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);break;
 40:           default: break;
 41:         }
 42:       }
 43:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
 44:     }
 45:   }
 46:   /* Cell heights */
 47:   PetscMalloc1(n, &bl->h);
 48:   if (bl->r > 1.) {
 49:     /* initial height h_0 = (r - 1) / (r^{n+1} - 1)
 50:        cell height    h_i = r^i h_0
 51:        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
 52:     */
 53:     PetscReal d = (bl->r - 1.)/(PetscPowRealInt(bl->r, n+1) - 1.);

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

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

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

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

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

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

424:   if (rt) *rt = -1;
425:   if (tr->trType && p >= 0) {
426:     PetscInt val;

428:     DMLabelGetValue(tr->trType, p, &val);
429:     if (val < 0) {
430:       DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
431:       return 0;
432:     }
433:     if (rt) *rt = val;
434:   }
435:   if (bl->Nt[source] < 0) {
436:     DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
437:   } else {
438:     *Nt     = bl->Nt[source];
439:     *target = bl->target[source];
440:     *size   = bl->size[source];
441:     *cone   = bl->cone[source];
442:     *ornt   = bl->ornt[source];
443:   }
444:   return 0;
445: }

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

453:   switch (pct) {
454:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
458:       for (d = 0; d < dE; ++d) out[d] = in[d] + bl->h[r] * (in[d + dE] - in[d]);
459:       break;
460:     default: DMPlexTransformMapCoordinatesBarycenter_Internal(tr, pct, ct, p, r, Nv, dE, in, out);
461:   }
462:   return 0;
463: }

465: static PetscErrorCode DMPlexTransformSetFromOptions_BL(PetscOptionItems *PetscOptionsObject, DMPlexTransform tr)
466: {
467:   DMPlexRefine_BL *bl = (DMPlexRefine_BL *) tr->data;
468:   PetscInt         cells[256], n = 256, i;
469:   PetscBool        flg;

472:   PetscOptionsHead(PetscOptionsObject,"DMPlexTransform Boundary Layer Options");
473:   PetscOptionsBoundedInt("-dm_plex_transform_bl_splits", "Number of divisions of a cell", "", bl->n, &bl->n, NULL, 1);
474:   PetscOptionsReal("-dm_plex_transform_bl_height_factor", "Factor increase for height at each division", "", bl->r, &bl->r, NULL);
475:   PetscOptionsIntArray("-dm_plex_transform_bl_ref_cell", "Mark cells for refinement", "", cells, &n, &flg);
476:   if (flg) {
477:     DMLabel active;

479:     DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active);
480:     for (i = 0; i < n; ++i) DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE);
481:     DMPlexTransformSetActive(tr, active);
482:     DMLabelDestroy(&active);
483:   }
484:   PetscOptionsTail();
485:   return 0;
486: }

488: static PetscErrorCode DMPlexTransformView_BL(DMPlexTransform tr, PetscViewer viewer)
489: {
490:   PetscBool      isascii;

494:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
495:   if (isascii) {
496:     PetscViewerFormat format;
497:     const char       *name;

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

511: static PetscErrorCode DMPlexTransformDestroy_BL(DMPlexTransform tr)
512: {
513:   DMPlexRefine_BL *bl = (DMPlexRefine_BL *) tr->data;
514:   PetscInt         ict;

516:   for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
517:     PetscFree4(bl->target[ict], bl->size[ict], bl->cone[ict], bl->ornt[ict]);
518:   }
519:   PetscFree5(bl->Nt, bl->target, bl->size, bl->cone, bl->ornt);
520:   PetscFree(bl->h);
521:   PetscFree(tr->data);
522:   return 0;
523: }

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

537: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform tr)
538: {
539:   DMPlexRefine_BL *bl;

542:   PetscNewLog(tr, &bl);
543:   tr->data = bl;

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

548:   DMPlexTransformInitialize_BL(tr);
549:   return 0;
550: }