xref: /petsc/src/dm/impls/plex/tests/ex4.c (revision 70a7d78aacfbd24b2e31399a3d8e056944bb7de3) !
1 static char help[] = "Tests for refinement of meshes created by hand\n\n";
2 
3 #include <petscdmplex.h>
4 
5 typedef struct {
6   PetscInt  debug;          /* The debugging level */
7   PetscInt  dim;            /* The topological mesh dimension */
8   PetscBool cellHybrid;     /* Use a hybrid mesh */
9   PetscBool cellSimplex;    /* Use simplices or hexes */
10   PetscBool testPartition;  /* Use a fixed partitioning for testing */
11   PetscInt  testNum;        /* The particular mesh to test */
12   PetscBool uninterpolate;  /* Uninterpolate the mesh at the end */
13   PetscBool reinterpolate;  /* Reinterpolate the mesh at the end */
14 } AppCtx;
15 
16 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
17 {
18   PetscErrorCode ierr;
19 
20   PetscFunctionBegin;
21   options->debug          = 0;
22   options->dim            = 2;
23   options->cellHybrid     = PETSC_TRUE;
24   options->cellSimplex    = PETSC_TRUE;
25   options->testPartition  = PETSC_TRUE;
26   options->testNum        = 0;
27   options->uninterpolate  = PETSC_FALSE;
28   options->reinterpolate  = PETSC_FALSE;
29 
30   ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr);
31   ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex4.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr);
32   ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr);
33   ierr = PetscOptionsBool("-cell_hybrid", "Use a hybrid mesh", "ex4.c", options->cellHybrid, &options->cellHybrid, NULL);CHKERRQ(ierr);
34   ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex4.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr);
35   ierr = PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex4.c", options->testPartition, &options->testPartition, NULL);CHKERRQ(ierr);
36   ierr = PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex4.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr);
37   ierr = PetscOptionsBool("-uninterpolate", "Uninterpolate the mesh at the end", "ex4.c", options->uninterpolate, &options->uninterpolate, NULL);CHKERRQ(ierr);
38   ierr = PetscOptionsBool("-reinterpolate", "Reinterpolate the mesh at the end", "ex4.c", options->reinterpolate, &options->reinterpolate, NULL);CHKERRQ(ierr);
39   ierr = PetscOptionsEnd();
40   PetscFunctionReturn(0);
41 }
42 
43 /* Two segments
44 
45   2-------0-------3-------1-------4
46 
47 become
48 
49   4---0---7---1---5---2---8---3---6
50 
51 */
52 PetscErrorCode CreateSimplex_1D(MPI_Comm comm, DM *dm)
53 {
54   PetscInt       depth = 1;
55   PetscMPIInt    rank;
56   PetscErrorCode ierr;
57 
58   PetscFunctionBegin;
59   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
60   if (!rank) {
61     PetscInt    numPoints[2]         = {3, 2};
62     PetscInt    coneSize[5]          = {2, 2, 0, 0, 0};
63     PetscInt    cones[4]             = {2, 3,  3, 4};
64     PetscInt    coneOrientations[16] = {0, 0,  0, 0};
65     PetscScalar vertexCoords[3]      = {-1.0, 0.0, 1.0};
66 
67     ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
68   } else {
69     PetscInt numPoints[2] = {0, 0};
70 
71     ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
72   }
73   PetscFunctionReturn(0);
74 }
75 
76 /* Two triangles
77         4
78       / | \
79      8  |  10
80     /   |   \
81    2  0 7  1 5
82     \   |   /
83      6  |  9
84       \ | /
85         3
86 
87 Becomes
88            10
89           / | \
90         21  |  26
91         /   |   \
92       14 2 20 4  16
93       /|\   |   /|\
94     22 | 28 | 32 | 25
95     /  |  \ | /  | 6\
96    8  29 3 13  7 31  11
97     \0 |  / | \  |  /
98     17 | 27 | 30 | 24
99       \|/   |   \|/
100       12 1 19 5  15
101         \   |   /
102         18  |  23
103           \ | /
104             9
105 */
106 PetscErrorCode CreateSimplex_2D(MPI_Comm comm, DM *dm)
107 {
108   PetscInt       depth = 2;
109   PetscMPIInt    rank;
110   PetscErrorCode ierr;
111 
112   PetscFunctionBegin;
113   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
114   if (!rank) {
115     PetscInt    numPoints[3]         = {4, 5, 2};
116     PetscInt    coneSize[11]         = {3, 3, 0, 0, 0, 0, 2, 2, 2, 2, 2};
117     PetscInt    cones[16]            = {6, 7, 8,  7, 9, 10,  2, 3,  3, 4,  4, 2,  3, 5,  5, 4};
118     PetscInt    coneOrientations[16] = {0, 0, 0, -2, 0,  0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0};
119     PetscScalar vertexCoords[8]      = {-0.5, 0.0,  0.0, -0.5,  0.0, 0.5,  0.5, 0.0};
120 
121     ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
122   } else {
123     PetscInt numPoints[3] = {0, 0, 0};
124 
125     ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
126   }
127   PetscFunctionReturn(0);
128 }
129 
130 /* Two triangles separated by a zero-volume cell with 4 vertices/2 edges
131         5--16--8
132       / |      | \
133     11  |      |  12
134     /   |      |   \
135    3  0 10  2 14 1  6
136     \   |      |   /
137      9  |      |  13
138       \ |      | /
139         4--15--7
140 */
141 PetscErrorCode CreateSimplexHybrid_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
142 {
143   DM             idm, hdm = NULL;
144   DMLabel        faultLabel, hybridLabel;
145   PetscInt       p;
146   PetscMPIInt    rank;
147   PetscErrorCode ierr;
148 
149   PetscFunctionBegin;
150   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
151   if (!rank) {
152     switch (testNum) {
153     case 0:
154     {
155       PetscInt    numPoints[2]        = {4, 2};
156       PetscInt    coneSize[6]         = {3, 3, 0, 0, 0, 0};
157       PetscInt    cones[6]            = {2, 3, 4,  5, 4, 3};
158       PetscInt    coneOrientations[6] = {0, 0, 0,  0, 0, 0};
159       PetscScalar vertexCoords[8]     = {-1.0, -0.5,  0.0, -0.5,  0.0, 0.5,  1.0, 0.5};
160       PetscInt    faultPoints[2]      = {3, 4};
161 
162       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
163       for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
164     }
165     break;
166     case 1:
167     {
168       PetscInt    numPoints[2]         = {5, 4};
169       PetscInt    coneSize[9]          = {3, 3, 3, 3, 0, 0, 0, 0, 0};
170       PetscInt    cones[12]            = {4, 5, 6,  6, 7, 4,  6, 5, 8,  6, 8, 7};
171       PetscInt    coneOrientations[12] = {0, 0, 0,  0, 0, 0,  0, 0, 0,  0, 0, 0};
172       PetscScalar vertexCoords[10]     = {-1.0, 0.0,  0.0, -1.0,  0.0, 0.0,  0.0, 1.0,  1.0, 0.0};
173       PetscInt    faultPoints[3]       = {5, 6, 7};
174 
175       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
176       for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
177     }
178     break;
179     default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
180     }
181     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
182     ierr = PetscObjectSetOptionsPrefix((PetscObject) idm, "in_");CHKERRQ(ierr);
183     ierr = DMSetFromOptions(idm);CHKERRQ(ierr);
184     ierr = DMViewFromOptions(idm, NULL, "-dm_view");CHKERRQ(ierr);
185     ierr = DMGetLabel(*dm, "fault", &faultLabel);CHKERRQ(ierr);
186     ierr = DMPlexCreateHybridMesh(idm, faultLabel, NULL, &hybridLabel, NULL, NULL, &hdm);CHKERRQ(ierr);
187     ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr);
188     ierr = DMDestroy(&idm);CHKERRQ(ierr);
189     ierr = DMDestroy(dm);CHKERRQ(ierr);
190     *dm  = hdm;
191   } else {
192     PetscInt numPoints[2] = {0, 0};
193 
194     ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
195     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
196     ierr = PetscObjectSetOptionsPrefix((PetscObject) idm, "in_");CHKERRQ(ierr);
197     ierr = DMSetFromOptions(idm);CHKERRQ(ierr);
198     ierr = DMViewFromOptions(idm, NULL, "-dm_view");CHKERRQ(ierr);
199     ierr = DMPlexCreateHybridMesh(idm, NULL, NULL, NULL, NULL, NULL, &hdm);CHKERRQ(ierr);
200     ierr = DMDestroy(&idm);CHKERRQ(ierr);
201     ierr = DMDestroy(dm);CHKERRQ(ierr);
202     *dm  = hdm;
203   }
204   PetscFunctionReturn(0);
205 }
206 
207 /* Two quadrilaterals
208 
209   5----10-----4----14-----7
210   |           |           |
211   |           |           |
212   |           |           |
213  11     0     9     1     13
214   |           |           |
215   |           |           |
216   |           |           |
217   2-----8-----3----12-----6
218 */
219 PetscErrorCode CreateTensorProduct_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
220 {
221   PetscInt       depth = 2;
222   PetscMPIInt    rank;
223   PetscErrorCode ierr;
224 
225   PetscFunctionBegin;
226   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
227   if (!rank) {
228     PetscInt    numPoints[3]         = {6, 7, 2};
229     PetscInt    coneSize[15]         = {4, 4, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2};
230     PetscInt    cones[22]            = {8, 9, 10, 11,  12, 13, 14,  9,  2, 3,  3, 4,  4, 5,  5, 2,  3, 6,  6, 7,  7, 4};
231     PetscInt    coneOrientations[22] = {0, 0,  0,  0,   0,  0,  0, -2,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0};
232     PetscScalar vertexCoords[12]     = {-1.0, -0.5,  0.0, -0.5,  0.0, 0.5,  -1.0, 0.5,  1.0, -0.5,  1.0, 0.5};
233 
234     ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
235   } else {
236     PetscInt numPoints[3] = {0, 0, 0};
237 
238     ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
239   }
240   PetscFunctionReturn(0);
241 }
242 
243 PetscErrorCode CreateTensorProductHybrid_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
244 {
245   DM             idm, hdm = NULL;
246   DMLabel        faultLabel, hybridLabel;
247   PetscInt       p;
248   PetscMPIInt    rank;
249   PetscErrorCode ierr;
250 
251   PetscFunctionBegin;
252   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
253   if (!rank) {
254     PetscInt    numPoints[2]        = {6, 2};
255     PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
256     PetscInt    cones[8]            = {2, 3, 4, 5,  3, 6, 7, 4,};
257     PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
258     PetscScalar vertexCoords[12]    = {-1.0, -0.5,  0.0, -0.5,  0.0, 0.5,  -1.0, 0.5,  1.0, -0.5,  1.0, 0.5};
259     PetscInt    faultPoints[2]      = {3, 4};
260 
261     ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
262     for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
263     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
264     ierr = PetscObjectSetOptionsPrefix((PetscObject) idm, "in_");CHKERRQ(ierr);
265     ierr = DMSetFromOptions(idm);CHKERRQ(ierr);
266     ierr = DMViewFromOptions(idm, NULL, "-dm_view");CHKERRQ(ierr);
267     ierr = DMGetLabel(*dm, "fault", &faultLabel);CHKERRQ(ierr);
268     ierr = DMPlexCreateHybridMesh(idm, faultLabel, NULL, &hybridLabel, NULL, NULL, &hdm);CHKERRQ(ierr);
269     ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr);
270   } else {
271     PetscInt numPoints[3] = {0, 0, 0};
272 
273     ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
274     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
275     ierr = PetscObjectSetOptionsPrefix((PetscObject) idm, "in_");CHKERRQ(ierr);
276     ierr = DMSetFromOptions(idm);CHKERRQ(ierr);
277     ierr = DMViewFromOptions(idm, NULL, "-dm_view");CHKERRQ(ierr);
278     ierr = DMPlexCreateHybridMesh(idm, NULL, NULL, NULL, NULL, NULL, &hdm);CHKERRQ(ierr);
279   }
280   ierr = DMDestroy(&idm);CHKERRQ(ierr);
281   ierr = DMDestroy(dm);CHKERRQ(ierr);
282   *dm  = hdm;
283   PetscFunctionReturn(0);
284 }
285 
286 /* Two tetrahedrons
287 
288  cell   5          5______    cell
289  0    / | \        |\      \     1
290     17  |  18      | 18 13  21
291     /8 19 10\     19  \      \
292    2-14-|----4     |   4--22--6
293     \ 9 | 7 /      |10 /      /
294     16  |  15      | 15  12 20
295       \ | /        |/      /
296         3          3------
297 */
298 PetscErrorCode CreateSimplex_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
299 {
300   DM             idm;
301   PetscInt       depth = 3;
302   PetscMPIInt    rank;
303   PetscErrorCode ierr;
304 
305   PetscFunctionBegin;
306   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
307   if (!rank) {
308     switch (testNum) {
309     case 0:
310     {
311       PetscInt    numPoints[4]         = {5, 9, 7, 2};
312       PetscInt    coneSize[23]         = {4, 4, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2};
313       PetscInt    cones[47]            = { 7,  8,  9, 10,  10, 11, 12, 13,  14, 15, 16,  17, 18, 14,  16, 19, 17,  15, 18, 19,  20, 21, 19,  15, 22, 20,  18, 21, 22,  2, 4,  4, 3,  3, 2,  2, 5,  5, 4,  3, 5,  3, 6,  6, 5,  4, 6};
314       PetscInt    coneOrientations[47] = { 0,  0,  0,  0,  -3,  0,  0,  0,   0,  0,  0,   0,  0, -2,  -2,  0, -2,  -2, -2, -2,   0,  0, -2,  -2,  0, -2,  -2, -2, -2,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0};
315       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};
316 
317       ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
318     }
319     break;
320     case 1:
321     {
322       PetscInt    numPoints[2]        = {5, 2};
323       PetscInt    coneSize[7]         = {4, 4, 0, 0, 0, 0, 0};
324       PetscInt    cones[8]            = {4, 3, 5, 2,  5, 3, 4, 6};
325       PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
326       PetscScalar vertexCoords[15]    = {-1.0, 0.0, 0.0,  0.0, -1.0, 0.0,  0.0, 0.0, 1.0,  0.0, 1.0, 0.0,  1.0, 0.0, 0.0};
327 
328       depth = 1;
329       ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
330       ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
331       ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr);
332       ierr = DMDestroy(dm);CHKERRQ(ierr);
333       *dm  = idm;
334     }
335     break;
336     case 2:
337     {
338       PetscInt    numPoints[2]        = {4, 1};
339       PetscInt    coneSize[5]         = {4, 0, 0, 0, 0};
340       PetscInt    cones[4]            = {2, 3, 4, 1};
341       PetscInt    coneOrientations[4] = {0, 0, 0, 0};
342       PetscScalar vertexCoords[12]    = {0.0, 0.0, 0.0,  1.0, 0.0, 0.0,  0.0, 1.0, 0.0,  0.0, 0.0, 1.0};
343 
344       depth = 1;
345       ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
346       ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
347       ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr);
348       ierr = DMDestroy(dm);CHKERRQ(ierr);
349       *dm  = idm;
350     }
351     break;
352     default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
353     }
354   } else {
355     PetscInt numPoints[4] = {0, 0, 0, 0};
356 
357     ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
358     switch (testNum) {
359     case 1:
360       ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
361       ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr);
362       ierr = DMDestroy(dm);CHKERRQ(ierr);
363       *dm  = idm;
364       break;
365     }
366   }
367   PetscFunctionReturn(0);
368 }
369 
370 /* Two tetrahedrons separated by a zero-volume cell with 6 vertices
371 
372  cell   6 ___33___10______    cell
373  0    / | \        |\      \     1
374     21  |  23      | 29     27
375     /12 24 14\    30  \      \
376    3-20-|----5--32-|---9--26--7
377     \ 13| 11/      |18 /      /
378     19  |  22      | 28     25
379       \ | /        |/      /
380         4----31----8------
381          cell 2
382 */
383 PetscErrorCode CreateSimplexHybrid_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
384 {
385   DM             idm, hdm = NULL;
386   DMLabel        faultLabel, hybridLabel;
387   PetscInt       p;
388   PetscMPIInt    rank;
389   PetscErrorCode ierr;
390 
391   PetscFunctionBegin;
392   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
393   if (!rank) {
394     switch (testNum) {
395     case 0:
396     {
397       PetscInt    numPoints[2]        = {5, 2};
398       PetscInt    coneSize[7]         = {4, 4, 0, 0, 0, 0, 0};
399       PetscInt    cones[8]            = {4, 3, 5, 2,  5, 3, 4, 6};
400       PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
401       PetscScalar vertexCoords[15]    = {-1.0, 0.0, 0.0,  0.0, -1.0, 0.0,  0.0, 0.0, 1.0,  0.0, 1.0, 0.0,  1.0, 0.0, 0.0};
402       PetscInt    faultPoints[3]      = {3, 4, 5};
403 
404       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
405       for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
406     }
407     break;
408     case 1:
409     {
410       /* Tets 0,3,5 and 1,2,4 */
411       PetscInt    numPoints[2]         = {9, 6};
412       PetscInt    coneSize[15]         = {4, 4, 4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0};
413       PetscInt    cones[24]            = { 7,  8,  9, 6,  11, 13,  9, 14,  10, 11, 13, 9,
414                                            9, 10, 11, 7,   9, 14, 13, 12,   7,  8, 11, 9};
415       PetscInt    coneOrientations[24] = { 0, 0,  0, 0,   0,  0,  0,  0,   0,  0,  0, 0,
416                                            0, 0,  0, 0,   0,  0,  0,  0,   0,  0,  0, 0};
417       PetscScalar vertexCoords[27]     = {-2.0, -1.0,  0.0,  -2.0,  0.0,  0.0,  -2.0,  0.0,  1.0,
418                                            0.0, -1.0,  0.0,   0.0,  0.0,  0.0,   0.0,  0.0,  1.0,
419                                            2.0, -1.0,  0.0,   2.0,  0.0,  0.0,   2.0,  0.0,  1.0};
420       PetscInt    faultPoints[3]       = {9, 10, 11};
421 
422       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
423       for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
424     }
425     break;
426     default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
427     }
428     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
429     ierr = PetscObjectSetOptionsPrefix((PetscObject) idm, "in_");CHKERRQ(ierr);
430     ierr = DMSetFromOptions(idm);CHKERRQ(ierr);
431     ierr = DMViewFromOptions(idm, NULL, "-dm_view");CHKERRQ(ierr);
432     ierr = DMGetLabel(*dm, "fault", &faultLabel);CHKERRQ(ierr);
433     ierr = DMPlexCreateHybridMesh(idm, faultLabel, NULL, &hybridLabel, NULL, NULL, &hdm);CHKERRQ(ierr);
434     ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr);
435     ierr = DMDestroy(&idm);CHKERRQ(ierr);
436     ierr = DMDestroy(dm);CHKERRQ(ierr);
437     *dm  = hdm;
438   } else {
439     PetscInt numPoints[4] = {0, 0, 0, 0};
440 
441     ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
442     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
443     ierr = PetscObjectSetOptionsPrefix((PetscObject) idm, "in_");CHKERRQ(ierr);
444     ierr = DMSetFromOptions(idm);CHKERRQ(ierr);
445     ierr = DMViewFromOptions(idm, NULL, "-dm_view");CHKERRQ(ierr);
446     ierr = DMPlexCreateHybridMesh(idm, NULL, NULL, NULL, NULL, NULL, &hdm);CHKERRQ(ierr);
447     ierr = DMDestroy(&idm);CHKERRQ(ierr);
448     ierr = DMDestroy(dm);CHKERRQ(ierr);
449     *dm  = hdm;
450   }
451   PetscFunctionReturn(0);
452 }
453 
454 PetscErrorCode CreateTensorProduct_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
455 {
456   DM             idm;
457   PetscMPIInt    rank;
458   PetscErrorCode ierr;
459 
460   PetscFunctionBegin;
461   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
462   if (!rank) {
463     switch (testNum) {
464     case 0:
465     {
466       PetscInt    numPoints[2]         = {12, 2};
467       PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
468       PetscInt    cones[16]            = {2, 3, 4, 5, 6, 7, 8, 9,  5, 4, 10, 11, 7, 12, 13, 8};
469       PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0,  0, 0,  0,  0, 0,  0,  0, 0};
470       PetscScalar vertexCoords[36]     = {-1.0, -0.5, -0.5,  -1.0,  0.5, -0.5,  0.0,  0.5, -0.5,   0.0, -0.5, -0.5,
471                                           -1.0, -0.5,  0.5,   0.0, -0.5,  0.5,  0.0,  0.5,  0.5,  -1.0,  0.5,  0.5,
472                                           1.0,  0.5, -0.5,   1.0, -0.5, -0.5,  1.0, -0.5,  0.5,   1.0,  0.5,  0.5};
473 
474       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
475     }
476     break;
477     case 1:
478     {
479       PetscInt    numPoints[2]        = {8, 1};
480       PetscInt    coneSize[9]         = {8, 0, 0, 0, 0, 0, 0, 0, 0};
481       PetscInt    cones[8]            = {1, 2, 3, 4, 5, 6, 7, 8};
482       PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
483       PetscScalar vertexCoords[24]    = {-1.0, -1.0, -1.0,  -1.0,  1.0, -1.0,  1.0,  1.0, -1.0,   1.0, -1.0, -1.0,
484                                          -1.0, -1.0,  1.0,   1.0, -1.0,  1.0,  1.0,  1.0,  1.0,  -1.0,  1.0,  1.0};
485 
486       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
487     }
488     break;
489     default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
490     }
491   } else {
492     PetscInt numPoints[4] = {0, 0, 0, 0};
493 
494     ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
495   }
496   ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
497   ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr);
498   ierr = DMDestroy(dm);CHKERRQ(ierr);
499   *dm  = idm;
500   PetscFunctionReturn(0);
501 }
502 
503 PetscErrorCode CreateTensorProductHybrid_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
504 {
505   DM             idm, hdm = NULL;
506   DMLabel        faultLabel;
507   PetscInt       p;
508   PetscMPIInt    rank;
509   PetscErrorCode ierr;
510 
511   PetscFunctionBegin;
512   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
513   if (!rank) {
514     switch (testNum) {
515     case 0:
516     {
517       PetscInt    numPoints[2]         = {12, 2};
518       PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
519       PetscInt    cones[16]            = {2, 3, 4, 5, 6, 7, 8, 9,  5, 4, 10, 11, 7, 12, 13, 8};
520       PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0,  0, 0,  0,  0, 0,  0,  0, 0};
521       PetscScalar vertexCoords[36]     = {-1.0, -0.5, -0.5,  -1.0,  0.5, -0.5,  0.0,  0.5, -0.5,   0.0, -0.5, -0.5,
522                                           -1.0, -0.5,  0.5,   0.0, -0.5,  0.5,  0.0,  0.5,  0.5,  -1.0,  0.5,  0.5,
523                                           1.0,  0.5, -0.5,   1.0, -0.5, -0.5,  1.0, -0.5,  0.5,   1.0,  0.5,  0.5};
524       PetscInt    faultPoints[4]       = {2, 3, 5, 6};
525 
526       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
527       for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
528     }
529     break;
530     case 1:
531     {
532       PetscInt    numPoints[2]         = {30, 7};
533       PetscInt    coneSize[37]         = {8,8,8,8,8,8,8, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
534       PetscInt    cones[56]            = { 8, 21, 20,  7, 13, 12, 23, 24,
535                                           14, 15, 10,  9, 13,  8, 21, 24,
536                                           15, 16, 11, 10, 24, 21, 22, 25,
537                                           30, 29, 28, 21, 35, 24, 33, 34,
538                                           24, 21, 30, 35, 25, 36, 31, 22,
539                                           27, 20, 21, 28, 32, 33, 24, 23,
540                                           15, 24, 13, 14, 19, 18, 17, 26};
541       PetscInt    coneOrientations[56] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
542       PetscScalar vertexCoords[90]     = {-2.0, -2.0, -2.0,  -2.0, -1.0, -2.0,  -3.0,  0.0, -2.0,  -2.0,  1.0, -2.0,  -2.0,  2.0, -2.0,  -2.0, -2.0,  0.0,
543                                           -2.0, -1.0,  0.0,  -3.0,  0.0,  0.0,  -2.0,  1.0,  0.0,  -2.0,  2.0,  0.0,  -2.0, -1.0,  2.0,  -3.0,  0.0,  2.0,
544                                           -2.0,  1.0,  2.0,   0.0, -2.0, -2.0,   0.0,  0.0, -2.0,   0.0,  2.0, -2.0,   0.0, -2.0,  0.0,   0.0,  0.0,  0.0,
545                                            0.0,  2.0,  0.0,   0.0,  0.0,  2.0,   2.0, -2.0, -2.0,   2.0, -1.0, -2.0,   3.0,  0.0, -2.0,   2.0,  1.0, -2.0,
546                                            2.0,  2.0, -2.0,   2.0, -2.0,  0.0,   2.0, -1.0,  0.0,   3.0,  0.0,  0.0,   2.0,  1.0,  0.0,   2.0,  2.0,  0.0};
547       PetscInt    faultPoints[6]       = {20, 21, 22, 23, 24, 25};
548 
549       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
550       for (p = 0; p < 6; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
551     }
552     break;
553     default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
554     }
555     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
556     ierr = PetscObjectSetOptionsPrefix((PetscObject) idm, "in_");CHKERRQ(ierr);
557     ierr = DMSetFromOptions(idm);CHKERRQ(ierr);
558     ierr = DMViewFromOptions(idm, NULL, "-dm_view");CHKERRQ(ierr);
559     ierr = DMGetLabel(*dm, "fault", &faultLabel);CHKERRQ(ierr);
560     ierr = DMPlexCreateHybridMesh(idm, faultLabel, NULL, NULL, NULL, NULL, &hdm);CHKERRQ(ierr);
561     ierr = DMDestroy(&idm);CHKERRQ(ierr);
562     ierr = DMDestroy(dm);CHKERRQ(ierr);
563     *dm  = hdm;
564   } else {
565     PetscInt numPoints[4] = {0, 0, 0, 0};
566 
567     ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
568     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
569     ierr = PetscObjectSetOptionsPrefix((PetscObject) idm, "in_");CHKERRQ(ierr);
570     ierr = DMSetFromOptions(idm);CHKERRQ(ierr);
571     ierr = DMViewFromOptions(idm, NULL, "-dm_view");CHKERRQ(ierr);
572     ierr = DMPlexCreateHybridMesh(idm, NULL, NULL, NULL, NULL, NULL, &hdm);CHKERRQ(ierr);
573     ierr = DMDestroy(&idm);CHKERRQ(ierr);
574     ierr = DMDestroy(dm);CHKERRQ(ierr);
575     *dm  = hdm;
576   }
577   PetscFunctionReturn(0);
578 }
579 
580 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
581 {
582   PetscInt       dim         = user->dim;
583   PetscBool      cellHybrid  = user->cellHybrid;
584   PetscBool      cellSimplex = user->cellSimplex;
585   PetscMPIInt    rank, size;
586   PetscErrorCode ierr;
587 
588   PetscFunctionBegin;
589   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
590   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
591   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
592   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
593   ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
594   switch (dim) {
595   case 1:
596     if (cellHybrid) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %d", dim);
597     ierr = CreateSimplex_1D(comm, dm);CHKERRQ(ierr);
598     break;
599   case 2:
600     if (cellSimplex) {
601       if (cellHybrid) {
602         ierr = CreateSimplexHybrid_2D(comm, user->testNum, dm);CHKERRQ(ierr);
603       } else {
604         ierr = CreateSimplex_2D(comm, dm);CHKERRQ(ierr);
605       }
606     } else {
607       if (cellHybrid) {
608         ierr = CreateTensorProductHybrid_2D(comm, user->testNum, dm);CHKERRQ(ierr);
609       } else {
610         ierr = CreateTensorProduct_2D(comm, user->testNum, dm);CHKERRQ(ierr);
611       }
612     }
613     break;
614   case 3:
615     if (cellSimplex) {
616       if (cellHybrid) {
617         ierr = CreateSimplexHybrid_3D(comm, user->testNum, dm);CHKERRQ(ierr);
618       } else {
619         ierr = CreateSimplex_3D(comm, user->testNum, dm);CHKERRQ(ierr);
620       }
621     } else {
622       if (cellHybrid) {
623         ierr = CreateTensorProductHybrid_3D(comm, user->testNum, dm);CHKERRQ(ierr);
624       } else {
625         ierr = CreateTensorProduct_3D(comm, user->testNum, dm);CHKERRQ(ierr);
626       }
627     }
628     break;
629   default:
630     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
631   }
632   if (user->testPartition && size > 1) {
633     PetscPartitioner part;
634     PetscInt  *sizes  = NULL;
635     PetscInt  *points = NULL;
636 
637     if (!rank) {
638       if (dim == 2 && cellSimplex && !cellHybrid && size == 2) {
639         switch (user->testNum) {
640         case 0: {
641           PetscInt triSizes_p2[2]  = {1, 1};
642           PetscInt triPoints_p2[2] = {0, 1};
643 
644           ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr);
645           ierr = PetscArraycpy(sizes,  triSizes_p2, 2);CHKERRQ(ierr);
646           ierr = PetscArraycpy(points, triPoints_p2, 2);CHKERRQ(ierr);break;}
647         default:
648           SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum);
649         }
650       } else if (dim == 2 && cellSimplex && cellHybrid && size == 2) {
651         switch (user->testNum) {
652         case 0: {
653           PetscInt triSizes_p2[2]  = {1, 2};
654           PetscInt triPoints_p2[3] = {0, 1, 2};
655 
656           ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr);
657           ierr = PetscArraycpy(sizes,  triSizes_p2, 2);CHKERRQ(ierr);
658           ierr = PetscArraycpy(points, triPoints_p2, 3);CHKERRQ(ierr);break;}
659         default:
660           SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular hybrid mesh on 2 procs", user->testNum);
661         }
662       } else if (dim == 2 && !cellSimplex && !cellHybrid && size == 2) {
663         switch (user->testNum) {
664         case 0: {
665           PetscInt quadSizes_p2[2]  = {1, 1};
666           PetscInt quadPoints_p2[2] = {0, 1};
667 
668           ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr);
669           ierr = PetscArraycpy(sizes,  quadSizes_p2, 2);CHKERRQ(ierr);
670           ierr = PetscArraycpy(points, quadPoints_p2, 2);CHKERRQ(ierr);break;}
671         default:
672           SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for quadrilateral mesh on 2 procs", user->testNum);
673         }
674       } else if (dim == 2 && !cellSimplex && cellHybrid && size == 2) {
675         switch (user->testNum) {
676         case 0: {
677           PetscInt quadSizes_p2[2]  = {1, 2};
678           PetscInt quadPoints_p2[3] = {0, 1, 2};
679 
680           ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr);
681           ierr = PetscArraycpy(sizes,  quadSizes_p2, 2);CHKERRQ(ierr);
682           ierr = PetscArraycpy(points, quadPoints_p2, 3);CHKERRQ(ierr);break;}
683         default:
684           SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for quadrilateral hybrid mesh on 2 procs", user->testNum);
685         }
686       } else if (dim == 3 && cellSimplex && !cellHybrid && size == 2) {
687         switch (user->testNum) {
688         case 0: {
689           PetscInt tetSizes_p2[2]  = {1, 1};
690           PetscInt tetPoints_p2[2] = {0, 1};
691 
692           ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr);
693           ierr = PetscArraycpy(sizes,  tetSizes_p2, 2);CHKERRQ(ierr);
694           ierr = PetscArraycpy(points, tetPoints_p2, 2);CHKERRQ(ierr);break;}
695         case 1: {
696           PetscInt tetSizes_p2[2]  = {1, 1};
697           PetscInt tetPoints_p2[2] = {0, 1};
698 
699           ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr);
700           ierr = PetscArraycpy(sizes,  tetSizes_p2, 2);CHKERRQ(ierr);
701           ierr = PetscArraycpy(points, tetPoints_p2, 2);CHKERRQ(ierr);break;}
702         default:
703           SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for tetrahedral mesh on 2 procs", user->testNum);
704         }
705       } else if (dim == 3 && cellSimplex && cellHybrid && size == 2) {
706         switch (user->testNum) {
707         case 0: {
708           PetscInt tetSizes_p2[2]  = {1, 2};
709           PetscInt tetPoints_p2[3] = {0, 1, 2};
710 
711           ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr);
712           ierr = PetscArraycpy(sizes,  tetSizes_p2, 2);CHKERRQ(ierr);
713           ierr = PetscArraycpy(points, tetPoints_p2, 3);CHKERRQ(ierr);break;}
714         case 1: {
715           PetscInt tetSizes_p2[2]  = {3, 4};
716           PetscInt tetPoints_p2[7] = {0, 3, 5, 1, 2, 4, 6};
717 
718           ierr = PetscMalloc2(2, &sizes, 7, &points);CHKERRQ(ierr);
719           ierr = PetscArraycpy(sizes,  tetSizes_p2, 2);CHKERRQ(ierr);
720           ierr = PetscArraycpy(points, tetPoints_p2, 7);CHKERRQ(ierr);break;}
721         default:
722           SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for tetrahedral hybrid mesh on 2 procs", user->testNum);
723         }
724       } else if (dim == 3 && !cellSimplex && !cellHybrid && size == 2) {
725         switch (user->testNum) {
726         case 0: {
727           PetscInt hexSizes_p2[2]  = {1, 1};
728           PetscInt hexPoints_p2[2] = {0, 1};
729 
730           ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr);
731           ierr = PetscArraycpy(sizes,  hexSizes_p2, 2);CHKERRQ(ierr);
732           ierr = PetscArraycpy(points, hexPoints_p2, 2);CHKERRQ(ierr);break;}
733         default:
734           SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for hexahedral mesh on 2 procs", user->testNum);
735         }
736       } else if (dim == 3 && !cellSimplex && cellHybrid && size == 2) {
737         switch (user->testNum) {
738         case 0: {
739           PetscInt hexSizes_p2[2]  = {1, 1};
740           PetscInt hexPoints_p2[2] = {0, 1};
741 
742           ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr);
743           ierr = PetscArraycpy(sizes,  hexSizes_p2, 2);CHKERRQ(ierr);
744           ierr = PetscArraycpy(points, hexPoints_p2, 2);CHKERRQ(ierr);break;}
745         case 1: {
746           PetscInt hexSizes_p2[2]  = {5, 4};
747           PetscInt hexPoints_p2[9] = {3, 4, 5, 7, 8, 0, 1, 2, 6};
748 
749           ierr = PetscMalloc2(2, &sizes, 9, &points);CHKERRQ(ierr);
750           ierr = PetscArraycpy(sizes,  hexSizes_p2, 2);CHKERRQ(ierr);
751           ierr = PetscArraycpy(points, hexPoints_p2, 9);CHKERRQ(ierr);break;}
752         default:
753           SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for hexahedral hybrid mesh on 2 procs", user->testNum);
754         }
755       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition");
756     }
757     ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr);
758     ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);CHKERRQ(ierr);
759     ierr = PetscPartitionerShellSetPartition(part, size, sizes, points);CHKERRQ(ierr);
760     ierr = PetscFree2(sizes, points);CHKERRQ(ierr);
761   } else {
762     PetscPartitioner part;
763 
764     ierr = DMPlexGetPartitioner(*dm,&part);CHKERRQ(ierr);
765     ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
766   }
767   {
768     DM pdm = NULL;
769 
770     ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr);
771     if (pdm) {
772       ierr = DMViewFromOptions(pdm, NULL, "-dm_view");CHKERRQ(ierr);
773       ierr = DMDestroy(dm);CHKERRQ(ierr);
774       *dm  = pdm;
775     }
776   }
777   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
778   if (user->uninterpolate || user->reinterpolate) {
779     DM udm = NULL;
780 
781     ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr);
782     ierr = DMPlexCopyCoordinates(*dm, udm);CHKERRQ(ierr);
783     ierr = DMDestroy(dm);CHKERRQ(ierr);
784     *dm  = udm;
785   }
786   if (user->reinterpolate) {
787     DM idm = NULL;
788 
789     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
790     ierr = DMPlexCopyCoordinates(*dm, idm);CHKERRQ(ierr);
791     ierr = DMDestroy(dm);CHKERRQ(ierr);
792     *dm  = idm;
793   }
794   ierr = PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh");CHKERRQ(ierr);
795   ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
796   ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "hyb_");CHKERRQ(ierr);
797   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
798   PetscFunctionReturn(0);
799 }
800 
801 int main(int argc, char **argv)
802 {
803   DM             dm;
804   AppCtx         user;                 /* user-defined work context */
805   PetscErrorCode ierr;
806 
807   ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
808   ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
809   ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
810   ierr = DMDestroy(&dm);CHKERRQ(ierr);
811   ierr = PetscFinalize();
812   return ierr;
813 }
814 
815 /*TEST
816 
817   # 1D Simplex 29-31
818   testset:
819     args: -dim 1 -cell_hybrid 0 -hyb_dm_plex_check_all -dm_plex_check_all
820     test:
821       suffix: 29
822     test:
823       suffix: 30
824       args: -dm_refine 1
825     test:
826       suffix: 31
827       args: -dm_refine 5
828 
829   # 2D Simplex 0-3
830   testset:
831     args: -dim 2 -cell_hybrid 0 -hyb_dm_plex_check_all -dm_plex_check_all
832     test:
833       suffix: 0
834       args: -hyb_dm_plex_check_faces
835     test:
836       suffix: 1
837       args: -dm_refine 1 -hyb_dm_plex_check_faces
838     test:
839       suffix: 2
840       nsize: 2
841       args: -hyb_dm_plex_check_faces
842     test:
843       suffix: 3
844       nsize: 2
845       args: -dm_refine 1 -hyb_dm_plex_check_faces
846     test:
847       suffix: 32
848       args: -dm_refine 1 -uninterpolate
849     test:
850       suffix: 33
851       nsize: 2
852       args: -dm_refine 1 -uninterpolate
853     test:
854       suffix: 34
855       nsize: 2
856       args: -dm_refine 3 -uninterpolate
857 
858   # 2D Hybrid Simplex 4-7
859   testset:
860     args: -dim 2 -hyb_dm_plex_check_all -in_dm_plex_check_all -dm_plex_check_all
861     test:
862       suffix: 4
863     test:
864       suffix: 5
865       args: -dm_refine 1
866     test:
867       suffix: 6
868       nsize: 2
869     test:
870       suffix: 7
871       nsize: 2
872       args: -dm_refine 1
873     test:
874       suffix: 24
875       args: -test_num 1 -dm_refine 1
876 
877   # 2D Quad 12-13
878   testset:
879     args: -dim 2 -cell_simplex 0 -cell_hybrid 0 -hyb_dm_plex_check_all -dm_plex_check_all
880     test:
881       suffix: 12
882       args: -dm_refine 1
883     test:
884       suffix: 13
885       nsize: 2
886       args: -dm_refine 1
887 
888   # 2D Hybrid Quad 27-28
889   testset:
890     args: -dim 2 -cell_simplex 0 -hyb_dm_plex_check_all -in_dm_plex_check_all -dm_plex_check_all
891     test:
892       suffix: 27
893       args: -dm_refine 1
894     test:
895       suffix: 28
896       nsize: 2
897       args: -dm_refine 1
898 
899   # 3D Simplex 8-11
900   testset:
901     args: -dim 3 -cell_hybrid 0 -hyb_dm_plex_check_all -dm_plex_check_all
902     test:
903       suffix: 8
904       args: -dm_refine 1
905     test:
906       suffix: 9
907       nsize: 2
908       args: -dm_refine 1
909     test:
910       suffix: 10
911       args: -test_num 1 -dm_refine 1
912     test:
913       suffix: 11
914       nsize: 2
915       args: -test_num 1 -dm_refine 1
916     test:
917       suffix: 25
918       args: -test_num 2 -dm_refine 2
919 
920   # 3D Hybrid Simplex 16-19
921   testset:
922     args: -dim 3 -hyb_dm_plex_check_all -in_dm_plex_check_all -dm_plex_check_all
923     test:
924       suffix: 16
925       args: -dm_refine 1
926     test:
927       suffix: 17
928       nsize: 2
929       args: -dm_refine 1
930     test:
931       suffix: 18
932       args: -test_num 1 -dm_refine 1
933     test:
934       suffix: 19
935       nsize: 2
936       args: -test_num 1 -dm_refine 1
937 
938   # 3D Hex 14-15
939   testset:
940     args: -dim 3 -cell_simplex 0 -cell_hybrid 0 -hyb_dm_plex_check_all -dm_plex_check_all
941     test:
942       suffix: 14
943       args: -dm_refine 1
944     test:
945       suffix: 15
946       nsize: 2
947      args: -dm_refine 1
948     test:
949       suffix: 26
950       args: -test_num 1 -dm_refine 2
951 
952   # 3D Hybrid Hex 20-23
953   testset:
954     args: -dim 3 -cell_simplex 0 -hyb_dm_plex_check_all -in_dm_plex_check_all -dm_plex_check_all
955     test:
956       suffix: 20
957       args: -dm_refine 1
958     test:
959       suffix: 21
960       nsize: 2
961       args: -dm_refine 1
962     test:
963       suffix: 22
964       args: -test_num 1 -dm_refine 1
965     test:
966       suffix: 23
967       nsize: 2
968       args: -test_num 1 -dm_refine 1
969 
970   # Hybrid interpolation
971   #   TODO Setup new tests (like -reinterpolate) that interpolate hybrid cells
972   testset:
973     nsize: 2
974     args: -test_partition 0 -petscpartitioner_type simple -dm_view -hyb_dm_plex_check_all -in_dm_plex_check_all -dm_plex_check_all
975     test:
976       suffix: hybint_2d_0
977       args: -dim 2 -dm_refine 2
978     test:
979       suffix: hybint_2d_1
980       args: -dim 2 -dm_refine 2 -test_num 1
981     test:
982       suffix: hybint_3d_0
983       args: -dim 3 -dm_refine 1
984     test:
985       suffix: hybint_3d_1
986       args: -dim 3 -dm_refine 1 -test_num 1
987 
988 TEST*/
989