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: %D\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(PetscOptionItems *PetscOptionsObject, DMPlexTransform tr)
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];
33: PetscOptionsHead(PetscOptionsObject, "DMPlexTransform Extrusion Options");
34: PetscOptionsBoundedInt("-dm_plex_transform_extrude_layers", "Number of layers to extrude", "", ex->layers, &nl, &flg, 1);
35: if (flg) DMPlexTransformExtrudeSetLayers(tr, nl);
36: PetscOptionsReal("-dm_plex_transform_extrude_thickness", "Total thickness of extruded layers", "", ex->thickness, &th, &flg);
37: if (flg) DMPlexTransformExtrudeSetThickness(tr, th);
38: PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg);
39: if (flg) DMPlexTransformExtrudeSetTensor(tr, tensor);
40: PetscOptionsBool("-dm_plex_transform_extrude_symmetric", "Extrude layers symmetrically about the surface", "", ex->symmetric, &sym, &flg);
41: if (flg) DMPlexTransformExtrudeSetSymmetric(tr, sym);
42: Nc = 3;
43: PetscOptionsRealArray("-dm_plex_transform_extrude_normal", "Input normal vector for extrusion", "DMPlexTransformExtrudeSetNormal", normal, &Nc, &flg);
44: if (flg) {
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: PetscOptionsTail();
64: return 0;
65: }
67: static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm)
68: {
69: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
70: PetscInt dim;
72: DMGetDimension(dm, &dim);
73: DMSetDimension(tdm, dim+1);
74: DMSetCoordinateDim(tdm, ex->cdimEx);
75: return 0;
76: }
78: /*
79: The refine types for extrusion are:
81: ct: For any normally extruded point
82: ct + 100: For any point which should just return itself
83: */
84: static PetscErrorCode DMPlexTransformSetUp_Extrude(DMPlexTransform tr)
85: {
86: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
87: DM dm;
88: DMLabel active;
89: DMPolytopeType ct;
90: PetscInt Nl = ex->layers, l, i, ict, Nc, No, coff, ooff;
92: DMPlexTransformGetDM(tr, &dm);
93: DMPlexTransformGetActive(tr, &active);
94: if (active) {
95: DMLabel celltype;
96: PetscInt pStart, pEnd, p;
98: DMPlexGetCellTypeLabel(dm, &celltype);
99: DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType);
100: DMPlexGetChart(dm, &pStart, &pEnd);
101: for (p = pStart; p < pEnd; ++p) {
102: PetscInt ct, val;
104: DMLabelGetValue(celltype, p, &ct);
105: DMLabelGetValue(active, p, &val);
106: if (val < 0) {DMLabelSetValue(tr->trType, p, ct + 100);}
107: else {DMLabelSetValue(tr->trType, p, ct);}
108: }
109: }
110: 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);
111: for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
112: ex->Nt[ict] = -1;
113: ex->target[ict] = NULL;
114: ex->size[ict] = NULL;
115: ex->cone[ict] = NULL;
116: ex->ornt[ict] = NULL;
117: }
118: /* DM_POLYTOPE_POINT produces Nl+1 points and Nl segments/tensor segments */
119: ct = DM_POLYTOPE_POINT;
120: ex->Nt[ct] = 2;
121: Nc = 6*Nl;
122: No = 2*Nl;
123: PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);
124: ex->target[ct][0] = DM_POLYTOPE_POINT;
125: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
126: ex->size[ct][0] = Nl+1;
127: ex->size[ct][1] = Nl;
128: /* cones for segments/tensor segments */
129: for (i = 0; i < Nl; ++i) {
130: ex->cone[ct][6*i+0] = DM_POLYTOPE_POINT;
131: ex->cone[ct][6*i+1] = 0;
132: ex->cone[ct][6*i+2] = i;
133: ex->cone[ct][6*i+3] = DM_POLYTOPE_POINT;
134: ex->cone[ct][6*i+4] = 0;
135: ex->cone[ct][6*i+5] = i+1;
136: }
137: for (i = 0; i < No; ++i) ex->ornt[ct][i] = 0;
138: /* DM_POLYTOPE_SEGMENT produces Nl+1 segments and Nl quads/tensor quads */
139: ct = DM_POLYTOPE_SEGMENT;
140: ex->Nt[ct] = 2;
141: Nc = 8*(Nl+1) + 14*Nl;
142: No = 2*(Nl+1) + 4*Nl;
143: PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);
144: ex->target[ct][0] = DM_POLYTOPE_SEGMENT;
145: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
146: ex->size[ct][0] = Nl+1;
147: ex->size[ct][1] = Nl;
148: /* cones for segments */
149: for (i = 0; i < Nl+1; ++i) {
150: ex->cone[ct][8*i+0] = DM_POLYTOPE_POINT;
151: ex->cone[ct][8*i+1] = 1;
152: ex->cone[ct][8*i+2] = 0;
153: ex->cone[ct][8*i+3] = i;
154: ex->cone[ct][8*i+4] = DM_POLYTOPE_POINT;
155: ex->cone[ct][8*i+5] = 1;
156: ex->cone[ct][8*i+6] = 1;
157: ex->cone[ct][8*i+7] = i;
158: }
159: for (i = 0; i < 2*(Nl+1); ++i) ex->ornt[ct][i] = 0;
160: /* cones for quads/tensor quads */
161: coff = 8*(Nl+1);
162: ooff = 2*(Nl+1);
163: for (i = 0; i < Nl; ++i) {
164: if (ex->useTensor) {
165: ex->cone[ct][coff+14*i+0] = DM_POLYTOPE_SEGMENT;
166: ex->cone[ct][coff+14*i+1] = 0;
167: ex->cone[ct][coff+14*i+2] = i;
168: ex->cone[ct][coff+14*i+3] = DM_POLYTOPE_SEGMENT;
169: ex->cone[ct][coff+14*i+4] = 0;
170: ex->cone[ct][coff+14*i+5] = i+1;
171: ex->cone[ct][coff+14*i+6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
172: ex->cone[ct][coff+14*i+7] = 1;
173: ex->cone[ct][coff+14*i+8] = 0;
174: ex->cone[ct][coff+14*i+9] = i;
175: ex->cone[ct][coff+14*i+10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
176: ex->cone[ct][coff+14*i+11] = 1;
177: ex->cone[ct][coff+14*i+12] = 1;
178: ex->cone[ct][coff+14*i+13] = i;
179: ex->ornt[ct][ooff+4*i+0] = 0;
180: ex->ornt[ct][ooff+4*i+1] = 0;
181: ex->ornt[ct][ooff+4*i+2] = 0;
182: ex->ornt[ct][ooff+4*i+3] = 0;
183: } else {
184: ex->cone[ct][coff+14*i+0] = DM_POLYTOPE_SEGMENT;
185: ex->cone[ct][coff+14*i+1] = 0;
186: ex->cone[ct][coff+14*i+2] = i;
187: ex->cone[ct][coff+14*i+3] = DM_POLYTOPE_SEGMENT;
188: ex->cone[ct][coff+14*i+4] = 1;
189: ex->cone[ct][coff+14*i+5] = 1;
190: ex->cone[ct][coff+14*i+6] = i;
191: ex->cone[ct][coff+14*i+7] = DM_POLYTOPE_SEGMENT;
192: ex->cone[ct][coff+14*i+8] = 0;
193: ex->cone[ct][coff+14*i+9] = i+1;
194: ex->cone[ct][coff+14*i+10] = DM_POLYTOPE_SEGMENT;
195: ex->cone[ct][coff+14*i+11] = 1;
196: ex->cone[ct][coff+14*i+12] = 0;
197: ex->cone[ct][coff+14*i+13] = i;
198: ex->ornt[ct][ooff+4*i+0] = 0;
199: ex->ornt[ct][ooff+4*i+1] = 0;
200: ex->ornt[ct][ooff+4*i+2] = -1;
201: ex->ornt[ct][ooff+4*i+3] = -1;
202: }
203: }
204: /* DM_POLYTOPE_TRIANGLE produces Nl+1 triangles and Nl triangular prisms/tensor triangular prisms */
205: ct = DM_POLYTOPE_TRIANGLE;
206: ex->Nt[ct] = 2;
207: Nc = 12*(Nl+1) + 18*Nl;
208: No = 3*(Nl+1) + 5*Nl;
209: PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);
210: ex->target[ct][0] = DM_POLYTOPE_TRIANGLE;
211: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM;
212: ex->size[ct][0] = Nl+1;
213: ex->size[ct][1] = Nl;
214: /* cones for triangles */
215: for (i = 0; i < Nl+1; ++i) {
216: ex->cone[ct][12*i+0] = DM_POLYTOPE_SEGMENT;
217: ex->cone[ct][12*i+1] = 1;
218: ex->cone[ct][12*i+2] = 0;
219: ex->cone[ct][12*i+3] = i;
220: ex->cone[ct][12*i+4] = DM_POLYTOPE_SEGMENT;
221: ex->cone[ct][12*i+5] = 1;
222: ex->cone[ct][12*i+6] = 1;
223: ex->cone[ct][12*i+7] = i;
224: ex->cone[ct][12*i+8] = DM_POLYTOPE_SEGMENT;
225: ex->cone[ct][12*i+9] = 1;
226: ex->cone[ct][12*i+10] = 2;
227: ex->cone[ct][12*i+11] = i;
228: }
229: for (i = 0; i < 3*(Nl+1); ++i) ex->ornt[ct][i] = 0;
230: /* cones for triangular prisms/tensor triangular prisms */
231: coff = 12*(Nl+1);
232: ooff = 3*(Nl+1);
233: for (i = 0; i < Nl; ++i) {
234: if (ex->useTensor) {
235: ex->cone[ct][coff+18*i+0] = DM_POLYTOPE_TRIANGLE;
236: ex->cone[ct][coff+18*i+1] = 0;
237: ex->cone[ct][coff+18*i+2] = i;
238: ex->cone[ct][coff+18*i+3] = DM_POLYTOPE_TRIANGLE;
239: ex->cone[ct][coff+18*i+4] = 0;
240: ex->cone[ct][coff+18*i+5] = i+1;
241: ex->cone[ct][coff+18*i+6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
242: ex->cone[ct][coff+18*i+7] = 1;
243: ex->cone[ct][coff+18*i+8] = 0;
244: ex->cone[ct][coff+18*i+9] = i;
245: ex->cone[ct][coff+18*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
246: ex->cone[ct][coff+18*i+11] = 1;
247: ex->cone[ct][coff+18*i+12] = 1;
248: ex->cone[ct][coff+18*i+13] = i;
249: ex->cone[ct][coff+18*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
250: ex->cone[ct][coff+18*i+15] = 1;
251: ex->cone[ct][coff+18*i+16] = 2;
252: ex->cone[ct][coff+18*i+17] = i;
253: ex->ornt[ct][ooff+5*i+0] = 0;
254: ex->ornt[ct][ooff+5*i+1] = 0;
255: ex->ornt[ct][ooff+5*i+2] = 0;
256: ex->ornt[ct][ooff+5*i+3] = 0;
257: ex->ornt[ct][ooff+5*i+4] = 0;
258: } else {
259: ex->cone[ct][coff+18*i+0] = DM_POLYTOPE_TRIANGLE;
260: ex->cone[ct][coff+18*i+1] = 0;
261: ex->cone[ct][coff+18*i+2] = i;
262: ex->cone[ct][coff+18*i+3] = DM_POLYTOPE_TRIANGLE;
263: ex->cone[ct][coff+18*i+4] = 0;
264: ex->cone[ct][coff+18*i+5] = i+1;
265: ex->cone[ct][coff+18*i+6] = DM_POLYTOPE_QUADRILATERAL;
266: ex->cone[ct][coff+18*i+7] = 1;
267: ex->cone[ct][coff+18*i+8] = 0;
268: ex->cone[ct][coff+18*i+9] = i;
269: ex->cone[ct][coff+18*i+10] = DM_POLYTOPE_QUADRILATERAL;
270: ex->cone[ct][coff+18*i+11] = 1;
271: ex->cone[ct][coff+18*i+12] = 1;
272: ex->cone[ct][coff+18*i+13] = i;
273: ex->cone[ct][coff+18*i+14] = DM_POLYTOPE_QUADRILATERAL;
274: ex->cone[ct][coff+18*i+15] = 1;
275: ex->cone[ct][coff+18*i+16] = 2;
276: ex->cone[ct][coff+18*i+17] = i;
277: ex->ornt[ct][ooff+5*i+0] = -2;
278: ex->ornt[ct][ooff+5*i+1] = 0;
279: ex->ornt[ct][ooff+5*i+2] = 0;
280: ex->ornt[ct][ooff+5*i+3] = 0;
281: ex->ornt[ct][ooff+5*i+4] = 0;
282: }
283: }
284: /* DM_POLYTOPE_QUADRILATERAL produces Nl+1 quads and Nl hexes/tensor hexes */
285: ct = DM_POLYTOPE_QUADRILATERAL;
286: ex->Nt[ct] = 2;
287: Nc = 16*(Nl+1) + 22*Nl;
288: No = 4*(Nl+1) + 6*Nl;
289: PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);
290: ex->target[ct][0] = DM_POLYTOPE_QUADRILATERAL;
291: ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON;
292: ex->size[ct][0] = Nl+1;
293: ex->size[ct][1] = Nl;
294: /* cones for quads */
295: for (i = 0; i < Nl+1; ++i) {
296: ex->cone[ct][16*i+0] = DM_POLYTOPE_SEGMENT;
297: ex->cone[ct][16*i+1] = 1;
298: ex->cone[ct][16*i+2] = 0;
299: ex->cone[ct][16*i+3] = i;
300: ex->cone[ct][16*i+4] = DM_POLYTOPE_SEGMENT;
301: ex->cone[ct][16*i+5] = 1;
302: ex->cone[ct][16*i+6] = 1;
303: ex->cone[ct][16*i+7] = i;
304: ex->cone[ct][16*i+8] = DM_POLYTOPE_SEGMENT;
305: ex->cone[ct][16*i+9] = 1;
306: ex->cone[ct][16*i+10] = 2;
307: ex->cone[ct][16*i+11] = i;
308: ex->cone[ct][16*i+12] = DM_POLYTOPE_SEGMENT;
309: ex->cone[ct][16*i+13] = 1;
310: ex->cone[ct][16*i+14] = 3;
311: ex->cone[ct][16*i+15] = i;
312: }
313: for (i = 0; i < 4*(Nl+1); ++i) ex->ornt[ct][i] = 0;
314: /* cones for hexes/tensor hexes */
315: coff = 16*(Nl+1);
316: ooff = 4*(Nl+1);
317: for (i = 0; i < Nl; ++i) {
318: if (ex->useTensor) {
319: ex->cone[ct][coff+22*i+0] = DM_POLYTOPE_QUADRILATERAL;
320: ex->cone[ct][coff+22*i+1] = 0;
321: ex->cone[ct][coff+22*i+2] = i;
322: ex->cone[ct][coff+22*i+3] = DM_POLYTOPE_QUADRILATERAL;
323: ex->cone[ct][coff+22*i+4] = 0;
324: ex->cone[ct][coff+22*i+5] = i+1;
325: ex->cone[ct][coff+22*i+6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
326: ex->cone[ct][coff+22*i+7] = 1;
327: ex->cone[ct][coff+22*i+8] = 0;
328: ex->cone[ct][coff+22*i+9] = i;
329: ex->cone[ct][coff+22*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
330: ex->cone[ct][coff+22*i+11] = 1;
331: ex->cone[ct][coff+22*i+12] = 1;
332: ex->cone[ct][coff+22*i+13] = i;
333: ex->cone[ct][coff+22*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
334: ex->cone[ct][coff+22*i+15] = 1;
335: ex->cone[ct][coff+22*i+16] = 2;
336: ex->cone[ct][coff+22*i+17] = i;
337: ex->cone[ct][coff+22*i+18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
338: ex->cone[ct][coff+22*i+19] = 1;
339: ex->cone[ct][coff+22*i+20] = 3;
340: ex->cone[ct][coff+22*i+21] = i;
341: ex->ornt[ct][ooff+6*i+0] = 0;
342: ex->ornt[ct][ooff+6*i+1] = 0;
343: ex->ornt[ct][ooff+6*i+2] = 0;
344: ex->ornt[ct][ooff+6*i+3] = 0;
345: ex->ornt[ct][ooff+6*i+4] = 0;
346: ex->ornt[ct][ooff+6*i+5] = 0;
347: } else {
348: ex->cone[ct][coff+22*i+0] = DM_POLYTOPE_QUADRILATERAL;
349: ex->cone[ct][coff+22*i+1] = 0;
350: ex->cone[ct][coff+22*i+2] = i;
351: ex->cone[ct][coff+22*i+3] = DM_POLYTOPE_QUADRILATERAL;
352: ex->cone[ct][coff+22*i+4] = 0;
353: ex->cone[ct][coff+22*i+5] = i+1;
354: ex->cone[ct][coff+22*i+6] = DM_POLYTOPE_QUADRILATERAL;
355: ex->cone[ct][coff+22*i+7] = 1;
356: ex->cone[ct][coff+22*i+8] = 0;
357: ex->cone[ct][coff+22*i+9] = i;
358: ex->cone[ct][coff+22*i+10] = DM_POLYTOPE_QUADRILATERAL;
359: ex->cone[ct][coff+22*i+11] = 1;
360: ex->cone[ct][coff+22*i+12] = 2;
361: ex->cone[ct][coff+22*i+13] = i;
362: ex->cone[ct][coff+22*i+14] = DM_POLYTOPE_QUADRILATERAL;
363: ex->cone[ct][coff+22*i+15] = 1;
364: ex->cone[ct][coff+22*i+16] = 1;
365: ex->cone[ct][coff+22*i+17] = i;
366: ex->cone[ct][coff+22*i+18] = DM_POLYTOPE_QUADRILATERAL;
367: ex->cone[ct][coff+22*i+19] = 1;
368: ex->cone[ct][coff+22*i+20] = 3;
369: ex->cone[ct][coff+22*i+21] = i;
370: ex->ornt[ct][ooff+6*i+0] = -2;
371: ex->ornt[ct][ooff+6*i+1] = 0;
372: ex->ornt[ct][ooff+6*i+2] = 0;
373: ex->ornt[ct][ooff+6*i+3] = 0;
374: ex->ornt[ct][ooff+6*i+4] = 0;
375: ex->ornt[ct][ooff+6*i+5] = 1;
376: }
377: }
378: /* Layers positions */
379: if (!ex->Nth) {
380: if (ex->symmetric) for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness*l)/ex->layers - ex->thickness/2;
381: else for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness*l)/ex->layers;
382: } else {
383: ex->thickness = 0.;
384: ex->layerPos[0] = 0.;
385: for (l = 0; l < ex->layers; ++l) {
386: const PetscReal t = ex->thicknesses[PetscMin(l, ex->Nth-1)];
387: ex->layerPos[l+1] = ex->layerPos[l] + t;
388: ex->thickness += t;
389: }
390: if (ex->symmetric) for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] -= ex->thickness/2.;
391: }
392: return 0;
393: }
395: static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr)
396: {
397: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
398: PetscInt ct;
400: for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {
401: PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct]);
402: }
403: PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt);
404: PetscFree(ex->layerPos);
405: PetscFree(ex);
406: return 0;
407: }
409: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
410: {
411: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
412: DMLabel trType = tr->trType;
413: PetscInt rt;
416: *rnew = r;
417: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
418: if (!so) return 0;
419: if (trType) {
420: DMLabelGetValue(tr->trType, sp, &rt);
421: if (rt >= 100) return 0;
422: }
423: if (ex->useTensor) {
424: switch (sct) {
425: case DM_POLYTOPE_POINT: break;
426: case DM_POLYTOPE_SEGMENT:
427: switch (tct) {
428: case DM_POLYTOPE_SEGMENT: break;
429: case DM_POLYTOPE_SEG_PRISM_TENSOR:
430: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
431: break;
432: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
433: }
434: break;
435: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
436: }
437: } else {
438: switch (sct) {
439: case DM_POLYTOPE_POINT: break;
440: case DM_POLYTOPE_SEGMENT:
441: switch (tct) {
442: case DM_POLYTOPE_SEGMENT: break;
443: case DM_POLYTOPE_QUADRILATERAL:
444: *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0);
445: break;
446: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
447: }
448: break;
449: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
450: }
451: }
452: return 0;
453: }
455: static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
456: {
457: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
458: DMLabel trType = tr->trType;
459: PetscBool ignore = PETSC_FALSE, identity = PETSC_FALSE;
460: PetscInt val = 0;
462: if (trType) {
463: DMLabelGetValue(trType, p, &val);
464: identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE;
465: } else {
466: ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE;
467: }
468: if (rt) *rt = val;
469: if (ignore) {
470: /* Ignore cells that cannot be extruded */
471: *Nt = 0;
472: *target = NULL;
473: *size = NULL;
474: *cone = NULL;
475: *ornt = NULL;
476: } else if (identity) {
477: DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
478: } else {
479: *Nt = ex->Nt[source];
480: *target = ex->target[source];
481: *size = ex->size[source];
482: *cone = ex->cone[source];
483: *ornt = ex->ornt[source];
484: }
485: return 0;
486: }
488: /* Computes new vertex along normal */
489: static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
490: {
491: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
492: DM dm;
493: PetscReal ones2[2] = {0., 1.}, ones3[3] = { 0., 0., 1.};
494: PetscReal normal[3] = {0., 0., 0.}, norm;
495: PetscBool computeNormal;
496: PetscInt dim, dEx = ex->cdimEx, cStart, cEnd, d;
504: DMPlexTransformGetDM(tr, &dm);
505: DMGetDimension(dm, &dim);
506: computeNormal = dim != ex->cdim && !ex->useNormal ? PETSC_TRUE : PETSC_FALSE;
507: if (computeNormal) {
508: PetscInt *closure = NULL;
509: PetscInt closureSize, cl;
511: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
512: DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);
513: for (cl = 0; cl < closureSize*2; cl += 2) {
514: if ((closure[cl] >= cStart) && (closure[cl] < cEnd)) {
515: PetscReal cnormal[3] = {0, 0, 0};
517: DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal);
518: for (d = 0; d < dEx; ++d) normal[d] += cnormal[d];
519: }
520: }
521: DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);
522: } else if (ex->useNormal) {
523: for (d = 0; d < dEx; ++d) normal[d] = ex->normal[d];
524: } else if (ex->cdimEx == 2) {
525: for (d = 0; d < dEx; ++d) normal[d] = ones2[d];
526: } else if (ex->cdimEx == 3) {
527: for (d = 0; d < dEx; ++d) normal[d] = ones3[d];
528: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
529: if (ex->normalFunc) {
530: PetscScalar n[3];
531: PetscReal x[3], dot;
533: for (d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]);
534: (*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL);
535: for (dot=0,d=0; d < dEx; d++) dot += PetscRealPart(n[d]) * normal[d];
536: for (d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]);
537: }
539: for (d = 0, norm = 0.0; d < dEx; ++d) norm += PetscSqr(normal[d]);
540: for (d = 0; d < dEx; ++d) normal[d] *= 1./PetscSqrtReal(norm);
541: for (d = 0; d < dEx; ++d) out[d] = normal[d]*ex->layerPos[r];
542: for (d = 0; d < ex->cdim; ++d) out[d] += in[d];
543: return 0;
544: }
546: static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr)
547: {
548: tr->ops->view = DMPlexTransformView_Extrude;
549: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude;
550: tr->ops->setup = DMPlexTransformSetUp_Extrude;
551: tr->ops->destroy = DMPlexTransformDestroy_Extrude;
552: tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude;
553: tr->ops->celltransform = DMPlexTransformCellTransform_Extrude;
554: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude;
555: tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude;
556: return 0;
557: }
559: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr)
560: {
561: DMPlexTransform_Extrude *ex;
562: DM dm;
563: PetscInt dim;
566: PetscNewLog(tr, &ex);
567: tr->data = ex;
568: ex->thickness = 1.;
569: ex->useTensor = PETSC_TRUE;
570: ex->symmetric = PETSC_FALSE;
571: ex->useNormal = PETSC_FALSE;
572: ex->layerPos = NULL;
573: DMPlexTransformGetDM(tr, &dm);
574: DMGetDimension(dm, &dim);
575: DMGetCoordinateDim(dm, &ex->cdim);
576: ex->cdimEx = ex->cdim == dim ? ex->cdim+1 : ex->cdim;
577: DMPlexTransformExtrudeSetLayers(tr, 1);
578: DMPlexTransformInitialize_Extrude(tr);
579: return 0;
580: }
582: /*@
583: DMPlexTransformExtrudeGetLayers - Get the number of extruded layers.
585: Not collective
587: Input Parameter:
588: . tr - The DMPlexTransform
590: Output Parameter:
591: . layers - The number of layers
593: Level: intermediate
595: .seealso: DMPlexTransformExtrudeSetLayers()
596: @*/
597: PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers)
598: {
599: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
603: *layers = ex->layers;
604: return 0;
605: }
607: /*@
608: DMPlexTransformExtrudeSetLayers - Set the number of extruded layers.
610: Not collective
612: Input Parameters:
613: + tr - The DMPlexTransform
614: - layers - The number of layers
616: Level: intermediate
618: .seealso: DMPlexTransformExtrudeGetLayers()
619: @*/
620: PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers)
621: {
622: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
625: ex->layers = layers;
626: PetscFree(ex->layerPos);
627: PetscCalloc1(ex->layers+1, &ex->layerPos);
628: return 0;
629: }
631: /*@
632: DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers
634: Not collective
636: Input Parameter:
637: . tr - The DMPlexTransform
639: Output Parameter:
640: . thickness - The total thickness of the layers
642: Level: intermediate
644: .seealso: DMPlexTransformExtrudeSetThickness()
645: @*/
646: PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness)
647: {
648: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
652: *thickness = ex->thickness;
653: return 0;
654: }
656: /*@
657: DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers
659: Not collective
661: Input Parameters:
662: + tr - The DMPlexTransform
663: - thickness - The total thickness of the layers
665: Level: intermediate
667: .seealso: DMPlexTransformExtrudeGetThickness()
668: @*/
669: PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness)
670: {
671: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
675: ex->thickness = thickness;
676: return 0;
677: }
679: /*@
680: DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells
682: Not collective
684: Input Parameter:
685: . tr - The DMPlexTransform
687: Output Parameter:
688: . useTensor - The flag to use tensor cells
690: 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.
692: Level: intermediate
694: .seealso: DMPlexTransformExtrudeSetTensor()
695: @*/
696: PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
697: {
698: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
702: *useTensor = ex->useTensor;
703: return 0;
704: }
706: /*@
707: DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells
709: Not collective
711: Input Parameters:
712: + tr - The DMPlexTransform
713: - useTensor - The flag for tensor cells
715: 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.
717: Level: intermediate
719: .seealso: DMPlexTransformExtrudeGetTensor()
720: @*/
721: PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
722: {
723: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
726: ex->useTensor = useTensor;
727: return 0;
728: }
730: /*@
731: DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface
733: Not collective
735: Input Parameter:
736: . tr - The DMPlexTransform
738: Output Parameter:
739: . symmetric - The flag to extrude symmetrically
741: Level: intermediate
743: .seealso: DMPlexTransformExtrudeSetSymmetric()
744: @*/
745: PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric)
746: {
747: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
751: *symmetric = ex->symmetric;
752: return 0;
753: }
755: /*@
756: DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface
758: Not collective
760: Input Parameters:
761: + tr - The DMPlexTransform
762: - symmetric - The flag to extrude symmetrically
764: Level: intermediate
766: .seealso: DMPlexTransformExtrudeGetSymmetric()
767: @*/
768: PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric)
769: {
770: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
773: ex->symmetric = symmetric;
774: return 0;
775: }
777: /*@
778: DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector
780: Not collective
782: Input Parameter:
783: . tr - The DMPlexTransform
785: Output Parameter:
786: . normal - The extrusion direction
788: Note: The user passes in an array, which is filled by the library.
790: Level: intermediate
792: .seealso: DMPlexTransformExtrudeSetNormal()
793: @*/
794: PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[])
795: {
796: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
797: PetscInt d;
800: if (ex->useNormal) {for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d];}
801: else {for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.;}
802: return 0;
803: }
805: /*@
806: DMPlexTransformExtrudeSetNormal - Set the extrusion normal
808: Not collective
810: Input Parameters:
811: + tr - The DMPlexTransform
812: - normal - The extrusion direction
814: Level: intermediate
816: .seealso: DMPlexTransformExtrudeGetNormal()
817: @*/
818: PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[])
819: {
820: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
821: PetscInt d;
824: ex->useNormal = PETSC_TRUE;
825: for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d];
826: return 0;
827: }
829: /*@C
830: DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal
832: Not collective
834: Input Parameters:
835: + tr - The DMPlexTransform
836: - normalFunc - A function determining the extrusion direction
838: Note: The calling sequence for the function is normalFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx)
839: $ dim - The coordinate dimension of the original mesh (usually a surface)
840: $ time - The current time, or 0.
841: $ x - The location of the current normal, in the coordinate space of the original mesh
842: $ r - The extrusion replica number (layer number) of this point
843: $ u - The user provides the computed normal on output; the sign and magnitude is not significant
844: $ ctx - An optional user context
846: Level: intermediate
848: .seealso: DMPlexTransformExtrudeGetNormal()
849: @*/
850: PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFunc normalFunc)
851: {
852: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
855: ex->normalFunc = normalFunc;
856: return 0;
857: }
859: /*@
860: DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer
862: Not collective
864: Input Parameters:
865: + tr - The DMPlexTransform
866: . Nth - The number of thicknesses
867: - thickness - The array of thicknesses
869: Level: intermediate
871: .seealso: DMPlexTransformExtrudeSetThickness(), DMPlexTransformExtrudeGetThickness()
872: @*/
873: PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[])
874: {
875: DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data;
876: PetscInt t;
880: ex->Nth = PetscMin(Nth, ex->layers);
881: PetscFree(ex->thicknesses);
882: PetscMalloc1(ex->Nth, &ex->thicknesses);
883: for (t = 0; t < ex->Nth; ++t) {
885: ex->thicknesses[t] = thicknesses[t];
886: }
887: return 0;
888: }