xref: /petsc/src/dm/impls/plex/tests/ex5.c (revision f97672e55eacc8688507b9471cd7ec2664d7f203)
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) {
869     PetscCall(DMPlexLabelComplete(*dm, matLabel));
870   }
871   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
872   PetscCall(DMHasLabel(*dm, "fault", &hasFault));
873   if (hasFault) {
874     DM      dmHybrid = NULL, dmInterface = NULL;
875     DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel;
876 
877     PetscCall(DMGetLabel(*dm, "fault", &faultLabel));
878     PetscCall(DMGetLabel(*dm, "faultBd", &faultBdLabel));
879     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid));
880     PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD));
881     PetscCall(DMLabelDestroy(&hybridLabel));
882     PetscCall(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD));
883     PetscCall(DMLabelDestroy(&splitLabel));
884     PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_interface_view"));
885     PetscCall(DMDestroy(&dmInterface));
886     PetscCall(DMDestroy(dm));
887     *dm  = dmHybrid;
888   }
889   PetscCall(DMHasLabel(*dm, "fault2", &hasFault2));
890   if (hasFault2) {
891     DM      dmHybrid = NULL;
892     DMLabel faultLabel, faultBdLabel, hybridLabel;
893 
894     PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_"));
895     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre"));
896     PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
897     PetscCall(DMSetFromOptions(*dm));
898     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
899     PetscCall(DMGetLabel(*dm, "fault2", &faultLabel));
900     PetscCall(DMGetLabel(*dm, "fault2Bd", &faultBdLabel));
901     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, NULL, &dmHybrid));
902     PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD));
903     PetscCall(DMLabelDestroy(&hybridLabel));
904     PetscCall(DMDestroy(dm));
905     *dm  = dmHybrid;
906   }
907   if (user->testPartition && size > 1) {
908     PetscPartitioner part;
909     PetscInt *sizes  = NULL;
910     PetscInt *points = NULL;
911 
912     if (rank == 0) {
913       if (dim == 2 && cellSimplex && size == 2) {
914         switch (user->testNum) {
915         case 0: {
916           PetscInt triSizes_p2[2]  = {1, 2};
917           PetscInt triPoints_p2[3] = {0, 1, 2};
918 
919           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
920           PetscCall(PetscArraycpy(sizes,  triSizes_p2, 2));
921           PetscCall(PetscArraycpy(points, triPoints_p2, 3));break;}
922         case 3: {
923           PetscInt triSizes_p2[2]  = {3, 3};
924           PetscInt triPoints_p2[6] = {0, 1, 2,  3, 4, 5};
925 
926           PetscCall(PetscMalloc2(2, &sizes, 6, &points));
927           PetscCall(PetscArraycpy(sizes,  triSizes_p2, 2));
928           PetscCall(PetscArraycpy(points, triPoints_p2, 6));break;}
929         default:
930           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 2 procs", user->testNum);
931         }
932       } else if (dim == 2 && cellSimplex && size == 6) {
933         switch (user->testNum) {
934         case 4: {
935           PetscInt triSizes_p6[6]  = {1, 1, 1, 1, 1, 1};
936           PetscInt triPoints_p6[6] = {0, 1, 2, 3, 4, 5};
937 
938           PetscCall(PetscMalloc2(6, &sizes, 6, &points));
939           PetscCall(PetscArraycpy(sizes,  triSizes_p6, 6));
940           PetscCall(PetscArraycpy(points, triPoints_p6, 6));break;}
941         default:
942           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 6 procs", user->testNum);
943         }
944       } else if (dim == 2 && !cellSimplex && size == 2) {
945         switch (user->testNum) {
946         case 0: {
947           PetscInt quadSizes_p2[2]  = {1, 2};
948           PetscInt quadPoints_p2[3] = {0, 1, 2};
949 
950           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
951           PetscCall(PetscArraycpy(sizes,  quadSizes_p2, 2));
952           PetscCall(PetscArraycpy(points, quadPoints_p2, 3));break;}
953         case 2: {
954           PetscInt quadSizes_p2[2]  = {1, 1};
955           PetscInt quadPoints_p2[2] = {0, 1};
956 
957           PetscCall(PetscMalloc2(2, &sizes, 2, &points));
958           PetscCall(PetscArraycpy(sizes,  quadSizes_p2, 2));
959           PetscCall(PetscArraycpy(points, quadPoints_p2, 2));break;}
960         case 3: {
961           PetscInt quadSizes_p2[2]  = {1, 1};
962           PetscInt quadPoints_p2[2] = {1, 0};
963 
964           PetscCall(PetscMalloc2(2, &sizes, 2, &points));
965           PetscCall(PetscArraycpy(sizes,  quadSizes_p2, 2));
966           PetscCall(PetscArraycpy(points, quadPoints_p2, 2));break;}
967         default:
968           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for quadrilateral mesh on 2 procs", user->testNum);
969         }
970       } else if (dim == 3 && cellSimplex && size == 2) {
971         switch (user->testNum) {
972         case 0: {
973           PetscInt tetSizes_p2[2]  = {1, 2};
974           PetscInt tetPoints_p2[3] = {0, 1, 2};
975 
976           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
977           PetscCall(PetscArraycpy(sizes,  tetSizes_p2, 2));
978           PetscCall(PetscArraycpy(points, tetPoints_p2, 3));break;}
979         default:
980           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for teterehedral mesh on 2 procs", user->testNum);
981         }
982       } else if (dim == 3 && !cellSimplex && size == 2) {
983         switch (user->testNum) {
984         case 0: {
985           PetscInt hexSizes_p2[2]  = {1, 2};
986           PetscInt hexPoints_p2[3] = {0, 1, 2};
987 
988           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
989           PetscCall(PetscArraycpy(sizes,  hexSizes_p2, 2));
990           PetscCall(PetscArraycpy(points, hexPoints_p2, 3));break;}
991         default:
992           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for hexahedral mesh on 2 procs", user->testNum);
993         }
994       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition");
995     }
996     PetscCall(DMPlexGetPartitioner(*dm, &part));
997     PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL));
998     PetscCall(PetscPartitionerShellSetPartition(part, size, sizes, points));
999     PetscCall(PetscFree2(sizes, points));
1000   }
1001   {
1002     DM pdm = NULL;
1003 
1004     /* Distribute mesh over processes */
1005     PetscCall(DMPlexDistribute(*dm, 0, NULL, &pdm));
1006     if (pdm) {
1007       PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view"));
1008       PetscCall(DMDestroy(dm));
1009       *dm  = pdm;
1010     }
1011   }
1012   PetscCall(DMHasLabel(*dm, "pfault", &hasParallelFault));
1013   if (hasParallelFault) {
1014     DM      dmHybrid = NULL, dmInterface;
1015     DMLabel faultLabel, faultBdLabel, hybridLabel;
1016 
1017     PetscCall(DMGetLabel(*dm, "pfault", &faultLabel));
1018     PetscCall(DMGetLabel(*dm, "pfaultBd", &faultBdLabel));
1019     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, &dmInterface, &dmHybrid));
1020     PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_fault_view"));
1021     {
1022       PetscViewer viewer;
1023       PetscMPIInt rank;
1024 
1025       PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank));
1026       PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer));
1027       PetscCall(PetscViewerASCIIPrintf(viewer, "Rank %d\n", rank));
1028       PetscCall(DMLabelView(hybridLabel, viewer));
1029       PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer));
1030       PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1031     }
1032     PetscCall(DMLabelDestroy(&hybridLabel));
1033     PetscCall(DMDestroy(&dmInterface));
1034     PetscCall(DMDestroy(dm));
1035     *dm  = dmHybrid;
1036   }
1037   PetscCall(PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh"));
1038   PetscCall(CreateFaultLabel(*dm));
1039   PetscCall(CreateDiscretization(*dm, user));
1040   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre"));
1041   PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
1042   PetscCall(DMSetFromOptions(*dm));
1043   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1044   PetscFunctionReturn(0);
1045 }
1046 
1047 static PetscErrorCode TestMesh(DM dm, AppCtx *user)
1048 {
1049   PetscFunctionBegin;
1050   PetscCall(DMPlexCheckSymmetry(dm));
1051   PetscCall(DMPlexCheckSkeleton(dm, 0));
1052   PetscCall(DMPlexCheckFaces(dm, 0));
1053   PetscFunctionReturn(0);
1054 }
1055 
1056 static PetscErrorCode TestDiscretization(DM dm, AppCtx *user)
1057 {
1058   PetscSection   s;
1059 
1060   PetscFunctionBegin;
1061   PetscCall(DMGetSection(dm, &s));
1062   PetscCall(PetscObjectViewFromOptions((PetscObject) s, NULL, "-local_section_view"));
1063   PetscFunctionReturn(0);
1064 }
1065 
1066 static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1067 {
1068   PetscInt d;
1069   for (d = 0; d < dim; ++d) u[d] = x[d];
1070   return 0;
1071 }
1072 
1073 static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1074 {
1075   PetscInt d;
1076   for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0);
1077   return 0;
1078 }
1079 
1080 static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1081 {
1082   PetscInt d;
1083   u[0] = -x[1];
1084   u[1] =  x[0];
1085   for (d = 2; d < dim; ++d) u[d] = x[d];
1086   return 0;
1087 }
1088 
1089 /* \lambda \cdot (\psi_u^- - \psi_u^+) */
1090 static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1091                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1092                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1093                     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1094 {
1095   const PetscInt Nc = uOff[1]-uOff[0];
1096   PetscInt       c;
1097   for (c = 0;  c < Nc;   ++c) {
1098     f0[c]    =   u[Nc*2+c] + x[Nc-c-1];
1099     f0[Nc+c] = -(u[Nc*2+c] + x[Nc-c-1]);
1100   }
1101 }
1102 
1103 /* (d - u^+ + u^-) \cdot \psi_\lambda */
1104 static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1105                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1106                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1107                     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1108 {
1109   const PetscInt Nc = uOff[2]-uOff[1];
1110   PetscInt       c;
1111 
1112   for (c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc+c];
1113 }
1114 
1115 /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */
1116 static void g0_bd_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1117                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1118                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1119                      PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1120 {
1121   const PetscInt Nc = uOff[1]-uOff[0];
1122   PetscInt       c;
1123 
1124   for (c = 0; c < Nc; ++c) {
1125     g0[(0 +c)*Nc+c] =  1.0;
1126     g0[(Nc+c)*Nc+c] = -1.0;
1127   }
1128 }
1129 
1130 /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */
1131 static void g0_bd_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1132                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1133                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1134                      PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1135 {
1136   const PetscInt Nc = uOff[2]-uOff[1];
1137   PetscInt       c;
1138 
1139   for (c = 0; c < Nc; ++c) {
1140     g0[c*Nc*2+c]    =  1.0;
1141     g0[c*Nc*2+Nc+c] = -1.0;
1142   }
1143 }
1144 
1145 static PetscErrorCode TestAssembly(DM dm, AppCtx *user)
1146 {
1147   Mat              J;
1148   Vec              locX, locF;
1149   PetscDS          probh;
1150   DMLabel          fault, material;
1151   IS               cohesiveCells;
1152   PetscWeakForm    wf;
1153   PetscFormKey     keys[3];
1154   PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx);
1155   PetscInt         dim, Nf, cMax, cEnd, id;
1156   PetscMPIInt      rank;
1157 
1158   PetscFunctionBegin;
1159   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank));
1160   PetscCall(DMGetDimension(dm, &dim));
1161   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax));
1162   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
1163   PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells));
1164   PetscCall(DMGetLabel(dm, "cohesive", &fault));
1165   PetscCall(DMGetLocalVector(dm, &locX));
1166   PetscCall(PetscObjectSetName((PetscObject) locX, "Local Solution"));
1167   PetscCall(DMGetLocalVector(dm, &locF));
1168   PetscCall(PetscObjectSetName((PetscObject) locF, "Local Residual"));
1169   PetscCall(DMCreateMatrix(dm, &J));
1170   PetscCall(PetscObjectSetName((PetscObject) J, "Jacobian"));
1171 
1172   /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */
1173   PetscCall(DMGetLabel(dm, "material", &material));
1174   id   = 1;
1175   initialGuess[0] = r;
1176   initialGuess[1] = NULL;
1177   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1178   id   = 2;
1179   initialGuess[0] = rp1;
1180   initialGuess[1] = NULL;
1181   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1182   id   = 1;
1183   initialGuess[0] = NULL;
1184   initialGuess[1] = phi;
1185   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1186   PetscCall(VecViewFromOptions(locX, NULL, "-local_solution_view"));
1187 
1188   PetscCall(DMGetCellDS(dm, cMax, &probh));
1189   PetscCall(PetscDSGetWeakForm(probh, &wf));
1190   PetscCall(PetscDSGetNumFields(probh, &Nf));
1191   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u, 0, NULL));
1192   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u, 0, NULL));
1193   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL));
1194   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL));
1195   if (Nf > 1) {
1196     PetscCall(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL));
1197     PetscCall(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL));
1198   }
1199   if (!rank) PetscCall(PetscDSView(probh, NULL));
1200 
1201   keys[0].label = NULL;
1202   keys[0].value = 0;
1203   keys[0].field = 0;
1204   keys[0].part  = 0;
1205   keys[1].label = material;
1206   keys[1].value = 2;
1207   keys[1].field = 0;
1208   keys[1].part  = 0;
1209   keys[2].label = fault;
1210   keys[2].value = 1;
1211   keys[2].field = 1;
1212   keys[2].part  = 0;
1213   PetscCall(VecSet(locF, 0.));
1214   PetscCall(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user));
1215   PetscCall(VecViewFromOptions(locF, NULL, "-local_residual_view"));
1216   PetscCall(MatZeroEntries(J));
1217   PetscCall(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user));
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   PetscCall(PetscInitialize(&argc, &argv, NULL,help));
1233   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1234   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1235   PetscCall(TestMesh(dm, &user));
1236   PetscCall(TestDiscretization(dm, &user));
1237   PetscCall(TestAssembly(dm, &user));
1238   PetscCall(DMDestroy(&dm));
1239   PetscCall(PetscFinalize());
1240   return 0;
1241 }
1242 
1243 /*TEST
1244   testset:
1245     args: -orig_dm_plex_check_all -dm_plex_check_all \
1246           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \
1247           -local_solution_view -local_residual_view -local_jacobian_view
1248     test:
1249       suffix: tri_0
1250       args: -dim 2
1251       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"
1252     test:
1253       suffix: tri_t1_0
1254       args: -dim 2 -test_num 1
1255       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"
1256     test:
1257       suffix: tri_t2_0
1258       args: -dim 2 -test_num 2
1259       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"
1260     test:
1261       suffix: tri_t5_0
1262       args: -dim 2 -test_num 5
1263       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"
1264     test:
1265       suffix: tet_0
1266       args: -dim 3
1267       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"
1268     test:
1269       suffix: tet_t1_0
1270       args: -dim 3 -test_num 1
1271       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"
1272 
1273   testset:
1274     args: -orig_dm_plex_check_all -dm_plex_check_all \
1275           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1276     test:
1277       suffix: tet_1
1278       nsize: 2
1279       args: -dim 3
1280       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"
1281     test:
1282       suffix: tri_1
1283       nsize: 2
1284       args: -dim 2
1285       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"
1286     test:
1287       suffix: tri_t3_0
1288       nsize: 2
1289       args: -dim 2 -test_num 3
1290       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"
1291     test:
1292       suffix: tri_t4_0
1293       nsize: 6
1294       args: -dim 2 -test_num 4
1295       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"
1296 
1297   testset:
1298     args: -orig_dm_plex_check_all -dm_plex_check_all \
1299           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1300     # 2D Quads
1301     test:
1302       suffix: quad_0
1303       args: -dim 2 -cell_simplex 0
1304       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"
1305     test:
1306       suffix: quad_1
1307       nsize: 2
1308       args: -dim 2 -cell_simplex 0
1309       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"
1310     test:
1311       suffix: quad_t1_0
1312       args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all
1313       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"
1314     test:
1315       suffix: quad_t2_0
1316       nsize: 2
1317       args: -dim 2 -cell_simplex 0 -test_num 2
1318       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"
1319     test:
1320       # TODO: The PetscSF is wrong here (connects to wrong side of split)
1321       suffix: quad_t3_0
1322       nsize: 2
1323       args: -dim 2 -cell_simplex 0 -test_num 3
1324       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"
1325     # 3D Hex
1326     test:
1327       suffix: hex_0
1328       args: -dim 3 -cell_simplex 0
1329       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"
1330     test:
1331       suffix: hex_1
1332       nsize: 2
1333       args: -dim 3 -cell_simplex 0
1334       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"
1335     test:
1336       suffix: hex_t1_0
1337       args: -dim 3 -cell_simplex 0 -test_num 1
1338       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"
1339     test:
1340       suffix: hex_t2_0
1341       args: -dim 3 -cell_simplex 0 -test_num 2
1342       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"
1343 
1344 TEST*/
1345