Actual source code: plextrextrude.c
1: #include <petsc/private/dmplextransformimpl.h>
3: static PetscErrorCode DMPlexTransformView_Extrude(DMPlexTransform tr, PetscViewer viewer)
4: {
5: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
6: PetscBool isascii;
8: PetscFunctionBegin;
11: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
12: if (isascii) {
13: const char *name;
15: PetscCall(PetscObjectGetName((PetscObject)tr, &name));
16: PetscCall(PetscViewerASCIIPrintf(viewer, "Extrusion transformation %s\n", name ? name : ""));
17: PetscCall(PetscViewerASCIIPrintf(viewer, " number of layers: %" PetscInt_FMT "\n", ex->layers));
18: PetscCall(PetscViewerASCIIPrintf(viewer, " create tensor cells: %s\n", ex->useTensor ? "YES" : "NO"));
19: if (ex->periodic) PetscCall(PetscViewerASCIIPrintf(viewer, " periodic\n"));
20: } else {
21: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
22: }
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
26: static PetscErrorCode DMPlexTransformSetFromOptions_Extrude(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
27: {
28: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
29: PetscReal th, normal[3], *thicknesses;
30: PetscInt nl, Nc;
31: PetscBool tensor, sym, per, flg;
32: char funcname[PETSC_MAX_PATH_LEN];
34: PetscFunctionBegin;
35: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Extrusion Options");
36: PetscCall(PetscOptionsBoundedInt("-dm_plex_transform_extrude_layers", "Number of layers to extrude", "", ex->layers, &nl, &flg, 1));
37: if (flg) PetscCall(DMPlexTransformExtrudeSetLayers(tr, nl));
38: PetscCall(PetscOptionsReal("-dm_plex_transform_extrude_thickness", "Total thickness of extruded layers", "", ex->thickness, &th, &flg));
39: if (flg) PetscCall(DMPlexTransformExtrudeSetThickness(tr, th));
40: PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg));
41: if (flg) PetscCall(DMPlexTransformExtrudeSetTensor(tr, tensor));
42: PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_symmetric", "Extrude layers symmetrically about the surface", "", ex->symmetric, &sym, &flg));
43: if (flg) PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, sym));
44: PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_periodic", "Extrude layers periodically about the surface", "", ex->periodic, &per, &flg));
45: if (flg) PetscCall(DMPlexTransformExtrudeSetPeriodic(tr, per));
46: Nc = 3;
47: PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_normal", "Input normal vector for extrusion", "DMPlexTransformExtrudeSetNormal", normal, &Nc, &flg));
48: if (flg) {
49: // Extrusion dimension might not yet be determined
50: 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);
51: PetscCall(DMPlexTransformExtrudeSetNormal(tr, normal));
52: }
53: PetscCall(PetscOptionsString("-dm_plex_transform_extrude_normal_function", "Function to determine normal vector", "DMPlexTransformExtrudeSetNormalFunction", funcname, funcname, sizeof(funcname), &flg));
54: if (flg) {
55: PetscSimplePointFn *normalFunc;
57: PetscCall(PetscDLSym(NULL, funcname, (void **)&normalFunc));
58: PetscCall(DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc));
59: }
60: nl = ex->layers;
61: PetscCall(PetscMalloc1(nl, &thicknesses));
62: PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_thicknesses", "Thickness of each individual extruded layer", "", thicknesses, &nl, &flg));
63: if (flg) {
64: PetscCheck(nl, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one thickness for -dm_plex_transform_extrude_thicknesses");
65: PetscCall(DMPlexTransformExtrudeSetThicknesses(tr, nl, thicknesses));
66: }
67: PetscCall(PetscFree(thicknesses));
68: PetscOptionsHeadEnd();
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: /* Determine the implicit dimension pre-extrusion (either the implicit dimension of the DM or of a point in the active set for the transform).
73: If that dimension is the same as the current coordinate dimension (ex->dim), the extruded mesh will have a coordinate dimension one greater;
74: Otherwise the coordinate dimension will be kept. */
75: static PetscErrorCode DMPlexTransformExtrudeComputeExtrusionDim(DMPlexTransform tr)
76: {
77: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
78: DM dm;
79: DMLabel active;
80: PetscInt dim, dimExtPoint, dimExtPointG;
82: PetscFunctionBegin;
83: PetscCall(DMPlexTransformGetDM(tr, &dm));
84: PetscCall(DMGetDimension(dm, &dim));
85: PetscCall(DMPlexTransformGetActive(tr, &active));
86: if (active) {
87: IS valueIS, pointIS;
88: const PetscInt *values, *points;
89: DMPolytopeType ct;
90: PetscInt Nv, Np;
92: dimExtPoint = 0;
93: PetscCall(DMLabelGetValueIS(active, &valueIS));
94: PetscCall(ISGetLocalSize(valueIS, &Nv));
95: PetscCall(ISGetIndices(valueIS, &values));
96: for (PetscInt v = 0; v < Nv; ++v) {
97: PetscCall(DMLabelGetStratumIS(active, values[v], &pointIS));
98: PetscCall(ISGetLocalSize(pointIS, &Np));
99: PetscCall(ISGetIndices(pointIS, &points));
100: for (PetscInt p = 0; p < Np; ++p) {
101: PetscCall(DMPlexGetCellType(dm, points[p], &ct));
102: dimExtPoint = PetscMax(dimExtPoint, DMPolytopeTypeGetDim(ct));
103: }
104: PetscCall(ISRestoreIndices(pointIS, &points));
105: PetscCall(ISDestroy(&pointIS));
106: }
107: PetscCall(ISRestoreIndices(valueIS, &values));
108: PetscCall(ISDestroy(&valueIS));
109: } else dimExtPoint = dim;
110: PetscCall(MPIU_Allreduce(&dimExtPoint, &dimExtPointG, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)tr)));
111: ex->dimEx = PetscMax(dim, dimExtPointG + 1);
112: ex->cdimEx = ex->cdim == dimExtPointG ? ex->cdim + 1 : ex->cdim;
113: PetscCheck(ex->dimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Topological dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->dimEx);
114: PetscCheck(ex->cdimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->cdimEx);
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm)
119: {
120: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
121: PetscInt dim;
123: PetscFunctionBegin;
124: PetscCall(DMGetDimension(dm, &dim));
125: PetscCall(DMSetDimension(tdm, ex->dimEx));
126: PetscCall(DMSetCoordinateDim(tdm, ex->cdimEx));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: /* DM_POLYTOPE_POINT produces
131: Nl+1 points, or Nl points periodically and
132: Nl segments, or tensor segments
133: */
134: static PetscErrorCode DMPlexTransformExtrudeSetUp_Point(DMPlexTransform_Extrude *ex, PetscInt Nl)
135: {
136: const DMPolytopeType ct = DM_POLYTOPE_POINT;
137: const PetscInt Np = ex->periodic ? Nl : Nl + 1;
138: PetscInt Nc, No;
140: PetscFunctionBegin;
141: ex->Nt[ct] = 2;
142: Nc = 6 * Nl;
143: No = 2 * Nl;
144: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
145: ex->target[ct][0] = DM_POLYTOPE_POINT;
146: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
147: ex->size[ct][0] = Np;
148: ex->size[ct][1] = Nl;
149: /* cones for segments/tensor segments */
150: for (PetscInt i = 0; i < Nl; ++i) {
151: ex->cone[ct][6 * i + 0] = DM_POLYTOPE_POINT;
152: ex->cone[ct][6 * i + 1] = 0;
153: ex->cone[ct][6 * i + 2] = i;
154: ex->cone[ct][6 * i + 3] = DM_POLYTOPE_POINT;
155: ex->cone[ct][6 * i + 4] = 0;
156: ex->cone[ct][6 * i + 5] = (i + 1) % Np;
157: }
158: for (PetscInt i = 0; i < No; ++i) ex->ornt[ct][i] = 0;
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: /* DM_POLYTOPE_SEGMENT produces
163: Nl+1 segments, or Nl segments periodically and
164: Nl quads, or tensor quads
165: */
166: static PetscErrorCode DMPlexTransformExtrudeSetUp_Segment(DMPlexTransform_Extrude *ex, PetscInt Nl)
167: {
168: const DMPolytopeType ct = DM_POLYTOPE_SEGMENT;
169: const PetscInt Np = ex->periodic ? Nl : Nl + 1;
170: PetscInt Nc, No, coff, ooff;
172: PetscFunctionBegin;
173: ex->Nt[ct] = 2;
174: Nc = 8 * Np + 14 * Nl;
175: No = 2 * Np + 4 * Nl;
176: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
177: ex->target[ct][0] = DM_POLYTOPE_SEGMENT;
178: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
179: ex->size[ct][0] = Np;
180: ex->size[ct][1] = Nl;
181: /* cones for segments */
182: for (PetscInt i = 0; i < Np; ++i) {
183: ex->cone[ct][8 * i + 0] = DM_POLYTOPE_POINT;
184: ex->cone[ct][8 * i + 1] = 1;
185: ex->cone[ct][8 * i + 2] = 0;
186: ex->cone[ct][8 * i + 3] = i;
187: ex->cone[ct][8 * i + 4] = DM_POLYTOPE_POINT;
188: ex->cone[ct][8 * i + 5] = 1;
189: ex->cone[ct][8 * i + 6] = 1;
190: ex->cone[ct][8 * i + 7] = i;
191: }
192: for (PetscInt i = 0; i < 2 * Np; ++i) ex->ornt[ct][i] = 0;
193: /* cones for quads/tensor quads */
194: coff = 8 * Np;
195: ooff = 2 * Np;
196: for (PetscInt i = 0; i < Nl; ++i) {
197: if (ex->useTensor) {
198: ex->cone[ct][coff + 14 * i + 0] = DM_POLYTOPE_SEGMENT;
199: ex->cone[ct][coff + 14 * i + 1] = 0;
200: ex->cone[ct][coff + 14 * i + 2] = i;
201: ex->cone[ct][coff + 14 * i + 3] = DM_POLYTOPE_SEGMENT;
202: ex->cone[ct][coff + 14 * i + 4] = 0;
203: ex->cone[ct][coff + 14 * i + 5] = (i + 1) % Np;
204: ex->cone[ct][coff + 14 * i + 6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
205: ex->cone[ct][coff + 14 * i + 7] = 1;
206: ex->cone[ct][coff + 14 * i + 8] = 0;
207: ex->cone[ct][coff + 14 * i + 9] = i;
208: ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
209: ex->cone[ct][coff + 14 * i + 11] = 1;
210: ex->cone[ct][coff + 14 * i + 12] = 1;
211: ex->cone[ct][coff + 14 * i + 13] = i;
212: ex->ornt[ct][ooff + 4 * i + 0] = 0;
213: ex->ornt[ct][ooff + 4 * i + 1] = 0;
214: ex->ornt[ct][ooff + 4 * i + 2] = 0;
215: ex->ornt[ct][ooff + 4 * i + 3] = 0;
216: } else {
217: ex->cone[ct][coff + 14 * i + 0] = DM_POLYTOPE_SEGMENT;
218: ex->cone[ct][coff + 14 * i + 1] = 0;
219: ex->cone[ct][coff + 14 * i + 2] = i;
220: ex->cone[ct][coff + 14 * i + 3] = DM_POLYTOPE_SEGMENT;
221: ex->cone[ct][coff + 14 * i + 4] = 1;
222: ex->cone[ct][coff + 14 * i + 5] = 1;
223: ex->cone[ct][coff + 14 * i + 6] = i;
224: ex->cone[ct][coff + 14 * i + 7] = DM_POLYTOPE_SEGMENT;
225: ex->cone[ct][coff + 14 * i + 8] = 0;
226: ex->cone[ct][coff + 14 * i + 9] = (i + 1) % Np;
227: ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_SEGMENT;
228: ex->cone[ct][coff + 14 * i + 11] = 1;
229: ex->cone[ct][coff + 14 * i + 12] = 0;
230: ex->cone[ct][coff + 14 * i + 13] = i;
231: ex->ornt[ct][ooff + 4 * i + 0] = 0;
232: ex->ornt[ct][ooff + 4 * i + 1] = 0;
233: ex->ornt[ct][ooff + 4 * i + 2] = -1;
234: ex->ornt[ct][ooff + 4 * i + 3] = -1;
235: }
236: }
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /* DM_POLYTOPE_TRIANGLE produces
241: Nl+1 triangles, or Nl triangles periodically and
242: Nl triangular prisms/tensor triangular prisms
243: */
244: static PetscErrorCode DMPlexTransformExtrudeSetUp_Triangle(DMPlexTransform_Extrude *ex, PetscInt Nl)
245: {
246: const DMPolytopeType ct = DM_POLYTOPE_TRIANGLE;
247: const PetscInt Np = ex->periodic ? Nl : Nl + 1;
248: PetscInt Nc, No, coff, ooff;
250: PetscFunctionBegin;
251: ex->Nt[ct] = 2;
252: Nc = 12 * Np + 18 * Nl;
253: No = 3 * Np + 5 * Nl;
254: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
255: ex->target[ct][0] = DM_POLYTOPE_TRIANGLE;
256: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM;
257: ex->size[ct][0] = Np;
258: ex->size[ct][1] = Nl;
259: /* cones for triangles */
260: for (PetscInt i = 0; i < Np; ++i) {
261: ex->cone[ct][12 * i + 0] = DM_POLYTOPE_SEGMENT;
262: ex->cone[ct][12 * i + 1] = 1;
263: ex->cone[ct][12 * i + 2] = 0;
264: ex->cone[ct][12 * i + 3] = i;
265: ex->cone[ct][12 * i + 4] = DM_POLYTOPE_SEGMENT;
266: ex->cone[ct][12 * i + 5] = 1;
267: ex->cone[ct][12 * i + 6] = 1;
268: ex->cone[ct][12 * i + 7] = i;
269: ex->cone[ct][12 * i + 8] = DM_POLYTOPE_SEGMENT;
270: ex->cone[ct][12 * i + 9] = 1;
271: ex->cone[ct][12 * i + 10] = 2;
272: ex->cone[ct][12 * i + 11] = i;
273: }
274: for (PetscInt i = 0; i < 3 * Np; ++i) ex->ornt[ct][i] = 0;
275: /* cones for triangular prisms/tensor triangular prisms */
276: coff = 12 * Np;
277: ooff = 3 * Np;
278: for (PetscInt i = 0; i < Nl; ++i) {
279: if (ex->useTensor) {
280: ex->cone[ct][coff + 18 * i + 0] = DM_POLYTOPE_TRIANGLE;
281: ex->cone[ct][coff + 18 * i + 1] = 0;
282: ex->cone[ct][coff + 18 * i + 2] = i;
283: ex->cone[ct][coff + 18 * i + 3] = DM_POLYTOPE_TRIANGLE;
284: ex->cone[ct][coff + 18 * i + 4] = 0;
285: ex->cone[ct][coff + 18 * i + 5] = (i + 1) % Np;
286: ex->cone[ct][coff + 18 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
287: ex->cone[ct][coff + 18 * i + 7] = 1;
288: ex->cone[ct][coff + 18 * i + 8] = 0;
289: ex->cone[ct][coff + 18 * i + 9] = i;
290: ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
291: ex->cone[ct][coff + 18 * i + 11] = 1;
292: ex->cone[ct][coff + 18 * i + 12] = 1;
293: ex->cone[ct][coff + 18 * i + 13] = i;
294: ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
295: ex->cone[ct][coff + 18 * i + 15] = 1;
296: ex->cone[ct][coff + 18 * i + 16] = 2;
297: ex->cone[ct][coff + 18 * i + 17] = i;
298: ex->ornt[ct][ooff + 5 * i + 0] = 0;
299: ex->ornt[ct][ooff + 5 * i + 1] = 0;
300: ex->ornt[ct][ooff + 5 * i + 2] = 0;
301: ex->ornt[ct][ooff + 5 * i + 3] = 0;
302: ex->ornt[ct][ooff + 5 * i + 4] = 0;
303: } else {
304: ex->cone[ct][coff + 18 * i + 0] = DM_POLYTOPE_TRIANGLE;
305: ex->cone[ct][coff + 18 * i + 1] = 0;
306: ex->cone[ct][coff + 18 * i + 2] = i;
307: ex->cone[ct][coff + 18 * i + 3] = DM_POLYTOPE_TRIANGLE;
308: ex->cone[ct][coff + 18 * i + 4] = 0;
309: ex->cone[ct][coff + 18 * i + 5] = (i + 1) % Np;
310: ex->cone[ct][coff + 18 * i + 6] = DM_POLYTOPE_QUADRILATERAL;
311: ex->cone[ct][coff + 18 * i + 7] = 1;
312: ex->cone[ct][coff + 18 * i + 8] = 0;
313: ex->cone[ct][coff + 18 * i + 9] = i;
314: ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
315: ex->cone[ct][coff + 18 * i + 11] = 1;
316: ex->cone[ct][coff + 18 * i + 12] = 1;
317: ex->cone[ct][coff + 18 * i + 13] = i;
318: ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
319: ex->cone[ct][coff + 18 * i + 15] = 1;
320: ex->cone[ct][coff + 18 * i + 16] = 2;
321: ex->cone[ct][coff + 18 * i + 17] = i;
322: ex->ornt[ct][ooff + 5 * i + 0] = -2;
323: ex->ornt[ct][ooff + 5 * i + 1] = 0;
324: ex->ornt[ct][ooff + 5 * i + 2] = 0;
325: ex->ornt[ct][ooff + 5 * i + 3] = 0;
326: ex->ornt[ct][ooff + 5 * i + 4] = 0;
327: }
328: }
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: /* DM_POLYTOPE_QUADRILATERAL produces
333: Nl+1 quads, or Nl quads periodically and
334: Nl hexes/tensor hexes
335: */
336: static PetscErrorCode DMPlexTransformExtrudeSetUp_Quadrilateral(DMPlexTransform_Extrude *ex, PetscInt Nl)
337: {
338: const DMPolytopeType ct = DM_POLYTOPE_QUADRILATERAL;
339: const PetscInt Np = ex->periodic ? Nl : Nl + 1;
340: PetscInt Nc, No, coff, ooff;
342: PetscFunctionBegin;
343: ex->Nt[ct] = 2;
344: Nc = 16 * Np + 22 * Nl;
345: No = 4 * Np + 6 * Nl;
346: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
347: ex->target[ct][0] = DM_POLYTOPE_QUADRILATERAL;
348: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON;
349: ex->size[ct][0] = Np;
350: ex->size[ct][1] = Nl;
351: /* cones for quads */
352: for (PetscInt i = 0; i < Np; ++i) {
353: ex->cone[ct][16 * i + 0] = DM_POLYTOPE_SEGMENT;
354: ex->cone[ct][16 * i + 1] = 1;
355: ex->cone[ct][16 * i + 2] = 0;
356: ex->cone[ct][16 * i + 3] = i;
357: ex->cone[ct][16 * i + 4] = DM_POLYTOPE_SEGMENT;
358: ex->cone[ct][16 * i + 5] = 1;
359: ex->cone[ct][16 * i + 6] = 1;
360: ex->cone[ct][16 * i + 7] = i;
361: ex->cone[ct][16 * i + 8] = DM_POLYTOPE_SEGMENT;
362: ex->cone[ct][16 * i + 9] = 1;
363: ex->cone[ct][16 * i + 10] = 2;
364: ex->cone[ct][16 * i + 11] = i;
365: ex->cone[ct][16 * i + 12] = DM_POLYTOPE_SEGMENT;
366: ex->cone[ct][16 * i + 13] = 1;
367: ex->cone[ct][16 * i + 14] = 3;
368: ex->cone[ct][16 * i + 15] = i;
369: }
370: for (PetscInt i = 0; i < 4 * Np; ++i) ex->ornt[ct][i] = 0;
371: /* cones for hexes/tensor hexes */
372: coff = 16 * Np;
373: ooff = 4 * Np;
374: for (PetscInt i = 0; i < Nl; ++i) {
375: if (ex->useTensor) {
376: ex->cone[ct][coff + 22 * i + 0] = DM_POLYTOPE_QUADRILATERAL;
377: ex->cone[ct][coff + 22 * i + 1] = 0;
378: ex->cone[ct][coff + 22 * i + 2] = i;
379: ex->cone[ct][coff + 22 * i + 3] = DM_POLYTOPE_QUADRILATERAL;
380: ex->cone[ct][coff + 22 * i + 4] = 0;
381: ex->cone[ct][coff + 22 * i + 5] = (i + 1) % Np;
382: ex->cone[ct][coff + 22 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
383: ex->cone[ct][coff + 22 * i + 7] = 1;
384: ex->cone[ct][coff + 22 * i + 8] = 0;
385: ex->cone[ct][coff + 22 * i + 9] = i;
386: ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
387: ex->cone[ct][coff + 22 * i + 11] = 1;
388: ex->cone[ct][coff + 22 * i + 12] = 1;
389: ex->cone[ct][coff + 22 * i + 13] = i;
390: ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
391: ex->cone[ct][coff + 22 * i + 15] = 1;
392: ex->cone[ct][coff + 22 * i + 16] = 2;
393: ex->cone[ct][coff + 22 * i + 17] = i;
394: ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
395: ex->cone[ct][coff + 22 * i + 19] = 1;
396: ex->cone[ct][coff + 22 * i + 20] = 3;
397: ex->cone[ct][coff + 22 * i + 21] = i;
398: ex->ornt[ct][ooff + 6 * i + 0] = 0;
399: ex->ornt[ct][ooff + 6 * i + 1] = 0;
400: ex->ornt[ct][ooff + 6 * i + 2] = 0;
401: ex->ornt[ct][ooff + 6 * i + 3] = 0;
402: ex->ornt[ct][ooff + 6 * i + 4] = 0;
403: ex->ornt[ct][ooff + 6 * i + 5] = 0;
404: } else {
405: ex->cone[ct][coff + 22 * i + 0] = DM_POLYTOPE_QUADRILATERAL;
406: ex->cone[ct][coff + 22 * i + 1] = 0;
407: ex->cone[ct][coff + 22 * i + 2] = i;
408: ex->cone[ct][coff + 22 * i + 3] = DM_POLYTOPE_QUADRILATERAL;
409: ex->cone[ct][coff + 22 * i + 4] = 0;
410: ex->cone[ct][coff + 22 * i + 5] = (i + 1) % Np;
411: ex->cone[ct][coff + 22 * i + 6] = DM_POLYTOPE_QUADRILATERAL;
412: ex->cone[ct][coff + 22 * i + 7] = 1;
413: ex->cone[ct][coff + 22 * i + 8] = 0;
414: ex->cone[ct][coff + 22 * i + 9] = i;
415: ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
416: ex->cone[ct][coff + 22 * i + 11] = 1;
417: ex->cone[ct][coff + 22 * i + 12] = 2;
418: ex->cone[ct][coff + 22 * i + 13] = i;
419: ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
420: ex->cone[ct][coff + 22 * i + 15] = 1;
421: ex->cone[ct][coff + 22 * i + 16] = 1;
422: ex->cone[ct][coff + 22 * i + 17] = i;
423: ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_QUADRILATERAL;
424: ex->cone[ct][coff + 22 * i + 19] = 1;
425: ex->cone[ct][coff + 22 * i + 20] = 3;
426: ex->cone[ct][coff + 22 * i + 21] = i;
427: ex->ornt[ct][ooff + 6 * i + 0] = -2;
428: ex->ornt[ct][ooff + 6 * i + 1] = 0;
429: ex->ornt[ct][ooff + 6 * i + 2] = 0;
430: ex->ornt[ct][ooff + 6 * i + 3] = 0;
431: ex->ornt[ct][ooff + 6 * i + 4] = 0;
432: ex->ornt[ct][ooff + 6 * i + 5] = 1;
433: }
434: }
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*
439: The refine types for extrusion are:
441: ct: For any normally extruded point
442: ct + 100: For any point which should just return itself
443: */
444: static PetscErrorCode DMPlexTransformSetUp_Extrude(DMPlexTransform tr)
445: {
446: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
447: DM dm;
448: DMLabel active;
449: PetscInt Nl = ex->layers, l, ict;
451: PetscFunctionBegin;
452: PetscCall(DMPlexTransformExtrudeComputeExtrusionDim(tr));
453: PetscCall(DMPlexTransformGetDM(tr, &dm));
454: PetscCall(DMPlexTransformGetActive(tr, &active));
455: if (active) {
456: DMLabel celltype;
457: PetscInt pStart, pEnd, p;
459: PetscCall(DMPlexGetCellTypeLabel(dm, &celltype));
460: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
461: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
462: for (p = pStart; p < pEnd; ++p) {
463: PetscInt ct, val;
465: PetscCall(DMLabelGetValue(celltype, p, &ct));
466: PetscCall(DMLabelGetValue(active, p, &val));
467: if (val < 0) {
468: PetscCall(DMLabelSetValue(tr->trType, p, ct + 100));
469: } else {
470: PetscCall(DMLabelSetValue(tr->trType, p, ct));
471: }
472: }
473: }
474: 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));
475: for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
476: ex->Nt[ict] = -1;
477: ex->target[ict] = NULL;
478: ex->size[ict] = NULL;
479: ex->cone[ict] = NULL;
480: ex->ornt[ict] = NULL;
481: }
482: PetscCall(DMPlexTransformExtrudeSetUp_Point(ex, Nl));
483: PetscCall(DMPlexTransformExtrudeSetUp_Segment(ex, Nl));
484: PetscCall(DMPlexTransformExtrudeSetUp_Triangle(ex, Nl));
485: PetscCall(DMPlexTransformExtrudeSetUp_Quadrilateral(ex, Nl));
486: /* Layers positions */
487: if (!ex->Nth) {
488: if (ex->symmetric)
489: for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers - ex->thickness / 2;
490: else
491: for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers;
492: } else {
493: ex->thickness = 0.;
494: ex->layerPos[0] = 0.;
495: for (l = 0; l < ex->layers; ++l) {
496: const PetscReal t = ex->thicknesses[PetscMin(l, ex->Nth - 1)];
497: ex->layerPos[l + 1] = ex->layerPos[l] + t;
498: ex->thickness += t;
499: }
500: if (ex->symmetric)
501: for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] -= ex->thickness / 2.;
502: }
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr)
507: {
508: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
509: PetscInt ct;
511: PetscFunctionBegin;
512: if (ex->target) {
513: for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) PetscCall(PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct]));
514: }
515: PetscCall(PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt));
516: PetscCall(PetscFree(ex->layerPos));
517: PetscCall(PetscFree(ex));
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
522: {
523: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
524: DMLabel trType = tr->trType;
525: PetscInt rt;
527: PetscFunctionBeginHot;
528: *rnew = r;
529: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
530: if (!so) PetscFunctionReturn(PETSC_SUCCESS);
531: if (trType) {
532: PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
533: if (rt >= 100) PetscFunctionReturn(PETSC_SUCCESS);
534: }
535: if (ex->useTensor) {
536: switch (sct) {
537: case DM_POLYTOPE_POINT:
538: break;
539: case DM_POLYTOPE_SEGMENT:
540: switch (tct) {
541: case DM_POLYTOPE_SEGMENT:
542: break;
543: case DM_POLYTOPE_SEG_PRISM_TENSOR:
544: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
545: break;
546: default:
547: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
548: }
549: break;
550: // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
551: case DM_POLYTOPE_TRIANGLE:
552: break;
553: case DM_POLYTOPE_QUADRILATERAL:
554: break;
555: default:
556: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
557: }
558: } else {
559: switch (sct) {
560: case DM_POLYTOPE_POINT:
561: break;
562: case DM_POLYTOPE_SEGMENT:
563: switch (tct) {
564: case DM_POLYTOPE_SEGMENT:
565: break;
566: case DM_POLYTOPE_QUADRILATERAL:
567: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0);
568: break;
569: default:
570: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
571: }
572: break;
573: // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
574: case DM_POLYTOPE_TRIANGLE:
575: break;
576: case DM_POLYTOPE_QUADRILATERAL:
577: break;
578: default:
579: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
580: }
581: }
582: PetscFunctionReturn(PETSC_SUCCESS);
583: }
585: static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
586: {
587: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
588: DMLabel trType = tr->trType;
589: PetscBool ignore = PETSC_FALSE, identity = PETSC_FALSE;
590: PetscInt val = 0;
592: PetscFunctionBegin;
593: if (trType) {
594: PetscCall(DMLabelGetValue(trType, p, &val));
595: identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE;
596: } else {
597: ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE;
598: }
599: if (rt) *rt = val;
600: if (ignore) {
601: /* Ignore cells that cannot be extruded */
602: *Nt = 0;
603: *target = NULL;
604: *size = NULL;
605: *cone = NULL;
606: *ornt = NULL;
607: } else if (identity) {
608: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
609: } else {
610: *Nt = ex->Nt[source];
611: *target = ex->target[source];
612: *size = ex->size[source];
613: *cone = ex->cone[source];
614: *ornt = ex->ornt[source];
615: }
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /* Computes new vertex along normal */
620: static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
621: {
622: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
623: DM dm;
624: DMLabel active;
625: PetscReal ones2[2] = {0., 1.}, ones3[3] = {0., 0., 1.};
626: PetscReal normal[3] = {0., 0., 0.}, norm;
627: PetscBool computeNormal;
628: PetscInt dim, dEx = ex->cdimEx, cStart, cEnd, d;
630: PetscFunctionBeginHot;
631: PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]);
632: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
633: PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv);
634: PetscCheck(dE == ex->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %" PetscInt_FMT " != %" PetscInt_FMT " original dimension", dE, ex->cdim);
635: PetscCheck(dEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", dEx);
637: PetscCall(DMPlexTransformGetDM(tr, &dm));
638: PetscCall(DMGetDimension(dm, &dim));
639: PetscCall(DMPlexTransformGetActive(tr, &active));
640: computeNormal = dim != ex->cdim && !ex->useNormal ? PETSC_TRUE : PETSC_FALSE;
641: if (computeNormal) {
642: PetscInt *closure = NULL;
643: PetscInt closureSize, cl;
645: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
646: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure));
647: for (cl = 0; cl < closureSize * 2; cl += 2) {
648: if ((closure[cl] >= cStart) && (closure[cl] < cEnd)) {
649: PetscReal cnormal[3] = {0, 0, 0};
651: PetscCall(DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal));
652: for (d = 0; d < dEx; ++d) normal[d] += cnormal[d];
653: }
654: }
655: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure));
656: } else if (ex->useNormal) {
657: for (d = 0; d < dEx; ++d) normal[d] = ex->normal[d];
658: } else if (active) { // find an active point in the closure of p and use its coordinate normal as the extrusion direction
659: PetscInt *closure = NULL;
660: PetscInt closureSize, cl, pStart, pEnd;
662: PetscCall(DMPlexGetDepthStratum(dm, ex->cdimEx - 1, &pStart, &pEnd));
663: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure));
664: for (cl = 0; cl < closureSize * 2; cl += 2) {
665: if ((closure[cl] >= pStart) && (closure[cl] < pEnd)) {
666: PetscReal cnormal[3] = {0, 0, 0};
667: const PetscInt *supp;
668: PetscInt suppSize;
670: PetscCall(DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal));
671: PetscCall(DMPlexGetSupportSize(dm, closure[cl], &suppSize));
672: PetscCall(DMPlexGetSupport(dm, closure[cl], &supp));
673: // Only use external faces, so I can get the orientation from any cell
674: if (suppSize == 1) {
675: const PetscInt *cone, *ornt;
676: PetscInt coneSize, c;
678: PetscCall(DMPlexGetConeSize(dm, supp[0], &coneSize));
679: PetscCall(DMPlexGetCone(dm, supp[0], &cone));
680: PetscCall(DMPlexGetConeOrientation(dm, supp[0], &ornt));
681: for (c = 0; c < coneSize; ++c)
682: if (cone[c] == closure[cl]) break;
683: PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Asymmetry in cone/support");
684: if (ornt[c] < 0)
685: for (d = 0; d < dEx; ++d) cnormal[d] *= -1.;
686: for (d = 0; d < dEx; ++d) normal[d] += cnormal[d];
687: }
688: }
689: }
690: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure));
691: } else if (ex->cdimEx == 2) {
692: for (d = 0; d < dEx; ++d) normal[d] = ones2[d];
693: } else if (ex->cdimEx == 3) {
694: for (d = 0; d < dEx; ++d) normal[d] = ones3[d];
695: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
696: if (ex->normalFunc) {
697: PetscScalar n[3];
698: PetscReal x[3], dot;
700: for (d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]);
701: PetscCall((*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL));
702: for (dot = 0, d = 0; d < dEx; d++) dot += PetscRealPart(n[d]) * normal[d];
703: for (d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]);
704: }
706: for (d = 0, norm = 0.0; d < dEx; ++d) norm += PetscSqr(normal[d]);
707: for (d = 0; d < dEx; ++d) normal[d] *= norm == 0.0 ? 1.0 : 1. / PetscSqrtReal(norm);
708: for (d = 0; d < dEx; ++d) out[d] = normal[d] * ex->layerPos[r];
709: for (d = 0; d < ex->cdim; ++d) out[d] += in[d];
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr)
714: {
715: PetscFunctionBegin;
716: tr->ops->view = DMPlexTransformView_Extrude;
717: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude;
718: tr->ops->setup = DMPlexTransformSetUp_Extrude;
719: tr->ops->destroy = DMPlexTransformDestroy_Extrude;
720: tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude;
721: tr->ops->celltransform = DMPlexTransformCellTransform_Extrude;
722: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude;
723: tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude;
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr)
728: {
729: DMPlexTransform_Extrude *ex;
730: DM dm;
731: PetscInt dim;
733: PetscFunctionBegin;
735: PetscCall(PetscNew(&ex));
736: tr->data = ex;
737: ex->thickness = 1.;
738: ex->useTensor = PETSC_TRUE;
739: ex->symmetric = PETSC_FALSE;
740: ex->periodic = PETSC_FALSE;
741: ex->useNormal = PETSC_FALSE;
742: ex->layerPos = NULL;
743: PetscCall(DMPlexTransformGetDM(tr, &dm));
744: PetscCall(DMGetDimension(dm, &dim));
745: PetscCall(DMGetCoordinateDim(dm, &ex->cdim));
746: PetscCall(DMPlexTransformExtrudeSetLayers(tr, 1));
747: PetscCall(DMPlexTransformInitialize_Extrude(tr));
748: PetscFunctionReturn(PETSC_SUCCESS);
749: }
751: /*@
752: DMPlexTransformExtrudeGetLayers - Get the number of extruded layers.
754: Not Collective
756: Input Parameter:
757: . tr - The `DMPlexTransform`
759: Output Parameter:
760: . layers - The number of layers
762: Level: intermediate
764: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetLayers()`
765: @*/
766: PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers)
767: {
768: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
770: PetscFunctionBegin;
772: PetscAssertPointer(layers, 2);
773: *layers = ex->layers;
774: PetscFunctionReturn(PETSC_SUCCESS);
775: }
777: /*@
778: DMPlexTransformExtrudeSetLayers - Set the number of extruded layers.
780: Not Collective
782: Input Parameters:
783: + tr - The `DMPlexTransform`
784: - layers - The number of layers
786: Level: intermediate
788: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetLayers()`
789: @*/
790: PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers)
791: {
792: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
794: PetscFunctionBegin;
796: ex->layers = layers;
797: PetscCall(PetscFree(ex->layerPos));
798: PetscCall(PetscCalloc1(ex->layers + 1, &ex->layerPos));
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: /*@
803: DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers
805: Not Collective
807: Input Parameter:
808: . tr - The `DMPlexTransform`
810: Output Parameter:
811: . thickness - The total thickness of the layers
813: Level: intermediate
815: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`
816: @*/
817: PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness)
818: {
819: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
821: PetscFunctionBegin;
823: PetscAssertPointer(thickness, 2);
824: *thickness = ex->thickness;
825: PetscFunctionReturn(PETSC_SUCCESS);
826: }
828: /*@
829: DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers
831: Not Collective
833: Input Parameters:
834: + tr - The `DMPlexTransform`
835: - thickness - The total thickness of the layers
837: Level: intermediate
839: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetThickness()`
840: @*/
841: PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness)
842: {
843: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
845: PetscFunctionBegin;
847: PetscCheck(thickness > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double)thickness);
848: ex->thickness = thickness;
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: /*@
853: DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells
855: Not Collective
857: Input Parameter:
858: . tr - The `DMPlexTransform`
860: Output Parameter:
861: . useTensor - The flag to use tensor cells
863: Note:
864: This flag determines the orientation behavior of the created points.
866: For example, if tensor is `PETSC_TRUE`, then
867: .vb
868: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
869: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
870: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
871: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
872: .ve
874: Level: intermediate
876: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetTensor()`
877: @*/
878: PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
879: {
880: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
882: PetscFunctionBegin;
884: PetscAssertPointer(useTensor, 2);
885: *useTensor = ex->useTensor;
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
889: /*@
890: DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells
892: Not Collective
894: Input Parameters:
895: + tr - The `DMPlexTransform`
896: - useTensor - The flag for tensor cells
898: Note:
899: This flag determines the orientation behavior of the created points
900: For example, if tensor is `PETSC_TRUE`, then
901: .vb
902: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
903: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
904: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
905: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
906: .ve
908: Level: intermediate
910: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetTensor()`
911: @*/
912: PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
913: {
914: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
916: PetscFunctionBegin;
918: ex->useTensor = useTensor;
919: PetscFunctionReturn(PETSC_SUCCESS);
920: }
922: /*@
923: DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface
925: Not Collective
927: Input Parameter:
928: . tr - The `DMPlexTransform`
930: Output Parameter:
931: . symmetric - The flag to extrude symmetrically
933: Level: intermediate
935: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetSymmetric()`, `DMPlexTransformExtrudeGetPeriodic()`
936: @*/
937: PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric)
938: {
939: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
941: PetscFunctionBegin;
943: PetscAssertPointer(symmetric, 2);
944: *symmetric = ex->symmetric;
945: PetscFunctionReturn(PETSC_SUCCESS);
946: }
948: /*@
949: DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface
951: Not Collective
953: Input Parameters:
954: + tr - The `DMPlexTransform`
955: - symmetric - The flag to extrude symmetrically
957: Level: intermediate
959: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetSymmetric()`, `DMPlexTransformExtrudeSetPeriodic()`
960: @*/
961: PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric)
962: {
963: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
965: PetscFunctionBegin;
967: ex->symmetric = symmetric;
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@
972: DMPlexTransformExtrudeGetPeriodic - Get the flag to extrude periodically from the initial surface
974: Not Collective
976: Input Parameter:
977: . tr - The `DMPlexTransform`
979: Output Parameter:
980: . periodic - The flag to extrude periodically
982: Level: intermediate
984: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetPeriodic()`, `DMPlexTransformExtrudeGetSymmetric()`
985: @*/
986: PetscErrorCode DMPlexTransformExtrudeGetPeriodic(DMPlexTransform tr, PetscBool *periodic)
987: {
988: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
990: PetscFunctionBegin;
992: PetscAssertPointer(periodic, 2);
993: *periodic = ex->periodic;
994: PetscFunctionReturn(PETSC_SUCCESS);
995: }
997: /*@
998: DMPlexTransformExtrudeSetPeriodic - Set the flag to extrude periodically from the initial surface
1000: Not Collective
1002: Input Parameters:
1003: + tr - The `DMPlexTransform`
1004: - periodic - The flag to extrude periodically
1006: Level: intermediate
1008: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetPeriodic()`, `DMPlexTransformExtrudeSetSymmetric()`
1009: @*/
1010: PetscErrorCode DMPlexTransformExtrudeSetPeriodic(DMPlexTransform tr, PetscBool periodic)
1011: {
1012: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1014: PetscFunctionBegin;
1016: ex->periodic = periodic;
1017: PetscFunctionReturn(PETSC_SUCCESS);
1018: }
1020: /*@
1021: DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector
1023: Not Collective
1025: Input Parameter:
1026: . tr - The `DMPlexTransform`
1028: Output Parameter:
1029: . normal - The extrusion direction
1031: Note:
1032: The user passes in an array, which is filled by the library.
1034: Level: intermediate
1036: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetNormal()`
1037: @*/
1038: PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[])
1039: {
1040: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1041: PetscInt d;
1043: PetscFunctionBegin;
1045: if (ex->useNormal) {
1046: for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d];
1047: } else {
1048: for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.;
1049: }
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: /*@
1054: DMPlexTransformExtrudeSetNormal - Set the extrusion normal
1056: Not Collective
1058: Input Parameters:
1059: + tr - The `DMPlexTransform`
1060: - normal - The extrusion direction
1062: Level: intermediate
1064: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`
1065: @*/
1066: PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[])
1067: {
1068: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1069: PetscInt d;
1071: PetscFunctionBegin;
1073: ex->useNormal = PETSC_TRUE;
1074: for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d];
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: /*@C
1079: DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal
1081: Not Collective
1083: Input Parameters:
1084: + tr - The `DMPlexTransform`
1085: - normalFunc - A function determining the extrusion direction, see `PetscSimplePointFn` for the calling sequence
1087: Level: intermediate
1089: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`, `PetscSimplePointFn`
1090: @*/
1091: PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFn *normalFunc)
1092: {
1093: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1095: PetscFunctionBegin;
1097: ex->normalFunc = normalFunc;
1098: PetscFunctionReturn(PETSC_SUCCESS);
1099: }
1101: /*@
1102: DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer
1104: Not Collective
1106: Input Parameters:
1107: + tr - The `DMPlexTransform`
1108: . Nth - The number of thicknesses
1109: - thicknesses - The array of thicknesses
1111: Level: intermediate
1113: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeGetThickness()`
1114: @*/
1115: PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[])
1116: {
1117: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1118: PetscInt t;
1120: PetscFunctionBegin;
1122: PetscCheck(Nth > 0, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %" PetscInt_FMT " must be positive", Nth);
1123: ex->Nth = PetscMin(Nth, ex->layers);
1124: PetscCall(PetscFree(ex->thicknesses));
1125: PetscCall(PetscMalloc1(ex->Nth, &ex->thicknesses));
1126: for (t = 0; t < ex->Nth; ++t) {
1127: PetscCheck(thicknesses[t] > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %" PetscInt_FMT " must be positive", (double)thicknesses[t], t);
1128: ex->thicknesses[t] = thicknesses[t];
1129: }
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }