Actual source code: ex99.c

  1: static char help[] = "Tests DMPlex Gmsh reader.\n\n";

  3: #include <petscdmplex.h>

  5: #if !defined(PETSC_GMSH_EXE)
  6:   #define PETSC_GMSH_EXE "gmsh"
  7: #endif

  9: #include <petscds.h>

 11: static void one(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar value[])
 12: {
 13:   value[0] = (PetscReal)1;
 14: }

 16: static PetscErrorCode CreateFE(DM dm)
 17: {
 18:   DM             cdm;
 19:   PetscSpace     P;
 20:   PetscDualSpace Q;
 21:   DM             K;
 22:   PetscFE        fe;
 23:   DMPolytopeType ptype;

 25:   PetscInt  dim, k;
 26:   PetscBool isSimplex;

 28:   PetscDS ds;

 31:   DMGetCoordinateDM(dm, &cdm);
 32:   DMGetField(cdm, 0, NULL, (PetscObject *)&fe);
 33:   PetscFEGetBasisSpace(fe, &P);
 34:   PetscFEGetDualSpace(fe, &Q);
 35:   PetscDualSpaceGetDM(Q, &K);
 36:   DMGetDimension(K, &dim);
 37:   PetscSpaceGetDegree(P, &k, NULL);
 38:   DMPlexGetCellType(K, 0, &ptype);
 39:   switch (ptype) {
 40:   case DM_POLYTOPE_QUADRILATERAL:
 41:   case DM_POLYTOPE_HEXAHEDRON:
 42:     isSimplex = PETSC_FALSE;
 43:     break;
 44:   default:
 45:     isSimplex = PETSC_TRUE;
 46:     break;
 47:   }

 49:   PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, isSimplex, k, PETSC_DETERMINE, &fe);
 50:   PetscFESetName(fe, "scalar");
 51:   DMAddField(dm, NULL, (PetscObject)fe);
 52:   PetscFEDestroy(&fe);
 53:   DMCreateDS(dm);

 55:   DMGetDS(dm, &ds);
 56:   PetscDSSetObjective(ds, 0, one);
 57:   return 0;
 58: }

 60: static PetscErrorCode CheckIntegral(DM dm, PetscReal integral, PetscReal tol)
 61: {
 62:   Vec         u;
 63:   PetscReal   rval;
 64:   PetscScalar result;

 67:   DMGetGlobalVector(dm, &u);
 68:   DMPlexComputeIntegralFEM(dm, u, &result, NULL);
 69:   DMRestoreGlobalVector(dm, &u);
 70:   rval = PetscRealPart(result);
 71:   if (integral > 0 && PetscAbsReal(integral - rval) > tol) {
 72:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Calculated value %g != %g actual value (error %g > %g tol)\n", (double)rval, (double)integral, (double)PetscAbsReal(integral - rval), (double)tol);
 73:   }
 74:   return 0;
 75: }

 77: int main(int argc, char **argv)
 78: {
 79:   DM                dm;
 80:   const char *const mshlist[] = {"seg", "tri", "qua", "tet", "wed", "hex", "vtx", "B2tri", "B2qua", "B3tet", "B3hex"};
 81:   const char *const fmtlist[] = {"msh22", "msh40", "msh41"};
 82:   PetscInt          msh       = 5;
 83:   PetscInt          fmt       = 2;
 84:   PetscBool         bin       = PETSC_TRUE;
 85:   PetscInt          dim       = 3;
 86:   PetscInt          order     = 2;

 88:   const char cmdtemplate[]            = "%s -format %s %s -%d -order %d %s -o %s";
 89:   char       gmsh[PETSC_MAX_PATH_LEN] = PETSC_GMSH_EXE;
 90:   char       tag[PETSC_MAX_PATH_LEN], path[PETSC_MAX_PATH_LEN];
 91:   char       geo[PETSC_MAX_PATH_LEN], geodir[PETSC_MAX_PATH_LEN] = ".";
 92:   char       out[PETSC_MAX_PATH_LEN], outdir[PETSC_MAX_PATH_LEN] = ".";
 93:   char       cmd[PETSC_MAX_PATH_LEN * 4];
 94:   PetscBool  set, flg;
 95:   FILE      *fp;

 98:   PetscInitialize(&argc, &argv, NULL, help);

100:   PetscStrncpy(geodir, "${PETSC_DIR}/share/petsc/datafiles/meshes", sizeof(geodir));
101:   PetscOptionsGetenv(PETSC_COMM_SELF, "GMSH", path, sizeof(path), &set);
102:   if (set) PetscStrncpy(gmsh, path, sizeof(gmsh));
103:   PetscOptionsGetString(NULL, NULL, "-gmsh", gmsh, sizeof(gmsh), NULL);
104:   PetscOptionsGetString(NULL, NULL, "-dir", geodir, sizeof(geodir), NULL);
105:   PetscOptionsGetString(NULL, NULL, "-out", outdir, sizeof(outdir), NULL);
106:   PetscOptionsGetEList(NULL, NULL, "-msh", mshlist, PETSC_STATIC_ARRAY_LENGTH(mshlist), &msh, NULL);
107:   PetscOptionsGetEList(NULL, NULL, "-fmt", fmtlist, PETSC_STATIC_ARRAY_LENGTH(fmtlist), &fmt, NULL);
108:   PetscOptionsGetBool(NULL, NULL, "-bin", &bin, NULL);
109:   PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);
110:   PetscOptionsGetInt(NULL, NULL, "-order", &order, NULL);
111:   if (fmt == 1) bin = PETSC_FALSE; /* Recent Gmsh releases cannot generate msh40+binary format*/

113:   { /* This test requires Gmsh >= 4.2.0 */
114:     char space[PETSC_MAX_PATH_LEN];
115:     int  inum = 0, major = 0, minor = 0, micro = 0;
116:     PetscSNPrintf(cmd, sizeof(cmd), "%s -info", gmsh);
117:     PetscPOpen(PETSC_COMM_SELF, NULL, cmd, "r", &fp);
118:     if (fp) inum = fscanf(fp, "Version %s %d.%d.%d", space, &major, &minor, &micro);
119:     PetscPClose(PETSC_COMM_SELF, fp);
120:     if (inum != 4 || major < 4 || (major == 4 && minor < 2)) {
121:       PetscPrintf(PETSC_COMM_SELF, "Gmsh>=4.2.0 not available\n");
122:       goto finish;
123:     }
124:   }

126:   PetscSNPrintf(tag, sizeof(tag), "%s-%d-%d-%s%s", mshlist[msh], (int)dim, (int)order, fmtlist[fmt], bin ? "-bin" : "");
127:   PetscSNPrintf(geo, sizeof(geo), "%s/gmsh-%s.geo", geodir, mshlist[msh]);
128:   PetscSNPrintf(out, sizeof(out), "%s/mesh-%s.msh", outdir, tag);
129:   PetscStrreplace(PETSC_COMM_SELF, geo, path, sizeof(path));
130:   PetscFixFilename(path, geo);
131:   PetscStrreplace(PETSC_COMM_SELF, out, path, sizeof(path));
132:   PetscFixFilename(path, out);
133:   PetscTestFile(geo, 'r', &flg);

136:   PetscSNPrintf(cmd, sizeof(cmd), cmdtemplate, gmsh, fmtlist[fmt], bin ? "-bin" : "", (int)dim, (int)order, geo, out);
137:   PetscPOpen(PETSC_COMM_SELF, NULL, cmd, "r", &fp);
138:   PetscPClose(PETSC_COMM_SELF, fp);

140:   DMPlexCreateFromFile(PETSC_COMM_SELF, out, "ex99_plex", PETSC_TRUE, &dm);
141:   PetscSNPrintf(tag, sizeof(tag), "mesh-%s", mshlist[msh]);
142:   PetscObjectSetName((PetscObject)dm, tag);
143:   DMSetFromOptions(dm);
144:   DMViewFromOptions(dm, NULL, "-dm_view");
145:   {
146:     PetscBool check;
147:     PetscReal integral = 0, tol = (PetscReal)1.0e-4;
148:     PetscOptionsGetReal(NULL, NULL, "-integral", &integral, &check);
149:     PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL);
150:     if (check) {
151:       CreateFE(dm);
152:       CheckIntegral(dm, integral, tol);
153:     }
154:   }
155:   DMDestroy(&dm);

157: finish:
158:   PetscFinalize();
159:   return 0;
160: }

162: /*TEST

164:   build:
165:     requires: defined(PETSC_HAVE_POPEN)

167:   test:
168:     requires: defined(PETSC_GMSH_EXE)
169:     args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes
170:     args: -msh {{vtx}separate_output}
171:     args: -order 1
172:     args: -fmt {{msh22 msh40 msh41}} -bin {{0 1}}
173:     args: -dm_view ::ascii_info_detail
174:     args: -dm_plex_check_all
175:     args: -dm_plex_gmsh_highorder false
176:     args: -dm_plex_boundary_label marker
177:     args: -dm_plex_gmsh_spacedim 3

179:   test:
180:     requires: defined(PETSC_GMSH_EXE)
181:     args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes
182:     args: -msh {{seg tri qua tet wed hex}separate_output}
183:     args: -order {{1 2 3 7}}
184:     args: -fmt {{msh22 msh40 msh41}} -bin {{0 1}}
185:     args: -dm_view ::ascii_info_detail
186:     args: -dm_plex_check_all
187:     args: -dm_plex_gmsh_highorder false
188:     args: -dm_plex_boundary_label marker

190:   testset:
191:     suffix: B2 # 2D ball
192:     requires: defined(PETSC_GMSH_EXE)
193:     args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes
194:     args: -msh {{B2tri B2qua}}
195:     args: -dim 2 -integral 3.141592653589793 # pi
196:     args: -order {{2 3 4 5 6 7 8 9}} -tol 0.05

198:   testset:
199:     suffix: B2_bnd # 2D ball boundary
200:     requires: defined(PETSC_GMSH_EXE)
201:     args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes
202:     args: -dm_plex_gmsh_spacedim 2
203:     args: -msh {{B2tri B2qua}}
204:     args: -dim 1 -integral 6.283185307179586 # 2*pi
205:     args: -order {{2 3 4 5 6 7 8 9}} -tol 0.05

207:   testset:
208:     suffix: B3 # 3D ball
209:     requires: defined(PETSC_GMSH_EXE)
210:     args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes
211:     args: -msh {{B3tet B3hex}}
212:     args: -dim 3 -integral 4.1887902047863905 # 4/3*pi
213:     args: -order {{2 3 4 5}} -tol 0.20

215:   testset:
216:     suffix: B3_bnd # 3D ball boundary
217:     requires: defined(PETSC_GMSH_EXE)
218:     args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes
219:     args: -dm_plex_gmsh_spacedim 3
220:     args: -msh {{B3tet B3hex}}
221:     args: -dim 2 -integral 12.566370614359172 # 4*pi
222:     args: -order {{2 3 4 5 6 7 8 9}} -tol 0.25

224:   testset:
225:     suffix: B_lin # linear discretizations
226:     requires: defined(PETSC_GMSH_EXE)
227:     args: -dir ${wPETSC_DIR}/share/petsc/datafiles/meshes
228:     args: -dm_plex_gmsh_highorder true
229:     args: -dm_plex_gmsh_project true
230:     args: -dm_plex_gmsh_project_petscspace_degree {{1 2 3}separate_output}
231:     args: -dm_plex_gmsh_fe_view
232:     args: -dm_plex_gmsh_project_fe_view
233:     args: -order 1 -tol 1e-4
234:     test:
235:       suffix: dim-1
236:       args: -dm_plex_gmsh_spacedim 2
237:       args: -msh {{B2tri B2qua}separate_output}
238:       args: -dim 1 -integral 5.656854249492381 # 4*sqrt(2)
239:     test:
240:       suffix: dim-2
241:       args: -dm_plex_gmsh_spacedim 2
242:       args: -msh {{B2tri B2qua}separate_output}
243:       args: -dim 2 -integral 2.000000000000000 # 2
244:     test:
245:       suffix: dim-2_msh-B3tet
246:       args: -dm_plex_gmsh_spacedim 3
247:       args: -msh B3tet -dim 2 -integral 9.914478
248:     test:
249:       suffix: dim-2_msh-B3hex
250:       args: -dm_plex_gmsh_spacedim 3
251:       args: -msh B3hex -dim 2 -integral 8.000000
252:     test:
253:       suffix: dim-3_msh-B3tet
254:       args: -dm_plex_gmsh_spacedim 3
255:       args: -msh B3tet -dim 3 -integral 2.666649
256:     test:
257:       suffix: dim-3_msh-B3hex
258:       args: -dm_plex_gmsh_spacedim 3
259:       args: -msh B3hex -dim 3 -integral 1.539600

261: TEST*/