Actual source code: plexrefbl.c
1: #include <petsc/private/dmplextransformimpl.h>
3: /*
4: When using the active label, refine types are as follows:
6: Two types are allocated for
7: DM_POLYTOPE_POINT_PRISM_TENSOR
8: DM_POLYTOPE_SEG_PRISM_TENSOR
9: DM_POLYTOPE_TRI_PRISM_TENSOR
10: DM_POLYTOPE_QUAD_PRISM_TENSOR
11: the first is identity, and the second refines.
12: All other types are identity. All the refine types are numbered in order.
13: */
14: static PetscErrorCode DMPlexTransformSetUp_BL(DMPlexTransform tr)
15: {
16: DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
17: const PetscInt n = bl->n;
18: DMPolytopeType ct;
19: DM dm;
20: DMLabel active;
21: PetscInt Nc, No, coff, i, ict;
23: PetscFunctionBegin;
24: /* If no label is given, split all tensor cells */
25: PetscCall(DMPlexTransformGetDM(tr, &dm));
26: PetscCall(DMPlexTransformGetActive(tr, &active));
27: if (active) {
28: IS refineIS;
29: const PetscInt *refineCells;
30: PetscInt pStart, pEnd, p, c;
32: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
33: PetscCall(DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS));
34: PetscCall(DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc));
35: if (refineIS) PetscCall(ISGetIndices(refineIS, &refineCells));
36: for (c = 0; c < Nc; ++c) {
37: const PetscInt cell = refineCells[c];
38: PetscInt *closure = NULL;
39: PetscInt Ncl, cl;
41: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
42: for (cl = 0; cl < Ncl * 2; cl += 2) {
43: const PetscInt point = closure[cl];
44: PetscInt val;
46: PetscCall(DMPlexGetCellType(dm, point, &ct));
47: val = (PetscInt)ct;
48: if (ct > DM_POLYTOPE_POINT_PRISM_TENSOR) ++val;
49: if (ct > DM_POLYTOPE_SEG_PRISM_TENSOR) ++val;
50: if (ct > DM_POLYTOPE_TRI_PRISM_TENSOR) ++val;
51: if (ct > DM_POLYTOPE_QUAD_PRISM_TENSOR) ++val;
52: switch (ct) {
53: case DM_POLYTOPE_POINT_PRISM_TENSOR:
54: case DM_POLYTOPE_SEG_PRISM_TENSOR:
55: case DM_POLYTOPE_TRI_PRISM_TENSOR:
56: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
57: PetscCall(DMLabelSetValue(tr->trType, point, val + 1));
58: break;
59: default:
60: break;
61: }
62: }
63: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
64: }
65: if (refineIS) PetscCall(ISRestoreIndices(refineIS, &refineCells));
66: PetscCall(ISDestroy(&refineIS));
67: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
68: for (p = pStart; p < pEnd; ++p) {
69: PetscInt val;
71: PetscCall(DMLabelGetValue(tr->trType, p, &val));
72: if (val < 0) {
73: PetscCall(DMPlexGetCellType(dm, p, &ct));
74: val = (PetscInt)ct;
75: if (ct > DM_POLYTOPE_POINT_PRISM_TENSOR) ++val;
76: if (ct > DM_POLYTOPE_SEG_PRISM_TENSOR) ++val;
77: if (ct > DM_POLYTOPE_TRI_PRISM_TENSOR) ++val;
78: if (ct > DM_POLYTOPE_QUAD_PRISM_TENSOR) ++val;
79: PetscCall(DMLabelSetValue(tr->trType, p, val));
80: }
81: }
82: }
83: /* Cell heights */
84: PetscCall(PetscMalloc1(n, &bl->h));
85: if (bl->r > 1.) {
86: /* initial height h_0 = (r - 1) / (r^{n+1} - 1)
87: cell height h_i = r^i h_0
88: 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
89: */
90: PetscReal d = (bl->r - 1.) / (PetscPowRealInt(bl->r, n + 1) - 1.);
92: bl->h[0] = d;
93: for (i = 1; i < n; ++i) {
94: d *= bl->r;
95: bl->h[i] = bl->h[i - 1] + d;
96: }
97: } else {
98: /* equal division */
99: for (i = 0; i < n; ++i) bl->h[i] = (i + 1.) / (n + 1);
100: }
102: PetscCall(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));
103: for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) {
104: bl->Nt[ict] = -1;
105: bl->target[ict] = NULL;
106: bl->size[ict] = NULL;
107: bl->cone[ict] = NULL;
108: bl->ornt[ict] = NULL;
109: }
110: /* DM_POLYTOPE_POINT_PRISM_TENSOR produces n points and n+1 tensor segments */
111: ct = DM_POLYTOPE_POINT_PRISM_TENSOR;
112: bl->Nt[ct] = 2;
113: Nc = 7 * 2 + 6 * (n - 1);
114: No = 2 * (n + 1);
115: PetscCall(PetscMalloc4(bl->Nt[ct], &bl->target[ct], bl->Nt[ct], &bl->size[ct], Nc, &bl->cone[ct], No, &bl->ornt[ct]));
116: bl->target[ct][0] = DM_POLYTOPE_POINT;
117: bl->target[ct][1] = DM_POLYTOPE_POINT_PRISM_TENSOR;
118: bl->size[ct][0] = n;
119: bl->size[ct][1] = n + 1;
120: /* cones for tensor segments */
121: bl->cone[ct][0] = DM_POLYTOPE_POINT;
122: bl->cone[ct][1] = 1;
123: bl->cone[ct][2] = 0;
124: bl->cone[ct][3] = 0;
125: bl->cone[ct][4] = DM_POLYTOPE_POINT;
126: bl->cone[ct][5] = 0;
127: bl->cone[ct][6] = 0;
128: for (i = 0; i < n - 1; ++i) {
129: bl->cone[ct][7 + 6 * i + 0] = DM_POLYTOPE_POINT;
130: bl->cone[ct][7 + 6 * i + 1] = 0;
131: bl->cone[ct][7 + 6 * i + 2] = i;
132: bl->cone[ct][7 + 6 * i + 3] = DM_POLYTOPE_POINT;
133: bl->cone[ct][7 + 6 * i + 4] = 0;
134: bl->cone[ct][7 + 6 * i + 5] = i + 1;
135: }
136: bl->cone[ct][7 + 6 * (n - 1) + 0] = DM_POLYTOPE_POINT;
137: bl->cone[ct][7 + 6 * (n - 1) + 1] = 0;
138: bl->cone[ct][7 + 6 * (n - 1) + 2] = n - 1;
139: bl->cone[ct][7 + 6 * (n - 1) + 3] = DM_POLYTOPE_POINT;
140: bl->cone[ct][7 + 6 * (n - 1) + 4] = 1;
141: bl->cone[ct][7 + 6 * (n - 1) + 5] = 1;
142: bl->cone[ct][7 + 6 * (n - 1) + 6] = 0;
143: for (i = 0; i < No; i++) bl->ornt[ct][i] = 0;
144: /* DM_POLYTOPE_SEG_PRISM_TENSOR produces n segments and n+1 tensor quads */
145: ct = DM_POLYTOPE_SEG_PRISM_TENSOR;
146: bl->Nt[ct] = 2;
147: Nc = 8 * n + 15 * 2 + 14 * (n - 1);
148: No = 2 * n + 4 * (n + 1);
149: PetscCall(PetscMalloc4(bl->Nt[ct], &bl->target[ct], bl->Nt[ct], &bl->size[ct], Nc, &bl->cone[ct], No, &bl->ornt[ct]));
150: bl->target[ct][0] = DM_POLYTOPE_SEGMENT;
151: bl->target[ct][1] = DM_POLYTOPE_SEG_PRISM_TENSOR;
152: bl->size[ct][0] = n;
153: bl->size[ct][1] = n + 1;
154: /* cones for segments */
155: for (i = 0; i < n; ++i) {
156: bl->cone[ct][8 * i + 0] = DM_POLYTOPE_POINT;
157: bl->cone[ct][8 * i + 1] = 1;
158: bl->cone[ct][8 * i + 2] = 2;
159: bl->cone[ct][8 * i + 3] = i;
160: bl->cone[ct][8 * i + 4] = DM_POLYTOPE_POINT;
161: bl->cone[ct][8 * i + 5] = 1;
162: bl->cone[ct][8 * i + 6] = 3;
163: bl->cone[ct][8 * i + 7] = i;
164: }
165: /* cones for tensor quads */
166: coff = 8 * n;
167: bl->cone[ct][coff + 0] = DM_POLYTOPE_SEGMENT;
168: bl->cone[ct][coff + 1] = 1;
169: bl->cone[ct][coff + 2] = 0;
170: bl->cone[ct][coff + 3] = 0;
171: bl->cone[ct][coff + 4] = DM_POLYTOPE_SEGMENT;
172: bl->cone[ct][coff + 5] = 0;
173: bl->cone[ct][coff + 6] = 0;
174: bl->cone[ct][coff + 7] = DM_POLYTOPE_POINT_PRISM_TENSOR;
175: bl->cone[ct][coff + 8] = 1;
176: bl->cone[ct][coff + 9] = 2;
177: bl->cone[ct][coff + 10] = 0;
178: bl->cone[ct][coff + 11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
179: bl->cone[ct][coff + 12] = 1;
180: bl->cone[ct][coff + 13] = 3;
181: bl->cone[ct][coff + 14] = 0;
182: for (i = 0; i < n - 1; ++i) {
183: bl->cone[ct][coff + 15 + 14 * i + 0] = DM_POLYTOPE_SEGMENT;
184: bl->cone[ct][coff + 15 + 14 * i + 1] = 0;
185: bl->cone[ct][coff + 15 + 14 * i + 2] = i;
186: bl->cone[ct][coff + 15 + 14 * i + 3] = DM_POLYTOPE_SEGMENT;
187: bl->cone[ct][coff + 15 + 14 * i + 4] = 0;
188: bl->cone[ct][coff + 15 + 14 * i + 5] = i + 1;
189: bl->cone[ct][coff + 15 + 14 * i + 6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
190: bl->cone[ct][coff + 15 + 14 * i + 7] = 1;
191: bl->cone[ct][coff + 15 + 14 * i + 8] = 2;
192: bl->cone[ct][coff + 15 + 14 * i + 9] = i + 1;
193: bl->cone[ct][coff + 15 + 14 * i + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
194: bl->cone[ct][coff + 15 + 14 * i + 11] = 1;
195: bl->cone[ct][coff + 15 + 14 * i + 12] = 3;
196: bl->cone[ct][coff + 15 + 14 * i + 13] = i + 1;
197: }
198: bl->cone[ct][coff + 15 + 14 * (n - 1) + 0] = DM_POLYTOPE_SEGMENT;
199: bl->cone[ct][coff + 15 + 14 * (n - 1) + 1] = 0;
200: bl->cone[ct][coff + 15 + 14 * (n - 1) + 2] = n - 1;
201: bl->cone[ct][coff + 15 + 14 * (n - 1) + 3] = DM_POLYTOPE_SEGMENT;
202: bl->cone[ct][coff + 15 + 14 * (n - 1) + 4] = 1;
203: bl->cone[ct][coff + 15 + 14 * (n - 1) + 5] = 1;
204: bl->cone[ct][coff + 15 + 14 * (n - 1) + 6] = 0;
205: bl->cone[ct][coff + 15 + 14 * (n - 1) + 7] = DM_POLYTOPE_POINT_PRISM_TENSOR;
206: bl->cone[ct][coff + 15 + 14 * (n - 1) + 8] = 1;
207: bl->cone[ct][coff + 15 + 14 * (n - 1) + 9] = 2;
208: bl->cone[ct][coff + 15 + 14 * (n - 1) + 10] = n;
209: bl->cone[ct][coff + 15 + 14 * (n - 1) + 11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
210: bl->cone[ct][coff + 15 + 14 * (n - 1) + 12] = 1;
211: bl->cone[ct][coff + 15 + 14 * (n - 1) + 13] = 3;
212: bl->cone[ct][coff + 15 + 14 * (n - 1) + 14] = n;
213: for (i = 0; i < No; i++) bl->ornt[ct][i] = 0;
214: /* DM_POLYTOPE_TRI_PRISM_TENSOR produces n triangles and n+1 tensor triangular prisms */
215: ct = DM_POLYTOPE_TRI_PRISM_TENSOR;
216: bl->Nt[ct] = 2;
217: Nc = 12 * n + 19 * 2 + 18 * (n - 1);
218: No = 3 * n + 5 * (n + 1);
219: PetscCall(PetscMalloc4(bl->Nt[ct], &bl->target[ct], bl->Nt[ct], &bl->size[ct], Nc, &bl->cone[ct], No, &bl->ornt[ct]));
220: bl->target[ct][0] = DM_POLYTOPE_TRIANGLE;
221: bl->target[ct][1] = DM_POLYTOPE_TRI_PRISM_TENSOR;
222: bl->size[ct][0] = n;
223: bl->size[ct][1] = n + 1;
224: /* cones for triangles */
225: for (i = 0; i < n; ++i) {
226: bl->cone[ct][12 * i + 0] = DM_POLYTOPE_SEGMENT;
227: bl->cone[ct][12 * i + 1] = 1;
228: bl->cone[ct][12 * i + 2] = 2;
229: bl->cone[ct][12 * i + 3] = i;
230: bl->cone[ct][12 * i + 4] = DM_POLYTOPE_SEGMENT;
231: bl->cone[ct][12 * i + 5] = 1;
232: bl->cone[ct][12 * i + 6] = 3;
233: bl->cone[ct][12 * i + 7] = i;
234: bl->cone[ct][12 * i + 8] = DM_POLYTOPE_SEGMENT;
235: bl->cone[ct][12 * i + 9] = 1;
236: bl->cone[ct][12 * i + 10] = 4;
237: bl->cone[ct][12 * i + 11] = i;
238: }
239: /* cones for triangular prisms */
240: coff = 12 * n;
241: bl->cone[ct][coff + 0] = DM_POLYTOPE_TRIANGLE;
242: bl->cone[ct][coff + 1] = 1;
243: bl->cone[ct][coff + 2] = 0;
244: bl->cone[ct][coff + 3] = 0;
245: bl->cone[ct][coff + 4] = DM_POLYTOPE_TRIANGLE;
246: bl->cone[ct][coff + 5] = 0;
247: bl->cone[ct][coff + 6] = 0;
248: bl->cone[ct][coff + 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
249: bl->cone[ct][coff + 8] = 1;
250: bl->cone[ct][coff + 9] = 2;
251: bl->cone[ct][coff + 10] = 0;
252: bl->cone[ct][coff + 11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
253: bl->cone[ct][coff + 12] = 1;
254: bl->cone[ct][coff + 13] = 3;
255: bl->cone[ct][coff + 14] = 0;
256: bl->cone[ct][coff + 15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
257: bl->cone[ct][coff + 16] = 1;
258: bl->cone[ct][coff + 17] = 4;
259: bl->cone[ct][coff + 18] = 0;
260: for (i = 0; i < n - 1; ++i) {
261: bl->cone[ct][coff + 19 + 18 * i + 0] = DM_POLYTOPE_TRIANGLE;
262: bl->cone[ct][coff + 19 + 18 * i + 1] = 0;
263: bl->cone[ct][coff + 19 + 18 * i + 2] = i;
264: bl->cone[ct][coff + 19 + 18 * i + 3] = DM_POLYTOPE_TRIANGLE;
265: bl->cone[ct][coff + 19 + 18 * i + 4] = 0;
266: bl->cone[ct][coff + 19 + 18 * i + 5] = i + 1;
267: bl->cone[ct][coff + 19 + 18 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
268: bl->cone[ct][coff + 19 + 18 * i + 7] = 1;
269: bl->cone[ct][coff + 19 + 18 * i + 8] = 2;
270: bl->cone[ct][coff + 19 + 18 * i + 9] = i + 1;
271: bl->cone[ct][coff + 19 + 18 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
272: bl->cone[ct][coff + 19 + 18 * i + 11] = 1;
273: bl->cone[ct][coff + 19 + 18 * i + 12] = 3;
274: bl->cone[ct][coff + 19 + 18 * i + 13] = i + 1;
275: bl->cone[ct][coff + 19 + 18 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
276: bl->cone[ct][coff + 19 + 18 * i + 15] = 1;
277: bl->cone[ct][coff + 19 + 18 * i + 16] = 4;
278: bl->cone[ct][coff + 19 + 18 * i + 17] = i + 1;
279: }
280: bl->cone[ct][coff + 19 + 18 * (n - 1) + 0] = DM_POLYTOPE_TRIANGLE;
281: bl->cone[ct][coff + 19 + 18 * (n - 1) + 1] = 0;
282: bl->cone[ct][coff + 19 + 18 * (n - 1) + 2] = n - 1;
283: bl->cone[ct][coff + 19 + 18 * (n - 1) + 3] = DM_POLYTOPE_TRIANGLE;
284: bl->cone[ct][coff + 19 + 18 * (n - 1) + 4] = 1;
285: bl->cone[ct][coff + 19 + 18 * (n - 1) + 5] = 1;
286: bl->cone[ct][coff + 19 + 18 * (n - 1) + 6] = 0;
287: bl->cone[ct][coff + 19 + 18 * (n - 1) + 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
288: bl->cone[ct][coff + 19 + 18 * (n - 1) + 8] = 1;
289: bl->cone[ct][coff + 19 + 18 * (n - 1) + 9] = 2;
290: bl->cone[ct][coff + 19 + 18 * (n - 1) + 10] = n;
291: bl->cone[ct][coff + 19 + 18 * (n - 1) + 11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
292: bl->cone[ct][coff + 19 + 18 * (n - 1) + 12] = 1;
293: bl->cone[ct][coff + 19 + 18 * (n - 1) + 13] = 3;
294: bl->cone[ct][coff + 19 + 18 * (n - 1) + 14] = n;
295: bl->cone[ct][coff + 19 + 18 * (n - 1) + 15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
296: bl->cone[ct][coff + 19 + 18 * (n - 1) + 16] = 1;
297: bl->cone[ct][coff + 19 + 18 * (n - 1) + 17] = 4;
298: bl->cone[ct][coff + 19 + 18 * (n - 1) + 18] = n;
299: for (i = 0; i < No; ++i) bl->ornt[ct][i] = 0;
300: /* DM_POLYTOPE_QUAD_PRISM_TENSOR produces n quads and n+1 tensor quad prisms */
301: ct = DM_POLYTOPE_QUAD_PRISM_TENSOR;
302: bl->Nt[ct] = 2;
303: Nc = 16 * n + 23 * 2 + 22 * (n - 1);
304: No = 4 * n + 6 * (n + 1);
305: PetscCall(PetscMalloc4(bl->Nt[ct], &bl->target[ct], bl->Nt[ct], &bl->size[ct], Nc, &bl->cone[ct], No, &bl->ornt[ct]));
306: bl->target[ct][0] = DM_POLYTOPE_QUADRILATERAL;
307: bl->target[ct][1] = DM_POLYTOPE_QUAD_PRISM_TENSOR;
308: bl->size[ct][0] = n;
309: bl->size[ct][1] = n + 1;
310: /* cones for quads */
311: for (i = 0; i < n; ++i) {
312: bl->cone[ct][16 * i + 0] = DM_POLYTOPE_SEGMENT;
313: bl->cone[ct][16 * i + 1] = 1;
314: bl->cone[ct][16 * i + 2] = 2;
315: bl->cone[ct][16 * i + 3] = i;
316: bl->cone[ct][16 * i + 4] = DM_POLYTOPE_SEGMENT;
317: bl->cone[ct][16 * i + 5] = 1;
318: bl->cone[ct][16 * i + 6] = 3;
319: bl->cone[ct][16 * i + 7] = i;
320: bl->cone[ct][16 * i + 8] = DM_POLYTOPE_SEGMENT;
321: bl->cone[ct][16 * i + 9] = 1;
322: bl->cone[ct][16 * i + 10] = 4;
323: bl->cone[ct][16 * i + 11] = i;
324: bl->cone[ct][16 * i + 12] = DM_POLYTOPE_SEGMENT;
325: bl->cone[ct][16 * i + 13] = 1;
326: bl->cone[ct][16 * i + 14] = 5;
327: bl->cone[ct][16 * i + 15] = i;
328: }
329: /* cones for quad prisms */
330: coff = 16 * n;
331: bl->cone[ct][coff + 0] = DM_POLYTOPE_QUADRILATERAL;
332: bl->cone[ct][coff + 1] = 1;
333: bl->cone[ct][coff + 2] = 0;
334: bl->cone[ct][coff + 3] = 0;
335: bl->cone[ct][coff + 4] = DM_POLYTOPE_QUADRILATERAL;
336: bl->cone[ct][coff + 5] = 0;
337: bl->cone[ct][coff + 6] = 0;
338: bl->cone[ct][coff + 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
339: bl->cone[ct][coff + 8] = 1;
340: bl->cone[ct][coff + 9] = 2;
341: bl->cone[ct][coff + 10] = 0;
342: bl->cone[ct][coff + 11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
343: bl->cone[ct][coff + 12] = 1;
344: bl->cone[ct][coff + 13] = 3;
345: bl->cone[ct][coff + 14] = 0;
346: bl->cone[ct][coff + 15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
347: bl->cone[ct][coff + 16] = 1;
348: bl->cone[ct][coff + 17] = 4;
349: bl->cone[ct][coff + 18] = 0;
350: bl->cone[ct][coff + 19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
351: bl->cone[ct][coff + 20] = 1;
352: bl->cone[ct][coff + 21] = 5;
353: bl->cone[ct][coff + 22] = 0;
354: for (i = 0; i < n - 1; ++i) {
355: bl->cone[ct][coff + 23 + 22 * i + 0] = DM_POLYTOPE_QUADRILATERAL;
356: bl->cone[ct][coff + 23 + 22 * i + 1] = 0;
357: bl->cone[ct][coff + 23 + 22 * i + 2] = i;
358: bl->cone[ct][coff + 23 + 22 * i + 3] = DM_POLYTOPE_QUADRILATERAL;
359: bl->cone[ct][coff + 23 + 22 * i + 4] = 0;
360: bl->cone[ct][coff + 23 + 22 * i + 5] = i + 1;
361: bl->cone[ct][coff + 23 + 22 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
362: bl->cone[ct][coff + 23 + 22 * i + 7] = 1;
363: bl->cone[ct][coff + 23 + 22 * i + 8] = 2;
364: bl->cone[ct][coff + 23 + 22 * i + 9] = i + 1;
365: bl->cone[ct][coff + 23 + 22 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
366: bl->cone[ct][coff + 23 + 22 * i + 11] = 1;
367: bl->cone[ct][coff + 23 + 22 * i + 12] = 3;
368: bl->cone[ct][coff + 23 + 22 * i + 13] = i + 1;
369: bl->cone[ct][coff + 23 + 22 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
370: bl->cone[ct][coff + 23 + 22 * i + 15] = 1;
371: bl->cone[ct][coff + 23 + 22 * i + 16] = 4;
372: bl->cone[ct][coff + 23 + 22 * i + 17] = i + 1;
373: bl->cone[ct][coff + 23 + 22 * i + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
374: bl->cone[ct][coff + 23 + 22 * i + 19] = 1;
375: bl->cone[ct][coff + 23 + 22 * i + 20] = 5;
376: bl->cone[ct][coff + 23 + 22 * i + 21] = i + 1;
377: }
378: bl->cone[ct][coff + 23 + 22 * (n - 1) + 0] = DM_POLYTOPE_QUADRILATERAL;
379: bl->cone[ct][coff + 23 + 22 * (n - 1) + 1] = 0;
380: bl->cone[ct][coff + 23 + 22 * (n - 1) + 2] = n - 1;
381: bl->cone[ct][coff + 23 + 22 * (n - 1) + 3] = DM_POLYTOPE_QUADRILATERAL;
382: bl->cone[ct][coff + 23 + 22 * (n - 1) + 4] = 1;
383: bl->cone[ct][coff + 23 + 22 * (n - 1) + 5] = 1;
384: bl->cone[ct][coff + 23 + 22 * (n - 1) + 6] = 0;
385: bl->cone[ct][coff + 23 + 22 * (n - 1) + 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
386: bl->cone[ct][coff + 23 + 22 * (n - 1) + 8] = 1;
387: bl->cone[ct][coff + 23 + 22 * (n - 1) + 9] = 2;
388: bl->cone[ct][coff + 23 + 22 * (n - 1) + 10] = n;
389: bl->cone[ct][coff + 23 + 22 * (n - 1) + 11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
390: bl->cone[ct][coff + 23 + 22 * (n - 1) + 12] = 1;
391: bl->cone[ct][coff + 23 + 22 * (n - 1) + 13] = 3;
392: bl->cone[ct][coff + 23 + 22 * (n - 1) + 14] = n;
393: bl->cone[ct][coff + 23 + 22 * (n - 1) + 15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
394: bl->cone[ct][coff + 23 + 22 * (n - 1) + 16] = 1;
395: bl->cone[ct][coff + 23 + 22 * (n - 1) + 17] = 4;
396: bl->cone[ct][coff + 23 + 22 * (n - 1) + 18] = n;
397: bl->cone[ct][coff + 23 + 22 * (n - 1) + 19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
398: bl->cone[ct][coff + 23 + 22 * (n - 1) + 20] = 1;
399: bl->cone[ct][coff + 23 + 22 * (n - 1) + 21] = 5;
400: bl->cone[ct][coff + 23 + 22 * (n - 1) + 22] = n;
401: for (i = 0; i < No; ++i) bl->ornt[ct][i] = 0;
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: static PetscErrorCode DMPlexTransformGetSubcellOrientation_BL(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
406: {
407: DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
408: const PetscInt n = bl->n;
409: PetscInt tquad_tquad_o[] = {0, 1, -2, -1, 1, 0, -1, -2, -2, -1, 0, 1, -1, -2, 1, 0};
411: PetscFunctionBeginHot;
412: *rnew = r;
413: *onew = o;
414: if (tr->trType) {
415: PetscInt rt;
417: PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
418: if ((rt != DM_POLYTOPE_POINT_PRISM_TENSOR + 1) && (rt != DM_POLYTOPE_SEG_PRISM_TENSOR + 2) && (rt != DM_POLYTOPE_TRI_PRISM_TENSOR + 3) && (rt != DM_POLYTOPE_QUAD_PRISM_TENSOR + 4)) {
419: PetscCall(DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew));
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
422: }
423: switch (sct) {
424: case DM_POLYTOPE_POINT_PRISM_TENSOR:
425: switch (tct) {
426: case DM_POLYTOPE_POINT:
427: *rnew = !so ? r : n - 1 - r;
428: break;
429: case DM_POLYTOPE_POINT_PRISM_TENSOR:
430: if (!so) {
431: *rnew = r;
432: *onew = o;
433: } else {
434: *rnew = n - r;
435: *onew = -(o + 1);
436: }
437: default:
438: break;
439: }
440: break;
441: case DM_POLYTOPE_SEG_PRISM_TENSOR:
442: switch (tct) {
443: case DM_POLYTOPE_SEGMENT:
444: *onew = (so == 0) || (so == 1) ? o : -(o + 1);
445: *rnew = (so == 0) || (so == -1) ? r : n - 1 - r;
446: break;
447: case DM_POLYTOPE_SEG_PRISM_TENSOR:
448: *onew = tquad_tquad_o[(so + 2) * 4 + o + 2];
449: *rnew = (so == 0) || (so == -1) ? r : n - r;
450: break;
451: default:
452: break;
453: }
454: break;
455: default:
456: PetscCall(DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew));
457: }
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: static PetscErrorCode DMPlexTransformCellTransform_BL(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
462: {
463: DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
465: PetscFunctionBeginHot;
466: if (rt) *rt = -1;
467: if (tr->trType && p >= 0) {
468: PetscInt val;
470: PetscCall(DMLabelGetValue(tr->trType, p, &val));
471: if (rt) *rt = val;
472: if ((val != DM_POLYTOPE_POINT_PRISM_TENSOR + 1) && (val != DM_POLYTOPE_SEG_PRISM_TENSOR + 2) && (val != DM_POLYTOPE_TRI_PRISM_TENSOR + 3) && (val != DM_POLYTOPE_QUAD_PRISM_TENSOR + 4)) {
473: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
476: }
477: if (bl->Nt[source] < 0) {
478: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
479: } else {
480: *Nt = bl->Nt[source];
481: *target = bl->target[source];
482: *size = bl->size[source];
483: *cone = bl->cone[source];
484: *ornt = bl->ornt[source];
485: }
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: static PetscErrorCode DMPlexTransformMapCoordinates_BL(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
490: {
491: DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
492: PetscInt d;
494: PetscFunctionBeginHot;
495: switch (pct) {
496: case DM_POLYTOPE_POINT_PRISM_TENSOR:
497: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for target point type %s", DMPolytopeTypes[ct]);
498: PetscCheck(Nv == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of parent vertices %" PetscInt_FMT " != 2", Nv);
499: PetscCheck(r < bl->n && r >= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Invalid replica %" PetscInt_FMT ", must be in [0, %" PetscInt_FMT ")", r, bl->n);
500: for (d = 0; d < dE; ++d) out[d] = in[d] + bl->h[r] * (in[d + dE] - in[d]);
501: break;
502: default:
503: PetscCall(DMPlexTransformMapCoordinatesBarycenter_Internal(tr, pct, ct, p, r, Nv, dE, in, out));
504: }
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode DMPlexTransformSetFromOptions_BL(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
509: {
510: DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
511: PetscInt cells[256], n = 256, i;
512: PetscBool flg;
514: PetscFunctionBegin;
515: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Boundary Layer Options");
516: PetscCall(PetscOptionsBoundedInt("-dm_plex_transform_bl_splits", "Number of divisions of a cell", "", bl->n, &bl->n, NULL, 1));
517: PetscCall(PetscOptionsReal("-dm_plex_transform_bl_height_factor", "Factor increase for height at each division", "", bl->r, &bl->r, NULL));
518: PetscCall(PetscOptionsIntArray("-dm_plex_transform_bl_ref_cell", "Mark cells for refinement", "", cells, &n, &flg));
519: if (flg) {
520: DMLabel active;
522: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active));
523: for (i = 0; i < n; ++i) PetscCall(DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE));
524: PetscCall(DMPlexTransformSetActive(tr, active));
525: PetscCall(DMLabelDestroy(&active));
526: }
527: PetscOptionsHeadEnd();
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: static PetscErrorCode DMPlexTransformView_BL(DMPlexTransform tr, PetscViewer viewer)
532: {
533: PetscBool isascii;
535: PetscFunctionBegin;
538: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
539: if (isascii) {
540: PetscViewerFormat format;
541: const char *name;
543: PetscCall(PetscObjectGetName((PetscObject)tr, &name));
544: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary Layer refinement %s\n", name ? name : ""));
545: PetscCall(PetscViewerGetFormat(viewer, &format));
546: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(DMLabelView(tr->trType, viewer));
547: } else {
548: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
549: }
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: static PetscErrorCode DMPlexTransformDestroy_BL(DMPlexTransform tr)
554: {
555: DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
556: PetscInt ict;
558: PetscFunctionBegin;
559: for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) PetscCall(PetscFree4(bl->target[ict], bl->size[ict], bl->cone[ict], bl->ornt[ict]));
560: PetscCall(PetscFree5(bl->Nt, bl->target, bl->size, bl->cone, bl->ornt));
561: PetscCall(PetscFree(bl->h));
562: PetscCall(PetscFree(tr->data));
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: static PetscErrorCode DMPlexTransformInitialize_BL(DMPlexTransform tr)
567: {
568: PetscFunctionBegin;
569: tr->ops->view = DMPlexTransformView_BL;
570: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_BL;
571: tr->ops->setup = DMPlexTransformSetUp_BL;
572: tr->ops->destroy = DMPlexTransformDestroy_BL;
573: tr->ops->setdimensions = DMPlexTransformSetDimensions_Internal;
574: tr->ops->celltransform = DMPlexTransformCellTransform_BL;
575: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_BL;
576: tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_BL;
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform tr)
581: {
582: DMPlexRefine_BL *bl;
584: PetscFunctionBegin;
586: PetscCall(PetscNew(&bl));
587: tr->data = bl;
589: bl->n = 1; /* 1 split -> 2 new cells */
590: bl->r = 1; /* linear coordinate progression */
592: PetscCall(DMPlexTransformInitialize_BL(tr));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }