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