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;
634: PetscInt rt;
636: PetscFunctionBeginHot;
637: *rnew = r;
638: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
639: if (!so) PetscFunctionReturn(PETSC_SUCCESS);
640: if (trType) {
641: PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
642: if (rt >= 100) PetscFunctionReturn(PETSC_SUCCESS);
643: }
644: if (ex->useTensor) {
645: switch (sct) {
646: case DM_POLYTOPE_POINT:
647: break;
648: case DM_POLYTOPE_SEGMENT:
649: switch (tct) {
650: case DM_POLYTOPE_SEGMENT:
651: break;
652: case DM_POLYTOPE_SEG_PRISM_TENSOR:
653: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
654: break;
655: default:
656: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
657: }
658: break;
659: // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
660: case DM_POLYTOPE_TRIANGLE:
661: break;
662: case DM_POLYTOPE_QUADRILATERAL:
663: break;
664: default:
665: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
666: }
667: } else {
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_QUADRILATERAL:
676: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 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: break;
685: case DM_POLYTOPE_QUADRILATERAL:
686: break;
687: default:
688: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
689: }
690: }
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
695: {
696: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
697: DMLabel trType = tr->trType;
698: PetscBool ignore = PETSC_FALSE, identity = PETSC_FALSE;
699: PetscInt val = 0;
701: PetscFunctionBegin;
702: if (trType) {
703: PetscCall(DMLabelGetValue(trType, p, &val));
704: identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE;
705: } else {
706: ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE;
707: }
708: if (rt) *rt = val;
709: if (ignore) {
710: /* Ignore cells that cannot be extruded */
711: *Nt = 0;
712: *target = NULL;
713: *size = NULL;
714: *cone = NULL;
715: *ornt = NULL;
716: } else if (identity) {
717: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
718: } else {
719: *Nt = ex->Nt[source];
720: *target = ex->target[source];
721: *size = ex->size[source];
722: *cone = ex->cone[source];
723: *ornt = ex->ornt[source];
724: }
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: /* Computes new vertex along normal */
729: static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
730: {
731: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
732: DM dm;
733: DMLabel active;
734: PetscReal ones2[2] = {0., 1.};
735: PetscReal ones3[3] = {0., 0., 1.};
736: PetscReal normal[3] = {0., 0., 0.};
737: PetscReal norm = 0.;
738: PetscInt dEx = ex->cdimEx;
739: PetscInt dim, cStart, cEnd;
741: PetscFunctionBeginHot;
742: PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]);
743: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
744: PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv);
745: PetscCheck(dE == ex->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %" PetscInt_FMT " != %" PetscInt_FMT " original dimension", dE, ex->cdim);
746: PetscCheck(dEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", dEx);
748: PetscCall(DMPlexTransformGetDM(tr, &dm));
749: PetscCall(DMGetDimension(dm, &dim));
750: PetscCall(DMPlexTransformGetActive(tr, &active));
751: switch (ex->normalAlg) {
752: case NORMAL_DEFAULT:
753: switch (ex->cdimEx) {
754: case 2:
755: for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones2[d];
756: break;
757: case 3:
758: for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones3[d];
759: break;
760: default:
761: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No default normal for dimension %" PetscInt_FMT, ex->cdimEx);
762: }
763: break;
764: case NORMAL_INPUT:
765: for (PetscInt d = 0; d < dEx; ++d) normal[d] = ex->normal[d];
766: break;
767: case NORMAL_COMPUTE: {
768: PetscInt *star = NULL;
769: PetscInt starSize;
771: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
772: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star));
773: for (PetscInt st = 0; st < starSize * 2; st += 2) {
774: if ((star[st] >= cStart) && (star[st] < cEnd)) {
775: PetscReal cnormal[3] = {0, 0, 0};
777: PetscCall(DMPlexComputeCellGeometryFVM(dm, star[st], NULL, NULL, cnormal));
778: for (PetscInt d = 0; d < dEx; ++d) normal[d] += cnormal[d];
779: }
780: }
781: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star));
782: } break;
783: case NORMAL_COMPUTE_BD: {
784: const PetscScalar *a;
785: PetscScalar *vnormal;
787: PetscCall(VecGetArrayRead(ex->vecNormal, &a));
788: PetscCall(DMPlexPointLocalRead(ex->dmNormal, p, a, (void *)&vnormal));
789: for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscRealPart(vnormal[d]);
790: PetscCall(VecRestoreArrayRead(ex->vecNormal, &a));
791: } break;
792: default:
793: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
794: }
795: if (ex->normalFunc) {
796: PetscScalar n[3];
797: PetscReal x[3], dot = 0.;
799: for (PetscInt d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]);
800: PetscCall((*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL));
801: for (PetscInt d = 0; d < dEx; ++d) dot += PetscRealPart(n[d]) * normal[d];
802: for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]);
803: }
804: for (PetscInt d = 0; d < dEx; ++d) norm += PetscSqr(normal[d]);
805: for (PetscInt d = 0; d < dEx; ++d) normal[d] *= norm == 0.0 ? 1.0 : 1. / PetscSqrtReal(norm);
806: for (PetscInt d = 0; d < dEx; ++d) out[d] = normal[d] * ex->layerPos[r];
807: for (PetscInt d = 0; d < ex->cdim; ++d) out[d] += in[d];
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr)
812: {
813: PetscFunctionBegin;
814: tr->ops->view = DMPlexTransformView_Extrude;
815: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude;
816: tr->ops->setup = DMPlexTransformSetUp_Extrude;
817: tr->ops->destroy = DMPlexTransformDestroy_Extrude;
818: tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude;
819: tr->ops->celltransform = DMPlexTransformCellTransform_Extrude;
820: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude;
821: tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude;
822: PetscFunctionReturn(PETSC_SUCCESS);
823: }
825: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr)
826: {
827: DMPlexTransform_Extrude *ex;
828: DM dm;
829: PetscInt dim;
831: PetscFunctionBegin;
833: PetscCall(PetscNew(&ex));
834: tr->data = ex;
835: ex->thickness = 1.;
836: ex->useTensor = PETSC_TRUE;
837: ex->symmetric = PETSC_FALSE;
838: ex->periodic = PETSC_FALSE;
839: ex->normalAlg = NORMAL_DEFAULT;
840: ex->layerPos = NULL;
841: PetscCall(DMPlexTransformGetDM(tr, &dm));
842: PetscCall(DMGetDimension(dm, &dim));
843: PetscCall(DMGetCoordinateDim(dm, &ex->cdim));
844: PetscCall(DMPlexTransformExtrudeSetLayers(tr, 1));
845: PetscCall(DMPlexTransformInitialize_Extrude(tr));
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: /*@
850: DMPlexTransformExtrudeGetLayers - Get the number of extruded layers.
852: Not Collective
854: Input Parameter:
855: . tr - The `DMPlexTransform`
857: Output Parameter:
858: . layers - The number of layers
860: Level: intermediate
862: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetLayers()`
863: @*/
864: PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers)
865: {
866: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
868: PetscFunctionBegin;
870: PetscAssertPointer(layers, 2);
871: *layers = ex->layers;
872: PetscFunctionReturn(PETSC_SUCCESS);
873: }
875: /*@
876: DMPlexTransformExtrudeSetLayers - Set the number of extruded layers.
878: Not Collective
880: Input Parameters:
881: + tr - The `DMPlexTransform`
882: - layers - The number of layers
884: Level: intermediate
886: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetLayers()`
887: @*/
888: PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers)
889: {
890: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
892: PetscFunctionBegin;
894: ex->layers = layers;
895: PetscCall(PetscFree(ex->layerPos));
896: PetscCall(PetscCalloc1(ex->layers + 1, &ex->layerPos));
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: /*@
901: DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers
903: Not Collective
905: Input Parameter:
906: . tr - The `DMPlexTransform`
908: Output Parameter:
909: . thickness - The total thickness of the layers
911: Level: intermediate
913: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`
914: @*/
915: PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness)
916: {
917: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
919: PetscFunctionBegin;
921: PetscAssertPointer(thickness, 2);
922: *thickness = ex->thickness;
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: /*@
927: DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers
929: Not Collective
931: Input Parameters:
932: + tr - The `DMPlexTransform`
933: - thickness - The total thickness of the layers
935: Level: intermediate
937: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetThickness()`
938: @*/
939: PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness)
940: {
941: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
943: PetscFunctionBegin;
945: PetscCheck(thickness > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double)thickness);
946: ex->thickness = thickness;
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
950: /*@
951: DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells
953: Not Collective
955: Input Parameter:
956: . tr - The `DMPlexTransform`
958: Output Parameter:
959: . useTensor - The flag to use tensor cells
961: Note:
962: This flag determines the orientation behavior of the created points.
964: For example, if tensor is `PETSC_TRUE`, then
965: .vb
966: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
967: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
968: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
969: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
970: .ve
972: Level: intermediate
974: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetTensor()`
975: @*/
976: PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
977: {
978: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
980: PetscFunctionBegin;
982: PetscAssertPointer(useTensor, 2);
983: *useTensor = ex->useTensor;
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: /*@
988: DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells
990: Not Collective
992: Input Parameters:
993: + tr - The `DMPlexTransform`
994: - useTensor - The flag for tensor cells
996: Note:
997: This flag determines the orientation behavior of the created points
998: For example, if tensor is `PETSC_TRUE`, then
999: .vb
1000: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1001: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1002: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1003: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1004: .ve
1006: Level: intermediate
1008: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetTensor()`
1009: @*/
1010: PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
1011: {
1012: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1014: PetscFunctionBegin;
1016: ex->useTensor = useTensor;
1017: PetscFunctionReturn(PETSC_SUCCESS);
1018: }
1020: /*@
1021: DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface
1023: Not Collective
1025: Input Parameter:
1026: . tr - The `DMPlexTransform`
1028: Output Parameter:
1029: . symmetric - The flag to extrude symmetrically
1031: Level: intermediate
1033: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetSymmetric()`, `DMPlexTransformExtrudeGetPeriodic()`
1034: @*/
1035: PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric)
1036: {
1037: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1039: PetscFunctionBegin;
1041: PetscAssertPointer(symmetric, 2);
1042: *symmetric = ex->symmetric;
1043: PetscFunctionReturn(PETSC_SUCCESS);
1044: }
1046: /*@
1047: DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface
1049: Not Collective
1051: Input Parameters:
1052: + tr - The `DMPlexTransform`
1053: - symmetric - The flag to extrude symmetrically
1055: Level: intermediate
1057: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetSymmetric()`, `DMPlexTransformExtrudeSetPeriodic()`
1058: @*/
1059: PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric)
1060: {
1061: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1063: PetscFunctionBegin;
1065: ex->symmetric = symmetric;
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }
1069: /*@
1070: DMPlexTransformExtrudeGetPeriodic - Get the flag to extrude periodically from the initial surface
1072: Not Collective
1074: Input Parameter:
1075: . tr - The `DMPlexTransform`
1077: Output Parameter:
1078: . periodic - The flag to extrude periodically
1080: Level: intermediate
1082: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetPeriodic()`, `DMPlexTransformExtrudeGetSymmetric()`
1083: @*/
1084: PetscErrorCode DMPlexTransformExtrudeGetPeriodic(DMPlexTransform tr, PetscBool *periodic)
1085: {
1086: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1088: PetscFunctionBegin;
1090: PetscAssertPointer(periodic, 2);
1091: *periodic = ex->periodic;
1092: PetscFunctionReturn(PETSC_SUCCESS);
1093: }
1095: /*@
1096: DMPlexTransformExtrudeSetPeriodic - Set the flag to extrude periodically from the initial surface
1098: Not Collective
1100: Input Parameters:
1101: + tr - The `DMPlexTransform`
1102: - periodic - The flag to extrude periodically
1104: Level: intermediate
1106: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetPeriodic()`, `DMPlexTransformExtrudeSetSymmetric()`
1107: @*/
1108: PetscErrorCode DMPlexTransformExtrudeSetPeriodic(DMPlexTransform tr, PetscBool periodic)
1109: {
1110: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1112: PetscFunctionBegin;
1114: ex->periodic = periodic;
1115: PetscFunctionReturn(PETSC_SUCCESS);
1116: }
1118: /*@
1119: DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector
1121: Not Collective
1123: Input Parameter:
1124: . tr - The `DMPlexTransform`
1126: Output Parameter:
1127: . normal - The extrusion direction
1129: Note:
1130: The user passes in an array, which is filled by the library.
1132: Level: intermediate
1134: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetNormal()`
1135: @*/
1136: PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[])
1137: {
1138: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1139: PetscInt d;
1141: PetscFunctionBegin;
1143: if (ex->normalAlg == NORMAL_INPUT) {
1144: for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d];
1145: } else {
1146: for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.;
1147: }
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: /*@
1152: DMPlexTransformExtrudeSetNormal - Set the extrusion normal
1154: Not Collective
1156: Input Parameters:
1157: + tr - The `DMPlexTransform`
1158: - normal - The extrusion direction
1160: Level: intermediate
1162: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`
1163: @*/
1164: PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[])
1165: {
1166: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1167: PetscInt d;
1169: PetscFunctionBegin;
1171: ex->normalAlg = NORMAL_INPUT;
1172: for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d];
1173: PetscFunctionReturn(PETSC_SUCCESS);
1174: }
1176: /*@C
1177: DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal
1179: Not Collective
1181: Input Parameters:
1182: + tr - The `DMPlexTransform`
1183: - normalFunc - A function determining the extrusion direction, see `PetscSimplePointFn` for the calling sequence
1185: Level: intermediate
1187: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`, `PetscSimplePointFn`
1188: @*/
1189: PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFn *normalFunc)
1190: {
1191: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1193: PetscFunctionBegin;
1195: ex->normalFunc = normalFunc;
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: /*@
1200: DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer
1202: Not Collective
1204: Input Parameters:
1205: + tr - The `DMPlexTransform`
1206: . Nth - The number of thicknesses
1207: - thicknesses - The array of thicknesses
1209: Level: intermediate
1211: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeGetThickness()`
1212: @*/
1213: PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[])
1214: {
1215: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1216: PetscInt t;
1218: PetscFunctionBegin;
1220: PetscCheck(Nth > 0, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %" PetscInt_FMT " must be positive", Nth);
1221: ex->Nth = PetscMin(Nth, ex->layers);
1222: PetscCall(PetscFree(ex->thicknesses));
1223: PetscCall(PetscMalloc1(ex->Nth, &ex->thicknesses));
1224: for (t = 0; t < ex->Nth; ++t) {
1225: PetscCheck(thicknesses[t] > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %" PetscInt_FMT " must be positive", (double)thicknesses[t], t);
1226: ex->thicknesses[t] = thicknesses[t];
1227: }
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }