Actual source code: plexrefbl.c

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

  3: static PetscErrorCode DMPlexTransformSetUp_BL(DMPlexTransform tr)
  4: {
  5:   DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
  6:   const PetscInt   n  = bl->n;
  7:   DMPolytopeType   ct;
  8:   DM               dm;
  9:   DMLabel          active;
 10:   PetscInt         Nc, No, coff, i, ict;

 12:   /* If no label is given, split all tensor cells */
 13:   DMPlexTransformGetDM(tr, &dm);
 14:   DMPlexTransformGetActive(tr, &active);
 15:   if (active) {
 16:     IS              refineIS;
 17:     const PetscInt *refineCells;
 18:     PetscInt        c;

 20:     DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType);
 21:     DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS);
 22:     DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc);
 23:     if (refineIS) ISGetIndices(refineIS, &refineCells);
 24:     for (c = 0; c < Nc; ++c) {
 25:       const PetscInt cell    = refineCells[c];
 26:       PetscInt      *closure = NULL;
 27:       PetscInt       Ncl, cl;

 29:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
 30:       for (cl = 0; cl < Ncl; cl += 2) {
 31:         const PetscInt point = closure[cl];

 33:         DMPlexGetCellType(dm, point, &ct);
 34:         switch (ct) {
 35:         case DM_POLYTOPE_POINT_PRISM_TENSOR:
 36:         case DM_POLYTOPE_SEG_PRISM_TENSOR:
 37:         case DM_POLYTOPE_TRI_PRISM_TENSOR:
 38:         case DM_POLYTOPE_QUAD_PRISM_TENSOR:
 39:           DMLabelSetValue(tr->trType, point, 1);
 40:           break;
 41:         default:
 42:           break;
 43:         }
 44:       }
 45:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
 46:     }
 47:   }
 48:   /* Cell heights */
 49:   PetscMalloc1(n, &bl->h);
 50:   if (bl->r > 1.) {
 51:     /* initial height h_0 = (r - 1) / (r^{n+1} - 1)
 52:        cell height    h_i = r^i h_0
 53:        so that \sum^n_{i = 0} h_i = h_0 \sum^n_{i=0} r^i = h_0 (r^{n+1} - 1) / (r - 1) = 1
 54:     */
 55:     PetscReal d = (bl->r - 1.) / (PetscPowRealInt(bl->r, n + 1) - 1.);

 57:     bl->h[0] = d;
 58:     for (i = 1; i < n; ++i) {
 59:       d *= bl->r;
 60:       bl->h[i] = bl->h[i - 1] + d;
 61:     }
 62:   } else {
 63:     /* equal division */
 64:     for (i = 0; i < n; ++i) bl->h[i] = (i + 1.) / (n + 1);
 65:   }

 67:   PetscMalloc5(DM_NUM_POLYTOPES, &bl->Nt, DM_NUM_POLYTOPES, &bl->target, DM_NUM_POLYTOPES, &bl->size, DM_NUM_POLYTOPES, &bl->cone, DM_NUM_POLYTOPES, &bl->ornt);
 68:   for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
 69:     bl->Nt[ict]     = -1;
 70:     bl->target[ict] = NULL;
 71:     bl->size[ict]   = NULL;
 72:     bl->cone[ict]   = NULL;
 73:     bl->ornt[ict]   = NULL;
 74:   }
 75:   /* DM_POLYTOPE_POINT_PRISM_TENSOR produces n points and n+1 tensor segments */
 76:   ct         = DM_POLYTOPE_POINT_PRISM_TENSOR;
 77:   bl->Nt[ct] = 2;
 78:   Nc         = 7 * 2 + 6 * (n - 1);
 79:   No         = 2 * (n + 1);
 80:   PetscMalloc4(bl->Nt[ct], &bl->target[ct], bl->Nt[ct], &bl->size[ct], Nc, &bl->cone[ct], No, &bl->ornt[ct]);
 81:   bl->target[ct][0] = DM_POLYTOPE_POINT;
 82:   bl->target[ct][1] = DM_POLYTOPE_POINT_PRISM_TENSOR;
 83:   bl->size[ct][0]   = n;
 84:   bl->size[ct][1]   = n + 1;
 85:   /*   cones for tensor segments */
 86:   bl->cone[ct][0] = DM_POLYTOPE_POINT;
 87:   bl->cone[ct][1] = 1;
 88:   bl->cone[ct][2] = 0;
 89:   bl->cone[ct][3] = 0;
 90:   bl->cone[ct][4] = DM_POLYTOPE_POINT;
 91:   bl->cone[ct][5] = 0;
 92:   bl->cone[ct][6] = 0;
 93:   for (i = 0; i < n - 1; ++i) {
 94:     bl->cone[ct][7 + 6 * i + 0] = DM_POLYTOPE_POINT;
 95:     bl->cone[ct][7 + 6 * i + 1] = 0;
 96:     bl->cone[ct][7 + 6 * i + 2] = i;
 97:     bl->cone[ct][7 + 6 * i + 3] = DM_POLYTOPE_POINT;
 98:     bl->cone[ct][7 + 6 * i + 4] = 0;
 99:     bl->cone[ct][7 + 6 * i + 5] = i + 1;
100:   }
101:   bl->cone[ct][7 + 6 * (n - 1) + 0] = DM_POLYTOPE_POINT;
102:   bl->cone[ct][7 + 6 * (n - 1) + 1] = 0;
103:   bl->cone[ct][7 + 6 * (n - 1) + 2] = n - 1;
104:   bl->cone[ct][7 + 6 * (n - 1) + 3] = DM_POLYTOPE_POINT;
105:   bl->cone[ct][7 + 6 * (n - 1) + 4] = 1;
106:   bl->cone[ct][7 + 6 * (n - 1) + 5] = 1;
107:   bl->cone[ct][7 + 6 * (n - 1) + 6] = 0;
108:   for (i = 0; i < No; i++) bl->ornt[ct][i] = 0;
109:   /* DM_POLYTOPE_SEG_PRISM_TENSOR produces n segments and n+1 tensor quads */
110:   ct         = DM_POLYTOPE_SEG_PRISM_TENSOR;
111:   bl->Nt[ct] = 2;
112:   Nc         = 8 * n + 15 * 2 + 14 * (n - 1);
113:   No         = 2 * n + 4 * (n + 1);
114:   PetscMalloc4(bl->Nt[ct], &bl->target[ct], bl->Nt[ct], &bl->size[ct], Nc, &bl->cone[ct], No, &bl->ornt[ct]);
115:   bl->target[ct][0] = DM_POLYTOPE_SEGMENT;
116:   bl->target[ct][1] = DM_POLYTOPE_SEG_PRISM_TENSOR;
117:   bl->size[ct][0]   = n;
118:   bl->size[ct][1]   = n + 1;
119:   /*   cones for segments */
120:   for (i = 0; i < n; ++i) {
121:     bl->cone[ct][8 * i + 0] = DM_POLYTOPE_POINT;
122:     bl->cone[ct][8 * i + 1] = 1;
123:     bl->cone[ct][8 * i + 2] = 2;
124:     bl->cone[ct][8 * i + 3] = i;
125:     bl->cone[ct][8 * i + 4] = DM_POLYTOPE_POINT;
126:     bl->cone[ct][8 * i + 5] = 1;
127:     bl->cone[ct][8 * i + 6] = 3;
128:     bl->cone[ct][8 * i + 7] = i;
129:   }
130:   /*   cones for tensor quads */
131:   coff                    = 8 * n;
132:   bl->cone[ct][coff + 0]  = DM_POLYTOPE_SEGMENT;
133:   bl->cone[ct][coff + 1]  = 1;
134:   bl->cone[ct][coff + 2]  = 0;
135:   bl->cone[ct][coff + 3]  = 0;
136:   bl->cone[ct][coff + 4]  = DM_POLYTOPE_SEGMENT;
137:   bl->cone[ct][coff + 5]  = 0;
138:   bl->cone[ct][coff + 6]  = 0;
139:   bl->cone[ct][coff + 7]  = DM_POLYTOPE_POINT_PRISM_TENSOR;
140:   bl->cone[ct][coff + 8]  = 1;
141:   bl->cone[ct][coff + 9]  = 2;
142:   bl->cone[ct][coff + 10] = 0;
143:   bl->cone[ct][coff + 11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
144:   bl->cone[ct][coff + 12] = 1;
145:   bl->cone[ct][coff + 13] = 3;
146:   bl->cone[ct][coff + 14] = 0;
147:   for (i = 0; i < n - 1; ++i) {
148:     bl->cone[ct][coff + 15 + 14 * i + 0]  = DM_POLYTOPE_SEGMENT;
149:     bl->cone[ct][coff + 15 + 14 * i + 1]  = 0;
150:     bl->cone[ct][coff + 15 + 14 * i + 2]  = i;
151:     bl->cone[ct][coff + 15 + 14 * i + 3]  = DM_POLYTOPE_SEGMENT;
152:     bl->cone[ct][coff + 15 + 14 * i + 4]  = 0;
153:     bl->cone[ct][coff + 15 + 14 * i + 5]  = i + 1;
154:     bl->cone[ct][coff + 15 + 14 * i + 6]  = DM_POLYTOPE_POINT_PRISM_TENSOR;
155:     bl->cone[ct][coff + 15 + 14 * i + 7]  = 1;
156:     bl->cone[ct][coff + 15 + 14 * i + 8]  = 2;
157:     bl->cone[ct][coff + 15 + 14 * i + 9]  = i + 1;
158:     bl->cone[ct][coff + 15 + 14 * i + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
159:     bl->cone[ct][coff + 15 + 14 * i + 11] = 1;
160:     bl->cone[ct][coff + 15 + 14 * i + 12] = 3;
161:     bl->cone[ct][coff + 15 + 14 * i + 13] = i + 1;
162:   }
163:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 0]  = DM_POLYTOPE_SEGMENT;
164:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 1]  = 0;
165:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 2]  = n - 1;
166:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 3]  = DM_POLYTOPE_SEGMENT;
167:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 4]  = 1;
168:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 5]  = 1;
169:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 6]  = 0;
170:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 7]  = DM_POLYTOPE_POINT_PRISM_TENSOR;
171:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 8]  = 1;
172:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 9]  = 2;
173:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 10] = n;
174:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
175:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 12] = 1;
176:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 13] = 3;
177:   bl->cone[ct][coff + 15 + 14 * (n - 1) + 14] = n;
178:   for (i = 0; i < No; i++) bl->ornt[ct][i] = 0;
179:   /* DM_POLYTOPE_TRI_PRISM_TENSOR produces n triangles and n+1 tensor triangular prisms */
180:   ct         = DM_POLYTOPE_TRI_PRISM_TENSOR;
181:   bl->Nt[ct] = 2;
182:   Nc         = 12 * n + 19 * 2 + 18 * (n - 1);
183:   No         = 3 * n + 5 * (n + 1);
184:   PetscMalloc4(bl->Nt[ct], &bl->target[ct], bl->Nt[ct], &bl->size[ct], Nc, &bl->cone[ct], No, &bl->ornt[ct]);
185:   bl->target[ct][0] = DM_POLYTOPE_TRIANGLE;
186:   bl->target[ct][1] = DM_POLYTOPE_TRI_PRISM_TENSOR;
187:   bl->size[ct][0]   = n;
188:   bl->size[ct][1]   = n + 1;
189:   /*   cones for triangles */
190:   for (i = 0; i < n; ++i) {
191:     bl->cone[ct][12 * i + 0]  = DM_POLYTOPE_SEGMENT;
192:     bl->cone[ct][12 * i + 1]  = 1;
193:     bl->cone[ct][12 * i + 2]  = 2;
194:     bl->cone[ct][12 * i + 3]  = i;
195:     bl->cone[ct][12 * i + 4]  = DM_POLYTOPE_SEGMENT;
196:     bl->cone[ct][12 * i + 5]  = 1;
197:     bl->cone[ct][12 * i + 6]  = 3;
198:     bl->cone[ct][12 * i + 7]  = i;
199:     bl->cone[ct][12 * i + 8]  = DM_POLYTOPE_SEGMENT;
200:     bl->cone[ct][12 * i + 9]  = 1;
201:     bl->cone[ct][12 * i + 10] = 4;
202:     bl->cone[ct][12 * i + 11] = i;
203:   }
204:   /*   cones for triangular prisms */
205:   coff                    = 12 * n;
206:   bl->cone[ct][coff + 0]  = DM_POLYTOPE_TRIANGLE;
207:   bl->cone[ct][coff + 1]  = 1;
208:   bl->cone[ct][coff + 2]  = 0;
209:   bl->cone[ct][coff + 3]  = 0;
210:   bl->cone[ct][coff + 4]  = DM_POLYTOPE_TRIANGLE;
211:   bl->cone[ct][coff + 5]  = 0;
212:   bl->cone[ct][coff + 6]  = 0;
213:   bl->cone[ct][coff + 7]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
214:   bl->cone[ct][coff + 8]  = 1;
215:   bl->cone[ct][coff + 9]  = 2;
216:   bl->cone[ct][coff + 10] = 0;
217:   bl->cone[ct][coff + 11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
218:   bl->cone[ct][coff + 12] = 1;
219:   bl->cone[ct][coff + 13] = 3;
220:   bl->cone[ct][coff + 14] = 0;
221:   bl->cone[ct][coff + 15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
222:   bl->cone[ct][coff + 16] = 1;
223:   bl->cone[ct][coff + 17] = 4;
224:   bl->cone[ct][coff + 18] = 0;
225:   for (i = 0; i < n - 1; ++i) {
226:     bl->cone[ct][coff + 19 + 18 * i + 0]  = DM_POLYTOPE_TRIANGLE;
227:     bl->cone[ct][coff + 19 + 18 * i + 1]  = 0;
228:     bl->cone[ct][coff + 19 + 18 * i + 2]  = i;
229:     bl->cone[ct][coff + 19 + 18 * i + 3]  = DM_POLYTOPE_TRIANGLE;
230:     bl->cone[ct][coff + 19 + 18 * i + 4]  = 0;
231:     bl->cone[ct][coff + 19 + 18 * i + 5]  = i + 1;
232:     bl->cone[ct][coff + 19 + 18 * i + 6]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
233:     bl->cone[ct][coff + 19 + 18 * i + 7]  = 1;
234:     bl->cone[ct][coff + 19 + 18 * i + 8]  = 2;
235:     bl->cone[ct][coff + 19 + 18 * i + 9]  = i + 1;
236:     bl->cone[ct][coff + 19 + 18 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
237:     bl->cone[ct][coff + 19 + 18 * i + 11] = 1;
238:     bl->cone[ct][coff + 19 + 18 * i + 12] = 3;
239:     bl->cone[ct][coff + 19 + 18 * i + 13] = i + 1;
240:     bl->cone[ct][coff + 19 + 18 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
241:     bl->cone[ct][coff + 19 + 18 * i + 15] = 1;
242:     bl->cone[ct][coff + 19 + 18 * i + 16] = 4;
243:     bl->cone[ct][coff + 19 + 18 * i + 17] = i + 1;
244:   }
245:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 0]  = DM_POLYTOPE_TRIANGLE;
246:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 1]  = 0;
247:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 2]  = n - 1;
248:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 3]  = DM_POLYTOPE_TRIANGLE;
249:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 4]  = 1;
250:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 5]  = 1;
251:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 6]  = 0;
252:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 7]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
253:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 8]  = 1;
254:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 9]  = 2;
255:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 10] = n;
256:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
257:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 12] = 1;
258:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 13] = 3;
259:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 14] = n;
260:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
261:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 16] = 1;
262:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 17] = 4;
263:   bl->cone[ct][coff + 19 + 18 * (n - 1) + 18] = n;
264:   for (i = 0; i < No; ++i) bl->ornt[ct][i] = 0;
265:   /* DM_POLYTOPE_QUAD_PRISM_TENSOR produces n quads and n+1 tensor quad prisms */
266:   ct         = DM_POLYTOPE_QUAD_PRISM_TENSOR;
267:   bl->Nt[ct] = 2;
268:   Nc         = 16 * n + 23 * 2 + 22 * (n - 1);
269:   No         = 4 * n + 6 * (n + 1);
270:   PetscMalloc4(bl->Nt[ct], &bl->target[ct], bl->Nt[ct], &bl->size[ct], Nc, &bl->cone[ct], No, &bl->ornt[ct]);
271:   bl->target[ct][0] = DM_POLYTOPE_QUADRILATERAL;
272:   bl->target[ct][1] = DM_POLYTOPE_QUAD_PRISM_TENSOR;
273:   bl->size[ct][0]   = n;
274:   bl->size[ct][1]   = n + 1;
275:   /*  cones for quads */
276:   for (i = 0; i < n; ++i) {
277:     bl->cone[ct][16 * i + 0]  = DM_POLYTOPE_SEGMENT;
278:     bl->cone[ct][16 * i + 1]  = 1;
279:     bl->cone[ct][16 * i + 2]  = 2;
280:     bl->cone[ct][16 * i + 3]  = i;
281:     bl->cone[ct][16 * i + 4]  = DM_POLYTOPE_SEGMENT;
282:     bl->cone[ct][16 * i + 5]  = 1;
283:     bl->cone[ct][16 * i + 6]  = 3;
284:     bl->cone[ct][16 * i + 7]  = i;
285:     bl->cone[ct][16 * i + 8]  = DM_POLYTOPE_SEGMENT;
286:     bl->cone[ct][16 * i + 9]  = 1;
287:     bl->cone[ct][16 * i + 10] = 4;
288:     bl->cone[ct][16 * i + 11] = i;
289:     bl->cone[ct][16 * i + 12] = DM_POLYTOPE_SEGMENT;
290:     bl->cone[ct][16 * i + 13] = 1;
291:     bl->cone[ct][16 * i + 14] = 5;
292:     bl->cone[ct][16 * i + 15] = i;
293:   }
294:   /*   cones for quad prisms */
295:   coff                    = 16 * n;
296:   bl->cone[ct][coff + 0]  = DM_POLYTOPE_QUADRILATERAL;
297:   bl->cone[ct][coff + 1]  = 1;
298:   bl->cone[ct][coff + 2]  = 0;
299:   bl->cone[ct][coff + 3]  = 0;
300:   bl->cone[ct][coff + 4]  = DM_POLYTOPE_QUADRILATERAL;
301:   bl->cone[ct][coff + 5]  = 0;
302:   bl->cone[ct][coff + 6]  = 0;
303:   bl->cone[ct][coff + 7]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
304:   bl->cone[ct][coff + 8]  = 1;
305:   bl->cone[ct][coff + 9]  = 2;
306:   bl->cone[ct][coff + 10] = 0;
307:   bl->cone[ct][coff + 11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
308:   bl->cone[ct][coff + 12] = 1;
309:   bl->cone[ct][coff + 13] = 3;
310:   bl->cone[ct][coff + 14] = 0;
311:   bl->cone[ct][coff + 15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
312:   bl->cone[ct][coff + 16] = 1;
313:   bl->cone[ct][coff + 17] = 4;
314:   bl->cone[ct][coff + 18] = 0;
315:   bl->cone[ct][coff + 19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
316:   bl->cone[ct][coff + 20] = 1;
317:   bl->cone[ct][coff + 21] = 5;
318:   bl->cone[ct][coff + 22] = 0;
319:   for (i = 0; i < n - 1; ++i) {
320:     bl->cone[ct][coff + 23 + 22 * i + 0]  = DM_POLYTOPE_QUADRILATERAL;
321:     bl->cone[ct][coff + 23 + 22 * i + 1]  = 0;
322:     bl->cone[ct][coff + 23 + 22 * i + 2]  = i;
323:     bl->cone[ct][coff + 23 + 22 * i + 3]  = DM_POLYTOPE_QUADRILATERAL;
324:     bl->cone[ct][coff + 23 + 22 * i + 4]  = 0;
325:     bl->cone[ct][coff + 23 + 22 * i + 5]  = i + 1;
326:     bl->cone[ct][coff + 23 + 22 * i + 6]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
327:     bl->cone[ct][coff + 23 + 22 * i + 7]  = 1;
328:     bl->cone[ct][coff + 23 + 22 * i + 8]  = 2;
329:     bl->cone[ct][coff + 23 + 22 * i + 9]  = i + 1;
330:     bl->cone[ct][coff + 23 + 22 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
331:     bl->cone[ct][coff + 23 + 22 * i + 11] = 1;
332:     bl->cone[ct][coff + 23 + 22 * i + 12] = 3;
333:     bl->cone[ct][coff + 23 + 22 * i + 13] = i + 1;
334:     bl->cone[ct][coff + 23 + 22 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
335:     bl->cone[ct][coff + 23 + 22 * i + 15] = 1;
336:     bl->cone[ct][coff + 23 + 22 * i + 16] = 4;
337:     bl->cone[ct][coff + 23 + 22 * i + 17] = i + 1;
338:     bl->cone[ct][coff + 23 + 22 * i + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
339:     bl->cone[ct][coff + 23 + 22 * i + 19] = 1;
340:     bl->cone[ct][coff + 23 + 22 * i + 20] = 5;
341:     bl->cone[ct][coff + 23 + 22 * i + 21] = i + 1;
342:   }
343:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 0]  = DM_POLYTOPE_QUADRILATERAL;
344:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 1]  = 0;
345:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 2]  = n - 1;
346:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 3]  = DM_POLYTOPE_QUADRILATERAL;
347:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 4]  = 1;
348:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 5]  = 1;
349:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 6]  = 0;
350:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 7]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
351:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 8]  = 1;
352:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 9]  = 2;
353:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 10] = n;
354:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
355:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 12] = 1;
356:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 13] = 3;
357:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 14] = n;
358:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
359:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 16] = 1;
360:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 17] = 4;
361:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 18] = n;
362:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
363:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 20] = 1;
364:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 21] = 5;
365:   bl->cone[ct][coff + 23 + 22 * (n - 1) + 22] = n;
366:   for (i = 0; i < No; ++i) bl->ornt[ct][i] = 0;
367:   return 0;
368: }

370: static PetscErrorCode DMPlexTransformGetSubcellOrientation_BL(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
371: {
372:   DMPlexRefine_BL *bl              = (DMPlexRefine_BL *)tr->data;
373:   const PetscInt   n               = bl->n;
374:   PetscInt         tquad_tquad_o[] = {0, 1, -2, -1, 1, 0, -1, -2, -2, -1, 0, 1, -1, -2, 1, 0};

377:   *rnew = r;
378:   *onew = o;
379:   if (tr->trType) {
380:     PetscInt rt;

382:     DMLabelGetValue(tr->trType, sp, &rt);
383:     if (rt < 0) {
384:       DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
385:       return 0;
386:     }
387:   }
388:   switch (sct) {
389:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
390:     switch (tct) {
391:     case DM_POLYTOPE_POINT:
392:       *rnew = !so ? r : n - 1 - r;
393:       break;
394:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
395:       if (!so) {
396:         *rnew = r;
397:         *onew = o;
398:       } else {
399:         *rnew = n - r;
400:         *onew = -(o + 1);
401:       }
402:     default:
403:       break;
404:     }
405:     break;
406:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
407:     switch (tct) {
408:     case DM_POLYTOPE_SEGMENT:
409:       *onew = (so == 0) || (so == 1) ? o : -(o + 1);
410:       *rnew = (so == 0) || (so == -1) ? r : n - 1 - r;
411:       break;
412:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
413:       *onew = tquad_tquad_o[(so + 2) * 4 + o + 2];
414:       *rnew = (so == 0) || (so == -1) ? r : n - r;
415:       break;
416:     default:
417:       break;
418:     }
419:     break;
420:   default:
421:     DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
422:   }
423:   return 0;
424: }

426: static PetscErrorCode DMPlexTransformCellTransform_BL(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
427: {
428:   DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;

431:   if (rt) *rt = -1;
432:   if (tr->trType && p >= 0) {
433:     PetscInt val;

435:     DMLabelGetValue(tr->trType, p, &val);
436:     if (val < 0) {
437:       DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
438:       return 0;
439:     }
440:     if (rt) *rt = val;
441:   }
442:   if (bl->Nt[source] < 0) {
443:     DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
444:   } else {
445:     *Nt     = bl->Nt[source];
446:     *target = bl->target[source];
447:     *size   = bl->size[source];
448:     *cone   = bl->cone[source];
449:     *ornt   = bl->ornt[source];
450:   }
451:   return 0;
452: }

454: static PetscErrorCode DMPlexTransformMapCoordinates_BL(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
455: {
456:   DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
457:   PetscInt         d;

460:   switch (pct) {
461:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
465:     for (d = 0; d < dE; ++d) out[d] = in[d] + bl->h[r] * (in[d + dE] - in[d]);
466:     break;
467:   default:
468:     DMPlexTransformMapCoordinatesBarycenter_Internal(tr, pct, ct, p, r, Nv, dE, in, out);
469:   }
470:   return 0;
471: }

473: static PetscErrorCode DMPlexTransformSetFromOptions_BL(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
474: {
475:   DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
476:   PetscInt         cells[256], n = 256, i;
477:   PetscBool        flg;

479:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Boundary Layer Options");
480:   PetscOptionsBoundedInt("-dm_plex_transform_bl_splits", "Number of divisions of a cell", "", bl->n, &bl->n, NULL, 1);
481:   PetscOptionsReal("-dm_plex_transform_bl_height_factor", "Factor increase for height at each division", "", bl->r, &bl->r, NULL);
482:   PetscOptionsIntArray("-dm_plex_transform_bl_ref_cell", "Mark cells for refinement", "", cells, &n, &flg);
483:   if (flg) {
484:     DMLabel active;

486:     DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active);
487:     for (i = 0; i < n; ++i) DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE);
488:     DMPlexTransformSetActive(tr, active);
489:     DMLabelDestroy(&active);
490:   }
491:   PetscOptionsHeadEnd();
492:   return 0;
493: }

495: static PetscErrorCode DMPlexTransformView_BL(DMPlexTransform tr, PetscViewer viewer)
496: {
497:   PetscBool isascii;

501:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
502:   if (isascii) {
503:     PetscViewerFormat format;
504:     const char       *name;

506:     PetscObjectGetName((PetscObject)tr, &name);
507:     PetscViewerASCIIPrintf(viewer, "Boundary Layer refinement %s\n", name ? name : "");
508:     PetscViewerGetFormat(viewer, &format);
509:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) DMLabelView(tr->trType, viewer);
510:   } else {
511:     SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
512:   }
513:   return 0;
514: }

516: static PetscErrorCode DMPlexTransformDestroy_BL(DMPlexTransform tr)
517: {
518:   DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
519:   PetscInt         ict;

521:   for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) PetscFree4(bl->target[ict], bl->size[ict], bl->cone[ict], bl->ornt[ict]);
522:   PetscFree5(bl->Nt, bl->target, bl->size, bl->cone, bl->ornt);
523:   PetscFree(bl->h);
524:   PetscFree(tr->data);
525:   return 0;
526: }

528: static PetscErrorCode DMPlexTransformInitialize_BL(DMPlexTransform tr)
529: {
530:   tr->ops->view                  = DMPlexTransformView_BL;
531:   tr->ops->setfromoptions        = DMPlexTransformSetFromOptions_BL;
532:   tr->ops->setup                 = DMPlexTransformSetUp_BL;
533:   tr->ops->destroy               = DMPlexTransformDestroy_BL;
534:   tr->ops->celltransform         = DMPlexTransformCellTransform_BL;
535:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_BL;
536:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinates_BL;
537:   return 0;
538: }

540: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform tr)
541: {
542:   DMPlexRefine_BL *bl;

545:   PetscNew(&bl);
546:   tr->data = bl;

548:   bl->n = 1; /* 1 split -> 2 new cells */
549:   bl->r = 1; /* linear coordinate progression */

551:   DMPlexTransformInitialize_BL(tr);
552:   return 0;
553: }