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:   PetscBool testAssembly;   // Flag for assembly test
398:   PetscInt  testNum;        /* The particular mesh to test */
399:   PetscInt  cohesiveFields; /* The number of cohesive fields */
400: } AppCtx;

402: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
403: {
404:   PetscFunctionBegin;
405:   options->debug          = 0;
406:   options->dim            = 2;
407:   options->cellSimplex    = PETSC_TRUE;
408:   options->testPartition  = PETSC_TRUE;
409:   options->testAssembly   = PETSC_TRUE;
410:   options->testNum        = 0;
411:   options->cohesiveFields = 1;

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

425: static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
426: {
427:   DM          idm;
428:   PetscInt    p;
429:   PetscMPIInt rank;

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

444:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
445:       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
446:       for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
447:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
448:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
449:     } break;
450:     case 1: {
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};

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));
463:       PetscCall(DMSetLabelValue(*dm, "material", 3, 1));
464:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
465:       PetscCall(DMSetLabelValue(*dm, "material", 2, 2));
466:     } break;
467:     case 2:
468:     case 3:
469:     case 4: {
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, 7, 8, 10, 10, 13, 12, 7, 10, 12};
473:       PetscInt    coneOrientations[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
474:       PetscScalar vertexCoords[16]     = {
475:         -1., -1., 0., -1., 1., -1., -1., 0., 1., 0., -1., 1., 0., 1., 1., 1.,
476:       };
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};

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)
489:         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
490:       if (testNum == 3 || testNum == 4)
491:         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1));
492:     } break;
493:     case 5: {
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};

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));
506:       PetscCall(DMSetLabelValue(*dm, "material", 2, 0));
507:       PetscCall(DMSetLabelValue(*dm, "material", 4, 0));
508:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
509:       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
510:       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
511:     } break;
512:     default:
513:       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %" PetscInt_FMT, testNum);
514:     }
515:   } else {
516:     PetscInt numPoints[3] = {0, 0, 0};

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

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

534:   PetscFunctionBegin;
535:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
536:   if (rank == 0) {
537:     switch (testNum) {
538:     case 0: {
539:       PetscInt    numPoints[4]         = {5, 7, 9, 2};
540:       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};
541:       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};
542:       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};
543:       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};
544:       PetscInt    markerPoints[20]     = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1};
545:       PetscInt    faultPoints[3]       = {3, 4, 5};

547:       PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
548:       for (p = 0; p < 10; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
549:       for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1));
550:       PetscCall(DMSetLabelValue(dm, "material", 0, 1));
551:       PetscCall(DMSetLabelValue(dm, "material", 1, 2));
552:     } break;
553:     case 1: {
554:       PetscInt    numPoints[4]         = {6, 13, 12, 4};
555:       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};
556:       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,
557:                                           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};
558:       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,
559:                                           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};
560:       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};
561:       PetscInt    markerPoints[14]     = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1};
562:       PetscInt    faultPoints[4]       = {5, 6, 7, 8};

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

576:     PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL));
577:     PetscCall(DMCreateLabel(dm, "fault"));
578:   }
579:   PetscFunctionReturn(PETSC_SUCCESS);
580: }

582: static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
583: {
584:   DM          idm;
585:   PetscInt    p;
586:   PetscMPIInt rank;

588:   PetscFunctionBegin;
589:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
590:   if (rank == 0) {
591:     switch (testNum) {
592:     case 0:
593:     case 2:
594:     case 3: {
595:       PetscInt    numPoints[2]        = {6, 2};
596:       PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
597:       PetscInt    cones[8]            = {2, 3, 4, 5, 3, 6, 7, 4};
598:       PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
599:       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};
600:       PetscInt    markerPoints[12]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1};
601:       PetscInt    faultPoints[2]      = {3, 4};

603:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
604:       for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
605:       if (testNum == 0)
606:         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1));
607:       if (testNum == 2 || testNum == 3)
608:         for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1));
609:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
610:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
611:     } break;
612:     case 1: {
613:       PetscInt    numPoints[2]         = {16, 9};
614:       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};
615:       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};
616:       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};
617:       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};
618:       PetscInt    faultPoints[4]       = {11, 15, 19, 20};
619:       PetscInt    fault2Points[2]      = {17, 18};

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

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

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

658:   PetscFunctionBegin;
659:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
660:   if (rank == 0) {
661:     switch (testNum) {
662:     case 0: {
663:       PetscInt    numPoints[2]         = {12, 2};
664:       PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
665:       PetscInt    cones[16]            = {2, 5, 4, 3, 6, 7, 8, 9, 3, 4, 11, 10, 7, 12, 13, 8};
666:       PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
667:       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};
668:       PetscInt    markerPoints[52]     = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1, 8, 1, 9, 1};
669:       PetscInt    faultPoints[4]       = {3, 4, 7, 8};

671:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
672:       PetscCall(DMPlexInterpolate(*dm, &idm));
673:       for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(idm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
674:       for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
675:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
676:       PetscCall(DMSetLabelValue(*dm, "material", 1, 2));
677:     } break;
678:     case 1: {
679:       /* Cell Adjacency Graph:
680:         0 -- { 8, 13, 21, 24} --> 1
681:         0 -- {20, 21, 23, 24} --> 5 F
682:         1 -- {10, 15, 21, 24} --> 2
683:         1 -- {13, 14, 15, 24} --> 6
684:         2 -- {21, 22, 24, 25} --> 4 F
685:         3 -- {21, 24, 30, 35} --> 4
686:         3 -- {21, 24, 28, 33} --> 5
687:        */
688:       PetscInt numPoints[2] = {30, 7};
689:       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};
690:       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};
691:       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};
692:       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,
693:                                        -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,
694:                                        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};
695:       PetscInt    faultPoints[6]    = {20, 21, 22, 23, 24, 25};

697:       PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords));
698:       PetscCall(DMPlexInterpolate(*dm, &idm));
699:       for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1));
700:       PetscCall(DMSetLabelValue(*dm, "material", 0, 1));
701:       PetscCall(DMSetLabelValue(*dm, "material", 1, 1));
702:       PetscCall(DMSetLabelValue(*dm, "material", 2, 1));
703:       PetscCall(DMSetLabelValue(*dm, "material", 3, 2));
704:       PetscCall(DMSetLabelValue(*dm, "material", 4, 2));
705:       PetscCall(DMSetLabelValue(*dm, "material", 5, 2));
706:       PetscCall(DMSetLabelValue(*dm, "material", 6, 2));
707:     } break;
708:     case 2: {
709:       /* Buried fault edge */
710:       PetscInt    numPoints[2]         = {18, 4};
711:       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};
712:       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};
713:       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};
714:       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,
715:                                           -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};
716:       PetscInt    faultPoints[4]       = {7, 8, 16, 17};

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

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

742: static PetscErrorCode CreateFaultLabel(DM dm)
743: {
744:   DMLabel  label;
745:   PetscInt dim, h, pStart, pEnd, pMax, p;

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

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

765:   PetscFunctionBegin;
766:   PetscCall(DMGetDimension(dm, &dim));
767:   PetscCall(DMGetLabel(dm, "cohesive", &fault));
768:   PetscCall(DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD));

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

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

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

792:   PetscCall(DMCreateDS(dm));
793:   PetscFunctionReturn(PETSC_SUCCESS);
794: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1036:   PetscFunctionBegin;
1037:   PetscCall(DMGetLocalSection(dm, &s));
1038:   PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-local_section_view"));
1039:   PetscFunctionReturn(PETSC_SUCCESS);
1040: }

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

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

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

1065: static void add_fields(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 f[])
1066: {
1067:   PetscInt       d;
1068:   const PetscInt offN = 0;
1069:   const PetscInt offP = dim;
1070:   for (d = 0; d < dim; ++d) f[d] = u[offN + d] + u[offP + d];
1071: }

1073: static void normal_field(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 f[])
1074: {
1075:   PetscInt d;
1076:   for (d = 0; d < dim; ++d) f[d] = n[d];
1077: }

1079: /* \lambda \cdot (\psi_u^- - \psi_u^+) */
1080: 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[])
1081: {
1082:   const PetscInt Nc = dim + 1;
1083:   for (PetscInt c = 0; c < Nc; ++c) f0[c] = -u[uOff[1] + c];
1084: }

1086: 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[])
1087: {
1088:   const PetscInt Nc = dim + 1;
1089:   for (PetscInt c = 0; c < Nc; ++c) f0[c] = u[uOff[1] + c];
1090: }

1092: /* (d - u^+ + u^-) \cdot \psi_\lambda */
1093: 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[])
1094: {
1095:   const PetscInt Nc = uOff[2] - uOff[1];

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

1100: /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */
1101: 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[])
1102: {
1103:   const PetscInt Nc = dim + 1;
1104:   for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = -1.0;
1105: }

1107: 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[])
1108: {
1109:   const PetscInt Nc = dim + 1;
1110:   for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = 1.0;
1111: }

1113: /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */
1114: 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[])
1115: {
1116:   const PetscInt Nc = uOff[2] - uOff[1];

1118:   for (PetscInt c = 0; c < Nc; ++c) {
1119:     g0[c * Nc + c]           = -1.0;
1120:     g0[Nc * Nc + c * Nc + c] = 1.0;
1121:   }
1122: }

1124: static PetscErrorCode TestAssembly(DM dm, AppCtx *user)
1125: {
1126:   Mat           J;
1127:   Vec           locX, locF, locW;
1128:   PetscDS       probh;
1129:   DMLabel       fault, material;
1130:   DM            dmFault;
1131:   IS            cohesiveCells;
1132:   PetscFE       fe;
1133:   PetscWeakForm wf;
1134:   PetscFormKey  keys[3];
1135:   PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx);
1136:   PetscInt    dim, Nf, cMax, cEnd, id;
1137:   PetscMPIInt rank;

1139:   PetscFunctionBegin;
1140:   if (!user->testAssembly) PetscFunctionReturn(PETSC_SUCCESS);
1141:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1142:   PetscCall(DMGetDimension(dm, &dim));
1143:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax));
1144:   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
1145:   PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells));
1146:   PetscCall(DMGetLabel(dm, "cohesive", &fault));
1147:   PetscCall(DMGetLocalVector(dm, &locX));
1148:   PetscCall(PetscObjectSetName((PetscObject)locX, "Local Solution"));
1149:   PetscCall(DMGetLocalVector(dm, &locF));
1150:   PetscCall(PetscObjectSetName((PetscObject)locF, "Local Residual"));
1151:   PetscCall(DMCreateMatrix(dm, &J));
1152:   PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));

1154:   /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */
1155:   PetscCall(DMGetLabel(dm, "material", &material));
1156:   id              = 1;
1157:   initialGuess[0] = r;
1158:   initialGuess[1] = NULL;
1159:   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1160:   id              = 2;
1161:   initialGuess[0] = rp1;
1162:   initialGuess[1] = NULL;
1163:   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1164:   id              = 1;
1165:   initialGuess[0] = NULL;
1166:   initialGuess[1] = phi;
1167:   PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX));
1168:   PetscCall(VecViewFromOptions(locX, NULL, "-local_solution_view"));

1170:   /* Test projection to fault mesh */
1171:   PetscCall(DMPlexCreateCohesiveSubmesh(dm, PETSC_FALSE, NULL, 0, &dmFault));
1172:   PetscCall(DMViewFromOptions(dmFault, NULL, "-fault_view"));
1173:   PetscCall(DMPlexOrient(dmFault));
1174:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, "fault_field_", PETSC_DETERMINE, &fe));
1175:   PetscCall(PetscFESetName(fe, "fault_field"));
1176:   PetscCall(DMAddField(dmFault, NULL, (PetscObject)fe));
1177:   PetscCall(PetscFEDestroy(&fe));
1178:   PetscCall(DMCreateDS(dmFault));
1179:   PetscCall(DMGetLocalVector(dmFault, &locW));
1180:   PetscCall(DMViewFromOptions(dmFault, NULL, "-cohesive_view"));
1181:   void (*faultFuncs[1])(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 f[]);

1183:   DMLabel  depthLabel;
1184:   PetscInt depth;
1185:   PetscCall(DMPlexGetDepthLabel(dmFault, &depthLabel));
1186:   PetscCall(DMPlexGetDepth(dmFault, &depth));
1187:   id = depth - 1;
1188:   /* w = r + rp1 */
1189:   faultFuncs[0] = add_fields;
1190:   PetscCall(DMProjectBdFieldLabelLocal(dmFault, 0.0, depthLabel, 1, &id, PETSC_DETERMINE, NULL, locX, faultFuncs, INSERT_VALUES, locW));
1191:   PetscCall(VecViewFromOptions(locW, NULL, "-local_projection_view"));

1193:   /* w = fault_normal */
1194:   faultFuncs[0] = normal_field;
1195:   PetscCall(DMProjectBdFieldLabelLocal(dmFault, 0.0, depthLabel, 1, &id, PETSC_DETERMINE, NULL, locX, faultFuncs, INSERT_VALUES, locW));
1196:   PetscCall(VecViewFromOptions(locW, NULL, "-local_projection_view"));
1197:   PetscCall(DMRestoreLocalVector(dmFault, &locW));
1198:   PetscCall(DMDestroy(&dmFault));

1200:   PetscCall(DMGetCellDS(dm, cMax, &probh, NULL));
1201:   PetscCall(PetscDSGetWeakForm(probh, &wf));
1202:   PetscCall(PetscDSGetNumFields(probh, &Nf));
1203:   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u_neg, 0, NULL));
1204:   PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u_pos, 0, NULL));
1205:   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul_neg, 0, NULL, 0, NULL, 0, NULL));
1206:   PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul_pos, 0, NULL, 0, NULL, 0, NULL));
1207:   if (Nf > 1) {
1208:     PetscCall(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL));
1209:     PetscCall(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL));
1210:   }
1211:   if (rank == 0) PetscCall(PetscDSView(probh, NULL));

1213:   keys[0].label = material;
1214:   keys[0].value = 1;
1215:   keys[0].field = 0;
1216:   keys[0].part  = 0;
1217:   keys[1].label = material;
1218:   keys[1].value = 2;
1219:   keys[1].field = 0;
1220:   keys[1].part  = 0;
1221:   keys[2].label = fault;
1222:   keys[2].value = 1;
1223:   keys[2].field = 1;
1224:   keys[2].part  = 0;
1225:   PetscCall(VecSet(locF, 0.));
1226:   PetscCall(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user));
1227:   PetscCall(VecViewFromOptions(locF, NULL, "-local_residual_view"));
1228:   PetscCall(MatZeroEntries(J));
1229:   PetscCall(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user));
1230:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1231:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1232:   PetscCall(MatViewFromOptions(J, NULL, "-local_jacobian_view"));

1234:   PetscCall(DMRestoreLocalVector(dm, &locX));
1235:   PetscCall(DMRestoreLocalVector(dm, &locF));
1236:   PetscCall(MatDestroy(&J));
1237:   PetscCall(ISDestroy(&cohesiveCells));

1239:   if (cMax < cEnd) {
1240:     PetscDS         ds;
1241:     PetscFE         fe;
1242:     PetscQuadrature quad;
1243:     IS             *perm;
1244:     const PetscInt *cone;
1245:     PetscInt        Na, a;

1247:     PetscCall(DMPlexGetCone(dm, cMax, &cone));
1248:     PetscCall(DMGetCellDS(dm, cMax, &ds, NULL));
1249:     PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
1250:     PetscCall(PetscFEGetQuadrature(fe, &quad));
1251:     PetscCall(PetscQuadratureComputePermutations(quad, &Na, &perm));
1252:     for (a = 0; a < Na; ++a) PetscCall(ISDestroy(&perm[a]));
1253:     PetscCall(PetscFree(perm));
1254:   }
1255:   PetscFunctionReturn(PETSC_SUCCESS);
1256: }

1258: int main(int argc, char **argv)
1259: {
1260:   DM     dm;
1261:   AppCtx user; /* user-defined work context */

1263:   PetscFunctionBeginUser;
1264:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1265:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1266:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1267:   PetscCall(TestMesh(dm, &user));
1268:   PetscCall(TestDiscretization(dm, &user));
1269:   PetscCall(TestAssembly(dm, &user));
1270:   PetscCall(DMDestroy(&dm));
1271:   PetscCall(PetscFinalize());
1272:   return 0;
1273: }

1275: /*TEST
1276:   testset:
1277:     args: -orig_dm_plex_check_all -dm_plex_check_all \
1278:           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \
1279:           -local_solution_view -local_residual_view -local_jacobian_view
1280:     test:
1281:       suffix: tri_0
1282:       args: -dim 2
1283:       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"
1284:     test:
1285:       suffix: tri_0_perm
1286:       args: -dim 2 -dm_reorder_section -dm_reorder_section_type cohesive
1287:       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" -e "s/_neg//g" -e "s/_pos//g" -e "s~_ZL.*~~g"
1288:     test:
1289:       suffix: tri_t1_0
1290:       args: -dim 2 -test_num 1
1291:       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"
1292:     test:
1293:       suffix: tri_t1_0_perm
1294:       args: -dim 2 -test_num 1 -dm_reorder_section -dm_reorder_section_type cohesive
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" -e "s/_neg//g" -e "s/_pos//g" -e "s~_ZL.*~~g"
1296:     test:
1297:       suffix: tri_t2_0
1298:       args: -dim 2 -test_num 2
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: tri_t5_0
1303:       args: -dim 2 -test_num 5
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: tet_0
1308:       args: -dim 3
1309:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1311:     test:
1312:       suffix: tet_t1_0
1313:       args: -dim 3 -test_num 1
1314:       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"

1316:   testset:
1317:     args: -orig_dm_plex_check_all -dm_plex_check_all \
1318:           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1319:     test:
1320:       suffix: tet_1
1321:       nsize: 2
1322:       args: -dim 3
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"

1325:     test:
1326:       suffix: tri_1
1327:       nsize: 2
1328:       args: -dim 2
1329:       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"

1331:     test:
1332:       suffix: tri_t3_0
1333:       nsize: 2
1334:       args: -dim 2 -test_num 3
1335:       filter: sed -e "s/_start//g" -e "s/f0_bd_u_neg//g" -e "s/f0_bd_u_pos//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul_neg//g" -e "s/g0_bd_ul_pos//g" -e "s/g0_bd_lu//g" -e "s~_ZL.*~~g"

1337:     test:
1338:       suffix: tri_t4_0
1339:       nsize: 6
1340:       args: -dim 2 -test_num 4
1341:       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"

1343:   testset:
1344:     args: -orig_dm_plex_check_all -dm_plex_check_all \
1345:           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1346:     # 2D Quads
1347:     test:
1348:       suffix: quad_0
1349:       args: -dim 2 -cell_simplex 0
1350:       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"

1352:     test:
1353:       suffix: quad_1
1354:       nsize: 2
1355:       args: -dim 2 -cell_simplex 0
1356:       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"

1358:     # Caanot run the assembly test because we cannot orient a fault mesh in a T-shape
1359:     test:
1360:       suffix: quad_t1_0
1361:       args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all -test_assembly 0
1362:       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"

1364:     test:
1365:       suffix: quad_t2_0
1366:       nsize: 2
1367:       args: -dim 2 -cell_simplex 0 -test_num 2
1368:       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"

1370:     test:
1371:       # TODO: The PetscSF is wrong here (connects to wrong side of split)
1372:       suffix: quad_t3_0
1373:       nsize: 2
1374:       args: -dim 2 -cell_simplex 0 -test_num 3
1375:       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"

1377:     # 3D Hex
1378:     test:
1379:       suffix: hex_0
1380:       args: -dim 3 -cell_simplex 0
1381:       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"
1382:     test:
1383:       suffix: hex_1
1384:       nsize: 2
1385:       args: -dim 3 -cell_simplex 0
1386:       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"
1387:     test:
1388:       suffix: hex_t1_0
1389:       args: -dim 3 -cell_simplex 0 -test_num 1
1390:       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"
1391:     test:
1392:       suffix: hex_t2_0
1393:       args: -dim 3 -cell_simplex 0 -test_num 2
1394:       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"

1396: TEST*/