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