xref: /petsc/src/dm/impls/plex/tests/ex5.c (revision fbf9dbe564678ed6eff1806adbc4c4f01b9743f4)
1 static char help[] = "Tests for creation of hybrid meshes\n\n";
2 
3 /* TODO
4  - Propagate hybridSize with distribution
5  - Test with multiple fault segments
6  - Test with embedded fault
7  - Test with multiple faults
8  - Move over all PyLith tests?
9 */
10 
11 #include <petscdmplex.h>
12 #include <petscds.h>
13 #include <petsc/private/dmpleximpl.h>
14 /* List of test meshes
15 
16 Triangle
17 --------
18 Test 0:
19 Two triangles sharing a face
20 
21         4
22       / | \
23      8  |  9
24     /   |   \
25    2  0 7 1  5
26     \   |   /
27      6  |  10
28       \ | /
29         3
30 
31 should become two triangles separated by a zero-volume cell with 4 vertices
32 
33         5--16--8              4--12--6 3
34       / |      | \          / |      | | \
35     11  |      |  12       9  |      | |  4
36     /   |      |   \      /   |      | |   \
37    3  0 10  2 14 1  6    2  0 8  1  10 6 0  1
38     \   |      |   /      \   |      | |   /
39      9  |      |  13       7  |      | |  5
40       \ |      | /          \ |      | | /
41         4--15--7              3--11--5 2
42 
43 Test 1:
44 Four triangles sharing two faces which are oriented against each other
45 
46           9
47          / \
48         /   \
49       17  2  16
50       /       \
51      /         \
52     8-----15----5
53      \         /|\
54       \       / | \
55       18  3  12 |  14
56         \   /   |   \
57          \ /    |    \
58           4  0 11  1  7
59            \    |    /
60             \   |   /
61             10  |  13
62               \ | /
63                \|/
64                 6
65 
66 Fault mesh
67 
68 0 --> 0
69 1 --> 1
70 2 --> 2
71 3 --> 3
72 4 --> 5
73 5 --> 6
74 6 --> 8
75 7 --> 11
76 8 --> 15
77 
78        2
79        |
80   6----8----4
81        |    |
82        3    |
83           0-7-1
84             |
85             |
86             5
87 
88 should become four triangles separated by two zero-volume cells with 4 vertices
89 
90           11
91           / \
92          /   \
93         /     \
94       22   2   21
95       /         \
96      /           \
97    10-----20------7
98 28  |     5    26/ \
99    14----25----12   \
100      \         /|   |\
101       \       / |   | \
102       23  3  17 |   |  19
103         \   /   |   |   \
104          \ /    |   |    \
105           6  0 24 4 16 1  9
106            \    |   |    /
107             \   |   |   /
108             15  |   |  18
109               \ |   | /
110                \|   |/
111                13---8
112                  27
113 
114 Test 2:
115 Six triangles sharing one face
116 
117 11-----12------13
118  |     /|\     |
119  | 1  / | \ 4  |
120  |   /  |  \   |
121  |  /   |   \  |
122  | /    |    \ |
123  |/     |     \|
124  9  2   |   5  10
125  |\     |     /|
126  | \    |    / |
127  |  \   |   /  |
128  |   \  |  /   |
129  | 0  \ | / 3  |
130  |     \|/     |
131  6------7------8
132 
133 Test 3:
134 This is Test 2 on two processes. After the fault, we have
135 
136  6--12--7    7--20-10--16-8
137  |     /|    |     |\     |
138  | 1  / |    |     | \  1 |
139 13  11  |    |     |  17  15
140  |  /   |    |     |   \  |
141  | /    |    |     |    \ |
142  |/     |    |     |     \|
143  5   2  14  11  3 18  2   6
144  |\     |    |     |     /|
145  | \    |    |     |    / |
146  |  \   |    |     |   /  |
147 10   9  |    |     |  14  13
148  | 0  \ |    |     | /  0 |
149  |     \|    |     |/     |
150  3---8--4    4--19-9--12--5
151 
152 Test 4:
153 This is Test 2 on six processes. After the fault, we have
154 
155 Test 5:
156 
157   Fault only on points 2 and 5:
158 
159         6
160       / | \
161     13  |  17
162     /  15   \
163    7  0 | 1  9
164    |\   |   /|
165    | 14 | 16 |
166    |  \ | /  |
167  18| 2  8  3 |21
168    |  / | \  |
169    | 19 | 20 |
170    |/   |   \|
171   10  4 | 5  12
172     \  23   /
173     22  |  24
174       \ | /
175        11
176 
177 Tetrahedron
178 -----------
179 Test 0:
180 Two tets sharing a face
181 
182  cell   5 _______    cell
183  0    / | \      \       1
184     16  |  18     22
185     /8 19 10\      \
186    2-15-|----4--21--6
187     \  9| 7 /      /
188     14  |  17     20
189       \ | /      /
190         3-------
191 
192 should become two tetrahedrons separated by a zero-volume cell with 3 faces/3 edges/6 vertices
193 
194  cell   6 ___36___10______    cell
195  0    / | \        |\      \     1
196     24  |  26      | 32     30
197     /12 27 14\    33  \      \
198    3-23-|----5--35-|---9--29--7
199     \ 13| 11/      |18 /      /
200     22  |  25      | 31     28
201       \ | /        |/      /
202         4----34----8------
203          cell 2
204 
205 In parallel,
206 
207  cell   5 ___28____8      4______    cell
208  0    / | \        |\     |\      \     0
209     19  |   21     | 24   | 13  6  11
210     /10 22 12\    25  \   |8 \      \
211    2-18-|----4--27-|---7  14  3--10--1
212     \ 11| 9 /      |13 /  |  /      /
213     17  |  20      | 23   | 12  5  9
214       \ | /        |/     |/      /
215         3----26----6      2------
216          cell 1
217 
218 Test 1:
219 Four tets sharing two faces
220 
221 Cells:    0-3,4-5
222 Vertices: 6-15
223 Faces:    16-29,30-34
224 Edges:    35-52,53-56
225 
226 Quadrilateral
227 -------------
228 Test 0:
229 Two quads sharing a face
230 
231    5--10---4--14---7
232    |       |       |
233   11   0   9   1  13
234    |       |       |
235    2---8---3--12---6
236 
237 should become two quads separated by a zero-volume cell with 4 vertices
238 
239    6--13---5-20-10--17---8    5--10---4-14--7  4---7---2
240    |       |     |       |    |       |     |  |       |
241   14   0  12  2 18   1  16   11   0   9  1 12  8   0   6
242    |       |     |       |    |       |     |  |       |
243    3--11---4-19--9--15---7    2---8---3-13--6  3---5---1
244 
245 Test 1:
246 
247 Original mesh with 9 cells,
248 
249   9-----10-----11-----12
250   |      |     ||      |
251   |      |     ||      |
252   |   0  |  1  ||  2   |
253   |      |     ||      |
254  13-----14-----15-----16
255   |      |     ||      |
256   |      |     ||      |
257   |  3   |  4  ||  5   |
258   |      |     ||      |
259  17-----18-----19=====20
260   |      |      |      |
261   |      |      |      |
262   |  6   |  7   |  8   |
263   |      |      |      |
264  21-----22-----23-----24
265 
266 After first fault,
267 
268  12 ----13 ----14-28 ----15
269   |      |      |  |      |
270   |  0   |  1   | 9|  2   |
271   |      |      |  |      |
272   |      |      |  |      |
273  16 ----17 ----18-29 ----19
274   |      |      |  |      |
275   |  3   |  4   |10|  5   |
276   |      |      |  |      |
277   |      |      |  |      |
278  20 ----21-----22-30 ----23
279   |      |      |  \--11- |
280   |  6   |  7   |     8   |
281   |      |      |         |
282   |      |      |         |
283  24 ----25 ----26--------27
284 
285 After second fault,
286 
287  14 ----15 ----16-30 ----17
288   |      |      |  |      |
289   |  0   |  1   | 9|  2   |
290   |      |      |  |      |
291   |      |      |  |      |
292  18 ----19 ----20-31 ----21
293   |      |      |  |      |
294   |  3   |  4   |10|  5   |
295   |      |      |  |      |
296   |      |      |  |      |
297  33 ----34-----24-32 ----25
298   |  12  | 13 / |  \-11-- |
299  22 ----23---/  |         |
300   |      |      |         |
301   |  6   |   7  |     8   |
302   |      |      |         |
303   |      |      |         |
304  26 ----27 ----28--------29
305 
306  Test 2:
307  Two quads sharing a face in parallel
308 
309     4---7---3  2---8---4
310     |       |  |       |
311     8   0   6  5   0   7
312     |       |  |       |
313     1---5---2  1---6---3
314 
315  should become two quads separated by a zero-volume cell with 4 vertices
316 
317      4---7---3  3-14--7--11---5
318      |       |  |     |       |
319      8   0   6  8  1  12  0   10
320      |       |  |     |       |
321      1---5---2  2-13--6---9---4
322 
323  Test 3:
324  Like Test 2, but with different partition
325 
326      5--10---4-14--7   2---8---4
327      |       |     |   |       |
328     11   0   9  1  12  5   0   7
329      |       |     |   |       |
330      2---8---3-13--6   1---6---3
331 
332 Hexahedron
333 ----------
334 Test 0:
335 Two hexes sharing a face
336 
337 cell   9-----31------8-----42------13 cell
338 0     /|            /|            /|     1
339     32 |   15      30|   21      41|
340     /  |          /  |          /  |
341    6-----29------7-----40------12  |
342    |   |     18  |   |     24  |   |
343    |  36         |  35         |   44
344    |19 |         |17 |         |23 |
345   33   |  16    34   |   22   43   |
346    |   5-----27--|---4-----39--|---11
347    |  /          |  /          |  /
348    | 28   14     | 26    20    | 38
349    |/            |/            |/
350    2-----25------3-----37------10
351 
352 should become two hexes separated by a zero-volume cell with 8 vertices
353 
354                          cell 2
355 cell  10-----41------9-----62------18----52------14 cell
356 0     /|            /|            /|            /|     1
357     42 |   20      40|  32       56|   26      51|
358     /  |          /  |          /  |          /  |
359    7-----39------8-----61------17--|-50------13  |
360    |   |     23  |   |         |   |     29  |   |
361    |  46         |  45         |   58        |   54
362    |24 |         |22 |         |30 |         |28 |
363   43   |  21    44   |        57   |   27   53   |
364    |   6-----37--|---5-----60--|---16----49--|---12
365    |  /          |  /          |  /          |  /
366    | 38   19     | 36   31     | 55    25    | 48
367    |/            |/            |/            |/
368    3-----35------4-----59------15----47------11
369 
370 In parallel,
371 
372                          cell 2
373 cell   9-----31------8-----44------13     8----20------4  cell
374 0     /|            /|            /|     /|           /|     1
375     32 |   15      30|  22       38|   24 |  10      19|
376     /  |          /  |          /  |   /  |         /  |
377    6-----29------7-----43------12  |  7----18------3   |
378    |   |     18  |   |         |   |  |   |    13  |   |
379    |  36         |  35         |   40 |  26        |   22
380    |19 |         |17 |         |20 |  |14 |        |12 |
381   33   |  16    34   |        39   |  25  |  11   21   |
382    |   5-----27--|---4-----42--|---11 |   6----17--|---2
383    |  /          |  /          |  /   |  /         |  /
384    | 28   14     | 26   21     | 37   |23     9    | 16
385    |/            |/            |/     |/           |/
386    2-----25------3-----41------10     5----15------1
387 
388 Test 1:
389 
390 */
391 
392 typedef struct {
393   PetscInt  debug;          /* The debugging level */
394   PetscInt  dim;            /* The topological mesh dimension */
395   PetscBool cellSimplex;    /* Use simplices or hexes */
396   PetscBool testPartition;  /* Use a fixed partitioning for testing */
397   PetscInt  testNum;        /* The particular mesh to test */
398   PetscInt  cohesiveFields; /* The number of cohesive fields */
399 } AppCtx;
400 
401 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
402 {
403   PetscFunctionBegin;
404   options->debug          = 0;
405   options->dim            = 2;
406   options->cellSimplex    = PETSC_TRUE;
407   options->testPartition  = PETSC_TRUE;
408   options->testNum        = 0;
409   options->cohesiveFields = 1;
410 
411   PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
412   PetscCall(PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL, 0));
413   PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL, 1, 3));
414   PetscCall(PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL));
415   PetscCall(PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL));
416   PetscCall(PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL, 0));
417   PetscCall(PetscOptionsBoundedInt("-cohesive_fields", "The number of cohesive fields", "ex5.c", options->cohesiveFields, &options->cohesiveFields, NULL, 0));
418   PetscOptionsEnd();
419   PetscFunctionReturn(PETSC_SUCCESS);
420 }
421 
422 static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
423 {
424   DM          idm;
425   PetscInt    p;
426   PetscMPIInt rank;
427 
428   PetscFunctionBegin;
429   PetscCallMPI(MPI_Comm_rank(comm, &rank));
430   if (rank == 0) {
431     switch (testNum) {
432     case 0: {
433       PetscInt    numPoints[2]        = {4, 2};
434       PetscInt    coneSize[6]         = {3, 3, 0, 0, 0, 0};
435       PetscInt    cones[6]            = {2, 3, 4, 5, 4, 3};
436       PetscInt    coneOrientations[6] = {0, 0, 0, 0, 0, 0};
437       PetscScalar vertexCoords[8]     = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
438       PetscInt    markerPoints[8]     = {2, 1, 3, 1, 4, 1, 5, 1};
439       PetscInt    faultPoints[2]      = {3, 4};
440 
441       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
442       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
443       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
444       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
445       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
446     } break;
447     case 1: {
448       PetscInt    numPoints[2]         = {6, 4};
449       PetscInt    coneSize[10]         = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0};
450       PetscInt    cones[12]            = {4, 6, 5, 5, 6, 7, 8, 5, 9, 8, 4, 5};
451       PetscInt    coneOrientations[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
452       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};
453       PetscInt    markerPoints[6]      = {4, 1, 6, 1, 8, 1};
454       PetscInt    faultPoints[3]       = {5, 6, 8};
455 
456       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
457       for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
458       for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
459       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
460       PetscCall(DMSetLabelValue(*dm, "material", 3, 1));
461       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
462       PetscCall(DMSetLabelValue(*dm, "material", 2, 2));
463     } break;
464     case 2:
465     case 3:
466     case 4: {
467       PetscInt    numPoints[2]         = {8, 6};
468       PetscInt    coneSize[14]         = {3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0};
469       PetscInt    cones[18]            = {6, 7, 9, 9, 12, 11, 7, 12, 9, 7, 8, 10, 10, 13, 12, 7, 10, 12};
470       PetscInt    coneOrientations[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
471       PetscScalar vertexCoords[16]     = {
472         -1., -1., 0., -1., 1., -1., -1., 0., 1., 0., -1., 1., 0., 1., 1., 1.,
473       };
474       PetscInt markerPoints[16] = {6, 1, 7, 1, 8, 1, 9, 1, 10, 1, 11, 1, 12, 1, 13, 1};
475       PetscInt faultPoints[2]   = {7, 12};
476 
477       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
478       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
479       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
480       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
481       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
482       PetscCall(DMSetLabelValue(*dm, "material", 4, 2));
483       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
484       for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
485       if (testNum == 2)
486         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
487       if (testNum == 3 || testNum == 4)
488         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1));
489     } break;
490     case 5: {
491       PetscInt    numPoints[2]         = {7, 6};
492       PetscInt    coneSize[13]         = {3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0};
493       PetscInt    cones[18]            = {6, 7, 8, 8, 9, 6, 7, 10, 8, 9, 8, 12, 8, 10, 11, 11, 12, 8};
494       PetscInt    coneOrientations[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
495       PetscScalar vertexCoords[14]     = {0.0, 2.0, -1.0, 1.0, 0.0, 0.0, 1.0, 1.0, -1.0, -1.0, 0.0, -2.0, 1.0, -1.0};
496       PetscInt    faultPoints[2]       = {8, 11};
497       PetscInt    faultBdPoints[1]     = {8};
498 
499       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
500       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
501       for (p = 0; p < 1; ++p) PetscCall(DMSetLabelValue(*dm, "faultBd", faultBdPoints[p], 1));
502       PetscCall(DMSetLabelValue(*dm, "material", 0, 0));
503       PetscCall(DMSetLabelValue(*dm, "material", 2, 0));
504       PetscCall(DMSetLabelValue(*dm, "material", 4, 0));
505       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
506       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
507       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
508     } break;
509     default:
510       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
511     }
512   } else {
513     PetscInt numPoints[3] = {0, 0, 0};
514 
515     PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL));
516     if (testNum == 3 || testNum == 4) PetscCall(DMCreateLabel(*dm, "pfault"));
517     else PetscCall(DMCreateLabel(*dm, "fault"));
518   }
519   PetscCall(DMPlexInterpolate(*dm, &idm));
520   PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view"));
521   PetscCall(DMDestroy(dm));
522   *dm = idm;
523   PetscFunctionReturn(PETSC_SUCCESS);
524 }
525 
526 static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, AppCtx *user, DM dm)
527 {
528   PetscInt    depth = 3, testNum = user->testNum, p;
529   PetscMPIInt rank;
530 
531   PetscFunctionBegin;
532   PetscCallMPI(MPI_Comm_rank(comm, &rank));
533   if (rank == 0) {
534     switch (testNum) {
535     case 0: {
536       PetscInt    numPoints[4]         = {5, 7, 9, 2};
537       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};
538       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};
539       PetscInt    coneOrientations[47] = {0, 0, 0, 0, 0, -2, 2, 2, 0, -1, -1, 0, -1, -1, 0, -1, -1, 0, 0, 0, 0, 0, -1, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
540       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};
541       PetscInt    markerPoints[20]     = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1};
542       PetscInt    faultPoints[3]       = {3, 4, 5};
543 
544       PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
545       for (p = 0; p < 10; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
546       for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1));
547       PetscCall(DMSetLabelValue(dm, "material", 0, 1));
548       PetscCall(DMSetLabelValue(dm, "material", 1, 2));
549     } break;
550     case 1: {
551       PetscInt    numPoints[4]         = {6, 13, 12, 4};
552       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};
553       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,
554                                           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};
555       PetscInt    coneOrientations[78] = {0, 0, 0, 0,  -2, 1,  0, 2,  0, 0,  -3, 0,  0,  -3, -1, 0, 0, 0, 0, 0, 0, -1, -1, 0, -1, -1, -1, -1, 0, 0, 0, 0, 0, -1, 0, -1, -1, 0, 0,
556                                           0, 0, 0, -1, -1, -1, 0, -1, 0, -1, -1, -1, -1, 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};
557       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};
558       PetscInt    markerPoints[14]     = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1};
559       PetscInt    faultPoints[4]       = {5, 6, 7, 8};
560 
561       PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
562       for (p = 0; p < 7; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
563       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1));
564       PetscCall(DMSetLabelValue(dm, "material", 0, 1));
565       PetscCall(DMSetLabelValue(dm, "material", 1, 2));
566     } break;
567     default:
568       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
569     }
570   } else {
571     PetscInt numPoints[4] = {0, 0, 0, 0};
572 
573     PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL));
574     PetscCall(DMCreateLabel(dm, "fault"));
575   }
576   PetscFunctionReturn(PETSC_SUCCESS);
577 }
578 
579 static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
580 {
581   DM          idm;
582   PetscInt    p;
583   PetscMPIInt rank;
584 
585   PetscFunctionBegin;
586   PetscCallMPI(MPI_Comm_rank(comm, &rank));
587   if (rank == 0) {
588     switch (testNum) {
589     case 0:
590     case 2:
591     case 3: {
592       PetscInt    numPoints[2]        = {6, 2};
593       PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
594       PetscInt    cones[8]            = {2, 3, 4, 5, 3, 6, 7, 4};
595       PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
596       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};
597       PetscInt    markerPoints[12]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1};
598       PetscInt    faultPoints[2]      = {3, 4};
599 
600       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
601       for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
602       if (testNum == 0)
603         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
604       if (testNum == 2 || testNum == 3)
605         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1));
606       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
607       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
608     } break;
609     case 1: {
610       PetscInt    numPoints[2]         = {16, 9};
611       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};
612       PetscInt    cones[36]            = {9, 13, 14, 10, 10, 14, 15, 11, 11, 15, 16, 12, 13, 17, 18, 14, 14, 18, 19, 15, 15, 19, 20, 16, 17, 21, 22, 18, 18, 22, 23, 19, 19, 23, 24, 20};
613       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};
614       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, -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};
615       PetscInt    faultPoints[4]       = {11, 15, 19, 20};
616       PetscInt    fault2Points[2]      = {17, 18};
617 
618       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
619       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
620       for (p = 3; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "faultBd", faultPoints[p], 1));
621       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault2", fault2Points[p], 1));
622       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
623       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
624       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
625       PetscCall(DMSetLabelValue(*dm, "material", 3, 1));
626       PetscCall(DMSetLabelValue(*dm, "material", 4, 1));
627       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
628       PetscCall(DMSetLabelValue(*dm, "material", 6, 2));
629       PetscCall(DMSetLabelValue(*dm, "material", 7, 2));
630       PetscCall(DMSetLabelValue(*dm, "material", 8, 2));
631     } break;
632     default:
633       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
634     }
635   } else {
636     PetscInt numPoints[3] = {0, 0, 0};
637 
638     PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL));
639     if (testNum == 2 || testNum == 3) PetscCall(DMCreateLabel(*dm, "pfault"));
640     else PetscCall(DMCreateLabel(*dm, "fault"));
641   }
642   PetscCall(DMPlexInterpolate(*dm, &idm));
643   PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view"));
644   PetscCall(DMDestroy(dm));
645   *dm = idm;
646   PetscFunctionReturn(PETSC_SUCCESS);
647 }
648 
649 static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
650 {
651   DM          idm;
652   PetscInt    depth = 3, p;
653   PetscMPIInt rank;
654 
655   PetscFunctionBegin;
656   PetscCallMPI(MPI_Comm_rank(comm, &rank));
657   if (rank == 0) {
658     switch (testNum) {
659     case 0: {
660       PetscInt    numPoints[2]         = {12, 2};
661       PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
662       PetscInt    cones[16]            = {2, 5, 4, 3, 6, 7, 8, 9, 3, 4, 11, 10, 7, 12, 13, 8};
663       PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
664       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, -0.5, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, -0.5, 1.0, 1.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0};
665       PetscInt    markerPoints[52]     = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1, 8, 1, 9, 1};
666       PetscInt    faultPoints[4]       = {3, 4, 7, 8};
667 
668       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
669       PetscCall(DMPlexInterpolate(*dm, &idm));
670       for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(idm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
671       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
672       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
673       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
674     } break;
675     case 1: {
676       /* Cell Adjacency Graph:
677         0 -- { 8, 13, 21, 24} --> 1
678         0 -- {20, 21, 23, 24} --> 5 F
679         1 -- {10, 15, 21, 24} --> 2
680         1 -- {13, 14, 15, 24} --> 6
681         2 -- {21, 22, 24, 25} --> 4 F
682         3 -- {21, 24, 30, 35} --> 4
683         3 -- {21, 24, 28, 33} --> 5
684        */
685       PetscInt numPoints[2] = {30, 7};
686       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};
687       PetscInt cones[56] = {8, 21, 20, 7, 13, 12, 23, 24, 14, 15, 10, 9, 13, 8, 21, 24, 15, 16, 11, 10, 24, 21, 22, 25, 30, 29, 28, 21, 35, 24, 33, 34, 24, 21, 30, 35, 25, 36, 31, 22, 27, 20, 21, 28, 32, 33, 24, 23, 15, 24, 13, 14, 19, 18, 17, 26};
688       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};
689       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,  -2.0, -1.0, 0.0, -3.0, 0.0, 0.0, -2.0, 1.0, 0.0, -2.0, 2.0, 0.0,
690                                        -2.0, -1.0, 2.0,  -3.0, 0.0,  2.0,  -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, 0.0,  2.0, 0.0, 0.0,  0.0, 2.0,
691                                        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,  2.0,  -1.0, 0.0, 3.0,  0.0, 0.0, 2.0,  1.0, 0.0, 2.0,  2.0, 0.0};
692       PetscInt    faultPoints[6]    = {20, 21, 22, 23, 24, 25};
693 
694       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
695       PetscCall(DMPlexInterpolate(*dm, &idm));
696       for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
697       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
698       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
699       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
700       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
701       PetscCall(DMSetLabelValue(*dm, "material", 4, 2));
702       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
703       PetscCall(DMSetLabelValue(*dm, "material", 6, 2));
704     } break;
705     case 2: {
706       /* Buried fault edge */
707       PetscInt    numPoints[2]         = {18, 4};
708       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};
709       PetscInt    cones[32]            = {4, 5, 8, 7, 13, 16, 17, 14, 5, 6, 9, 8, 14, 17, 18, 15, 7, 8, 11, 10, 16, 19, 20, 17, 8, 9, 12, 11, 17, 20, 21, 18};
710       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};
711       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, 2.0, -2.0, 0.0, 2.0, 0.0, 0.0, 2.0, 2.0, 0.0,
712                                           -2.0, -2.0, 2.0, -2.0, 0.0, 2.0, -2.0, 2.0, 2.0, 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};
713       PetscInt    faultPoints[4]       = {7, 8, 16, 17};
714 
715       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
716       PetscCall(DMPlexInterpolate(*dm, &idm));
717       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
718       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
719       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
720       PetscCall(DMSetLabelValue(*dm, "material", 2, 2));
721       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
722     } break;
723     default:
724       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
725     }
726   } else {
727     PetscInt numPoints[4] = {0, 0, 0, 0};
728 
729     PetscCall(DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL));
730     PetscCall(DMPlexInterpolate(*dm, &idm));
731     PetscCall(DMCreateLabel(idm, "fault"));
732   }
733   PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view"));
734   PetscCall(DMDestroy(dm));
735   *dm = idm;
736   PetscFunctionReturn(PETSC_SUCCESS);
737 }
738 
739 static PetscErrorCode CreateFaultLabel(DM dm)
740 {
741   DMLabel  label;
742   PetscInt dim, h, pStart, pEnd, pMax, p;
743 
744   PetscFunctionBegin;
745   PetscCall(DMGetDimension(dm, &dim));
746   PetscCall(DMCreateLabel(dm, "cohesive"));
747   PetscCall(DMGetLabel(dm, "cohesive", &label));
748   for (h = 0; h <= dim; ++h) {
749     PetscCall(DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax));
750     PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd));
751     for (p = pMax; p < pEnd; ++p) PetscCall(DMLabelSetValue(label, p, 1));
752   }
753   PetscFunctionReturn(PETSC_SUCCESS);
754 }
755 
756 static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user)
757 {
758   PetscFE  fe;
759   DMLabel  fault;
760   PetscInt dim, Ncf = user->cohesiveFields, f;
761 
762   PetscFunctionBegin;
763   PetscCall(DMGetDimension(dm, &dim));
764   PetscCall(DMGetLabel(dm, "cohesive", &fault));
765   PetscCall(DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD));
766 
767   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe));
768   PetscCall(PetscFESetName(fe, "displacement"));
769   PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
770   PetscCall(PetscFEDestroy(&fe));
771 
772   if (Ncf > 0) {
773     PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe));
774     PetscCall(PetscFESetName(fe, "fault traction"));
775     PetscCall(DMAddField(dm, fault, (PetscObject)fe));
776     PetscCall(PetscFEDestroy(&fe));
777   }
778   for (f = 1; f < Ncf; ++f) {
779     char name[256], opt[256];
780 
781     PetscCall(PetscSNPrintf(name, 256, "fault field %" PetscInt_FMT, f));
782     PetscCall(PetscSNPrintf(opt, 256, "faultfield_%" PetscInt_FMT "_", f));
783     PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe));
784     PetscCall(PetscFESetName(fe, name));
785     PetscCall(DMAddField(dm, fault, (PetscObject)fe));
786     PetscCall(PetscFEDestroy(&fe));
787   }
788 
789   PetscCall(DMCreateDS(dm));
790   PetscFunctionReturn(PETSC_SUCCESS);
791 }
792 
793 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
794 {
795   PetscInt    dim         = user->dim;
796   PetscBool   cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault;
797   PetscMPIInt rank, size;
798   DMLabel     matLabel;
799 
800   PetscFunctionBegin;
801   PetscCallMPI(MPI_Comm_rank(comm, &rank));
802   PetscCallMPI(MPI_Comm_size(comm, &size));
803   PetscCall(DMCreate(comm, dm));
804   PetscCall(DMSetType(*dm, DMPLEX));
805   PetscCall(DMSetDimension(*dm, dim));
806   switch (dim) {
807   case 2:
808     if (cellSimplex) {
809       PetscCall(CreateSimplex_2D(comm, user->testNum, dm));
810     } else {
811       PetscCall(CreateQuad_2D(comm, user->testNum, dm));
812     }
813     break;
814   case 3:
815     if (cellSimplex) {
816       PetscCall(CreateSimplex_3D(comm, user, *dm));
817     } else {
818       PetscCall(CreateHex_3D(comm, user->testNum, dm));
819     }
820     break;
821   default:
822     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %" PetscInt_FMT, dim);
823   }
824   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "orig_"));
825   PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
826   PetscCall(DMSetFromOptions(*dm));
827   PetscCall(DMGetLabel(*dm, "material", &matLabel));
828   if (matLabel) PetscCall(DMPlexLabelComplete(*dm, matLabel));
829   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
830   PetscCall(DMHasLabel(*dm, "fault", &hasFault));
831   if (hasFault) {
832     DM      dmHybrid = NULL, dmInterface = NULL;
833     DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel;
834 
835     PetscCall(DMGetLabel(*dm, "fault", &faultLabel));
836     PetscCall(DMGetLabel(*dm, "faultBd", &faultBdLabel));
837     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid));
838     PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD));
839     PetscCall(DMLabelDestroy(&hybridLabel));
840     PetscCall(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD));
841     PetscCall(DMLabelDestroy(&splitLabel));
842     PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_interface_view"));
843     PetscCall(DMDestroy(&dmInterface));
844     PetscCall(DMDestroy(dm));
845     *dm = dmHybrid;
846   }
847   PetscCall(DMHasLabel(*dm, "fault2", &hasFault2));
848   if (hasFault2) {
849     DM      dmHybrid = NULL;
850     DMLabel faultLabel, faultBdLabel, hybridLabel;
851 
852     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "faulted_"));
853     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre"));
854     PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
855     PetscCall(DMSetFromOptions(*dm));
856     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
857     PetscCall(DMGetLabel(*dm, "fault2", &faultLabel));
858     PetscCall(DMGetLabel(*dm, "fault2Bd", &faultBdLabel));
859     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, NULL, &dmHybrid));
860     PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD));
861     PetscCall(DMLabelDestroy(&hybridLabel));
862     PetscCall(DMDestroy(dm));
863     *dm = dmHybrid;
864   }
865   if (user->testPartition && size > 1) {
866     PetscPartitioner part;
867     PetscInt        *sizes  = NULL;
868     PetscInt        *points = NULL;
869 
870     if (rank == 0) {
871       if (dim == 2 && cellSimplex && size == 2) {
872         switch (user->testNum) {
873         case 0: {
874           PetscInt triSizes_p2[2]  = {1, 2};
875           PetscInt triPoints_p2[3] = {0, 1, 2};
876 
877           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
878           PetscCall(PetscArraycpy(sizes, triSizes_p2, 2));
879           PetscCall(PetscArraycpy(points, triPoints_p2, 3));
880           break;
881         }
882         case 3: {
883           PetscInt triSizes_p2[2]  = {3, 3};
884           PetscInt triPoints_p2[6] = {0, 1, 2, 3, 4, 5};
885 
886           PetscCall(PetscMalloc2(2, &sizes, 6, &points));
887           PetscCall(PetscArraycpy(sizes, triSizes_p2, 2));
888           PetscCall(PetscArraycpy(points, triPoints_p2, 6));
889           break;
890         }
891         default:
892           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 2 procs", user->testNum);
893         }
894       } else if (dim == 2 && cellSimplex && size == 6) {
895         switch (user->testNum) {
896         case 4: {
897           PetscInt triSizes_p6[6]  = {1, 1, 1, 1, 1, 1};
898           PetscInt triPoints_p6[6] = {0, 1, 2, 3, 4, 5};
899 
900           PetscCall(PetscMalloc2(6, &sizes, 6, &points));
901           PetscCall(PetscArraycpy(sizes, triSizes_p6, 6));
902           PetscCall(PetscArraycpy(points, triPoints_p6, 6));
903           break;
904         }
905         default:
906           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 6 procs", user->testNum);
907         }
908       } else if (dim == 2 && !cellSimplex && size == 2) {
909         switch (user->testNum) {
910         case 0: {
911           PetscInt quadSizes_p2[2]  = {1, 2};
912           PetscInt quadPoints_p2[3] = {0, 1, 2};
913 
914           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
915           PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2));
916           PetscCall(PetscArraycpy(points, quadPoints_p2, 3));
917           break;
918         }
919         case 2: {
920           PetscInt quadSizes_p2[2]  = {1, 1};
921           PetscInt quadPoints_p2[2] = {0, 1};
922 
923           PetscCall(PetscMalloc2(2, &sizes, 2, &points));
924           PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2));
925           PetscCall(PetscArraycpy(points, quadPoints_p2, 2));
926           break;
927         }
928         case 3: {
929           PetscInt quadSizes_p2[2]  = {1, 1};
930           PetscInt quadPoints_p2[2] = {1, 0};
931 
932           PetscCall(PetscMalloc2(2, &sizes, 2, &points));
933           PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2));
934           PetscCall(PetscArraycpy(points, quadPoints_p2, 2));
935           break;
936         }
937         default:
938           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for quadrilateral mesh on 2 procs", user->testNum);
939         }
940       } else if (dim == 3 && cellSimplex && size == 2) {
941         switch (user->testNum) {
942         case 0: {
943           PetscInt tetSizes_p2[2]  = {1, 2};
944           PetscInt tetPoints_p2[3] = {0, 1, 2};
945 
946           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
947           PetscCall(PetscArraycpy(sizes, tetSizes_p2, 2));
948           PetscCall(PetscArraycpy(points, tetPoints_p2, 3));
949           break;
950         }
951         default:
952           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for teterehedral mesh on 2 procs", user->testNum);
953         }
954       } else if (dim == 3 && !cellSimplex && size == 2) {
955         switch (user->testNum) {
956         case 0: {
957           PetscInt hexSizes_p2[2]  = {1, 2};
958           PetscInt hexPoints_p2[3] = {0, 1, 2};
959 
960           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
961           PetscCall(PetscArraycpy(sizes, hexSizes_p2, 2));
962           PetscCall(PetscArraycpy(points, hexPoints_p2, 3));
963           break;
964         }
965         default:
966           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for hexahedral mesh on 2 procs", user->testNum);
967         }
968       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition");
969     }
970     PetscCall(DMPlexGetPartitioner(*dm, &part));
971     PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL));
972     PetscCall(PetscPartitionerShellSetPartition(part, size, sizes, points));
973     PetscCall(PetscFree2(sizes, points));
974   }
975   {
976     DM pdm = NULL;
977 
978     /* Distribute mesh over processes */
979     PetscCall(DMPlexDistribute(*dm, 0, NULL, &pdm));
980     if (pdm) {
981       PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view"));
982       PetscCall(DMDestroy(dm));
983       *dm = pdm;
984     }
985   }
986   PetscCall(DMHasLabel(*dm, "pfault", &hasParallelFault));
987   if (hasParallelFault) {
988     DM      dmHybrid = NULL, dmInterface;
989     DMLabel faultLabel, faultBdLabel, hybridLabel;
990 
991     PetscCall(DMGetLabel(*dm, "pfault", &faultLabel));
992     PetscCall(DMGetLabel(*dm, "pfaultBd", &faultBdLabel));
993     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, &dmInterface, &dmHybrid));
994     PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_fault_view"));
995     {
996       PetscViewer viewer;
997       PetscMPIInt rank;
998 
999       PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)*dm), &rank));
1000       PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer));
1001       PetscCall(PetscViewerASCIIPrintf(viewer, "Rank %d\n", rank));
1002       PetscCall(DMLabelView(hybridLabel, viewer));
1003       PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer));
1004       PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1005     }
1006     PetscCall(DMLabelDestroy(&hybridLabel));
1007     PetscCall(DMDestroy(&dmInterface));
1008     PetscCall(DMDestroy(dm));
1009     *dm = dmHybrid;
1010   }
1011   PetscCall(PetscObjectSetName((PetscObject)*dm, "Hybrid Mesh"));
1012   PetscCall(CreateFaultLabel(*dm));
1013   PetscCall(CreateDiscretization(*dm, user));
1014   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre"));
1015   PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
1016   PetscCall(DMSetFromOptions(*dm));
1017   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1018   PetscFunctionReturn(PETSC_SUCCESS);
1019 }
1020 
1021 static PetscErrorCode TestMesh(DM dm, AppCtx *user)
1022 {
1023   PetscFunctionBegin;
1024   PetscCall(DMPlexCheckSymmetry(dm));
1025   PetscCall(DMPlexCheckSkeleton(dm, 0));
1026   PetscCall(DMPlexCheckFaces(dm, 0));
1027   PetscFunctionReturn(PETSC_SUCCESS);
1028 }
1029 
1030 static PetscErrorCode TestDiscretization(DM dm, AppCtx *user)
1031 {
1032   PetscSection s;
1033 
1034   PetscFunctionBegin;
1035   PetscCall(DMGetSection(dm, &s));
1036   PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-local_section_view"));
1037   PetscFunctionReturn(PETSC_SUCCESS);
1038 }
1039 
1040 static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1041 {
1042   PetscInt d;
1043   for (d = 0; d < dim; ++d) u[d] = x[d];
1044   return PETSC_SUCCESS;
1045 }
1046 
1047 static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1048 {
1049   PetscInt d;
1050   for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0);
1051   return PETSC_SUCCESS;
1052 }
1053 
1054 static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1055 {
1056   PetscInt d;
1057   u[0] = -x[1];
1058   u[1] = x[0];
1059   for (d = 2; d < dim; ++d) u[d] = x[d];
1060   return PETSC_SUCCESS;
1061 }
1062 
1063 /* \lambda \cdot (\psi_u^- - \psi_u^+) */
1064 static void f0_bd_u_neg(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1065 {
1066   const PetscInt Nc = dim + 1;
1067   for (PetscInt c = 0; c < Nc; ++c) f0[c] = -u[uOff[1] + c];
1068 }
1069 
1070 static void f0_bd_u_pos(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1071 {
1072   const PetscInt Nc = dim + 1;
1073   for (PetscInt c = 0; c < Nc; ++c) f0[c] = u[uOff[1] + c];
1074 }
1075 
1076 /* (d - u^+ + u^-) \cdot \psi_\lambda */
1077 static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1078 {
1079   const PetscInt Nc = uOff[2] - uOff[1];
1080 
1081   for (PetscInt c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc + c];
1082 }
1083 
1084 /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */
1085 static void g0_bd_ul_neg(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1086 {
1087   const PetscInt Nc = dim + 1;
1088   for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = -1.0;
1089 }
1090 
1091 static void g0_bd_ul_pos(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1092 {
1093   const PetscInt Nc = dim + 1;
1094   for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = 1.0;
1095 }
1096 
1097 /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */
1098 static void g0_bd_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1099 {
1100   const PetscInt Nc = uOff[2] - uOff[1];
1101 
1102   for (PetscInt c = 0; c < Nc; ++c) {
1103     g0[c * Nc + c]           = -1.0;
1104     g0[Nc * Nc + c * Nc + c] = 1.0;
1105   }
1106 }
1107 
1108 static PetscErrorCode TestAssembly(DM dm, AppCtx *user)
1109 {
1110   Mat           J;
1111   Vec           locX, locF;
1112   PetscDS       probh;
1113   DMLabel       fault, material;
1114   IS            cohesiveCells;
1115   PetscWeakForm wf;
1116   PetscFormKey  keys[3];
1117   PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx);
1118   PetscInt    dim, Nf, cMax, cEnd, id;
1119   PetscMPIInt rank;
1120 
1121   PetscFunctionBegin;
1122   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1123   PetscCall(DMGetDimension(dm, &dim));
1124   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax));
1125   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
1126   PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells));
1127   PetscCall(DMGetLabel(dm, "cohesive", &fault));
1128   PetscCall(DMGetLocalVector(dm, &locX));
1129   PetscCall(PetscObjectSetName((PetscObject)locX, "Local Solution"));
1130   PetscCall(DMGetLocalVector(dm, &locF));
1131   PetscCall(PetscObjectSetName((PetscObject)locF, "Local Residual"));
1132   PetscCall(DMCreateMatrix(dm, &J));
1133   PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
1134 
1135   /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */
1136   PetscCall(DMGetLabel(dm, "material", &material));
1137   id              = 1;
1138   initialGuess[0] = r;
1139   initialGuess[1] = NULL;
1140   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1141   id              = 2;
1142   initialGuess[0] = rp1;
1143   initialGuess[1] = NULL;
1144   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1145   id              = 1;
1146   initialGuess[0] = NULL;
1147   initialGuess[1] = phi;
1148   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1149   PetscCall(VecViewFromOptions(locX, NULL, "-local_solution_view"));
1150 
1151   PetscCall(DMGetCellDS(dm, cMax, &probh, NULL));
1152   PetscCall(PetscDSGetWeakForm(probh, &wf));
1153   PetscCall(PetscDSGetNumFields(probh, &Nf));
1154   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u_neg, 0, NULL));
1155   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u_pos, 0, NULL));
1156   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul_neg, 0, NULL, 0, NULL, 0, NULL));
1157   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul_pos, 0, NULL, 0, NULL, 0, NULL));
1158   if (Nf > 1) {
1159     PetscCall(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL));
1160     PetscCall(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL));
1161   }
1162   if (rank == 0) PetscCall(PetscDSView(probh, NULL));
1163 
1164   keys[0].label = material;
1165   keys[0].value = 1;
1166   keys[0].field = 0;
1167   keys[0].part  = 0;
1168   keys[1].label = material;
1169   keys[1].value = 2;
1170   keys[1].field = 0;
1171   keys[1].part  = 0;
1172   keys[2].label = fault;
1173   keys[2].value = 1;
1174   keys[2].field = 1;
1175   keys[2].part  = 0;
1176   PetscCall(VecSet(locF, 0.));
1177   PetscCall(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user));
1178   PetscCall(VecViewFromOptions(locF, NULL, "-local_residual_view"));
1179   PetscCall(MatZeroEntries(J));
1180   PetscCall(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user));
1181   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1182   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1183   PetscCall(MatViewFromOptions(J, NULL, "-local_jacobian_view"));
1184 
1185   PetscCall(DMRestoreLocalVector(dm, &locX));
1186   PetscCall(DMRestoreLocalVector(dm, &locF));
1187   PetscCall(MatDestroy(&J));
1188   PetscCall(ISDestroy(&cohesiveCells));
1189 
1190   if (cMax < cEnd) {
1191     PetscDS         ds;
1192     PetscFE         fe;
1193     PetscQuadrature quad;
1194     IS             *perm;
1195     const PetscInt *cone;
1196     PetscInt        Na, a;
1197 
1198     PetscCall(DMPlexGetCone(dm, cMax, &cone));
1199     PetscCall(DMGetCellDS(dm, cMax, &ds, NULL));
1200     PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
1201     PetscCall(PetscFEGetQuadrature(fe, &quad));
1202     PetscCall(PetscQuadratureComputePermutations(quad, &Na, &perm));
1203     for (a = 0; a < Na; ++a) PetscCall(ISDestroy(&perm[a]));
1204     PetscCall(PetscFree(perm));
1205   }
1206   PetscFunctionReturn(PETSC_SUCCESS);
1207 }
1208 
1209 int main(int argc, char **argv)
1210 {
1211   DM     dm;
1212   AppCtx user; /* user-defined work context */
1213 
1214   PetscFunctionBeginUser;
1215   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1216   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1217   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1218   PetscCall(TestMesh(dm, &user));
1219   PetscCall(TestDiscretization(dm, &user));
1220   PetscCall(TestAssembly(dm, &user));
1221   PetscCall(DMDestroy(&dm));
1222   PetscCall(PetscFinalize());
1223   return 0;
1224 }
1225 
1226 /*TEST
1227   testset:
1228     args: -orig_dm_plex_check_all -dm_plex_check_all \
1229           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \
1230           -local_solution_view -local_residual_view -local_jacobian_view
1231     test:
1232       suffix: tri_0
1233       args: -dim 2
1234       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1235     test:
1236       suffix: tri_t1_0
1237       args: -dim 2 -test_num 1
1238       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1239     test:
1240       suffix: tri_t2_0
1241       args: -dim 2 -test_num 2
1242       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1243 
1244     test:
1245       suffix: tri_t5_0
1246       args: -dim 2 -test_num 5
1247       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1248 
1249     test:
1250       suffix: tet_0
1251       args: -dim 3
1252       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1253 
1254     test:
1255       suffix: tet_t1_0
1256       args: -dim 3 -test_num 1
1257       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1258 
1259   testset:
1260     args: -orig_dm_plex_check_all -dm_plex_check_all \
1261           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1262     test:
1263       suffix: tet_1
1264       nsize: 2
1265       args: -dim 3
1266       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1267 
1268     test:
1269       suffix: tri_1
1270       nsize: 2
1271       args: -dim 2
1272       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1273 
1274     test:
1275       suffix: tri_t3_0
1276       nsize: 2
1277       args: -dim 2 -test_num 3
1278       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1279 
1280     test:
1281       suffix: tri_t4_0
1282       nsize: 6
1283       args: -dim 2 -test_num 4
1284       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1285 
1286   testset:
1287     args: -orig_dm_plex_check_all -dm_plex_check_all \
1288           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1289     # 2D Quads
1290     test:
1291       suffix: quad_0
1292       args: -dim 2 -cell_simplex 0
1293       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1294 
1295     test:
1296       suffix: quad_1
1297       nsize: 2
1298       args: -dim 2 -cell_simplex 0
1299       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1300 
1301     test:
1302       suffix: quad_t1_0
1303       args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all
1304       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1305 
1306     test:
1307       suffix: quad_t2_0
1308       nsize: 2
1309       args: -dim 2 -cell_simplex 0 -test_num 2
1310       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1311 
1312     test:
1313       # TODO: The PetscSF is wrong here (connects to wrong side of split)
1314       suffix: quad_t3_0
1315       nsize: 2
1316       args: -dim 2 -cell_simplex 0 -test_num 3
1317       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1318 
1319     # 3D Hex
1320     test:
1321       suffix: hex_0
1322       args: -dim 3 -cell_simplex 0
1323       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1324     test:
1325       suffix: hex_1
1326       nsize: 2
1327       args: -dim 3 -cell_simplex 0
1328       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1329     test:
1330       suffix: hex_t1_0
1331       args: -dim 3 -cell_simplex 0 -test_num 1
1332       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1333     test:
1334       suffix: hex_t2_0
1335       args: -dim 3 -cell_simplex 0 -test_num 2
1336       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1337 
1338 TEST*/
1339