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;
11: PetscErrorCode ierr;
14: /* If no label is given, split all tensor cells */
15: DMPlexTransformGetDM(tr, &dm);
16: DMPlexTransformGetActive(tr, &active);
17: if (active) {
18: IS refineIS;
19: const PetscInt *refineCells;
20: PetscInt c;
22: DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType);
23: DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS);
24: DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc);
25: if (refineIS) {ISGetIndices(refineIS, &refineCells);}
26: for (c = 0; c < Nc; ++c) {
27: const PetscInt cell = refineCells[c];
28: PetscInt *closure = NULL;
29: PetscInt Ncl, cl;
31: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
32: for (cl = 0; cl < Ncl; cl += 2) {
33: const PetscInt point = closure[cl];
35: DMPlexGetCellType(dm, point, &ct);
36: switch (ct) {
37: case DM_POLYTOPE_POINT_PRISM_TENSOR:
38: case DM_POLYTOPE_SEG_PRISM_TENSOR:
39: case DM_POLYTOPE_TRI_PRISM_TENSOR:
40: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
41: DMLabelSetValue(tr->trType, point, 1);break;
42: default: 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: PetscErrorCode ierr;
375: PetscInt tquad_tquad_o[] = { 0, 1, -2, -1,
376: 1, 0, -1, -2,
377: -2, -1, 0, 1,
378: -1, -2, 1, 0};
381: *rnew = r;
382: *onew = o;
383: if (tr->trType) {
384: PetscInt rt;
386: DMLabelGetValue(tr->trType, sp, &rt);
387: if (rt < 0) {
388: DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
389: return(0);
390: }
391: }
392: switch (sct) {
393: case DM_POLYTOPE_POINT_PRISM_TENSOR:
394: switch (tct) {
395: case DM_POLYTOPE_POINT:
396: *rnew = !so ? r : n-1 - r;
397: break;
398: case DM_POLYTOPE_POINT_PRISM_TENSOR:
399: if (!so) {*rnew = r; *onew = o;}
400: else {*rnew = n - r; *onew = -(o+1);}
401: default: break;
402: }
403: break;
404: case DM_POLYTOPE_SEG_PRISM_TENSOR:
405: switch (tct) {
406: case DM_POLYTOPE_SEGMENT:
407: *onew = (so == 0) || (so == 1) ? o : -(o+1);
408: *rnew = (so == 0) || (so == -1) ? r : n-1 - r;
409: break;
410: case DM_POLYTOPE_SEG_PRISM_TENSOR:
411: *onew = tquad_tquad_o[(so+2)*4+o+2];
412: *rnew = (so == 0) || (so == -1) ? r : n - r;
413: break;
414: default: break;
415: }
416: break;
417: default: DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
418: }
419: return(0);
420: }
422: static PetscErrorCode DMPlexTransformCellTransform_BL(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
423: {
424: DMPlexRefine_BL *bl = (DMPlexRefine_BL *) tr->data;
425: PetscErrorCode ierr;
428: if (rt) *rt = -1;
429: if (tr->trType && p >= 0) {
430: PetscInt val;
432: DMLabelGetValue(tr->trType, p, &val);
433: if (val < 0) {
434: DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
435: return(0);
436: }
437: if (rt) *rt = val;
438: }
439: if (bl->Nt[source] < 0) {
440: DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
441: } else {
442: *Nt = bl->Nt[source];
443: *target = bl->target[source];
444: *size = bl->size[source];
445: *cone = bl->cone[source];
446: *ornt = bl->ornt[source];
447: }
448: return(0);
449: }
451: static PetscErrorCode DMPlexTransformMapCoordinates_BL(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
452: {
453: DMPlexRefine_BL *bl = (DMPlexRefine_BL *) tr->data;
454: PetscInt d;
455: PetscErrorCode ierr;
458: switch (pct) {
459: case DM_POLYTOPE_POINT_PRISM_TENSOR:
460: if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for target point type %s", DMPolytopeTypes[ct]);
461: if (Nv != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of parent vertices %D != 2", Nv);
462: if (r >= bl->n || r < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Invalid replica %D, must be in [0, %D)", r, bl->n);
463: for (d = 0; d < dE; ++d) out[d] = in[d] + bl->h[r] * (in[d + dE] - in[d]);
464: break;
465: default: DMPlexTransformMapCoordinatesBarycenter_Internal(tr, pct, ct, r, Nv, dE, in, out);
466: }
467: return(0);
468: }
470: static PetscErrorCode DMPlexTransformSetFromOptions_BL(PetscOptionItems *PetscOptionsObject, DMPlexTransform tr)
471: {
472: DMPlexRefine_BL *bl = (DMPlexRefine_BL *) tr->data;
473: PetscInt cells[256], n = 256, i;
474: PetscBool flg;
475: PetscErrorCode ierr;
479: PetscOptionsHead(PetscOptionsObject,"DMPlex 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: PetscOptionsTail();
492: return(0);
493: }
495: static PetscErrorCode DMPlexTransformView_BL(DMPlexTransform tr, PetscViewer viewer)
496: {
497: PetscBool isascii;
503: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
504: if (isascii) {
505: PetscViewerFormat format;
506: const char *name;
508: PetscObjectGetName((PetscObject) tr, &name);
509: PetscViewerASCIIPrintf(viewer, "Boundary Layer refinement %s\n", name ? name : "");
510: PetscViewerGetFormat(viewer, &format);
511: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
512: DMLabelView(tr->trType, viewer);
513: }
514: } else {
515: SETERRQ1(PetscObjectComm((PetscObject) tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject) viewer)->type_name);
516: }
517: return(0);
518: }
520: static PetscErrorCode DMPlexTransformDestroy_BL(DMPlexTransform tr)
521: {
522: DMPlexRefine_BL *bl = (DMPlexRefine_BL *) tr->data;
523: PetscInt ict;
524: PetscErrorCode ierr;
527: for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
528: PetscFree4(bl->target[ict], bl->size[ict], bl->cone[ict], bl->ornt[ict]);
529: }
530: PetscFree5(bl->Nt, bl->target, bl->size, bl->cone, bl->ornt);
531: PetscFree(bl->h);
532: PetscFree(tr->data);
533: return(0);
534: }
536: static PetscErrorCode DMPlexTransformInitialize_BL(DMPlexTransform tr)
537: {
539: tr->ops->view = DMPlexTransformView_BL;
540: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_BL;
541: tr->ops->setup = DMPlexTransformSetUp_BL;
542: tr->ops->destroy = DMPlexTransformDestroy_BL;
543: tr->ops->celltransform = DMPlexTransformCellTransform_BL;
544: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_BL;
545: tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_BL;
546: return(0);
547: }
549: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform tr)
550: {
551: DMPlexRefine_BL *bl;
552: PetscErrorCode ierr;
556: PetscNewLog(tr, &bl);
557: tr->data = bl;
559: bl->n = 1; /* 1 split -> 2 new cells */
560: bl->r = 1; /* linear coordinate progression */
562: DMPlexTransformInitialize_BL(tr);
563: return(0);
564: }