Actual source code: plextrcohesive.c
1: #include <petsc/private/dmplextransformimpl.h>
3: #include <petsc/private/dmlabelimpl.h>
5: /*
6: The cohesive transformation extrudes cells into a mesh from faces along an internal boundary.
8: Orientation:
10: We will say that a face has a positive and negative side. The positive side is defined by the cell which attaches the face with a positive orientation, and the negative side cell attaches it with a negative orientation (a reflection). However, this means that the positive side is in the opposite direction of the face normal, and the negative side is in the direction of the face normal, since all cells have outward facing normals. For clarity, in 2D the cross product of the normal and the edge is in the positive z direction.
12: Labeling:
14: We require an active label on input, which marks all points on the internal surface. Each point is
15: labeled with its depth. This label is passed to DMPlexLabelCohesiveComplete(), which adds all points
16: which ``impinge'' on the surface, meaning a point has a face on the surface. These points are labeled
17: with celltype + 100 on the positive side, and -(celltype + 100) on the negative side.
19: Point Creation:
21: We split points on the fault surface, creating a new partner point for each one. The negative side
22: receives the old point, while the positive side receives the new partner. In addition, points are
23: created with the two split points as boundaries. For example, split vertices have a segment between
24: them, split edges a quadrilaterial, split triangles a prism, and split quads a hexahedron. By
25: default, these spanning points have tensor ordering, but the user can choose to have them use the
26: outward normal convention instead.
28: */
30: static PetscErrorCode DMPlexTransformView_Cohesive(DMPlexTransform tr, PetscViewer viewer)
31: {
32: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
33: PetscBool isascii;
35: PetscFunctionBegin;
38: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
39: if (isascii) {
40: const char *name;
42: PetscCall(PetscObjectGetName((PetscObject)tr, &name));
43: PetscCall(PetscViewerASCIIPrintf(viewer, "Cohesive extrusion transformation %s\n", name ? name : ""));
44: PetscCall(PetscViewerASCIIPrintf(viewer, " create tensor cells: %s\n", ex->useTensor ? "YES" : "NO"));
45: } else {
46: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
47: }
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode DMPlexTransformSetFromOptions_Cohesive(DMPlexTransform tr, PetscOptionItems PetscOptionsObject)
52: {
53: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
54: PetscReal width;
55: PetscBool tensor, flg;
57: PetscFunctionBegin;
58: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Cohesive Extrusion Options");
59: PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg));
60: if (flg) PetscCall(DMPlexTransformCohesiveExtrudeSetTensor(tr, tensor));
61: PetscCall(PetscOptionsReal("-dm_plex_transform_cohesive_width", "Width of a cohesive cell", "", ex->width, &width, &flg));
62: if (flg) PetscCall(DMPlexTransformCohesiveExtrudeSetWidth(tr, width));
63: PetscCall(PetscOptionsInt("-dm_plex_transform_cohesive_debug", "Det debugging level", "", ex->debug, &ex->debug, NULL));
64: PetscOptionsHeadEnd();
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: /*
69: ComputeSplitFaceNumber - Compute an encoding describing which faces of p are split by the surface
71: Not collective
73: Input Parameters:
74: + dm - The `DM`
75: . label - `DMLabel` marking the surface and adjacent points
76: - p - Impinging point, adjacent to the surface
78: Output Parameter:
79: . fsplit - A number encoding the faces which are split by the surface
81: Level: developer
83: Note: We will use a bit encoding, where bit k is 1 if face k is split.
85: .seealso: ComputeUnsplitFaceNumber()
86: */
87: static PetscErrorCode ComputeSplitFaceNumber(DM dm, DMLabel label, PetscInt p, PetscInt *fsplit)
88: {
89: const PetscInt *cone;
90: PetscInt coneSize, val;
92: PetscFunctionBegin;
93: *fsplit = 0;
94: PetscCall(DMPlexGetCone(dm, p, &cone));
95: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
96: PetscCheck(coneSize < (PetscInt)sizeof(*fsplit) * 8, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cone of size %" PetscInt_FMT " is too large to be contained in an integer", coneSize);
97: for (PetscInt c = 0; c < coneSize; ++c) {
98: PetscCall(DMLabelGetValue(label, cone[c], &val));
99: if (val >= 0 && val < 100) *fsplit |= 1 << c;
100: }
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /*
105: ComputeUnsplitFaceNumber - Compute an encoding describing which faces of p are unsplit by the surface
107: Not collective
109: Input Parameters:
110: + dm - The `DM`
111: . label - `DMLabel` marking the surface and adjacent points
112: - p - Split point, on the surface
114: Output Parameter:
115: . funsplit - A number encoding the faces which are split by the surface
117: Level: developer
119: Note: We will use a bit encoding, where bit k is 1 if face k is unsplit.
121: .seealso: ComputeSplitFaceNumber()
122: */
123: static PetscErrorCode ComputeUnsplitFaceNumber(DM dm, DMLabel label, PetscInt p, PetscInt *funsplit)
124: {
125: const PetscInt *cone;
126: PetscInt coneSize, val;
128: PetscFunctionBegin;
129: *funsplit = 0;
130: PetscCall(DMPlexGetCone(dm, p, &cone));
131: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
132: PetscCheck(coneSize < (PetscInt)sizeof(*funsplit) * 8, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cone of size %" PetscInt_FMT " is too large to be contained in an integer", coneSize);
133: for (PetscInt c = 0; c < coneSize; ++c) {
134: PetscCall(DMLabelGetValue(label, cone[c], &val));
135: if (val >= 200) *funsplit |= 1 << c;
136: }
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /* DM_POLYTOPE_POINT produces
141: 2 points when split, or 1 point when unsplit, and
142: 1 segment, or tensor segment
143: */
144: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Point(DMPlexTransform_Cohesive *ex)
145: {
146: PetscInt rt, Nc, No;
148: PetscFunctionBegin;
149: // Unsplit vertex
150: rt = DM_POLYTOPE_POINT * 2 + 1;
151: ex->Nt[rt] = 2;
152: Nc = 6;
153: No = 2;
154: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
155: ex->target[rt][0] = DM_POLYTOPE_POINT;
156: ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
157: ex->size[rt][0] = 1;
158: ex->size[rt][1] = 1;
159: // cone for segment/tensor segment
160: ex->cone[rt][0] = DM_POLYTOPE_POINT;
161: ex->cone[rt][1] = 0;
162: ex->cone[rt][2] = 0;
163: ex->cone[rt][3] = DM_POLYTOPE_POINT;
164: ex->cone[rt][4] = 0;
165: ex->cone[rt][5] = 0;
166: for (PetscInt i = 0; i < No; ++i) ex->ornt[rt][i] = 0;
167: // Split vertex
168: rt = (DM_POLYTOPE_POINT * 2 + 1) * 100 + 0;
169: ex->Nt[rt] = 2;
170: Nc = 6;
171: No = 2;
172: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
173: ex->target[rt][0] = DM_POLYTOPE_POINT;
174: ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
175: ex->size[rt][0] = 2;
176: ex->size[rt][1] = 1;
177: // cone for segment/tensor segment
178: ex->cone[rt][0] = DM_POLYTOPE_POINT;
179: ex->cone[rt][1] = 0;
180: ex->cone[rt][2] = 0;
181: ex->cone[rt][3] = DM_POLYTOPE_POINT;
182: ex->cone[rt][4] = 0;
183: ex->cone[rt][5] = 1;
184: for (PetscInt i = 0; i < No; ++i) ex->ornt[rt][i] = 0;
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: /* DM_POLYTOPE_SEGMENT produces
189: 2 segments when split, or 1 segment when unsplit, and
190: 1 quad, or tensor quad
191: */
192: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Segment(DMPlexTransform_Cohesive *ex)
193: {
194: PetscInt rt, Nc, No, coff, ooff;
196: PetscFunctionBegin;
197: // Unsplit segment
198: rt = DM_POLYTOPE_SEGMENT * 2 + 1;
199: ex->Nt[rt] = 2;
200: Nc = 8 + 14;
201: No = 2 + 4;
202: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
203: ex->target[rt][0] = DM_POLYTOPE_SEGMENT;
204: ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
205: ex->size[rt][0] = 1;
206: ex->size[rt][1] = 1;
207: // cones for segment
208: ex->cone[rt][0] = DM_POLYTOPE_POINT;
209: ex->cone[rt][1] = 1;
210: ex->cone[rt][2] = 0;
211: ex->cone[rt][3] = 0;
212: ex->cone[rt][4] = DM_POLYTOPE_POINT;
213: ex->cone[rt][5] = 1;
214: ex->cone[rt][6] = 1;
215: ex->cone[rt][7] = 0;
216: for (PetscInt i = 0; i < 2; ++i) ex->ornt[rt][i] = 0;
217: // cone for quad/tensor quad
218: coff = 8;
219: ooff = 2;
220: if (ex->useTensor) {
221: ex->cone[rt][coff + 0] = DM_POLYTOPE_SEGMENT;
222: ex->cone[rt][coff + 1] = 0;
223: ex->cone[rt][coff + 2] = 0;
224: ex->cone[rt][coff + 3] = DM_POLYTOPE_SEGMENT;
225: ex->cone[rt][coff + 4] = 0;
226: ex->cone[rt][coff + 5] = 0;
227: ex->cone[rt][coff + 6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
228: ex->cone[rt][coff + 7] = 1;
229: ex->cone[rt][coff + 8] = 0;
230: ex->cone[rt][coff + 9] = 0;
231: ex->cone[rt][coff + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
232: ex->cone[rt][coff + 11] = 1;
233: ex->cone[rt][coff + 12] = 1;
234: ex->cone[rt][coff + 13] = 0;
235: ex->ornt[rt][ooff + 0] = 0;
236: ex->ornt[rt][ooff + 1] = 0;
237: ex->ornt[rt][ooff + 2] = 0;
238: ex->ornt[rt][ooff + 3] = 0;
239: } else {
240: ex->cone[rt][coff + 0] = DM_POLYTOPE_SEGMENT;
241: ex->cone[rt][coff + 1] = 0;
242: ex->cone[rt][coff + 2] = 0;
243: ex->cone[rt][coff + 3] = DM_POLYTOPE_SEGMENT;
244: ex->cone[rt][coff + 4] = 1;
245: ex->cone[rt][coff + 5] = 1;
246: ex->cone[rt][coff + 6] = 0;
247: ex->cone[rt][coff + 7] = DM_POLYTOPE_SEGMENT;
248: ex->cone[rt][coff + 8] = 0;
249: ex->cone[rt][coff + 9] = 0;
250: ex->cone[rt][coff + 10] = DM_POLYTOPE_SEGMENT;
251: ex->cone[rt][coff + 11] = 1;
252: ex->cone[rt][coff + 12] = 0;
253: ex->cone[rt][coff + 13] = 0;
254: ex->ornt[rt][ooff + 0] = 0;
255: ex->ornt[rt][ooff + 1] = 0;
256: ex->ornt[rt][ooff + 2] = -1;
257: ex->ornt[rt][ooff + 3] = -1;
258: }
259: // Split segment
260: // 0: no unsplit vertex
261: // 1: unsplit vertex 0
262: // 2: unsplit vertex 1
263: // 3: both vertices unsplit (impossible)
264: for (PetscInt s = 0; s < 3; ++s) {
265: rt = (DM_POLYTOPE_SEGMENT * 2 + 1) * 100 + s;
266: ex->Nt[rt] = 2;
267: Nc = 8 * 2 + 14;
268: No = 2 * 2 + 4;
269: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
270: ex->target[rt][0] = DM_POLYTOPE_SEGMENT;
271: ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
272: ex->size[rt][0] = 2;
273: ex->size[rt][1] = 1;
274: // cones for segments
275: for (PetscInt i = 0; i < 2; ++i) {
276: ex->cone[rt][8 * i + 0] = DM_POLYTOPE_POINT;
277: ex->cone[rt][8 * i + 1] = 1;
278: ex->cone[rt][8 * i + 2] = 0;
279: ex->cone[rt][8 * i + 3] = s == 1 ? 0 : i;
280: ex->cone[rt][8 * i + 4] = DM_POLYTOPE_POINT;
281: ex->cone[rt][8 * i + 5] = 1;
282: ex->cone[rt][8 * i + 6] = 1;
283: ex->cone[rt][8 * i + 7] = s == 2 ? 0 : i;
284: }
285: for (PetscInt i = 0; i < 2 * 2; ++i) ex->ornt[rt][i] = 0;
286: // cone for quad/tensor quad
287: coff = 8 * 2;
288: ooff = 2 * 2;
289: if (ex->useTensor) {
290: ex->cone[rt][coff + 0] = DM_POLYTOPE_SEGMENT;
291: ex->cone[rt][coff + 1] = 0;
292: ex->cone[rt][coff + 2] = 0;
293: ex->cone[rt][coff + 3] = DM_POLYTOPE_SEGMENT;
294: ex->cone[rt][coff + 4] = 0;
295: ex->cone[rt][coff + 5] = 1;
296: ex->cone[rt][coff + 6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
297: ex->cone[rt][coff + 7] = 1;
298: ex->cone[rt][coff + 8] = 0;
299: ex->cone[rt][coff + 9] = 0;
300: ex->cone[rt][coff + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
301: ex->cone[rt][coff + 11] = 1;
302: ex->cone[rt][coff + 12] = 1;
303: ex->cone[rt][coff + 13] = 0;
304: ex->ornt[rt][ooff + 0] = 0;
305: ex->ornt[rt][ooff + 1] = 0;
306: ex->ornt[rt][ooff + 2] = 0;
307: ex->ornt[rt][ooff + 3] = 0;
308: } else {
309: ex->cone[rt][coff + 0] = DM_POLYTOPE_SEGMENT;
310: ex->cone[rt][coff + 1] = 0;
311: ex->cone[rt][coff + 2] = 0;
312: ex->cone[rt][coff + 3] = DM_POLYTOPE_SEGMENT;
313: ex->cone[rt][coff + 4] = 1;
314: ex->cone[rt][coff + 5] = 1;
315: ex->cone[rt][coff + 6] = 0;
316: ex->cone[rt][coff + 7] = DM_POLYTOPE_SEGMENT;
317: ex->cone[rt][coff + 8] = 0;
318: ex->cone[rt][coff + 9] = 1;
319: ex->cone[rt][coff + 10] = DM_POLYTOPE_SEGMENT;
320: ex->cone[rt][coff + 11] = 1;
321: ex->cone[rt][coff + 12] = 0;
322: ex->cone[rt][coff + 13] = 0;
323: ex->ornt[rt][ooff + 0] = 0;
324: ex->ornt[rt][ooff + 1] = 0;
325: ex->ornt[rt][ooff + 2] = -1;
326: ex->ornt[rt][ooff + 3] = -1;
327: }
328: }
329: // Impinging segment
330: // 0: no splits (impossible)
331: // 1: split vertex 0
332: // 2: split vertex 1
333: // 3: split both vertices (impossible)
334: for (PetscInt s = 1; s < 3; ++s) {
335: rt = (DM_POLYTOPE_SEGMENT * 2 + 0) * 100 + s;
336: ex->Nt[rt] = 1;
337: Nc = 8;
338: No = 2;
339: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
340: ex->target[rt][0] = DM_POLYTOPE_SEGMENT;
341: ex->size[rt][0] = 1;
342: // cone for segment
343: ex->cone[rt][0] = DM_POLYTOPE_POINT;
344: ex->cone[rt][1] = 1;
345: ex->cone[rt][2] = 0;
346: ex->cone[rt][3] = s == 1 ? 1 : 0;
347: ex->cone[rt][4] = DM_POLYTOPE_POINT;
348: ex->cone[rt][5] = 1;
349: ex->cone[rt][6] = 1;
350: ex->cone[rt][7] = s == 2 ? 1 : 0;
351: for (PetscInt i = 0; i < 2; ++i) ex->ornt[rt][i] = 0;
352: }
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /* DM_POLYTOPE_TRIANGLE produces
357: 2 triangles, and
358: 1 triangular prism/tensor triangular prism
359: */
360: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Triangle(DMPlexTransform_Cohesive *ex)
361: {
362: PetscInt rt, Nc, No, coff, ooff;
364: PetscFunctionBegin;
365: // No unsplit triangles
366: // Split triangles
367: // 0: no unsplit edge
368: // 1: unsplit edge 0
369: // 2: unsplit edge 1
370: // 3: unsplit edge 0 1
371: // 4: unsplit edge 2
372: // 5: unsplit edge 0 2
373: // 6: unsplit edge 1 2
374: // 7: all edges unsplit (impossible)
375: for (PetscInt s = 0; s < 7; ++s) {
376: rt = (DM_POLYTOPE_TRIANGLE * 2 + 1) * 100 + s;
377: ex->Nt[rt] = 2;
378: Nc = 12 * 2 + 18;
379: No = 3 * 2 + 5;
380: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
381: ex->target[rt][0] = DM_POLYTOPE_TRIANGLE;
382: ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM;
383: ex->size[rt][0] = 2;
384: ex->size[rt][1] = 1;
385: // cones for triangles
386: for (PetscInt i = 0; i < 2; ++i) {
387: ex->cone[rt][12 * i + 0] = DM_POLYTOPE_SEGMENT;
388: ex->cone[rt][12 * i + 1] = 1;
389: ex->cone[rt][12 * i + 2] = 0;
390: ex->cone[rt][12 * i + 3] = s & 1 ? 0 : i;
391: ex->cone[rt][12 * i + 4] = DM_POLYTOPE_SEGMENT;
392: ex->cone[rt][12 * i + 5] = 1;
393: ex->cone[rt][12 * i + 6] = 1;
394: ex->cone[rt][12 * i + 7] = s & 2 ? 0 : i;
395: ex->cone[rt][12 * i + 8] = DM_POLYTOPE_SEGMENT;
396: ex->cone[rt][12 * i + 9] = 1;
397: ex->cone[rt][12 * i + 10] = 2;
398: ex->cone[rt][12 * i + 11] = s & 4 ? 0 : i;
399: }
400: for (PetscInt i = 0; i < 3 * 2; ++i) ex->ornt[rt][i] = 0;
401: // cone for triangular prism/tensor triangular prism
402: coff = 12 * 2;
403: ooff = 3 * 2;
404: if (ex->useTensor) {
405: ex->cone[rt][coff + 0] = DM_POLYTOPE_TRIANGLE;
406: ex->cone[rt][coff + 1] = 0;
407: ex->cone[rt][coff + 2] = 0;
408: ex->cone[rt][coff + 3] = DM_POLYTOPE_TRIANGLE;
409: ex->cone[rt][coff + 4] = 0;
410: ex->cone[rt][coff + 5] = 1;
411: ex->cone[rt][coff + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
412: ex->cone[rt][coff + 7] = 1;
413: ex->cone[rt][coff + 8] = 0;
414: ex->cone[rt][coff + 9] = 0;
415: ex->cone[rt][coff + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
416: ex->cone[rt][coff + 11] = 1;
417: ex->cone[rt][coff + 12] = 1;
418: ex->cone[rt][coff + 13] = 0;
419: ex->cone[rt][coff + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
420: ex->cone[rt][coff + 15] = 1;
421: ex->cone[rt][coff + 16] = 2;
422: ex->cone[rt][coff + 17] = 0;
423: ex->ornt[rt][ooff + 0] = 0;
424: ex->ornt[rt][ooff + 1] = 0;
425: ex->ornt[rt][ooff + 2] = 0;
426: ex->ornt[rt][ooff + 3] = 0;
427: ex->ornt[rt][ooff + 4] = 0;
428: } else {
429: ex->cone[rt][coff + 0] = DM_POLYTOPE_TRIANGLE;
430: ex->cone[rt][coff + 1] = 0;
431: ex->cone[rt][coff + 2] = 0;
432: ex->cone[rt][coff + 3] = DM_POLYTOPE_TRIANGLE;
433: ex->cone[rt][coff + 4] = 0;
434: ex->cone[rt][coff + 5] = 1;
435: ex->cone[rt][coff + 6] = DM_POLYTOPE_QUADRILATERAL;
436: ex->cone[rt][coff + 7] = 1;
437: ex->cone[rt][coff + 8] = 0;
438: ex->cone[rt][coff + 9] = 0;
439: ex->cone[rt][coff + 10] = DM_POLYTOPE_QUADRILATERAL;
440: ex->cone[rt][coff + 11] = 1;
441: ex->cone[rt][coff + 12] = 1;
442: ex->cone[rt][coff + 13] = 0;
443: ex->cone[rt][coff + 14] = DM_POLYTOPE_QUADRILATERAL;
444: ex->cone[rt][coff + 15] = 1;
445: ex->cone[rt][coff + 16] = 2;
446: ex->cone[rt][coff + 17] = 0;
447: ex->ornt[rt][ooff + 0] = -2;
448: ex->ornt[rt][ooff + 1] = 0;
449: ex->ornt[rt][ooff + 2] = 0;
450: ex->ornt[rt][ooff + 3] = 0;
451: ex->ornt[rt][ooff + 4] = 0;
452: }
453: }
454: // Impinging triangles
455: // 0: no splits (impossible)
456: // 1: split edge 0
457: // 2: split edge 1
458: // 3: split edges 0 and 1
459: // 4: split edge 2
460: // 5: split edges 0 and 2
461: // 6: split edges 1 and 2
462: // 7: split all edges (impossible)
463: for (PetscInt s = 1; s < 7; ++s) {
464: rt = (DM_POLYTOPE_TRIANGLE * 2 + 0) * 100 + s;
465: ex->Nt[rt] = 1;
466: Nc = 12;
467: No = 3;
468: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
469: ex->target[rt][0] = DM_POLYTOPE_TRIANGLE;
470: ex->size[rt][0] = 1;
471: // cone for triangle
472: ex->cone[rt][0] = DM_POLYTOPE_SEGMENT;
473: ex->cone[rt][1] = 1;
474: ex->cone[rt][2] = 0;
475: ex->cone[rt][3] = s & 1 ? 1 : 0;
476: ex->cone[rt][4] = DM_POLYTOPE_SEGMENT;
477: ex->cone[rt][5] = 1;
478: ex->cone[rt][6] = 1;
479: ex->cone[rt][7] = s & 2 ? 1 : 0;
480: ex->cone[rt][8] = DM_POLYTOPE_SEGMENT;
481: ex->cone[rt][9] = 1;
482: ex->cone[rt][10] = 2;
483: ex->cone[rt][11] = s & 4 ? 1 : 0;
484: for (PetscInt i = 0; i < 3; ++i) ex->ornt[rt][i] = 0;
485: }
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: /* DM_POLYTOPE_QUADRILATERAL produces
490: 2 quads, and
491: 1 hex/tensor hex
492: */
493: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Quadrilateral(DMPlexTransform_Cohesive *ex)
494: {
495: PetscInt rt, Nc, No, coff, ooff;
497: PetscFunctionBegin;
498: // No unsplit quadrilaterals
499: // Split quadrilateral
500: // 0: no unsplit edge
501: // 1: unsplit edge 0
502: // 2: unsplit edge 1
503: // 3: unsplit edge 0 1
504: // 4: unsplit edge 2
505: // 5: unsplit edge 0 2
506: // 6: unsplit edge 1 2
507: // 7: unsplit edge 0 1 2
508: // 8: unsplit edge 3
509: // 9: unsplit edge 0 3
510: // 10: unsplit edge 1 3
511: // 11: unsplit edge 0 1 3
512: // 12: unsplit edge 2 3
513: // 13: unsplit edge 0 2 3
514: // 14: unsplit edge 1 2 3
515: // 15: all edges unsplit (impossible)
516: for (PetscInt s = 0; s < 15; ++s) {
517: rt = (DM_POLYTOPE_QUADRILATERAL * 2 + 1) * 100 + s;
518: ex->Nt[rt] = 2;
519: Nc = 16 * 2 + 22;
520: No = 4 * 2 + 6;
521: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
522: ex->target[rt][0] = DM_POLYTOPE_QUADRILATERAL;
523: ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON;
524: ex->size[rt][0] = 2;
525: ex->size[rt][1] = 1;
526: // cones for quads
527: for (PetscInt i = 0; i < 2; ++i) {
528: ex->cone[rt][16 * i + 0] = DM_POLYTOPE_SEGMENT;
529: ex->cone[rt][16 * i + 1] = 1;
530: ex->cone[rt][16 * i + 2] = 0;
531: ex->cone[rt][16 * i + 3] = s & 1 ? 0 : i;
532: ex->cone[rt][16 * i + 4] = DM_POLYTOPE_SEGMENT;
533: ex->cone[rt][16 * i + 5] = 1;
534: ex->cone[rt][16 * i + 6] = 1;
535: ex->cone[rt][16 * i + 7] = s & 2 ? 0 : i;
536: ex->cone[rt][16 * i + 8] = DM_POLYTOPE_SEGMENT;
537: ex->cone[rt][16 * i + 9] = 1;
538: ex->cone[rt][16 * i + 10] = 2;
539: ex->cone[rt][16 * i + 11] = s & 4 ? 0 : i;
540: ex->cone[rt][16 * i + 12] = DM_POLYTOPE_SEGMENT;
541: ex->cone[rt][16 * i + 13] = 1;
542: ex->cone[rt][16 * i + 14] = 3;
543: ex->cone[rt][16 * i + 15] = s & 8 ? 0 : i;
544: }
545: for (PetscInt i = 0; i < 4 * 2; ++i) ex->ornt[rt][i] = 0;
546: // cones for hexes/tensor hexes
547: coff = 16 * 2;
548: ooff = 4 * 2;
549: if (ex->useTensor) {
550: ex->cone[rt][coff + 0] = DM_POLYTOPE_QUADRILATERAL;
551: ex->cone[rt][coff + 1] = 0;
552: ex->cone[rt][coff + 2] = 0;
553: ex->cone[rt][coff + 3] = DM_POLYTOPE_QUADRILATERAL;
554: ex->cone[rt][coff + 4] = 0;
555: ex->cone[rt][coff + 5] = 1;
556: ex->cone[rt][coff + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
557: ex->cone[rt][coff + 7] = 1;
558: ex->cone[rt][coff + 8] = 0;
559: ex->cone[rt][coff + 9] = 0;
560: ex->cone[rt][coff + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
561: ex->cone[rt][coff + 11] = 1;
562: ex->cone[rt][coff + 12] = 1;
563: ex->cone[rt][coff + 13] = 0;
564: ex->cone[rt][coff + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
565: ex->cone[rt][coff + 15] = 1;
566: ex->cone[rt][coff + 16] = 2;
567: ex->cone[rt][coff + 17] = 0;
568: ex->cone[rt][coff + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
569: ex->cone[rt][coff + 19] = 1;
570: ex->cone[rt][coff + 20] = 3;
571: ex->cone[rt][coff + 21] = 0;
572: ex->ornt[rt][ooff + 0] = 0;
573: ex->ornt[rt][ooff + 1] = 0;
574: ex->ornt[rt][ooff + 2] = 0;
575: ex->ornt[rt][ooff + 3] = 0;
576: ex->ornt[rt][ooff + 4] = 0;
577: ex->ornt[rt][ooff + 5] = 0;
578: } else {
579: ex->cone[rt][coff + 0] = DM_POLYTOPE_QUADRILATERAL;
580: ex->cone[rt][coff + 1] = 0;
581: ex->cone[rt][coff + 2] = 0;
582: ex->cone[rt][coff + 3] = DM_POLYTOPE_QUADRILATERAL;
583: ex->cone[rt][coff + 4] = 0;
584: ex->cone[rt][coff + 5] = 1;
585: ex->cone[rt][coff + 6] = DM_POLYTOPE_QUADRILATERAL;
586: ex->cone[rt][coff + 7] = 1;
587: ex->cone[rt][coff + 8] = 0;
588: ex->cone[rt][coff + 9] = 0;
589: ex->cone[rt][coff + 10] = DM_POLYTOPE_QUADRILATERAL;
590: ex->cone[rt][coff + 11] = 1;
591: ex->cone[rt][coff + 12] = 2;
592: ex->cone[rt][coff + 13] = 0;
593: ex->cone[rt][coff + 14] = DM_POLYTOPE_QUADRILATERAL;
594: ex->cone[rt][coff + 15] = 1;
595: ex->cone[rt][coff + 16] = 1;
596: ex->cone[rt][coff + 17] = 0;
597: ex->cone[rt][coff + 18] = DM_POLYTOPE_QUADRILATERAL;
598: ex->cone[rt][coff + 19] = 1;
599: ex->cone[rt][coff + 20] = 3;
600: ex->cone[rt][coff + 21] = 0;
601: ex->ornt[rt][ooff + 0] = -2;
602: ex->ornt[rt][ooff + 1] = 0;
603: ex->ornt[rt][ooff + 2] = 0;
604: ex->ornt[rt][ooff + 3] = 0;
605: ex->ornt[rt][ooff + 4] = 0;
606: ex->ornt[rt][ooff + 5] = 1;
607: }
608: }
609: // Impinging quadrilaterals
610: // 0: no splits (impossible)
611: // 1: split edge 0
612: // 2: split edge 1
613: // 3: split edges 0 and 1
614: // 4: split edge 2
615: // 5: split edges 0 and 2
616: // 6: split edges 1 and 2
617: // 7: split edges 0, 1, and 2
618: // 8: split edge 3
619: // 9: split edges 0 and 3
620: // 10: split edges 1 and 3
621: // 11: split edges 0, 1, and 3
622: // 12: split edges 2 and 3
623: // 13: split edges 0, 2, and 3
624: // 14: split edges 1, 2, and 3
625: // 15: split all edges (impossible)
626: for (PetscInt s = 1; s < 15; ++s) {
627: rt = (DM_POLYTOPE_QUADRILATERAL * 2 + 0) * 100 + s;
628: ex->Nt[rt] = 1;
629: Nc = 16;
630: No = 4;
631: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
632: ex->target[rt][0] = DM_POLYTOPE_QUADRILATERAL;
633: ex->size[rt][0] = 1;
634: // cone for quadrilateral
635: ex->cone[rt][0] = DM_POLYTOPE_SEGMENT;
636: ex->cone[rt][1] = 1;
637: ex->cone[rt][2] = 0;
638: ex->cone[rt][3] = s & 1 ? 1 : 0;
639: ex->cone[rt][4] = DM_POLYTOPE_SEGMENT;
640: ex->cone[rt][5] = 1;
641: ex->cone[rt][6] = 1;
642: ex->cone[rt][7] = s & 2 ? 1 : 0;
643: ex->cone[rt][8] = DM_POLYTOPE_SEGMENT;
644: ex->cone[rt][9] = 1;
645: ex->cone[rt][10] = 2;
646: ex->cone[rt][11] = s & 4 ? 1 : 0;
647: ex->cone[rt][12] = DM_POLYTOPE_SEGMENT;
648: ex->cone[rt][13] = 1;
649: ex->cone[rt][14] = 3;
650: ex->cone[rt][15] = s & 8 ? 1 : 0;
651: for (PetscInt i = 0; i < 4; ++i) ex->ornt[rt][i] = 0;
652: }
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Tetrahedron(DMPlexTransform_Cohesive *ex)
657: {
658: PetscInt rt, Nc, No;
660: PetscFunctionBegin;
661: // Impinging tetrahedra
662: // 0: no splits (impossible)
663: // 1: split face 0
664: // 2: split face 1
665: // 3: split faces 0 and 1
666: // 4: split face 2
667: // 5: split faces 0 and 2
668: // 6: split faces 1 and 2
669: // 7: split faces 0, 1, and 2
670: // 8: split face 3
671: // 9: split faces 0 and 3
672: // 10: split faces 1 and 3
673: // 11: split faces 0, 1, and 3
674: // 12: split faces 2 and 3
675: // 13: split faces 0, 2, and 3
676: // 14: split faces 1, 2, and 3
677: // 15: split all faces (impossible)
678: for (PetscInt s = 1; s < 15; ++s) {
679: rt = (DM_POLYTOPE_TETRAHEDRON * 2 + 0) * 100 + s;
680: ex->Nt[rt] = 1;
681: Nc = 16;
682: No = 4;
683: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
684: ex->target[rt][0] = DM_POLYTOPE_TETRAHEDRON;
685: ex->size[rt][0] = 1;
686: // cone for triangle
687: ex->cone[rt][0] = DM_POLYTOPE_TRIANGLE;
688: ex->cone[rt][1] = 1;
689: ex->cone[rt][2] = 0;
690: ex->cone[rt][3] = s & 1 ? 1 : 0;
691: ex->cone[rt][4] = DM_POLYTOPE_TRIANGLE;
692: ex->cone[rt][5] = 1;
693: ex->cone[rt][6] = 1;
694: ex->cone[rt][7] = s & 2 ? 1 : 0;
695: ex->cone[rt][8] = DM_POLYTOPE_TRIANGLE;
696: ex->cone[rt][9] = 1;
697: ex->cone[rt][10] = 2;
698: ex->cone[rt][11] = s & 4 ? 1 : 0;
699: ex->cone[rt][12] = DM_POLYTOPE_TRIANGLE;
700: ex->cone[rt][13] = 1;
701: ex->cone[rt][14] = 3;
702: ex->cone[rt][15] = s & 8 ? 1 : 0;
703: for (PetscInt i = 0; i < 4; ++i) ex->ornt[rt][i] = 0;
704: }
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Hexahedron(DMPlexTransform_Cohesive *ex)
709: {
710: PetscInt rt, Nc, No;
712: PetscFunctionBegin;
713: // Impinging hexahedra
714: // 0: no splits (impossible)
715: // bit is set if the face is split
716: // 63: split all faces (impossible)
717: for (PetscInt s = 1; s < 63; ++s) {
718: rt = (DM_POLYTOPE_HEXAHEDRON * 2 + 0) * 100 + s;
719: ex->Nt[rt] = 1;
720: Nc = 24;
721: No = 6;
722: PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
723: ex->target[rt][0] = DM_POLYTOPE_HEXAHEDRON;
724: ex->size[rt][0] = 1;
725: // cone for hexahedron
726: ex->cone[rt][0] = DM_POLYTOPE_QUADRILATERAL;
727: ex->cone[rt][1] = 1;
728: ex->cone[rt][2] = 0;
729: ex->cone[rt][3] = s & 1 ? 1 : 0;
730: ex->cone[rt][4] = DM_POLYTOPE_QUADRILATERAL;
731: ex->cone[rt][5] = 1;
732: ex->cone[rt][6] = 1;
733: ex->cone[rt][7] = s & 2 ? 1 : 0;
734: ex->cone[rt][8] = DM_POLYTOPE_QUADRILATERAL;
735: ex->cone[rt][9] = 1;
736: ex->cone[rt][10] = 2;
737: ex->cone[rt][11] = s & 4 ? 1 : 0;
738: ex->cone[rt][12] = DM_POLYTOPE_QUADRILATERAL;
739: ex->cone[rt][13] = 1;
740: ex->cone[rt][14] = 3;
741: ex->cone[rt][15] = s & 8 ? 1 : 0;
742: ex->cone[rt][16] = DM_POLYTOPE_QUADRILATERAL;
743: ex->cone[rt][17] = 1;
744: ex->cone[rt][18] = 4;
745: ex->cone[rt][19] = s & 16 ? 1 : 0;
746: ex->cone[rt][20] = DM_POLYTOPE_QUADRILATERAL;
747: ex->cone[rt][21] = 1;
748: ex->cone[rt][22] = 5;
749: ex->cone[rt][23] = s & 32 ? 1 : 0;
750: for (PetscInt i = 0; i < 6; ++i) ex->ornt[rt][i] = 0;
751: }
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
755: /*
756: The refine types for cohesive extrusion are:
758: ct * 2 + 0: For any point which should just return itself
759: ct * 2 + 1: For unsplit points
760: (ct * 2 + 0) * 100 + fsplit: For impinging points, one type for each combination of split faces
761: (ct * 2 + 1) * 100 + fsplit: For split points, one type for each combination of unsplit faces
762: */
763: static PetscErrorCode DMPlexTransformSetUp_Cohesive(DMPlexTransform tr)
764: {
765: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
766: DM dm;
767: DMLabel active, celltype;
768: PetscInt numRt, pStart, pEnd, ict;
770: PetscFunctionBegin;
771: PetscCall(DMPlexTransformGetDM(tr, &dm));
772: PetscCall(DMPlexTransformGetActive(tr, &active));
773: PetscCheck(active, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cohesive extrusion requires an active label");
774: PetscCall(DMPlexGetCellTypeLabel(dm, &celltype));
775: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
776: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
777: PetscCall(DMLabelMakeAllInvalid_Internal(active));
778: for (PetscInt p = pStart; p < pEnd; ++p) {
779: PetscInt ct, val;
781: PetscCall(DMLabelGetValue(celltype, p, &ct));
782: PetscCall(DMLabelGetValue(active, p, &val));
783: if (val < 0) {
784: // Also negative size impinging points
785: // ct * 2 + 0 is the identity transform
786: PetscCall(DMLabelSetValue(tr->trType, p, ct * 2 + 0));
787: } else {
788: PetscInt fsplit = -1, funsplit = -1;
790: // Unsplit points ct * 2 + 1
791: if (val >= 200) {
792: // Cohesive cells cannot be unsplit
793: // This is faulty inheritance through the label
794: if (ct == DM_POLYTOPE_POINT_PRISM_TENSOR || ct == DM_POLYTOPE_SEG_PRISM_TENSOR || ct == DM_POLYTOPE_TRI_PRISM_TENSOR || ct == DM_POLYTOPE_QUAD_PRISM_TENSOR) PetscCall(DMLabelSetValue(tr->trType, p, ct * 2 + 0));
795: else PetscCall(DMLabelSetValue(tr->trType, p, ct * 2 + 1));
796: } else if (val >= 100) {
797: // Impinging points: (ct * 2 + 0) * 100 + fsplit
798: PetscCall(ComputeSplitFaceNumber(dm, active, p, &fsplit));
799: if (!fsplit) PetscCall(DMLabelSetValue(tr->trType, p, ct * 2 + 0));
800: else PetscCall(DMLabelSetValue(tr->trType, p, (ct * 2 + 0) * 100 + fsplit));
801: } else {
802: // Split points: (ct * 2 + 1) * 100 + funsplit
803: PetscCall(ComputeUnsplitFaceNumber(dm, active, p, &funsplit));
804: PetscCall(DMLabelSetValue(tr->trType, p, (ct * 2 + 1) * 100 + funsplit));
805: }
806: }
807: }
808: if (ex->debug) {
809: PetscCall(DMLabelView(active, NULL));
810: PetscCall(DMLabelView(tr->trType, NULL));
811: }
812: numRt = DM_NUM_POLYTOPES * 2 * 100;
813: PetscCall(PetscMalloc5(numRt, &ex->Nt, numRt, &ex->target, numRt, &ex->size, numRt, &ex->cone, numRt, &ex->ornt));
814: for (ict = 0; ict < numRt; ++ict) {
815: ex->Nt[ict] = -1;
816: ex->target[ict] = NULL;
817: ex->size[ict] = NULL;
818: ex->cone[ict] = NULL;
819: ex->ornt[ict] = NULL;
820: }
821: PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Point(ex));
822: PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Segment(ex));
823: PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Triangle(ex));
824: PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Quadrilateral(ex));
825: PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Tetrahedron(ex));
826: PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Hexahedron(ex));
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: static PetscErrorCode DMPlexTransformDestroy_Cohesive(DMPlexTransform tr)
831: {
832: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
833: PetscInt ct;
835: PetscFunctionBegin;
836: if (ex->target) {
837: for (ct = 0; ct < DM_NUM_POLYTOPES * 2 * 100; ++ct) PetscCall(PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct]));
838: }
839: PetscCall(PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt));
840: PetscCall(PetscFree(ex));
841: PetscFunctionReturn(PETSC_SUCCESS);
842: }
844: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Cohesive(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
845: {
846: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
847: DMLabel trType = tr->trType;
848: PetscInt rt;
850: PetscFunctionBeginHot;
851: *rnew = r;
852: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
853: if (!so) PetscFunctionReturn(PETSC_SUCCESS);
854: if (trType) {
855: PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
856: if (rt < 100 && !(rt % 2)) PetscFunctionReturn(PETSC_SUCCESS);
857: }
858: if (ex->useTensor) {
859: switch (sct) {
860: case DM_POLYTOPE_POINT:
861: break;
862: case DM_POLYTOPE_SEGMENT:
863: switch (tct) {
864: case DM_POLYTOPE_SEGMENT:
865: break;
866: case DM_POLYTOPE_SEG_PRISM_TENSOR:
867: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
868: break;
869: default:
870: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
871: }
872: break;
873: // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
874: case DM_POLYTOPE_TRIANGLE:
875: break;
876: case DM_POLYTOPE_QUADRILATERAL:
877: break;
878: default:
879: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
880: }
881: } else {
882: switch (sct) {
883: case DM_POLYTOPE_POINT:
884: break;
885: case DM_POLYTOPE_SEGMENT:
886: switch (tct) {
887: case DM_POLYTOPE_SEGMENT:
888: break;
889: case DM_POLYTOPE_QUADRILATERAL:
890: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0);
891: break;
892: default:
893: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
894: }
895: break;
896: // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
897: case DM_POLYTOPE_TRIANGLE:
898: break;
899: case DM_POLYTOPE_QUADRILATERAL:
900: break;
901: default:
902: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
903: }
904: }
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: static PetscErrorCode DMPlexTransformCellTransform_Cohesive(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
909: {
910: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
911: DMLabel trType = tr->trType;
912: PetscBool identity = PETSC_FALSE;
913: PetscInt val = 0;
915: PetscFunctionBegin;
916: PetscCheck(trType, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing transform type label");
917: PetscCall(DMLabelGetValue(trType, p, &val));
918: identity = val < 100 && !(val % 2) ? PETSC_TRUE : PETSC_FALSE;
919: if (rt) *rt = val;
920: if (identity) {
921: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
922: } else {
923: *Nt = ex->Nt[val];
924: *target = ex->target[val];
925: *size = ex->size[val];
926: *cone = ex->cone[val];
927: *ornt = ex->ornt[val];
928: }
929: PetscFunctionReturn(PETSC_SUCCESS);
930: }
932: /* New vertices have the same coordinates */
933: static PetscErrorCode DMPlexTransformMapCoordinates_Cohesive(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
934: {
935: PetscReal width;
936: PetscInt pval;
938: PetscFunctionBeginHot;
939: PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]);
940: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
941: PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv);
942: PetscCheck(r < 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should only have two replicas, not %" PetscInt_FMT, r);
944: PetscCall(DMPlexTransformCohesiveExtrudeGetWidth(tr, &width));
945: PetscCall(DMLabelGetValue(tr->trType, p, &pval));
946: if (width == 0. || pval < 100) {
947: for (PetscInt d = 0; d < dE; ++d) out[d] = in[d];
948: } else {
949: DM dm;
950: PetscReal avgNormal[3] = {0., 0., 0.}, norm = 0.;
951: PetscInt *star = NULL;
952: PetscInt Nst, fStart, fEnd, Nf = 0;
954: PetscCall(DMPlexTransformGetDM(tr, &dm));
955: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
956: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &Nst, &star));
957: // Get support faces that are split, refine type (ct * 2 + 1) * 100 + fsplit
958: for (PetscInt st = 0; st < Nst * 2; st += 2) {
959: DMPolytopeType ct;
960: PetscInt val;
962: if (star[st] < fStart || star[st] >= fEnd) continue;
963: PetscCall(DMPlexGetCellType(dm, star[st], &ct));
964: PetscCall(DMLabelGetValue(tr->trType, star[st], &val));
965: if (val < ((PetscInt)ct * 2 + 1) * 100) continue;
966: star[Nf++] = star[st];
967: }
968: PetscCheck(Nf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Split vertex %" PetscInt_FMT " must be connected to at least one split face", p);
969: // Average normals
970: for (PetscInt f = 0; f < Nf; ++f) {
971: PetscReal normal[3], vol;
973: PetscCall(DMPlexComputeCellGeometryFVM(dm, star[f], &vol, NULL, normal));
974: for (PetscInt d = 0; d < dE; ++d) avgNormal[d] += normal[d];
975: }
976: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &Nst, &star));
977: // Normalize normal
978: for (PetscInt d = 0; d < dE; ++d) norm += PetscSqr(avgNormal[d]);
979: norm = PetscSqrtReal(norm);
980: for (PetscInt d = 0; d < dE; ++d) avgNormal[d] /= norm;
981: // Symmetrically push vertices along normal
982: for (PetscInt d = 0; d < dE; ++d) out[d] = in[d] + width * avgNormal[d] * (r ? -0.5 : 0.5);
983: }
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: static PetscErrorCode DMPlexTransformInitialize_Cohesive(DMPlexTransform tr)
988: {
989: PetscFunctionBegin;
990: tr->ops->view = DMPlexTransformView_Cohesive;
991: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Cohesive;
992: tr->ops->setup = DMPlexTransformSetUp_Cohesive;
993: tr->ops->destroy = DMPlexTransformDestroy_Cohesive;
994: tr->ops->setdimensions = DMPlexTransformSetDimensions_Internal;
995: tr->ops->celltransform = DMPlexTransformCellTransform_Cohesive;
996: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Cohesive;
997: tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Cohesive;
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform tr)
1002: {
1003: DMPlexTransform_Cohesive *ex;
1005: PetscFunctionBegin;
1007: PetscCall(PetscNew(&ex));
1008: tr->data = ex;
1009: ex->useTensor = PETSC_TRUE;
1010: PetscCall(DMPlexTransformInitialize_Cohesive(tr));
1011: PetscFunctionReturn(PETSC_SUCCESS);
1012: }
1014: /*@
1015: DMPlexTransformCohesiveExtrudeGetTensor - Get the flag to use tensor cells
1017: Not Collective
1019: Input Parameter:
1020: . tr - The `DMPlexTransform`
1022: Output Parameter:
1023: . useTensor - The flag to use tensor cells
1025: Note:
1026: This flag determines the orientation behavior of the created points.
1028: For example, if tensor is `PETSC_TRUE`, then
1029: .vb
1030: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1031: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1032: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1033: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1034: .ve
1036: Level: intermediate
1038: .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeSetTensor()`, `DMPlexTransformExtrudeGetTensor()`
1039: @*/
1040: PetscErrorCode DMPlexTransformCohesiveExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
1041: {
1042: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
1044: PetscFunctionBegin;
1046: PetscAssertPointer(useTensor, 2);
1047: *useTensor = ex->useTensor;
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: /*@
1052: DMPlexTransformCohesiveExtrudeSetTensor - Set the flag to use tensor cells
1054: Not Collective
1056: Input Parameters:
1057: + tr - The `DMPlexTransform`
1058: - useTensor - The flag for tensor cells
1060: Note:
1061: This flag determines the orientation behavior of the created points
1062: For example, if tensor is `PETSC_TRUE`, then
1063: .vb
1064: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1065: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1066: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1067: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1068: .ve
1070: Level: intermediate
1072: .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeGetTensor()`, `DMPlexTransformExtrudeSetTensor()`
1073: @*/
1074: PetscErrorCode DMPlexTransformCohesiveExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
1075: {
1076: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
1078: PetscFunctionBegin;
1080: ex->useTensor = useTensor;
1081: PetscFunctionReturn(PETSC_SUCCESS);
1082: }
1084: /*@
1085: DMPlexTransformCohesiveExtrudeGetWidth - Get the width of extruded cells
1087: Not Collective
1089: Input Parameter:
1090: . tr - The `DMPlexTransform`
1092: Output Parameter:
1093: . width - The width of extruded cells, or 0.
1095: Level: intermediate
1097: .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeSetWidth()`
1098: @*/
1099: PetscErrorCode DMPlexTransformCohesiveExtrudeGetWidth(DMPlexTransform tr, PetscReal *width)
1100: {
1101: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
1103: PetscFunctionBegin;
1105: PetscAssertPointer(width, 2);
1106: *width = ex->width;
1107: PetscFunctionReturn(PETSC_SUCCESS);
1108: }
1110: /*@
1111: DMPlexTransformCohesiveExtrudeSetWidth - Set the width of extruded cells
1113: Not Collective
1115: Input Parameters:
1116: + tr - The `DMPlexTransform`
1117: - width - The width of the extruded cells, or 0.
1119: Level: intermediate
1121: .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeGetWidth()`
1122: @*/
1123: PetscErrorCode DMPlexTransformCohesiveExtrudeSetWidth(DMPlexTransform tr, PetscReal width)
1124: {
1125: DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
1127: PetscFunctionBegin;
1129: ex->width = width;
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }