xref: /petsc/src/dm/impls/plex/tests/ex4.c (revision 030f984af8d8bb4c203755d35bded3c05b3d83ce)
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();CHKERRQ(ierr);
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, -1, 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, -1,  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,  -2,  0,  0,  0,   0,  0,  0,   0,  0, -1,  -1,  0, -1,  -1, -1, -1,   0,  0, -1,  -1,  0, -1,  -1, -1, -1,  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 = DMViewFromOptions(*dm, NULL, "-dm_view_pre");CHKERRQ(ierr);
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