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: }