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