xref: /petsc/src/dm/impls/plex/tests/ex5.c (revision 76d901e46dda72c1afe96306c7cb4731c47d4e87)
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(0);
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     {
434       PetscInt    numPoints[2]        = {4, 2};
435       PetscInt    coneSize[6]         = {3, 3, 0, 0, 0, 0};
436       PetscInt    cones[6]            = {2, 3, 4,  5, 4, 3};
437       PetscInt    coneOrientations[6] = {0, 0, 0,  0, 0, 0};
438       PetscScalar vertexCoords[8]     = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
439       PetscInt    markerPoints[8]     = {2, 1, 3, 1, 4, 1, 5, 1};
440       PetscInt    faultPoints[2]      = {3, 4};
441 
442       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
443       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
444       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
445       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
446       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
447     }
448     break;
449     case 1:
450     {
451       PetscInt    numPoints[2]         = {6, 4};
452       PetscInt    coneSize[10]         = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0};
453       PetscInt    cones[12]            = {4, 6, 5,  5, 6, 7,  8, 5, 9,  8, 4, 5};
454       PetscInt    coneOrientations[12] = {0, 0, 0,  0, 0, 0,  0, 0, 0,  0, 0, 0};
455       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};
456       PetscInt    markerPoints[6]      = {4, 1, 6, 1, 8, 1};
457       PetscInt    faultPoints[3]       = {5, 6, 8};
458 
459       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
460       for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
461       for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
462       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));PetscCall(DMSetLabelValue(*dm, "material", 3, 1));
463       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));PetscCall(DMSetLabelValue(*dm, "material", 2, 2));
464     }
465     break;
466     case 2:
467     case 3:
468     case 4:
469     {
470       PetscInt    numPoints[2]         = {8, 6};
471       PetscInt    coneSize[14]         = {3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0};
472       PetscInt    cones[18]            = {6, 7,  9,   9, 12, 11,  7, 12,  9,
473                                           7, 8, 10,  10, 13, 12,  7, 10, 12};
474       PetscInt    coneOrientations[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
475       PetscScalar vertexCoords[16]     = {-1., -1.,  0., -1.,  1., -1.,  -1., 0.,  1., 0.,
476                                           -1.,  1.,  0.,  1.,  1.,  1.,};
477       PetscInt    markerPoints[16]     = {6, 1, 7, 1, 8, 1, 9, 1, 10, 1, 11, 1, 12, 1, 13, 1};
478       PetscInt    faultPoints[2]       = {7, 12};
479 
480       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
481       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
482       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
483       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
484       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
485       PetscCall(DMSetLabelValue(*dm, "material", 4, 2));
486       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
487       for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
488       if (testNum == 2) for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
489       if (testNum == 3 || testNum == 4) for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1));
490     }
491     break;
492     case 5:
493     {
494       PetscInt    numPoints[2]         = {7, 6};
495       PetscInt    coneSize[13]         = {3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0};
496       PetscInt    cones[18]            = {6, 7, 8,  8, 9, 6,  7, 10, 8,  9, 8, 12,  8, 10, 11,  11, 12, 8};
497       PetscInt    coneOrientations[18] = {0, 0, 0,  0, 0, 0,  0,  0, 0,  0, 0,  0,  0,  0,  0,   0,  0, 0};
498       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};
499       PetscInt    faultPoints[2]       = {8, 11};
500       PetscInt    faultBdPoints[1]     = {8};
501 
502       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
503       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
504       for (p = 0; p < 1; ++p) PetscCall(DMSetLabelValue(*dm, "faultBd", faultBdPoints[p], 1));
505       PetscCall(DMSetLabelValue(*dm, "material", 0, 0));PetscCall(DMSetLabelValue(*dm, "material", 2, 0));PetscCall(DMSetLabelValue(*dm, "material", 4, 0));
506       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));PetscCall(DMSetLabelValue(*dm, "material", 3, 2));PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
507     }
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(0);
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     {
537       PetscInt    numPoints[4]         = {5, 7, 9, 2};
538       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};
539       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};
540       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};
541       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};
542       PetscInt    markerPoints[20]     = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1};
543       PetscInt    faultPoints[3]      = {3, 4, 5};
544 
545       PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
546       for (p = 0; p < 10; ++p) {
547         PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
548       }
549       for (p = 0; p < 3; ++p) {
550         PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1));
551       }
552       PetscCall(DMSetLabelValue(dm, "material", 0, 1));
553       PetscCall(DMSetLabelValue(dm, "material", 1, 2));
554     }
555     break;
556     case 1:
557     {
558       PetscInt    numPoints[4]         = {6, 13, 12, 4};
559       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};
560       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};
561       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,  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};
562       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};
563       PetscInt    markerPoints[14]     = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1};
564       PetscInt    faultPoints[4]       = {5, 6, 7, 8};
565 
566       PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
567       for (p = 0; p < 7; ++p) {
568         PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
569       }
570       for (p = 0; p < 4; ++p) {
571         PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1));
572       }
573       PetscCall(DMSetLabelValue(dm, "material", 0, 1));
574       PetscCall(DMSetLabelValue(dm, "material", 1, 2));
575     }
576     break;
577     default:
578       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
579     }
580   } else {
581     PetscInt numPoints[4] = {0, 0, 0, 0};
582 
583     PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL));
584     PetscCall(DMCreateLabel(dm, "fault"));
585   }
586   PetscFunctionReturn(0);
587 }
588 
589 static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
590 {
591   DM             idm;
592   PetscInt       p;
593   PetscMPIInt    rank;
594 
595   PetscFunctionBegin;
596   PetscCallMPI(MPI_Comm_rank(comm, &rank));
597   if (rank == 0) {
598     switch (testNum) {
599     case 0:
600     case 2:
601     case 3:
602     {
603       PetscInt    numPoints[2]        = {6, 2};
604       PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
605       PetscInt    cones[8]            = {2, 3, 4, 5,  3, 6, 7, 4};
606       PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
607       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};
608       PetscInt    markerPoints[12]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1};
609       PetscInt    faultPoints[2]      = {3, 4};
610 
611       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
612       for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
613       if (testNum == 0) for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
614       if (testNum == 2 || testNum == 3) for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1));
615       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
616       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
617     }
618     break;
619     case 1:
620     {
621       PetscInt    numPoints[2]         = {16, 9};
622       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};
623       PetscInt    cones[36]            = {9,  13, 14, 10,
624                                           10, 14, 15, 11,
625                                           11, 15, 16, 12,
626                                           13, 17, 18, 14,
627                                           14, 18, 19, 15,
628                                           15, 19, 20, 16,
629                                           17, 21, 22, 18,
630                                           18, 22, 23, 19,
631                                           19, 23, 24, 20};
632       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};
633       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,
634                                           -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};
635       PetscInt    faultPoints[4]       = {11, 15, 19, 20};
636       PetscInt    fault2Points[2]      = {17, 18};
637 
638       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
639       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "fault",  faultPoints[p], 1));
640       for (p = 3; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "faultBd", faultPoints[p], 1));
641       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault2", fault2Points[p], 1));
642       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
643       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
644       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
645       PetscCall(DMSetLabelValue(*dm, "material", 3, 1));
646       PetscCall(DMSetLabelValue(*dm, "material", 4, 1));
647       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
648       PetscCall(DMSetLabelValue(*dm, "material", 6, 2));
649       PetscCall(DMSetLabelValue(*dm, "material", 7, 2));
650       PetscCall(DMSetLabelValue(*dm, "material", 8, 2));
651     }
652     break;
653     default:
654       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
655     }
656   } else {
657     PetscInt numPoints[3] = {0, 0, 0};
658 
659     PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL));
660     if (testNum == 2 || testNum == 3) PetscCall(DMCreateLabel(*dm, "pfault"));
661     else                              PetscCall(DMCreateLabel(*dm, "fault"));
662   }
663   PetscCall(DMPlexInterpolate(*dm, &idm));
664   PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view"));
665   PetscCall(DMDestroy(dm));
666   *dm  = idm;
667   PetscFunctionReturn(0);
668 }
669 
670 static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
671 {
672   DM             idm;
673   PetscInt       depth = 3, p;
674   PetscMPIInt    rank;
675 
676   PetscFunctionBegin;
677   PetscCallMPI(MPI_Comm_rank(comm, &rank));
678   if (rank == 0) {
679     switch (testNum) {
680     case 0:
681     {
682       PetscInt    numPoints[2]         = {12, 2};
683       PetscInt    coneSize[14]         = {8,8, 0,0,0,0,0,0,0,0,0,0,0,0};
684       PetscInt    cones[16]            = {2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8};
685       PetscInt    coneOrientations[16] = {0,0,0,0,0,0,0,0,  0,0, 0,0 ,0, 0, 0,0};
686       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,
687                                           -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0,
688                                            0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0,  0.5,1.0,1.0};
689       PetscInt    markerPoints[52]     = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1};
690       PetscInt    faultPoints[4]       = {3, 4, 7, 8};
691 
692       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
693       PetscCall(DMPlexInterpolate(*dm, &idm));
694       for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(idm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
695       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
696       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
697       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
698     }
699     break;
700     case 1:
701     {
702       /* Cell Adjacency Graph:
703         0 -- { 8, 13, 21, 24} --> 1
704         0 -- {20, 21, 23, 24} --> 5 F
705         1 -- {10, 15, 21, 24} --> 2
706         1 -- {13, 14, 15, 24} --> 6
707         2 -- {21, 22, 24, 25} --> 4 F
708         3 -- {21, 24, 30, 35} --> 4
709         3 -- {21, 24, 28, 33} --> 5
710        */
711       PetscInt    numPoints[2]         = {30, 7};
712       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};
713       PetscInt    cones[56]            = { 8, 21, 20,  7, 13, 12, 23, 24,
714                                           14, 15, 10,  9, 13,  8, 21, 24,
715                                           15, 16, 11, 10, 24, 21, 22, 25,
716                                           30, 29, 28, 21, 35, 24, 33, 34,
717                                           24, 21, 30, 35, 25, 36, 31, 22,
718                                           27, 20, 21, 28, 32, 33, 24, 23,
719                                           15, 24, 13, 14, 19, 18, 17, 26};
720       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};
721       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,
722                                           -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,
723                                           -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,
724                                            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,
725                                            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};
726       PetscInt    faultPoints[6]       = {20, 21, 22, 23, 24, 25};
727 
728       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
729       PetscCall(DMPlexInterpolate(*dm, &idm));
730       for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
731       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
732       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
733       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
734       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
735       PetscCall(DMSetLabelValue(*dm, "material", 4, 2));
736       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
737       PetscCall(DMSetLabelValue(*dm, "material", 6, 2));
738     }
739     break;
740     case 2:
741     {
742       /* Buried fault edge */
743       PetscInt    numPoints[2]         = {18, 4};
744       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};
745       PetscInt    cones[32]            = { 4,  5,  8,  7, 13, 16, 17, 14,
746                                            5,  6,  9,  8, 14, 17, 18, 15,
747                                            7,  8, 11, 10, 16, 19, 20, 17,
748                                            8,  9, 12, 11, 17, 20, 21, 18};
749       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};
750       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,
751                                            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,
752                                            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};
753       PetscInt    faultPoints[4]       = {7, 8, 16, 17};
754 
755       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
756       PetscCall(DMPlexInterpolate(*dm, &idm));
757       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
758       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
759       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
760       PetscCall(DMSetLabelValue(*dm, "material", 2, 2));
761       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
762     }
763     break;
764     default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
765     }
766   } else {
767     PetscInt numPoints[4] = {0, 0, 0, 0};
768 
769     PetscCall(DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL));
770     PetscCall(DMPlexInterpolate(*dm, &idm));
771     PetscCall(DMCreateLabel(idm, "fault"));
772   }
773   PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view"));
774   PetscCall(DMDestroy(dm));
775   *dm  = idm;
776   PetscFunctionReturn(0);
777 }
778 
779 static PetscErrorCode CreateFaultLabel(DM dm)
780 {
781   DMLabel        label;
782   PetscInt       dim, h, pStart, pEnd, pMax, p;
783 
784   PetscFunctionBegin;
785   PetscCall(DMGetDimension(dm, &dim));
786   PetscCall(DMCreateLabel(dm, "cohesive"));
787   PetscCall(DMGetLabel(dm, "cohesive", &label));
788   for (h = 0; h <= dim; ++h) {
789     PetscCall(DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax));
790     PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd));
791     for (p = pMax; p < pEnd; ++p) PetscCall(DMLabelSetValue(label, p, 1));
792   }
793   PetscFunctionReturn(0);
794 }
795 
796 static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user)
797 {
798   PetscFE        fe;
799   DMLabel        fault;
800   PetscInt       dim, Ncf = user->cohesiveFields, f;
801 
802   PetscFunctionBegin;
803   PetscCall(DMGetDimension(dm, &dim));
804   PetscCall(DMGetLabel(dm, "cohesive", &fault));
805   PetscCall(DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD));
806 
807   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe));
808   PetscCall(PetscFESetName(fe, "displacement"));
809   PetscCall(DMAddField(dm, NULL, (PetscObject) fe));
810   PetscCall(PetscFEDestroy(&fe));
811 
812   if (Ncf > 0) {
813     PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe));
814     PetscCall(PetscFESetName(fe, "fault traction"));
815     PetscCall(DMAddField(dm, fault, (PetscObject) fe));
816     PetscCall(PetscFEDestroy(&fe));
817   }
818   for (f = 1; f < Ncf; ++f) {
819     char name[256], opt[256];
820 
821     PetscCall(PetscSNPrintf(name, 256, "fault field %" PetscInt_FMT, f));
822     PetscCall(PetscSNPrintf(opt,  256, "faultfield_%" PetscInt_FMT "_", f));
823     PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe));
824     PetscCall(PetscFESetName(fe, name));
825     PetscCall(DMAddField(dm, fault, (PetscObject) fe));
826     PetscCall(PetscFEDestroy(&fe));
827   }
828 
829   PetscCall(DMCreateDS(dm));
830   PetscFunctionReturn(0);
831 }
832 
833 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
834 {
835   PetscInt       dim          = user->dim;
836   PetscBool      cellSimplex  = user->cellSimplex, hasFault, hasFault2, hasParallelFault;
837   PetscMPIInt    rank, size;
838   DMLabel        matLabel;
839 
840   PetscFunctionBegin;
841   PetscCallMPI(MPI_Comm_rank(comm, &rank));
842   PetscCallMPI(MPI_Comm_size(comm, &size));
843   PetscCall(DMCreate(comm, dm));
844   PetscCall(DMSetType(*dm, DMPLEX));
845   PetscCall(DMSetDimension(*dm, dim));
846   switch (dim) {
847   case 2:
848     if (cellSimplex) {
849       PetscCall(CreateSimplex_2D(comm, user->testNum, dm));
850     } else {
851       PetscCall(CreateQuad_2D(comm, user->testNum, dm));
852     }
853     break;
854   case 3:
855     if (cellSimplex) {
856       PetscCall(CreateSimplex_3D(comm, user, *dm));
857     } else {
858       PetscCall(CreateHex_3D(comm, user->testNum, dm));
859     }
860     break;
861   default:
862     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %" PetscInt_FMT, dim);
863   }
864   PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_"));
865   PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
866   PetscCall(DMSetFromOptions(*dm));
867   PetscCall(DMGetLabel(*dm, "material", &matLabel));
868   if (matLabel) PetscCall(DMPlexLabelComplete(*dm, matLabel));
869   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
870   PetscCall(DMHasLabel(*dm, "fault", &hasFault));
871   if (hasFault) {
872     DM      dmHybrid = NULL, dmInterface = NULL;
873     DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel;
874 
875     PetscCall(DMGetLabel(*dm, "fault", &faultLabel));
876     PetscCall(DMGetLabel(*dm, "faultBd", &faultBdLabel));
877     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid));
878     PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD));
879     PetscCall(DMLabelDestroy(&hybridLabel));
880     PetscCall(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD));
881     PetscCall(DMLabelDestroy(&splitLabel));
882     PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_interface_view"));
883     PetscCall(DMDestroy(&dmInterface));
884     PetscCall(DMDestroy(dm));
885     *dm  = dmHybrid;
886   }
887   PetscCall(DMHasLabel(*dm, "fault2", &hasFault2));
888   if (hasFault2) {
889     DM      dmHybrid = NULL;
890     DMLabel faultLabel, faultBdLabel, hybridLabel;
891 
892     PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_"));
893     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre"));
894     PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
895     PetscCall(DMSetFromOptions(*dm));
896     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
897     PetscCall(DMGetLabel(*dm, "fault2", &faultLabel));
898     PetscCall(DMGetLabel(*dm, "fault2Bd", &faultBdLabel));
899     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, NULL, &dmHybrid));
900     PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD));
901     PetscCall(DMLabelDestroy(&hybridLabel));
902     PetscCall(DMDestroy(dm));
903     *dm  = dmHybrid;
904   }
905   if (user->testPartition && size > 1) {
906     PetscPartitioner part;
907     PetscInt *sizes  = NULL;
908     PetscInt *points = NULL;
909 
910     if (rank == 0) {
911       if (dim == 2 && cellSimplex && size == 2) {
912         switch (user->testNum) {
913         case 0: {
914           PetscInt triSizes_p2[2]  = {1, 2};
915           PetscInt triPoints_p2[3] = {0, 1, 2};
916 
917           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
918           PetscCall(PetscArraycpy(sizes,  triSizes_p2, 2));
919           PetscCall(PetscArraycpy(points, triPoints_p2, 3));break;}
920         case 3: {
921           PetscInt triSizes_p2[2]  = {3, 3};
922           PetscInt triPoints_p2[6] = {0, 1, 2,  3, 4, 5};
923 
924           PetscCall(PetscMalloc2(2, &sizes, 6, &points));
925           PetscCall(PetscArraycpy(sizes,  triSizes_p2, 2));
926           PetscCall(PetscArraycpy(points, triPoints_p2, 6));break;}
927         default:
928           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 2 procs", user->testNum);
929         }
930       } else if (dim == 2 && cellSimplex && size == 6) {
931         switch (user->testNum) {
932         case 4: {
933           PetscInt triSizes_p6[6]  = {1, 1, 1, 1, 1, 1};
934           PetscInt triPoints_p6[6] = {0, 1, 2, 3, 4, 5};
935 
936           PetscCall(PetscMalloc2(6, &sizes, 6, &points));
937           PetscCall(PetscArraycpy(sizes,  triSizes_p6, 6));
938           PetscCall(PetscArraycpy(points, triPoints_p6, 6));break;}
939         default:
940           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 6 procs", user->testNum);
941         }
942       } else if (dim == 2 && !cellSimplex && size == 2) {
943         switch (user->testNum) {
944         case 0: {
945           PetscInt quadSizes_p2[2]  = {1, 2};
946           PetscInt quadPoints_p2[3] = {0, 1, 2};
947 
948           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
949           PetscCall(PetscArraycpy(sizes,  quadSizes_p2, 2));
950           PetscCall(PetscArraycpy(points, quadPoints_p2, 3));break;}
951         case 2: {
952           PetscInt quadSizes_p2[2]  = {1, 1};
953           PetscInt quadPoints_p2[2] = {0, 1};
954 
955           PetscCall(PetscMalloc2(2, &sizes, 2, &points));
956           PetscCall(PetscArraycpy(sizes,  quadSizes_p2, 2));
957           PetscCall(PetscArraycpy(points, quadPoints_p2, 2));break;}
958         case 3: {
959           PetscInt quadSizes_p2[2]  = {1, 1};
960           PetscInt quadPoints_p2[2] = {1, 0};
961 
962           PetscCall(PetscMalloc2(2, &sizes, 2, &points));
963           PetscCall(PetscArraycpy(sizes,  quadSizes_p2, 2));
964           PetscCall(PetscArraycpy(points, quadPoints_p2, 2));break;}
965         default:
966           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for quadrilateral mesh on 2 procs", user->testNum);
967         }
968       } else if (dim == 3 && cellSimplex && size == 2) {
969         switch (user->testNum) {
970         case 0: {
971           PetscInt tetSizes_p2[2]  = {1, 2};
972           PetscInt tetPoints_p2[3] = {0, 1, 2};
973 
974           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
975           PetscCall(PetscArraycpy(sizes,  tetSizes_p2, 2));
976           PetscCall(PetscArraycpy(points, tetPoints_p2, 3));break;}
977         default:
978           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for teterehedral mesh on 2 procs", user->testNum);
979         }
980       } else if (dim == 3 && !cellSimplex && size == 2) {
981         switch (user->testNum) {
982         case 0: {
983           PetscInt hexSizes_p2[2]  = {1, 2};
984           PetscInt hexPoints_p2[3] = {0, 1, 2};
985 
986           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
987           PetscCall(PetscArraycpy(sizes,  hexSizes_p2, 2));
988           PetscCall(PetscArraycpy(points, hexPoints_p2, 3));break;}
989         default:
990           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for hexahedral mesh on 2 procs", user->testNum);
991         }
992       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition");
993     }
994     PetscCall(DMPlexGetPartitioner(*dm, &part));
995     PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL));
996     PetscCall(PetscPartitionerShellSetPartition(part, size, sizes, points));
997     PetscCall(PetscFree2(sizes, points));
998   }
999   {
1000     DM pdm = NULL;
1001 
1002     /* Distribute mesh over processes */
1003     PetscCall(DMPlexDistribute(*dm, 0, NULL, &pdm));
1004     if (pdm) {
1005       PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view"));
1006       PetscCall(DMDestroy(dm));
1007       *dm  = pdm;
1008     }
1009   }
1010   PetscCall(DMHasLabel(*dm, "pfault", &hasParallelFault));
1011   if (hasParallelFault) {
1012     DM      dmHybrid = NULL, dmInterface;
1013     DMLabel faultLabel, faultBdLabel, hybridLabel;
1014 
1015     PetscCall(DMGetLabel(*dm, "pfault", &faultLabel));
1016     PetscCall(DMGetLabel(*dm, "pfaultBd", &faultBdLabel));
1017     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, &dmInterface, &dmHybrid));
1018     PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_fault_view"));
1019     {
1020       PetscViewer viewer;
1021       PetscMPIInt rank;
1022 
1023       PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank));
1024       PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer));
1025       PetscCall(PetscViewerASCIIPrintf(viewer, "Rank %d\n", rank));
1026       PetscCall(DMLabelView(hybridLabel, viewer));
1027       PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer));
1028       PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1029     }
1030     PetscCall(DMLabelDestroy(&hybridLabel));
1031     PetscCall(DMDestroy(&dmInterface));
1032     PetscCall(DMDestroy(dm));
1033     *dm  = dmHybrid;
1034   }
1035   PetscCall(PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh"));
1036   PetscCall(CreateFaultLabel(*dm));
1037   PetscCall(CreateDiscretization(*dm, user));
1038   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre"));
1039   PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
1040   PetscCall(DMSetFromOptions(*dm));
1041   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1042   PetscFunctionReturn(0);
1043 }
1044 
1045 static PetscErrorCode TestMesh(DM dm, AppCtx *user)
1046 {
1047   PetscFunctionBegin;
1048   PetscCall(DMPlexCheckSymmetry(dm));
1049   PetscCall(DMPlexCheckSkeleton(dm, 0));
1050   PetscCall(DMPlexCheckFaces(dm, 0));
1051   PetscFunctionReturn(0);
1052 }
1053 
1054 static PetscErrorCode TestDiscretization(DM dm, AppCtx *user)
1055 {
1056   PetscSection   s;
1057 
1058   PetscFunctionBegin;
1059   PetscCall(DMGetSection(dm, &s));
1060   PetscCall(PetscObjectViewFromOptions((PetscObject) s, NULL, "-local_section_view"));
1061   PetscFunctionReturn(0);
1062 }
1063 
1064 static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1065 {
1066   PetscInt d;
1067   for (d = 0; d < dim; ++d) u[d] = x[d];
1068   return 0;
1069 }
1070 
1071 static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1072 {
1073   PetscInt d;
1074   for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0);
1075   return 0;
1076 }
1077 
1078 static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1079 {
1080   PetscInt d;
1081   u[0] = -x[1];
1082   u[1] =  x[0];
1083   for (d = 2; d < dim; ++d) u[d] = x[d];
1084   return 0;
1085 }
1086 
1087 /* \lambda \cdot (\psi_u^- - \psi_u^+) */
1088 static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1089                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1090                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1091                     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1092 {
1093   const PetscInt Nc = uOff[1]-uOff[0];
1094   PetscInt       c;
1095   for (c = 0;  c < Nc;   ++c) {
1096     f0[c]    =   u[Nc*2+c] + x[Nc-c-1];
1097     f0[Nc+c] = -(u[Nc*2+c] + x[Nc-c-1]);
1098   }
1099 }
1100 
1101 /* (d - u^+ + u^-) \cdot \psi_\lambda */
1102 static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1103                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1104                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1105                     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1106 {
1107   const PetscInt Nc = uOff[2]-uOff[1];
1108   PetscInt       c;
1109 
1110   for (c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc+c];
1111 }
1112 
1113 /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */
1114 static void g0_bd_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1115                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1116                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1117                      PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1118 {
1119   const PetscInt Nc = uOff[1]-uOff[0];
1120   PetscInt       c;
1121 
1122   for (c = 0; c < Nc; ++c) {
1123     g0[(0 +c)*Nc+c] =  1.0;
1124     g0[(Nc+c)*Nc+c] = -1.0;
1125   }
1126 }
1127 
1128 /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */
1129 static void g0_bd_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1130                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1131                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1132                      PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1133 {
1134   const PetscInt Nc = uOff[2]-uOff[1];
1135   PetscInt       c;
1136 
1137   for (c = 0; c < Nc; ++c) {
1138     g0[c*Nc*2+c]    =  1.0;
1139     g0[c*Nc*2+Nc+c] = -1.0;
1140   }
1141 }
1142 
1143 static PetscErrorCode TestAssembly(DM dm, AppCtx *user)
1144 {
1145   Mat              J;
1146   Vec              locX, locF;
1147   PetscDS          probh;
1148   DMLabel          fault, material;
1149   IS               cohesiveCells;
1150   PetscWeakForm    wf;
1151   PetscFormKey     keys[3];
1152   PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx);
1153   PetscInt         dim, Nf, cMax, cEnd, id;
1154   PetscMPIInt      rank;
1155 
1156   PetscFunctionBegin;
1157   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank));
1158   PetscCall(DMGetDimension(dm, &dim));
1159   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax));
1160   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
1161   PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells));
1162   PetscCall(DMGetLabel(dm, "cohesive", &fault));
1163   PetscCall(DMGetLocalVector(dm, &locX));
1164   PetscCall(PetscObjectSetName((PetscObject) locX, "Local Solution"));
1165   PetscCall(DMGetLocalVector(dm, &locF));
1166   PetscCall(PetscObjectSetName((PetscObject) locF, "Local Residual"));
1167   PetscCall(DMCreateMatrix(dm, &J));
1168   PetscCall(PetscObjectSetName((PetscObject) J, "Jacobian"));
1169 
1170   /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */
1171   PetscCall(DMGetLabel(dm, "material", &material));
1172   id   = 1;
1173   initialGuess[0] = r;
1174   initialGuess[1] = NULL;
1175   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1176   id   = 2;
1177   initialGuess[0] = rp1;
1178   initialGuess[1] = NULL;
1179   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1180   id   = 1;
1181   initialGuess[0] = NULL;
1182   initialGuess[1] = phi;
1183   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1184   PetscCall(VecViewFromOptions(locX, NULL, "-local_solution_view"));
1185 
1186   PetscCall(DMGetCellDS(dm, cMax, &probh));
1187   PetscCall(PetscDSGetWeakForm(probh, &wf));
1188   PetscCall(PetscDSGetNumFields(probh, &Nf));
1189   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u, 0, NULL));
1190   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u, 0, NULL));
1191   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL));
1192   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL));
1193   if (Nf > 1) {
1194     PetscCall(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL));
1195     PetscCall(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL));
1196   }
1197   if (rank == 0) PetscCall(PetscDSView(probh, NULL));
1198 
1199   keys[0].label = NULL;
1200   keys[0].value = 0;
1201   keys[0].field = 0;
1202   keys[0].part  = 0;
1203   keys[1].label = material;
1204   keys[1].value = 2;
1205   keys[1].field = 0;
1206   keys[1].part  = 0;
1207   keys[2].label = fault;
1208   keys[2].value = 1;
1209   keys[2].field = 1;
1210   keys[2].part  = 0;
1211   PetscCall(VecSet(locF, 0.));
1212   PetscCall(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user));
1213   PetscCall(VecViewFromOptions(locF, NULL, "-local_residual_view"));
1214   PetscCall(MatZeroEntries(J));
1215   PetscCall(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user));
1216   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1217   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1218   PetscCall(MatViewFromOptions(J, NULL, "-local_jacobian_view"));
1219 
1220   PetscCall(DMRestoreLocalVector(dm, &locX));
1221   PetscCall(DMRestoreLocalVector(dm, &locF));
1222   PetscCall(MatDestroy(&J));
1223   PetscCall(ISDestroy(&cohesiveCells));
1224   PetscFunctionReturn(0);
1225 }
1226 
1227 int main(int argc, char **argv)
1228 {
1229   DM             dm;
1230   AppCtx         user;                 /* user-defined work context */
1231 
1232   PetscFunctionBeginUser;
1233   PetscCall(PetscInitialize(&argc, &argv, NULL,help));
1234   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1235   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1236   PetscCall(TestMesh(dm, &user));
1237   PetscCall(TestDiscretization(dm, &user));
1238   PetscCall(TestAssembly(dm, &user));
1239   PetscCall(DMDestroy(&dm));
1240   PetscCall(PetscFinalize());
1241   return 0;
1242 }
1243 
1244 /*TEST
1245   testset:
1246     args: -orig_dm_plex_check_all -dm_plex_check_all \
1247           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \
1248           -local_solution_view -local_residual_view -local_jacobian_view
1249     test:
1250       suffix: tri_0
1251       args: -dim 2
1252       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1253     test:
1254       suffix: tri_t1_0
1255       args: -dim 2 -test_num 1
1256       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1257     test:
1258       suffix: tri_t2_0
1259       args: -dim 2 -test_num 2
1260       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1261     test:
1262       suffix: tri_t5_0
1263       args: -dim 2 -test_num 5
1264       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1265     test:
1266       suffix: tet_0
1267       args: -dim 3
1268       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1269     test:
1270       suffix: tet_t1_0
1271       args: -dim 3 -test_num 1
1272       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1273 
1274   testset:
1275     args: -orig_dm_plex_check_all -dm_plex_check_all \
1276           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1277     test:
1278       suffix: tet_1
1279       nsize: 2
1280       args: -dim 3
1281       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1282     test:
1283       suffix: tri_1
1284       nsize: 2
1285       args: -dim 2
1286       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1287     test:
1288       suffix: tri_t3_0
1289       nsize: 2
1290       args: -dim 2 -test_num 3
1291       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1292     test:
1293       suffix: tri_t4_0
1294       nsize: 6
1295       args: -dim 2 -test_num 4
1296       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1297 
1298   testset:
1299     args: -orig_dm_plex_check_all -dm_plex_check_all \
1300           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1301     # 2D Quads
1302     test:
1303       suffix: quad_0
1304       args: -dim 2 -cell_simplex 0
1305       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1306     test:
1307       suffix: quad_1
1308       nsize: 2
1309       args: -dim 2 -cell_simplex 0
1310       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1311     test:
1312       suffix: quad_t1_0
1313       args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all
1314       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1315     test:
1316       suffix: quad_t2_0
1317       nsize: 2
1318       args: -dim 2 -cell_simplex 0 -test_num 2
1319       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1320     test:
1321       # TODO: The PetscSF is wrong here (connects to wrong side of split)
1322       suffix: quad_t3_0
1323       nsize: 2
1324       args: -dim 2 -cell_simplex 0 -test_num 3
1325       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1326     # 3D Hex
1327     test:
1328       suffix: hex_0
1329       args: -dim 3 -cell_simplex 0
1330       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1331     test:
1332       suffix: hex_1
1333       nsize: 2
1334       args: -dim 3 -cell_simplex 0
1335       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1336     test:
1337       suffix: hex_t1_0
1338       args: -dim 3 -cell_simplex 0 -test_num 1
1339       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1340     test:
1341       suffix: hex_t2_0
1342       args: -dim 3 -cell_simplex 0 -test_num 2
1343       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1344 
1345 TEST*/
1346