Actual source code: plextrextrude.c

  1: #include <petsc/private/dmplextransformimpl.h>

  3: const char *const PlexNormalAlgs[] = {"default", "input", "compute", "compute_bd"};

  5: static PetscErrorCode DMPlexTransformView_Extrude(DMPlexTransform tr, PetscViewer viewer)
  6: {
  7:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
  8:   PetscBool                isascii;

 10:   PetscFunctionBegin;
 13:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
 14:   if (isascii) {
 15:     DM          dm;
 16:     DMLabel     active;
 17:     PetscInt    dim;
 18:     const char *name;

 20:     PetscCall(PetscObjectGetName((PetscObject)tr, &name));
 21:     PetscCall(DMPlexTransformGetDM(tr, &dm));
 22:     PetscCall(DMGetDimension(dm, &dim));
 23:     PetscCall(DMPlexTransformGetActive(tr, &active));

 25:     PetscCall(PetscViewerASCIIPrintf(viewer, "Extrusion transformation %s\n", name ? name : ""));
 26:     PetscCall(PetscViewerASCIIPrintf(viewer, "  number of layers: %" PetscInt_FMT "\n", ex->layers));
 27:     PetscCall(PetscViewerASCIIPrintf(viewer, "  create tensor cells: %s\n", ex->useTensor ? "YES" : "NO"));
 28:     if (ex->periodic) PetscCall(PetscViewerASCIIPrintf(viewer, "  periodic\n"));
 29:     PetscCall(PetscViewerASCIIPrintf(viewer, "  normal algorithm: %s\n", PlexNormalAlgs[ex->normalAlg]));
 30:     if (ex->normalFunc) PetscCall(PetscViewerASCIIPrintf(viewer, "  normal modified by user function\n"));
 31:   } else {
 32:     SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
 33:   }
 34:   PetscFunctionReturn(PETSC_SUCCESS);
 35: }

 37: static PetscErrorCode DMPlexTransformSetFromOptions_Extrude(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
 38: {
 39:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
 40:   PetscReal                th, normal[3], *thicknesses;
 41:   PetscInt                 nl, Nc;
 42:   PetscBool                tensor, sym, per, flg;
 43:   char                     funcname[PETSC_MAX_PATH_LEN];

 45:   PetscFunctionBegin;
 46:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Extrusion Options");
 47:   PetscCall(PetscOptionsBoundedInt("-dm_plex_transform_extrude_layers", "Number of layers to extrude", "", ex->layers, &nl, &flg, 1));
 48:   if (flg) PetscCall(DMPlexTransformExtrudeSetLayers(tr, nl));
 49:   PetscCall(PetscOptionsReal("-dm_plex_transform_extrude_thickness", "Total thickness of extruded layers", "", ex->thickness, &th, &flg));
 50:   if (flg) PetscCall(DMPlexTransformExtrudeSetThickness(tr, th));
 51:   PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg));
 52:   if (flg) PetscCall(DMPlexTransformExtrudeSetTensor(tr, tensor));
 53:   PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_symmetric", "Extrude layers symmetrically about the surface", "", ex->symmetric, &sym, &flg));
 54:   if (flg) PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, sym));
 55:   PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_periodic", "Extrude layers periodically about the surface", "", ex->periodic, &per, &flg));
 56:   if (flg) PetscCall(DMPlexTransformExtrudeSetPeriodic(tr, per));
 57:   Nc = 3;
 58:   PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_normal", "Input normal vector for extrusion", "DMPlexTransformExtrudeSetNormal", normal, &Nc, &flg));
 59:   if (flg) {
 60:     // Extrusion dimension might not yet be determined
 61:     PetscCheck(!ex->cdimEx || Nc == ex->cdimEx, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_SIZ, "Input normal has size %" PetscInt_FMT " != %" PetscInt_FMT " extruded coordinate dimension", Nc, ex->cdimEx);
 62:     PetscCall(DMPlexTransformExtrudeSetNormal(tr, normal));
 63:   }
 64:   PetscCall(PetscOptionsString("-dm_plex_transform_extrude_normal_function", "Function to determine normal vector", "DMPlexTransformExtrudeSetNormalFunction", funcname, funcname, sizeof(funcname), &flg));
 65:   if (flg) {
 66:     PetscSimplePointFn *normalFunc;

 68:     PetscCall(PetscDLSym(NULL, funcname, (void **)&normalFunc));
 69:     PetscCall(DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc));
 70:   }
 71:   nl = ex->layers;
 72:   PetscCall(PetscMalloc1(nl, &thicknesses));
 73:   PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_thicknesses", "Thickness of each individual extruded layer", "", thicknesses, &nl, &flg));
 74:   if (flg) {
 75:     PetscCheck(nl, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one thickness for -dm_plex_transform_extrude_thicknesses");
 76:     PetscCall(DMPlexTransformExtrudeSetThicknesses(tr, nl, thicknesses));
 77:   }
 78:   PetscCall(PetscFree(thicknesses));
 79:   PetscOptionsHeadEnd();
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }

 83: /* Determine the implicit dimension pre-extrusion (either the implicit dimension of the DM or of a point in the active set for the transform).
 84:    If that dimension is the same as the current coordinate dimension (ex->dim), the extruded mesh will have a coordinate dimension one greater;
 85:    Otherwise the coordinate dimension will be kept. */
 86: static PetscErrorCode DMPlexTransformExtrudeComputeExtrusionDim(DMPlexTransform tr)
 87: {
 88:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
 89:   DM                       dm;
 90:   DMLabel                  active;
 91:   PetscInt                 dim, dimExtPoint, dimExtPointG;

 93:   PetscFunctionBegin;
 94:   PetscCall(DMPlexTransformGetDM(tr, &dm));
 95:   PetscCall(DMGetDimension(dm, &dim));
 96:   PetscCall(DMPlexTransformGetActive(tr, &active));
 97:   if (active) {
 98:     IS              valueIS, pointIS;
 99:     const PetscInt *values, *points;
100:     DMPolytopeType  ct;
101:     PetscInt        Nv, Np;

103:     dimExtPoint = 0;
104:     PetscCall(DMLabelGetValueIS(active, &valueIS));
105:     PetscCall(ISGetLocalSize(valueIS, &Nv));
106:     PetscCall(ISGetIndices(valueIS, &values));
107:     for (PetscInt v = 0; v < Nv; ++v) {
108:       PetscCall(DMLabelGetStratumIS(active, values[v], &pointIS));
109:       PetscCall(ISGetLocalSize(pointIS, &Np));
110:       PetscCall(ISGetIndices(pointIS, &points));
111:       for (PetscInt p = 0; p < Np; ++p) {
112:         PetscCall(DMPlexGetCellType(dm, points[p], &ct));
113:         dimExtPoint = PetscMax(dimExtPoint, DMPolytopeTypeGetDim(ct));
114:       }
115:       PetscCall(ISRestoreIndices(pointIS, &points));
116:       PetscCall(ISDestroy(&pointIS));
117:     }
118:     PetscCall(ISRestoreIndices(valueIS, &values));
119:     PetscCall(ISDestroy(&valueIS));
120:   } else dimExtPoint = dim;
121:   PetscCallMPI(MPIU_Allreduce(&dimExtPoint, &dimExtPointG, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)tr)));
122:   ex->dimEx  = PetscMax(dim, dimExtPointG + 1);
123:   ex->cdimEx = ex->cdim == dimExtPointG ? ex->cdim + 1 : ex->cdim;
124:   PetscCheck(ex->dimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Topological dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->dimEx);
125:   PetscCheck(ex->cdimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->cdimEx);
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm)
130: {
131:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
132:   PetscInt                 dim;

134:   PetscFunctionBegin;
135:   PetscCall(DMGetDimension(dm, &dim));
136:   PetscCall(DMSetDimension(tdm, ex->dimEx));
137:   PetscCall(DMSetCoordinateDim(tdm, ex->cdimEx));
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: /* DM_POLYTOPE_POINT produces
142:      Nl+1 points, or Nl points periodically and
143:      Nl segments, or tensor segments
144: */
145: static PetscErrorCode DMPlexTransformExtrudeSetUp_Point(DMPlexTransform_Extrude *ex, PetscInt Nl)
146: {
147:   const DMPolytopeType ct = DM_POLYTOPE_POINT;
148:   const PetscInt       Np = ex->periodic ? Nl : Nl + 1;
149:   PetscInt             Nc, No;

151:   PetscFunctionBegin;
152:   ex->Nt[ct] = 2;
153:   Nc         = 6 * Nl;
154:   No         = 2 * Nl;
155:   PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
156:   ex->target[ct][0] = DM_POLYTOPE_POINT;
157:   ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
158:   ex->size[ct][0]   = Np;
159:   ex->size[ct][1]   = Nl;
160:   /*   cones for segments/tensor segments */
161:   for (PetscInt i = 0; i < Nl; ++i) {
162:     ex->cone[ct][6 * i + 0] = DM_POLYTOPE_POINT;
163:     ex->cone[ct][6 * i + 1] = 0;
164:     ex->cone[ct][6 * i + 2] = i;
165:     ex->cone[ct][6 * i + 3] = DM_POLYTOPE_POINT;
166:     ex->cone[ct][6 * i + 4] = 0;
167:     ex->cone[ct][6 * i + 5] = (i + 1) % Np;
168:   }
169:   for (PetscInt i = 0; i < No; ++i) ex->ornt[ct][i] = 0;
170:   PetscFunctionReturn(PETSC_SUCCESS);
171: }

173: /* DM_POLYTOPE_SEGMENT produces
174:      Nl+1 segments, or Nl segments periodically and
175:      Nl quads, or tensor quads
176: */
177: static PetscErrorCode DMPlexTransformExtrudeSetUp_Segment(DMPlexTransform_Extrude *ex, PetscInt Nl)
178: {
179:   const DMPolytopeType ct = DM_POLYTOPE_SEGMENT;
180:   const PetscInt       Np = ex->periodic ? Nl : Nl + 1;
181:   PetscInt             Nc, No, coff, ooff;

183:   PetscFunctionBegin;
184:   ex->Nt[ct] = 2;
185:   Nc         = 8 * Np + 14 * Nl;
186:   No         = 2 * Np + 4 * Nl;
187:   PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
188:   ex->target[ct][0] = DM_POLYTOPE_SEGMENT;
189:   ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
190:   ex->size[ct][0]   = Np;
191:   ex->size[ct][1]   = Nl;
192:   /*   cones for segments */
193:   for (PetscInt i = 0; i < Np; ++i) {
194:     ex->cone[ct][8 * i + 0] = DM_POLYTOPE_POINT;
195:     ex->cone[ct][8 * i + 1] = 1;
196:     ex->cone[ct][8 * i + 2] = 0;
197:     ex->cone[ct][8 * i + 3] = i;
198:     ex->cone[ct][8 * i + 4] = DM_POLYTOPE_POINT;
199:     ex->cone[ct][8 * i + 5] = 1;
200:     ex->cone[ct][8 * i + 6] = 1;
201:     ex->cone[ct][8 * i + 7] = i;
202:   }
203:   for (PetscInt i = 0; i < 2 * Np; ++i) ex->ornt[ct][i] = 0;
204:   /*   cones for quads/tensor quads */
205:   coff = 8 * Np;
206:   ooff = 2 * Np;
207:   for (PetscInt i = 0; i < Nl; ++i) {
208:     if (ex->useTensor) {
209:       ex->cone[ct][coff + 14 * i + 0]  = DM_POLYTOPE_SEGMENT;
210:       ex->cone[ct][coff + 14 * i + 1]  = 0;
211:       ex->cone[ct][coff + 14 * i + 2]  = i;
212:       ex->cone[ct][coff + 14 * i + 3]  = DM_POLYTOPE_SEGMENT;
213:       ex->cone[ct][coff + 14 * i + 4]  = 0;
214:       ex->cone[ct][coff + 14 * i + 5]  = (i + 1) % Np;
215:       ex->cone[ct][coff + 14 * i + 6]  = DM_POLYTOPE_POINT_PRISM_TENSOR;
216:       ex->cone[ct][coff + 14 * i + 7]  = 1;
217:       ex->cone[ct][coff + 14 * i + 8]  = 0;
218:       ex->cone[ct][coff + 14 * i + 9]  = i;
219:       ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
220:       ex->cone[ct][coff + 14 * i + 11] = 1;
221:       ex->cone[ct][coff + 14 * i + 12] = 1;
222:       ex->cone[ct][coff + 14 * i + 13] = i;
223:       ex->ornt[ct][ooff + 4 * i + 0]   = 0;
224:       ex->ornt[ct][ooff + 4 * i + 1]   = 0;
225:       ex->ornt[ct][ooff + 4 * i + 2]   = 0;
226:       ex->ornt[ct][ooff + 4 * i + 3]   = 0;
227:     } else {
228:       ex->cone[ct][coff + 14 * i + 0]  = DM_POLYTOPE_SEGMENT;
229:       ex->cone[ct][coff + 14 * i + 1]  = 0;
230:       ex->cone[ct][coff + 14 * i + 2]  = i;
231:       ex->cone[ct][coff + 14 * i + 3]  = DM_POLYTOPE_SEGMENT;
232:       ex->cone[ct][coff + 14 * i + 4]  = 1;
233:       ex->cone[ct][coff + 14 * i + 5]  = 1;
234:       ex->cone[ct][coff + 14 * i + 6]  = i;
235:       ex->cone[ct][coff + 14 * i + 7]  = DM_POLYTOPE_SEGMENT;
236:       ex->cone[ct][coff + 14 * i + 8]  = 0;
237:       ex->cone[ct][coff + 14 * i + 9]  = (i + 1) % Np;
238:       ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_SEGMENT;
239:       ex->cone[ct][coff + 14 * i + 11] = 1;
240:       ex->cone[ct][coff + 14 * i + 12] = 0;
241:       ex->cone[ct][coff + 14 * i + 13] = i;
242:       ex->ornt[ct][ooff + 4 * i + 0]   = 0;
243:       ex->ornt[ct][ooff + 4 * i + 1]   = 0;
244:       ex->ornt[ct][ooff + 4 * i + 2]   = -1;
245:       ex->ornt[ct][ooff + 4 * i + 3]   = -1;
246:     }
247:   }
248:   PetscFunctionReturn(PETSC_SUCCESS);
249: }

251: /* DM_POLYTOPE_TRIANGLE produces
252:      Nl+1 triangles, or Nl triangles periodically and
253:      Nl triangular prisms/tensor triangular prisms
254: */
255: static PetscErrorCode DMPlexTransformExtrudeSetUp_Triangle(DMPlexTransform_Extrude *ex, PetscInt Nl)
256: {
257:   const DMPolytopeType ct = DM_POLYTOPE_TRIANGLE;
258:   const PetscInt       Np = ex->periodic ? Nl : Nl + 1;
259:   PetscInt             Nc, No, coff, ooff;

261:   PetscFunctionBegin;
262:   ex->Nt[ct] = 2;
263:   Nc         = 12 * Np + 18 * Nl;
264:   No         = 3 * Np + 5 * Nl;
265:   PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
266:   ex->target[ct][0] = DM_POLYTOPE_TRIANGLE;
267:   ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM;
268:   ex->size[ct][0]   = Np;
269:   ex->size[ct][1]   = Nl;
270:   /*   cones for triangles */
271:   for (PetscInt i = 0; i < Np; ++i) {
272:     ex->cone[ct][12 * i + 0]  = DM_POLYTOPE_SEGMENT;
273:     ex->cone[ct][12 * i + 1]  = 1;
274:     ex->cone[ct][12 * i + 2]  = 0;
275:     ex->cone[ct][12 * i + 3]  = i;
276:     ex->cone[ct][12 * i + 4]  = DM_POLYTOPE_SEGMENT;
277:     ex->cone[ct][12 * i + 5]  = 1;
278:     ex->cone[ct][12 * i + 6]  = 1;
279:     ex->cone[ct][12 * i + 7]  = i;
280:     ex->cone[ct][12 * i + 8]  = DM_POLYTOPE_SEGMENT;
281:     ex->cone[ct][12 * i + 9]  = 1;
282:     ex->cone[ct][12 * i + 10] = 2;
283:     ex->cone[ct][12 * i + 11] = i;
284:   }
285:   for (PetscInt i = 0; i < 3 * Np; ++i) ex->ornt[ct][i] = 0;
286:   /*   cones for triangular prisms/tensor triangular prisms */
287:   coff = 12 * Np;
288:   ooff = 3 * Np;
289:   for (PetscInt i = 0; i < Nl; ++i) {
290:     if (ex->useTensor) {
291:       ex->cone[ct][coff + 18 * i + 0]  = DM_POLYTOPE_TRIANGLE;
292:       ex->cone[ct][coff + 18 * i + 1]  = 0;
293:       ex->cone[ct][coff + 18 * i + 2]  = i;
294:       ex->cone[ct][coff + 18 * i + 3]  = DM_POLYTOPE_TRIANGLE;
295:       ex->cone[ct][coff + 18 * i + 4]  = 0;
296:       ex->cone[ct][coff + 18 * i + 5]  = (i + 1) % Np;
297:       ex->cone[ct][coff + 18 * i + 6]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
298:       ex->cone[ct][coff + 18 * i + 7]  = 1;
299:       ex->cone[ct][coff + 18 * i + 8]  = 0;
300:       ex->cone[ct][coff + 18 * i + 9]  = i;
301:       ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
302:       ex->cone[ct][coff + 18 * i + 11] = 1;
303:       ex->cone[ct][coff + 18 * i + 12] = 1;
304:       ex->cone[ct][coff + 18 * i + 13] = i;
305:       ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
306:       ex->cone[ct][coff + 18 * i + 15] = 1;
307:       ex->cone[ct][coff + 18 * i + 16] = 2;
308:       ex->cone[ct][coff + 18 * i + 17] = i;
309:       ex->ornt[ct][ooff + 5 * i + 0]   = 0;
310:       ex->ornt[ct][ooff + 5 * i + 1]   = 0;
311:       ex->ornt[ct][ooff + 5 * i + 2]   = 0;
312:       ex->ornt[ct][ooff + 5 * i + 3]   = 0;
313:       ex->ornt[ct][ooff + 5 * i + 4]   = 0;
314:     } else {
315:       ex->cone[ct][coff + 18 * i + 0]  = DM_POLYTOPE_TRIANGLE;
316:       ex->cone[ct][coff + 18 * i + 1]  = 0;
317:       ex->cone[ct][coff + 18 * i + 2]  = i;
318:       ex->cone[ct][coff + 18 * i + 3]  = DM_POLYTOPE_TRIANGLE;
319:       ex->cone[ct][coff + 18 * i + 4]  = 0;
320:       ex->cone[ct][coff + 18 * i + 5]  = (i + 1) % Np;
321:       ex->cone[ct][coff + 18 * i + 6]  = DM_POLYTOPE_QUADRILATERAL;
322:       ex->cone[ct][coff + 18 * i + 7]  = 1;
323:       ex->cone[ct][coff + 18 * i + 8]  = 0;
324:       ex->cone[ct][coff + 18 * i + 9]  = i;
325:       ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
326:       ex->cone[ct][coff + 18 * i + 11] = 1;
327:       ex->cone[ct][coff + 18 * i + 12] = 1;
328:       ex->cone[ct][coff + 18 * i + 13] = i;
329:       ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
330:       ex->cone[ct][coff + 18 * i + 15] = 1;
331:       ex->cone[ct][coff + 18 * i + 16] = 2;
332:       ex->cone[ct][coff + 18 * i + 17] = i;
333:       ex->ornt[ct][ooff + 5 * i + 0]   = -2;
334:       ex->ornt[ct][ooff + 5 * i + 1]   = 0;
335:       ex->ornt[ct][ooff + 5 * i + 2]   = 0;
336:       ex->ornt[ct][ooff + 5 * i + 3]   = 0;
337:       ex->ornt[ct][ooff + 5 * i + 4]   = 0;
338:     }
339:   }
340:   PetscFunctionReturn(PETSC_SUCCESS);
341: }

343: /* DM_POLYTOPE_QUADRILATERAL produces
344:      Nl+1 quads, or Nl quads periodically and
345:      Nl hexes/tensor hexes
346: */
347: static PetscErrorCode DMPlexTransformExtrudeSetUp_Quadrilateral(DMPlexTransform_Extrude *ex, PetscInt Nl)
348: {
349:   const DMPolytopeType ct = DM_POLYTOPE_QUADRILATERAL;
350:   const PetscInt       Np = ex->periodic ? Nl : Nl + 1;
351:   PetscInt             Nc, No, coff, ooff;

353:   PetscFunctionBegin;
354:   ex->Nt[ct] = 2;
355:   Nc         = 16 * Np + 22 * Nl;
356:   No         = 4 * Np + 6 * Nl;
357:   PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]));
358:   ex->target[ct][0] = DM_POLYTOPE_QUADRILATERAL;
359:   ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON;
360:   ex->size[ct][0]   = Np;
361:   ex->size[ct][1]   = Nl;
362:   /*   cones for quads */
363:   for (PetscInt i = 0; i < Np; ++i) {
364:     ex->cone[ct][16 * i + 0]  = DM_POLYTOPE_SEGMENT;
365:     ex->cone[ct][16 * i + 1]  = 1;
366:     ex->cone[ct][16 * i + 2]  = 0;
367:     ex->cone[ct][16 * i + 3]  = i;
368:     ex->cone[ct][16 * i + 4]  = DM_POLYTOPE_SEGMENT;
369:     ex->cone[ct][16 * i + 5]  = 1;
370:     ex->cone[ct][16 * i + 6]  = 1;
371:     ex->cone[ct][16 * i + 7]  = i;
372:     ex->cone[ct][16 * i + 8]  = DM_POLYTOPE_SEGMENT;
373:     ex->cone[ct][16 * i + 9]  = 1;
374:     ex->cone[ct][16 * i + 10] = 2;
375:     ex->cone[ct][16 * i + 11] = i;
376:     ex->cone[ct][16 * i + 12] = DM_POLYTOPE_SEGMENT;
377:     ex->cone[ct][16 * i + 13] = 1;
378:     ex->cone[ct][16 * i + 14] = 3;
379:     ex->cone[ct][16 * i + 15] = i;
380:   }
381:   for (PetscInt i = 0; i < 4 * Np; ++i) ex->ornt[ct][i] = 0;
382:   /*   cones for hexes/tensor hexes */
383:   coff = 16 * Np;
384:   ooff = 4 * Np;
385:   for (PetscInt i = 0; i < Nl; ++i) {
386:     if (ex->useTensor) {
387:       ex->cone[ct][coff + 22 * i + 0]  = DM_POLYTOPE_QUADRILATERAL;
388:       ex->cone[ct][coff + 22 * i + 1]  = 0;
389:       ex->cone[ct][coff + 22 * i + 2]  = i;
390:       ex->cone[ct][coff + 22 * i + 3]  = DM_POLYTOPE_QUADRILATERAL;
391:       ex->cone[ct][coff + 22 * i + 4]  = 0;
392:       ex->cone[ct][coff + 22 * i + 5]  = (i + 1) % Np;
393:       ex->cone[ct][coff + 22 * i + 6]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
394:       ex->cone[ct][coff + 22 * i + 7]  = 1;
395:       ex->cone[ct][coff + 22 * i + 8]  = 0;
396:       ex->cone[ct][coff + 22 * i + 9]  = i;
397:       ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
398:       ex->cone[ct][coff + 22 * i + 11] = 1;
399:       ex->cone[ct][coff + 22 * i + 12] = 1;
400:       ex->cone[ct][coff + 22 * i + 13] = i;
401:       ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
402:       ex->cone[ct][coff + 22 * i + 15] = 1;
403:       ex->cone[ct][coff + 22 * i + 16] = 2;
404:       ex->cone[ct][coff + 22 * i + 17] = i;
405:       ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
406:       ex->cone[ct][coff + 22 * i + 19] = 1;
407:       ex->cone[ct][coff + 22 * i + 20] = 3;
408:       ex->cone[ct][coff + 22 * i + 21] = i;
409:       ex->ornt[ct][ooff + 6 * i + 0]   = 0;
410:       ex->ornt[ct][ooff + 6 * i + 1]   = 0;
411:       ex->ornt[ct][ooff + 6 * i + 2]   = 0;
412:       ex->ornt[ct][ooff + 6 * i + 3]   = 0;
413:       ex->ornt[ct][ooff + 6 * i + 4]   = 0;
414:       ex->ornt[ct][ooff + 6 * i + 5]   = 0;
415:     } else {
416:       ex->cone[ct][coff + 22 * i + 0]  = DM_POLYTOPE_QUADRILATERAL;
417:       ex->cone[ct][coff + 22 * i + 1]  = 0;
418:       ex->cone[ct][coff + 22 * i + 2]  = i;
419:       ex->cone[ct][coff + 22 * i + 3]  = DM_POLYTOPE_QUADRILATERAL;
420:       ex->cone[ct][coff + 22 * i + 4]  = 0;
421:       ex->cone[ct][coff + 22 * i + 5]  = (i + 1) % Np;
422:       ex->cone[ct][coff + 22 * i + 6]  = DM_POLYTOPE_QUADRILATERAL;
423:       ex->cone[ct][coff + 22 * i + 7]  = 1;
424:       ex->cone[ct][coff + 22 * i + 8]  = 0;
425:       ex->cone[ct][coff + 22 * i + 9]  = i;
426:       ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
427:       ex->cone[ct][coff + 22 * i + 11] = 1;
428:       ex->cone[ct][coff + 22 * i + 12] = 2;
429:       ex->cone[ct][coff + 22 * i + 13] = i;
430:       ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
431:       ex->cone[ct][coff + 22 * i + 15] = 1;
432:       ex->cone[ct][coff + 22 * i + 16] = 1;
433:       ex->cone[ct][coff + 22 * i + 17] = i;
434:       ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_QUADRILATERAL;
435:       ex->cone[ct][coff + 22 * i + 19] = 1;
436:       ex->cone[ct][coff + 22 * i + 20] = 3;
437:       ex->cone[ct][coff + 22 * i + 21] = i;
438:       ex->ornt[ct][ooff + 6 * i + 0]   = -2;
439:       ex->ornt[ct][ooff + 6 * i + 1]   = 0;
440:       ex->ornt[ct][ooff + 6 * i + 2]   = 0;
441:       ex->ornt[ct][ooff + 6 * i + 3]   = 0;
442:       ex->ornt[ct][ooff + 6 * i + 4]   = 0;
443:       ex->ornt[ct][ooff + 6 * i + 5]   = 1;
444:     }
445:   }
446:   PetscFunctionReturn(PETSC_SUCCESS);
447: }

449: /*
450:   The refine types for extrusion are:

452:   ct:       For any normally extruded point
453:   ct + 100: For any point which should just return itself
454: */
455: static PetscErrorCode DMPlexTransformSetUp_Extrude(DMPlexTransform tr)
456: {
457:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
458:   DM                       dm;
459:   DMLabel                  active;
460:   PetscInt                 Nl = ex->layers, l, ict, dim;

462:   PetscFunctionBegin;
463:   PetscCall(DMPlexTransformExtrudeComputeExtrusionDim(tr));
464:   PetscCall(DMPlexTransformGetDM(tr, &dm));
465:   PetscCall(DMGetDimension(dm, &dim));
466:   PetscCall(DMPlexTransformGetActive(tr, &active));
467:   if (active) {
468:     DMLabel  celltype;
469:     PetscInt pStart, pEnd, p;

471:     PetscCall(DMPlexGetCellTypeLabel(dm, &celltype));
472:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
473:     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
474:     for (p = pStart; p < pEnd; ++p) {
475:       PetscInt ct, val;

477:       PetscCall(DMLabelGetValue(celltype, p, &ct));
478:       PetscCall(DMLabelGetValue(active, p, &val));
479:       if (val < 0) {
480:         PetscCall(DMLabelSetValue(tr->trType, p, ct + 100));
481:       } else {
482:         PetscCall(DMLabelSetValue(tr->trType, p, ct));
483:       }
484:     }
485:   }
486:   if (ex->normalAlg != NORMAL_INPUT) {
487:     if (dim != ex->cdim) ex->normalAlg = NORMAL_COMPUTE;
488:     else if (active) ex->normalAlg = NORMAL_COMPUTE_BD;
489:   }
490:   // Need this to determine face sharing
491:   PetscMPIInt size;

493:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
494:   if (size > 1) {
495:     PetscSF  sf;
496:     PetscInt Nr;

498:     PetscCall(DMGetPointSF(dm, &sf));
499:     PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
500:     if (Nr >= 0) {
501:       PetscCall(PetscSFComputeDegreeBegin(sf, &ex->degree));
502:       PetscCall(PetscSFComputeDegreeEnd(sf, &ex->degree));
503:     }
504:   }
505:   // Create normal field
506:   if (ex->normalAlg == NORMAL_COMPUTE_BD) {
507:     PetscSection s;
508:     PetscInt     vStart, vEnd;

510:     PetscMPIInt rank;
511:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

513:     PetscCall(DMClone(dm, &ex->dmNormal));
514:     PetscCall(DMGetLocalSection(ex->dmNormal, &s));
515:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
516:     PetscCall(PetscSectionSetNumFields(s, 1));
517:     PetscCall(PetscSectionSetChart(s, vStart, vEnd));
518:     for (PetscInt v = vStart; v < vEnd; ++v) PetscCall(PetscSectionSetDof(s, v, ex->cdimEx));
519:     PetscCall(PetscSectionSetUp(s));
520:     PetscCall(DMCreateLocalVector(ex->dmNormal, &ex->vecNormal));
521:     PetscCall(PetscObjectSetName((PetscObject)ex->vecNormal, "Normal Field"));

523:     // find an active point in the closure of v and use its coordinate normal as the extrusion direction
524:     PetscSF         sf;
525:     const PetscInt *leaves;
526:     PetscScalar    *a, *normal;
527:     PetscInt        Nl;

529:     PetscCall(DMGetPointSF(ex->dmNormal, &sf));
530:     PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
531:     PetscCall(VecGetArrayWrite(ex->vecNormal, &a));
532:     for (PetscInt v = vStart; v < vEnd; ++v) {
533:       PetscInt *star = NULL;
534:       PetscInt  starSize, pStart, pEnd;

536:       PetscCall(DMPlexGetDepthStratum(ex->dmNormal, ex->cdimEx - 1, &pStart, &pEnd));
537:       PetscCall(DMPlexGetTransitiveClosure(ex->dmNormal, v, PETSC_FALSE, &starSize, &star));
538:       PetscCall(DMPlexPointLocalRef(ex->dmNormal, v, a, &normal));
539:       for (PetscInt st = 0; st < starSize * 2; st += 2) {
540:         const PetscInt face = star[st];
541:         if ((face >= pStart) && (face < pEnd)) {
542:           PetscReal       cnormal[3] = {0, 0, 0};
543:           const PetscInt *supp;
544:           PetscInt        suppSize, floc = -1;
545:           PetscBool       shared;

547:           PetscCall(DMPlexComputeCellGeometryFVM(ex->dmNormal, face, NULL, NULL, cnormal));
548:           PetscCall(DMPlexGetSupportSize(ex->dmNormal, face, &suppSize));
549:           PetscCall(DMPlexGetSupport(ex->dmNormal, face, &supp));
550:           // Only use external faces, so I can get the orientation from any cell
551:           if (leaves) PetscCall(PetscFindInt(face, Nl, leaves, &floc));
552:           shared = floc >= 0 || (ex->degree && ex->degree[face]) ? PETSC_TRUE : PETSC_FALSE;
553:           if (suppSize == 1 && !shared) {
554:             const PetscInt *cone, *ornt;
555:             PetscInt        coneSize, c;

557:             PetscCall(DMPlexGetConeSize(ex->dmNormal, supp[0], &coneSize));
558:             PetscCall(DMPlexGetCone(ex->dmNormal, supp[0], &cone));
559:             PetscCall(DMPlexGetConeOrientation(ex->dmNormal, supp[0], &ornt));
560:             for (c = 0; c < coneSize; ++c)
561:               if (cone[c] == face) break;
562:             PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Asymmetry in cone/support");
563:             if (ornt[c] < 0)
564:               for (PetscInt d = 0; d < ex->cdimEx; ++d) cnormal[d] *= -1.;
565:             for (PetscInt d = 0; d < ex->cdimEx; ++d) normal[d] += cnormal[d];
566:           }
567:         }
568:       }
569:       PetscCall(DMPlexRestoreTransitiveClosure(ex->dmNormal, v, PETSC_FALSE, &starSize, &star));
570:     }
571:     PetscCall(VecRestoreArrayWrite(ex->vecNormal, &a));

573:     Vec g;

575:     PetscCall(DMGetGlobalVector(ex->dmNormal, &g));
576:     PetscCall(VecSet(g, 0.));
577:     PetscCall(DMLocalToGlobal(ex->dmNormal, ex->vecNormal, ADD_VALUES, g));
578:     PetscCall(DMGlobalToLocal(ex->dmNormal, g, INSERT_VALUES, ex->vecNormal));
579:     PetscCall(DMRestoreGlobalVector(ex->dmNormal, &g));
580:   }
581:   PetscCall(PetscMalloc5(DM_NUM_POLYTOPES, &ex->Nt, DM_NUM_POLYTOPES, &ex->target, DM_NUM_POLYTOPES, &ex->size, DM_NUM_POLYTOPES, &ex->cone, DM_NUM_POLYTOPES, &ex->ornt));
582:   for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
583:     ex->Nt[ict]     = -1;
584:     ex->target[ict] = NULL;
585:     ex->size[ict]   = NULL;
586:     ex->cone[ict]   = NULL;
587:     ex->ornt[ict]   = NULL;
588:   }
589:   PetscCall(DMPlexTransformExtrudeSetUp_Point(ex, Nl));
590:   PetscCall(DMPlexTransformExtrudeSetUp_Segment(ex, Nl));
591:   PetscCall(DMPlexTransformExtrudeSetUp_Triangle(ex, Nl));
592:   PetscCall(DMPlexTransformExtrudeSetUp_Quadrilateral(ex, Nl));
593:   /* Layers positions */
594:   if (!ex->Nth) {
595:     if (ex->symmetric)
596:       for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers - ex->thickness / 2;
597:     else
598:       for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers;
599:   } else {
600:     ex->thickness   = 0.;
601:     ex->layerPos[0] = 0.;
602:     for (l = 0; l < ex->layers; ++l) {
603:       const PetscReal t   = ex->thicknesses[PetscMin(l, ex->Nth - 1)];
604:       ex->layerPos[l + 1] = ex->layerPos[l] + t;
605:       ex->thickness += t;
606:     }
607:     if (ex->symmetric)
608:       for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] -= ex->thickness / 2.;
609:   }
610:   PetscFunctionReturn(PETSC_SUCCESS);
611: }

613: static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr)
614: {
615:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
616:   PetscInt                 ct;

618:   PetscFunctionBegin;
619:   if (ex->target) {
620:     for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) PetscCall(PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct]));
621:   }
622:   PetscCall(PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt));
623:   PetscCall(PetscFree(ex->layerPos));
624:   PetscCall(DMDestroy(&ex->dmNormal));
625:   PetscCall(VecDestroy(&ex->vecNormal));
626:   PetscCall(PetscFree(ex));
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
631: {
632:   DMPlexTransform_Extrude *ex     = (DMPlexTransform_Extrude *)tr->data;
633:   DMLabel                  trType = tr->trType;
634:   PetscInt                 rt;

636:   PetscFunctionBeginHot;
637:   *rnew = r;
638:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
639:   if (!so) PetscFunctionReturn(PETSC_SUCCESS);
640:   if (trType) {
641:     PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
642:     if (rt >= 100) PetscFunctionReturn(PETSC_SUCCESS);
643:   }
644:   if (ex->useTensor) {
645:     switch (sct) {
646:     case DM_POLYTOPE_POINT:
647:       break;
648:     case DM_POLYTOPE_SEGMENT:
649:       switch (tct) {
650:       case DM_POLYTOPE_SEGMENT:
651:         break;
652:       case DM_POLYTOPE_SEG_PRISM_TENSOR:
653:         *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
654:         break;
655:       default:
656:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
657:       }
658:       break;
659:     // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
660:     case DM_POLYTOPE_TRIANGLE:
661:       break;
662:     case DM_POLYTOPE_QUADRILATERAL:
663:       break;
664:     default:
665:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
666:     }
667:   } else {
668:     switch (sct) {
669:     case DM_POLYTOPE_POINT:
670:       break;
671:     case DM_POLYTOPE_SEGMENT:
672:       switch (tct) {
673:       case DM_POLYTOPE_SEGMENT:
674:         break;
675:       case DM_POLYTOPE_QUADRILATERAL:
676:         *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0);
677:         break;
678:       default:
679:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
680:       }
681:       break;
682:     // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
683:     case DM_POLYTOPE_TRIANGLE:
684:       break;
685:     case DM_POLYTOPE_QUADRILATERAL:
686:       break;
687:     default:
688:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
689:     }
690:   }
691:   PetscFunctionReturn(PETSC_SUCCESS);
692: }

694: static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
695: {
696:   DMPlexTransform_Extrude *ex     = (DMPlexTransform_Extrude *)tr->data;
697:   DMLabel                  trType = tr->trType;
698:   PetscBool                ignore = PETSC_FALSE, identity = PETSC_FALSE;
699:   PetscInt                 val = 0;

701:   PetscFunctionBegin;
702:   if (trType) {
703:     PetscCall(DMLabelGetValue(trType, p, &val));
704:     identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE;
705:   } else {
706:     ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE;
707:   }
708:   if (rt) *rt = val;
709:   if (ignore) {
710:     /* Ignore cells that cannot be extruded */
711:     *Nt     = 0;
712:     *target = NULL;
713:     *size   = NULL;
714:     *cone   = NULL;
715:     *ornt   = NULL;
716:   } else if (identity) {
717:     PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
718:   } else {
719:     *Nt     = ex->Nt[source];
720:     *target = ex->target[source];
721:     *size   = ex->size[source];
722:     *cone   = ex->cone[source];
723:     *ornt   = ex->ornt[source];
724:   }
725:   PetscFunctionReturn(PETSC_SUCCESS);
726: }

728: /* Computes new vertex along normal */
729: static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
730: {
731:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
732:   DM                       dm;
733:   DMLabel                  active;
734:   PetscReal                ones2[2]  = {0., 1.};
735:   PetscReal                ones3[3]  = {0., 0., 1.};
736:   PetscReal                normal[3] = {0., 0., 0.};
737:   PetscReal                norm      = 0.;
738:   PetscInt                 dEx       = ex->cdimEx;
739:   PetscInt                 dim, cStart, cEnd;

741:   PetscFunctionBeginHot;
742:   PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]);
743:   PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
744:   PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv);
745:   PetscCheck(dE == ex->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %" PetscInt_FMT " != %" PetscInt_FMT " original dimension", dE, ex->cdim);
746:   PetscCheck(dEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", dEx);

748:   PetscCall(DMPlexTransformGetDM(tr, &dm));
749:   PetscCall(DMGetDimension(dm, &dim));
750:   PetscCall(DMPlexTransformGetActive(tr, &active));
751:   switch (ex->normalAlg) {
752:   case NORMAL_DEFAULT:
753:     switch (ex->cdimEx) {
754:     case 2:
755:       for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones2[d];
756:       break;
757:     case 3:
758:       for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones3[d];
759:       break;
760:     default:
761:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No default normal for dimension %" PetscInt_FMT, ex->cdimEx);
762:     }
763:     break;
764:   case NORMAL_INPUT:
765:     for (PetscInt d = 0; d < dEx; ++d) normal[d] = ex->normal[d];
766:     break;
767:   case NORMAL_COMPUTE: {
768:     PetscInt *star = NULL;
769:     PetscInt  starSize;

771:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
772:     PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star));
773:     for (PetscInt st = 0; st < starSize * 2; st += 2) {
774:       if ((star[st] >= cStart) && (star[st] < cEnd)) {
775:         PetscReal cnormal[3] = {0, 0, 0};

777:         PetscCall(DMPlexComputeCellGeometryFVM(dm, star[st], NULL, NULL, cnormal));
778:         for (PetscInt d = 0; d < dEx; ++d) normal[d] += cnormal[d];
779:       }
780:     }
781:     PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star));
782:   } break;
783:   case NORMAL_COMPUTE_BD: {
784:     const PetscScalar *a;
785:     PetscScalar       *vnormal;

787:     PetscCall(VecGetArrayRead(ex->vecNormal, &a));
788:     PetscCall(DMPlexPointLocalRead(ex->dmNormal, p, a, (void *)&vnormal));
789:     for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscRealPart(vnormal[d]);
790:     PetscCall(VecRestoreArrayRead(ex->vecNormal, &a));
791:   } break;
792:   default:
793:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
794:   }
795:   if (ex->normalFunc) {
796:     PetscScalar n[3];
797:     PetscReal   x[3], dot = 0.;

799:     for (PetscInt d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]);
800:     PetscCall((*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL));
801:     for (PetscInt d = 0; d < dEx; ++d) dot += PetscRealPart(n[d]) * normal[d];
802:     for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]);
803:   }
804:   for (PetscInt d = 0; d < dEx; ++d) norm += PetscSqr(normal[d]);
805:   for (PetscInt d = 0; d < dEx; ++d) normal[d] *= norm == 0.0 ? 1.0 : 1. / PetscSqrtReal(norm);
806:   for (PetscInt d = 0; d < dEx; ++d) out[d] = normal[d] * ex->layerPos[r];
807:   for (PetscInt d = 0; d < ex->cdim; ++d) out[d] += in[d];
808:   PetscFunctionReturn(PETSC_SUCCESS);
809: }

811: static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr)
812: {
813:   PetscFunctionBegin;
814:   tr->ops->view                  = DMPlexTransformView_Extrude;
815:   tr->ops->setfromoptions        = DMPlexTransformSetFromOptions_Extrude;
816:   tr->ops->setup                 = DMPlexTransformSetUp_Extrude;
817:   tr->ops->destroy               = DMPlexTransformDestroy_Extrude;
818:   tr->ops->setdimensions         = DMPlexTransformSetDimensions_Extrude;
819:   tr->ops->celltransform         = DMPlexTransformCellTransform_Extrude;
820:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude;
821:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinates_Extrude;
822:   PetscFunctionReturn(PETSC_SUCCESS);
823: }

825: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr)
826: {
827:   DMPlexTransform_Extrude *ex;
828:   DM                       dm;
829:   PetscInt                 dim;

831:   PetscFunctionBegin;
833:   PetscCall(PetscNew(&ex));
834:   tr->data      = ex;
835:   ex->thickness = 1.;
836:   ex->useTensor = PETSC_TRUE;
837:   ex->symmetric = PETSC_FALSE;
838:   ex->periodic  = PETSC_FALSE;
839:   ex->normalAlg = NORMAL_DEFAULT;
840:   ex->layerPos  = NULL;
841:   PetscCall(DMPlexTransformGetDM(tr, &dm));
842:   PetscCall(DMGetDimension(dm, &dim));
843:   PetscCall(DMGetCoordinateDim(dm, &ex->cdim));
844:   PetscCall(DMPlexTransformExtrudeSetLayers(tr, 1));
845:   PetscCall(DMPlexTransformInitialize_Extrude(tr));
846:   PetscFunctionReturn(PETSC_SUCCESS);
847: }

849: /*@
850:   DMPlexTransformExtrudeGetLayers - Get the number of extruded layers.

852:   Not Collective

854:   Input Parameter:
855: . tr - The `DMPlexTransform`

857:   Output Parameter:
858: . layers - The number of layers

860:   Level: intermediate

862: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetLayers()`
863: @*/
864: PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers)
865: {
866:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

868:   PetscFunctionBegin;
870:   PetscAssertPointer(layers, 2);
871:   *layers = ex->layers;
872:   PetscFunctionReturn(PETSC_SUCCESS);
873: }

875: /*@
876:   DMPlexTransformExtrudeSetLayers - Set the number of extruded layers.

878:   Not Collective

880:   Input Parameters:
881: + tr     - The `DMPlexTransform`
882: - layers - The number of layers

884:   Level: intermediate

886: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetLayers()`
887: @*/
888: PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers)
889: {
890:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

892:   PetscFunctionBegin;
894:   ex->layers = layers;
895:   PetscCall(PetscFree(ex->layerPos));
896:   PetscCall(PetscCalloc1(ex->layers + 1, &ex->layerPos));
897:   PetscFunctionReturn(PETSC_SUCCESS);
898: }

900: /*@
901:   DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers

903:   Not Collective

905:   Input Parameter:
906: . tr - The `DMPlexTransform`

908:   Output Parameter:
909: . thickness - The total thickness of the layers

911:   Level: intermediate

913: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`
914: @*/
915: PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness)
916: {
917:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

919:   PetscFunctionBegin;
921:   PetscAssertPointer(thickness, 2);
922:   *thickness = ex->thickness;
923:   PetscFunctionReturn(PETSC_SUCCESS);
924: }

926: /*@
927:   DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers

929:   Not Collective

931:   Input Parameters:
932: + tr        - The `DMPlexTransform`
933: - thickness - The total thickness of the layers

935:   Level: intermediate

937: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetThickness()`
938: @*/
939: PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness)
940: {
941:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

943:   PetscFunctionBegin;
945:   PetscCheck(thickness > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double)thickness);
946:   ex->thickness = thickness;
947:   PetscFunctionReturn(PETSC_SUCCESS);
948: }

950: /*@
951:   DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells

953:   Not Collective

955:   Input Parameter:
956: . tr - The `DMPlexTransform`

958:   Output Parameter:
959: . useTensor - The flag to use tensor cells

961:   Note:
962:   This flag determines the orientation behavior of the created points.

964:   For example, if tensor is `PETSC_TRUE`, then
965: .vb
966:   DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
967:   DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
968:   DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
969:   DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
970: .ve

972:   Level: intermediate

974: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetTensor()`
975: @*/
976: PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
977: {
978:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

980:   PetscFunctionBegin;
982:   PetscAssertPointer(useTensor, 2);
983:   *useTensor = ex->useTensor;
984:   PetscFunctionReturn(PETSC_SUCCESS);
985: }

987: /*@
988:   DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells

990:   Not Collective

992:   Input Parameters:
993: + tr        - The `DMPlexTransform`
994: - useTensor - The flag for tensor cells

996:   Note:
997:   This flag determines the orientation behavior of the created points
998:   For example, if tensor is `PETSC_TRUE`, then
999: .vb
1000:   DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1001:   DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1002:   DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1003:   DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1004: .ve

1006:   Level: intermediate

1008: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetTensor()`
1009: @*/
1010: PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
1011: {
1012:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

1014:   PetscFunctionBegin;
1016:   ex->useTensor = useTensor;
1017:   PetscFunctionReturn(PETSC_SUCCESS);
1018: }

1020: /*@
1021:   DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface

1023:   Not Collective

1025:   Input Parameter:
1026: . tr - The `DMPlexTransform`

1028:   Output Parameter:
1029: . symmetric - The flag to extrude symmetrically

1031:   Level: intermediate

1033: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetSymmetric()`, `DMPlexTransformExtrudeGetPeriodic()`
1034: @*/
1035: PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric)
1036: {
1037:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

1039:   PetscFunctionBegin;
1041:   PetscAssertPointer(symmetric, 2);
1042:   *symmetric = ex->symmetric;
1043:   PetscFunctionReturn(PETSC_SUCCESS);
1044: }

1046: /*@
1047:   DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface

1049:   Not Collective

1051:   Input Parameters:
1052: + tr        - The `DMPlexTransform`
1053: - symmetric - The flag to extrude symmetrically

1055:   Level: intermediate

1057: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetSymmetric()`, `DMPlexTransformExtrudeSetPeriodic()`
1058: @*/
1059: PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric)
1060: {
1061:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

1063:   PetscFunctionBegin;
1065:   ex->symmetric = symmetric;
1066:   PetscFunctionReturn(PETSC_SUCCESS);
1067: }

1069: /*@
1070:   DMPlexTransformExtrudeGetPeriodic - Get the flag to extrude periodically from the initial surface

1072:   Not Collective

1074:   Input Parameter:
1075: . tr - The `DMPlexTransform`

1077:   Output Parameter:
1078: . periodic - The flag to extrude periodically

1080:   Level: intermediate

1082: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetPeriodic()`, `DMPlexTransformExtrudeGetSymmetric()`
1083: @*/
1084: PetscErrorCode DMPlexTransformExtrudeGetPeriodic(DMPlexTransform tr, PetscBool *periodic)
1085: {
1086:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

1088:   PetscFunctionBegin;
1090:   PetscAssertPointer(periodic, 2);
1091:   *periodic = ex->periodic;
1092:   PetscFunctionReturn(PETSC_SUCCESS);
1093: }

1095: /*@
1096:   DMPlexTransformExtrudeSetPeriodic - Set the flag to extrude periodically from the initial surface

1098:   Not Collective

1100:   Input Parameters:
1101: + tr       - The `DMPlexTransform`
1102: - periodic - The flag to extrude periodically

1104:   Level: intermediate

1106: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetPeriodic()`, `DMPlexTransformExtrudeSetSymmetric()`
1107: @*/
1108: PetscErrorCode DMPlexTransformExtrudeSetPeriodic(DMPlexTransform tr, PetscBool periodic)
1109: {
1110:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

1112:   PetscFunctionBegin;
1114:   ex->periodic = periodic;
1115:   PetscFunctionReturn(PETSC_SUCCESS);
1116: }

1118: /*@
1119:   DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector

1121:   Not Collective

1123:   Input Parameter:
1124: . tr - The `DMPlexTransform`

1126:   Output Parameter:
1127: . normal - The extrusion direction

1129:   Note:
1130:   The user passes in an array, which is filled by the library.

1132:   Level: intermediate

1134: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetNormal()`
1135: @*/
1136: PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[])
1137: {
1138:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1139:   PetscInt                 d;

1141:   PetscFunctionBegin;
1143:   if (ex->normalAlg == NORMAL_INPUT) {
1144:     for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d];
1145:   } else {
1146:     for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.;
1147:   }
1148:   PetscFunctionReturn(PETSC_SUCCESS);
1149: }

1151: /*@
1152:   DMPlexTransformExtrudeSetNormal - Set the extrusion normal

1154:   Not Collective

1156:   Input Parameters:
1157: + tr     - The `DMPlexTransform`
1158: - normal - The extrusion direction

1160:   Level: intermediate

1162: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`
1163: @*/
1164: PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[])
1165: {
1166:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1167:   PetscInt                 d;

1169:   PetscFunctionBegin;
1171:   ex->normalAlg = NORMAL_INPUT;
1172:   for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d];
1173:   PetscFunctionReturn(PETSC_SUCCESS);
1174: }

1176: /*@C
1177:   DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal

1179:   Not Collective

1181:   Input Parameters:
1182: + tr         - The `DMPlexTransform`
1183: - normalFunc - A function determining the extrusion direction, see `PetscSimplePointFn` for the calling sequence

1185:   Level: intermediate

1187: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`, `PetscSimplePointFn`
1188: @*/
1189: PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFn *normalFunc)
1190: {
1191:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

1193:   PetscFunctionBegin;
1195:   ex->normalFunc = normalFunc;
1196:   PetscFunctionReturn(PETSC_SUCCESS);
1197: }

1199: /*@
1200:   DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer

1202:   Not Collective

1204:   Input Parameters:
1205: + tr          - The `DMPlexTransform`
1206: . Nth         - The number of thicknesses
1207: - thicknesses - The array of thicknesses

1209:   Level: intermediate

1211: .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeGetThickness()`
1212: @*/
1213: PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[])
1214: {
1215:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1216:   PetscInt                 t;

1218:   PetscFunctionBegin;
1220:   PetscCheck(Nth > 0, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %" PetscInt_FMT " must be positive", Nth);
1221:   ex->Nth = PetscMin(Nth, ex->layers);
1222:   PetscCall(PetscFree(ex->thicknesses));
1223:   PetscCall(PetscMalloc1(ex->Nth, &ex->thicknesses));
1224:   for (t = 0; t < ex->Nth; ++t) {
1225:     PetscCheck(thicknesses[t] > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %" PetscInt_FMT " must be positive", (double)thicknesses[t], t);
1226:     ex->thicknesses[t] = thicknesses[t];
1227:   }
1228:   PetscFunctionReturn(PETSC_SUCCESS);
1229: }