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