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