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