Actual source code: ex5.c

  1: static char help[] = "Tests for creation of hybrid meshes\n\n";

  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: */

 11: #include <petscdmplex.h>
 12: #include <petscds.h>
 13: #include <petsc/private/dmpleximpl.h>
 14: /* List of test meshes

 16: Triangle
 17: --------
 18: Test 0:
 19: Two triangles sharing a face

 21:         4
 22:       / | \
 23:      8  |  9
 24:     /   |   \
 25:    2  0 7 1  5
 26:     \   |   /
 27:      6  |  10
 28:       \ | /
 29:         3

 31: should become two triangles separated by a zero-volume cell with 4 vertices

 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

 43: Test 1:
 44: Four triangles sharing two faces which are oriented against each other

 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

 66: Fault mesh

 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

 78:        2
 79:        |
 80:   6----8----4
 81:        |    |
 82:        3    |
 83:           0-7-1
 84:             |
 85:             |
 86:             5

 88: should become four triangles separated by two zero-volume cells with 4 vertices

 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

114: Test 2:
115: Six triangles sharing one face

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

133: Test 3:
134: This is Test 2 on two processes. After the fault, we have

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

152: Test 4:
153: This is Test 2 on six processes. After the fault, we have

155: Test 5:

157:   Fault only on points 2 and 5:

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

177: Tetrahedron
178: -----------
179: Test 0:
180: Two tets sharing a face

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-------

192: should become two tetrahedrons separated by a zero-volume cell with 3 faces/3 edges/6 vertices

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

205: In parallel,

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

218: Test 1:
219: Four tets sharing two faces

221: Cells:    0-3,4-5
222: Vertices: 6-15
223: Faces:    16-29,30-34
224: Edges:    35-52,53-56

226: Quadrilateral
227: -------------
228: Test 0:
229: Two quads sharing a face

231:    5--10---4--14---7
232:    |       |       |
233:   11   0   9   1  13
234:    |       |       |
235:    2---8---3--12---6

237: should become two quads separated by a zero-volume cell with 4 vertices

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

245: Test 1:

247: Original mesh with 9 cells,

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

266: After first fault,

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

285: After second fault,

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

306:  Test 2:
307:  Two quads sharing a face in parallel

309:     4---7---3  2---8---4
310:     |       |  |       |
311:     8   0   6  5   0   7
312:     |       |  |       |
313:     1---5---2  1---6---3

315:  should become two quads separated by a zero-volume cell with 4 vertices

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

323:  Test 3:
324:  Like Test 2, but with different partition

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

332: Hexahedron
333: ----------
334: Test 0:
335: Two hexes sharing a face

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

352: should become two hexes separated by a zero-volume cell with 8 vertices

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

370: In parallel,

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

388: Test 1:

390: */

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;

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;

411:   PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
412:   PetscCall(PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL, 0));
413:   PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL, 1, 3));
414:   PetscCall(PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL));
415:   PetscCall(PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL));
416:   PetscCall(PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL, 0));
417:   PetscCall(PetscOptionsBoundedInt("-cohesive_fields", "The number of cohesive fields", "ex5.c", options->cohesiveFields, &options->cohesiveFields, NULL, 0));
418:   PetscOptionsEnd();
419:   PetscFunctionReturn(PETSC_SUCCESS);
420: }

422: static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
423: {
424:   DM          idm;
425:   PetscInt    p;
426:   PetscMPIInt rank;

428:   PetscFunctionBegin;
429:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
430:   if (rank == 0) {
431:     switch (testNum) {
432:     case 0: {
433:       PetscInt    numPoints[2]        = {4, 2};
434:       PetscInt    coneSize[6]         = {3, 3, 0, 0, 0, 0};
435:       PetscInt    cones[6]            = {2, 3, 4, 5, 4, 3};
436:       PetscInt    coneOrientations[6] = {0, 0, 0, 0, 0, 0};
437:       PetscScalar vertexCoords[8]     = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
438:       PetscInt    markerPoints[8]     = {2, 1, 3, 1, 4, 1, 5, 1};
439:       PetscInt    faultPoints[2]      = {3, 4};

441:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
442:       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
443:       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
444:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
445:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
446:     } break;
447:     case 1: {
448:       PetscInt    numPoints[2]         = {6, 4};
449:       PetscInt    coneSize[10]         = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0};
450:       PetscInt    cones[12]            = {4, 6, 5, 5, 6, 7, 8, 5, 9, 8, 4, 5};
451:       PetscInt    coneOrientations[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
452:       PetscScalar vertexCoords[12]     = {-1.0, 0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 0.0, -2.0, 1.0, -1.0, 2.0};
453:       PetscInt    markerPoints[6]      = {4, 1, 6, 1, 8, 1};
454:       PetscInt    faultPoints[3]       = {5, 6, 8};

456:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
457:       for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
458:       for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
459:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
460:       PetscCall(DMSetLabelValue(*dm, "material", 3, 1));
461:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
462:       PetscCall(DMSetLabelValue(*dm, "material", 2, 2));
463:     } break;
464:     case 2:
465:     case 3:
466:     case 4: {
467:       PetscInt    numPoints[2]         = {8, 6};
468:       PetscInt    coneSize[14]         = {3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0};
469:       PetscInt    cones[18]            = {6, 7, 9, 9, 12, 11, 7, 12, 9, 7, 8, 10, 10, 13, 12, 7, 10, 12};
470:       PetscInt    coneOrientations[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
471:       PetscScalar vertexCoords[16]     = {
472:         -1., -1., 0., -1., 1., -1., -1., 0., 1., 0., -1., 1., 0., 1., 1., 1.,
473:       };
474:       PetscInt markerPoints[16] = {6, 1, 7, 1, 8, 1, 9, 1, 10, 1, 11, 1, 12, 1, 13, 1};
475:       PetscInt faultPoints[2]   = {7, 12};

477:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
478:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
479:       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
480:       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
481:       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
482:       PetscCall(DMSetLabelValue(*dm, "material", 4, 2));
483:       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
484:       for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
485:       if (testNum == 2)
486:         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
487:       if (testNum == 3 || testNum == 4)
488:         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1));
489:     } break;
490:     case 5: {
491:       PetscInt    numPoints[2]         = {7, 6};
492:       PetscInt    coneSize[13]         = {3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0};
493:       PetscInt    cones[18]            = {6, 7, 8, 8, 9, 6, 7, 10, 8, 9, 8, 12, 8, 10, 11, 11, 12, 8};
494:       PetscInt    coneOrientations[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
495:       PetscScalar vertexCoords[14]     = {0.0, 2.0, -1.0, 1.0, 0.0, 0.0, 1.0, 1.0, -1.0, -1.0, 0.0, -2.0, 1.0, -1.0};
496:       PetscInt    faultPoints[2]       = {8, 11};
497:       PetscInt    faultBdPoints[1]     = {8};

499:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
500:       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
501:       for (p = 0; p < 1; ++p) PetscCall(DMSetLabelValue(*dm, "faultBd", faultBdPoints[p], 1));
502:       PetscCall(DMSetLabelValue(*dm, "material", 0, 0));
503:       PetscCall(DMSetLabelValue(*dm, "material", 2, 0));
504:       PetscCall(DMSetLabelValue(*dm, "material", 4, 0));
505:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
506:       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
507:       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
508:     } break;
509:     default:
510:       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
511:     }
512:   } else {
513:     PetscInt numPoints[3] = {0, 0, 0};

515:     PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL));
516:     if (testNum == 3 || testNum == 4) PetscCall(DMCreateLabel(*dm, "pfault"));
517:     else PetscCall(DMCreateLabel(*dm, "fault"));
518:   }
519:   PetscCall(DMPlexInterpolate(*dm, &idm));
520:   PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view"));
521:   PetscCall(DMDestroy(dm));
522:   *dm = idm;
523:   PetscFunctionReturn(PETSC_SUCCESS);
524: }

526: static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, AppCtx *user, DM dm)
527: {
528:   PetscInt    depth = 3, testNum = user->testNum, p;
529:   PetscMPIInt rank;

531:   PetscFunctionBegin;
532:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
533:   if (rank == 0) {
534:     switch (testNum) {
535:     case 0: {
536:       PetscInt    numPoints[4]         = {5, 7, 9, 2};
537:       PetscInt    coneSize[23]         = {4, 4, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2};
538:       PetscInt    cones[47]            = {7, 8, 9, 10, 11, 10, 13, 12, 15, 17, 14, 16, 18, 15, 14, 19, 16, 17, 18, 19, 17, 21, 20, 18, 22, 21, 22, 19, 20, 2, 3, 2, 4, 2, 5, 3, 4, 4, 5, 5, 3, 3, 6, 4, 6, 5, 6};
539:       PetscInt    coneOrientations[47] = {0, 0, 0, 0, 0, -2, 2, 2, 0, -1, -1, 0, -1, -1, 0, -1, -1, 0, 0, 0, 0, 0, -1, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
540:       PetscScalar vertexCoords[15]     = {0.0, 0.0, -0.5, 0.0, -0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5};
541:       PetscInt    markerPoints[20]     = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1};
542:       PetscInt    faultPoints[3]       = {3, 4, 5};

544:       PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
545:       for (p = 0; p < 10; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
546:       for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1));
547:       PetscCall(DMSetLabelValue(dm, "material", 0, 1));
548:       PetscCall(DMSetLabelValue(dm, "material", 1, 2));
549:     } break;
550:     case 1: {
551:       PetscInt    numPoints[4]         = {6, 13, 12, 4};
552:       PetscInt    coneSize[35]         = {4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2};
553:       PetscInt    cones[78]            = {10, 11, 12, 13, 10, 15, 16, 14, 17, 18, 14, 19, 20, 13, 19, 21, 22, 23, 24, 25, 26, 22, 24, 27, 25, 23, 26, 27, 28, 29, 23, 24, 30, 28, 22, 29, 30, 31, 32,
554:                                           28, 29, 33, 31, 32, 33, 23, 26, 34, 33, 34, 27, 32, 6,  5,  5,  7,  7,  6,  6,  4,  4,  5,  7,  4,  7,  9,  9,  5,  6,  9,  9,  8,  8,  7,  5,  8,  4,  8};
555:       PetscInt    coneOrientations[78] = {0, 0, 0, 0,  -2, 1,  0, 2,  0, 0,  -3, 0,  0,  -3, -1, 0, 0, 0, 0, 0, 0, -1, -1, 0, -1, -1, -1, -1, 0, 0, 0, 0, 0, -1, 0, -1, -1, 0, 0,
556:                                           0, 0, 0, -1, -1, -1, 0, -1, 0, -1, -1, -1, -1, 0,  0,  0, 0, 0, 0, 0, 0, 0,  0,  0, 0,  0,  0,  0,  0, 0, 0, 0, 0, 0,  0, 0,  0,  0, 0};
557:       PetscScalar vertexCoords[18]     = {-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 1.0, 0.0, 0.0};
558:       PetscInt    markerPoints[14]     = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1};
559:       PetscInt    faultPoints[4]       = {5, 6, 7, 8};

561:       PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
562:       for (p = 0; p < 7; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
563:       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1));
564:       PetscCall(DMSetLabelValue(dm, "material", 0, 1));
565:       PetscCall(DMSetLabelValue(dm, "material", 1, 2));
566:     } break;
567:     default:
568:       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
569:     }
570:   } else {
571:     PetscInt numPoints[4] = {0, 0, 0, 0};

573:     PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL));
574:     PetscCall(DMCreateLabel(dm, "fault"));
575:   }
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }

579: static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
580: {
581:   DM          idm;
582:   PetscInt    p;
583:   PetscMPIInt rank;

585:   PetscFunctionBegin;
586:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
587:   if (rank == 0) {
588:     switch (testNum) {
589:     case 0:
590:     case 2:
591:     case 3: {
592:       PetscInt    numPoints[2]        = {6, 2};
593:       PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
594:       PetscInt    cones[8]            = {2, 3, 4, 5, 3, 6, 7, 4};
595:       PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
596:       PetscScalar vertexCoords[12]    = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0, -0.5, 1.0, 0.5, 0.0, 0.5, 1.0};
597:       PetscInt    markerPoints[12]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1};
598:       PetscInt    faultPoints[2]      = {3, 4};

600:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
601:       for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
602:       if (testNum == 0)
603:         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
604:       if (testNum == 2 || testNum == 3)
605:         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1));
606:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
607:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
608:     } break;
609:     case 1: {
610:       PetscInt    numPoints[2]         = {16, 9};
611:       PetscInt    coneSize[25]         = {4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
612:       PetscInt    cones[36]            = {9, 13, 14, 10, 10, 14, 15, 11, 11, 15, 16, 12, 13, 17, 18, 14, 14, 18, 19, 15, 15, 19, 20, 16, 17, 21, 22, 18, 18, 22, 23, 19, 19, 23, 24, 20};
613:       PetscInt    coneOrientations[36] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
614:       PetscScalar vertexCoords[32]     = {-3.0, 3.0, -1.0, 3.0, 1.0, 3.0, 3.0, 3.0, -3.0, 1.0, -1.0, 1.0, 1.0, 1.0, 3.0, 1.0, -3.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3.0, -1.0, -3.0, -3.0, -1.0, -3.0, 1.0, -3.0, 3.0, -3.0};
615:       PetscInt    faultPoints[4]       = {11, 15, 19, 20};
616:       PetscInt    fault2Points[2]      = {17, 18};

618:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
619:       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
620:       for (p = 3; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "faultBd", faultPoints[p], 1));
621:       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault2", fault2Points[p], 1));
622:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
623:       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
624:       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
625:       PetscCall(DMSetLabelValue(*dm, "material", 3, 1));
626:       PetscCall(DMSetLabelValue(*dm, "material", 4, 1));
627:       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
628:       PetscCall(DMSetLabelValue(*dm, "material", 6, 2));
629:       PetscCall(DMSetLabelValue(*dm, "material", 7, 2));
630:       PetscCall(DMSetLabelValue(*dm, "material", 8, 2));
631:     } break;
632:     default:
633:       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
634:     }
635:   } else {
636:     PetscInt numPoints[3] = {0, 0, 0};

638:     PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL));
639:     if (testNum == 2 || testNum == 3) PetscCall(DMCreateLabel(*dm, "pfault"));
640:     else PetscCall(DMCreateLabel(*dm, "fault"));
641:   }
642:   PetscCall(DMPlexInterpolate(*dm, &idm));
643:   PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view"));
644:   PetscCall(DMDestroy(dm));
645:   *dm = idm;
646:   PetscFunctionReturn(PETSC_SUCCESS);
647: }

649: static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
650: {
651:   DM          idm;
652:   PetscInt    depth = 3, p;
653:   PetscMPIInt rank;

655:   PetscFunctionBegin;
656:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
657:   if (rank == 0) {
658:     switch (testNum) {
659:     case 0: {
660:       PetscInt    numPoints[2]         = {12, 2};
661:       PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
662:       PetscInt    cones[16]            = {2, 5, 4, 3, 6, 7, 8, 9, 3, 4, 11, 10, 7, 12, 13, 8};
663:       PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
664:       PetscScalar vertexCoords[36]     = {-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -0.5, 1.0, 0.0, -0.5, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, -0.5, 1.0, 1.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0};
665:       PetscInt    markerPoints[52]     = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1, 8, 1, 9, 1};
666:       PetscInt    faultPoints[4]       = {3, 4, 7, 8};

668:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
669:       PetscCall(DMPlexInterpolate(*dm, &idm));
670:       for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(idm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
671:       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
672:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
673:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
674:     } break;
675:     case 1: {
676:       /* Cell Adjacency Graph:
677:         0 -- { 8, 13, 21, 24} --> 1
678:         0 -- {20, 21, 23, 24} --> 5 F
679:         1 -- {10, 15, 21, 24} --> 2
680:         1 -- {13, 14, 15, 24} --> 6
681:         2 -- {21, 22, 24, 25} --> 4 F
682:         3 -- {21, 24, 30, 35} --> 4
683:         3 -- {21, 24, 28, 33} --> 5
684:        */
685:       PetscInt numPoints[2] = {30, 7};
686:       PetscInt coneSize[37] = {8, 8, 8, 8, 8, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
687:       PetscInt cones[56] = {8, 21, 20, 7, 13, 12, 23, 24, 14, 15, 10, 9, 13, 8, 21, 24, 15, 16, 11, 10, 24, 21, 22, 25, 30, 29, 28, 21, 35, 24, 33, 34, 24, 21, 30, 35, 25, 36, 31, 22, 27, 20, 21, 28, 32, 33, 24, 23, 15, 24, 13, 14, 19, 18, 17, 26};
688:       PetscInt coneOrientations[56] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
689:       PetscScalar vertexCoords[90]  = {-2.0, -2.0, -2.0, -2.0, -1.0, -2.0, -3.0, 0.0, -2.0, -2.0, 1.0,  -2.0, -2.0, 2.0, -2.0, -2.0, -2.0, 0.0,  -2.0, -1.0, 0.0, -3.0, 0.0, 0.0, -2.0, 1.0, 0.0, -2.0, 2.0, 0.0,
690:                                        -2.0, -1.0, 2.0,  -3.0, 0.0,  2.0,  -2.0, 1.0, 2.0,  0.0,  -2.0, -2.0, 0.0,  0.0, -2.0, 0.0,  2.0,  -2.0, 0.0,  -2.0, 0.0, 0.0,  0.0, 0.0, 0.0,  2.0, 0.0, 0.0,  0.0, 2.0,
691:                                        2.0,  -2.0, -2.0, 2.0,  -1.0, -2.0, 3.0,  0.0, -2.0, 2.0,  1.0,  -2.0, 2.0,  2.0, -2.0, 2.0,  -2.0, 0.0,  2.0,  -1.0, 0.0, 3.0,  0.0, 0.0, 2.0,  1.0, 0.0, 2.0,  2.0, 0.0};
692:       PetscInt    faultPoints[6]    = {20, 21, 22, 23, 24, 25};

694:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
695:       PetscCall(DMPlexInterpolate(*dm, &idm));
696:       for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
697:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
698:       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
699:       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
700:       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
701:       PetscCall(DMSetLabelValue(*dm, "material", 4, 2));
702:       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
703:       PetscCall(DMSetLabelValue(*dm, "material", 6, 2));
704:     } break;
705:     case 2: {
706:       /* Buried fault edge */
707:       PetscInt    numPoints[2]         = {18, 4};
708:       PetscInt    coneSize[22]         = {8, 8, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
709:       PetscInt    cones[32]            = {4, 5, 8, 7, 13, 16, 17, 14, 5, 6, 9, 8, 14, 17, 18, 15, 7, 8, 11, 10, 16, 19, 20, 17, 8, 9, 12, 11, 17, 20, 21, 18};
710:       PetscInt    coneOrientations[32] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
711:       PetscScalar vertexCoords[54]     = {-2.0, -2.0, 0.0, -2.0, 0.0, 0.0, -2.0, 2.0, 0.0, 0.0, -2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 2.0, -2.0, 0.0, 2.0, 0.0, 0.0, 2.0, 2.0, 0.0,
712:                                           -2.0, -2.0, 2.0, -2.0, 0.0, 2.0, -2.0, 2.0, 2.0, 0.0, -2.0, 2.0, 0.0, 0.0, 2.0, 0.0, 2.0, 2.0, 2.0, -2.0, 2.0, 2.0, 0.0, 2.0, 2.0, 2.0, 2.0};
713:       PetscInt    faultPoints[4]       = {7, 8, 16, 17};

715:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
716:       PetscCall(DMPlexInterpolate(*dm, &idm));
717:       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
718:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
719:       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
720:       PetscCall(DMSetLabelValue(*dm, "material", 2, 2));
721:       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
722:     } break;
723:     default:
724:       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
725:     }
726:   } else {
727:     PetscInt numPoints[4] = {0, 0, 0, 0};

729:     PetscCall(DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL));
730:     PetscCall(DMPlexInterpolate(*dm, &idm));
731:     PetscCall(DMCreateLabel(idm, "fault"));
732:   }
733:   PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view"));
734:   PetscCall(DMDestroy(dm));
735:   *dm = idm;
736:   PetscFunctionReturn(PETSC_SUCCESS);
737: }

739: static PetscErrorCode CreateFaultLabel(DM dm)
740: {
741:   DMLabel  label;
742:   PetscInt dim, h, pStart, pEnd, pMax, p;

744:   PetscFunctionBegin;
745:   PetscCall(DMGetDimension(dm, &dim));
746:   PetscCall(DMCreateLabel(dm, "cohesive"));
747:   PetscCall(DMGetLabel(dm, "cohesive", &label));
748:   for (h = 0; h <= dim; ++h) {
749:     PetscCall(DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax));
750:     PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd));
751:     for (p = pMax; p < pEnd; ++p) PetscCall(DMLabelSetValue(label, p, 1));
752:   }
753:   PetscFunctionReturn(PETSC_SUCCESS);
754: }

756: static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user)
757: {
758:   PetscFE  fe;
759:   DMLabel  fault;
760:   PetscInt dim, Ncf = user->cohesiveFields, f;

762:   PetscFunctionBegin;
763:   PetscCall(DMGetDimension(dm, &dim));
764:   PetscCall(DMGetLabel(dm, "cohesive", &fault));
765:   PetscCall(DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD));

767:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe));
768:   PetscCall(PetscFESetName(fe, "displacement"));
769:   PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
770:   PetscCall(PetscFEDestroy(&fe));

772:   if (Ncf > 0) {
773:     PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe));
774:     PetscCall(PetscFESetName(fe, "fault traction"));
775:     PetscCall(DMAddField(dm, fault, (PetscObject)fe));
776:     PetscCall(PetscFEDestroy(&fe));
777:   }
778:   for (f = 1; f < Ncf; ++f) {
779:     char name[256], opt[256];

781:     PetscCall(PetscSNPrintf(name, 256, "fault field %" PetscInt_FMT, f));
782:     PetscCall(PetscSNPrintf(opt, 256, "faultfield_%" PetscInt_FMT "_", f));
783:     PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe));
784:     PetscCall(PetscFESetName(fe, name));
785:     PetscCall(DMAddField(dm, fault, (PetscObject)fe));
786:     PetscCall(PetscFEDestroy(&fe));
787:   }

789:   PetscCall(DMCreateDS(dm));
790:   PetscFunctionReturn(PETSC_SUCCESS);
791: }

793: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
794: {
795:   PetscInt    dim         = user->dim;
796:   PetscBool   cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault;
797:   PetscMPIInt rank, size;
798:   DMLabel     matLabel;

800:   PetscFunctionBegin;
801:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
802:   PetscCallMPI(MPI_Comm_size(comm, &size));
803:   PetscCall(DMCreate(comm, dm));
804:   PetscCall(DMSetType(*dm, DMPLEX));
805:   PetscCall(DMSetDimension(*dm, dim));
806:   switch (dim) {
807:   case 2:
808:     if (cellSimplex) {
809:       PetscCall(CreateSimplex_2D(comm, user->testNum, dm));
810:     } else {
811:       PetscCall(CreateQuad_2D(comm, user->testNum, dm));
812:     }
813:     break;
814:   case 3:
815:     if (cellSimplex) {
816:       PetscCall(CreateSimplex_3D(comm, user, *dm));
817:     } else {
818:       PetscCall(CreateHex_3D(comm, user->testNum, dm));
819:     }
820:     break;
821:   default:
822:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %" PetscInt_FMT, dim);
823:   }
824:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "orig_"));
825:   PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
826:   PetscCall(DMSetFromOptions(*dm));
827:   PetscCall(DMGetLabel(*dm, "material", &matLabel));
828:   if (matLabel) PetscCall(DMPlexLabelComplete(*dm, matLabel));
829:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
830:   PetscCall(DMHasLabel(*dm, "fault", &hasFault));
831:   if (hasFault) {
832:     DM      dmHybrid = NULL, dmInterface = NULL;
833:     DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel;

835:     PetscCall(DMGetLabel(*dm, "fault", &faultLabel));
836:     PetscCall(DMGetLabel(*dm, "faultBd", &faultBdLabel));
837:     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid));
838:     PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD));
839:     PetscCall(DMLabelDestroy(&hybridLabel));
840:     PetscCall(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD));
841:     PetscCall(DMLabelDestroy(&splitLabel));
842:     PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_interface_view"));
843:     PetscCall(DMDestroy(&dmInterface));
844:     PetscCall(DMDestroy(dm));
845:     *dm = dmHybrid;
846:   }
847:   PetscCall(DMHasLabel(*dm, "fault2", &hasFault2));
848:   if (hasFault2) {
849:     DM      dmHybrid = NULL;
850:     DMLabel faultLabel, faultBdLabel, hybridLabel;

852:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "faulted_"));
853:     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre"));
854:     PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
855:     PetscCall(DMSetFromOptions(*dm));
856:     PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
857:     PetscCall(DMGetLabel(*dm, "fault2", &faultLabel));
858:     PetscCall(DMGetLabel(*dm, "fault2Bd", &faultBdLabel));
859:     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, NULL, &dmHybrid));
860:     PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD));
861:     PetscCall(DMLabelDestroy(&hybridLabel));
862:     PetscCall(DMDestroy(dm));
863:     *dm = dmHybrid;
864:   }
865:   if (user->testPartition && size > 1) {
866:     PetscPartitioner part;
867:     PetscInt        *sizes  = NULL;
868:     PetscInt        *points = NULL;

870:     if (rank == 0) {
871:       if (dim == 2 && cellSimplex && size == 2) {
872:         switch (user->testNum) {
873:         case 0: {
874:           PetscInt triSizes_p2[2]  = {1, 2};
875:           PetscInt triPoints_p2[3] = {0, 1, 2};

877:           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
878:           PetscCall(PetscArraycpy(sizes, triSizes_p2, 2));
879:           PetscCall(PetscArraycpy(points, triPoints_p2, 3));
880:           break;
881:         }
882:         case 3: {
883:           PetscInt triSizes_p2[2]  = {3, 3};
884:           PetscInt triPoints_p2[6] = {0, 1, 2, 3, 4, 5};

886:           PetscCall(PetscMalloc2(2, &sizes, 6, &points));
887:           PetscCall(PetscArraycpy(sizes, triSizes_p2, 2));
888:           PetscCall(PetscArraycpy(points, triPoints_p2, 6));
889:           break;
890:         }
891:         default:
892:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 2 procs", user->testNum);
893:         }
894:       } else if (dim == 2 && cellSimplex && size == 6) {
895:         switch (user->testNum) {
896:         case 4: {
897:           PetscInt triSizes_p6[6]  = {1, 1, 1, 1, 1, 1};
898:           PetscInt triPoints_p6[6] = {0, 1, 2, 3, 4, 5};

900:           PetscCall(PetscMalloc2(6, &sizes, 6, &points));
901:           PetscCall(PetscArraycpy(sizes, triSizes_p6, 6));
902:           PetscCall(PetscArraycpy(points, triPoints_p6, 6));
903:           break;
904:         }
905:         default:
906:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 6 procs", user->testNum);
907:         }
908:       } else if (dim == 2 && !cellSimplex && size == 2) {
909:         switch (user->testNum) {
910:         case 0: {
911:           PetscInt quadSizes_p2[2]  = {1, 2};
912:           PetscInt quadPoints_p2[3] = {0, 1, 2};

914:           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
915:           PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2));
916:           PetscCall(PetscArraycpy(points, quadPoints_p2, 3));
917:           break;
918:         }
919:         case 2: {
920:           PetscInt quadSizes_p2[2]  = {1, 1};
921:           PetscInt quadPoints_p2[2] = {0, 1};

923:           PetscCall(PetscMalloc2(2, &sizes, 2, &points));
924:           PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2));
925:           PetscCall(PetscArraycpy(points, quadPoints_p2, 2));
926:           break;
927:         }
928:         case 3: {
929:           PetscInt quadSizes_p2[2]  = {1, 1};
930:           PetscInt quadPoints_p2[2] = {1, 0};

932:           PetscCall(PetscMalloc2(2, &sizes, 2, &points));
933:           PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2));
934:           PetscCall(PetscArraycpy(points, quadPoints_p2, 2));
935:           break;
936:         }
937:         default:
938:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for quadrilateral mesh on 2 procs", user->testNum);
939:         }
940:       } else if (dim == 3 && cellSimplex && size == 2) {
941:         switch (user->testNum) {
942:         case 0: {
943:           PetscInt tetSizes_p2[2]  = {1, 2};
944:           PetscInt tetPoints_p2[3] = {0, 1, 2};

946:           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
947:           PetscCall(PetscArraycpy(sizes, tetSizes_p2, 2));
948:           PetscCall(PetscArraycpy(points, tetPoints_p2, 3));
949:           break;
950:         }
951:         default:
952:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for teterehedral mesh on 2 procs", user->testNum);
953:         }
954:       } else if (dim == 3 && !cellSimplex && size == 2) {
955:         switch (user->testNum) {
956:         case 0: {
957:           PetscInt hexSizes_p2[2]  = {1, 2};
958:           PetscInt hexPoints_p2[3] = {0, 1, 2};

960:           PetscCall(PetscMalloc2(2, &sizes, 3, &points));
961:           PetscCall(PetscArraycpy(sizes, hexSizes_p2, 2));
962:           PetscCall(PetscArraycpy(points, hexPoints_p2, 3));
963:           break;
964:         }
965:         default:
966:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for hexahedral mesh on 2 procs", user->testNum);
967:         }
968:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition");
969:     }
970:     PetscCall(DMPlexGetPartitioner(*dm, &part));
971:     PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL));
972:     PetscCall(PetscPartitionerShellSetPartition(part, size, sizes, points));
973:     PetscCall(PetscFree2(sizes, points));
974:   }
975:   {
976:     DM pdm = NULL;

978:     /* Distribute mesh over processes */
979:     PetscCall(DMPlexDistribute(*dm, 0, NULL, &pdm));
980:     if (pdm) {
981:       PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view"));
982:       PetscCall(DMDestroy(dm));
983:       *dm = pdm;
984:     }
985:   }
986:   PetscCall(DMHasLabel(*dm, "pfault", &hasParallelFault));
987:   if (hasParallelFault) {
988:     DM      dmHybrid = NULL, dmInterface;
989:     DMLabel faultLabel, faultBdLabel, hybridLabel;

991:     PetscCall(DMGetLabel(*dm, "pfault", &faultLabel));
992:     PetscCall(DMGetLabel(*dm, "pfaultBd", &faultBdLabel));
993:     PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, &dmInterface, &dmHybrid));
994:     PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_fault_view"));
995:     {
996:       PetscViewer viewer;
997:       PetscMPIInt rank;

999:       PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)*dm), &rank));
1000:       PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer));
1001:       PetscCall(PetscViewerASCIIPrintf(viewer, "Rank %d\n", rank));
1002:       PetscCall(DMLabelView(hybridLabel, viewer));
1003:       PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer));
1004:       PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1005:     }
1006:     PetscCall(DMLabelDestroy(&hybridLabel));
1007:     PetscCall(DMDestroy(&dmInterface));
1008:     PetscCall(DMDestroy(dm));
1009:     *dm = dmHybrid;
1010:   }
1011:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Hybrid Mesh"));
1012:   PetscCall(CreateFaultLabel(*dm));
1013:   PetscCall(CreateDiscretization(*dm, user));
1014:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre"));
1015:   PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
1016:   PetscCall(DMSetFromOptions(*dm));
1017:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1018:   PetscFunctionReturn(PETSC_SUCCESS);
1019: }

1021: static PetscErrorCode TestMesh(DM dm, AppCtx *user)
1022: {
1023:   PetscFunctionBegin;
1024:   PetscCall(DMPlexCheckSymmetry(dm));
1025:   PetscCall(DMPlexCheckSkeleton(dm, 0));
1026:   PetscCall(DMPlexCheckFaces(dm, 0));
1027:   PetscFunctionReturn(PETSC_SUCCESS);
1028: }

1030: static PetscErrorCode TestDiscretization(DM dm, AppCtx *user)
1031: {
1032:   PetscSection s;

1034:   PetscFunctionBegin;
1035:   PetscCall(DMGetSection(dm, &s));
1036:   PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-local_section_view"));
1037:   PetscFunctionReturn(PETSC_SUCCESS);
1038: }

1040: static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1041: {
1042:   PetscInt d;
1043:   for (d = 0; d < dim; ++d) u[d] = x[d];
1044:   return PETSC_SUCCESS;
1045: }

1047: static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1048: {
1049:   PetscInt d;
1050:   for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0);
1051:   return PETSC_SUCCESS;
1052: }

1054: static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1055: {
1056:   PetscInt d;
1057:   u[0] = -x[1];
1058:   u[1] = x[0];
1059:   for (d = 2; d < dim; ++d) u[d] = x[d];
1060:   return PETSC_SUCCESS;
1061: }

1063: /* \lambda \cdot (\psi_u^- - \psi_u^+) */
1064: static void f0_bd_u_neg(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1065: {
1066:   const PetscInt Nc = dim + 1;
1067:   for (PetscInt c = 0; c < Nc; ++c) f0[c] = -u[uOff[1] + c];
1068: }

1070: static void f0_bd_u_pos(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1071: {
1072:   const PetscInt Nc = dim + 1;
1073:   for (PetscInt c = 0; c < Nc; ++c) f0[c] = u[uOff[1] + c];
1074: }

1076: /* (d - u^+ + u^-) \cdot \psi_\lambda */
1077: static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1078: {
1079:   const PetscInt Nc = uOff[2] - uOff[1];

1081:   for (PetscInt c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc + c];
1082: }

1084: /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */
1085: static void g0_bd_ul_neg(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1086: {
1087:   const PetscInt Nc = dim + 1;
1088:   for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = -1.0;
1089: }

1091: static void g0_bd_ul_pos(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1092: {
1093:   const PetscInt Nc = dim + 1;
1094:   for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = 1.0;
1095: }

1097: /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */
1098: static void g0_bd_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1099: {
1100:   const PetscInt Nc = uOff[2] - uOff[1];

1102:   for (PetscInt c = 0; c < Nc; ++c) {
1103:     g0[c * Nc + c]           = -1.0;
1104:     g0[Nc * Nc + c * Nc + c] = 1.0;
1105:   }
1106: }

1108: static PetscErrorCode TestAssembly(DM dm, AppCtx *user)
1109: {
1110:   Mat           J;
1111:   Vec           locX, locF;
1112:   PetscDS       probh;
1113:   DMLabel       fault, material;
1114:   IS            cohesiveCells;
1115:   PetscWeakForm wf;
1116:   PetscFormKey  keys[3];
1117:   PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx);
1118:   PetscInt    dim, Nf, cMax, cEnd, id;
1119:   PetscMPIInt rank;

1121:   PetscFunctionBegin;
1122:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1123:   PetscCall(DMGetDimension(dm, &dim));
1124:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax));
1125:   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
1126:   PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells));
1127:   PetscCall(DMGetLabel(dm, "cohesive", &fault));
1128:   PetscCall(DMGetLocalVector(dm, &locX));
1129:   PetscCall(PetscObjectSetName((PetscObject)locX, "Local Solution"));
1130:   PetscCall(DMGetLocalVector(dm, &locF));
1131:   PetscCall(PetscObjectSetName((PetscObject)locF, "Local Residual"));
1132:   PetscCall(DMCreateMatrix(dm, &J));
1133:   PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));

1135:   /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */
1136:   PetscCall(DMGetLabel(dm, "material", &material));
1137:   id              = 1;
1138:   initialGuess[0] = r;
1139:   initialGuess[1] = NULL;
1140:   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1141:   id              = 2;
1142:   initialGuess[0] = rp1;
1143:   initialGuess[1] = NULL;
1144:   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1145:   id              = 1;
1146:   initialGuess[0] = NULL;
1147:   initialGuess[1] = phi;
1148:   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1149:   PetscCall(VecViewFromOptions(locX, NULL, "-local_solution_view"));

1151:   PetscCall(DMGetCellDS(dm, cMax, &probh, NULL));
1152:   PetscCall(PetscDSGetWeakForm(probh, &wf));
1153:   PetscCall(PetscDSGetNumFields(probh, &Nf));
1154:   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u_neg, 0, NULL));
1155:   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u_pos, 0, NULL));
1156:   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul_neg, 0, NULL, 0, NULL, 0, NULL));
1157:   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul_pos, 0, NULL, 0, NULL, 0, NULL));
1158:   if (Nf > 1) {
1159:     PetscCall(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL));
1160:     PetscCall(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL));
1161:   }
1162:   if (rank == 0) PetscCall(PetscDSView(probh, NULL));

1164:   keys[0].label = material;
1165:   keys[0].value = 1;
1166:   keys[0].field = 0;
1167:   keys[0].part  = 0;
1168:   keys[1].label = material;
1169:   keys[1].value = 2;
1170:   keys[1].field = 0;
1171:   keys[1].part  = 0;
1172:   keys[2].label = fault;
1173:   keys[2].value = 1;
1174:   keys[2].field = 1;
1175:   keys[2].part  = 0;
1176:   PetscCall(VecSet(locF, 0.));
1177:   PetscCall(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user));
1178:   PetscCall(VecViewFromOptions(locF, NULL, "-local_residual_view"));
1179:   PetscCall(MatZeroEntries(J));
1180:   PetscCall(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user));
1181:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1182:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1183:   PetscCall(MatViewFromOptions(J, NULL, "-local_jacobian_view"));

1185:   PetscCall(DMRestoreLocalVector(dm, &locX));
1186:   PetscCall(DMRestoreLocalVector(dm, &locF));
1187:   PetscCall(MatDestroy(&J));
1188:   PetscCall(ISDestroy(&cohesiveCells));

1190:   if (cMax < cEnd) {
1191:     PetscDS         ds;
1192:     PetscFE         fe;
1193:     PetscQuadrature quad;
1194:     IS             *perm;
1195:     const PetscInt *cone;
1196:     PetscInt        Na, a;

1198:     PetscCall(DMPlexGetCone(dm, cMax, &cone));
1199:     PetscCall(DMGetCellDS(dm, cMax, &ds, NULL));
1200:     PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
1201:     PetscCall(PetscFEGetQuadrature(fe, &quad));
1202:     PetscCall(PetscQuadratureComputePermutations(quad, &Na, &perm));
1203:     for (a = 0; a < Na; ++a) PetscCall(ISDestroy(&perm[a]));
1204:     PetscCall(PetscFree(perm));
1205:   }
1206:   PetscFunctionReturn(PETSC_SUCCESS);
1207: }

1209: int main(int argc, char **argv)
1210: {
1211:   DM     dm;
1212:   AppCtx user; /* user-defined work context */

1214:   PetscFunctionBeginUser;
1215:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1216:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1217:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1218:   PetscCall(TestMesh(dm, &user));
1219:   PetscCall(TestDiscretization(dm, &user));
1220:   PetscCall(TestAssembly(dm, &user));
1221:   PetscCall(DMDestroy(&dm));
1222:   PetscCall(PetscFinalize());
1223:   return 0;
1224: }

1226: /*TEST
1227:   testset:
1228:     args: -orig_dm_plex_check_all -dm_plex_check_all \
1229:           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \
1230:           -local_solution_view -local_residual_view -local_jacobian_view
1231:     test:
1232:       suffix: tri_0
1233:       args: -dim 2
1234:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1235:     test:
1236:       suffix: tri_t1_0
1237:       args: -dim 2 -test_num 1
1238:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1239:     test:
1240:       suffix: tri_t2_0
1241:       args: -dim 2 -test_num 2
1242:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1244:     test:
1245:       suffix: tri_t5_0
1246:       args: -dim 2 -test_num 5
1247:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1249:     test:
1250:       suffix: tet_0
1251:       args: -dim 3
1252:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1254:     test:
1255:       suffix: tet_t1_0
1256:       args: -dim 3 -test_num 1
1257:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1259:   testset:
1260:     args: -orig_dm_plex_check_all -dm_plex_check_all \
1261:           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1262:     test:
1263:       suffix: tet_1
1264:       nsize: 2
1265:       args: -dim 3
1266:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1268:     test:
1269:       suffix: tri_1
1270:       nsize: 2
1271:       args: -dim 2
1272:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1274:     test:
1275:       suffix: tri_t3_0
1276:       nsize: 2
1277:       args: -dim 2 -test_num 3
1278:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1280:     test:
1281:       suffix: tri_t4_0
1282:       nsize: 6
1283:       args: -dim 2 -test_num 4
1284:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1286:   testset:
1287:     args: -orig_dm_plex_check_all -dm_plex_check_all \
1288:           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1289:     # 2D Quads
1290:     test:
1291:       suffix: quad_0
1292:       args: -dim 2 -cell_simplex 0
1293:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1295:     test:
1296:       suffix: quad_1
1297:       nsize: 2
1298:       args: -dim 2 -cell_simplex 0
1299:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1301:     test:
1302:       suffix: quad_t1_0
1303:       args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all
1304:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1306:     test:
1307:       suffix: quad_t2_0
1308:       nsize: 2
1309:       args: -dim 2 -cell_simplex 0 -test_num 2
1310:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1312:     test:
1313:       # TODO: The PetscSF is wrong here (connects to wrong side of split)
1314:       suffix: quad_t3_0
1315:       nsize: 2
1316:       args: -dim 2 -cell_simplex 0 -test_num 3
1317:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1319:     # 3D Hex
1320:     test:
1321:       suffix: hex_0
1322:       args: -dim 3 -cell_simplex 0
1323:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1324:     test:
1325:       suffix: hex_1
1326:       nsize: 2
1327:       args: -dim 3 -cell_simplex 0
1328:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1329:     test:
1330:       suffix: hex_t1_0
1331:       args: -dim 3 -cell_simplex 0 -test_num 1
1332:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"
1333:     test:
1334:       suffix: hex_t2_0
1335:       args: -dim 3 -cell_simplex 0 -test_num 2
1336:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1338: TEST*/