Actual source code: plextrextrude.c
1: #include <petsc/private/dmplextransformimpl.h>
3: const char *const PlexNormalAlgs[] = {"default", "input", "compute", "compute_bd"};
5: static PetscErrorCode DMPlexTransformView_Extrude(DMPlexTransform tr, PetscViewer viewer)
6: {
7: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
8: PetscBool isascii;
10: PetscFunctionBegin;
13: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
14: if (isascii) {
15: DM dm;
16: DMLabel active;
17: PetscInt dim;
18: const char *name;
20: PetscCall(PetscObjectGetName((PetscObject)tr, &name));
21: PetscCall(DMPlexTransformGetDM(tr, &dm));
22: PetscCall(DMGetDimension(dm, &dim));
23: PetscCall(DMPlexTransformGetActive(tr, &active));
25: PetscCall(PetscViewerASCIIPrintf(viewer, "Extrusion transformation %s\n", name ? name : ""));
26: PetscCall(PetscViewerASCIIPrintf(viewer, " number of layers: %" PetscInt_FMT "\n", ex->layers));
27: PetscCall(PetscViewerASCIIPrintf(viewer, " create tensor cells: %s\n", ex->useTensor ? "YES" : "NO"));
28: if (ex->periodic) PetscCall(PetscViewerASCIIPrintf(viewer, " periodic\n"));
29: PetscCall(PetscViewerASCIIPrintf(viewer, " normal algorithm: %s\n", PlexNormalAlgs[ex->normalAlg]));
30: if (ex->normalFunc) PetscCall(PetscViewerASCIIPrintf(viewer, " normal modified by user function\n"));
31: } else {
32: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
33: }
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode DMPlexTransformSetFromOptions_Extrude(DMPlexTransform tr, PetscOptionItems PetscOptionsObject)
38: {
39: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
40: PetscReal th, normal[3], *thicknesses;
41: PetscInt nl, Nc;
42: PetscBool tensor, sym, per, flg;
43: char funcname[PETSC_MAX_PATH_LEN];
45: PetscFunctionBegin;
46: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Extrusion Options");
47: PetscCall(PetscOptionsBoundedInt("-dm_plex_transform_extrude_layers", "Number of layers to extrude", "", ex->layers, &nl, &flg, 1));
48: if (flg) PetscCall(DMPlexTransformExtrudeSetLayers(tr, nl));
49: PetscCall(PetscOptionsReal("-dm_plex_transform_extrude_thickness", "Total thickness of extruded layers", "", ex->thickness, &th, &flg));
50: if (flg) PetscCall(DMPlexTransformExtrudeSetThickness(tr, th));
51: PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg));
52: if (flg) PetscCall(DMPlexTransformExtrudeSetTensor(tr, tensor));
53: PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_symmetric", "Extrude layers symmetrically about the surface", "", ex->symmetric, &sym, &flg));
54: if (flg) PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, sym));
55: PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_periodic", "Extrude layers periodically about the surface", "", ex->periodic, &per, &flg));
56: if (flg) PetscCall(DMPlexTransformExtrudeSetPeriodic(tr, per));
57: Nc = 3;
58: PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_normal", "Input normal vector for extrusion", "DMPlexTransformExtrudeSetNormal", normal, &Nc, &flg));
59: if (flg) {
60: // Extrusion dimension might not yet be determined
61: PetscCheck(!ex->cdimEx || Nc == ex->cdimEx, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_SIZ, "Input normal has size %" PetscInt_FMT " != %" PetscInt_FMT " extruded coordinate dimension", Nc, ex->cdimEx);
62: PetscCall(DMPlexTransformExtrudeSetNormal(tr, normal));
63: }
64: PetscCall(PetscOptionsString("-dm_plex_transform_extrude_normal_function", "Function to determine normal vector", "DMPlexTransformExtrudeSetNormalFunction", funcname, funcname, sizeof(funcname), &flg));
65: if (flg) {
66: PetscSimplePointFn *normalFunc;
68: PetscCall(PetscDLSym(NULL, funcname, (void **)&normalFunc));
69: PetscCall(DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc));
70: }
71: nl = ex->layers;
72: PetscCall(PetscMalloc1(nl, &thicknesses));
73: PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_thicknesses", "Thickness of each individual extruded layer", "", thicknesses, &nl, &flg));
74: if (flg) {
75: PetscCheck(nl, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one thickness for -dm_plex_transform_extrude_thicknesses");
76: PetscCall(DMPlexTransformExtrudeSetThicknesses(tr, nl, thicknesses));
77: }
78: PetscCall(PetscFree(thicknesses));
79: PetscOptionsHeadEnd();
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: /* Determine the implicit dimension pre-extrusion (either the implicit dimension of the DM or of a point in the active set for the transform).
84: If that dimension is the same as the current coordinate dimension (ex->dim), the extruded mesh will have a coordinate dimension one greater;
85: Otherwise the coordinate dimension will be kept. */
86: static PetscErrorCode DMPlexTransformExtrudeComputeExtrusionDim(DMPlexTransform tr)
87: {
88: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
89: DM dm;
90: DMLabel active;
91: PetscInt dim, dimExtPoint, dimExtPointG;
93: PetscFunctionBegin;
94: PetscCall(DMPlexTransformGetDM(tr, &dm));
95: PetscCall(DMGetDimension(dm, &dim));
96: PetscCall(DMPlexTransformGetActive(tr, &active));
97: if (active) {
98: IS valueIS, pointIS;
99: const PetscInt *values, *points;
100: DMPolytopeType ct;
101: PetscInt Nv, Np;
103: dimExtPoint = 0;
104: PetscCall(DMLabelGetValueIS(active, &valueIS));
105: PetscCall(ISGetLocalSize(valueIS, &Nv));
106: PetscCall(ISGetIndices(valueIS, &values));
107: for (PetscInt v = 0; v < Nv; ++v) {
108: PetscCall(DMLabelGetStratumIS(active, values[v], &pointIS));
109: PetscCall(ISGetLocalSize(pointIS, &Np));
110: PetscCall(ISGetIndices(pointIS, &points));
111: for (PetscInt p = 0; p < Np; ++p) {
112: PetscCall(DMPlexGetCellType(dm, points[p], &ct));
113: dimExtPoint = PetscMax(dimExtPoint, DMPolytopeTypeGetDim(ct));
114: }
115: PetscCall(ISRestoreIndices(pointIS, &points));
116: PetscCall(ISDestroy(&pointIS));
117: }
118: PetscCall(ISRestoreIndices(valueIS, &values));
119: PetscCall(ISDestroy(&valueIS));
120: } else dimExtPoint = dim;
121: PetscCallMPI(MPIU_Allreduce(&dimExtPoint, &dimExtPointG, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)tr)));
122: ex->dimEx = PetscMax(dim, dimExtPointG + 1);
123: ex->cdimEx = ex->cdim == dimExtPointG ? ex->cdim + 1 : ex->cdim;
124: PetscCheck(ex->dimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Topological dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->dimEx);
125: PetscCheck(ex->cdimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->cdimEx);
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm)
130: {
131: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
132: PetscInt dim;
134: PetscFunctionBegin;
135: PetscCall(DMGetDimension(dm, &dim));
136: PetscCall(DMSetDimension(tdm, ex->dimEx));
137: PetscCall(DMSetCoordinateDim(tdm, ex->cdimEx));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /* DM_POLYTOPE_POINT produces
142: Nl+1 points, or Nl points periodically and
143: Nl segments, or tensor segments
144: */
145: static PetscErrorCode DMPlexTransformExtrudeSetUp_Point(DMPlexTransform_Extrude *ex, PetscInt Nl)
146: {
147: const DMPolytopeType ct = DM_POLYTOPE_POINT;
148: const PetscInt Np = ex->periodic ? Nl : Nl + 1;
149: PetscInt Nc, No;
151: PetscFunctionBegin;
152: ex->Nt[ct] = 2;
153: Nc = 6 * Nl;
154: No = 2 * Nl;
155: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
156: ex->target[ct][0] = DM_POLYTOPE_POINT;
157: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
158: ex->size[ct][0] = Np;
159: ex->size[ct][1] = Nl;
160: /* cones for segments/tensor segments */
161: for (PetscInt i = 0; i < Nl; ++i) {
162: ex->cone[ct][6 * i + 0] = DM_POLYTOPE_POINT;
163: ex->cone[ct][6 * i + 1] = 0;
164: ex->cone[ct][6 * i + 2] = i;
165: ex->cone[ct][6 * i + 3] = DM_POLYTOPE_POINT;
166: ex->cone[ct][6 * i + 4] = 0;
167: ex->cone[ct][6 * i + 5] = (i + 1) % Np;
168: }
169: for (PetscInt i = 0; i < No; ++i) ex->ornt[ct][i] = 0;
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /* DM_POLYTOPE_SEGMENT produces
174: Nl+1 segments, or Nl segments periodically and
175: Nl quads, or tensor quads
176: */
177: static PetscErrorCode DMPlexTransformExtrudeSetUp_Segment(DMPlexTransform_Extrude *ex, PetscInt Nl)
178: {
179: const DMPolytopeType ct = DM_POLYTOPE_SEGMENT;
180: const PetscInt Np = ex->periodic ? Nl : Nl + 1;
181: PetscInt Nc, No, coff, ooff;
183: PetscFunctionBegin;
184: ex->Nt[ct] = 2;
185: Nc = 8 * Np + 14 * Nl;
186: No = 2 * Np + 4 * Nl;
187: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
188: ex->target[ct][0] = DM_POLYTOPE_SEGMENT;
189: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
190: ex->size[ct][0] = Np;
191: ex->size[ct][1] = Nl;
192: /* cones for segments */
193: for (PetscInt i = 0; i < Np; ++i) {
194: ex->cone[ct][8 * i + 0] = DM_POLYTOPE_POINT;
195: ex->cone[ct][8 * i + 1] = 1;
196: ex->cone[ct][8 * i + 2] = 0;
197: ex->cone[ct][8 * i + 3] = i;
198: ex->cone[ct][8 * i + 4] = DM_POLYTOPE_POINT;
199: ex->cone[ct][8 * i + 5] = 1;
200: ex->cone[ct][8 * i + 6] = 1;
201: ex->cone[ct][8 * i + 7] = i;
202: }
203: for (PetscInt i = 0; i < 2 * Np; ++i) ex->ornt[ct][i] = 0;
204: /* cones for quads/tensor quads */
205: coff = 8 * Np;
206: ooff = 2 * Np;
207: for (PetscInt i = 0; i < Nl; ++i) {
208: if (ex->useTensor) {
209: ex->cone[ct][coff + 14 * i + 0] = DM_POLYTOPE_SEGMENT;
210: ex->cone[ct][coff + 14 * i + 1] = 0;
211: ex->cone[ct][coff + 14 * i + 2] = i;
212: ex->cone[ct][coff + 14 * i + 3] = DM_POLYTOPE_SEGMENT;
213: ex->cone[ct][coff + 14 * i + 4] = 0;
214: ex->cone[ct][coff + 14 * i + 5] = (i + 1) % Np;
215: ex->cone[ct][coff + 14 * i + 6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
216: ex->cone[ct][coff + 14 * i + 7] = 1;
217: ex->cone[ct][coff + 14 * i + 8] = 0;
218: ex->cone[ct][coff + 14 * i + 9] = i;
219: ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
220: ex->cone[ct][coff + 14 * i + 11] = 1;
221: ex->cone[ct][coff + 14 * i + 12] = 1;
222: ex->cone[ct][coff + 14 * i + 13] = i;
223: ex->ornt[ct][ooff + 4 * i + 0] = 0;
224: ex->ornt[ct][ooff + 4 * i + 1] = 0;
225: ex->ornt[ct][ooff + 4 * i + 2] = 0;
226: ex->ornt[ct][ooff + 4 * i + 3] = 0;
227: } else {
228: ex->cone[ct][coff + 14 * i + 0] = DM_POLYTOPE_SEGMENT;
229: ex->cone[ct][coff + 14 * i + 1] = 0;
230: ex->cone[ct][coff + 14 * i + 2] = i;
231: ex->cone[ct][coff + 14 * i + 3] = DM_POLYTOPE_SEGMENT;
232: ex->cone[ct][coff + 14 * i + 4] = 1;
233: ex->cone[ct][coff + 14 * i + 5] = 1;
234: ex->cone[ct][coff + 14 * i + 6] = i;
235: ex->cone[ct][coff + 14 * i + 7] = DM_POLYTOPE_SEGMENT;
236: ex->cone[ct][coff + 14 * i + 8] = 0;
237: ex->cone[ct][coff + 14 * i + 9] = (i + 1) % Np;
238: ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_SEGMENT;
239: ex->cone[ct][coff + 14 * i + 11] = 1;
240: ex->cone[ct][coff + 14 * i + 12] = 0;
241: ex->cone[ct][coff + 14 * i + 13] = i;
242: ex->ornt[ct][ooff + 4 * i + 0] = 0;
243: ex->ornt[ct][ooff + 4 * i + 1] = 0;
244: ex->ornt[ct][ooff + 4 * i + 2] = -1;
245: ex->ornt[ct][ooff + 4 * i + 3] = -1;
246: }
247: }
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: /* DM_POLYTOPE_TRIANGLE produces
252: Nl+1 triangles, or Nl triangles periodically and
253: Nl triangular prisms/tensor triangular prisms
254: */
255: static PetscErrorCode DMPlexTransformExtrudeSetUp_Triangle(DMPlexTransform_Extrude *ex, PetscInt Nl)
256: {
257: const DMPolytopeType ct = DM_POLYTOPE_TRIANGLE;
258: const PetscInt Np = ex->periodic ? Nl : Nl + 1;
259: PetscInt Nc, No, coff, ooff;
261: PetscFunctionBegin;
262: ex->Nt[ct] = 2;
263: Nc = 12 * Np + 18 * Nl;
264: No = 3 * Np + 5 * Nl;
265: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
266: ex->target[ct][0] = DM_POLYTOPE_TRIANGLE;
267: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM;
268: ex->size[ct][0] = Np;
269: ex->size[ct][1] = Nl;
270: /* cones for triangles */
271: for (PetscInt i = 0; i < Np; ++i) {
272: ex->cone[ct][12 * i + 0] = DM_POLYTOPE_SEGMENT;
273: ex->cone[ct][12 * i + 1] = 1;
274: ex->cone[ct][12 * i + 2] = 0;
275: ex->cone[ct][12 * i + 3] = i;
276: ex->cone[ct][12 * i + 4] = DM_POLYTOPE_SEGMENT;
277: ex->cone[ct][12 * i + 5] = 1;
278: ex->cone[ct][12 * i + 6] = 1;
279: ex->cone[ct][12 * i + 7] = i;
280: ex->cone[ct][12 * i + 8] = DM_POLYTOPE_SEGMENT;
281: ex->cone[ct][12 * i + 9] = 1;
282: ex->cone[ct][12 * i + 10] = 2;
283: ex->cone[ct][12 * i + 11] = i;
284: }
285: for (PetscInt i = 0; i < 3 * Np; ++i) ex->ornt[ct][i] = 0;
286: /* cones for triangular prisms/tensor triangular prisms */
287: coff = 12 * Np;
288: ooff = 3 * Np;
289: for (PetscInt i = 0; i < Nl; ++i) {
290: if (ex->useTensor) {
291: ex->cone[ct][coff + 18 * i + 0] = DM_POLYTOPE_TRIANGLE;
292: ex->cone[ct][coff + 18 * i + 1] = 0;
293: ex->cone[ct][coff + 18 * i + 2] = i;
294: ex->cone[ct][coff + 18 * i + 3] = DM_POLYTOPE_TRIANGLE;
295: ex->cone[ct][coff + 18 * i + 4] = 0;
296: ex->cone[ct][coff + 18 * i + 5] = (i + 1) % Np;
297: ex->cone[ct][coff + 18 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
298: ex->cone[ct][coff + 18 * i + 7] = 1;
299: ex->cone[ct][coff + 18 * i + 8] = 0;
300: ex->cone[ct][coff + 18 * i + 9] = i;
301: ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
302: ex->cone[ct][coff + 18 * i + 11] = 1;
303: ex->cone[ct][coff + 18 * i + 12] = 1;
304: ex->cone[ct][coff + 18 * i + 13] = i;
305: ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
306: ex->cone[ct][coff + 18 * i + 15] = 1;
307: ex->cone[ct][coff + 18 * i + 16] = 2;
308: ex->cone[ct][coff + 18 * i + 17] = i;
309: ex->ornt[ct][ooff + 5 * i + 0] = 0;
310: ex->ornt[ct][ooff + 5 * i + 1] = 0;
311: ex->ornt[ct][ooff + 5 * i + 2] = 0;
312: ex->ornt[ct][ooff + 5 * i + 3] = 0;
313: ex->ornt[ct][ooff + 5 * i + 4] = 0;
314: } else {
315: ex->cone[ct][coff + 18 * i + 0] = DM_POLYTOPE_TRIANGLE;
316: ex->cone[ct][coff + 18 * i + 1] = 0;
317: ex->cone[ct][coff + 18 * i + 2] = i;
318: ex->cone[ct][coff + 18 * i + 3] = DM_POLYTOPE_TRIANGLE;
319: ex->cone[ct][coff + 18 * i + 4] = 0;
320: ex->cone[ct][coff + 18 * i + 5] = (i + 1) % Np;
321: ex->cone[ct][coff + 18 * i + 6] = DM_POLYTOPE_QUADRILATERAL;
322: ex->cone[ct][coff + 18 * i + 7] = 1;
323: ex->cone[ct][coff + 18 * i + 8] = 0;
324: ex->cone[ct][coff + 18 * i + 9] = i;
325: ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
326: ex->cone[ct][coff + 18 * i + 11] = 1;
327: ex->cone[ct][coff + 18 * i + 12] = 1;
328: ex->cone[ct][coff + 18 * i + 13] = i;
329: ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
330: ex->cone[ct][coff + 18 * i + 15] = 1;
331: ex->cone[ct][coff + 18 * i + 16] = 2;
332: ex->cone[ct][coff + 18 * i + 17] = i;
333: ex->ornt[ct][ooff + 5 * i + 0] = -2;
334: ex->ornt[ct][ooff + 5 * i + 1] = 0;
335: ex->ornt[ct][ooff + 5 * i + 2] = 0;
336: ex->ornt[ct][ooff + 5 * i + 3] = 0;
337: ex->ornt[ct][ooff + 5 * i + 4] = 0;
338: }
339: }
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /* DM_POLYTOPE_QUADRILATERAL produces
344: Nl+1 quads, or Nl quads periodically and
345: Nl hexes/tensor hexes
346: */
347: static PetscErrorCode DMPlexTransformExtrudeSetUp_Quadrilateral(DMPlexTransform_Extrude *ex, PetscInt Nl)
348: {
349: const DMPolytopeType ct = DM_POLYTOPE_QUADRILATERAL;
350: const PetscInt Np = ex->periodic ? Nl : Nl + 1;
351: PetscInt Nc, No, coff, ooff;
353: PetscFunctionBegin;
354: ex->Nt[ct] = 2;
355: Nc = 16 * Np + 22 * Nl;
356: No = 4 * Np + 6 * Nl;
357: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
358: ex->target[ct][0] = DM_POLYTOPE_QUADRILATERAL;
359: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON;
360: ex->size[ct][0] = Np;
361: ex->size[ct][1] = Nl;
362: /* cones for quads */
363: for (PetscInt i = 0; i < Np; ++i) {
364: ex->cone[ct][16 * i + 0] = DM_POLYTOPE_SEGMENT;
365: ex->cone[ct][16 * i + 1] = 1;
366: ex->cone[ct][16 * i + 2] = 0;
367: ex->cone[ct][16 * i + 3] = i;
368: ex->cone[ct][16 * i + 4] = DM_POLYTOPE_SEGMENT;
369: ex->cone[ct][16 * i + 5] = 1;
370: ex->cone[ct][16 * i + 6] = 1;
371: ex->cone[ct][16 * i + 7] = i;
372: ex->cone[ct][16 * i + 8] = DM_POLYTOPE_SEGMENT;
373: ex->cone[ct][16 * i + 9] = 1;
374: ex->cone[ct][16 * i + 10] = 2;
375: ex->cone[ct][16 * i + 11] = i;
376: ex->cone[ct][16 * i + 12] = DM_POLYTOPE_SEGMENT;
377: ex->cone[ct][16 * i + 13] = 1;
378: ex->cone[ct][16 * i + 14] = 3;
379: ex->cone[ct][16 * i + 15] = i;
380: }
381: for (PetscInt i = 0; i < 4 * Np; ++i) ex->ornt[ct][i] = 0;
382: /* cones for hexes/tensor hexes */
383: coff = 16 * Np;
384: ooff = 4 * Np;
385: for (PetscInt i = 0; i < Nl; ++i) {
386: if (ex->useTensor) {
387: ex->cone[ct][coff + 22 * i + 0] = DM_POLYTOPE_QUADRILATERAL;
388: ex->cone[ct][coff + 22 * i + 1] = 0;
389: ex->cone[ct][coff + 22 * i + 2] = i;
390: ex->cone[ct][coff + 22 * i + 3] = DM_POLYTOPE_QUADRILATERAL;
391: ex->cone[ct][coff + 22 * i + 4] = 0;
392: ex->cone[ct][coff + 22 * i + 5] = (i + 1) % Np;
393: ex->cone[ct][coff + 22 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
394: ex->cone[ct][coff + 22 * i + 7] = 1;
395: ex->cone[ct][coff + 22 * i + 8] = 0;
396: ex->cone[ct][coff + 22 * i + 9] = i;
397: ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
398: ex->cone[ct][coff + 22 * i + 11] = 1;
399: ex->cone[ct][coff + 22 * i + 12] = 1;
400: ex->cone[ct][coff + 22 * i + 13] = i;
401: ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
402: ex->cone[ct][coff + 22 * i + 15] = 1;
403: ex->cone[ct][coff + 22 * i + 16] = 2;
404: ex->cone[ct][coff + 22 * i + 17] = i;
405: ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
406: ex->cone[ct][coff + 22 * i + 19] = 1;
407: ex->cone[ct][coff + 22 * i + 20] = 3;
408: ex->cone[ct][coff + 22 * i + 21] = i;
409: ex->ornt[ct][ooff + 6 * i + 0] = 0;
410: ex->ornt[ct][ooff + 6 * i + 1] = 0;
411: ex->ornt[ct][ooff + 6 * i + 2] = 0;
412: ex->ornt[ct][ooff + 6 * i + 3] = 0;
413: ex->ornt[ct][ooff + 6 * i + 4] = 0;
414: ex->ornt[ct][ooff + 6 * i + 5] = 0;
415: } else {
416: ex->cone[ct][coff + 22 * i + 0] = DM_POLYTOPE_QUADRILATERAL;
417: ex->cone[ct][coff + 22 * i + 1] = 0;
418: ex->cone[ct][coff + 22 * i + 2] = i;
419: ex->cone[ct][coff + 22 * i + 3] = DM_POLYTOPE_QUADRILATERAL;
420: ex->cone[ct][coff + 22 * i + 4] = 0;
421: ex->cone[ct][coff + 22 * i + 5] = (i + 1) % Np;
422: ex->cone[ct][coff + 22 * i + 6] = DM_POLYTOPE_QUADRILATERAL;
423: ex->cone[ct][coff + 22 * i + 7] = 1;
424: ex->cone[ct][coff + 22 * i + 8] = 0;
425: ex->cone[ct][coff + 22 * i + 9] = i;
426: ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
427: ex->cone[ct][coff + 22 * i + 11] = 1;
428: ex->cone[ct][coff + 22 * i + 12] = 2;
429: ex->cone[ct][coff + 22 * i + 13] = i;
430: ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
431: ex->cone[ct][coff + 22 * i + 15] = 1;
432: ex->cone[ct][coff + 22 * i + 16] = 1;
433: ex->cone[ct][coff + 22 * i + 17] = i;
434: ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_QUADRILATERAL;
435: ex->cone[ct][coff + 22 * i + 19] = 1;
436: ex->cone[ct][coff + 22 * i + 20] = 3;
437: ex->cone[ct][coff + 22 * i + 21] = i;
438: ex->ornt[ct][ooff + 6 * i + 0] = -2;
439: ex->ornt[ct][ooff + 6 * i + 1] = 0;
440: ex->ornt[ct][ooff + 6 * i + 2] = 0;
441: ex->ornt[ct][ooff + 6 * i + 3] = 0;
442: ex->ornt[ct][ooff + 6 * i + 4] = 0;
443: ex->ornt[ct][ooff + 6 * i + 5] = 1;
444: }
445: }
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: /*
450: The refine types for extrusion are:
452: ct: For any normally extruded point
453: ct + 100: For any point which should just return itself
454: */
455: static PetscErrorCode DMPlexTransformSetUp_Extrude(DMPlexTransform tr)
456: {
457: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
458: DM dm;
459: DMLabel active;
460: PetscInt Nl = ex->layers, l, ict, dim;
462: PetscFunctionBegin;
463: PetscCall(DMPlexTransformExtrudeComputeExtrusionDim(tr));
464: PetscCall(DMPlexTransformGetDM(tr, &dm));
465: PetscCall(DMGetDimension(dm, &dim));
466: PetscCall(DMPlexTransformGetActive(tr, &active));
467: if (active) {
468: DMLabel celltype;
469: PetscInt pStart, pEnd, p;
471: PetscCall(DMPlexGetCellTypeLabel(dm, &celltype));
472: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
473: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
474: for (p = pStart; p < pEnd; ++p) {
475: PetscInt ct, val;
477: PetscCall(DMLabelGetValue(celltype, p, &ct));
478: PetscCall(DMLabelGetValue(active, p, &val));
479: if (val < 0) {
480: PetscCall(DMLabelSetValue(tr->trType, p, ct + 100));
481: } else {
482: PetscCall(DMLabelSetValue(tr->trType, p, ct));
483: }
484: }
485: }
486: if (ex->normalAlg != NORMAL_INPUT) {
487: if (dim != ex->cdim) ex->normalAlg = NORMAL_COMPUTE;
488: else if (active) ex->normalAlg = NORMAL_COMPUTE_BD;
489: }
490: // Need this to determine face sharing
491: PetscMPIInt size;
493: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
494: if (size > 1) {
495: PetscSF sf;
496: PetscInt Nr;
498: PetscCall(DMGetPointSF(dm, &sf));
499: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
500: if (Nr >= 0) {
501: PetscCall(PetscSFComputeDegreeBegin(sf, &ex->degree));
502: PetscCall(PetscSFComputeDegreeEnd(sf, &ex->degree));
503: }
504: }
505: // Create normal field
506: if (ex->normalAlg == NORMAL_COMPUTE_BD) {
507: PetscSection s;
508: PetscInt vStart, vEnd;
510: PetscMPIInt rank;
511: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
513: PetscCall(DMClone(dm, &ex->dmNormal));
514: PetscCall(DMGetLocalSection(ex->dmNormal, &s));
515: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
516: PetscCall(PetscSectionSetNumFields(s, 1));
517: PetscCall(PetscSectionSetChart(s, vStart, vEnd));
518: for (PetscInt v = vStart; v < vEnd; ++v) PetscCall(PetscSectionSetDof(s, v, ex->cdimEx));
519: PetscCall(PetscSectionSetUp(s));
520: PetscCall(DMCreateLocalVector(ex->dmNormal, &ex->vecNormal));
521: PetscCall(PetscObjectSetName((PetscObject)ex->vecNormal, "Normal Field"));
523: // find an active point in the closure of v and use its coordinate normal as the extrusion direction
524: PetscSF sf;
525: const PetscInt *leaves;
526: PetscScalar *a, *normal;
527: PetscInt Nl;
529: PetscCall(DMGetPointSF(ex->dmNormal, &sf));
530: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
531: PetscCall(VecGetArrayWrite(ex->vecNormal, &a));
532: for (PetscInt v = vStart; v < vEnd; ++v) {
533: PetscInt *star = NULL;
534: PetscInt starSize, pStart, pEnd;
536: PetscCall(DMPlexGetDepthStratum(ex->dmNormal, ex->cdimEx - 1, &pStart, &pEnd));
537: PetscCall(DMPlexGetTransitiveClosure(ex->dmNormal, v, PETSC_FALSE, &starSize, &star));
538: PetscCall(DMPlexPointLocalRef(ex->dmNormal, v, a, &normal));
539: for (PetscInt st = 0; st < starSize * 2; st += 2) {
540: const PetscInt face = star[st];
541: if ((face >= pStart) && (face < pEnd)) {
542: PetscReal cnormal[3] = {0, 0, 0};
543: const PetscInt *supp;
544: PetscInt suppSize, floc = -1;
545: PetscBool shared;
547: PetscCall(DMPlexComputeCellGeometryFVM(ex->dmNormal, face, NULL, NULL, cnormal));
548: PetscCall(DMPlexGetSupportSize(ex->dmNormal, face, &suppSize));
549: PetscCall(DMPlexGetSupport(ex->dmNormal, face, &supp));
550: // Only use external faces, so I can get the orientation from any cell
551: if (leaves) PetscCall(PetscFindInt(face, Nl, leaves, &floc));
552: shared = floc >= 0 || (ex->degree && ex->degree[face]) ? PETSC_TRUE : PETSC_FALSE;
553: if (suppSize == 1 && !shared) {
554: const PetscInt *cone, *ornt;
555: PetscInt coneSize, c;
557: PetscCall(DMPlexGetConeSize(ex->dmNormal, supp[0], &coneSize));
558: PetscCall(DMPlexGetCone(ex->dmNormal, supp[0], &cone));
559: PetscCall(DMPlexGetConeOrientation(ex->dmNormal, supp[0], &ornt));
560: for (c = 0; c < coneSize; ++c)
561: if (cone[c] == face) break;
562: PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Asymmetry in cone/support");
563: if (ornt[c] < 0)
564: for (PetscInt d = 0; d < ex->cdimEx; ++d) cnormal[d] *= -1.;
565: for (PetscInt d = 0; d < ex->cdimEx; ++d) normal[d] += cnormal[d];
566: }
567: }
568: }
569: PetscCall(DMPlexRestoreTransitiveClosure(ex->dmNormal, v, PETSC_FALSE, &starSize, &star));
570: }
571: PetscCall(VecRestoreArrayWrite(ex->vecNormal, &a));
573: Vec g;
575: PetscCall(DMGetGlobalVector(ex->dmNormal, &g));
576: PetscCall(VecSet(g, 0.));
577: PetscCall(DMLocalToGlobal(ex->dmNormal, ex->vecNormal, ADD_VALUES, g));
578: PetscCall(DMGlobalToLocal(ex->dmNormal, g, INSERT_VALUES, ex->vecNormal));
579: PetscCall(DMRestoreGlobalVector(ex->dmNormal, &g));
580: }
581: PetscCall(PetscMalloc5(DM_NUM_POLYTOPES, &ex->Nt, DM_NUM_POLYTOPES, &ex->target, DM_NUM_POLYTOPES, &ex->size, DM_NUM_POLYTOPES, &ex->cone, DM_NUM_POLYTOPES, &ex->ornt));
582: for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
583: ex->Nt[ict] = -1;
584: ex->target[ict] = NULL;
585: ex->size[ict] = NULL;
586: ex->cone[ict] = NULL;
587: ex->ornt[ict] = NULL;
588: }
589: PetscCall(DMPlexTransformExtrudeSetUp_Point(ex, Nl));
590: PetscCall(DMPlexTransformExtrudeSetUp_Segment(ex, Nl));
591: PetscCall(DMPlexTransformExtrudeSetUp_Triangle(ex, Nl));
592: PetscCall(DMPlexTransformExtrudeSetUp_Quadrilateral(ex, Nl));
593: /* Layers positions */
594: if (!ex->Nth) {
595: if (ex->symmetric)
596: for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers - ex->thickness / 2;
597: else
598: for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers;
599: } else {
600: ex->thickness = 0.;
601: ex->layerPos[0] = 0.;
602: for (l = 0; l < ex->layers; ++l) {
603: const PetscReal t = ex->thicknesses[PetscMin(l, ex->Nth - 1)];
604: ex->layerPos[l + 1] = ex->layerPos[l] + t;
605: ex->thickness += t;
606: }
607: if (ex->symmetric)
608: for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] -= ex->thickness / 2.;
609: }
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr)
614: {
615: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
616: PetscInt ct;
618: PetscFunctionBegin;
619: if (ex->target) {
620: for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) PetscCall(PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct]));
621: }
622: PetscCall(PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt));
623: PetscCall(PetscFree(ex->layerPos));
624: PetscCall(DMDestroy(&ex->dmNormal));
625: PetscCall(VecDestroy(&ex->vecNormal));
626: PetscCall(PetscFree(ex));
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
631: {
632: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
633: DMLabel trType = tr->trType, active;
634: PetscBool onBd = PETSC_FALSE;
635: PetscInt rt;
637: PetscFunctionBeginHot;
638: *rnew = r;
639: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
640: PetscCall(DMPlexTransformGetActive(tr, &active));
641: if (!so && !active) PetscFunctionReturn(PETSC_SUCCESS);
642: if (trType) {
643: PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
644: if (rt >= 100) PetscFunctionReturn(PETSC_SUCCESS);
645: }
646: if (active) {
647: // Get orientation of boundary face in cell
648: if (DMPolytopeTypeGetDim(sct) == ex->dimEx - 1) {
649: DM dm;
650: const PetscInt *supp, *cone, *ornt;
651: PetscInt suppSize, coneSize, c;
653: PetscCall(DMPlexTransformGetDM(tr, &dm));
654: PetscCall(DMPlexGetSupportSize(dm, sp, &suppSize));
655: PetscCall(DMPlexGetSupport(dm, sp, &supp));
656: PetscCheck(suppSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Source point %" PetscInt_FMT " is not a boundary face", sp);
657: PetscCall(DMPlexGetConeSize(dm, supp[0], &coneSize));
658: PetscCall(DMPlexGetOrientedCone(dm, supp[0], &cone, &ornt));
659: for (c = 0; c < coneSize; ++c)
660: if (cone[c] == sp) break;
661: PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Source point %" PetscInt_FMT " not found in cone of support %" PetscInt_FMT, sp, supp[0]);
662: o = ornt[c];
663: onBd = PETSC_TRUE;
664: PetscCall(DMPlexRestoreOrientedCone(dm, supp[0], &cone, &ornt));
665: }
666: }
667: if (ex->useTensor) {
668: switch (sct) {
669: case DM_POLYTOPE_POINT:
670: break;
671: case DM_POLYTOPE_SEGMENT:
672: switch (tct) {
673: case DM_POLYTOPE_SEGMENT:
674: break;
675: case DM_POLYTOPE_SEG_PRISM_TENSOR:
676: *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
677: break;
678: default:
679: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
680: }
681: break;
682: // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
683: case DM_POLYTOPE_TRIANGLE:
684: switch (tct) {
685: case DM_POLYTOPE_TRIANGLE:
686: break;
687: case DM_POLYTOPE_TRI_PRISM_TENSOR:
688: *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
689: break;
690: default:
691: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
692: }
693: break;
694: case DM_POLYTOPE_QUADRILATERAL:
695: switch (tct) {
696: case DM_POLYTOPE_QUADRILATERAL:
697: break;
698: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
699: *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
700: break;
701: default:
702: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
703: }
704: break;
705: default:
706: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
707: }
708: } else {
709: switch (sct) {
710: case DM_POLYTOPE_POINT:
711: break;
712: case DM_POLYTOPE_SEGMENT:
713: switch (tct) {
714: case DM_POLYTOPE_SEGMENT:
715: break;
716: case DM_POLYTOPE_QUADRILATERAL:
717: *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -3) : DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0);
718: break;
719: default:
720: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
721: }
722: break;
723: case DM_POLYTOPE_TRIANGLE:
724: switch (tct) {
725: case DM_POLYTOPE_TRIANGLE:
726: break;
727: case DM_POLYTOPE_TRI_PRISM:
728: *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
729: break;
730: default:
731: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
732: }
733: break;
734: case DM_POLYTOPE_QUADRILATERAL:
735: switch (tct) {
736: case DM_POLYTOPE_QUADRILATERAL:
737: break;
738: case DM_POLYTOPE_HEXAHEDRON:
739: *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
740: break;
741: default:
742: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
743: }
744: break;
745: default:
746: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
747: }
748: }
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
753: {
754: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
755: DMLabel trType = tr->trType;
756: PetscBool ignore = PETSC_FALSE, identity = PETSC_FALSE;
757: PetscInt val = 0;
759: PetscFunctionBegin;
760: if (trType) {
761: PetscCall(DMLabelGetValue(trType, p, &val));
762: identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE;
763: } else {
764: ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE;
765: }
766: if (rt) *rt = val;
767: if (ignore) {
768: /* Ignore cells that cannot be extruded */
769: *Nt = 0;
770: *target = NULL;
771: *size = NULL;
772: *cone = NULL;
773: *ornt = NULL;
774: } else if (identity) {
775: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
776: } else {
777: *Nt = ex->Nt[source];
778: *target = ex->target[source];
779: *size = ex->size[source];
780: *cone = ex->cone[source];
781: *ornt = ex->ornt[source];
782: }
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }
786: /* Computes new vertex along normal */
787: static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
788: {
789: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
790: DM dm;
791: DMLabel active;
792: PetscReal ones2[2] = {0., 1.};
793: PetscReal ones3[3] = {0., 0., 1.};
794: PetscReal normal[3] = {0., 0., 0.};
795: PetscReal norm = 0.;
796: PetscInt dEx = ex->cdimEx;
797: PetscInt dim, cStart, cEnd;
799: PetscFunctionBeginHot;
800: PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]);
801: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
802: PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv);
803: PetscCheck(dE == ex->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %" PetscInt_FMT " != %" PetscInt_FMT " original dimension", dE, ex->cdim);
804: PetscCheck(dEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", dEx);
806: PetscCall(DMPlexTransformGetDM(tr, &dm));
807: PetscCall(DMGetDimension(dm, &dim));
808: PetscCall(DMPlexTransformGetActive(tr, &active));
809: switch (ex->normalAlg) {
810: case NORMAL_DEFAULT:
811: switch (ex->cdimEx) {
812: case 2:
813: for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones2[d];
814: break;
815: case 3:
816: for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones3[d];
817: break;
818: default:
819: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No default normal for dimension %" PetscInt_FMT, ex->cdimEx);
820: }
821: break;
822: case NORMAL_INPUT:
823: for (PetscInt d = 0; d < dEx; ++d) normal[d] = ex->normal[d];
824: break;
825: case NORMAL_COMPUTE: {
826: PetscInt *star = NULL;
827: PetscInt starSize;
829: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
830: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star));
831: for (PetscInt st = 0; st < starSize * 2; st += 2) {
832: if ((star[st] >= cStart) && (star[st] < cEnd)) {
833: PetscReal cnormal[3] = {0, 0, 0};
835: PetscCall(DMPlexComputeCellGeometryFVM(dm, star[st], NULL, NULL, cnormal));
836: for (PetscInt d = 0; d < dEx; ++d) normal[d] += cnormal[d];
837: }
838: }
839: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star));
840: } break;
841: case NORMAL_COMPUTE_BD: {
842: const PetscScalar *a;
843: PetscScalar *vnormal;
845: PetscCall(VecGetArrayRead(ex->vecNormal, &a));
846: PetscCall(DMPlexPointLocalRead(ex->dmNormal, p, a, (void *)&vnormal));
847: for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscRealPart(vnormal[d]);
848: PetscCall(VecRestoreArrayRead(ex->vecNormal, &a));
849: } break;
850: default:
851: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
852: }
853: if (ex->normalFunc) {
854: PetscScalar n[3];
855: PetscReal x[3], dot = 0.;
857: for (PetscInt d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]);
858: PetscCall((*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL));
859: for (PetscInt d = 0; d < dEx; ++d) dot += PetscRealPart(n[d]) * normal[d];
860: for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]);
861: }
862: for (PetscInt d = 0; d < dEx; ++d) norm += PetscSqr(normal[d]);
863: for (PetscInt d = 0; d < dEx; ++d) normal[d] *= norm == 0.0 ? 1.0 : 1. / PetscSqrtReal(norm);
864: for (PetscInt d = 0; d < dEx; ++d) out[d] = normal[d] * ex->layerPos[r];
865: for (PetscInt d = 0; d < ex->cdim; ++d) out[d] += in[d];
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr)
870: {
871: PetscFunctionBegin;
872: tr->ops->view = DMPlexTransformView_Extrude;
873: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude;
874: tr->ops->setup = DMPlexTransformSetUp_Extrude;
875: tr->ops->destroy = DMPlexTransformDestroy_Extrude;
876: tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude;
877: tr->ops->celltransform = DMPlexTransformCellTransform_Extrude;
878: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude;
879: tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude;
880: PetscFunctionReturn(PETSC_SUCCESS);
881: }
883: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr)
884: {
885: DMPlexTransform_Extrude *ex;
886: DM dm;
887: PetscInt dim;
889: PetscFunctionBegin;
891: PetscCall(PetscNew(&ex));
892: tr->data = ex;
893: ex->thickness = 1.;
894: ex->useTensor = PETSC_TRUE;
895: ex->symmetric = PETSC_FALSE;
896: ex->periodic = PETSC_FALSE;
897: ex->normalAlg = NORMAL_DEFAULT;
898: ex->layerPos = NULL;
899: PetscCall(DMPlexTransformGetDM(tr, &dm));
900: PetscCall(DMGetDimension(dm, &dim));
901: PetscCall(DMGetCoordinateDim(dm, &ex->cdim));
902: PetscCall(DMPlexTransformExtrudeSetLayers(tr, 1));
903: PetscCall(DMPlexTransformInitialize_Extrude(tr));
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
907: /*@
908: DMPlexTransformExtrudeGetLayers - Get the number of extruded layers.
910: Not Collective
912: Input Parameter:
913: . tr - The `DMPlexTransform`
915: Output Parameter:
916: . layers - The number of layers
918: Level: intermediate
920: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetLayers()`
921: @*/
922: PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers)
923: {
924: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
926: PetscFunctionBegin;
928: PetscAssertPointer(layers, 2);
929: *layers = ex->layers;
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: /*@
934: DMPlexTransformExtrudeSetLayers - Set the number of extruded layers.
936: Not Collective
938: Input Parameters:
939: + tr - The `DMPlexTransform`
940: - layers - The number of layers
942: Level: intermediate
944: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetLayers()`
945: @*/
946: PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers)
947: {
948: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
950: PetscFunctionBegin;
952: ex->layers = layers;
953: PetscCall(PetscFree(ex->layerPos));
954: PetscCall(PetscCalloc1(ex->layers + 1, &ex->layerPos));
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: /*@
959: DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers
961: Not Collective
963: Input Parameter:
964: . tr - The `DMPlexTransform`
966: Output Parameter:
967: . thickness - The total thickness of the layers
969: Level: intermediate
971: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`
972: @*/
973: PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness)
974: {
975: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
977: PetscFunctionBegin;
979: PetscAssertPointer(thickness, 2);
980: *thickness = ex->thickness;
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }
984: /*@
985: DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers
987: Not Collective
989: Input Parameters:
990: + tr - The `DMPlexTransform`
991: - thickness - The total thickness of the layers
993: Level: intermediate
995: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetThickness()`
996: @*/
997: PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness)
998: {
999: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1001: PetscFunctionBegin;
1003: PetscCheck(thickness > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double)thickness);
1004: ex->thickness = thickness;
1005: PetscFunctionReturn(PETSC_SUCCESS);
1006: }
1008: /*@
1009: DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells
1011: Not Collective
1013: Input Parameter:
1014: . tr - The `DMPlexTransform`
1016: Output Parameter:
1017: . useTensor - The flag to use tensor cells
1019: Note:
1020: This flag determines the orientation behavior of the created points.
1022: For example, if tensor is `PETSC_TRUE`, then
1023: .vb
1024: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1025: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1026: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1027: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1028: .ve
1030: Level: intermediate
1032: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetTensor()`
1033: @*/
1034: PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
1035: {
1036: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1038: PetscFunctionBegin;
1040: PetscAssertPointer(useTensor, 2);
1041: *useTensor = ex->useTensor;
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: /*@
1046: DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells
1048: Not Collective
1050: Input Parameters:
1051: + tr - The `DMPlexTransform`
1052: - useTensor - The flag for tensor cells
1054: Note:
1055: This flag determines the orientation behavior of the created points
1056: For example, if tensor is `PETSC_TRUE`, then
1057: .vb
1058: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1059: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1060: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1061: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1062: .ve
1064: Level: intermediate
1066: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetTensor()`
1067: @*/
1068: PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
1069: {
1070: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1072: PetscFunctionBegin;
1074: ex->useTensor = useTensor;
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: /*@
1079: DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface
1081: Not Collective
1083: Input Parameter:
1084: . tr - The `DMPlexTransform`
1086: Output Parameter:
1087: . symmetric - The flag to extrude symmetrically
1089: Level: intermediate
1091: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetSymmetric()`, `DMPlexTransformExtrudeGetPeriodic()`
1092: @*/
1093: PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric)
1094: {
1095: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1097: PetscFunctionBegin;
1099: PetscAssertPointer(symmetric, 2);
1100: *symmetric = ex->symmetric;
1101: PetscFunctionReturn(PETSC_SUCCESS);
1102: }
1104: /*@
1105: DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface
1107: Not Collective
1109: Input Parameters:
1110: + tr - The `DMPlexTransform`
1111: - symmetric - The flag to extrude symmetrically
1113: Level: intermediate
1115: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetSymmetric()`, `DMPlexTransformExtrudeSetPeriodic()`
1116: @*/
1117: PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric)
1118: {
1119: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1121: PetscFunctionBegin;
1123: ex->symmetric = symmetric;
1124: PetscFunctionReturn(PETSC_SUCCESS);
1125: }
1127: /*@
1128: DMPlexTransformExtrudeGetPeriodic - Get the flag to extrude periodically from the initial surface
1130: Not Collective
1132: Input Parameter:
1133: . tr - The `DMPlexTransform`
1135: Output Parameter:
1136: . periodic - The flag to extrude periodically
1138: Level: intermediate
1140: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetPeriodic()`, `DMPlexTransformExtrudeGetSymmetric()`
1141: @*/
1142: PetscErrorCode DMPlexTransformExtrudeGetPeriodic(DMPlexTransform tr, PetscBool *periodic)
1143: {
1144: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1146: PetscFunctionBegin;
1148: PetscAssertPointer(periodic, 2);
1149: *periodic = ex->periodic;
1150: PetscFunctionReturn(PETSC_SUCCESS);
1151: }
1153: /*@
1154: DMPlexTransformExtrudeSetPeriodic - Set the flag to extrude periodically from the initial surface
1156: Not Collective
1158: Input Parameters:
1159: + tr - The `DMPlexTransform`
1160: - periodic - The flag to extrude periodically
1162: Level: intermediate
1164: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetPeriodic()`, `DMPlexTransformExtrudeSetSymmetric()`
1165: @*/
1166: PetscErrorCode DMPlexTransformExtrudeSetPeriodic(DMPlexTransform tr, PetscBool periodic)
1167: {
1168: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1170: PetscFunctionBegin;
1172: ex->periodic = periodic;
1173: PetscFunctionReturn(PETSC_SUCCESS);
1174: }
1176: /*@
1177: DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector
1179: Not Collective
1181: Input Parameter:
1182: . tr - The `DMPlexTransform`
1184: Output Parameter:
1185: . normal - The extrusion direction
1187: Note:
1188: The user passes in an array, which is filled by the library.
1190: Level: intermediate
1192: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetNormal()`
1193: @*/
1194: PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[])
1195: {
1196: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1197: PetscInt d;
1199: PetscFunctionBegin;
1201: if (ex->normalAlg == NORMAL_INPUT) {
1202: for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d];
1203: } else {
1204: for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.;
1205: }
1206: PetscFunctionReturn(PETSC_SUCCESS);
1207: }
1209: /*@
1210: DMPlexTransformExtrudeSetNormal - Set the extrusion normal
1212: Not Collective
1214: Input Parameters:
1215: + tr - The `DMPlexTransform`
1216: - normal - The extrusion direction
1218: Level: intermediate
1220: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`
1221: @*/
1222: PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[])
1223: {
1224: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1225: PetscInt d;
1227: PetscFunctionBegin;
1229: ex->normalAlg = NORMAL_INPUT;
1230: if (!ex->cdimEx) PetscCall(DMPlexTransformExtrudeComputeExtrusionDim(tr));
1231: for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d];
1232: PetscFunctionReturn(PETSC_SUCCESS);
1233: }
1235: /*@C
1236: DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal
1238: Not Collective
1240: Input Parameters:
1241: + tr - The `DMPlexTransform`
1242: - normalFunc - A function determining the extrusion direction, see `PetscSimplePointFn` for the calling sequence
1244: Level: intermediate
1246: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`, `PetscSimplePointFn`
1247: @*/
1248: PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFn *normalFunc)
1249: {
1250: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1252: PetscFunctionBegin;
1254: ex->normalFunc = normalFunc;
1255: PetscFunctionReturn(PETSC_SUCCESS);
1256: }
1258: /*@
1259: DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer
1261: Not Collective
1263: Input Parameters:
1264: + tr - The `DMPlexTransform`
1265: . Nth - The number of thicknesses
1266: - thicknesses - The array of thicknesses
1268: Level: intermediate
1270: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeGetThickness()`
1271: @*/
1272: PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[])
1273: {
1274: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1275: PetscInt t;
1277: PetscFunctionBegin;
1279: PetscCheck(Nth > 0, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %" PetscInt_FMT " must be positive", Nth);
1280: ex->Nth = PetscMin(Nth, ex->layers);
1281: PetscCall(PetscFree(ex->thicknesses));
1282: PetscCall(PetscMalloc1(ex->Nth, &ex->thicknesses));
1283: for (t = 0; t < ex->Nth; ++t) {
1284: PetscCheck(thicknesses[t] > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %" PetscInt_FMT " must be positive", (double)thicknesses[t], t);
1285: ex->thicknesses[t] = thicknesses[t];
1286: }
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }