1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/
2
3 #include <petsc/private/dmlabelimpl.h> // For DMLabelMakeAllInvalid_Internal()
4
5 /*
6 The cohesive transformation extrudes cells into a mesh from faces along an internal boundary.
7
8 Orientation:
9
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.
11
12 Labeling:
13
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.
18
19 Point Creation:
20
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.
27
28 */
29
DMPlexTransformView_Cohesive(DMPlexTransform tr,PetscViewer viewer)30 static PetscErrorCode DMPlexTransformView_Cohesive(DMPlexTransform tr, PetscViewer viewer)
31 {
32 DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
33 PetscBool isascii;
34
35 PetscFunctionBegin;
36 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
37 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
38 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
39 if (isascii) {
40 const char *name;
41
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 }
50
DMPlexTransformSetFromOptions_Cohesive(DMPlexTransform tr,PetscOptionItems PetscOptionsObject)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;
56
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 }
67
68 /*
69 ComputeSplitFaceNumber - Compute an encoding describing which faces of p are split by the surface
70
71 Not collective
72
73 Input Parameters:
74 + dm - The `DM`
75 . label - `DMLabel` marking the surface and adjacent points
76 - p - Impinging point, adjacent to the surface
77
78 Output Parameter:
79 . fsplit - A number encoding the faces which are split by the surface
80
81 Level: developer
82
83 Note: We will use a bit encoding, where bit k is 1 if face k is split.
84
85 .seealso: ComputeUnsplitFaceNumber()
86 */
ComputeSplitFaceNumber(DM dm,DMLabel label,PetscInt p,PetscInt * fsplit)87 static PetscErrorCode ComputeSplitFaceNumber(DM dm, DMLabel label, PetscInt p, PetscInt *fsplit)
88 {
89 const PetscInt *cone;
90 PetscInt coneSize, val;
91
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 }
103
104 /*
105 ComputeUnsplitFaceNumber - Compute an encoding describing which faces of p are unsplit by the surface
106
107 Not collective
108
109 Input Parameters:
110 + dm - The `DM`
111 . label - `DMLabel` marking the surface and adjacent points
112 - p - Split point, on the surface
113
114 Output Parameter:
115 . funsplit - A number encoding the faces which are split by the surface
116
117 Level: developer
118
119 Note: We will use a bit encoding, where bit k is 1 if face k is unsplit.
120
121 .seealso: ComputeSplitFaceNumber()
122 */
ComputeUnsplitFaceNumber(DM dm,DMLabel label,PetscInt p,PetscInt * funsplit)123 static PetscErrorCode ComputeUnsplitFaceNumber(DM dm, DMLabel label, PetscInt p, PetscInt *funsplit)
124 {
125 const PetscInt *cone;
126 PetscInt coneSize, val;
127
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 }
139
140 /* DM_POLYTOPE_POINT produces
141 2 points when split, or 1 point when unsplit, and
142 1 segment, or tensor segment
143 */
DMPlexTransformCohesiveExtrudeSetUp_Point(DMPlexTransform_Cohesive * ex)144 static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Point(DMPlexTransform_Cohesive *ex)
145 {
146 PetscInt rt, Nc, No;
147
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 }
187
188 /* DM_POLYTOPE_SEGMENT produces
189 2 segments when split, or 1 segment when unsplit, and
190 1 quad, or tensor quad
191 */
DMPlexTransformCohesiveExtrudeSetUp_Segment(DMPlexTransform_Cohesive * ex)192 static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Segment(DMPlexTransform_Cohesive *ex)
193 {
194 PetscInt rt, Nc, No, coff, ooff;
195
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 }
355
356 /* DM_POLYTOPE_TRIANGLE produces
357 2 triangles, and
358 1 triangular prism/tensor triangular prism
359 */
DMPlexTransformCohesiveExtrudeSetUp_Triangle(DMPlexTransform_Cohesive * ex)360 static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Triangle(DMPlexTransform_Cohesive *ex)
361 {
362 PetscInt rt, Nc, No, coff, ooff;
363
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 }
488
489 /* DM_POLYTOPE_QUADRILATERAL produces
490 2 quads, and
491 1 hex/tensor hex
492 */
DMPlexTransformCohesiveExtrudeSetUp_Quadrilateral(DMPlexTransform_Cohesive * ex)493 static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Quadrilateral(DMPlexTransform_Cohesive *ex)
494 {
495 PetscInt rt, Nc, No, coff, ooff;
496
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 }
655
DMPlexTransformCohesiveExtrudeSetUp_Tetrahedron(DMPlexTransform_Cohesive * ex)656 static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Tetrahedron(DMPlexTransform_Cohesive *ex)
657 {
658 PetscInt rt, Nc, No;
659
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 }
707
DMPlexTransformCohesiveExtrudeSetUp_Hexahedron(DMPlexTransform_Cohesive * ex)708 static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Hexahedron(DMPlexTransform_Cohesive *ex)
709 {
710 PetscInt rt, Nc, No;
711
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 }
754
755 /*
756 The refine types for cohesive extrusion are:
757
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 */
DMPlexTransformSetUp_Cohesive(DMPlexTransform tr)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;
769
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;
780
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;
789
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 }
829
DMPlexTransformDestroy_Cohesive(DMPlexTransform tr)830 static PetscErrorCode DMPlexTransformDestroy_Cohesive(DMPlexTransform tr)
831 {
832 DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
833 PetscInt ct;
834
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 }
843
DMPlexTransformGetSubcellOrientation_Cohesive(DMPlexTransform tr,DMPolytopeType sct,PetscInt sp,PetscInt so,DMPolytopeType tct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)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;
849
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 }
907
DMPlexTransformCellTransform_Cohesive(DMPlexTransform tr,DMPolytopeType source,PetscInt p,PetscInt * rt,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])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;
914
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 }
931
OrderCohesiveSupport_Private(DM dm,PetscInt p)932 static PetscErrorCode OrderCohesiveSupport_Private(DM dm, PetscInt p)
933 {
934 const PetscInt *cone;
935 PetscInt csize;
936
937 PetscFunctionBegin;
938 PetscCall(DMPlexGetCone(dm, p, &cone));
939 PetscCall(DMPlexGetConeSize(dm, p, &csize));
940 PetscCheck(csize > 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cone size for cohesive cell should be > 2, not %" PetscInt_FMT, csize);
941 for (PetscInt s = 0; s < 2; ++s) {
942 const PetscInt *supp;
943 PetscInt Ns, neighbor;
944
945 PetscCall(DMPlexGetSupport(dm, cone[s], &supp));
946 PetscCall(DMPlexGetSupportSize(dm, cone[s], &Ns));
947 // Could check here that the face is in the pointSF
948 if (Ns == 1) continue;
949 PetscCheck(Ns == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cohesive cell endcap should have support of size 2, not %" PetscInt_FMT, Ns);
950 PetscCheck((supp[0] == p) != (supp[1] == p), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cohesive cell %" PetscInt_FMT " endcap must have cell in support once", p);
951 neighbor = supp[s] == p ? supp[1 - s] : -1;
952 if (neighbor >= 0) {
953 PetscCall(DMPlexInsertSupport(dm, cone[s], s, neighbor));
954 PetscCall(DMPlexInsertSupport(dm, cone[s], 1 - s, p));
955 }
956 }
957 PetscFunctionReturn(PETSC_SUCCESS);
958 }
959
960 // We need the supports of endcap faces on cohesive cells to have the same orientation
961 // We cannot just fix split points, since we destroy support ordering with DMPLexSymmetrize()
DMPlexTransformOrderSupports_Cohesive(DMPlexTransform tr,DM dm,DM tdm)962 static PetscErrorCode DMPlexTransformOrderSupports_Cohesive(DMPlexTransform tr, DM dm, DM tdm)
963 {
964 PetscInt dim, pStart, pEnd;
965
966 PetscFunctionBegin;
967 PetscCall(DMGetDimension(dm, &dim));
968 PetscCall(DMPlexGetChart(tdm, &pStart, &pEnd));
969 for (PetscInt p = pStart; p < pEnd; ++p) {
970 DMPolytopeType ct;
971
972 PetscCall(DMPlexGetCellType(tdm, p, &ct));
973 switch (dim) {
974 case 2:
975 if (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) PetscCall(OrderCohesiveSupport_Private(tdm, p));
976 break;
977 case 3:
978 if (ct == DM_POLYTOPE_TRI_PRISM_TENSOR || ct == DM_POLYTOPE_QUAD_PRISM_TENSOR) PetscCall(OrderCohesiveSupport_Private(tdm, p));
979 break;
980 default:
981 break;
982 }
983 }
984 PetscFunctionReturn(PETSC_SUCCESS);
985 }
986
987 /* New vertices have the same coordinates */
DMPlexTransformMapCoordinates_Cohesive(DMPlexTransform tr,DMPolytopeType pct,DMPolytopeType ct,PetscInt p,PetscInt r,PetscInt Nv,PetscInt dE,const PetscScalar in[],PetscScalar out[])988 static PetscErrorCode DMPlexTransformMapCoordinates_Cohesive(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
989 {
990 PetscReal width;
991 PetscInt pval;
992
993 PetscFunctionBeginHot;
994 PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]);
995 PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
996 PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv);
997 PetscCheck(r < 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should only have two replicas, not %" PetscInt_FMT, r);
998
999 PetscCall(DMPlexTransformCohesiveExtrudeGetWidth(tr, &width));
1000 PetscCall(DMLabelGetValue(tr->trType, p, &pval));
1001 if (width == 0. || pval < 100) {
1002 for (PetscInt d = 0; d < dE; ++d) out[d] = in[d];
1003 } else {
1004 DM dm;
1005 PetscReal avgNormal[3] = {0., 0., 0.}, norm = 0.;
1006 PetscInt *star = NULL;
1007 PetscInt Nst, fStart, fEnd, Nf = 0;
1008
1009 PetscCall(DMPlexTransformGetDM(tr, &dm));
1010 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
1011 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &Nst, &star));
1012 // Get support faces that are split, refine type (ct * 2 + 1) * 100 + fsplit
1013 for (PetscInt st = 0; st < Nst * 2; st += 2) {
1014 DMPolytopeType ct;
1015 PetscInt val;
1016
1017 if (star[st] < fStart || star[st] >= fEnd) continue;
1018 PetscCall(DMPlexGetCellType(dm, star[st], &ct));
1019 PetscCall(DMLabelGetValue(tr->trType, star[st], &val));
1020 if (val < ((PetscInt)ct * 2 + 1) * 100) continue;
1021 star[Nf++] = star[st];
1022 }
1023 PetscCheck(Nf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Split vertex %" PetscInt_FMT " must be connected to at least one split face", p);
1024 // Average normals
1025 for (PetscInt f = 0; f < Nf; ++f) {
1026 PetscReal normal[3], vol;
1027
1028 PetscCall(DMPlexComputeCellGeometryFVM(dm, star[f], &vol, NULL, normal));
1029 for (PetscInt d = 0; d < dE; ++d) avgNormal[d] += normal[d];
1030 }
1031 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &Nst, &star));
1032 // Normalize normal
1033 for (PetscInt d = 0; d < dE; ++d) norm += PetscSqr(avgNormal[d]);
1034 norm = PetscSqrtReal(norm);
1035 for (PetscInt d = 0; d < dE; ++d) avgNormal[d] /= norm;
1036 // Symmetrically push vertices along normal
1037 for (PetscInt d = 0; d < dE; ++d) out[d] = in[d] + width * avgNormal[d] * (r ? -0.5 : 0.5);
1038 }
1039 PetscFunctionReturn(PETSC_SUCCESS);
1040 }
1041
DMPlexTransformInitialize_Cohesive(DMPlexTransform tr)1042 static PetscErrorCode DMPlexTransformInitialize_Cohesive(DMPlexTransform tr)
1043 {
1044 PetscFunctionBegin;
1045 tr->ops->view = DMPlexTransformView_Cohesive;
1046 tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Cohesive;
1047 tr->ops->setup = DMPlexTransformSetUp_Cohesive;
1048 tr->ops->destroy = DMPlexTransformDestroy_Cohesive;
1049 tr->ops->setdimensions = DMPlexTransformSetDimensions_Internal;
1050 tr->ops->celltransform = DMPlexTransformCellTransform_Cohesive;
1051 tr->ops->ordersupports = DMPlexTransformOrderSupports_Cohesive;
1052 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Cohesive;
1053 tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Cohesive;
1054 PetscFunctionReturn(PETSC_SUCCESS);
1055 }
1056
DMPlexTransformCreate_Cohesive(DMPlexTransform tr)1057 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform tr)
1058 {
1059 DMPlexTransform_Cohesive *ex;
1060
1061 PetscFunctionBegin;
1062 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1063 PetscCall(PetscNew(&ex));
1064 tr->data = ex;
1065 ex->useTensor = PETSC_TRUE;
1066 PetscCall(DMPlexTransformInitialize_Cohesive(tr));
1067 PetscFunctionReturn(PETSC_SUCCESS);
1068 }
1069
1070 /*@
1071 DMPlexTransformCohesiveExtrudeGetTensor - Get the flag to use tensor cells
1072
1073 Not Collective
1074
1075 Input Parameter:
1076 . tr - The `DMPlexTransform`
1077
1078 Output Parameter:
1079 . useTensor - The flag to use tensor cells
1080
1081 Note:
1082 This flag determines the orientation behavior of the created points.
1083
1084 For example, if tensor is `PETSC_TRUE`, then
1085 .vb
1086 DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1087 DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1088 DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1089 DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1090 .ve
1091
1092 Level: intermediate
1093
1094 .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeSetTensor()`, `DMPlexTransformExtrudeGetTensor()`
1095 @*/
DMPlexTransformCohesiveExtrudeGetTensor(DMPlexTransform tr,PetscBool * useTensor)1096 PetscErrorCode DMPlexTransformCohesiveExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
1097 {
1098 DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
1099
1100 PetscFunctionBegin;
1101 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1102 PetscAssertPointer(useTensor, 2);
1103 *useTensor = ex->useTensor;
1104 PetscFunctionReturn(PETSC_SUCCESS);
1105 }
1106
1107 /*@
1108 DMPlexTransformCohesiveExtrudeSetTensor - Set the flag to use tensor cells
1109
1110 Not Collective
1111
1112 Input Parameters:
1113 + tr - The `DMPlexTransform`
1114 - useTensor - The flag for tensor cells
1115
1116 Note:
1117 This flag determines the orientation behavior of the created points
1118 For example, if tensor is `PETSC_TRUE`, then
1119 .vb
1120 DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1121 DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1122 DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1123 DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1124 .ve
1125
1126 Level: intermediate
1127
1128 .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeGetTensor()`, `DMPlexTransformExtrudeSetTensor()`
1129 @*/
DMPlexTransformCohesiveExtrudeSetTensor(DMPlexTransform tr,PetscBool useTensor)1130 PetscErrorCode DMPlexTransformCohesiveExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
1131 {
1132 DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
1133
1134 PetscFunctionBegin;
1135 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1136 ex->useTensor = useTensor;
1137 PetscFunctionReturn(PETSC_SUCCESS);
1138 }
1139
1140 /*@
1141 DMPlexTransformCohesiveExtrudeGetWidth - Get the width of extruded cells
1142
1143 Not Collective
1144
1145 Input Parameter:
1146 . tr - The `DMPlexTransform`
1147
1148 Output Parameter:
1149 . width - The width of extruded cells, or 0.
1150
1151 Level: intermediate
1152
1153 .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeSetWidth()`
1154 @*/
DMPlexTransformCohesiveExtrudeGetWidth(DMPlexTransform tr,PetscReal * width)1155 PetscErrorCode DMPlexTransformCohesiveExtrudeGetWidth(DMPlexTransform tr, PetscReal *width)
1156 {
1157 DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
1158
1159 PetscFunctionBegin;
1160 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1161 PetscAssertPointer(width, 2);
1162 *width = ex->width;
1163 PetscFunctionReturn(PETSC_SUCCESS);
1164 }
1165
1166 /*@
1167 DMPlexTransformCohesiveExtrudeSetWidth - Set the width of extruded cells
1168
1169 Not Collective
1170
1171 Input Parameters:
1172 + tr - The `DMPlexTransform`
1173 - width - The width of the extruded cells, or 0.
1174
1175 Level: intermediate
1176
1177 .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeGetWidth()`
1178 @*/
DMPlexTransformCohesiveExtrudeSetWidth(DMPlexTransform tr,PetscReal width)1179 PetscErrorCode DMPlexTransformCohesiveExtrudeSetWidth(DMPlexTransform tr, PetscReal width)
1180 {
1181 DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
1182
1183 PetscFunctionBegin;
1184 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1185 ex->width = width;
1186 PetscFunctionReturn(PETSC_SUCCESS);
1187 }
1188
1189 /*@
1190 DMPlexTransformCohesiveExtrudeGetUnsplit - Get a new label marking the unsplit points in the transformed mesh
1191
1192 Not Collective
1193
1194 Input Parameter:
1195 . tr - The `DMPlexTransform`
1196
1197 Output Parameter:
1198 . unsplit - A new `DMLabel` marking the unsplit points in the transformed mesh
1199
1200 Level: intermediate
1201
1202 Note:
1203 This label should be destroyed by the caller.
1204
1205 .seealso: `DMPlexTransform`, `DMPlexTransformGetTransformTypes()`
1206 @*/
DMPlexTransformCohesiveExtrudeGetUnsplit(DMPlexTransform tr,DMLabel * unsplit)1207 PetscErrorCode DMPlexTransformCohesiveExtrudeGetUnsplit(DMPlexTransform tr, DMLabel *unsplit)
1208 {
1209 DM dm;
1210 DMLabel trTypes;
1211 IS valueIS;
1212 const PetscInt *values;
1213 PetscInt Nv;
1214
1215 PetscFunctionBegin;
1216 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1217 PetscAssertPointer(unsplit, 2);
1218 PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)tr), "Unsplit Points", unsplit));
1219 PetscCall(DMPlexTransformGetDM(tr, &dm));
1220 PetscCall(DMPlexTransformGetTransformTypes(tr, &trTypes));
1221 PetscCall(DMLabelGetValueIS(trTypes, &valueIS));
1222 PetscCall(ISGetLocalSize(valueIS, &Nv));
1223 PetscCall(ISGetIndices(valueIS, &values));
1224 for (PetscInt v = 0; v < Nv; ++v) {
1225 const PetscInt val = values[v];
1226 IS pointIS;
1227 const PetscInt *points;
1228 PetscInt Np;
1229
1230 if (val > 2 * DM_NUM_POLYTOPES || !(val % 2)) continue;
1231 PetscCall(DMLabelGetStratumIS(trTypes, val, &pointIS));
1232 PetscCall(ISGetLocalSize(pointIS, &Np));
1233 PetscCall(ISGetIndices(pointIS, &points));
1234 for (PetscInt p = 0; p < Np; ++p) {
1235 const PetscInt point = points[p];
1236 DMPolytopeType ct;
1237 DMPolytopeType *rct;
1238 PetscInt *rsize, *rcone, *rornt;
1239 PetscInt Nct, pNew = 0;
1240
1241 PetscCall(DMPlexGetCellType(dm, point, &ct));
1242 PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1243 for (PetscInt n = 0; n < Nct; ++n) {
1244 for (PetscInt r = 0; r < rsize[n]; ++r) {
1245 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
1246 PetscCall(DMLabelSetValue(*unsplit, pNew, val + tr->labelReplicaInc * r));
1247 }
1248 }
1249 }
1250 PetscCall(ISRestoreIndices(pointIS, &points));
1251 PetscCall(ISDestroy(&pointIS));
1252 }
1253 PetscCall(ISRestoreIndices(valueIS, &values));
1254 PetscCall(ISDestroy(&valueIS));
1255 PetscFunctionReturn(PETSC_SUCCESS);
1256 }
1257