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