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:   options->debug          = 0;
404:   options->dim            = 2;
405:   options->cellSimplex    = PETSC_TRUE;
406:   options->testPartition  = PETSC_TRUE;
407:   options->testNum        = 0;
408:   options->cohesiveFields = 1;

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

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

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

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

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

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

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

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

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

529:   MPI_Comm_rank(comm, &rank);
530:   if (rank == 0) {
531:     switch (testNum) {
532:     case 0: {
533:       PetscInt    numPoints[4]         = {5, 7, 9, 2};
534:       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};
535:       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};
536:       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};
537:       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};
538:       PetscInt    markerPoints[20]     = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1};
539:       PetscInt    faultPoints[3]       = {3, 4, 5};

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

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

570:     DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);
571:     DMCreateLabel(dm, "fault");
572:   }
573:   return 0;
574: }

576: static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
577: {
578:   DM          idm;
579:   PetscInt    p;
580:   PetscMPIInt rank;

582:   MPI_Comm_rank(comm, &rank);
583:   if (rank == 0) {
584:     switch (testNum) {
585:     case 0:
586:     case 2:
587:     case 3: {
588:       PetscInt    numPoints[2]        = {6, 2};
589:       PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
590:       PetscInt    cones[8]            = {2, 3, 4, 5, 3, 6, 7, 4};
591:       PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
592:       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};
593:       PetscInt    markerPoints[12]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1};
594:       PetscInt    faultPoints[2]      = {3, 4};

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

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

634:     DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);
635:     if (testNum == 2 || testNum == 3) DMCreateLabel(*dm, "pfault");
636:     else DMCreateLabel(*dm, "fault");
637:   }
638:   DMPlexInterpolate(*dm, &idm);
639:   DMViewFromOptions(idm, NULL, "-in_dm_view");
640:   DMDestroy(dm);
641:   *dm = idm;
642:   return 0;
643: }

645: static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
646: {
647:   DM          idm;
648:   PetscInt    depth = 3, p;
649:   PetscMPIInt rank;

651:   MPI_Comm_rank(comm, &rank);
652:   if (rank == 0) {
653:     switch (testNum) {
654:     case 0: {
655:       PetscInt    numPoints[2]         = {12, 2};
656:       PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
657:       PetscInt    cones[16]            = {2, 5, 4, 3, 6, 7, 8, 9, 3, 4, 11, 10, 7, 12, 13, 8};
658:       PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
659:       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};
660:       PetscInt    markerPoints[52]     = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1, 8, 1, 9, 1};
661:       PetscInt    faultPoints[4]       = {3, 4, 7, 8};

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

689:       DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
690:       DMPlexInterpolate(*dm, &idm);
691:       for (p = 0; p < 6; ++p) DMSetLabelValue(idm, "fault", faultPoints[p], 1);
692:       DMSetLabelValue(*dm, "material", 0, 1);
693:       DMSetLabelValue(*dm, "material", 1, 1);
694:       DMSetLabelValue(*dm, "material", 2, 1);
695:       DMSetLabelValue(*dm, "material", 3, 2);
696:       DMSetLabelValue(*dm, "material", 4, 2);
697:       DMSetLabelValue(*dm, "material", 5, 2);
698:       DMSetLabelValue(*dm, "material", 6, 2);
699:     } break;
700:     case 2: {
701:       /* Buried fault edge */
702:       PetscInt    numPoints[2]         = {18, 4};
703:       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};
704:       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};
705:       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};
706:       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,
707:                                           -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};
708:       PetscInt    faultPoints[4]       = {7, 8, 16, 17};

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

724:     DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL);
725:     DMPlexInterpolate(*dm, &idm);
726:     DMCreateLabel(idm, "fault");
727:   }
728:   DMViewFromOptions(idm, NULL, "-in_dm_view");
729:   DMDestroy(dm);
730:   *dm = idm;
731:   return 0;
732: }

734: static PetscErrorCode CreateFaultLabel(DM dm)
735: {
736:   DMLabel  label;
737:   PetscInt dim, h, pStart, pEnd, pMax, p;

739:   DMGetDimension(dm, &dim);
740:   DMCreateLabel(dm, "cohesive");
741:   DMGetLabel(dm, "cohesive", &label);
742:   for (h = 0; h <= dim; ++h) {
743:     DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax);
744:     DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
745:     for (p = pMax; p < pEnd; ++p) DMLabelSetValue(label, p, 1);
746:   }
747:   return 0;
748: }

750: static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user)
751: {
752:   PetscFE  fe;
753:   DMLabel  fault;
754:   PetscInt dim, Ncf = user->cohesiveFields, f;

756:   DMGetDimension(dm, &dim);
757:   DMGetLabel(dm, "cohesive", &fault);
758:   DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD);

760:   PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe);
761:   PetscFESetName(fe, "displacement");
762:   DMAddField(dm, NULL, (PetscObject)fe);
763:   PetscFEDestroy(&fe);

765:   if (Ncf > 0) {
766:     PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe);
767:     PetscFESetName(fe, "fault traction");
768:     DMAddField(dm, fault, (PetscObject)fe);
769:     PetscFEDestroy(&fe);
770:   }
771:   for (f = 1; f < Ncf; ++f) {
772:     char name[256], opt[256];

774:     PetscSNPrintf(name, 256, "fault field %" PetscInt_FMT, f);
775:     PetscSNPrintf(opt, 256, "faultfield_%" PetscInt_FMT "_", f);
776:     PetscFECreateDefault(PETSC_COMM_SELF, dim - 1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe);
777:     PetscFESetName(fe, name);
778:     DMAddField(dm, fault, (PetscObject)fe);
779:     PetscFEDestroy(&fe);
780:   }

782:   DMCreateDS(dm);
783:   return 0;
784: }

786: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
787: {
788:   PetscInt    dim         = user->dim;
789:   PetscBool   cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault;
790:   PetscMPIInt rank, size;
791:   DMLabel     matLabel;

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

827:     DMGetLabel(*dm, "fault", &faultLabel);
828:     DMGetLabel(*dm, "faultBd", &faultBdLabel);
829:     DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid);
830:     DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);
831:     DMLabelDestroy(&hybridLabel);
832:     DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD);
833:     DMLabelDestroy(&splitLabel);
834:     DMViewFromOptions(dmInterface, NULL, "-dm_interface_view");
835:     DMDestroy(&dmInterface);
836:     DMDestroy(dm);
837:     *dm = dmHybrid;
838:   }
839:   DMHasLabel(*dm, "fault2", &hasFault2);
840:   if (hasFault2) {
841:     DM      dmHybrid = NULL;
842:     DMLabel faultLabel, faultBdLabel, hybridLabel;

844:     PetscObjectSetOptionsPrefix((PetscObject)*dm, "faulted_");
845:     DMViewFromOptions(*dm, NULL, "-dm_view_pre");
846:     DMPlexDistributeSetDefault(*dm, PETSC_FALSE);
847:     DMSetFromOptions(*dm);
848:     DMViewFromOptions(*dm, NULL, "-dm_view");
849:     DMGetLabel(*dm, "fault2", &faultLabel);
850:     DMGetLabel(*dm, "fault2Bd", &faultBdLabel);
851:     DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, NULL, &dmHybrid);
852:     DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);
853:     DMLabelDestroy(&hybridLabel);
854:     DMDestroy(dm);
855:     *dm = dmHybrid;
856:   }
857:   if (user->testPartition && size > 1) {
858:     PetscPartitioner part;
859:     PetscInt        *sizes  = NULL;
860:     PetscInt        *points = NULL;

862:     if (rank == 0) {
863:       if (dim == 2 && cellSimplex && size == 2) {
864:         switch (user->testNum) {
865:         case 0: {
866:           PetscInt triSizes_p2[2]  = {1, 2};
867:           PetscInt triPoints_p2[3] = {0, 1, 2};

869:           PetscMalloc2(2, &sizes, 3, &points);
870:           PetscArraycpy(sizes, triSizes_p2, 2);
871:           PetscArraycpy(points, triPoints_p2, 3);
872:           break;
873:         }
874:         case 3: {
875:           PetscInt triSizes_p2[2]  = {3, 3};
876:           PetscInt triPoints_p2[6] = {0, 1, 2, 3, 4, 5};

878:           PetscMalloc2(2, &sizes, 6, &points);
879:           PetscArraycpy(sizes, triSizes_p2, 2);
880:           PetscArraycpy(points, triPoints_p2, 6);
881:           break;
882:         }
883:         default:
884:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 2 procs", user->testNum);
885:         }
886:       } else if (dim == 2 && cellSimplex && size == 6) {
887:         switch (user->testNum) {
888:         case 4: {
889:           PetscInt triSizes_p6[6]  = {1, 1, 1, 1, 1, 1};
890:           PetscInt triPoints_p6[6] = {0, 1, 2, 3, 4, 5};

892:           PetscMalloc2(6, &sizes, 6, &points);
893:           PetscArraycpy(sizes, triSizes_p6, 6);
894:           PetscArraycpy(points, triPoints_p6, 6);
895:           break;
896:         }
897:         default:
898:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for triangular mesh on 6 procs", user->testNum);
899:         }
900:       } else if (dim == 2 && !cellSimplex && size == 2) {
901:         switch (user->testNum) {
902:         case 0: {
903:           PetscInt quadSizes_p2[2]  = {1, 2};
904:           PetscInt quadPoints_p2[3] = {0, 1, 2};

906:           PetscMalloc2(2, &sizes, 3, &points);
907:           PetscArraycpy(sizes, quadSizes_p2, 2);
908:           PetscArraycpy(points, quadPoints_p2, 3);
909:           break;
910:         }
911:         case 2: {
912:           PetscInt quadSizes_p2[2]  = {1, 1};
913:           PetscInt quadPoints_p2[2] = {0, 1};

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

924:           PetscMalloc2(2, &sizes, 2, &points);
925:           PetscArraycpy(sizes, quadSizes_p2, 2);
926:           PetscArraycpy(points, quadPoints_p2, 2);
927:           break;
928:         }
929:         default:
930:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for quadrilateral mesh on 2 procs", user->testNum);
931:         }
932:       } else if (dim == 3 && cellSimplex && size == 2) {
933:         switch (user->testNum) {
934:         case 0: {
935:           PetscInt tetSizes_p2[2]  = {1, 2};
936:           PetscInt tetPoints_p2[3] = {0, 1, 2};

938:           PetscMalloc2(2, &sizes, 3, &points);
939:           PetscArraycpy(sizes, tetSizes_p2, 2);
940:           PetscArraycpy(points, tetPoints_p2, 3);
941:           break;
942:         }
943:         default:
944:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for teterehedral mesh on 2 procs", user->testNum);
945:         }
946:       } else if (dim == 3 && !cellSimplex && size == 2) {
947:         switch (user->testNum) {
948:         case 0: {
949:           PetscInt hexSizes_p2[2]  = {1, 2};
950:           PetscInt hexPoints_p2[3] = {0, 1, 2};

952:           PetscMalloc2(2, &sizes, 3, &points);
953:           PetscArraycpy(sizes, hexSizes_p2, 2);
954:           PetscArraycpy(points, hexPoints_p2, 3);
955:           break;
956:         }
957:         default:
958:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %" PetscInt_FMT " for hexahedral mesh on 2 procs", user->testNum);
959:         }
960:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition");
961:     }
962:     DMPlexGetPartitioner(*dm, &part);
963:     PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);
964:     PetscPartitionerShellSetPartition(part, size, sizes, points);
965:     PetscFree2(sizes, points);
966:   }
967:   {
968:     DM pdm = NULL;

970:     /* Distribute mesh over processes */
971:     DMPlexDistribute(*dm, 0, NULL, &pdm);
972:     if (pdm) {
973:       DMViewFromOptions(pdm, NULL, "-dm_view");
974:       DMDestroy(dm);
975:       *dm = pdm;
976:     }
977:   }
978:   DMHasLabel(*dm, "pfault", &hasParallelFault);
979:   if (hasParallelFault) {
980:     DM      dmHybrid = NULL, dmInterface;
981:     DMLabel faultLabel, faultBdLabel, hybridLabel;

983:     DMGetLabel(*dm, "pfault", &faultLabel);
984:     DMGetLabel(*dm, "pfaultBd", &faultBdLabel);
985:     DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, 1, &hybridLabel, NULL, &dmInterface, &dmHybrid);
986:     DMViewFromOptions(dmInterface, NULL, "-dm_fault_view");
987:     {
988:       PetscViewer viewer;
989:       PetscMPIInt rank;

991:       MPI_Comm_rank(PetscObjectComm((PetscObject)*dm), &rank);
992:       PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer);
993:       PetscViewerASCIIPrintf(viewer, "Rank %d\n", rank);
994:       DMLabelView(hybridLabel, viewer);
995:       PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &viewer);
996:       PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
997:     }
998:     DMLabelDestroy(&hybridLabel);
999:     DMDestroy(&dmInterface);
1000:     DMDestroy(dm);
1001:     *dm = dmHybrid;
1002:   }
1003:   PetscObjectSetName((PetscObject)*dm, "Hybrid Mesh");
1004:   CreateFaultLabel(*dm);
1005:   CreateDiscretization(*dm, user);
1006:   DMViewFromOptions(*dm, NULL, "-dm_view_pre");
1007:   DMPlexDistributeSetDefault(*dm, PETSC_FALSE);
1008:   DMSetFromOptions(*dm);
1009:   DMViewFromOptions(*dm, NULL, "-dm_view");
1010:   return 0;
1011: }

1013: static PetscErrorCode TestMesh(DM dm, AppCtx *user)
1014: {
1015:   DMPlexCheckSymmetry(dm);
1016:   DMPlexCheckSkeleton(dm, 0);
1017:   DMPlexCheckFaces(dm, 0);
1018:   return 0;
1019: }

1021: static PetscErrorCode TestDiscretization(DM dm, AppCtx *user)
1022: {
1023:   PetscSection s;

1025:   DMGetSection(dm, &s);
1026:   PetscObjectViewFromOptions((PetscObject)s, NULL, "-local_section_view");
1027:   return 0;
1028: }

1030: static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1031: {
1032:   PetscInt d;
1033:   for (d = 0; d < dim; ++d) u[d] = x[d];
1034:   return 0;
1035: }

1037: static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1038: {
1039:   PetscInt d;
1040:   for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0);
1041:   return 0;
1042: }

1044: static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1045: {
1046:   PetscInt d;
1047:   u[0] = -x[1];
1048:   u[1] = x[0];
1049:   for (d = 2; d < dim; ++d) u[d] = x[d];
1050:   return 0;
1051: }

1053: /* \lambda \cdot (\psi_u^- - \psi_u^+) */
1054: static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1055: {
1056:   const PetscInt Nc = uOff[1] - uOff[0];
1057:   PetscInt       c;
1058:   for (c = 0; c < Nc; ++c) {
1059:     f0[c]      = u[Nc * 2 + c] + x[Nc - c - 1];
1060:     f0[Nc + c] = -(u[Nc * 2 + c] + x[Nc - c - 1]);
1061:   }
1062: }

1064: /* (d - u^+ + u^-) \cdot \psi_\lambda */
1065: 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[])
1066: {
1067:   const PetscInt Nc = uOff[2] - uOff[1];
1068:   PetscInt       c;

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

1073: /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */
1074: static void g0_bd_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1075: {
1076:   const PetscInt Nc = uOff[1] - uOff[0];
1077:   PetscInt       c;

1079:   for (c = 0; c < Nc; ++c) {
1080:     g0[(0 + c) * Nc + c]  = 1.0;
1081:     g0[(Nc + c) * Nc + c] = -1.0;
1082:   }
1083: }

1085: /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */
1086: 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[])
1087: {
1088:   const PetscInt Nc = uOff[2] - uOff[1];
1089:   PetscInt       c;

1091:   for (c = 0; c < Nc; ++c) {
1092:     g0[c * Nc * 2 + c]      = 1.0;
1093:     g0[c * Nc * 2 + Nc + c] = -1.0;
1094:   }
1095: }

1097: static PetscErrorCode TestAssembly(DM dm, AppCtx *user)
1098: {
1099:   Mat           J;
1100:   Vec           locX, locF;
1101:   PetscDS       probh;
1102:   DMLabel       fault, material;
1103:   IS            cohesiveCells;
1104:   PetscWeakForm wf;
1105:   PetscFormKey  keys[3];
1106:   PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx);
1107:   PetscInt    dim, Nf, cMax, cEnd, id;
1108:   PetscMPIInt rank;

1110:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1111:   DMGetDimension(dm, &dim);
1112:   DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax);
1113:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1114:   ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells);
1115:   DMGetLabel(dm, "cohesive", &fault);
1116:   DMGetLocalVector(dm, &locX);
1117:   PetscObjectSetName((PetscObject)locX, "Local Solution");
1118:   DMGetLocalVector(dm, &locF);
1119:   PetscObjectSetName((PetscObject)locF, "Local Residual");
1120:   DMCreateMatrix(dm, &J);
1121:   PetscObjectSetName((PetscObject)J, "Jacobian");

1123:   /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */
1124:   DMGetLabel(dm, "material", &material);
1125:   id              = 1;
1126:   initialGuess[0] = r;
1127:   initialGuess[1] = NULL;
1128:   DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);
1129:   id              = 2;
1130:   initialGuess[0] = rp1;
1131:   initialGuess[1] = NULL;
1132:   DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);
1133:   id              = 1;
1134:   initialGuess[0] = NULL;
1135:   initialGuess[1] = phi;
1136:   DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);
1137:   VecViewFromOptions(locX, NULL, "-local_solution_view");

1139:   DMGetCellDS(dm, cMax, &probh);
1140:   PetscDSGetWeakForm(probh, &wf);
1141:   PetscDSGetNumFields(probh, &Nf);
1142:   PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u, 0, NULL);
1143:   PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u, 0, NULL);
1144:   PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL);
1145:   PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL);
1146:   if (Nf > 1) {
1147:     PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL);
1148:     PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL);
1149:   }
1150:   if (rank == 0) PetscDSView(probh, NULL);

1152:   keys[0].label = NULL;
1153:   keys[0].value = 0;
1154:   keys[0].field = 0;
1155:   keys[0].part  = 0;
1156:   keys[1].label = material;
1157:   keys[1].value = 2;
1158:   keys[1].field = 0;
1159:   keys[1].part  = 0;
1160:   keys[2].label = fault;
1161:   keys[2].value = 1;
1162:   keys[2].field = 1;
1163:   keys[2].part  = 0;
1164:   VecSet(locF, 0.);
1165:   DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user);
1166:   VecViewFromOptions(locF, NULL, "-local_residual_view");
1167:   MatZeroEntries(J);
1168:   DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user);
1169:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
1170:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
1171:   MatViewFromOptions(J, NULL, "-local_jacobian_view");

1173:   DMRestoreLocalVector(dm, &locX);
1174:   DMRestoreLocalVector(dm, &locF);
1175:   MatDestroy(&J);
1176:   ISDestroy(&cohesiveCells);
1177:   return 0;
1178: }

1180: int main(int argc, char **argv)
1181: {
1182:   DM     dm;
1183:   AppCtx user; /* user-defined work context */

1186:   PetscInitialize(&argc, &argv, NULL, help);
1187:   ProcessOptions(PETSC_COMM_WORLD, &user);
1188:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
1189:   TestMesh(dm, &user);
1190:   TestDiscretization(dm, &user);
1191:   TestAssembly(dm, &user);
1192:   DMDestroy(&dm);
1193:   PetscFinalize();
1194:   return 0;
1195: }

1197: /*TEST
1198:   testset:
1199:     args: -orig_dm_plex_check_all -dm_plex_check_all \
1200:           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \
1201:           -local_solution_view -local_residual_view -local_jacobian_view
1202:     test:
1203:       suffix: tri_0
1204:       args: -dim 2
1205:       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"
1206:     test:
1207:       suffix: tri_t1_0
1208:       args: -dim 2 -test_num 1
1209:       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"
1210:     test:
1211:       suffix: tri_t2_0
1212:       args: -dim 2 -test_num 2
1213:       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"
1214:     test:
1215:       suffix: tri_t5_0
1216:       args: -dim 2 -test_num 5
1217:       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"
1218:     test:
1219:       suffix: tet_0
1220:       args: -dim 3
1221:       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1222:     test:
1223:       suffix: tet_t1_0
1224:       args: -dim 3 -test_num 1
1225:       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"

1227:   testset:
1228:     args: -orig_dm_plex_check_all -dm_plex_check_all \
1229:           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1230:     test:
1231:       suffix: tet_1
1232:       nsize: 2
1233:       args: -dim 3
1234:       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"
1235:     test:
1236:       suffix: tri_1
1237:       nsize: 2
1238:       args: -dim 2
1239:       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"
1240:     test:
1241:       suffix: tri_t3_0
1242:       nsize: 2
1243:       args: -dim 2 -test_num 3
1244:       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"
1245:     test:
1246:       suffix: tri_t4_0
1247:       nsize: 6
1248:       args: -dim 2 -test_num 4
1249:       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"

1251:   testset:
1252:     args: -orig_dm_plex_check_all -dm_plex_check_all \
1253:           -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1
1254:     # 2D Quads
1255:     test:
1256:       suffix: quad_0
1257:       args: -dim 2 -cell_simplex 0
1258:       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"
1259:     test:
1260:       suffix: quad_1
1261:       nsize: 2
1262:       args: -dim 2 -cell_simplex 0
1263:       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1264:     test:
1265:       suffix: quad_t1_0
1266:       args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all
1267:       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1268:     test:
1269:       suffix: quad_t2_0
1270:       nsize: 2
1271:       args: -dim 2 -cell_simplex 0 -test_num 2
1272:       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1273:     test:
1274:       # TODO: The PetscSF is wrong here (connects to wrong side of split)
1275:       suffix: quad_t3_0
1276:       nsize: 2
1277:       args: -dim 2 -cell_simplex 0 -test_num 3
1278:       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"
1279:     # 3D Hex
1280:     test:
1281:       suffix: hex_0
1282:       args: -dim 3 -cell_simplex 0
1283:       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"
1284:     test:
1285:       suffix: hex_1
1286:       nsize: 2
1287:       args: -dim 3 -cell_simplex 0
1288:       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"
1289:     test:
1290:       suffix: hex_t1_0
1291:       args: -dim 3 -cell_simplex 0 -test_num 1
1292:       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"
1293:     test:
1294:       suffix: hex_t2_0
1295:       args: -dim 3 -cell_simplex 0 -test_num 2
1296:       filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g"

1298: TEST*/