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: } else {
20: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
21: }
22: PetscFunctionReturn(PETSC_SUCCESS);
23: }
25: static PetscErrorCode DMPlexTransformSetFromOptions_Extrude(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
26: {
27: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
28: PetscReal th, normal[3], *thicknesses;
29: PetscInt nl, Nc;
30: PetscBool tensor, sym, flg;
31: char funcname[PETSC_MAX_PATH_LEN];
33: PetscFunctionBegin;
34: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Extrusion Options");
35: PetscCall(PetscOptionsBoundedInt("-dm_plex_transform_extrude_layers", "Number of layers to extrude", "", ex->layers, &nl, &flg, 1));
36: if (flg) PetscCall(DMPlexTransformExtrudeSetLayers(tr, nl));
37: PetscCall(PetscOptionsReal("-dm_plex_transform_extrude_thickness", "Total thickness of extruded layers", "", ex->thickness, &th, &flg));
38: if (flg) PetscCall(DMPlexTransformExtrudeSetThickness(tr, th));
39: PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg));
40: if (flg) PetscCall(DMPlexTransformExtrudeSetTensor(tr, tensor));
41: PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_symmetric", "Extrude layers symmetrically about the surface", "", ex->symmetric, &sym, &flg));
42: if (flg) PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, sym));
43: Nc = 3;
44: PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_normal", "Input normal vector for extrusion", "DMPlexTransformExtrudeSetNormal", normal, &Nc, &flg));
45: if (flg) {
46: // Extrusion dimension might not yet be determined
47: 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);
48: PetscCall(DMPlexTransformExtrudeSetNormal(tr, normal));
49: }
50: PetscCall(PetscOptionsString("-dm_plex_transform_extrude_normal_function", "Function to determine normal vector", "DMPlexTransformExtrudeSetNormalFunction", funcname, funcname, sizeof(funcname), &flg));
51: if (flg) {
52: PetscSimplePointFunc normalFunc;
54: PetscCall(PetscDLSym(NULL, funcname, (void **)&normalFunc));
55: PetscCall(DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc));
56: }
57: nl = ex->layers;
58: PetscCall(PetscMalloc1(nl, &thicknesses));
59: PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_thicknesses", "Thickness of each individual extruded layer", "", thicknesses, &nl, &flg));
60: if (flg) {
61: PetscCheck(nl, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one thickness for -dm_plex_transform_extrude_thicknesses");
62: PetscCall(DMPlexTransformExtrudeSetThicknesses(tr, nl, thicknesses));
63: }
64: PetscCall(PetscFree(thicknesses));
65: PetscOptionsHeadEnd();
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: /* Determine the implicit dimension pre-extrusion (either the implicit dimension of the DM or of a point in the active set for the transform).
70: If that dimension is the same as the current coordinate dimension (ex->dim), the extruded mesh will have a coordinate dimension one greater;
71: Otherwise the coordinate dimension will be kept. */
72: static PetscErrorCode DMPlexTransformExtrudeComputeExtrusionDim(DMPlexTransform tr)
73: {
74: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
75: DM dm;
76: DMLabel active;
77: PetscInt dim, dimExtPoint, dimExtPointG;
79: PetscFunctionBegin;
80: PetscCall(DMPlexTransformGetDM(tr, &dm));
81: PetscCall(DMGetDimension(dm, &dim));
82: PetscCall(DMPlexTransformGetActive(tr, &active));
83: if (active) {
84: IS valueIS, pointIS;
85: const PetscInt *values, *points;
86: DMPolytopeType ct;
87: PetscInt Nv, Np;
89: dimExtPoint = 0;
90: PetscCall(DMLabelGetValueIS(active, &valueIS));
91: PetscCall(ISGetLocalSize(valueIS, &Nv));
92: PetscCall(ISGetIndices(valueIS, &values));
93: for (PetscInt v = 0; v < Nv; ++v) {
94: PetscCall(DMLabelGetStratumIS(active, values[v], &pointIS));
95: PetscCall(ISGetLocalSize(pointIS, &Np));
96: PetscCall(ISGetIndices(pointIS, &points));
97: for (PetscInt p = 0; p < Np; ++p) {
98: PetscCall(DMPlexGetCellType(dm, points[p], &ct));
99: dimExtPoint = PetscMax(dimExtPoint, DMPolytopeTypeGetDim(ct));
100: }
101: PetscCall(ISRestoreIndices(pointIS, &points));
102: PetscCall(ISDestroy(&pointIS));
103: }
104: PetscCall(ISRestoreIndices(valueIS, &values));
105: PetscCall(ISDestroy(&valueIS));
106: } else dimExtPoint = dim;
107: PetscCall(MPIU_Allreduce(&dimExtPoint, &dimExtPointG, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)tr)));
108: ex->dimEx = PetscMax(dim, dimExtPointG + 1);
109: ex->cdimEx = ex->cdim == dimExtPointG ? ex->cdim + 1 : ex->cdim;
110: PetscCheck(ex->dimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Topological dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->dimEx);
111: PetscCheck(ex->cdimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->cdimEx);
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm)
116: {
117: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
118: PetscInt dim;
120: PetscFunctionBegin;
121: PetscCall(DMGetDimension(dm, &dim));
122: PetscCall(DMSetDimension(tdm, ex->dimEx));
123: PetscCall(DMSetCoordinateDim(tdm, ex->cdimEx));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: /*
128: The refine types for extrusion are:
130: ct: For any normally extruded point
131: ct + 100: For any point which should just return itself
132: */
133: static PetscErrorCode DMPlexTransformSetUp_Extrude(DMPlexTransform tr)
134: {
135: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
136: DM dm;
137: DMLabel active;
138: DMPolytopeType ct;
139: PetscInt Nl = ex->layers, l, i, ict, Nc, No, coff, ooff;
141: PetscFunctionBegin;
142: PetscCall(DMPlexTransformExtrudeComputeExtrusionDim(tr));
143: PetscCall(DMPlexTransformGetDM(tr, &dm));
144: PetscCall(DMPlexTransformGetActive(tr, &active));
145: if (active) {
146: DMLabel celltype;
147: PetscInt pStart, pEnd, p;
149: PetscCall(DMPlexGetCellTypeLabel(dm, &celltype));
150: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
151: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
152: for (p = pStart; p < pEnd; ++p) {
153: PetscInt ct, val;
155: PetscCall(DMLabelGetValue(celltype, p, &ct));
156: PetscCall(DMLabelGetValue(active, p, &val));
157: if (val < 0) {
158: PetscCall(DMLabelSetValue(tr->trType, p, ct + 100));
159: } else {
160: PetscCall(DMLabelSetValue(tr->trType, p, ct));
161: }
162: }
163: }
164: 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));
165: for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
166: ex->Nt[ict] = -1;
167: ex->target[ict] = NULL;
168: ex->size[ict] = NULL;
169: ex->cone[ict] = NULL;
170: ex->ornt[ict] = NULL;
171: }
172: /* DM_POLYTOPE_POINT produces Nl+1 points and Nl segments/tensor segments */
173: ct = DM_POLYTOPE_POINT;
174: ex->Nt[ct] = 2;
175: Nc = 6 * Nl;
176: No = 2 * Nl;
177: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
178: ex->target[ct][0] = DM_POLYTOPE_POINT;
179: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
180: ex->size[ct][0] = Nl + 1;
181: ex->size[ct][1] = Nl;
182: /* cones for segments/tensor segments */
183: for (i = 0; i < Nl; ++i) {
184: ex->cone[ct][6 * i + 0] = DM_POLYTOPE_POINT;
185: ex->cone[ct][6 * i + 1] = 0;
186: ex->cone[ct][6 * i + 2] = i;
187: ex->cone[ct][6 * i + 3] = DM_POLYTOPE_POINT;
188: ex->cone[ct][6 * i + 4] = 0;
189: ex->cone[ct][6 * i + 5] = i + 1;
190: }
191: for (i = 0; i < No; ++i) ex->ornt[ct][i] = 0;
192: /* DM_POLYTOPE_SEGMENT produces Nl+1 segments and Nl quads/tensor quads */
193: ct = DM_POLYTOPE_SEGMENT;
194: ex->Nt[ct] = 2;
195: Nc = 8 * (Nl + 1) + 14 * Nl;
196: No = 2 * (Nl + 1) + 4 * Nl;
197: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
198: ex->target[ct][0] = DM_POLYTOPE_SEGMENT;
199: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
200: ex->size[ct][0] = Nl + 1;
201: ex->size[ct][1] = Nl;
202: /* cones for segments */
203: for (i = 0; i < Nl + 1; ++i) {
204: ex->cone[ct][8 * i + 0] = DM_POLYTOPE_POINT;
205: ex->cone[ct][8 * i + 1] = 1;
206: ex->cone[ct][8 * i + 2] = 0;
207: ex->cone[ct][8 * i + 3] = i;
208: ex->cone[ct][8 * i + 4] = DM_POLYTOPE_POINT;
209: ex->cone[ct][8 * i + 5] = 1;
210: ex->cone[ct][8 * i + 6] = 1;
211: ex->cone[ct][8 * i + 7] = i;
212: }
213: for (i = 0; i < 2 * (Nl + 1); ++i) ex->ornt[ct][i] = 0;
214: /* cones for quads/tensor quads */
215: coff = 8 * (Nl + 1);
216: ooff = 2 * (Nl + 1);
217: for (i = 0; i < Nl; ++i) {
218: if (ex->useTensor) {
219: ex->cone[ct][coff + 14 * i + 0] = DM_POLYTOPE_SEGMENT;
220: ex->cone[ct][coff + 14 * i + 1] = 0;
221: ex->cone[ct][coff + 14 * i + 2] = i;
222: ex->cone[ct][coff + 14 * i + 3] = DM_POLYTOPE_SEGMENT;
223: ex->cone[ct][coff + 14 * i + 4] = 0;
224: ex->cone[ct][coff + 14 * i + 5] = i + 1;
225: ex->cone[ct][coff + 14 * i + 6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
226: ex->cone[ct][coff + 14 * i + 7] = 1;
227: ex->cone[ct][coff + 14 * i + 8] = 0;
228: ex->cone[ct][coff + 14 * i + 9] = i;
229: ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
230: ex->cone[ct][coff + 14 * i + 11] = 1;
231: ex->cone[ct][coff + 14 * i + 12] = 1;
232: ex->cone[ct][coff + 14 * i + 13] = i;
233: ex->ornt[ct][ooff + 4 * i + 0] = 0;
234: ex->ornt[ct][ooff + 4 * i + 1] = 0;
235: ex->ornt[ct][ooff + 4 * i + 2] = 0;
236: ex->ornt[ct][ooff + 4 * i + 3] = 0;
237: } else {
238: ex->cone[ct][coff + 14 * i + 0] = DM_POLYTOPE_SEGMENT;
239: ex->cone[ct][coff + 14 * i + 1] = 0;
240: ex->cone[ct][coff + 14 * i + 2] = i;
241: ex->cone[ct][coff + 14 * i + 3] = DM_POLYTOPE_SEGMENT;
242: ex->cone[ct][coff + 14 * i + 4] = 1;
243: ex->cone[ct][coff + 14 * i + 5] = 1;
244: ex->cone[ct][coff + 14 * i + 6] = i;
245: ex->cone[ct][coff + 14 * i + 7] = DM_POLYTOPE_SEGMENT;
246: ex->cone[ct][coff + 14 * i + 8] = 0;
247: ex->cone[ct][coff + 14 * i + 9] = i + 1;
248: ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_SEGMENT;
249: ex->cone[ct][coff + 14 * i + 11] = 1;
250: ex->cone[ct][coff + 14 * i + 12] = 0;
251: ex->cone[ct][coff + 14 * i + 13] = i;
252: ex->ornt[ct][ooff + 4 * i + 0] = 0;
253: ex->ornt[ct][ooff + 4 * i + 1] = 0;
254: ex->ornt[ct][ooff + 4 * i + 2] = -1;
255: ex->ornt[ct][ooff + 4 * i + 3] = -1;
256: }
257: }
258: /* DM_POLYTOPE_TRIANGLE produces Nl+1 triangles and Nl triangular prisms/tensor triangular prisms */
259: ct = DM_POLYTOPE_TRIANGLE;
260: ex->Nt[ct] = 2;
261: Nc = 12 * (Nl + 1) + 18 * Nl;
262: No = 3 * (Nl + 1) + 5 * Nl;
263: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
264: ex->target[ct][0] = DM_POLYTOPE_TRIANGLE;
265: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM;
266: ex->size[ct][0] = Nl + 1;
267: ex->size[ct][1] = Nl;
268: /* cones for triangles */
269: for (i = 0; i < Nl + 1; ++i) {
270: ex->cone[ct][12 * i + 0] = DM_POLYTOPE_SEGMENT;
271: ex->cone[ct][12 * i + 1] = 1;
272: ex->cone[ct][12 * i + 2] = 0;
273: ex->cone[ct][12 * i + 3] = i;
274: ex->cone[ct][12 * i + 4] = DM_POLYTOPE_SEGMENT;
275: ex->cone[ct][12 * i + 5] = 1;
276: ex->cone[ct][12 * i + 6] = 1;
277: ex->cone[ct][12 * i + 7] = i;
278: ex->cone[ct][12 * i + 8] = DM_POLYTOPE_SEGMENT;
279: ex->cone[ct][12 * i + 9] = 1;
280: ex->cone[ct][12 * i + 10] = 2;
281: ex->cone[ct][12 * i + 11] = i;
282: }
283: for (i = 0; i < 3 * (Nl + 1); ++i) ex->ornt[ct][i] = 0;
284: /* cones for triangular prisms/tensor triangular prisms */
285: coff = 12 * (Nl + 1);
286: ooff = 3 * (Nl + 1);
287: for (i = 0; i < Nl; ++i) {
288: if (ex->useTensor) {
289: ex->cone[ct][coff + 18 * i + 0] = DM_POLYTOPE_TRIANGLE;
290: ex->cone[ct][coff + 18 * i + 1] = 0;
291: ex->cone[ct][coff + 18 * i + 2] = i;
292: ex->cone[ct][coff + 18 * i + 3] = DM_POLYTOPE_TRIANGLE;
293: ex->cone[ct][coff + 18 * i + 4] = 0;
294: ex->cone[ct][coff + 18 * i + 5] = i + 1;
295: ex->cone[ct][coff + 18 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
296: ex->cone[ct][coff + 18 * i + 7] = 1;
297: ex->cone[ct][coff + 18 * i + 8] = 0;
298: ex->cone[ct][coff + 18 * i + 9] = i;
299: ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
300: ex->cone[ct][coff + 18 * i + 11] = 1;
301: ex->cone[ct][coff + 18 * i + 12] = 1;
302: ex->cone[ct][coff + 18 * i + 13] = i;
303: ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
304: ex->cone[ct][coff + 18 * i + 15] = 1;
305: ex->cone[ct][coff + 18 * i + 16] = 2;
306: ex->cone[ct][coff + 18 * i + 17] = i;
307: ex->ornt[ct][ooff + 5 * i + 0] = 0;
308: ex->ornt[ct][ooff + 5 * i + 1] = 0;
309: ex->ornt[ct][ooff + 5 * i + 2] = 0;
310: ex->ornt[ct][ooff + 5 * i + 3] = 0;
311: ex->ornt[ct][ooff + 5 * i + 4] = 0;
312: } else {
313: ex->cone[ct][coff + 18 * i + 0] = DM_POLYTOPE_TRIANGLE;
314: ex->cone[ct][coff + 18 * i + 1] = 0;
315: ex->cone[ct][coff + 18 * i + 2] = i;
316: ex->cone[ct][coff + 18 * i + 3] = DM_POLYTOPE_TRIANGLE;
317: ex->cone[ct][coff + 18 * i + 4] = 0;
318: ex->cone[ct][coff + 18 * i + 5] = i + 1;
319: ex->cone[ct][coff + 18 * i + 6] = DM_POLYTOPE_QUADRILATERAL;
320: ex->cone[ct][coff + 18 * i + 7] = 1;
321: ex->cone[ct][coff + 18 * i + 8] = 0;
322: ex->cone[ct][coff + 18 * i + 9] = i;
323: ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
324: ex->cone[ct][coff + 18 * i + 11] = 1;
325: ex->cone[ct][coff + 18 * i + 12] = 1;
326: ex->cone[ct][coff + 18 * i + 13] = i;
327: ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
328: ex->cone[ct][coff + 18 * i + 15] = 1;
329: ex->cone[ct][coff + 18 * i + 16] = 2;
330: ex->cone[ct][coff + 18 * i + 17] = i;
331: ex->ornt[ct][ooff + 5 * i + 0] = -2;
332: ex->ornt[ct][ooff + 5 * i + 1] = 0;
333: ex->ornt[ct][ooff + 5 * i + 2] = 0;
334: ex->ornt[ct][ooff + 5 * i + 3] = 0;
335: ex->ornt[ct][ooff + 5 * i + 4] = 0;
336: }
337: }
338: /* DM_POLYTOPE_QUADRILATERAL produces Nl+1 quads and Nl hexes/tensor hexes */
339: ct = DM_POLYTOPE_QUADRILATERAL;
340: ex->Nt[ct] = 2;
341: Nc = 16 * (Nl + 1) + 22 * Nl;
342: No = 4 * (Nl + 1) + 6 * Nl;
343: PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
344: ex->target[ct][0] = DM_POLYTOPE_QUADRILATERAL;
345: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON;
346: ex->size[ct][0] = Nl + 1;
347: ex->size[ct][1] = Nl;
348: /* cones for quads */
349: for (i = 0; i < Nl + 1; ++i) {
350: ex->cone[ct][16 * i + 0] = DM_POLYTOPE_SEGMENT;
351: ex->cone[ct][16 * i + 1] = 1;
352: ex->cone[ct][16 * i + 2] = 0;
353: ex->cone[ct][16 * i + 3] = i;
354: ex->cone[ct][16 * i + 4] = DM_POLYTOPE_SEGMENT;
355: ex->cone[ct][16 * i + 5] = 1;
356: ex->cone[ct][16 * i + 6] = 1;
357: ex->cone[ct][16 * i + 7] = i;
358: ex->cone[ct][16 * i + 8] = DM_POLYTOPE_SEGMENT;
359: ex->cone[ct][16 * i + 9] = 1;
360: ex->cone[ct][16 * i + 10] = 2;
361: ex->cone[ct][16 * i + 11] = i;
362: ex->cone[ct][16 * i + 12] = DM_POLYTOPE_SEGMENT;
363: ex->cone[ct][16 * i + 13] = 1;
364: ex->cone[ct][16 * i + 14] = 3;
365: ex->cone[ct][16 * i + 15] = i;
366: }
367: for (i = 0; i < 4 * (Nl + 1); ++i) ex->ornt[ct][i] = 0;
368: /* cones for hexes/tensor hexes */
369: coff = 16 * (Nl + 1);
370: ooff = 4 * (Nl + 1);
371: for (i = 0; i < Nl; ++i) {
372: if (ex->useTensor) {
373: ex->cone[ct][coff + 22 * i + 0] = DM_POLYTOPE_QUADRILATERAL;
374: ex->cone[ct][coff + 22 * i + 1] = 0;
375: ex->cone[ct][coff + 22 * i + 2] = i;
376: ex->cone[ct][coff + 22 * i + 3] = DM_POLYTOPE_QUADRILATERAL;
377: ex->cone[ct][coff + 22 * i + 4] = 0;
378: ex->cone[ct][coff + 22 * i + 5] = i + 1;
379: ex->cone[ct][coff + 22 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
380: ex->cone[ct][coff + 22 * i + 7] = 1;
381: ex->cone[ct][coff + 22 * i + 8] = 0;
382: ex->cone[ct][coff + 22 * i + 9] = i;
383: ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
384: ex->cone[ct][coff + 22 * i + 11] = 1;
385: ex->cone[ct][coff + 22 * i + 12] = 1;
386: ex->cone[ct][coff + 22 * i + 13] = i;
387: ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
388: ex->cone[ct][coff + 22 * i + 15] = 1;
389: ex->cone[ct][coff + 22 * i + 16] = 2;
390: ex->cone[ct][coff + 22 * i + 17] = i;
391: ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
392: ex->cone[ct][coff + 22 * i + 19] = 1;
393: ex->cone[ct][coff + 22 * i + 20] = 3;
394: ex->cone[ct][coff + 22 * i + 21] = i;
395: ex->ornt[ct][ooff + 6 * i + 0] = 0;
396: ex->ornt[ct][ooff + 6 * i + 1] = 0;
397: ex->ornt[ct][ooff + 6 * i + 2] = 0;
398: ex->ornt[ct][ooff + 6 * i + 3] = 0;
399: ex->ornt[ct][ooff + 6 * i + 4] = 0;
400: ex->ornt[ct][ooff + 6 * i + 5] = 0;
401: } else {
402: ex->cone[ct][coff + 22 * i + 0] = DM_POLYTOPE_QUADRILATERAL;
403: ex->cone[ct][coff + 22 * i + 1] = 0;
404: ex->cone[ct][coff + 22 * i + 2] = i;
405: ex->cone[ct][coff + 22 * i + 3] = DM_POLYTOPE_QUADRILATERAL;
406: ex->cone[ct][coff + 22 * i + 4] = 0;
407: ex->cone[ct][coff + 22 * i + 5] = i + 1;
408: ex->cone[ct][coff + 22 * i + 6] = DM_POLYTOPE_QUADRILATERAL;
409: ex->cone[ct][coff + 22 * i + 7] = 1;
410: ex->cone[ct][coff + 22 * i + 8] = 0;
411: ex->cone[ct][coff + 22 * i + 9] = i;
412: ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
413: ex->cone[ct][coff + 22 * i + 11] = 1;
414: ex->cone[ct][coff + 22 * i + 12] = 2;
415: ex->cone[ct][coff + 22 * i + 13] = i;
416: ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
417: ex->cone[ct][coff + 22 * i + 15] = 1;
418: ex->cone[ct][coff + 22 * i + 16] = 1;
419: ex->cone[ct][coff + 22 * i + 17] = i;
420: ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_QUADRILATERAL;
421: ex->cone[ct][coff + 22 * i + 19] = 1;
422: ex->cone[ct][coff + 22 * i + 20] = 3;
423: ex->cone[ct][coff + 22 * i + 21] = i;
424: ex->ornt[ct][ooff + 6 * i + 0] = -2;
425: ex->ornt[ct][ooff + 6 * i + 1] = 0;
426: ex->ornt[ct][ooff + 6 * i + 2] = 0;
427: ex->ornt[ct][ooff + 6 * i + 3] = 0;
428: ex->ornt[ct][ooff + 6 * i + 4] = 0;
429: ex->ornt[ct][ooff + 6 * i + 5] = 1;
430: }
431: }
432: /* Layers positions */
433: if (!ex->Nth) {
434: if (ex->symmetric)
435: for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers - ex->thickness / 2;
436: else
437: for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers;
438: } else {
439: ex->thickness = 0.;
440: ex->layerPos[0] = 0.;
441: for (l = 0; l < ex->layers; ++l) {
442: const PetscReal t = ex->thicknesses[PetscMin(l, ex->Nth - 1)];
443: ex->layerPos[l + 1] = ex->layerPos[l] + t;
444: ex->thickness += t;
445: }
446: if (ex->symmetric)
447: for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] -= ex->thickness / 2.;
448: }
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr)
453: {
454: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
455: PetscInt ct;
457: PetscFunctionBegin;
458: if (ex->target) {
459: for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) PetscCall(PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct]));
460: }
461: PetscCall(PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt));
462: PetscCall(PetscFree(ex->layerPos));
463: PetscCall(PetscFree(ex));
464: PetscFunctionReturn(PETSC_SUCCESS);
465: }
467: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
468: {
469: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
470: DMLabel trType = tr->trType;
471: PetscInt rt;
473: PetscFunctionBeginHot;
474: *rnew = r;
475: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
476: if (!so) PetscFunctionReturn(PETSC_SUCCESS);
477: if (trType) {
478: PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
479: if (rt >= 100) PetscFunctionReturn(PETSC_SUCCESS);
480: }
481: if (ex->useTensor) {
482: switch (sct) {
483: case DM_POLYTOPE_POINT:
484: break;
485: case DM_POLYTOPE_SEGMENT:
486: switch (tct) {
487: case DM_POLYTOPE_SEGMENT:
488: break;
489: case DM_POLYTOPE_SEG_PRISM_TENSOR:
490: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
491: break;
492: default:
493: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
494: }
495: break;
496: // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
497: case DM_POLYTOPE_TRIANGLE:
498: break;
499: case DM_POLYTOPE_QUADRILATERAL:
500: break;
501: default:
502: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
503: }
504: } else {
505: switch (sct) {
506: case DM_POLYTOPE_POINT:
507: break;
508: case DM_POLYTOPE_SEGMENT:
509: switch (tct) {
510: case DM_POLYTOPE_SEGMENT:
511: break;
512: case DM_POLYTOPE_QUADRILATERAL:
513: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0);
514: break;
515: default:
516: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
517: }
518: break;
519: // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
520: case DM_POLYTOPE_TRIANGLE:
521: break;
522: case DM_POLYTOPE_QUADRILATERAL:
523: break;
524: default:
525: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
526: }
527: }
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
532: {
533: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
534: DMLabel trType = tr->trType;
535: PetscBool ignore = PETSC_FALSE, identity = PETSC_FALSE;
536: PetscInt val = 0;
538: PetscFunctionBegin;
539: if (trType) {
540: PetscCall(DMLabelGetValue(trType, p, &val));
541: identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE;
542: } else {
543: ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE;
544: }
545: if (rt) *rt = val;
546: if (ignore) {
547: /* Ignore cells that cannot be extruded */
548: *Nt = 0;
549: *target = NULL;
550: *size = NULL;
551: *cone = NULL;
552: *ornt = NULL;
553: } else if (identity) {
554: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
555: } else {
556: *Nt = ex->Nt[source];
557: *target = ex->target[source];
558: *size = ex->size[source];
559: *cone = ex->cone[source];
560: *ornt = ex->ornt[source];
561: }
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: /* Computes new vertex along normal */
566: static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
567: {
568: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
569: DM dm;
570: DMLabel active;
571: PetscReal ones2[2] = {0., 1.}, ones3[3] = {0., 0., 1.};
572: PetscReal normal[3] = {0., 0., 0.}, norm;
573: PetscBool computeNormal;
574: PetscInt dim, dEx = ex->cdimEx, cStart, cEnd, d;
576: PetscFunctionBeginHot;
577: PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]);
578: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
579: PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv);
580: PetscCheck(dE == ex->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %" PetscInt_FMT " != %" PetscInt_FMT " original dimension", dE, ex->cdim);
581: PetscCheck(dEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", dEx);
583: PetscCall(DMPlexTransformGetDM(tr, &dm));
584: PetscCall(DMGetDimension(dm, &dim));
585: PetscCall(DMPlexTransformGetActive(tr, &active));
586: computeNormal = dim != ex->cdim && !ex->useNormal ? PETSC_TRUE : PETSC_FALSE;
587: if (computeNormal) {
588: PetscInt *closure = NULL;
589: PetscInt closureSize, cl;
591: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
592: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure));
593: for (cl = 0; cl < closureSize * 2; cl += 2) {
594: if ((closure[cl] >= cStart) && (closure[cl] < cEnd)) {
595: PetscReal cnormal[3] = {0, 0, 0};
597: PetscCall(DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal));
598: for (d = 0; d < dEx; ++d) normal[d] += cnormal[d];
599: }
600: }
601: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure));
602: } else if (ex->useNormal) {
603: for (d = 0; d < dEx; ++d) normal[d] = ex->normal[d];
604: } else if (active) { // find an active point in the closure of p and use its coordinate normal as the extrusion direction
605: PetscInt *closure = NULL;
606: PetscInt closureSize, cl, pStart, pEnd;
608: PetscCall(DMPlexGetDepthStratum(dm, ex->cdimEx - 1, &pStart, &pEnd));
609: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure));
610: for (cl = 0; cl < closureSize * 2; cl += 2) {
611: if ((closure[cl] >= pStart) && (closure[cl] < pEnd)) {
612: PetscReal cnormal[3] = {0, 0, 0};
613: const PetscInt *supp;
614: PetscInt suppSize;
616: PetscCall(DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal));
617: PetscCall(DMPlexGetSupportSize(dm, closure[cl], &suppSize));
618: PetscCall(DMPlexGetSupport(dm, closure[cl], &supp));
619: // Only use external faces, so I can get the orientation from any cell
620: if (suppSize == 1) {
621: const PetscInt *cone, *ornt;
622: PetscInt coneSize, c;
624: PetscCall(DMPlexGetConeSize(dm, supp[0], &coneSize));
625: PetscCall(DMPlexGetCone(dm, supp[0], &cone));
626: PetscCall(DMPlexGetConeOrientation(dm, supp[0], &ornt));
627: for (c = 0; c < coneSize; ++c)
628: if (cone[c] == closure[cl]) break;
629: PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Asymmetry in cone/support");
630: if (ornt[c] < 0)
631: for (d = 0; d < dEx; ++d) cnormal[d] *= -1.;
632: for (d = 0; d < dEx; ++d) normal[d] += cnormal[d];
633: }
634: }
635: }
636: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure));
637: } else if (ex->cdimEx == 2) {
638: for (d = 0; d < dEx; ++d) normal[d] = ones2[d];
639: } else if (ex->cdimEx == 3) {
640: for (d = 0; d < dEx; ++d) normal[d] = ones3[d];
641: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
642: if (ex->normalFunc) {
643: PetscScalar n[3];
644: PetscReal x[3], dot;
646: for (d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]);
647: PetscCall((*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL));
648: for (dot = 0, d = 0; d < dEx; d++) dot += PetscRealPart(n[d]) * normal[d];
649: for (d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]);
650: }
652: for (d = 0, norm = 0.0; d < dEx; ++d) norm += PetscSqr(normal[d]);
653: for (d = 0; d < dEx; ++d) normal[d] *= norm == 0.0 ? 1.0 : 1. / PetscSqrtReal(norm);
654: for (d = 0; d < dEx; ++d) out[d] = normal[d] * ex->layerPos[r];
655: for (d = 0; d < ex->cdim; ++d) out[d] += in[d];
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr)
660: {
661: PetscFunctionBegin;
662: tr->ops->view = DMPlexTransformView_Extrude;
663: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude;
664: tr->ops->setup = DMPlexTransformSetUp_Extrude;
665: tr->ops->destroy = DMPlexTransformDestroy_Extrude;
666: tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude;
667: tr->ops->celltransform = DMPlexTransformCellTransform_Extrude;
668: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude;
669: tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude;
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr)
674: {
675: DMPlexTransform_Extrude *ex;
676: DM dm;
677: PetscInt dim;
679: PetscFunctionBegin;
681: PetscCall(PetscNew(&ex));
682: tr->data = ex;
683: ex->thickness = 1.;
684: ex->useTensor = PETSC_TRUE;
685: ex->symmetric = PETSC_FALSE;
686: ex->useNormal = PETSC_FALSE;
687: ex->layerPos = NULL;
688: PetscCall(DMPlexTransformGetDM(tr, &dm));
689: PetscCall(DMGetDimension(dm, &dim));
690: PetscCall(DMGetCoordinateDim(dm, &ex->cdim));
691: PetscCall(DMPlexTransformExtrudeSetLayers(tr, 1));
692: PetscCall(DMPlexTransformInitialize_Extrude(tr));
693: PetscFunctionReturn(PETSC_SUCCESS);
694: }
696: /*@
697: DMPlexTransformExtrudeGetLayers - Get the number of extruded layers.
699: Not Collective
701: Input Parameter:
702: . tr - The `DMPlexTransform`
704: Output Parameter:
705: . layers - The number of layers
707: Level: intermediate
709: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetLayers()`
710: @*/
711: PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers)
712: {
713: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
715: PetscFunctionBegin;
718: *layers = ex->layers;
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }
722: /*@
723: DMPlexTransformExtrudeSetLayers - Set the number of extruded layers.
725: Not Collective
727: Input Parameters:
728: + tr - The `DMPlexTransform`
729: - layers - The number of layers
731: Level: intermediate
733: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetLayers()`
734: @*/
735: PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers)
736: {
737: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
739: PetscFunctionBegin;
741: ex->layers = layers;
742: PetscCall(PetscFree(ex->layerPos));
743: PetscCall(PetscCalloc1(ex->layers + 1, &ex->layerPos));
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: /*@
748: DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers
750: Not Collective
752: Input Parameter:
753: . tr - The `DMPlexTransform`
755: Output Parameter:
756: . thickness - The total thickness of the layers
758: Level: intermediate
760: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`
761: @*/
762: PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness)
763: {
764: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
766: PetscFunctionBegin;
769: *thickness = ex->thickness;
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: /*@
774: DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers
776: Not Collective
778: Input Parameters:
779: + tr - The `DMPlexTransform`
780: - thickness - The total thickness of the layers
782: Level: intermediate
784: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetThickness()`
785: @*/
786: PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness)
787: {
788: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
790: PetscFunctionBegin;
792: PetscCheck(thickness > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double)thickness);
793: ex->thickness = thickness;
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: /*@
798: DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells
800: Not Collective
802: Input Parameter:
803: . tr - The `DMPlexTransform`
805: Output Parameter:
806: . useTensor - The flag to use tensor cells
808: Note:
809: This flag determines the orientation behavior of the created points.
811: For example, if tensor is `PETSC_TRUE`, then
812: .vb
813: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
814: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
815: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
816: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
817: .ve
819: Level: intermediate
821: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetTensor()`
822: @*/
823: PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
824: {
825: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
827: PetscFunctionBegin;
830: *useTensor = ex->useTensor;
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: /*@
835: DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells
837: Not Collective
839: Input Parameters:
840: + tr - The `DMPlexTransform`
841: - useTensor - The flag for tensor cells
843: Note:
844: This flag determines the orientation behavior of the created points
845: For example, if tensor is `PETSC_TRUE`, then
846: .vb
847: DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
848: DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
849: DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
850: DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
851: .ve
853: Level: intermediate
855: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetTensor()`
856: @*/
857: PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
858: {
859: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
861: PetscFunctionBegin;
863: ex->useTensor = useTensor;
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface
870: Not Collective
872: Input Parameter:
873: . tr - The `DMPlexTransform`
875: Output Parameter:
876: . symmetric - The flag to extrude symmetrically
878: Level: intermediate
880: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetSymmetric()`
881: @*/
882: PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric)
883: {
884: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
886: PetscFunctionBegin;
889: *symmetric = ex->symmetric;
890: PetscFunctionReturn(PETSC_SUCCESS);
891: }
893: /*@
894: DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface
896: Not Collective
898: Input Parameters:
899: + tr - The `DMPlexTransform`
900: - symmetric - The flag to extrude symmetrically
902: Level: intermediate
904: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetSymmetric()`
905: @*/
906: PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric)
907: {
908: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
910: PetscFunctionBegin;
912: ex->symmetric = symmetric;
913: PetscFunctionReturn(PETSC_SUCCESS);
914: }
916: /*@
917: DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector
919: Not Collective
921: Input Parameter:
922: . tr - The `DMPlexTransform`
924: Output Parameter:
925: . normal - The extrusion direction
927: Note:
928: The user passes in an array, which is filled by the library.
930: Level: intermediate
932: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetNormal()`
933: @*/
934: PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[])
935: {
936: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
937: PetscInt d;
939: PetscFunctionBegin;
941: if (ex->useNormal) {
942: for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d];
943: } else {
944: for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.;
945: }
946: PetscFunctionReturn(PETSC_SUCCESS);
947: }
949: /*@
950: DMPlexTransformExtrudeSetNormal - Set the extrusion normal
952: Not Collective
954: Input Parameters:
955: + tr - The `DMPlexTransform`
956: - normal - The extrusion direction
958: Level: intermediate
960: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`
961: @*/
962: PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[])
963: {
964: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
965: PetscInt d;
967: PetscFunctionBegin;
969: ex->useNormal = PETSC_TRUE;
970: for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d];
971: PetscFunctionReturn(PETSC_SUCCESS);
972: }
974: /*@C
975: DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal
977: Not Collective
979: Input Parameters:
980: + tr - The `DMPlexTransform`
981: - normalFunc - A function determining the extrusion direction
983: Calling sequence of `normalFunc`:
984: $ PetscErrorCode normalFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx)
985: + dim - The coordinate dimension of the original mesh (usually a surface)
986: . time - The current time, or 0.
987: . x - The location of the current normal, in the coordinate space of the original mesh
988: . r - The extrusion replica number (layer number) of this point
989: . u - The user provides the computed normal on output; the sign and magnitude is not significant
990: - ctx - An optional user context
992: Level: intermediate
994: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`
995: @*/
996: PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFunc normalFunc)
997: {
998: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1000: PetscFunctionBegin;
1002: ex->normalFunc = normalFunc;
1003: PetscFunctionReturn(PETSC_SUCCESS);
1004: }
1006: /*@
1007: DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer
1009: Not Collective
1011: Input Parameters:
1012: + tr - The `DMPlexTransform`
1013: . Nth - The number of thicknesses
1014: - thickness - The array of thicknesses
1016: Level: intermediate
1018: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeGetThickness()`
1019: @*/
1020: PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[])
1021: {
1022: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1023: PetscInt t;
1025: PetscFunctionBegin;
1027: PetscCheck(Nth > 0, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %" PetscInt_FMT " must be positive", Nth);
1028: ex->Nth = PetscMin(Nth, ex->layers);
1029: PetscCall(PetscFree(ex->thicknesses));
1030: PetscCall(PetscMalloc1(ex->Nth, &ex->thicknesses));
1031: for (t = 0; t < ex->Nth; ++t) {
1032: PetscCheck(thicknesses[t] > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %" PetscInt_FMT " must be positive", (double)thicknesses[t], t);
1033: ex->thicknesses[t] = thicknesses[t];
1034: }
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }