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