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;
30: PetscFunctionBeginUser;
31: PetscCall(DMGetCoordinateDM(dm, &cdm));
32: PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
33: PetscCall(PetscFEGetBasisSpace(fe, &P));
34: PetscCall(PetscFEGetDualSpace(fe, &Q));
35: PetscCall(PetscDualSpaceGetDM(Q, &K));
36: PetscCall(DMGetDimension(K, &dim));
37: PetscCall(PetscSpaceGetDegree(P, &k, NULL));
38: PetscCall(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: PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, isSimplex, k, PETSC_DETERMINE, &fe));
50: PetscCall(PetscFESetName(fe, "scalar"));
51: PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
52: PetscCall(PetscFEDestroy(&fe));
53: PetscCall(DMCreateDS(dm));
55: PetscCall(DMGetDS(dm, &ds));
56: PetscCall(PetscDSSetObjective(ds, 0, one));
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode CheckIntegral(DM dm, PetscReal integral, PetscReal tol)
61: {
62: Vec u;
63: PetscReal rval;
64: PetscScalar result;
66: PetscFunctionBeginUser;
67: PetscCall(DMGetGlobalVector(dm, &u));
68: PetscCall(DMPlexComputeIntegralFEM(dm, u, &result, NULL));
69: PetscCall(DMRestoreGlobalVector(dm, &u));
70: rval = PetscRealPart(result);
71: if (integral > 0 && PetscAbsReal(integral - rval) > tol) {
72: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
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;
97: PetscFunctionBeginUser;
98: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
100: PetscCall(PetscStrncpy(geodir, "${PETSC_DIR}/share/petsc/datafiles/meshes", sizeof(geodir)));
101: PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "GMSH", path, sizeof(path), &set));
102: if (set) PetscCall(PetscStrncpy(gmsh, path, sizeof(gmsh)));
103: PetscCall(PetscOptionsGetString(NULL, NULL, "-gmsh", gmsh, sizeof(gmsh), NULL));
104: PetscCall(PetscOptionsGetString(NULL, NULL, "-dir", geodir, sizeof(geodir), NULL));
105: PetscCall(PetscOptionsGetString(NULL, NULL, "-out", outdir, sizeof(outdir), NULL));
106: PetscCall(PetscOptionsGetEList(NULL, NULL, "-msh", mshlist, PETSC_STATIC_ARRAY_LENGTH(mshlist), &msh, NULL));
107: PetscCall(PetscOptionsGetEList(NULL, NULL, "-fmt", fmtlist, PETSC_STATIC_ARRAY_LENGTH(fmtlist), &fmt, NULL));
108: PetscCall(PetscOptionsGetBool(NULL, NULL, "-bin", &bin, NULL));
109: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
110: PetscCall(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: PetscCall(PetscSNPrintf(cmd, sizeof(cmd), "%s -info", gmsh));
117: PetscCall(PetscPOpen(PETSC_COMM_SELF, NULL, cmd, "r", &fp));
118: if (fp) inum = fscanf(fp, "Version %s %d.%d.%d", space, &major, &minor, µ);
119: PetscCall(PetscPClose(PETSC_COMM_SELF, fp));
120: if (inum != 4 || major < 4 || (major == 4 && minor < 2)) {
121: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Gmsh>=4.2.0 not available\n"));
122: goto finish;
123: }
124: }
126: PetscCall(PetscSNPrintf(tag, sizeof(tag), "%s-%d-%d-%s%s", mshlist[msh], (int)dim, (int)order, fmtlist[fmt], bin ? "-bin" : ""));
127: PetscCall(PetscSNPrintf(geo, sizeof(geo), "%s/gmsh-%s.geo", geodir, mshlist[msh]));
128: PetscCall(PetscSNPrintf(out, sizeof(out), "%s/mesh-%s.msh", outdir, tag));
129: PetscCall(PetscStrreplace(PETSC_COMM_SELF, geo, path, sizeof(path)));
130: PetscCall(PetscFixFilename(path, geo));
131: PetscCall(PetscStrreplace(PETSC_COMM_SELF, out, path, sizeof(path)));
132: PetscCall(PetscFixFilename(path, out));
133: PetscCall(PetscTestFile(geo, 'r', &flg));
134: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "File not found: %s", geo);
136: PetscCall(PetscSNPrintf(cmd, sizeof(cmd), cmdtemplate, gmsh, fmtlist[fmt], bin ? "-bin" : "", (int)dim, (int)order, geo, out));
137: PetscCall(PetscPOpen(PETSC_COMM_SELF, NULL, cmd, "r", &fp));
138: PetscCall(PetscPClose(PETSC_COMM_SELF, fp));
140: PetscCall(DMPlexCreateFromFile(PETSC_COMM_SELF, out, "ex99_plex", PETSC_TRUE, &dm));
141: PetscCall(PetscSNPrintf(tag, sizeof(tag), "mesh-%s", mshlist[msh]));
142: PetscCall(PetscObjectSetName((PetscObject)dm, tag));
143: PetscCall(DMSetFromOptions(dm));
144: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
145: {
146: PetscBool check;
147: PetscReal integral = 0, tol = (PetscReal)1.0e-4;
148: PetscCall(PetscOptionsGetReal(NULL, NULL, "-integral", &integral, &check));
149: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));
150: if (check) {
151: PetscCall(CreateFE(dm));
152: PetscCall(CheckIntegral(dm, integral, tol));
153: }
154: }
155: PetscCall(DMDestroy(&dm));
157: finish:
158: PetscCall(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*/