1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/
2
3 const char *const PlexNormalAlgs[] = {"default", "input", "compute", "compute_bd"};
4
DMPlexTransformView_Extrude(DMPlexTransform tr,PetscViewer viewer)5 static PetscErrorCode DMPlexTransformView_Extrude(DMPlexTransform tr, PetscViewer viewer)
6 {
7 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
8 PetscBool isascii;
9
10 PetscFunctionBegin;
11 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
12 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
13 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
14 if (isascii) {
15 DM dm;
16 DMLabel active;
17 PetscInt dim;
18 const char *name;
19
20 PetscCall(PetscObjectGetName((PetscObject)tr, &name));
21 PetscCall(DMPlexTransformGetDM(tr, &dm));
22 PetscCall(DMGetDimension(dm, &dim));
23 PetscCall(DMPlexTransformGetActive(tr, &active));
24
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 }
36
DMPlexTransformSetFromOptions_Extrude(DMPlexTransform tr,PetscOptionItems PetscOptionsObject)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];
44
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", NULL, funcname, sizeof(funcname), &flg));
65 if (flg) {
66 PetscSimplePointFn *normalFunc;
67
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 }
82
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. */
DMPlexTransformExtrudeComputeExtrusionDim(DMPlexTransform tr)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;
92
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;
102
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 }
128
DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr,DM dm,DM tdm)129 static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm)
130 {
131 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
132 PetscInt dim;
133
134 PetscFunctionBegin;
135 PetscCall(DMGetDimension(dm, &dim));
136 PetscCall(DMSetDimension(tdm, ex->dimEx));
137 PetscCall(DMSetCoordinateDim(tdm, ex->cdimEx));
138 PetscFunctionReturn(PETSC_SUCCESS);
139 }
140
141 /* DM_POLYTOPE_POINT produces
142 Nl+1 points, or Nl points periodically and
143 Nl segments, or tensor segments
144 */
DMPlexTransformExtrudeSetUp_Point(DMPlexTransform_Extrude * ex,PetscInt Nl)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;
150
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 }
172
173 /* DM_POLYTOPE_SEGMENT produces
174 Nl+1 segments, or Nl segments periodically and
175 Nl quads, or tensor quads
176 */
DMPlexTransformExtrudeSetUp_Segment(DMPlexTransform_Extrude * ex,PetscInt Nl)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;
182
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 }
250
251 /* DM_POLYTOPE_TRIANGLE produces
252 Nl+1 triangles, or Nl triangles periodically and
253 Nl triangular prisms/tensor triangular prisms
254 */
DMPlexTransformExtrudeSetUp_Triangle(DMPlexTransform_Extrude * ex,PetscInt Nl)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;
260
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 }
342
343 /* DM_POLYTOPE_QUADRILATERAL produces
344 Nl+1 quads, or Nl quads periodically and
345 Nl hexes/tensor hexes
346 */
DMPlexTransformExtrudeSetUp_Quadrilateral(DMPlexTransform_Extrude * ex,PetscInt Nl)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;
352
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 }
448
449 /*
450 The refine types for extrusion are:
451
452 ct: For any normally extruded point
453 ct + 100: For any point which should just return itself
454 */
DMPlexTransformSetUp_Extrude(DMPlexTransform tr)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;
461
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;
470
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;
476
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;
492
493 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
494 if (size > 1) {
495 PetscSF sf;
496 PetscInt Nr;
497
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;
509
510 PetscMPIInt rank;
511 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
512
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"));
522
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;
528
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;
535
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;
546
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;
556
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));
572
573 Vec g;
574
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 }
612
DMPlexTransformDestroy_Extrude(DMPlexTransform tr)613 static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr)
614 {
615 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
616 PetscInt ct;
617
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 }
629
DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr,DMPolytopeType sct,PetscInt sp,PetscInt so,DMPolytopeType tct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)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, active;
634 PetscBool onBd = PETSC_FALSE;
635 PetscInt rt;
636
637 PetscFunctionBeginHot;
638 *rnew = r;
639 *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
640 PetscCall(DMPlexTransformGetActive(tr, &active));
641 if (!so && !active) PetscFunctionReturn(PETSC_SUCCESS);
642 if (trType) {
643 PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
644 if (rt >= 100) PetscFunctionReturn(PETSC_SUCCESS);
645 }
646 if (active) {
647 // Get orientation of boundary face in cell
648 if (DMPolytopeTypeGetDim(sct) == ex->dimEx - 1) {
649 DM dm;
650 const PetscInt *supp, *cone, *ornt;
651 PetscInt suppSize, coneSize, c;
652
653 PetscCall(DMPlexTransformGetDM(tr, &dm));
654 PetscCall(DMPlexGetSupportSize(dm, sp, &suppSize));
655 PetscCall(DMPlexGetSupport(dm, sp, &supp));
656 PetscCheck(suppSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Source point %" PetscInt_FMT " is not a boundary face", sp);
657 PetscCall(DMPlexGetConeSize(dm, supp[0], &coneSize));
658 PetscCall(DMPlexGetOrientedCone(dm, supp[0], &cone, &ornt));
659 for (c = 0; c < coneSize; ++c)
660 if (cone[c] == sp) break;
661 PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Source point %" PetscInt_FMT " not found in cone of support %" PetscInt_FMT, sp, supp[0]);
662 o = ornt[c];
663 onBd = PETSC_TRUE;
664 PetscCall(DMPlexRestoreOrientedCone(dm, supp[0], &cone, &ornt));
665 }
666 }
667 if (ex->useTensor) {
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_SEG_PRISM_TENSOR:
676 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 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 switch (tct) {
685 case DM_POLYTOPE_TRIANGLE:
686 break;
687 case DM_POLYTOPE_TRI_PRISM_TENSOR:
688 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
689 break;
690 default:
691 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
692 }
693 break;
694 case DM_POLYTOPE_QUADRILATERAL:
695 switch (tct) {
696 case DM_POLYTOPE_QUADRILATERAL:
697 break;
698 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
699 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
700 break;
701 default:
702 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
703 }
704 break;
705 default:
706 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
707 }
708 } else {
709 switch (sct) {
710 case DM_POLYTOPE_POINT:
711 break;
712 case DM_POLYTOPE_SEGMENT:
713 switch (tct) {
714 case DM_POLYTOPE_SEGMENT:
715 break;
716 case DM_POLYTOPE_QUADRILATERAL:
717 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -3) : DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0);
718 break;
719 default:
720 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
721 }
722 break;
723 case DM_POLYTOPE_TRIANGLE:
724 switch (tct) {
725 case DM_POLYTOPE_TRIANGLE:
726 break;
727 case DM_POLYTOPE_TRI_PRISM:
728 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
729 break;
730 default:
731 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
732 }
733 break;
734 case DM_POLYTOPE_QUADRILATERAL:
735 switch (tct) {
736 case DM_POLYTOPE_QUADRILATERAL:
737 break;
738 case DM_POLYTOPE_HEXAHEDRON:
739 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
740 break;
741 default:
742 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
743 }
744 break;
745 default:
746 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
747 }
748 }
749 PetscFunctionReturn(PETSC_SUCCESS);
750 }
751
DMPlexTransformCellTransform_Extrude(DMPlexTransform tr,DMPolytopeType source,PetscInt p,PetscInt * rt,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])752 static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
753 {
754 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
755 DMLabel trType = tr->trType;
756 PetscBool ignore = PETSC_FALSE, identity = PETSC_FALSE;
757 PetscInt val = 0;
758
759 PetscFunctionBegin;
760 if (trType) {
761 PetscCall(DMLabelGetValue(trType, p, &val));
762 identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE;
763 } else {
764 ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE;
765 }
766 if (rt) *rt = val;
767 if (ignore) {
768 /* Ignore cells that cannot be extruded */
769 *Nt = 0;
770 *target = NULL;
771 *size = NULL;
772 *cone = NULL;
773 *ornt = NULL;
774 } else if (identity) {
775 PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
776 } else {
777 *Nt = ex->Nt[source];
778 *target = ex->target[source];
779 *size = ex->size[source];
780 *cone = ex->cone[source];
781 *ornt = ex->ornt[source];
782 }
783 PetscFunctionReturn(PETSC_SUCCESS);
784 }
785
786 /* Computes new vertex along normal */
DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr,DMPolytopeType pct,DMPolytopeType ct,PetscInt p,PetscInt r,PetscInt Nv,PetscInt dE,const PetscScalar in[],PetscScalar out[])787 static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
788 {
789 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
790 DM dm;
791 DMLabel active;
792 PetscReal ones2[2] = {0., 1.};
793 PetscReal ones3[3] = {0., 0., 1.};
794 PetscReal normal[3] = {0., 0., 0.};
795 PetscReal norm = 0.;
796 PetscInt dEx = ex->cdimEx;
797 PetscInt dim, cStart, cEnd;
798
799 PetscFunctionBeginHot;
800 PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]);
801 PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
802 PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv);
803 PetscCheck(dE == ex->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %" PetscInt_FMT " != %" PetscInt_FMT " original dimension", dE, ex->cdim);
804 PetscCheck(dEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", dEx);
805
806 PetscCall(DMPlexTransformGetDM(tr, &dm));
807 PetscCall(DMGetDimension(dm, &dim));
808 PetscCall(DMPlexTransformGetActive(tr, &active));
809 switch (ex->normalAlg) {
810 case NORMAL_DEFAULT:
811 switch (ex->cdimEx) {
812 case 2:
813 for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones2[d];
814 break;
815 case 3:
816 for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones3[d];
817 break;
818 default:
819 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No default normal for dimension %" PetscInt_FMT, ex->cdimEx);
820 }
821 break;
822 case NORMAL_INPUT:
823 for (PetscInt d = 0; d < dEx; ++d) normal[d] = ex->normal[d];
824 break;
825 case NORMAL_COMPUTE: {
826 PetscInt *star = NULL;
827 PetscInt starSize;
828
829 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
830 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star));
831 for (PetscInt st = 0; st < starSize * 2; st += 2) {
832 if ((star[st] >= cStart) && (star[st] < cEnd)) {
833 PetscReal cnormal[3] = {0, 0, 0};
834
835 PetscCall(DMPlexComputeCellGeometryFVM(dm, star[st], NULL, NULL, cnormal));
836 for (PetscInt d = 0; d < dEx; ++d) normal[d] += cnormal[d];
837 }
838 }
839 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star));
840 } break;
841 case NORMAL_COMPUTE_BD: {
842 const PetscScalar *a;
843 PetscScalar *vnormal;
844
845 PetscCall(VecGetArrayRead(ex->vecNormal, &a));
846 PetscCall(DMPlexPointLocalRead(ex->dmNormal, p, a, (void *)&vnormal));
847 for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscRealPart(vnormal[d]);
848 PetscCall(VecRestoreArrayRead(ex->vecNormal, &a));
849 } break;
850 default:
851 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
852 }
853 if (ex->normalFunc) {
854 PetscScalar n[3];
855 PetscReal x[3], dot = 0.;
856
857 for (PetscInt d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]);
858 PetscCall((*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL));
859 for (PetscInt d = 0; d < dEx; ++d) dot += PetscRealPart(n[d]) * normal[d];
860 for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]);
861 }
862 for (PetscInt d = 0; d < dEx; ++d) norm += PetscSqr(normal[d]);
863 for (PetscInt d = 0; d < dEx; ++d) normal[d] *= norm == 0.0 ? 1.0 : 1. / PetscSqrtReal(norm);
864 for (PetscInt d = 0; d < dEx; ++d) out[d] = normal[d] * ex->layerPos[r];
865 for (PetscInt d = 0; d < ex->cdim; ++d) out[d] += in[d];
866 PetscFunctionReturn(PETSC_SUCCESS);
867 }
868
DMPlexTransformInitialize_Extrude(DMPlexTransform tr)869 static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr)
870 {
871 PetscFunctionBegin;
872 tr->ops->view = DMPlexTransformView_Extrude;
873 tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude;
874 tr->ops->setup = DMPlexTransformSetUp_Extrude;
875 tr->ops->destroy = DMPlexTransformDestroy_Extrude;
876 tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude;
877 tr->ops->celltransform = DMPlexTransformCellTransform_Extrude;
878 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude;
879 tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude;
880 PetscFunctionReturn(PETSC_SUCCESS);
881 }
882
DMPlexTransformCreate_Extrude(DMPlexTransform tr)883 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr)
884 {
885 DMPlexTransform_Extrude *ex;
886 DM dm;
887 PetscInt dim;
888
889 PetscFunctionBegin;
890 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
891 PetscCall(PetscNew(&ex));
892 tr->data = ex;
893 ex->thickness = 1.;
894 ex->useTensor = PETSC_TRUE;
895 ex->symmetric = PETSC_FALSE;
896 ex->periodic = PETSC_FALSE;
897 ex->normalAlg = NORMAL_DEFAULT;
898 ex->layerPos = NULL;
899 PetscCall(DMPlexTransformGetDM(tr, &dm));
900 PetscCall(DMGetDimension(dm, &dim));
901 PetscCall(DMGetCoordinateDim(dm, &ex->cdim));
902 PetscCall(DMPlexTransformExtrudeSetLayers(tr, 1));
903 PetscCall(DMPlexTransformInitialize_Extrude(tr));
904 PetscFunctionReturn(PETSC_SUCCESS);
905 }
906
907 /*@
908 DMPlexTransformExtrudeGetLayers - Get the number of extruded layers.
909
910 Not Collective
911
912 Input Parameter:
913 . tr - The `DMPlexTransform`
914
915 Output Parameter:
916 . layers - The number of layers
917
918 Level: intermediate
919
920 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetLayers()`
921 @*/
DMPlexTransformExtrudeGetLayers(DMPlexTransform tr,PetscInt * layers)922 PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers)
923 {
924 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
925
926 PetscFunctionBegin;
927 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
928 PetscAssertPointer(layers, 2);
929 *layers = ex->layers;
930 PetscFunctionReturn(PETSC_SUCCESS);
931 }
932
933 /*@
934 DMPlexTransformExtrudeSetLayers - Set the number of extruded layers.
935
936 Not Collective
937
938 Input Parameters:
939 + tr - The `DMPlexTransform`
940 - layers - The number of layers
941
942 Level: intermediate
943
944 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetLayers()`
945 @*/
DMPlexTransformExtrudeSetLayers(DMPlexTransform tr,PetscInt layers)946 PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers)
947 {
948 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
949
950 PetscFunctionBegin;
951 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
952 ex->layers = layers;
953 PetscCall(PetscFree(ex->layerPos));
954 PetscCall(PetscCalloc1(ex->layers + 1, &ex->layerPos));
955 PetscFunctionReturn(PETSC_SUCCESS);
956 }
957
958 /*@
959 DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers
960
961 Not Collective
962
963 Input Parameter:
964 . tr - The `DMPlexTransform`
965
966 Output Parameter:
967 . thickness - The total thickness of the layers
968
969 Level: intermediate
970
971 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`
972 @*/
DMPlexTransformExtrudeGetThickness(DMPlexTransform tr,PetscReal * thickness)973 PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness)
974 {
975 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
976
977 PetscFunctionBegin;
978 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
979 PetscAssertPointer(thickness, 2);
980 *thickness = ex->thickness;
981 PetscFunctionReturn(PETSC_SUCCESS);
982 }
983
984 /*@
985 DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers
986
987 Not Collective
988
989 Input Parameters:
990 + tr - The `DMPlexTransform`
991 - thickness - The total thickness of the layers
992
993 Level: intermediate
994
995 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetThickness()`
996 @*/
DMPlexTransformExtrudeSetThickness(DMPlexTransform tr,PetscReal thickness)997 PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness)
998 {
999 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1000
1001 PetscFunctionBegin;
1002 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1003 PetscCheck(thickness > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double)thickness);
1004 ex->thickness = thickness;
1005 PetscFunctionReturn(PETSC_SUCCESS);
1006 }
1007
1008 /*@
1009 DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells
1010
1011 Not Collective
1012
1013 Input Parameter:
1014 . tr - The `DMPlexTransform`
1015
1016 Output Parameter:
1017 . useTensor - The flag to use tensor cells
1018
1019 Note:
1020 This flag determines the orientation behavior of the created points.
1021
1022 For example, if tensor is `PETSC_TRUE`, then
1023 .vb
1024 DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1025 DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1026 DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1027 DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1028 .ve
1029
1030 Level: intermediate
1031
1032 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetTensor()`
1033 @*/
DMPlexTransformExtrudeGetTensor(DMPlexTransform tr,PetscBool * useTensor)1034 PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
1035 {
1036 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1037
1038 PetscFunctionBegin;
1039 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1040 PetscAssertPointer(useTensor, 2);
1041 *useTensor = ex->useTensor;
1042 PetscFunctionReturn(PETSC_SUCCESS);
1043 }
1044
1045 /*@
1046 DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells
1047
1048 Not Collective
1049
1050 Input Parameters:
1051 + tr - The `DMPlexTransform`
1052 - useTensor - The flag for tensor cells
1053
1054 Note:
1055 This flag determines the orientation behavior of the created points
1056 For example, if tensor is `PETSC_TRUE`, then
1057 .vb
1058 DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
1059 DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
1060 DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
1061 DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
1062 .ve
1063
1064 Level: intermediate
1065
1066 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetTensor()`
1067 @*/
DMPlexTransformExtrudeSetTensor(DMPlexTransform tr,PetscBool useTensor)1068 PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
1069 {
1070 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1071
1072 PetscFunctionBegin;
1073 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1074 ex->useTensor = useTensor;
1075 PetscFunctionReturn(PETSC_SUCCESS);
1076 }
1077
1078 /*@
1079 DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface
1080
1081 Not Collective
1082
1083 Input Parameter:
1084 . tr - The `DMPlexTransform`
1085
1086 Output Parameter:
1087 . symmetric - The flag to extrude symmetrically
1088
1089 Level: intermediate
1090
1091 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetSymmetric()`, `DMPlexTransformExtrudeGetPeriodic()`
1092 @*/
DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr,PetscBool * symmetric)1093 PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric)
1094 {
1095 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1096
1097 PetscFunctionBegin;
1098 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1099 PetscAssertPointer(symmetric, 2);
1100 *symmetric = ex->symmetric;
1101 PetscFunctionReturn(PETSC_SUCCESS);
1102 }
1103
1104 /*@
1105 DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface
1106
1107 Not Collective
1108
1109 Input Parameters:
1110 + tr - The `DMPlexTransform`
1111 - symmetric - The flag to extrude symmetrically
1112
1113 Level: intermediate
1114
1115 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetSymmetric()`, `DMPlexTransformExtrudeSetPeriodic()`
1116 @*/
DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr,PetscBool symmetric)1117 PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric)
1118 {
1119 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1120
1121 PetscFunctionBegin;
1122 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1123 ex->symmetric = symmetric;
1124 PetscFunctionReturn(PETSC_SUCCESS);
1125 }
1126
1127 /*@
1128 DMPlexTransformExtrudeGetPeriodic - Get the flag to extrude periodically from the initial surface
1129
1130 Not Collective
1131
1132 Input Parameter:
1133 . tr - The `DMPlexTransform`
1134
1135 Output Parameter:
1136 . periodic - The flag to extrude periodically
1137
1138 Level: intermediate
1139
1140 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetPeriodic()`, `DMPlexTransformExtrudeGetSymmetric()`
1141 @*/
DMPlexTransformExtrudeGetPeriodic(DMPlexTransform tr,PetscBool * periodic)1142 PetscErrorCode DMPlexTransformExtrudeGetPeriodic(DMPlexTransform tr, PetscBool *periodic)
1143 {
1144 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1145
1146 PetscFunctionBegin;
1147 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1148 PetscAssertPointer(periodic, 2);
1149 *periodic = ex->periodic;
1150 PetscFunctionReturn(PETSC_SUCCESS);
1151 }
1152
1153 /*@
1154 DMPlexTransformExtrudeSetPeriodic - Set the flag to extrude periodically from the initial surface
1155
1156 Not Collective
1157
1158 Input Parameters:
1159 + tr - The `DMPlexTransform`
1160 - periodic - The flag to extrude periodically
1161
1162 Level: intermediate
1163
1164 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetPeriodic()`, `DMPlexTransformExtrudeSetSymmetric()`
1165 @*/
DMPlexTransformExtrudeSetPeriodic(DMPlexTransform tr,PetscBool periodic)1166 PetscErrorCode DMPlexTransformExtrudeSetPeriodic(DMPlexTransform tr, PetscBool periodic)
1167 {
1168 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1169
1170 PetscFunctionBegin;
1171 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1172 ex->periodic = periodic;
1173 PetscFunctionReturn(PETSC_SUCCESS);
1174 }
1175
1176 /*@
1177 DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector
1178
1179 Not Collective
1180
1181 Input Parameter:
1182 . tr - The `DMPlexTransform`
1183
1184 Output Parameter:
1185 . normal - The extrusion direction
1186
1187 Note:
1188 The user passes in an array, which is filled by the library.
1189
1190 Level: intermediate
1191
1192 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetNormal()`
1193 @*/
DMPlexTransformExtrudeGetNormal(DMPlexTransform tr,PetscReal normal[])1194 PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[])
1195 {
1196 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1197 PetscInt d;
1198
1199 PetscFunctionBegin;
1200 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1201 if (ex->normalAlg == NORMAL_INPUT) {
1202 for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d];
1203 } else {
1204 for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.;
1205 }
1206 PetscFunctionReturn(PETSC_SUCCESS);
1207 }
1208
1209 /*@
1210 DMPlexTransformExtrudeSetNormal - Set the extrusion normal
1211
1212 Not Collective
1213
1214 Input Parameters:
1215 + tr - The `DMPlexTransform`
1216 - normal - The extrusion direction
1217
1218 Level: intermediate
1219
1220 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`
1221 @*/
DMPlexTransformExtrudeSetNormal(DMPlexTransform tr,const PetscReal normal[])1222 PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[])
1223 {
1224 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1225 PetscInt d;
1226
1227 PetscFunctionBegin;
1228 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1229 ex->normalAlg = NORMAL_INPUT;
1230 if (!ex->cdimEx) PetscCall(DMPlexTransformExtrudeComputeExtrusionDim(tr));
1231 for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d];
1232 PetscFunctionReturn(PETSC_SUCCESS);
1233 }
1234
1235 /*@C
1236 DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal
1237
1238 Not Collective
1239
1240 Input Parameters:
1241 + tr - The `DMPlexTransform`
1242 - normalFunc - A function determining the extrusion direction, see `PetscSimplePointFn` for the calling sequence
1243
1244 Level: intermediate
1245
1246 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`, `PetscSimplePointFn`
1247 @*/
DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr,PetscSimplePointFn * normalFunc)1248 PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFn *normalFunc)
1249 {
1250 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1251
1252 PetscFunctionBegin;
1253 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1254 ex->normalFunc = normalFunc;
1255 PetscFunctionReturn(PETSC_SUCCESS);
1256 }
1257
1258 /*@
1259 DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer
1260
1261 Not Collective
1262
1263 Input Parameters:
1264 + tr - The `DMPlexTransform`
1265 . Nth - The number of thicknesses
1266 - thicknesses - The array of thicknesses
1267
1268 Level: intermediate
1269
1270 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeGetThickness()`
1271 @*/
DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr,PetscInt Nth,const PetscReal thicknesses[])1272 PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[])
1273 {
1274 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data;
1275 PetscInt t;
1276
1277 PetscFunctionBegin;
1278 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
1279 PetscCheck(Nth > 0, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %" PetscInt_FMT " must be positive", Nth);
1280 ex->Nth = PetscMin(Nth, ex->layers);
1281 PetscCall(PetscFree(ex->thicknesses));
1282 PetscCall(PetscMalloc1(ex->Nth, &ex->thicknesses));
1283 for (t = 0; t < ex->Nth; ++t) {
1284 PetscCheck(thicknesses[t] > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %" PetscInt_FMT " must be positive", (double)thicknesses[t], t);
1285 ex->thicknesses[t] = thicknesses[t];
1286 }
1287 PetscFunctionReturn(PETSC_SUCCESS);
1288 }
1289