xref: /petsc/src/dm/impls/plex/transform/impls/refine/bl/plexrefbl.c (revision 6d8694c4fbab79f9439f1ad13c0386ba7ee1ca4b) 
1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/
2 
3 /*
4   When using the active label, refine types are as follows:
5 
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 */
DMPlexTransformSetUp_BL(DMPlexTransform tr)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;
22 
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;
31 
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;
40 
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;
45 
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;
70 
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.);
91 
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   }
101 
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 }
404 
DMPlexTransformGetSubcellOrientation_BL(DMPlexTransform tr,DMPolytopeType sct,PetscInt sp,PetscInt so,DMPolytopeType tct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)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};
410 
411   PetscFunctionBeginHot;
412   *rnew = r;
413   *onew = o;
414   if (tr->trType) {
415     PetscInt rt;
416 
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 }
460 
DMPlexTransformCellTransform_BL(DMPlexTransform tr,DMPolytopeType source,PetscInt p,PetscInt * rt,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])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;
464 
465   PetscFunctionBeginHot;
466   if (rt) *rt = -1;
467   if (tr->trType && p >= 0) {
468     PetscInt val;
469 
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 }
488 
DMPlexTransformMapCoordinates_BL(DMPlexTransform tr,DMPolytopeType pct,DMPolytopeType ct,PetscInt p,PetscInt r,PetscInt Nv,PetscInt dE,const PetscScalar in[],PetscScalar out[])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;
493 
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 }
507 
DMPlexTransformSetFromOptions_BL(DMPlexTransform tr,PetscOptionItems PetscOptionsObject)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;
513 
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;
521 
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 }
530 
DMPlexTransformView_BL(DMPlexTransform tr,PetscViewer viewer)531 static PetscErrorCode DMPlexTransformView_BL(DMPlexTransform tr, PetscViewer viewer)
532 {
533   PetscBool isascii;
534 
535   PetscFunctionBegin;
536   PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
537   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
538   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
539   if (isascii) {
540     PetscViewerFormat format;
541     const char       *name;
542 
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 }
552 
DMPlexTransformDestroy_BL(DMPlexTransform tr)553 static PetscErrorCode DMPlexTransformDestroy_BL(DMPlexTransform tr)
554 {
555   DMPlexRefine_BL *bl = (DMPlexRefine_BL *)tr->data;
556   PetscInt         ict;
557 
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 }
565 
DMPlexTransformInitialize_BL(DMPlexTransform tr)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 }
579 
DMPlexTransformCreate_BL(DMPlexTransform tr)580 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform tr)
581 {
582   DMPlexRefine_BL *bl;
583 
584   PetscFunctionBegin;
585   PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
586   PetscCall(PetscNew(&bl));
587   tr->data = bl;
588 
589   bl->n = 1; /* 1 split -> 2 new cells */
590   bl->r = 1; /* linear coordinate progression */
591 
592   PetscCall(DMPlexTransformInitialize_BL(tr));
593   PetscFunctionReturn(PETSC_SUCCESS);
594 }
595