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