1 static char help[] = "Tests for mesh interpolation\n\n";
2
3 /* TODO
4 */
5
6 #include <petscdmplex.h>
7
8 /* List of test meshes
9
10 Triangle
11 --------
12 Test 0:
13 Two triangles sharing a face
14
15 4
16 / | \
17 / | \
18 / | \
19 2 0 | 1 5
20 \ | /
21 \ | /
22 \ | /
23 3
24
25 should become
26
27 4
28 / | \
29 8 | 9
30 / | \
31 2 0 7 1 5
32 \ | /
33 6 | 10
34 \ | /
35 3
36
37 Tetrahedron
38 -----------
39 Test 0:
40 Two tets sharing a face
41
42 cell 5 _______ cell
43 0 / | \ \ 1
44 / | \ \
45 / | \ \
46 2----|----4-----6
47 \ | / /
48 \ | / /
49 \ | / /
50 3-------
51
52 should become
53
54 cell 5 _______ cell
55 0 / | \ \ 1
56 / | \ \
57 17 | 18 13 22
58 / 8 19 10 \ \
59 / | \ \
60 2---14-|------4--20--6
61 \ | / /
62 \ 9 | 7 / /
63 16 | 15 11 21
64 \ | / /
65 \ | / /
66 3-------
67
68
69 Quadrilateral
70 -------------
71 Test 0:
72 Two quads sharing a face
73
74 5-------4-------7
75 | | |
76 | 0 | 1 |
77 | | |
78 2-------3-------6
79
80 should become
81
82 5--10---4--14---7
83 | | |
84 11 0 9 1 13
85 | | |
86 2---8---3--12---6
87
88 Test 1:
89 A quad and a triangle sharing a face
90
91 5-------4
92 | | \
93 | 0 | \
94 | | 1 \
95 2-------3----6
96
97 should become
98
99 5---9---4
100 | | \
101 10 0 8 12
102 | | 1 \
103 2---7---3-11-6
104
105 Hexahedron
106 ----------
107 Test 0:
108 Two hexes sharing a face
109
110 cell 9-------------8-------------13 cell
111 0 /| /| /| 1
112 / | 15 / | 21 / |
113 / | / | / |
114 6-------------7-------------12 |
115 | | 18 | | 24 | |
116 | | | | | |
117 |19 | |17 | |23 |
118 | | 16 | | 22 | |
119 | 5---------|---4---------|---11
120 | / | / | /
121 | / 14 | / 20 | /
122 |/ |/ |/
123 2-------------3-------------10
124
125 should become
126
127 cell 9-----31------8-----42------13 cell
128 0 /| /| /| 1
129 32 | 15 30| 21 41|
130 / | / | / |
131 6-----29------7-----40------12 |
132 | | 17 | | 23 | |
133 | 35 | 36 | 44
134 |19 | |18 | |24 |
135 34 | 16 33 | 22 43 |
136 | 5-----26--|---4-----37--|---11
137 | / | / | /
138 | 25 14 | 27 20 | 38
139 |/ |/ |/
140 2-----28------3-----39------10
141
142 */
143
144 typedef struct {
145 DM dm;
146 PetscInt debug; /* The debugging level */
147 PetscInt testNum; /* Indicates the mesh to create */
148 PetscInt dim; /* The topological mesh dimension */
149 PetscBool cellSimplex; /* Use simplices or hexes */
150 PetscBool useGenerator; /* Construct mesh with a mesh generator */
151 char filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
152 } AppCtx;
153
ProcessOptions(MPI_Comm comm,AppCtx * options)154 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
155 {
156 PetscErrorCode ierr;
157
158 PetscFunctionBegin;
159 options->debug = 0;
160 options->testNum = 0;
161 options->dim = 2;
162 options->cellSimplex = PETSC_TRUE;
163 options->useGenerator = PETSC_FALSE;
164 options->filename[0] = '\0';
165
166 ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr);
167 ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex7.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr);
168 ierr = PetscOptionsBoundedInt("-testnum", "The mesh to create", "ex7.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr);
169 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex7.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr);
170 ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex7.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr);
171 ierr = PetscOptionsBool("-use_generator", "Use a mesh generator to build the mesh", "ex7.c", options->useGenerator, &options->useGenerator, NULL);CHKERRQ(ierr);
172 ierr = PetscOptionsString("-filename", "The mesh file", "ex7.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr);
173 ierr = PetscOptionsEnd();
174 PetscFunctionReturn(0);
175 }
176
CreateSimplex_2D(MPI_Comm comm,DM dm)177 PetscErrorCode CreateSimplex_2D(MPI_Comm comm, DM dm)
178 {
179 PetscInt depth = 1, testNum = 0, p;
180 PetscMPIInt rank;
181 PetscErrorCode ierr;
182
183 PetscFunctionBegin;
184 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
185 if (!rank) {
186 switch (testNum) {
187 case 0:
188 {
189 PetscInt numPoints[2] = {4, 2};
190 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0};
191 PetscInt cones[6] = {2, 3, 4, 5, 4, 3};
192 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
193 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
194 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1};
195
196 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
197 for (p = 0; p < 4; ++p) {
198 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);
199 }
200 }
201 break;
202 default:
203 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
204 }
205 } else {
206 PetscInt numPoints[2] = {0, 0};
207
208 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
209 }
210 PetscFunctionReturn(0);
211 }
212
CreateSimplex_3D(MPI_Comm comm,DM dm)213 PetscErrorCode CreateSimplex_3D(MPI_Comm comm, DM dm)
214 {
215 PetscInt depth = 1, testNum = 0, p;
216 PetscMPIInt rank;
217 PetscErrorCode ierr;
218
219 PetscFunctionBegin;
220 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
221 if (!rank) {
222 switch (testNum) {
223 case 0:
224 {
225 PetscInt numPoints[2] = {5, 2};
226 PetscInt coneSize[23] = {4, 4, 0, 0, 0, 0, 0};
227 PetscInt cones[8] = {2, 4, 3, 5, 3, 4, 6, 5};
228 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
229 PetscScalar vertexCoords[15] = {0.0, 0.0, -0.5, 0.0, -0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5};
230 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1};
231
232 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
233 for (p = 0; p < 4; ++p) {
234 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);
235 }
236 }
237 break;
238 default:
239 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
240 }
241 } else {
242 PetscInt numPoints[2] = {0, 0};
243
244 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
245 }
246 PetscFunctionReturn(0);
247 }
248
CreateQuad_2D(MPI_Comm comm,PetscInt testNum,DM dm)249 PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM dm)
250 {
251 PetscInt depth = 1, p;
252 PetscMPIInt rank;
253 PetscErrorCode ierr;
254
255 PetscFunctionBegin;
256 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
257 if (!rank) {
258 switch (testNum) {
259 case 0:
260 {
261 PetscInt numPoints[2] = {6, 2};
262 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0};
263 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4};
264 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
265 PetscScalar vertexCoords[12] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0, -0.5, 1.0, 0.5, 0.0, 0.5, 1.0};
266 PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1};
267
268 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
269 for (p = 0; p < 6; ++p) {
270 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);
271 }
272 }
273 break;
274 case 1:
275 {
276 PetscInt numPoints[2] = {5, 2};
277 PetscInt coneSize[7] = {4, 3, 0, 0, 0, 0, 0};
278 PetscInt cones[7] = {2, 3, 4, 5, 3, 6, 4};
279 PetscInt coneOrientations[7] = {0, 0, 0, 0, 0, 0, 0};
280 PetscScalar vertexCoords[10] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0, -0.5, 1.0, 0.5, 0.0};
281 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
282
283 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
284 for (p = 0; p < 5; ++p) {
285 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);
286 }
287 }
288 break;
289 default:
290 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
291 }
292 } else {
293 PetscInt numPoints[2] = {0, 0};
294
295 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
296 }
297 PetscFunctionReturn(0);
298 }
299
CreateHex_3D(MPI_Comm comm,DM dm)300 PetscErrorCode CreateHex_3D(MPI_Comm comm, DM dm)
301 {
302 PetscInt depth = 1, testNum = 0, p;
303 PetscMPIInt rank;
304 PetscErrorCode ierr;
305
306 PetscFunctionBegin;
307 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
308 if (!rank) {
309 switch (testNum) {
310 case 0:
311 {
312 PetscInt numPoints[2] = {12, 2};
313 PetscInt coneSize[14] = {8, 8, 0,0,0,0,0,0,0,0,0,0,0,0};
314 PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8};
315 PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0, 0,0, 0, 0,0};
316 PetscScalar vertexCoords[36] = {-0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0,
317 -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0,
318 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0};
319 PetscInt markerPoints[16] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1};
320
321 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
322 for (p = 0; p < 8; ++p) {
323 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);
324 }
325 }
326 break;
327 default:
328 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
329 }
330 } else {
331 PetscInt numPoints[2] = {0, 0};
332
333 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
334 }
335 PetscFunctionReturn(0);
336 }
337
CheckMesh(DM dm,AppCtx * user)338 PetscErrorCode CheckMesh(DM dm, AppCtx *user)
339 {
340 PetscReal detJ, J[9], refVol = 1.0;
341 PetscReal vol;
342 PetscInt dim, depth, d, cStart, cEnd, c;
343 PetscErrorCode ierr;
344
345 PetscFunctionBegin;
346 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
347 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
348 for (d = 0; d < dim; ++d) {
349 refVol *= 2.0;
350 }
351 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
352 for (c = cStart; c < cEnd; ++c) {
353 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, NULL, J, NULL, &detJ);CHKERRQ(ierr);
354 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh cell %d is inverted, |J| = %g", c, detJ);
355 if (user->debug) {PetscPrintf(PETSC_COMM_SELF, "FEM Volume: %g\n", detJ*refVol);CHKERRQ(ierr);}
356 if (depth > 1) {
357 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr);
358 if (vol <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh cell %d is inverted, vol = %g", c, vol);
359 if (user->debug) {PetscPrintf(PETSC_COMM_SELF, "FVM Volume: %g\n", vol);CHKERRQ(ierr);}
360 }
361 }
362 PetscFunctionReturn(0);
363 }
364
CompareCones(DM dm,DM idm)365 PetscErrorCode CompareCones(DM dm, DM idm)
366 {
367 PetscInt cStart, cEnd, c, vStart, vEnd, v;
368 PetscErrorCode ierr;
369
370 PetscFunctionBegin;
371 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
372 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
373 for (c = cStart; c < cEnd; ++c) {
374 const PetscInt *cone;
375 PetscInt *points = NULL, numPoints, p, numVertices = 0, coneSize;
376
377 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
378 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr);
379 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
380 for (p = 0; p < numPoints*2; p += 2) {
381 const PetscInt point = points[p];
382 if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
383 }
384 if (numVertices != coneSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In cell %d, cone size %d != %d vertices in closure", c, coneSize, numVertices);
385 for (v = 0; v < numVertices; ++v) {
386 if (cone[v] != points[v]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In cell %d, cone point %d is %d != %d vertex in closure", c, v, cone[v], points[v]);
387 }
388 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
389 }
390 PetscFunctionReturn(0);
391 }
392
CreateMesh(MPI_Comm comm,PetscInt testNum,AppCtx * user,DM * dm)393 PetscErrorCode CreateMesh(MPI_Comm comm, PetscInt testNum, AppCtx *user, DM *dm)
394 {
395 PetscInt dim = user->dim;
396 PetscBool cellSimplex = user->cellSimplex;
397 PetscBool useGenerator = user->useGenerator;
398 const char *filename = user->filename;
399 size_t len;
400 PetscMPIInt rank;
401 PetscErrorCode ierr;
402
403 PetscFunctionBegin;
404 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
405 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
406 if (len) {
407 ierr = DMPlexCreateFromFile(comm, filename, PETSC_FALSE, dm);CHKERRQ(ierr);
408 ierr = DMGetDimension(*dm, &dim);CHKERRQ(ierr);
409 } else if (useGenerator) {
410 ierr = DMPlexCreateBoxMesh(comm, dim, cellSimplex, NULL, NULL, NULL, NULL, PETSC_FALSE, dm);CHKERRQ(ierr);
411 } else {
412 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
413 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
414 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
415 switch (dim) {
416 case 2:
417 if (cellSimplex) {
418 ierr = CreateSimplex_2D(comm, *dm);CHKERRQ(ierr);
419 } else {
420 ierr = CreateQuad_2D(comm, testNum, *dm);CHKERRQ(ierr);
421 }
422 break;
423 case 3:
424 if (cellSimplex) {
425 ierr = CreateSimplex_3D(comm, *dm);CHKERRQ(ierr);
426 } else {
427 ierr = CreateHex_3D(comm, *dm);CHKERRQ(ierr);
428 }
429 break;
430 default:
431 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
432 }
433 }
434 {
435 DM idm;
436
437 ierr = CheckMesh(*dm, user);CHKERRQ(ierr);
438 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
439 ierr = CompareCones(*dm, idm);CHKERRQ(ierr);
440 ierr = DMDestroy(dm);CHKERRQ(ierr);
441 *dm = idm;
442 }
443 {
444 DM distributedMesh = NULL;
445 PetscPartitioner part;
446
447 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr);
448 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
449 /* Distribute mesh over processes */
450 ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr);
451 if (distributedMesh) {
452 ierr = DMViewFromOptions(distributedMesh, NULL, "-dm_view");CHKERRQ(ierr);
453 ierr = DMDestroy(dm);CHKERRQ(ierr);
454 *dm = distributedMesh;
455 }
456 }
457 ierr = PetscObjectSetName((PetscObject) *dm, "Interpolated Mesh");CHKERRQ(ierr);
458 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
459 user->dm = *dm;
460 PetscFunctionReturn(0);
461 }
462
main(int argc,char ** argv)463 int main(int argc, char **argv)
464 {
465 AppCtx user; /* user-defined work context */
466 PetscErrorCode ierr;
467
468 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
469 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
470 ierr = CreateMesh(PETSC_COMM_WORLD, user.testNum, &user, &user.dm);CHKERRQ(ierr);
471 ierr = CheckMesh(user.dm, &user);CHKERRQ(ierr);
472 ierr = DMDestroy(&user.dm);CHKERRQ(ierr);
473 ierr = PetscFinalize();
474 return ierr;
475 }
476
477 /*TEST
478
479 # Two cell test meshes 0-7
480 test:
481 suffix: 0
482 args: -dim 2 -dm_view ascii::ascii_info_detail
483 test:
484 suffix: 1
485 nsize: 2
486 args: -petscpartitioner_type simple -dim 2 -dm_view ascii::ascii_info_detail
487 test:
488 suffix: 2
489 args: -dim 2 -cell_simplex 0 -dm_view ascii::ascii_info_detail
490 test:
491 suffix: 3
492 nsize: 2
493 args: -petscpartitioner_type simple -dim 2 -cell_simplex 0 -dm_view ascii::ascii_info_detail
494 test:
495 suffix: 4
496 args: -dim 3 -dm_view ascii::ascii_info_detail
497 test:
498 suffix: 5
499 nsize: 2
500 args: -petscpartitioner_type simple -dim 3 -dm_view ascii::ascii_info_detail
501 test:
502 suffix: 6
503 args: -dim 3 -cell_simplex 0 -dm_view ascii::ascii_info_detail
504 test:
505 suffix: 7
506 nsize: 2
507 args: -petscpartitioner_type simple -dim 3 -cell_simplex 0 -dm_view ascii::ascii_info_detail
508 # 2D Hybrid Mesh 8
509 test:
510 suffix: 8
511 args: -dim 2 -cell_simplex 0 -testnum 1 -dm_view ascii::ascii_info_detail
512 # TetGen meshes 9-10
513 test:
514 suffix: 9
515 requires: triangle
516 args: -dim 2 -use_generator -dm_view ascii::ascii_info_detail
517 test:
518 suffix: 10
519 requires: ctetgen
520 args: -dim 3 -use_generator -dm_view ascii::ascii_info_detail
521 # Cubit meshes 11
522 test:
523 suffix: 11
524 requires: exodusii
525 args: -dim 3 -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo -dm_view ascii::ascii_info_detail
526
527 TEST*/
528