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