Actual source code: plextrextrude.c

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

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

 10:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
 11:   if (isascii) {
 12:     const char *name;

 14:     PetscObjectGetName((PetscObject)tr, &name);
 15:     PetscViewerASCIIPrintf(viewer, "Extrusion transformation %s\n", name ? name : "");
 16:     PetscViewerASCIIPrintf(viewer, "  number of layers: %" PetscInt_FMT "\n", ex->layers);
 17:     PetscViewerASCIIPrintf(viewer, "  create tensor cells: %s\n", ex->useTensor ? "YES" : "NO");
 18:   } else {
 19:     SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
 20:   }
 21:   return 0;
 22: }

 24: static PetscErrorCode DMPlexTransformSetFromOptions_Extrude(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
 25: {
 26:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
 27:   PetscReal                th, normal[3], *thicknesses;
 28:   PetscInt                 nl, Nc;
 29:   PetscBool                tensor, sym, flg;
 30:   char                     funcname[PETSC_MAX_PATH_LEN];

 32:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Extrusion Options");
 33:   PetscOptionsBoundedInt("-dm_plex_transform_extrude_layers", "Number of layers to extrude", "", ex->layers, &nl, &flg, 1);
 34:   if (flg) DMPlexTransformExtrudeSetLayers(tr, nl);
 35:   PetscOptionsReal("-dm_plex_transform_extrude_thickness", "Total thickness of extruded layers", "", ex->thickness, &th, &flg);
 36:   if (flg) DMPlexTransformExtrudeSetThickness(tr, th);
 37:   PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg);
 38:   if (flg) DMPlexTransformExtrudeSetTensor(tr, tensor);
 39:   PetscOptionsBool("-dm_plex_transform_extrude_symmetric", "Extrude layers symmetrically about the surface", "", ex->symmetric, &sym, &flg);
 40:   if (flg) DMPlexTransformExtrudeSetSymmetric(tr, sym);
 41:   Nc = 3;
 42:   PetscOptionsRealArray("-dm_plex_transform_extrude_normal", "Input normal vector for extrusion", "DMPlexTransformExtrudeSetNormal", normal, &Nc, &flg);
 43:   if (flg) {
 44:     // Extrusion dimension might not yet be determined
 46:     DMPlexTransformExtrudeSetNormal(tr, normal);
 47:   }
 48:   PetscOptionsString("-dm_plex_transform_extrude_normal_function", "Function to determine normal vector", "DMPlexTransformExtrudeSetNormalFunction", funcname, funcname, sizeof(funcname), &flg);
 49:   if (flg) {
 50:     PetscSimplePointFunc normalFunc;

 52:     PetscDLSym(NULL, funcname, (void **)&normalFunc);
 53:     DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc);
 54:   }
 55:   nl = ex->layers;
 56:   PetscMalloc1(nl, &thicknesses);
 57:   PetscOptionsRealArray("-dm_plex_transform_extrude_thicknesses", "Thickness of each individual extruded layer", "", thicknesses, &nl, &flg);
 58:   if (flg) {
 60:     DMPlexTransformExtrudeSetThicknesses(tr, nl, thicknesses);
 61:   }
 62:   PetscFree(thicknesses);
 63:   PetscOptionsHeadEnd();
 64:   return 0;
 65: }

 67: /* Determine the implicit dimension pre-extrusion (either the implicit dimension of the DM or of a point in the active set for the transform).
 68:    If that dimension is the same as the current coordinate dimension (ex->dim), the extruded mesh will have a coordinate dimension one greater;
 69:    Otherwise the coordinate dimension will be kept. */
 70: static PetscErrorCode DMPlexTransformExtrudeComputeExtrusionDim(DMPlexTransform tr)
 71: {
 72:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
 73:   DM                       dm;
 74:   DMLabel                  active;
 75:   PetscInt                 dim, dimExtPoint, dimExtPointG;

 77:   DMPlexTransformGetDM(tr, &dm);
 78:   DMGetDimension(dm, &dim);
 79:   DMPlexTransformGetActive(tr, &active);
 80:   if (active) {
 81:     IS              valueIS, pointIS;
 82:     const PetscInt *values, *points;
 83:     DMPolytopeType  ct;
 84:     PetscInt        Nv, Np;

 86:     dimExtPoint = 0;
 87:     DMLabelGetValueIS(active, &valueIS);
 88:     ISGetLocalSize(valueIS, &Nv);
 89:     ISGetIndices(valueIS, &values);
 90:     for (PetscInt v = 0; v < Nv; ++v) {
 91:       DMLabelGetStratumIS(active, values[v], &pointIS);
 92:       ISGetLocalSize(pointIS, &Np);
 93:       ISGetIndices(pointIS, &points);
 94:       for (PetscInt p = 0; p < Np; ++p) {
 95:         DMPlexGetCellType(dm, points[p], &ct);
 96:         dimExtPoint = PetscMax(dimExtPoint, DMPolytopeTypeGetDim(ct));
 97:       }
 98:       ISRestoreIndices(pointIS, &points);
 99:       ISDestroy(&pointIS);
100:     }
101:     ISRestoreIndices(valueIS, &values);
102:     ISDestroy(&valueIS);
103:   } else dimExtPoint = dim;
104:   MPI_Allreduce(&dimExtPoint, &dimExtPointG, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)tr));
105:   ex->dimEx  = PetscMax(dim, dimExtPointG + 1);
106:   ex->cdimEx = ex->cdim == dimExtPointG ? ex->cdim + 1 : ex->cdim;
109:   return 0;
110: }

112: static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm)
113: {
114:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
115:   PetscInt                 dim;

117:   DMGetDimension(dm, &dim);
118:   DMSetDimension(tdm, ex->dimEx);
119:   DMSetCoordinateDim(tdm, ex->cdimEx);
120:   return 0;
121: }

123: /*
124:   The refine types for extrusion are:

126:   ct:       For any normally extruded point
127:   ct + 100: For any point which should just return itself
128: */
129: static PetscErrorCode DMPlexTransformSetUp_Extrude(DMPlexTransform tr)
130: {
131:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
132:   DM                       dm;
133:   DMLabel                  active;
134:   DMPolytopeType           ct;
135:   PetscInt                 Nl = ex->layers, l, i, ict, Nc, No, coff, ooff;

137:   DMPlexTransformExtrudeComputeExtrusionDim(tr);
138:   DMPlexTransformGetDM(tr, &dm);
139:   DMPlexTransformGetActive(tr, &active);
140:   if (active) {
141:     DMLabel  celltype;
142:     PetscInt pStart, pEnd, p;

144:     DMPlexGetCellTypeLabel(dm, &celltype);
145:     DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType);
146:     DMPlexGetChart(dm, &pStart, &pEnd);
147:     for (p = pStart; p < pEnd; ++p) {
148:       PetscInt ct, val;

150:       DMLabelGetValue(celltype, p, &ct);
151:       DMLabelGetValue(active, p, &val);
152:       if (val < 0) {
153:         DMLabelSetValue(tr->trType, p, ct + 100);
154:       } else {
155:         DMLabelSetValue(tr->trType, p, ct);
156:       }
157:     }
158:   }
159:   PetscMalloc5(DM_NUM_POLYTOPES, &ex->Nt, DM_NUM_POLYTOPES, &ex->target, DM_NUM_POLYTOPES, &ex->size, DM_NUM_POLYTOPES, &ex->cone, DM_NUM_POLYTOPES, &ex->ornt);
160:   for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
161:     ex->Nt[ict]     = -1;
162:     ex->target[ict] = NULL;
163:     ex->size[ict]   = NULL;
164:     ex->cone[ict]   = NULL;
165:     ex->ornt[ict]   = NULL;
166:   }
167:   /* DM_POLYTOPE_POINT produces Nl+1 points and Nl segments/tensor segments */
168:   ct         = DM_POLYTOPE_POINT;
169:   ex->Nt[ct] = 2;
170:   Nc         = 6 * Nl;
171:   No         = 2 * Nl;
172:   PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);
173:   ex->target[ct][0] = DM_POLYTOPE_POINT;
174:   ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
175:   ex->size[ct][0]   = Nl + 1;
176:   ex->size[ct][1]   = Nl;
177:   /*   cones for segments/tensor segments */
178:   for (i = 0; i < Nl; ++i) {
179:     ex->cone[ct][6 * i + 0] = DM_POLYTOPE_POINT;
180:     ex->cone[ct][6 * i + 1] = 0;
181:     ex->cone[ct][6 * i + 2] = i;
182:     ex->cone[ct][6 * i + 3] = DM_POLYTOPE_POINT;
183:     ex->cone[ct][6 * i + 4] = 0;
184:     ex->cone[ct][6 * i + 5] = i + 1;
185:   }
186:   for (i = 0; i < No; ++i) ex->ornt[ct][i] = 0;
187:   /* DM_POLYTOPE_SEGMENT produces Nl+1 segments and Nl quads/tensor quads */
188:   ct         = DM_POLYTOPE_SEGMENT;
189:   ex->Nt[ct] = 2;
190:   Nc         = 8 * (Nl + 1) + 14 * Nl;
191:   No         = 2 * (Nl + 1) + 4 * Nl;
192:   PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);
193:   ex->target[ct][0] = DM_POLYTOPE_SEGMENT;
194:   ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
195:   ex->size[ct][0]   = Nl + 1;
196:   ex->size[ct][1]   = Nl;
197:   /*   cones for segments */
198:   for (i = 0; i < Nl + 1; ++i) {
199:     ex->cone[ct][8 * i + 0] = DM_POLYTOPE_POINT;
200:     ex->cone[ct][8 * i + 1] = 1;
201:     ex->cone[ct][8 * i + 2] = 0;
202:     ex->cone[ct][8 * i + 3] = i;
203:     ex->cone[ct][8 * i + 4] = DM_POLYTOPE_POINT;
204:     ex->cone[ct][8 * i + 5] = 1;
205:     ex->cone[ct][8 * i + 6] = 1;
206:     ex->cone[ct][8 * i + 7] = i;
207:   }
208:   for (i = 0; i < 2 * (Nl + 1); ++i) ex->ornt[ct][i] = 0;
209:   /*   cones for quads/tensor quads */
210:   coff = 8 * (Nl + 1);
211:   ooff = 2 * (Nl + 1);
212:   for (i = 0; i < Nl; ++i) {
213:     if (ex->useTensor) {
214:       ex->cone[ct][coff + 14 * i + 0]  = DM_POLYTOPE_SEGMENT;
215:       ex->cone[ct][coff + 14 * i + 1]  = 0;
216:       ex->cone[ct][coff + 14 * i + 2]  = i;
217:       ex->cone[ct][coff + 14 * i + 3]  = DM_POLYTOPE_SEGMENT;
218:       ex->cone[ct][coff + 14 * i + 4]  = 0;
219:       ex->cone[ct][coff + 14 * i + 5]  = i + 1;
220:       ex->cone[ct][coff + 14 * i + 6]  = DM_POLYTOPE_POINT_PRISM_TENSOR;
221:       ex->cone[ct][coff + 14 * i + 7]  = 1;
222:       ex->cone[ct][coff + 14 * i + 8]  = 0;
223:       ex->cone[ct][coff + 14 * i + 9]  = i;
224:       ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
225:       ex->cone[ct][coff + 14 * i + 11] = 1;
226:       ex->cone[ct][coff + 14 * i + 12] = 1;
227:       ex->cone[ct][coff + 14 * i + 13] = i;
228:       ex->ornt[ct][ooff + 4 * i + 0]   = 0;
229:       ex->ornt[ct][ooff + 4 * i + 1]   = 0;
230:       ex->ornt[ct][ooff + 4 * i + 2]   = 0;
231:       ex->ornt[ct][ooff + 4 * i + 3]   = 0;
232:     } else {
233:       ex->cone[ct][coff + 14 * i + 0]  = DM_POLYTOPE_SEGMENT;
234:       ex->cone[ct][coff + 14 * i + 1]  = 0;
235:       ex->cone[ct][coff + 14 * i + 2]  = i;
236:       ex->cone[ct][coff + 14 * i + 3]  = DM_POLYTOPE_SEGMENT;
237:       ex->cone[ct][coff + 14 * i + 4]  = 1;
238:       ex->cone[ct][coff + 14 * i + 5]  = 1;
239:       ex->cone[ct][coff + 14 * i + 6]  = i;
240:       ex->cone[ct][coff + 14 * i + 7]  = DM_POLYTOPE_SEGMENT;
241:       ex->cone[ct][coff + 14 * i + 8]  = 0;
242:       ex->cone[ct][coff + 14 * i + 9]  = i + 1;
243:       ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_SEGMENT;
244:       ex->cone[ct][coff + 14 * i + 11] = 1;
245:       ex->cone[ct][coff + 14 * i + 12] = 0;
246:       ex->cone[ct][coff + 14 * i + 13] = i;
247:       ex->ornt[ct][ooff + 4 * i + 0]   = 0;
248:       ex->ornt[ct][ooff + 4 * i + 1]   = 0;
249:       ex->ornt[ct][ooff + 4 * i + 2]   = -1;
250:       ex->ornt[ct][ooff + 4 * i + 3]   = -1;
251:     }
252:   }
253:   /* DM_POLYTOPE_TRIANGLE produces Nl+1 triangles and Nl triangular prisms/tensor triangular prisms */
254:   ct         = DM_POLYTOPE_TRIANGLE;
255:   ex->Nt[ct] = 2;
256:   Nc         = 12 * (Nl + 1) + 18 * Nl;
257:   No         = 3 * (Nl + 1) + 5 * Nl;
258:   PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);
259:   ex->target[ct][0] = DM_POLYTOPE_TRIANGLE;
260:   ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM;
261:   ex->size[ct][0]   = Nl + 1;
262:   ex->size[ct][1]   = Nl;
263:   /*   cones for triangles */
264:   for (i = 0; i < Nl + 1; ++i) {
265:     ex->cone[ct][12 * i + 0]  = DM_POLYTOPE_SEGMENT;
266:     ex->cone[ct][12 * i + 1]  = 1;
267:     ex->cone[ct][12 * i + 2]  = 0;
268:     ex->cone[ct][12 * i + 3]  = i;
269:     ex->cone[ct][12 * i + 4]  = DM_POLYTOPE_SEGMENT;
270:     ex->cone[ct][12 * i + 5]  = 1;
271:     ex->cone[ct][12 * i + 6]  = 1;
272:     ex->cone[ct][12 * i + 7]  = i;
273:     ex->cone[ct][12 * i + 8]  = DM_POLYTOPE_SEGMENT;
274:     ex->cone[ct][12 * i + 9]  = 1;
275:     ex->cone[ct][12 * i + 10] = 2;
276:     ex->cone[ct][12 * i + 11] = i;
277:   }
278:   for (i = 0; i < 3 * (Nl + 1); ++i) ex->ornt[ct][i] = 0;
279:   /*   cones for triangular prisms/tensor triangular prisms */
280:   coff = 12 * (Nl + 1);
281:   ooff = 3 * (Nl + 1);
282:   for (i = 0; i < Nl; ++i) {
283:     if (ex->useTensor) {
284:       ex->cone[ct][coff + 18 * i + 0]  = DM_POLYTOPE_TRIANGLE;
285:       ex->cone[ct][coff + 18 * i + 1]  = 0;
286:       ex->cone[ct][coff + 18 * i + 2]  = i;
287:       ex->cone[ct][coff + 18 * i + 3]  = DM_POLYTOPE_TRIANGLE;
288:       ex->cone[ct][coff + 18 * i + 4]  = 0;
289:       ex->cone[ct][coff + 18 * i + 5]  = i + 1;
290:       ex->cone[ct][coff + 18 * i + 6]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
291:       ex->cone[ct][coff + 18 * i + 7]  = 1;
292:       ex->cone[ct][coff + 18 * i + 8]  = 0;
293:       ex->cone[ct][coff + 18 * i + 9]  = i;
294:       ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
295:       ex->cone[ct][coff + 18 * i + 11] = 1;
296:       ex->cone[ct][coff + 18 * i + 12] = 1;
297:       ex->cone[ct][coff + 18 * i + 13] = i;
298:       ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
299:       ex->cone[ct][coff + 18 * i + 15] = 1;
300:       ex->cone[ct][coff + 18 * i + 16] = 2;
301:       ex->cone[ct][coff + 18 * i + 17] = i;
302:       ex->ornt[ct][ooff + 5 * i + 0]   = 0;
303:       ex->ornt[ct][ooff + 5 * i + 1]   = 0;
304:       ex->ornt[ct][ooff + 5 * i + 2]   = 0;
305:       ex->ornt[ct][ooff + 5 * i + 3]   = 0;
306:       ex->ornt[ct][ooff + 5 * i + 4]   = 0;
307:     } else {
308:       ex->cone[ct][coff + 18 * i + 0]  = DM_POLYTOPE_TRIANGLE;
309:       ex->cone[ct][coff + 18 * i + 1]  = 0;
310:       ex->cone[ct][coff + 18 * i + 2]  = i;
311:       ex->cone[ct][coff + 18 * i + 3]  = DM_POLYTOPE_TRIANGLE;
312:       ex->cone[ct][coff + 18 * i + 4]  = 0;
313:       ex->cone[ct][coff + 18 * i + 5]  = i + 1;
314:       ex->cone[ct][coff + 18 * i + 6]  = DM_POLYTOPE_QUADRILATERAL;
315:       ex->cone[ct][coff + 18 * i + 7]  = 1;
316:       ex->cone[ct][coff + 18 * i + 8]  = 0;
317:       ex->cone[ct][coff + 18 * i + 9]  = i;
318:       ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
319:       ex->cone[ct][coff + 18 * i + 11] = 1;
320:       ex->cone[ct][coff + 18 * i + 12] = 1;
321:       ex->cone[ct][coff + 18 * i + 13] = i;
322:       ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
323:       ex->cone[ct][coff + 18 * i + 15] = 1;
324:       ex->cone[ct][coff + 18 * i + 16] = 2;
325:       ex->cone[ct][coff + 18 * i + 17] = i;
326:       ex->ornt[ct][ooff + 5 * i + 0]   = -2;
327:       ex->ornt[ct][ooff + 5 * i + 1]   = 0;
328:       ex->ornt[ct][ooff + 5 * i + 2]   = 0;
329:       ex->ornt[ct][ooff + 5 * i + 3]   = 0;
330:       ex->ornt[ct][ooff + 5 * i + 4]   = 0;
331:     }
332:   }
333:   /* DM_POLYTOPE_QUADRILATERAL produces Nl+1 quads and Nl hexes/tensor hexes */
334:   ct         = DM_POLYTOPE_QUADRILATERAL;
335:   ex->Nt[ct] = 2;
336:   Nc         = 16 * (Nl + 1) + 22 * Nl;
337:   No         = 4 * (Nl + 1) + 6 * Nl;
338:   PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);
339:   ex->target[ct][0] = DM_POLYTOPE_QUADRILATERAL;
340:   ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON;
341:   ex->size[ct][0]   = Nl + 1;
342:   ex->size[ct][1]   = Nl;
343:   /*   cones for quads */
344:   for (i = 0; i < Nl + 1; ++i) {
345:     ex->cone[ct][16 * i + 0]  = DM_POLYTOPE_SEGMENT;
346:     ex->cone[ct][16 * i + 1]  = 1;
347:     ex->cone[ct][16 * i + 2]  = 0;
348:     ex->cone[ct][16 * i + 3]  = i;
349:     ex->cone[ct][16 * i + 4]  = DM_POLYTOPE_SEGMENT;
350:     ex->cone[ct][16 * i + 5]  = 1;
351:     ex->cone[ct][16 * i + 6]  = 1;
352:     ex->cone[ct][16 * i + 7]  = i;
353:     ex->cone[ct][16 * i + 8]  = DM_POLYTOPE_SEGMENT;
354:     ex->cone[ct][16 * i + 9]  = 1;
355:     ex->cone[ct][16 * i + 10] = 2;
356:     ex->cone[ct][16 * i + 11] = i;
357:     ex->cone[ct][16 * i + 12] = DM_POLYTOPE_SEGMENT;
358:     ex->cone[ct][16 * i + 13] = 1;
359:     ex->cone[ct][16 * i + 14] = 3;
360:     ex->cone[ct][16 * i + 15] = i;
361:   }
362:   for (i = 0; i < 4 * (Nl + 1); ++i) ex->ornt[ct][i] = 0;
363:   /*   cones for hexes/tensor hexes */
364:   coff = 16 * (Nl + 1);
365:   ooff = 4 * (Nl + 1);
366:   for (i = 0; i < Nl; ++i) {
367:     if (ex->useTensor) {
368:       ex->cone[ct][coff + 22 * i + 0]  = DM_POLYTOPE_QUADRILATERAL;
369:       ex->cone[ct][coff + 22 * i + 1]  = 0;
370:       ex->cone[ct][coff + 22 * i + 2]  = i;
371:       ex->cone[ct][coff + 22 * i + 3]  = DM_POLYTOPE_QUADRILATERAL;
372:       ex->cone[ct][coff + 22 * i + 4]  = 0;
373:       ex->cone[ct][coff + 22 * i + 5]  = i + 1;
374:       ex->cone[ct][coff + 22 * i + 6]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
375:       ex->cone[ct][coff + 22 * i + 7]  = 1;
376:       ex->cone[ct][coff + 22 * i + 8]  = 0;
377:       ex->cone[ct][coff + 22 * i + 9]  = i;
378:       ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
379:       ex->cone[ct][coff + 22 * i + 11] = 1;
380:       ex->cone[ct][coff + 22 * i + 12] = 1;
381:       ex->cone[ct][coff + 22 * i + 13] = i;
382:       ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
383:       ex->cone[ct][coff + 22 * i + 15] = 1;
384:       ex->cone[ct][coff + 22 * i + 16] = 2;
385:       ex->cone[ct][coff + 22 * i + 17] = i;
386:       ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
387:       ex->cone[ct][coff + 22 * i + 19] = 1;
388:       ex->cone[ct][coff + 22 * i + 20] = 3;
389:       ex->cone[ct][coff + 22 * i + 21] = i;
390:       ex->ornt[ct][ooff + 6 * i + 0]   = 0;
391:       ex->ornt[ct][ooff + 6 * i + 1]   = 0;
392:       ex->ornt[ct][ooff + 6 * i + 2]   = 0;
393:       ex->ornt[ct][ooff + 6 * i + 3]   = 0;
394:       ex->ornt[ct][ooff + 6 * i + 4]   = 0;
395:       ex->ornt[ct][ooff + 6 * i + 5]   = 0;
396:     } else {
397:       ex->cone[ct][coff + 22 * i + 0]  = DM_POLYTOPE_QUADRILATERAL;
398:       ex->cone[ct][coff + 22 * i + 1]  = 0;
399:       ex->cone[ct][coff + 22 * i + 2]  = i;
400:       ex->cone[ct][coff + 22 * i + 3]  = DM_POLYTOPE_QUADRILATERAL;
401:       ex->cone[ct][coff + 22 * i + 4]  = 0;
402:       ex->cone[ct][coff + 22 * i + 5]  = i + 1;
403:       ex->cone[ct][coff + 22 * i + 6]  = DM_POLYTOPE_QUADRILATERAL;
404:       ex->cone[ct][coff + 22 * i + 7]  = 1;
405:       ex->cone[ct][coff + 22 * i + 8]  = 0;
406:       ex->cone[ct][coff + 22 * i + 9]  = i;
407:       ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_QUADRILATERAL;
408:       ex->cone[ct][coff + 22 * i + 11] = 1;
409:       ex->cone[ct][coff + 22 * i + 12] = 2;
410:       ex->cone[ct][coff + 22 * i + 13] = i;
411:       ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_QUADRILATERAL;
412:       ex->cone[ct][coff + 22 * i + 15] = 1;
413:       ex->cone[ct][coff + 22 * i + 16] = 1;
414:       ex->cone[ct][coff + 22 * i + 17] = i;
415:       ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_QUADRILATERAL;
416:       ex->cone[ct][coff + 22 * i + 19] = 1;
417:       ex->cone[ct][coff + 22 * i + 20] = 3;
418:       ex->cone[ct][coff + 22 * i + 21] = i;
419:       ex->ornt[ct][ooff + 6 * i + 0]   = -2;
420:       ex->ornt[ct][ooff + 6 * i + 1]   = 0;
421:       ex->ornt[ct][ooff + 6 * i + 2]   = 0;
422:       ex->ornt[ct][ooff + 6 * i + 3]   = 0;
423:       ex->ornt[ct][ooff + 6 * i + 4]   = 0;
424:       ex->ornt[ct][ooff + 6 * i + 5]   = 1;
425:     }
426:   }
427:   /* Layers positions */
428:   if (!ex->Nth) {
429:     if (ex->symmetric)
430:       for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers - ex->thickness / 2;
431:     else
432:       for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers;
433:   } else {
434:     ex->thickness   = 0.;
435:     ex->layerPos[0] = 0.;
436:     for (l = 0; l < ex->layers; ++l) {
437:       const PetscReal t   = ex->thicknesses[PetscMin(l, ex->Nth - 1)];
438:       ex->layerPos[l + 1] = ex->layerPos[l] + t;
439:       ex->thickness += t;
440:     }
441:     if (ex->symmetric)
442:       for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] -= ex->thickness / 2.;
443:   }
444:   return 0;
445: }

447: static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr)
448: {
449:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
450:   PetscInt                 ct;

452:   for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct]);
453:   PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt);
454:   PetscFree(ex->layerPos);
455:   PetscFree(ex);
456:   return 0;
457: }

459: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
460: {
461:   DMPlexTransform_Extrude *ex     = (DMPlexTransform_Extrude *)tr->data;
462:   DMLabel                  trType = tr->trType;
463:   PetscInt                 rt;

466:   *rnew = r;
467:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
468:   if (!so) return 0;
469:   if (trType) {
470:     DMLabelGetValue(tr->trType, sp, &rt);
471:     if (rt >= 100) return 0;
472:   }
473:   if (ex->useTensor) {
474:     switch (sct) {
475:     case DM_POLYTOPE_POINT:
476:       break;
477:     case DM_POLYTOPE_SEGMENT:
478:       switch (tct) {
479:       case DM_POLYTOPE_SEGMENT:
480:         break;
481:       case DM_POLYTOPE_SEG_PRISM_TENSOR:
482:         *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
483:         break;
484:       default:
485:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
486:       }
487:       break;
488:     // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
489:     case DM_POLYTOPE_TRIANGLE:
490:       break;
491:     case DM_POLYTOPE_QUADRILATERAL:
492:       break;
493:     default:
494:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
495:     }
496:   } else {
497:     switch (sct) {
498:     case DM_POLYTOPE_POINT:
499:       break;
500:     case DM_POLYTOPE_SEGMENT:
501:       switch (tct) {
502:       case DM_POLYTOPE_SEGMENT:
503:         break;
504:       case DM_POLYTOPE_QUADRILATERAL:
505:         *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0);
506:         break;
507:       default:
508:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
509:       }
510:       break;
511:     // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
512:     case DM_POLYTOPE_TRIANGLE:
513:       break;
514:     case DM_POLYTOPE_QUADRILATERAL:
515:       break;
516:     default:
517:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
518:     }
519:   }
520:   return 0;
521: }

523: static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
524: {
525:   DMPlexTransform_Extrude *ex     = (DMPlexTransform_Extrude *)tr->data;
526:   DMLabel                  trType = tr->trType;
527:   PetscBool                ignore = PETSC_FALSE, identity = PETSC_FALSE;
528:   PetscInt                 val = 0;

530:   if (trType) {
531:     DMLabelGetValue(trType, p, &val);
532:     identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE;
533:   } else {
534:     ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE;
535:   }
536:   if (rt) *rt = val;
537:   if (ignore) {
538:     /* Ignore cells that cannot be extruded */
539:     *Nt     = 0;
540:     *target = NULL;
541:     *size   = NULL;
542:     *cone   = NULL;
543:     *ornt   = NULL;
544:   } else if (identity) {
545:     DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
546:   } else {
547:     *Nt     = ex->Nt[source];
548:     *target = ex->target[source];
549:     *size   = ex->size[source];
550:     *cone   = ex->cone[source];
551:     *ornt   = ex->ornt[source];
552:   }
553:   return 0;
554: }

556: /* Computes new vertex along normal */
557: static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
558: {
559:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
560:   DM                       dm;
561:   DMLabel                  active;
562:   PetscReal                ones2[2] = {0., 1.}, ones3[3] = {0., 0., 1.};
563:   PetscReal                normal[3] = {0., 0., 0.}, norm;
564:   PetscBool                computeNormal;
565:   PetscInt                 dim, dEx = ex->cdimEx, cStart, cEnd, d;


574:   DMPlexTransformGetDM(tr, &dm);
575:   DMGetDimension(dm, &dim);
576:   DMPlexTransformGetActive(tr, &active);
577:   computeNormal = dim != ex->cdim && !ex->useNormal ? PETSC_TRUE : PETSC_FALSE;
578:   if (computeNormal) {
579:     PetscInt *closure = NULL;
580:     PetscInt  closureSize, cl;

582:     DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
583:     DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);
584:     for (cl = 0; cl < closureSize * 2; cl += 2) {
585:       if ((closure[cl] >= cStart) && (closure[cl] < cEnd)) {
586:         PetscReal cnormal[3] = {0, 0, 0};

588:         DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal);
589:         for (d = 0; d < dEx; ++d) normal[d] += cnormal[d];
590:       }
591:     }
592:     DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);
593:   } else if (ex->useNormal) {
594:     for (d = 0; d < dEx; ++d) normal[d] = ex->normal[d];
595:   } else if (active) { // find an active point in the closure of p and use its coordinate normal as the extrusion direction
596:     PetscInt *closure = NULL;
597:     PetscInt  closureSize, cl, pStart, pEnd;

599:     DMPlexGetDepthStratum(dm, ex->cdimEx - 1, &pStart, &pEnd);
600:     DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);
601:     for (cl = 0; cl < closureSize * 2; cl += 2) {
602:       if ((closure[cl] >= pStart) && (closure[cl] < pEnd)) {
603:         PetscReal       cnormal[3] = {0, 0, 0};
604:         const PetscInt *supp;
605:         PetscInt        suppSize;

607:         DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal);
608:         DMPlexGetSupportSize(dm, closure[cl], &suppSize);
609:         DMPlexGetSupport(dm, closure[cl], &supp);
610:         // Only use external faces, so I can get the orientation from any cell
611:         if (suppSize == 1) {
612:           const PetscInt *cone, *ornt;
613:           PetscInt        coneSize, c;

615:           DMPlexGetConeSize(dm, supp[0], &coneSize);
616:           DMPlexGetCone(dm, supp[0], &cone);
617:           DMPlexGetConeOrientation(dm, supp[0], &ornt);
618:           for (c = 0; c < coneSize; ++c)
619:             if (cone[c] == closure[cl]) break;
621:           if (ornt[c] < 0)
622:             for (d = 0; d < dEx; ++d) cnormal[d] *= -1.;
623:           for (d = 0; d < dEx; ++d) normal[d] += cnormal[d];
624:         }
625:       }
626:     }
627:     DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);
628:   } else if (ex->cdimEx == 2) {
629:     for (d = 0; d < dEx; ++d) normal[d] = ones2[d];
630:   } else if (ex->cdimEx == 3) {
631:     for (d = 0; d < dEx; ++d) normal[d] = ones3[d];
632:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
633:   if (ex->normalFunc) {
634:     PetscScalar n[3];
635:     PetscReal   x[3], dot;

637:     for (d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]);
638:     (*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL);
639:     for (dot = 0, d = 0; d < dEx; d++) dot += PetscRealPart(n[d]) * normal[d];
640:     for (d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]);
641:   }

643:   for (d = 0, norm = 0.0; d < dEx; ++d) norm += PetscSqr(normal[d]);
644:   for (d = 0; d < dEx; ++d) normal[d] *= norm == 0.0 ? 1.0 : 1. / PetscSqrtReal(norm);
645:   for (d = 0; d < dEx; ++d) out[d] = normal[d] * ex->layerPos[r];
646:   for (d = 0; d < ex->cdim; ++d) out[d] += in[d];
647:   return 0;
648: }

650: static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr)
651: {
652:   tr->ops->view                  = DMPlexTransformView_Extrude;
653:   tr->ops->setfromoptions        = DMPlexTransformSetFromOptions_Extrude;
654:   tr->ops->setup                 = DMPlexTransformSetUp_Extrude;
655:   tr->ops->destroy               = DMPlexTransformDestroy_Extrude;
656:   tr->ops->setdimensions         = DMPlexTransformSetDimensions_Extrude;
657:   tr->ops->celltransform         = DMPlexTransformCellTransform_Extrude;
658:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude;
659:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinates_Extrude;
660:   return 0;
661: }

663: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr)
664: {
665:   DMPlexTransform_Extrude *ex;
666:   DM                       dm;
667:   PetscInt                 dim;

670:   PetscNew(&ex);
671:   tr->data      = ex;
672:   ex->thickness = 1.;
673:   ex->useTensor = PETSC_TRUE;
674:   ex->symmetric = PETSC_FALSE;
675:   ex->useNormal = PETSC_FALSE;
676:   ex->layerPos  = NULL;
677:   DMPlexTransformGetDM(tr, &dm);
678:   DMGetDimension(dm, &dim);
679:   DMGetCoordinateDim(dm, &ex->cdim);
680:   DMPlexTransformExtrudeSetLayers(tr, 1);
681:   DMPlexTransformInitialize_Extrude(tr);
682:   return 0;
683: }

685: /*@
686:   DMPlexTransformExtrudeGetLayers - Get the number of extruded layers.

688:   Not collective

690:   Input Parameter:
691: . tr  - The DMPlexTransform

693:   Output Parameter:
694: . layers - The number of layers

696:   Level: intermediate

698: .seealso: `DMPlexTransformExtrudeSetLayers()`
699: @*/
700: PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers)
701: {
702:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

706:   *layers = ex->layers;
707:   return 0;
708: }

710: /*@
711:   DMPlexTransformExtrudeSetLayers - Set the number of extruded layers.

713:   Not collective

715:   Input Parameters:
716: + tr  - The DMPlexTransform
717: - layers - The number of layers

719:   Level: intermediate

721: .seealso: `DMPlexTransformExtrudeGetLayers()`
722: @*/
723: PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers)
724: {
725:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

728:   ex->layers = layers;
729:   PetscFree(ex->layerPos);
730:   PetscCalloc1(ex->layers + 1, &ex->layerPos);
731:   return 0;
732: }

734: /*@
735:   DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers

737:   Not collective

739:   Input Parameter:
740: . tr  - The DMPlexTransform

742:   Output Parameter:
743: . thickness - The total thickness of the layers

745:   Level: intermediate

747: .seealso: `DMPlexTransformExtrudeSetThickness()`
748: @*/
749: PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness)
750: {
751:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

755:   *thickness = ex->thickness;
756:   return 0;
757: }

759: /*@
760:   DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers

762:   Not collective

764:   Input Parameters:
765: + tr  - The DMPlexTransform
766: - thickness - The total thickness of the layers

768:   Level: intermediate

770: .seealso: `DMPlexTransformExtrudeGetThickness()`
771: @*/
772: PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness)
773: {
774:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

778:   ex->thickness = thickness;
779:   return 0;
780: }

782: /*@
783:   DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells

785:   Not collective

787:   Input Parameter:
788: . tr  - The DMPlexTransform

790:   Output Parameter:
791: . useTensor - The flag to use tensor cells

793:   Note: This flag determines the orientation behavior of the created points. For example, if tensor is PETSC_TRUE, then DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT, DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.

795:   Level: intermediate

797: .seealso: `DMPlexTransformExtrudeSetTensor()`
798: @*/
799: PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
800: {
801:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

805:   *useTensor = ex->useTensor;
806:   return 0;
807: }

809: /*@
810:   DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells

812:   Not collective

814:   Input Parameters:
815: + tr  - The DMPlexTransform
816: - useTensor - The flag for tensor cells

818:   Note: This flag determines the orientation behavior of the created points. For example, if tensor is PETSC_TRUE, then DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT, DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.

820:   Level: intermediate

822: .seealso: `DMPlexTransformExtrudeGetTensor()`
823: @*/
824: PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
825: {
826:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

829:   ex->useTensor = useTensor;
830:   return 0;
831: }

833: /*@
834:   DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface

836:   Not collective

838:   Input Parameter:
839: . tr  - The DMPlexTransform

841:   Output Parameter:
842: . symmetric - The flag to extrude symmetrically

844:   Level: intermediate

846: .seealso: `DMPlexTransformExtrudeSetSymmetric()`
847: @*/
848: PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric)
849: {
850:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

854:   *symmetric = ex->symmetric;
855:   return 0;
856: }

858: /*@
859:   DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface

861:   Not collective

863:   Input Parameters:
864: + tr  - The DMPlexTransform
865: - symmetric - The flag to extrude symmetrically

867:   Level: intermediate

869: .seealso: `DMPlexTransformExtrudeGetSymmetric()`
870: @*/
871: PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric)
872: {
873:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

876:   ex->symmetric = symmetric;
877:   return 0;
878: }

880: /*@
881:   DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector

883:   Not collective

885:   Input Parameter:
886: . tr  - The DMPlexTransform

888:   Output Parameter:
889: . normal - The extrusion direction

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

893:   Level: intermediate

895: .seealso: `DMPlexTransformExtrudeSetNormal()`
896: @*/
897: PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[])
898: {
899:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
900:   PetscInt                 d;

903:   if (ex->useNormal) {
904:     for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d];
905:   } else {
906:     for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.;
907:   }
908:   return 0;
909: }

911: /*@
912:   DMPlexTransformExtrudeSetNormal - Set the extrusion normal

914:   Not collective

916:   Input Parameters:
917: + tr     - The DMPlexTransform
918: - normal - The extrusion direction

920:   Level: intermediate

922: .seealso: `DMPlexTransformExtrudeGetNormal()`
923: @*/
924: PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[])
925: {
926:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
927:   PetscInt                 d;

930:   ex->useNormal = PETSC_TRUE;
931:   for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d];
932:   return 0;
933: }

935: /*@C
936:   DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal

938:   Not collective

940:   Input Parameters:
941: + tr     - The DMPlexTransform
942: - normalFunc - A function determining the extrusion direction

944:   Note: The calling sequence for the function is normalFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx)
945: $ dim  - The coordinate dimension of the original mesh (usually a surface)
946: $ time - The current time, or 0.
947: $ x    - The location of the current normal, in the coordinate space of the original mesh
948: $ r    - The extrusion replica number (layer number) of this point
949: $ u    - The user provides the computed normal on output; the sign and magnitude is not significant
950: $ ctx  - An optional user context

952:   Level: intermediate

954: .seealso: `DMPlexTransformExtrudeGetNormal()`
955: @*/
956: PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFunc normalFunc)
957: {
958:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;

961:   ex->normalFunc = normalFunc;
962:   return 0;
963: }

965: /*@
966:   DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer

968:   Not collective

970:   Input Parameters:
971: + tr  - The DMPlexTransform
972: . Nth - The number of thicknesses
973: - thickness - The array of thicknesses

975:   Level: intermediate

977: .seealso: `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeGetThickness()`
978: @*/
979: PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[])
980: {
981:   DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
982:   PetscInt                 t;

986:   ex->Nth = PetscMin(Nth, ex->layers);
987:   PetscFree(ex->thicknesses);
988:   PetscMalloc1(ex->Nth, &ex->thicknesses);
989:   for (t = 0; t < ex->Nth; ++t) {
991:     ex->thicknesses[t] = thicknesses[t];
992:   }
993:   return 0;
994: }