xref: /petsc/src/dm/impls/plex/tests/ex5.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94)
1*c4762a1bSJed Brown static char help[] = "Tests for creation of hybrid meshes\n\n";
2*c4762a1bSJed Brown 
3*c4762a1bSJed Brown /* TODO
4*c4762a1bSJed Brown  - Propagate hybridSize with distribution
5*c4762a1bSJed Brown  - Test with multiple fault segments
6*c4762a1bSJed Brown  - Test with embedded fault
7*c4762a1bSJed Brown  - Test with multiple faults
8*c4762a1bSJed Brown  - Move over all PyLith tests?
9*c4762a1bSJed Brown */
10*c4762a1bSJed Brown 
11*c4762a1bSJed Brown #include <petscdmplex.h>
12*c4762a1bSJed Brown /* List of test meshes
13*c4762a1bSJed Brown 
14*c4762a1bSJed Brown Triangle
15*c4762a1bSJed Brown --------
16*c4762a1bSJed Brown Test 0:
17*c4762a1bSJed Brown Two triangles sharing a face
18*c4762a1bSJed Brown 
19*c4762a1bSJed Brown         4
20*c4762a1bSJed Brown       / | \
21*c4762a1bSJed Brown      8  |  9
22*c4762a1bSJed Brown     /   |   \
23*c4762a1bSJed Brown    2  0 7 1  5
24*c4762a1bSJed Brown     \   |   /
25*c4762a1bSJed Brown      6  |  10
26*c4762a1bSJed Brown       \ | /
27*c4762a1bSJed Brown         3
28*c4762a1bSJed Brown 
29*c4762a1bSJed Brown should become two triangles separated by a zero-volume cell with 4 vertices
30*c4762a1bSJed Brown 
31*c4762a1bSJed Brown         5--16--8              4--12--6 3
32*c4762a1bSJed Brown       / |      | \          / |      | | \
33*c4762a1bSJed Brown     11  |      |  12       9  |      | |  4
34*c4762a1bSJed Brown     /   |      |   \      /   |      | |   \
35*c4762a1bSJed Brown    3  0 10  2 14 1  6    2  0 8  1  10 6 0  1
36*c4762a1bSJed Brown     \   |      |   /      \   |      | |   /
37*c4762a1bSJed Brown      9  |      |  13       7  |      | |  5
38*c4762a1bSJed Brown       \ |      | /          \ |      | | /
39*c4762a1bSJed Brown         4--15--7              3--11--5 2
40*c4762a1bSJed Brown 
41*c4762a1bSJed Brown Test 1:
42*c4762a1bSJed Brown Four triangles sharing two faces which are oriented against each other
43*c4762a1bSJed Brown 
44*c4762a1bSJed Brown           9
45*c4762a1bSJed Brown          / \
46*c4762a1bSJed Brown         /   \
47*c4762a1bSJed Brown       17  2  16
48*c4762a1bSJed Brown       /       \
49*c4762a1bSJed Brown      /         \
50*c4762a1bSJed Brown     8-----15----5
51*c4762a1bSJed Brown      \         /|\
52*c4762a1bSJed Brown       \       / | \
53*c4762a1bSJed Brown       18  3  12 |  14
54*c4762a1bSJed Brown         \   /   |   \
55*c4762a1bSJed Brown          \ /    |    \
56*c4762a1bSJed Brown           4  0 11  1  7
57*c4762a1bSJed Brown            \    |    /
58*c4762a1bSJed Brown             \   |   /
59*c4762a1bSJed Brown             10  |  13
60*c4762a1bSJed Brown               \ | /
61*c4762a1bSJed Brown                \|/
62*c4762a1bSJed Brown                 6
63*c4762a1bSJed Brown 
64*c4762a1bSJed Brown Fault mesh
65*c4762a1bSJed Brown 
66*c4762a1bSJed Brown 0 --> 0
67*c4762a1bSJed Brown 1 --> 1
68*c4762a1bSJed Brown 2 --> 2
69*c4762a1bSJed Brown 3 --> 3
70*c4762a1bSJed Brown 4 --> 5
71*c4762a1bSJed Brown 5 --> 6
72*c4762a1bSJed Brown 6 --> 8
73*c4762a1bSJed Brown 7 --> 11
74*c4762a1bSJed Brown 8 --> 15
75*c4762a1bSJed Brown 
76*c4762a1bSJed Brown        2
77*c4762a1bSJed Brown        |
78*c4762a1bSJed Brown   6----8----4
79*c4762a1bSJed Brown        |    |
80*c4762a1bSJed Brown        3    |
81*c4762a1bSJed Brown           0-7-1
82*c4762a1bSJed Brown             |
83*c4762a1bSJed Brown             |
84*c4762a1bSJed Brown             5
85*c4762a1bSJed Brown 
86*c4762a1bSJed Brown should become four triangles separated by two zero-volume cells with 4 vertices
87*c4762a1bSJed Brown 
88*c4762a1bSJed Brown           11
89*c4762a1bSJed Brown           / \
90*c4762a1bSJed Brown          /   \
91*c4762a1bSJed Brown         /     \
92*c4762a1bSJed Brown       22   2   21
93*c4762a1bSJed Brown       /         \
94*c4762a1bSJed Brown      /           \
95*c4762a1bSJed Brown    10-----20------7
96*c4762a1bSJed Brown 28  |     5    26/ \
97*c4762a1bSJed Brown    14----25----12   \
98*c4762a1bSJed Brown      \         /|   |\
99*c4762a1bSJed Brown       \       / |   | \
100*c4762a1bSJed Brown       23  3  17 |   |  19
101*c4762a1bSJed Brown         \   /   |   |   \
102*c4762a1bSJed Brown          \ /    |   |    \
103*c4762a1bSJed Brown           6  0 24 4 16 1  9
104*c4762a1bSJed Brown            \    |   |    /
105*c4762a1bSJed Brown             \   |   |   /
106*c4762a1bSJed Brown             15  |   |  18
107*c4762a1bSJed Brown               \ |   | /
108*c4762a1bSJed Brown                \|   |/
109*c4762a1bSJed Brown                13---8
110*c4762a1bSJed Brown                  27
111*c4762a1bSJed Brown 
112*c4762a1bSJed Brown Tetrahedron
113*c4762a1bSJed Brown -----------
114*c4762a1bSJed Brown Test 0:
115*c4762a1bSJed Brown Two tets sharing a face
116*c4762a1bSJed Brown 
117*c4762a1bSJed Brown  cell   5 _______    cell
118*c4762a1bSJed Brown  0    / | \      \       1
119*c4762a1bSJed Brown     16  |  18     22
120*c4762a1bSJed Brown     /8 19 10\      \
121*c4762a1bSJed Brown    2-15-|----4--21--6
122*c4762a1bSJed Brown     \  9| 7 /      /
123*c4762a1bSJed Brown     14  |  17     20
124*c4762a1bSJed Brown       \ | /      /
125*c4762a1bSJed Brown         3-------
126*c4762a1bSJed Brown 
127*c4762a1bSJed Brown should become two tetrahedrons separated by a zero-volume cell with 3 faces/3 edges/6 vertices
128*c4762a1bSJed Brown 
129*c4762a1bSJed Brown  cell   6 ___36___10______    cell
130*c4762a1bSJed Brown  0    / | \        |\      \     1
131*c4762a1bSJed Brown     24  |  26      | 32     30
132*c4762a1bSJed Brown     /12 27 14\    33  \      \
133*c4762a1bSJed Brown    3-23-|----5--35-|---9--29--7
134*c4762a1bSJed Brown     \ 13| 11/      |18 /      /
135*c4762a1bSJed Brown     22  |  25      | 31     28
136*c4762a1bSJed Brown       \ | /        |/      /
137*c4762a1bSJed Brown         4----34----8------
138*c4762a1bSJed Brown          cell 2
139*c4762a1bSJed Brown 
140*c4762a1bSJed Brown In parallel,
141*c4762a1bSJed Brown 
142*c4762a1bSJed Brown  cell   5 ___28____8      4______    cell
143*c4762a1bSJed Brown  0    / | \        |\     |\      \     0
144*c4762a1bSJed Brown     19  |   21     | 24   | 13  6  11
145*c4762a1bSJed Brown     /10 22 12\    25  \   |8 \      \
146*c4762a1bSJed Brown    2-18-|----4--27-|---7  14  3--10--1
147*c4762a1bSJed Brown     \ 11| 9 /      |13 /  |  /      /
148*c4762a1bSJed Brown     17  |  20      | 23   | 12  5  9
149*c4762a1bSJed Brown       \ | /        |/     |/      /
150*c4762a1bSJed Brown         3----26----6      2------
151*c4762a1bSJed Brown          cell 1
152*c4762a1bSJed Brown 
153*c4762a1bSJed Brown Test 1:
154*c4762a1bSJed Brown Four tets sharing two faces
155*c4762a1bSJed Brown 
156*c4762a1bSJed Brown Cells:    0-3,4-5
157*c4762a1bSJed Brown Vertices: 6-15
158*c4762a1bSJed Brown Faces:    16-29,30-34
159*c4762a1bSJed Brown Edges:    35-52,53-56
160*c4762a1bSJed Brown 
161*c4762a1bSJed Brown Quadrilateral
162*c4762a1bSJed Brown -------------
163*c4762a1bSJed Brown Test 0:
164*c4762a1bSJed Brown Two quads sharing a face
165*c4762a1bSJed Brown 
166*c4762a1bSJed Brown    5--10---4--14---7
167*c4762a1bSJed Brown    |       |       |
168*c4762a1bSJed Brown   11   0   9   1  13
169*c4762a1bSJed Brown    |       |       |
170*c4762a1bSJed Brown    2---8---3--12---6
171*c4762a1bSJed Brown 
172*c4762a1bSJed Brown should become two quads separated by a zero-volume cell with 4 vertices
173*c4762a1bSJed Brown 
174*c4762a1bSJed Brown    6--13---5-20-10--17---8    5--10---4-14--7  4---7---2
175*c4762a1bSJed Brown    |       |     |       |    |       |     |  |       |
176*c4762a1bSJed Brown   14   0  12  2 18   1  16   11   0   9  1 12  8   0   6
177*c4762a1bSJed Brown    |       |     |       |    |       |     |  |       |
178*c4762a1bSJed Brown    3--11---4-19--9--15---7    2---8---3-13--6  3---5---1
179*c4762a1bSJed Brown 
180*c4762a1bSJed Brown Test 1:
181*c4762a1bSJed Brown 
182*c4762a1bSJed Brown Original mesh with 9 cells,
183*c4762a1bSJed Brown 
184*c4762a1bSJed Brown   9 ----10 ----11 ----12
185*c4762a1bSJed Brown   |      |      |      |
186*c4762a1bSJed Brown   |      |      |      |
187*c4762a1bSJed Brown   |      |      |      |
188*c4762a1bSJed Brown   |      |      |      |
189*c4762a1bSJed Brown  13 ----14 ----15 ----16
190*c4762a1bSJed Brown   |      |      |      |
191*c4762a1bSJed Brown   |      |      |      |
192*c4762a1bSJed Brown   |      |      |      |
193*c4762a1bSJed Brown   |      |      |      |
194*c4762a1bSJed Brown  17 ----18 ----19 ----20
195*c4762a1bSJed Brown   |      |      |      |
196*c4762a1bSJed Brown   |      |      |      |
197*c4762a1bSJed Brown   |      |      |      |
198*c4762a1bSJed Brown   |      |      |      |
199*c4762a1bSJed Brown  21 ----22 ----23 ----24
200*c4762a1bSJed Brown 
201*c4762a1bSJed Brown After first fault,
202*c4762a1bSJed Brown 
203*c4762a1bSJed Brown  12 ----13 ----14-28 ----15
204*c4762a1bSJed Brown   |      |      |  |      |
205*c4762a1bSJed Brown   |  0   |  1   | 9|  2   |
206*c4762a1bSJed Brown   |      |      |  |      |
207*c4762a1bSJed Brown   |      |      |  |      |
208*c4762a1bSJed Brown  16 ----17 ----18-29 ----19
209*c4762a1bSJed Brown   |      |      |  |      |
210*c4762a1bSJed Brown   |  3   |  4   |10|  5   |
211*c4762a1bSJed Brown   |      |      |  |      |
212*c4762a1bSJed Brown   |      |      |  |      |
213*c4762a1bSJed Brown  20 ----21-----22-30 ----23
214*c4762a1bSJed Brown   |      |      |  \--11- |
215*c4762a1bSJed Brown   |  6   |  7   |     8   |
216*c4762a1bSJed Brown   |      |      |         |
217*c4762a1bSJed Brown   |      |      |         |
218*c4762a1bSJed Brown  24 ----25 ----26--------27
219*c4762a1bSJed Brown 
220*c4762a1bSJed Brown After second fault,
221*c4762a1bSJed Brown 
222*c4762a1bSJed Brown  14 ----15 ----16-30 ----17
223*c4762a1bSJed Brown   |      |      |  |      |
224*c4762a1bSJed Brown   |  0   |  1   | 9|  2   |
225*c4762a1bSJed Brown   |      |      |  |      |
226*c4762a1bSJed Brown   |      |      |  |      |
227*c4762a1bSJed Brown  18 ----19 ----20-31 ----21
228*c4762a1bSJed Brown   |      |      |  |      |
229*c4762a1bSJed Brown   |  3   |  4   |10|  5   |
230*c4762a1bSJed Brown   |      |      |  |      |
231*c4762a1bSJed Brown   |      |      |  |      |
232*c4762a1bSJed Brown  33 ----34-----24-32 ----25
233*c4762a1bSJed Brown   |  12  | 13 / |  \-11-- |
234*c4762a1bSJed Brown  22 ----23---/  |         |
235*c4762a1bSJed Brown   |      |   7  |     8   |
236*c4762a1bSJed Brown   |  6   |      |         |
237*c4762a1bSJed Brown   |      |      |         |
238*c4762a1bSJed Brown   |      |      |         |
239*c4762a1bSJed Brown  26 ----27 ----28--------29
240*c4762a1bSJed Brown 
241*c4762a1bSJed Brown Hexahedron
242*c4762a1bSJed Brown ----------
243*c4762a1bSJed Brown Test 0:
244*c4762a1bSJed Brown Two hexes sharing a face
245*c4762a1bSJed Brown 
246*c4762a1bSJed Brown cell   9-----31------8-----42------13 cell
247*c4762a1bSJed Brown 0     /|            /|            /|     1
248*c4762a1bSJed Brown     32 |   15      30|   21      41|
249*c4762a1bSJed Brown     /  |          /  |          /  |
250*c4762a1bSJed Brown    6-----29------7-----40------12  |
251*c4762a1bSJed Brown    |   |     18  |   |     24  |   |
252*c4762a1bSJed Brown    |  36         |  35         |   44
253*c4762a1bSJed Brown    |19 |         |17 |         |23 |
254*c4762a1bSJed Brown   33   |  16    34   |   22   43   |
255*c4762a1bSJed Brown    |   5-----27--|---4-----39--|---11
256*c4762a1bSJed Brown    |  /          |  /          |  /
257*c4762a1bSJed Brown    | 28   14     | 26    20    | 38
258*c4762a1bSJed Brown    |/            |/            |/
259*c4762a1bSJed Brown    2-----25------3-----37------10
260*c4762a1bSJed Brown 
261*c4762a1bSJed Brown should become two hexes separated by a zero-volume cell with 8 vertices
262*c4762a1bSJed Brown 
263*c4762a1bSJed Brown                          cell 2
264*c4762a1bSJed Brown cell  10-----41------9-----62------18----52------14 cell
265*c4762a1bSJed Brown 0     /|            /|            /|            /|     1
266*c4762a1bSJed Brown     42 |   20      40|  32       56|   26      51|
267*c4762a1bSJed Brown     /  |          /  |          /  |          /  |
268*c4762a1bSJed Brown    7-----39------8-----61------17--|-50------13  |
269*c4762a1bSJed Brown    |   |     23  |   |         |   |     29  |   |
270*c4762a1bSJed Brown    |  46         |  45         |   58        |   54
271*c4762a1bSJed Brown    |24 |         |22 |         |30 |         |28 |
272*c4762a1bSJed Brown   43   |  21    44   |        57   |   27   53   |
273*c4762a1bSJed Brown    |   6-----37--|---5-----60--|---16----49--|---12
274*c4762a1bSJed Brown    |  /          |  /          |  /          |  /
275*c4762a1bSJed Brown    | 38   19     | 36   31     | 55    25    | 48
276*c4762a1bSJed Brown    |/            |/            |/            |/
277*c4762a1bSJed Brown    3-----35------4-----59------15----47------11
278*c4762a1bSJed Brown 
279*c4762a1bSJed Brown In parallel,
280*c4762a1bSJed Brown 
281*c4762a1bSJed Brown                          cell 2
282*c4762a1bSJed Brown cell   9-----31------8-----44------13     8----20------4  cell
283*c4762a1bSJed Brown 0     /|            /|            /|     /|           /|     1
284*c4762a1bSJed Brown     32 |   15      30|  22       38|   24 |  10      19|
285*c4762a1bSJed Brown     /  |          /  |          /  |   /  |         /  |
286*c4762a1bSJed Brown    6-----29------7-----43------12  |  7----18------3   |
287*c4762a1bSJed Brown    |   |     18  |   |         |   |  |   |    13  |   |
288*c4762a1bSJed Brown    |  36         |  35         |   40 |  26        |   22
289*c4762a1bSJed Brown    |19 |         |17 |         |20 |  |14 |        |12 |
290*c4762a1bSJed Brown   33   |  16    34   |        39   |  25  |  11   21   |
291*c4762a1bSJed Brown    |   5-----27--|---4-----42--|---11 |   6----17--|---2
292*c4762a1bSJed Brown    |  /          |  /          |  /   |  /         |  /
293*c4762a1bSJed Brown    | 28   14     | 26   21     | 37   |23     9    | 16
294*c4762a1bSJed Brown    |/            |/            |/     |/           |/
295*c4762a1bSJed Brown    2-----25------3-----41------10     5----15------1
296*c4762a1bSJed Brown 
297*c4762a1bSJed Brown Test 1:
298*c4762a1bSJed Brown 
299*c4762a1bSJed Brown */
300*c4762a1bSJed Brown 
301*c4762a1bSJed Brown typedef struct {
302*c4762a1bSJed Brown   PetscInt  debug;         /* The debugging level */
303*c4762a1bSJed Brown   PetscInt  dim;           /* The topological mesh dimension */
304*c4762a1bSJed Brown   PetscBool cellSimplex;   /* Use simplices or hexes */
305*c4762a1bSJed Brown   PetscBool testPartition; /* Use a fixed partitioning for testing */
306*c4762a1bSJed Brown   PetscInt  testNum;       /* The particular mesh to test */
307*c4762a1bSJed Brown } AppCtx;
308*c4762a1bSJed Brown 
309*c4762a1bSJed Brown PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
310*c4762a1bSJed Brown {
311*c4762a1bSJed Brown   PetscErrorCode ierr;
312*c4762a1bSJed Brown 
313*c4762a1bSJed Brown   PetscFunctionBegin;
314*c4762a1bSJed Brown   options->debug         = 0;
315*c4762a1bSJed Brown   options->dim           = 2;
316*c4762a1bSJed Brown   options->cellSimplex   = PETSC_TRUE;
317*c4762a1bSJed Brown   options->testPartition = PETSC_TRUE;
318*c4762a1bSJed Brown   options->testNum       = 0;
319*c4762a1bSJed Brown 
320*c4762a1bSJed Brown   ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr);
321*c4762a1bSJed Brown   ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr);
322*c4762a1bSJed Brown   ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr);
323*c4762a1bSJed Brown   ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr);
324*c4762a1bSJed Brown   ierr = PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL);CHKERRQ(ierr);
325*c4762a1bSJed Brown   ierr = PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr);
326*c4762a1bSJed Brown   ierr = PetscOptionsEnd();
327*c4762a1bSJed Brown   PetscFunctionReturn(0);
328*c4762a1bSJed Brown }
329*c4762a1bSJed Brown 
330*c4762a1bSJed Brown PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
331*c4762a1bSJed Brown {
332*c4762a1bSJed Brown   DM             idm;
333*c4762a1bSJed Brown   PetscInt       p;
334*c4762a1bSJed Brown   PetscMPIInt    rank;
335*c4762a1bSJed Brown   PetscErrorCode ierr;
336*c4762a1bSJed Brown 
337*c4762a1bSJed Brown   PetscFunctionBegin;
338*c4762a1bSJed Brown   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
339*c4762a1bSJed Brown   if (!rank) {
340*c4762a1bSJed Brown     switch (testNum) {
341*c4762a1bSJed Brown     case 0:
342*c4762a1bSJed Brown     {
343*c4762a1bSJed Brown       PetscInt    numPoints[2]        = {4, 2};
344*c4762a1bSJed Brown       PetscInt    coneSize[6]         = {3, 3, 0, 0, 0, 0};
345*c4762a1bSJed Brown       PetscInt    cones[6]            = {2, 3, 4,  5, 4, 3};
346*c4762a1bSJed Brown       PetscInt    coneOrientations[6] = {0, 0, 0,  0, 0, 0};
347*c4762a1bSJed Brown       PetscScalar vertexCoords[8]     = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
348*c4762a1bSJed Brown       PetscInt    markerPoints[8]     = {2, 1, 3, 1, 4, 1, 5, 1};
349*c4762a1bSJed Brown       PetscInt    faultPoints[2]      = {3, 4};
350*c4762a1bSJed Brown 
351*c4762a1bSJed Brown       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
352*c4762a1bSJed Brown       for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
353*c4762a1bSJed Brown       for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
354*c4762a1bSJed Brown     }
355*c4762a1bSJed Brown     break;
356*c4762a1bSJed Brown     case 1:
357*c4762a1bSJed Brown     {
358*c4762a1bSJed Brown       PetscInt    numPoints[2]         = {6, 4};
359*c4762a1bSJed Brown       PetscInt    coneSize[10]         = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0};
360*c4762a1bSJed Brown       PetscInt    cones[12]            = {4, 6, 5,  5, 6, 7,  8, 5, 9,  8, 4, 5};
361*c4762a1bSJed Brown       PetscInt    coneOrientations[12] = {0, 0, 0,  0, 0, 0,  0, 0, 0,  0, 0, 0};
362*c4762a1bSJed Brown       PetscScalar vertexCoords[12]     = {-1.0, 0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 0.0, -2.0, 1.0, -1.0, 2.0};
363*c4762a1bSJed Brown       PetscInt    markerPoints[6]      = {4, 1, 6, 1, 8, 1};
364*c4762a1bSJed Brown       PetscInt    faultPoints[3]       = {5, 6, 8};
365*c4762a1bSJed Brown 
366*c4762a1bSJed Brown       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
367*c4762a1bSJed Brown       for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
368*c4762a1bSJed Brown       for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
369*c4762a1bSJed Brown     }
370*c4762a1bSJed Brown     break;
371*c4762a1bSJed Brown     default:
372*c4762a1bSJed Brown       SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
373*c4762a1bSJed Brown     }
374*c4762a1bSJed Brown   } else {
375*c4762a1bSJed Brown     PetscInt numPoints[3] = {0, 0, 0};
376*c4762a1bSJed Brown 
377*c4762a1bSJed Brown     ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
378*c4762a1bSJed Brown     ierr = DMCreateLabel(*dm, "fault");CHKERRQ(ierr);
379*c4762a1bSJed Brown   }
380*c4762a1bSJed Brown   ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
381*c4762a1bSJed Brown   ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr);
382*c4762a1bSJed Brown   ierr = DMDestroy(dm);CHKERRQ(ierr);
383*c4762a1bSJed Brown   *dm  = idm;
384*c4762a1bSJed Brown   PetscFunctionReturn(0);
385*c4762a1bSJed Brown }
386*c4762a1bSJed Brown 
387*c4762a1bSJed Brown PetscErrorCode CreateSimplex_3D(MPI_Comm comm, AppCtx *user, DM dm)
388*c4762a1bSJed Brown {
389*c4762a1bSJed Brown   PetscInt       depth = 3, testNum  = user->testNum, p;
390*c4762a1bSJed Brown   PetscMPIInt    rank;
391*c4762a1bSJed Brown   PetscErrorCode ierr;
392*c4762a1bSJed Brown 
393*c4762a1bSJed Brown   PetscFunctionBegin;
394*c4762a1bSJed Brown   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
395*c4762a1bSJed Brown   if (!rank) {
396*c4762a1bSJed Brown     switch (testNum) {
397*c4762a1bSJed Brown     case 0:
398*c4762a1bSJed Brown     {
399*c4762a1bSJed Brown       PetscInt    numPoints[4]         = {5, 7, 9, 2};
400*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};
401*c4762a1bSJed Brown       PetscInt    cones[47]            = {7, 8, 9, 10,  11, 10, 13, 12,  15, 17, 14,  16, 18, 15,  14, 19, 16,  17, 18, 19,  17, 21, 20,  18, 22, 21,  22, 19, 20,   2, 3,  2, 4,  2, 5,  3, 4,  4, 5,  5, 3,  3, 6,  4, 6,  5, 6};
402*c4762a1bSJed Brown       PetscInt    coneOrientations[47] = {0, 0, 0,  0,   0, -3,  2,  2,   0, -2, -2,   0, -2, -2,   0, -2, -2,   0,  0,  0,   0,  0, -2,   0,  0, -2,  -2,  0,  0,   0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0};
403*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};
404*c4762a1bSJed Brown       PetscInt    markerPoints[20]     = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1};
405*c4762a1bSJed Brown       PetscInt    faultPoints[3]      = {3, 4, 5};
406*c4762a1bSJed Brown 
407*c4762a1bSJed Brown       ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
408*c4762a1bSJed Brown       for (p = 0; p < 10; ++p) {
409*c4762a1bSJed Brown         ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);
410*c4762a1bSJed Brown       }
411*c4762a1bSJed Brown       for(p = 0; p < 3; ++p) {
412*c4762a1bSJed Brown         ierr = DMSetLabelValue(dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);
413*c4762a1bSJed Brown       }
414*c4762a1bSJed Brown     }
415*c4762a1bSJed Brown     break;
416*c4762a1bSJed Brown     case 1:
417*c4762a1bSJed Brown     {
418*c4762a1bSJed Brown       PetscInt    numPoints[4]         = {6, 13, 12, 4};
419*c4762a1bSJed Brown       PetscInt    coneSize[35]         = {4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2};
420*c4762a1bSJed Brown       PetscInt    cones[78]            = {10, 11, 12, 13,  10, 15, 16, 14,  17, 18, 14, 19,  20, 13, 19, 21,  22, 23, 24,  25, 26, 22,  24, 27, 25,  23, 26, 27,  28, 29, 23,  24, 30, 28,  22, 29, 30,   31, 32, 28,  29, 33, 31,  32, 33, 23,  26, 34, 33,  34, 27, 32,  6, 5,  5, 7,  7, 6,  6, 4,  4, 5,  7, 4,  7, 9,  9, 5,  6, 9,  9, 8,  8, 7,  5, 8,  4, 8};
421*c4762a1bSJed Brown       PetscInt    coneOrientations[78] = { 0,  0,  0,  0,  -3,  1,  0,  2,   0,  0, -1,  0,   0, -1, -2,  0,   0,  0,  0,   0,  0, -2,  -2,  0, -2,  -2, -2, -2,   0,  0,  0,   0,  0, -2,   0, -2, -2,    0,  0,  0,   0,  0, -2,  -2, -2,  0,  -2,  0, -2,  -2, -2, -2,  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};
422*c4762a1bSJed Brown       PetscScalar vertexCoords[18]     = {-1.0, 0.0, 0.0,  0.0, -1.0, 0.0,  0.0, 0.0, 1.0,  0.0, 1.0, 0.0,  0.0, 0.0, -1.0,  1.0, 0.0, 0.0};
423*c4762a1bSJed Brown       PetscInt    markerPoints[14]     = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1};
424*c4762a1bSJed Brown       PetscInt    faultPoints[4]       = {5, 6, 7, 8};
425*c4762a1bSJed Brown 
426*c4762a1bSJed Brown       ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
427*c4762a1bSJed Brown       for (p = 0; p < 7; ++p) {
428*c4762a1bSJed Brown         ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);
429*c4762a1bSJed Brown       }
430*c4762a1bSJed Brown       for(p = 0; p < 4; ++p) {
431*c4762a1bSJed Brown         ierr = DMSetLabelValue(dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);
432*c4762a1bSJed Brown       }
433*c4762a1bSJed Brown     }
434*c4762a1bSJed Brown     break;
435*c4762a1bSJed Brown     default:
436*c4762a1bSJed Brown       SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
437*c4762a1bSJed Brown     }
438*c4762a1bSJed Brown   } else {
439*c4762a1bSJed Brown     PetscInt numPoints[4] = {0, 0, 0, 0};
440*c4762a1bSJed Brown 
441*c4762a1bSJed Brown     ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
442*c4762a1bSJed Brown     ierr = DMCreateLabel(dm, "fault");CHKERRQ(ierr);
443*c4762a1bSJed Brown   }
444*c4762a1bSJed Brown   PetscFunctionReturn(0);
445*c4762a1bSJed Brown }
446*c4762a1bSJed Brown 
447*c4762a1bSJed Brown PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
448*c4762a1bSJed Brown {
449*c4762a1bSJed Brown   DM             idm;
450*c4762a1bSJed Brown   PetscInt       p;
451*c4762a1bSJed Brown   PetscMPIInt    rank;
452*c4762a1bSJed Brown   PetscErrorCode ierr;
453*c4762a1bSJed Brown 
454*c4762a1bSJed Brown   PetscFunctionBegin;
455*c4762a1bSJed Brown   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
456*c4762a1bSJed Brown   if (!rank) {
457*c4762a1bSJed Brown     switch (testNum) {
458*c4762a1bSJed Brown     case 0:
459*c4762a1bSJed Brown     case 2:
460*c4762a1bSJed Brown     {
461*c4762a1bSJed Brown       PetscInt    numPoints[2]        = {6, 2};
462*c4762a1bSJed Brown       PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
463*c4762a1bSJed Brown       PetscInt    cones[8]            = {2, 3, 4, 5,  3, 6, 7, 4};
464*c4762a1bSJed Brown       PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
465*c4762a1bSJed Brown       PetscScalar vertexCoords[12]    = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0, -0.5, 1.0, 0.5, 0.0, 0.5, 1.0};
466*c4762a1bSJed Brown       PetscInt    markerPoints[12]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1};
467*c4762a1bSJed Brown       PetscInt    faultPoints[2]      = {3, 4};
468*c4762a1bSJed Brown 
469*c4762a1bSJed Brown       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
470*c4762a1bSJed Brown       for (p = 0; p < 6; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
471*c4762a1bSJed Brown       if (testNum == 0) for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
472*c4762a1bSJed Brown       if (testNum == 2) for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "pfault", faultPoints[p], 1);CHKERRQ(ierr);}
473*c4762a1bSJed Brown     }
474*c4762a1bSJed Brown     break;
475*c4762a1bSJed Brown     case 1:
476*c4762a1bSJed Brown     {
477*c4762a1bSJed Brown       PetscInt    numPoints[2]         = {16, 9};
478*c4762a1bSJed Brown       PetscInt    coneSize[25]         = {4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
479*c4762a1bSJed Brown       PetscInt    cones[36]            = {9,  13, 14, 10,
480*c4762a1bSJed Brown                                           10, 14, 15, 11,
481*c4762a1bSJed Brown                                           11, 15, 16, 12,
482*c4762a1bSJed Brown                                           13, 17, 18, 14,
483*c4762a1bSJed Brown                                           14, 18, 19, 15,
484*c4762a1bSJed Brown                                           15, 19, 20, 16,
485*c4762a1bSJed Brown                                           17, 21, 22, 18,
486*c4762a1bSJed Brown                                           18, 22, 23, 19,
487*c4762a1bSJed Brown                                           19, 23, 24, 20};
488*c4762a1bSJed Brown       PetscInt    coneOrientations[36] = {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};
489*c4762a1bSJed Brown       PetscScalar vertexCoords[32]     = {-3.0,  3.0,  -1.0,  3.0,  1.0,  3.0,  3.0,  3.0,  -3.0,  1.0,  -1.0,  1.0,  1.0,  1.0,  3.0,  1.0,
490*c4762a1bSJed Brown                                           -3.0, -1.0,  -1.0, -1.0,  1.0, -1.0,  3.0, -1.0,  -3.0, -3.0,  -1.0, -3.0,  1.0, -3.0,  3.0, -3.0};
491*c4762a1bSJed Brown       PetscInt    faultPoints[3]       = {11, 15, 19};
492*c4762a1bSJed Brown       PetscInt    fault2Points[2]      = {17, 18};
493*c4762a1bSJed Brown 
494*c4762a1bSJed Brown       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
495*c4762a1bSJed Brown       for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault",  faultPoints[p], 1);CHKERRQ(ierr);}
496*c4762a1bSJed Brown       for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault2", fault2Points[p], 1);CHKERRQ(ierr);}
497*c4762a1bSJed Brown     }
498*c4762a1bSJed Brown     break;
499*c4762a1bSJed Brown     default:
500*c4762a1bSJed Brown       SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
501*c4762a1bSJed Brown     }
502*c4762a1bSJed Brown   } else {
503*c4762a1bSJed Brown     PetscInt numPoints[3] = {0, 0, 0};
504*c4762a1bSJed Brown 
505*c4762a1bSJed Brown     ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
506*c4762a1bSJed Brown     if (testNum == 2) {ierr = DMCreateLabel(*dm, "pfault");CHKERRQ(ierr);}
507*c4762a1bSJed Brown     else              {ierr = DMCreateLabel(*dm, "fault");CHKERRQ(ierr);}
508*c4762a1bSJed Brown   }
509*c4762a1bSJed Brown   ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
510*c4762a1bSJed Brown   ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr);
511*c4762a1bSJed Brown   ierr = DMDestroy(dm);CHKERRQ(ierr);
512*c4762a1bSJed Brown   *dm  = idm;
513*c4762a1bSJed Brown   PetscFunctionReturn(0);
514*c4762a1bSJed Brown }
515*c4762a1bSJed Brown 
516*c4762a1bSJed Brown PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
517*c4762a1bSJed Brown {
518*c4762a1bSJed Brown   DM             idm;
519*c4762a1bSJed Brown   PetscInt       depth = 3, p;
520*c4762a1bSJed Brown   PetscMPIInt    rank;
521*c4762a1bSJed Brown   PetscErrorCode ierr;
522*c4762a1bSJed Brown 
523*c4762a1bSJed Brown   PetscFunctionBegin;
524*c4762a1bSJed Brown   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
525*c4762a1bSJed Brown   if (!rank) {
526*c4762a1bSJed Brown     switch (testNum) {
527*c4762a1bSJed Brown     case 0:
528*c4762a1bSJed Brown     {
529*c4762a1bSJed Brown       PetscInt    numPoints[2]         = {12, 2};
530*c4762a1bSJed Brown       PetscInt    coneSize[14]         = {8,8, 0,0,0,0,0,0,0,0,0,0,0,0};
531*c4762a1bSJed Brown       PetscInt    cones[16]            = {2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8};
532*c4762a1bSJed Brown       PetscInt    coneOrientations[16] = {0,0,0,0,0,0,0,0,  0,0, 0,0 ,0, 0, 0,0};
533*c4762a1bSJed Brown       PetscScalar vertexCoords[36]     = {-0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0,
534*c4762a1bSJed Brown                                           -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0,
535*c4762a1bSJed Brown                                            0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0,  0.5,1.0,1.0};
536*c4762a1bSJed Brown       PetscInt    markerPoints[52]     = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1};
537*c4762a1bSJed Brown       PetscInt    faultPoints[4]       = {3, 4, 7, 8};
538*c4762a1bSJed Brown 
539*c4762a1bSJed Brown       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
540*c4762a1bSJed Brown       ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
541*c4762a1bSJed Brown       for(p = 0; p < 8; ++p) {ierr = DMSetLabelValue(idm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
542*c4762a1bSJed Brown       for(p = 0; p < 4; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
543*c4762a1bSJed Brown     }
544*c4762a1bSJed Brown     break;
545*c4762a1bSJed Brown     case 1:
546*c4762a1bSJed Brown     {
547*c4762a1bSJed Brown       /* Cell Adjacency Graph:
548*c4762a1bSJed Brown         0 -- { 8, 13, 21, 24} --> 1
549*c4762a1bSJed Brown         0 -- {20, 21, 23, 24} --> 5 F
550*c4762a1bSJed Brown         1 -- {10, 15, 21, 24} --> 2
551*c4762a1bSJed Brown         1 -- {13, 14, 15, 24} --> 6
552*c4762a1bSJed Brown         2 -- {21, 22, 24, 25} --> 4 F
553*c4762a1bSJed Brown         3 -- {21, 24, 30, 35} --> 4
554*c4762a1bSJed Brown         3 -- {21, 24, 28, 33} --> 5
555*c4762a1bSJed Brown        */
556*c4762a1bSJed Brown       PetscInt    numPoints[2]         = {30, 7};
557*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};
558*c4762a1bSJed Brown       PetscInt    cones[56]            = { 8, 21, 20,  7, 13, 12, 23, 24,
559*c4762a1bSJed Brown                                           14, 15, 10,  9, 13,  8, 21, 24,
560*c4762a1bSJed Brown                                           15, 16, 11, 10, 24, 21, 22, 25,
561*c4762a1bSJed Brown                                           30, 29, 28, 21, 35, 24, 33, 34,
562*c4762a1bSJed Brown                                           24, 21, 30, 35, 25, 36, 31, 22,
563*c4762a1bSJed Brown                                           27, 20, 21, 28, 32, 33, 24, 23,
564*c4762a1bSJed Brown                                           15, 24, 13, 14, 19, 18, 17, 26};
565*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};
566*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,
567*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,
568*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,
569*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,
570*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};
571*c4762a1bSJed Brown       PetscInt    faultPoints[6]       = {20, 21, 22, 23, 24, 25};
572*c4762a1bSJed Brown 
573*c4762a1bSJed Brown       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
574*c4762a1bSJed Brown       ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
575*c4762a1bSJed Brown       for(p = 0; p < 6; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
576*c4762a1bSJed Brown     }
577*c4762a1bSJed Brown     break;
578*c4762a1bSJed Brown     case 2:
579*c4762a1bSJed Brown     {
580*c4762a1bSJed Brown       /* Buried fault edge */
581*c4762a1bSJed Brown       PetscInt    numPoints[2]         = {18, 4};
582*c4762a1bSJed Brown       PetscInt    coneSize[22]         = {8,8,8,8, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
583*c4762a1bSJed Brown       PetscInt    cones[32]            = { 4,  5,  8,  7, 13, 16, 17, 14,
584*c4762a1bSJed Brown                                            5,  6,  9,  8, 14, 17, 18, 15,
585*c4762a1bSJed Brown                                            7,  8, 11, 10, 16, 19, 20, 17,
586*c4762a1bSJed Brown                                            8,  9, 12, 11, 17, 20, 21, 18};
587*c4762a1bSJed Brown       PetscInt    coneOrientations[32] = {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};
588*c4762a1bSJed Brown       PetscScalar vertexCoords[54]     = {-2.0, -2.0,  0.0,  -2.0,  0.0,  0.0,  -2.0,  2.0,  0.0,   0.0, -2.0,  0.0,   0.0,  0.0,  0.0,   0.0,  2.0,  0.0,
589*c4762a1bSJed Brown                                            2.0, -2.0,  0.0,   2.0,  0.0,  0.0,   2.0,  2.0,  0.0,  -2.0, -2.0,  2.0,  -2.0,  0.0,  2.0,  -2.0,  2.0,  2.0,
590*c4762a1bSJed Brown                                            0.0, -2.0,  2.0,   0.0,  0.0,  2.0,   0.0,  2.0,  2.0,   2.0, -2.0,  2.0,   2.0,  0.0,  2.0,   2.0,  2.0,  2.0};
591*c4762a1bSJed Brown       PetscInt    faultPoints[4]       = {7, 8, 16, 17};
592*c4762a1bSJed Brown 
593*c4762a1bSJed Brown       ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
594*c4762a1bSJed Brown       ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
595*c4762a1bSJed Brown       for(p = 0; p < 4; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);}
596*c4762a1bSJed Brown     }
597*c4762a1bSJed Brown     break;
598*c4762a1bSJed Brown     default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum);
599*c4762a1bSJed Brown     }
600*c4762a1bSJed Brown   } else {
601*c4762a1bSJed Brown     PetscInt numPoints[4] = {0, 0, 0, 0};
602*c4762a1bSJed Brown 
603*c4762a1bSJed Brown     ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
604*c4762a1bSJed Brown     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
605*c4762a1bSJed Brown     ierr = DMCreateLabel(idm, "fault");CHKERRQ(ierr);
606*c4762a1bSJed Brown   }
607*c4762a1bSJed Brown   ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr);
608*c4762a1bSJed Brown   ierr = DMDestroy(dm);CHKERRQ(ierr);
609*c4762a1bSJed Brown   *dm  = idm;
610*c4762a1bSJed Brown   PetscFunctionReturn(0);
611*c4762a1bSJed Brown }
612*c4762a1bSJed Brown 
613*c4762a1bSJed Brown PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
614*c4762a1bSJed Brown {
615*c4762a1bSJed Brown   PetscInt       dim          = user->dim;
616*c4762a1bSJed Brown   PetscBool      cellSimplex  = user->cellSimplex, hasFault, hasFault2, hasParallelFault;
617*c4762a1bSJed Brown   PetscMPIInt    rank, size;
618*c4762a1bSJed Brown   PetscErrorCode ierr;
619*c4762a1bSJed Brown 
620*c4762a1bSJed Brown   PetscFunctionBegin;
621*c4762a1bSJed Brown   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
622*c4762a1bSJed Brown   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
623*c4762a1bSJed Brown   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
624*c4762a1bSJed Brown   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
625*c4762a1bSJed Brown   ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
626*c4762a1bSJed Brown   switch (dim) {
627*c4762a1bSJed Brown   case 2:
628*c4762a1bSJed Brown     if (cellSimplex) {
629*c4762a1bSJed Brown       ierr = CreateSimplex_2D(comm, user->testNum, dm);CHKERRQ(ierr);
630*c4762a1bSJed Brown     } else {
631*c4762a1bSJed Brown       ierr = CreateQuad_2D(comm, user->testNum, dm);CHKERRQ(ierr);
632*c4762a1bSJed Brown     }
633*c4762a1bSJed Brown     break;
634*c4762a1bSJed Brown   case 3:
635*c4762a1bSJed Brown     if (cellSimplex) {
636*c4762a1bSJed Brown       ierr = CreateSimplex_3D(comm, user, *dm);CHKERRQ(ierr);
637*c4762a1bSJed Brown     } else {
638*c4762a1bSJed Brown       ierr = CreateHex_3D(comm, user->testNum, dm);CHKERRQ(ierr);
639*c4762a1bSJed Brown     }
640*c4762a1bSJed Brown     break;
641*c4762a1bSJed Brown   default:
642*c4762a1bSJed Brown     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %d", dim);
643*c4762a1bSJed Brown   }
644*c4762a1bSJed Brown   ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_");CHKERRQ(ierr);
645*c4762a1bSJed Brown   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
646*c4762a1bSJed Brown   ierr = DMHasLabel(*dm, "fault", &hasFault);CHKERRQ(ierr);
647*c4762a1bSJed Brown   if (hasFault) {
648*c4762a1bSJed Brown     DM      dmHybrid = NULL;
649*c4762a1bSJed Brown     DMLabel faultLabel, faultBdLabel, hybridLabel;
650*c4762a1bSJed Brown 
651*c4762a1bSJed Brown     ierr = DMGetLabel(*dm, "fault", &faultLabel);CHKERRQ(ierr);
652*c4762a1bSJed Brown     ierr = DMGetLabel(*dm, "faultBd", &faultBdLabel);CHKERRQ(ierr);
653*c4762a1bSJed Brown     ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr);
654*c4762a1bSJed Brown     ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
655*c4762a1bSJed Brown     ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr);
656*c4762a1bSJed Brown     ierr = DMDestroy(dm);CHKERRQ(ierr);
657*c4762a1bSJed Brown     *dm  = dmHybrid;
658*c4762a1bSJed Brown   }
659*c4762a1bSJed Brown   ierr = DMHasLabel(*dm, "fault2", &hasFault2);CHKERRQ(ierr);
660*c4762a1bSJed Brown   if (hasFault2) {
661*c4762a1bSJed Brown     DM      dmHybrid = NULL;
662*c4762a1bSJed Brown     DMLabel faultLabel, faultBdLabel, hybridLabel;
663*c4762a1bSJed Brown 
664*c4762a1bSJed Brown     ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_");CHKERRQ(ierr);
665*c4762a1bSJed Brown     ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
666*c4762a1bSJed Brown     ierr = DMViewFromOptions(*dm, NULL, "-faulted_dm_view");CHKERRQ(ierr);
667*c4762a1bSJed Brown     ierr = DMGetLabel(*dm, "fault2", &faultLabel);CHKERRQ(ierr);
668*c4762a1bSJed Brown     ierr = DMGetLabel(*dm, "fault2Bd", &faultBdLabel);CHKERRQ(ierr);
669*c4762a1bSJed Brown     ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr);
670*c4762a1bSJed Brown     ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
671*c4762a1bSJed Brown     ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr);
672*c4762a1bSJed Brown     ierr = DMDestroy(dm);CHKERRQ(ierr);
673*c4762a1bSJed Brown     *dm  = dmHybrid;
674*c4762a1bSJed Brown   }
675*c4762a1bSJed Brown   if (user->testPartition && size > 1) {
676*c4762a1bSJed Brown     PetscPartitioner part;
677*c4762a1bSJed Brown     PetscInt *sizes  = NULL;
678*c4762a1bSJed Brown     PetscInt *points = NULL;
679*c4762a1bSJed Brown 
680*c4762a1bSJed Brown     if (!rank) {
681*c4762a1bSJed Brown       if (dim == 2 && cellSimplex && size == 2) {
682*c4762a1bSJed Brown         switch (user->testNum) {
683*c4762a1bSJed Brown         case 0: {
684*c4762a1bSJed Brown           PetscInt triSizes_p2[2]  = {1, 2};
685*c4762a1bSJed Brown           PetscInt triPoints_p2[3] = {0, 1, 2};
686*c4762a1bSJed Brown 
687*c4762a1bSJed Brown           ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr);
688*c4762a1bSJed Brown           ierr = PetscArraycpy(sizes,  triSizes_p2, 2);CHKERRQ(ierr);
689*c4762a1bSJed Brown           ierr = PetscArraycpy(points, triPoints_p2, 3);CHKERRQ(ierr);break;}
690*c4762a1bSJed Brown         default:
691*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);
692*c4762a1bSJed Brown         }
693*c4762a1bSJed Brown       } else if (dim == 2 && !cellSimplex && size == 2) {
694*c4762a1bSJed Brown         switch (user->testNum) {
695*c4762a1bSJed Brown         case 0: {
696*c4762a1bSJed Brown           PetscInt quadSizes_p2[2]  = {1, 2};
697*c4762a1bSJed Brown           PetscInt quadPoints_p2[3] = {0, 1, 2};
698*c4762a1bSJed Brown 
699*c4762a1bSJed Brown           ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr);
700*c4762a1bSJed Brown           ierr = PetscArraycpy(sizes,  quadSizes_p2, 2);CHKERRQ(ierr);
701*c4762a1bSJed Brown           ierr = PetscArraycpy(points, quadPoints_p2, 3);CHKERRQ(ierr);break;}
702*c4762a1bSJed Brown         case 2: {
703*c4762a1bSJed Brown           PetscInt quadSizes_p2[2]  = {1, 1};
704*c4762a1bSJed Brown           PetscInt quadPoints_p2[2] = {0, 1};
705*c4762a1bSJed Brown 
706*c4762a1bSJed Brown           ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr);
707*c4762a1bSJed Brown           ierr = PetscArraycpy(sizes,  quadSizes_p2, 2);CHKERRQ(ierr);
708*c4762a1bSJed Brown           ierr = PetscArraycpy(points, quadPoints_p2, 2);CHKERRQ(ierr);break;}
709*c4762a1bSJed Brown         default:
710*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);
711*c4762a1bSJed Brown         }
712*c4762a1bSJed Brown       } else if (dim == 3 && cellSimplex && size == 2) {
713*c4762a1bSJed Brown         switch (user->testNum) {
714*c4762a1bSJed Brown         case 0: {
715*c4762a1bSJed Brown           PetscInt tetSizes_p2[2]  = {1, 2};
716*c4762a1bSJed Brown           PetscInt tetPoints_p2[3] = {0, 1, 2};
717*c4762a1bSJed Brown 
718*c4762a1bSJed Brown           ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr);
719*c4762a1bSJed Brown           ierr = PetscArraycpy(sizes,  tetSizes_p2, 2);CHKERRQ(ierr);
720*c4762a1bSJed Brown           ierr = PetscArraycpy(points, tetPoints_p2, 3);CHKERRQ(ierr);break;}
721*c4762a1bSJed Brown         default:
722*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);
723*c4762a1bSJed Brown         }
724*c4762a1bSJed Brown       } else if (dim == 3 && !cellSimplex && size == 2) {
725*c4762a1bSJed Brown         switch (user->testNum) {
726*c4762a1bSJed Brown         case 0: {
727*c4762a1bSJed Brown           PetscInt hexSizes_p2[2]  = {1, 2};
728*c4762a1bSJed Brown           PetscInt hexPoints_p2[3] = {0, 1, 2};
729*c4762a1bSJed Brown 
730*c4762a1bSJed Brown           ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr);
731*c4762a1bSJed Brown           ierr = PetscArraycpy(sizes,  hexSizes_p2, 2);CHKERRQ(ierr);
732*c4762a1bSJed Brown           ierr = PetscArraycpy(points, hexPoints_p2, 3);CHKERRQ(ierr);break;}
733*c4762a1bSJed Brown         default:
734*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);
735*c4762a1bSJed Brown         }
736*c4762a1bSJed Brown       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition");
737*c4762a1bSJed Brown     }
738*c4762a1bSJed Brown     ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr);
739*c4762a1bSJed Brown     ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);CHKERRQ(ierr);
740*c4762a1bSJed Brown     ierr = PetscPartitionerShellSetPartition(part, size, sizes, points);CHKERRQ(ierr);
741*c4762a1bSJed Brown     ierr = PetscFree2(sizes, points);CHKERRQ(ierr);
742*c4762a1bSJed Brown   }
743*c4762a1bSJed Brown   {
744*c4762a1bSJed Brown     DM pdm = NULL;
745*c4762a1bSJed Brown 
746*c4762a1bSJed Brown     /* Distribute mesh over processes */
747*c4762a1bSJed Brown     ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr);
748*c4762a1bSJed Brown     if (pdm) {
749*c4762a1bSJed Brown       ierr = DMViewFromOptions(pdm, NULL, "-dm_view");CHKERRQ(ierr);
750*c4762a1bSJed Brown       ierr = DMDestroy(dm);CHKERRQ(ierr);
751*c4762a1bSJed Brown       *dm  = pdm;
752*c4762a1bSJed Brown     }
753*c4762a1bSJed Brown   }
754*c4762a1bSJed Brown   ierr = DMHasLabel(*dm, "pfault", &hasParallelFault);CHKERRQ(ierr);
755*c4762a1bSJed Brown   if (hasParallelFault) {
756*c4762a1bSJed Brown     DM      dmHybrid = NULL;
757*c4762a1bSJed Brown     DMLabel faultLabel, faultBdLabel, hybridLabel;
758*c4762a1bSJed Brown 
759*c4762a1bSJed Brown     ierr = DMGetLabel(*dm, "pfault", &faultLabel);CHKERRQ(ierr);
760*c4762a1bSJed Brown     ierr = DMGetLabel(*dm, "pfaultBd", &faultBdLabel);CHKERRQ(ierr);
761*c4762a1bSJed Brown     ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr);
762*c4762a1bSJed Brown     ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
763*c4762a1bSJed Brown     ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr);
764*c4762a1bSJed Brown     ierr = DMDestroy(dm);CHKERRQ(ierr);
765*c4762a1bSJed Brown     *dm  = dmHybrid;
766*c4762a1bSJed Brown   }
767*c4762a1bSJed Brown   ierr = PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh");CHKERRQ(ierr);
768*c4762a1bSJed Brown   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
769*c4762a1bSJed Brown   ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
770*c4762a1bSJed Brown   PetscFunctionReturn(0);
771*c4762a1bSJed Brown }
772*c4762a1bSJed Brown 
773*c4762a1bSJed Brown int main(int argc, char **argv)
774*c4762a1bSJed Brown {
775*c4762a1bSJed Brown   DM             dm;
776*c4762a1bSJed Brown   AppCtx         user;                 /* user-defined work context */
777*c4762a1bSJed Brown   PetscErrorCode ierr;
778*c4762a1bSJed Brown 
779*c4762a1bSJed Brown   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
780*c4762a1bSJed Brown   ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
781*c4762a1bSJed Brown   ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
782*c4762a1bSJed Brown   ierr = DMDestroy(&dm);CHKERRQ(ierr);
783*c4762a1bSJed Brown   ierr = PetscFinalize();
784*c4762a1bSJed Brown   return ierr;
785*c4762a1bSJed Brown }
786*c4762a1bSJed Brown 
787*c4762a1bSJed Brown /*TEST
788*c4762a1bSJed Brown   testset:
789*c4762a1bSJed Brown     args: -orig_dm_plex_check_symmetry -orig_dm_plex_check_skeleton -orig_dm_plex_check_faces \
790*c4762a1bSJed Brown           -dm_view ascii::ascii_info_detail -dm_plex_check_symmetry -dm_plex_check_skeleton -dm_plex_check_faces
791*c4762a1bSJed Brown     # 2D Simplex
792*c4762a1bSJed Brown     test:
793*c4762a1bSJed Brown       suffix: tri_0
794*c4762a1bSJed Brown       args: -dim 2
795*c4762a1bSJed Brown     test:
796*c4762a1bSJed Brown       suffix: tri_1
797*c4762a1bSJed Brown       nsize: 2
798*c4762a1bSJed Brown       args: -dim 2
799*c4762a1bSJed Brown     test:
800*c4762a1bSJed Brown       suffix: tri_t1_0
801*c4762a1bSJed Brown       args: -dim 2 -test_num 1
802*c4762a1bSJed Brown     # 3D Simplex
803*c4762a1bSJed Brown     test:
804*c4762a1bSJed Brown       suffix: tet_0
805*c4762a1bSJed Brown       args: -dim 3
806*c4762a1bSJed Brown     test:
807*c4762a1bSJed Brown       suffix: tet_1
808*c4762a1bSJed Brown       nsize: 2
809*c4762a1bSJed Brown       args: -dim 3
810*c4762a1bSJed Brown     test:
811*c4762a1bSJed Brown       suffix: tet_t1_0
812*c4762a1bSJed Brown       args: -dim 3 -test_num 1
813*c4762a1bSJed Brown 
814*c4762a1bSJed Brown   testset:
815*c4762a1bSJed Brown     args: -orig_dm_plex_check_symmetry -orig_dm_plex_check_skeleton -orig_dm_plex_check_faces \
816*c4762a1bSJed Brown           -dm_view ascii::ascii_info_detail -dm_plex_check_symmetry -dm_plex_check_skeleton -dm_plex_check_faces
817*c4762a1bSJed Brown     # 2D Quads
818*c4762a1bSJed Brown     test:
819*c4762a1bSJed Brown       suffix: quad_0
820*c4762a1bSJed Brown       args: -dim 2 -cell_simplex 0
821*c4762a1bSJed Brown     test:
822*c4762a1bSJed Brown       suffix: quad_1
823*c4762a1bSJed Brown       nsize: 2
824*c4762a1bSJed Brown       args: -dim 2 -cell_simplex 0
825*c4762a1bSJed Brown     test:
826*c4762a1bSJed Brown       suffix: quad_t1_0
827*c4762a1bSJed Brown       args: -dim 2 -cell_simplex 0 -test_num 1 \
828*c4762a1bSJed Brown             -faulted_dm_plex_check_symmetry -faulted_dm_plex_check_skeleton -faulted_dm_plex_check_faces
829*c4762a1bSJed Brown     # 3D Hex
830*c4762a1bSJed Brown     test:
831*c4762a1bSJed Brown       suffix: hex_0
832*c4762a1bSJed Brown       args: -dim 3 -cell_simplex 0
833*c4762a1bSJed Brown     test:
834*c4762a1bSJed Brown       suffix: hex_1
835*c4762a1bSJed Brown       nsize: 2
836*c4762a1bSJed Brown       args: -dim 3 -cell_simplex 0
837*c4762a1bSJed Brown     test:
838*c4762a1bSJed Brown       suffix: hex_t1_0
839*c4762a1bSJed Brown       args: -dim 3 -cell_simplex 0 -test_num 1
840*c4762a1bSJed Brown     test:
841*c4762a1bSJed Brown       suffix: hex_t2_0
842*c4762a1bSJed Brown       args: -dim 3 -cell_simplex 0 -test_num 2
843*c4762a1bSJed Brown 
844*c4762a1bSJed Brown TEST*/
845