Actual source code: plexrefsbr.c
1: #include <petsc/private/dmplextransformimpl.h>
2: #include <petscsf.h>
4: PetscBool SBRcite = PETSC_FALSE;
5: const char SBRCitation[] = "@article{PlazaCarey2000,\n"
6: " title = {Local refinement of simplicial grids based on the skeleton},\n"
7: " journal = {Applied Numerical Mathematics},\n"
8: " author = {A. Plaza and Graham F. Carey},\n"
9: " volume = {32},\n"
10: " number = {3},\n"
11: " pages = {195--218},\n"
12: " doi = {10.1016/S0168-9274(99)00022-7},\n"
13: " year = {2000}\n}\n";
15: /*
16: Tetrahedron combinatorics, in the canonical local numbering used throughout this file (matching
17: the convention documented in DMPlexTransformCellRefine_Regular() for DM_POLYTOPE_TETRAHEDRON):
19: e0 = (v0,v1), e1 = (v1,v2), e2 = (v2,v0), e3 = (v0,v3), e4 = (v1,v3), e5 = (v2,v3)
20: f0 = (v0,v1,v2), f1 = (v0,v3,v1), f2 = (v0,v2,v3), f3 = (v2,v1,v3)
22: tetEdgeVert[e] - the 2 vertices of edge e, in its canonical direction
23: tetFaceVert[f] - the 3 vertices of face f, in its canonical order
24: tetFaceEdge[f] - the edges of face f, at local indices 0,1,2 (i.e. (v_0,v_1),(v_1,v_2),(v_2,v_0) of that face)
25: tetEdgeFaceLoc[e] - a (face, local edge index) path that reaches edge e in its canonical direction
26: tetVertPath[v] - a (face, local edge index, local vertex index) path that reaches vertex v
27: tetOppFace[x] - the 3 local vertex positions (not values) of the face opposite position x, i.e. f[3-x]
28: */
29: static const PetscInt tetEdgeVert[6][2] = {
30: {0, 1},
31: {1, 2},
32: {2, 0},
33: {0, 3},
34: {1, 3},
35: {2, 3}
36: };
37: static const PetscInt tetFaceVert[4][3] = {
38: {0, 1, 2},
39: {0, 3, 1},
40: {0, 2, 3},
41: {2, 1, 3}
42: };
43: static const PetscInt tetFaceEdge[4][3] = {
44: {0, 1, 2},
45: {3, 4, 0},
46: {2, 5, 3},
47: {1, 4, 5}
48: };
49: static const PetscInt tetEdgeFaceLoc[6][2] = {
50: {0, 0},
51: {0, 1},
52: {0, 2},
53: {1, 0},
54: {3, 1},
55: {2, 1}
56: };
57: static const PetscInt tetVertPath[4][3] = {
58: {0, 0, 0},
59: {0, 0, 1},
60: {0, 2, 0},
61: {1, 0, 1}
62: };
63: static const PetscInt tetOppFace[4][3] = {
64: {2, 1, 3},
65: {0, 2, 3},
66: {0, 3, 1},
67: {0, 1, 2}
68: };
69: /* tetVertFaces[v], tetEdgeFaces[e] - bitmasks over f0..f3 of the faces containing vertex v or edge e */
70: static const PetscInt tetVertFaces[4] = {7, 11, 13, 14};
71: static const PetscInt tetEdgeFaces[6] = {3, 9, 5, 6, 10, 12};
72: /* The edges of the triangle, triEdgeVert[e] = (v_e, v_{e+1}) */
73: static const PetscInt triEdgeVert[3][2] = {
74: {0, 1},
75: {1, 2},
76: {2, 0}
77: };
79: /* Find the tetrahedron edge index (0-5) connecting the two given (unordered) vertices */
80: static PetscInt SBREdgeIndexFromVerts_Private(PetscInt v0, PetscInt v1)
81: {
82: for (PetscInt k = 0; k < 6; ++k) {
83: if ((tetEdgeVert[k][0] == v0 && tetEdgeVert[k][1] == v1) || (tetEdgeVert[k][0] == v1 && tetEdgeVert[k][1] == v0)) return k;
84: }
85: return -1;
86: }
88: static PetscErrorCode SBRGetEdgeLen_Private(DMPlexTransform tr, PetscInt edge, PetscReal *len)
89: {
90: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
91: DM dm;
92: PetscInt off;
94: PetscFunctionBeginHot;
95: PetscCall(DMPlexTransformGetDM(tr, &dm));
96: PetscCall(PetscSectionGetOffset(sbr->secEdgeLen, edge, &off));
97: if (sbr->edgeLen[off] <= 0.0) {
98: DM cdm;
99: Vec coordsLocal;
100: const PetscScalar *coords;
101: const PetscInt *cone;
102: PetscScalar *cA, *cB;
103: PetscInt coneSize, cdim;
105: PetscCall(DMGetCoordinateDM(dm, &cdm));
106: PetscCall(DMPlexGetCone(dm, edge, &cone));
107: PetscCall(DMPlexGetConeSize(dm, edge, &coneSize));
108: PetscCheck(coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Edge %" PetscInt_FMT " cone size must be 2, not %" PetscInt_FMT, edge, coneSize);
109: PetscCall(DMGetCoordinateDim(dm, &cdim));
110: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordsLocal));
111: PetscCall(VecGetArrayRead(coordsLocal, &coords));
112: PetscCall(DMPlexPointLocalRead(cdm, cone[0], coords, &cA));
113: PetscCall(DMPlexPointLocalRead(cdm, cone[1], coords, &cB));
114: sbr->edgeLen[off] = DMPlex_DistD_Internal(cdim, cA, cB);
115: PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
116: }
117: *len = sbr->edgeLen[off];
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: static PetscErrorCode SBRGetEdgeMid_Private(DMPlexTransform tr, PetscInt edge, PetscReal mid[3])
122: {
123: DM dm, cdm;
124: Vec coordsLocal;
125: const PetscScalar *coords;
126: const PetscInt *cone;
127: PetscScalar *cA, *cB;
128: PetscInt cdim;
130: PetscFunctionBeginHot;
131: PetscCall(DMPlexTransformGetDM(tr, &dm));
132: PetscCall(DMGetCoordinateDM(dm, &cdm));
133: PetscCall(DMPlexGetCone(dm, edge, &cone));
134: PetscCall(DMGetCoordinateDim(dm, &cdim));
135: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordsLocal));
136: PetscCall(VecGetArrayRead(coordsLocal, &coords));
137: PetscCall(DMPlexPointLocalRead(cdm, cone[0], coords, &cA));
138: PetscCall(DMPlexPointLocalRead(cdm, cone[1], coords, &cB));
139: for (PetscInt d = 0; d < 3; ++d) mid[d] = d < cdim ? 0.5 * PetscRealPart(cA[d] + cB[d]) : 0.0;
140: PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /*
145: Determine whether edge eA is bisected before edge eB: longer edges come first, and exact length
146: ties are broken in favor of the lexicographically greater edge midpoint. The tie-breaker depends
147: only on the geometry, not on mesh point numbers or cone traversal order, so every face and every
148: cell containing the two edges resolves the tie the same way, on every process. This keeps the
149: subdivision each face chooses for itself consistent with the subdivision induced on that face by
150: the bisection of the cells above it (see the discussion at RefinementType).
151: */
152: static PetscErrorCode SBREdgePrecedes_Private(DMPlexTransform tr, PetscInt eA, PetscInt eB, PetscBool *precedes)
153: {
154: PetscReal lenA, lenB, midA[3], midB[3];
156: PetscFunctionBeginHot;
157: PetscCall(SBRGetEdgeLen_Private(tr, eA, &lenA));
158: PetscCall(SBRGetEdgeLen_Private(tr, eB, &lenB));
159: if (lenA != lenB) {
160: *precedes = lenA > lenB ? PETSC_TRUE : PETSC_FALSE;
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
163: PetscCall(SBRGetEdgeMid_Private(tr, eA, midA));
164: PetscCall(SBRGetEdgeMid_Private(tr, eB, midB));
165: *precedes = PETSC_FALSE;
166: for (PetscInt d = 0; d < 3; ++d) {
167: if (midA[d] != midB[d]) {
168: *precedes = midA[d] > midB[d] ? PETSC_TRUE : PETSC_FALSE;
169: break;
170: }
171: }
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: /*
176: Get the 6 edges of a tetrahedron, in the canonical local order used throughout this file:
177: e0 = (v0,v1), e1 = (v1,v2), e2 = (v2,v0), e3 = (v0,v3), e4 = (v1,v3), e5 = (v2,v3), where
178: v0..v3 are the tetrahedron's vertices in the order given by its own transitive closure. This
179: matches the convention documented in DMPlexTransformCellRefine_Regular() for DM_POLYTOPE_TETRAHEDRON.
180: The faces of the tetrahedron, in DMPlexGetCone() order, are f0 = (v0,v1,v2), f1 = (v0,v3,v1),
181: f2 = (v0,v2,v3), f3 = (v2,v1,v3), so each e_k below is read off face fIdx[k] at local edge index
182: lIdx[k], reoriented for that face's actual orientation as seen from the tetrahedron.
183: */
184: static PetscErrorCode SBRGetTetEdges_Private(DM dm, PetscInt tet, PetscInt edges[6])
185: {
186: const PetscInt *fcone, *forient;
188: PetscFunctionBeginHot;
189: PetscCall(DMPlexGetCone(dm, tet, &fcone));
190: PetscCall(DMPlexGetConeOrientation(dm, tet, &forient));
191: for (PetscInt k = 0; k < 6; ++k) {
192: const PetscInt floc = tetEdgeFaceLoc[k][0];
193: const PetscInt face = fcone[floc];
194: const PetscInt *arr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, forient[floc]);
195: const PetscInt li = arr[tetEdgeFaceLoc[k][1] * 2];
196: const PetscInt *econe;
198: PetscCall(DMPlexGetCone(dm, face, &econe));
199: edges[k] = econe[li];
200: }
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /* Find the first edge of a tetrahedron in the bisection order, that is its longest edge */
205: static PetscErrorCode SBRGetTetMaxEdge_Private(DMPlexTransform tr, PetscInt tet, PetscInt *maxedge)
206: {
207: DM dm;
208: PetscInt edges[6];
209: PetscBool prec;
211: PetscFunctionBeginHot;
212: PetscCall(DMPlexTransformGetDM(tr, &dm));
213: PetscCall(SBRGetTetEdges_Private(dm, tet, edges));
214: *maxedge = edges[0];
215: for (PetscInt k = 1; k < 6; ++k) {
216: PetscCall(SBREdgePrecedes_Private(tr, edges[k], *maxedge, &prec));
217: if (prec) *maxedge = edges[k];
218: }
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*
223: Mark the longest edge of a face-level cell (a triangle of a 3D mesh, or a maximal cell of a 2D
224: mesh), and the longest edge of every unprocessed cell above it. This is idempotent, so it is
225: safe to process a face again when it was marked both locally and on another process.
226: */
227: static PetscErrorCode SBRSplitFace_Private(DMPlexTransform tr, DMPlexPointQueue queue, PetscInt cell)
228: {
229: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
230: DM dm;
231: const PetscInt *cone, *tsupport;
232: PetscInt coneSize, tsupportSize, c, eval, maxedge;
233: PetscBool prec;
235: PetscFunctionBegin;
236: PetscCall(DMPlexTransformGetDM(tr, &dm));
237: PetscCall(DMPlexGetCone(dm, cell, &cone));
238: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
239: maxedge = cone[0];
240: for (c = 1; c < coneSize; ++c) {
241: PetscCall(SBREdgePrecedes_Private(tr, cone[c], maxedge, &prec));
242: if (prec) maxedge = cone[c];
243: }
244: PetscCall(DMLabelGetValue(sbr->splitPoints, maxedge, &eval));
245: if (eval != 1) {
246: PetscCall(DMLabelSetValue(sbr->splitPoints, maxedge, 1));
247: PetscCall(DMPlexPointQueueEnqueue(queue, maxedge));
248: }
249: PetscCall(DMLabelSetValue(sbr->splitPoints, cell, 2));
250: /* Propagate to the tetrahedra above this face, if any (empty in a 2D mesh) */
251: PetscCall(DMPlexGetSupport(dm, cell, &tsupport));
252: PetscCall(DMPlexGetSupportSize(dm, cell, &tsupportSize));
253: for (PetscInt ts = 0; ts < tsupportSize; ++ts) {
254: const PetscInt tet = tsupport[ts];
255: PetscInt tval, tmaxedge = -1, teval;
257: PetscCall(DMLabelGetValue(sbr->splitPoints, tet, &tval));
258: if (tval == 3) continue;
259: PetscCall(SBRGetTetMaxEdge_Private(tr, tet, &tmaxedge));
260: PetscCall(DMLabelGetValue(sbr->splitPoints, tmaxedge, &teval));
261: if (teval != 1) {
262: PetscCall(DMLabelSetValue(sbr->splitPoints, tmaxedge, 1));
263: PetscCall(DMPlexPointQueueEnqueue(queue, tmaxedge));
264: }
265: PetscCall(DMLabelSetValue(sbr->splitPoints, tet, 3));
266: }
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: /*
271: Mark local edges that should be split, ensuring conformity of the mesh skeleton.
273: This implements the closure step of Plaza & Carey, Section 3.1, generalized from 2D to 3D:
274: whenever an edge is marked, every non-conforming face containing it has its own longest edge
275: marked (as in the original 2D algorithm), and every tetrahedron touching such a face has its
276: own longest edge marked in turn (the "2.1/2.2" steps of the paper's 3D outline). In a 2D mesh,
277: the tetrahedron support of a face is empty and this reduces exactly to the original algorithm.
279: The queue also receives the points marked on other processes from the label propagation. A
280: face-level cell arriving this way is processed directly: the process that marked it has no
281: access to the cells above it on this process, so the local closure must be completed here.
282: Points of any other depth are ignored: vertices are never marked, and a cell arriving through
283: an overlapped point SF already had its closure ensured by its owner, whose edge marks arrive
284: through the same propagation.
285: */
286: static PetscErrorCode SBRSplitLocalEdges_Private(DMPlexTransform tr, DMPlexPointQueue queue)
287: {
288: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
289: DM dm;
290: PetscInt eStart, eEnd, fStart, fEnd;
292: PetscFunctionBegin;
293: PetscCall(DMPlexTransformGetDM(tr, &dm));
294: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
295: PetscCall(DMPlexGetDepthStratum(dm, 2, &fStart, &fEnd));
296: while (!DMPlexPointQueueEmpty(queue)) {
297: PetscInt p = -1;
299: PetscCall(DMPlexPointQueueDequeue(queue, &p));
300: if (p >= eStart && p < eEnd) {
301: const PetscInt *support;
302: PetscInt supportSize;
304: PetscCall(DMPlexGetSupport(dm, p, &support));
305: PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
306: for (PetscInt s = 0; s < supportSize; ++s) {
307: PetscInt cval;
309: PetscCall(DMLabelGetValue(sbr->splitPoints, support[s], &cval));
310: if (cval == 2) continue;
311: PetscCall(SBRSplitFace_Private(tr, queue, support[s]));
312: }
313: } else if (p >= fStart && p < fEnd) PetscCall(SBRSplitFace_Private(tr, queue, p));
314: }
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: static PetscErrorCode splitPoint(PETSC_UNUSED DMLabel label, PetscInt p, PETSC_UNUSED PetscInt val, PetscCtx ctx)
319: {
320: DMPlexPointQueue queue = (DMPlexPointQueue)ctx;
322: PetscFunctionBegin;
323: PetscCall(DMPlexPointQueueEnqueue(queue, p));
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /*
328: The 'splitPoints' label marks mesh points to be divided. It marks edges with 1, triangles with 2, and tetrahedra with 3.
329: Then the refinement type is calculated as follows:
331: RT_VERTEX: vertex
332: RT_EDGE: edge unsplit
333: RT_EDGE_SPLIT: edge split
334: RT_TRIANGLE: triangle unsplit
335: RT_TRIANGLE_SPLIT: maximal (2D) triangle with all edges split, subdivided regularly
336: RT_TRIANGLE_SPLIT_ij: triangle with edges i and j split, i bisected first
337: RT_TRIANGLE_SPLIT_k: triangle with only edge k split
338: RT_TRIANGLE_SPLIT_FAN_k: triangle with all edges split under a 3D cell, fanned around the midpoint of edge k
339: RT_TET: tetrahedron unsplit
340: RT_TET_SPLIT_BASE + code: tetrahedron with marked edges bisected in a given order
342: Following Plaza & Carey (2000), a marked cell is subdivided by bisecting its marked edges one at a
343: time, longest first: each bisection splits every subsimplex containing that edge, so the ordered
344: list of marked edges determines the whole subdivision, and the subdivision it induces on a face is
345: the bisection of the face's own marked edges in the same relative order. The classification below
346: therefore records the bisection order: for a pair of split triangle edges as the type
347: RT_TRIANGLE_SPLIT_ij, and for a tetrahedron as code = sum_i (e_i + 1) 7^i, the little-endian
348: base-7 encoding of its marked edges e_0 > e_1 > ... in bisection order, with edges e0..e5 numbered
349: as documented at SBRGetTetEdges_Private(). All length comparisons use SBREdgePrecedes_Private(),
350: whose geometric tie-breaking makes the order consistent between each face and the cells above it.
352: A fully marked triangle under a tetrahedron is fanned around the midpoint of its first-bisected
353: edge, matching the bisection cascade (the 4T partition of Rivara, Fig. 3 of the paper), and a
354: fully marked tetrahedron is bisected 6 times like any other marked tetrahedron. The regular
355: subdivisions are kept only for maximal (2D) triangles, where no compatibility with higher cells
356: constrains the interior of the subdivision, preserving the original 2D behavior of this transform.
357: */
358: typedef enum {
359: RT_VERTEX,
360: RT_EDGE,
361: RT_EDGE_SPLIT,
362: RT_TRIANGLE,
363: RT_TRIANGLE_SPLIT,
364: RT_TRIANGLE_SPLIT_01,
365: RT_TRIANGLE_SPLIT_10,
366: RT_TRIANGLE_SPLIT_12,
367: RT_TRIANGLE_SPLIT_21,
368: RT_TRIANGLE_SPLIT_20,
369: RT_TRIANGLE_SPLIT_02,
370: RT_TRIANGLE_SPLIT_0,
371: RT_TRIANGLE_SPLIT_1,
372: RT_TRIANGLE_SPLIT_2,
373: RT_TRIANGLE_SPLIT_FAN_0,
374: RT_TRIANGLE_SPLIT_FAN_1,
375: RT_TRIANGLE_SPLIT_FAN_2,
376: RT_TET,
377: RT_TET_SPLIT_BASE
378: } RefinementType;
380: static PetscErrorCode DMPlexTransformSetUp_SBR(DMPlexTransform tr)
381: {
382: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
383: DM dm;
384: DMLabel active;
385: PetscSF pointSF;
386: DMPlexPointQueue queue = NULL;
387: IS refineIS;
388: const PetscInt *refineCells;
389: PetscInt pStart, pEnd, p, eStart, eEnd, e, edgeLenSize, Nc, c;
390: PetscBool empty;
392: PetscFunctionBegin;
393: PetscCall(DMPlexTransformGetDM(tr, &dm));
394: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Split Points", &sbr->splitPoints));
395: /* Create edge lengths */
396: PetscCall(DMGetCoordinatesLocalSetUp(dm));
397: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
398: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen));
399: PetscCall(PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd));
400: for (e = eStart; e < eEnd; ++e) PetscCall(PetscSectionSetDof(sbr->secEdgeLen, e, 1));
401: PetscCall(PetscSectionSetUp(sbr->secEdgeLen));
402: PetscCall(PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize));
403: PetscCall(PetscCalloc1(edgeLenSize, &sbr->edgeLen));
404: /* Add edges of cells that are marked for refinement to edge queue */
405: PetscCall(DMPlexTransformGetActive(tr, &active));
406: PetscCheck(active, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONGSTATE, "DMPlexTransform must have an adaptation label in order to use SBR algorithm");
407: PetscCall(DMPlexPointQueueCreate(1024, &queue));
408: PetscCall(DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS));
409: PetscCall(DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc));
410: if (refineIS) PetscCall(ISGetIndices(refineIS, &refineCells));
411: for (c = 0; c < Nc; ++c) {
412: const PetscInt cell = refineCells[c];
413: PetscInt depth;
415: PetscCall(DMPlexGetPointDepth(dm, cell, &depth));
416: if (depth == 1) {
417: PetscCall(DMLabelSetValue(sbr->splitPoints, cell, 1));
418: PetscCall(DMPlexPointQueueEnqueue(queue, cell));
419: } else {
420: PetscInt *closure = NULL;
421: PetscInt Ncl;
423: PetscCall(DMLabelSetValue(sbr->splitPoints, cell, depth));
424: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
425: for (PetscInt cl = 0; cl < 2 * Ncl; cl += 2) {
426: const PetscInt edge = closure[cl];
428: if (edge >= eStart && edge < eEnd) {
429: PetscCall(DMLabelSetValue(sbr->splitPoints, edge, 1));
430: PetscCall(DMPlexPointQueueEnqueue(queue, edge));
431: }
432: }
433: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
434: }
435: }
436: if (refineIS) PetscCall(ISRestoreIndices(refineIS, &refineCells));
437: PetscCall(ISDestroy(&refineIS));
438: /* Setup communication */
439: PetscCall(DMGetPointSF(dm, &pointSF));
440: PetscCall(DMLabelPropagateBegin(sbr->splitPoints, pointSF));
441: /* While edge queue is not empty: */
442: PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty));
443: while (!empty) {
444: PetscCall(SBRSplitLocalEdges_Private(tr, queue));
445: /* Communicate marked edges
446: An easy implementation is to allocate an array the size of the number of points. We put the splitPoints marks into the
447: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
449: TODO: We could use in-place communication with a different SF
450: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
451: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
453: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
454: values will have 1+0=1 and old values will have 1+1=2. Loop over these, resetting the values to 1, and adding any new
455: edge to the queue.
456: */
457: PetscCall(DMLabelPropagatePush(sbr->splitPoints, pointSF, MPI_MAX, splitPoint, queue));
458: PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty));
459: }
460: PetscCall(DMLabelPropagateEnd(sbr->splitPoints, pointSF));
461: /* Calculate refineType for each cell */
462: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
463: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
464: for (p = pStart; p < pEnd; ++p) {
465: DMLabel trType = tr->trType;
466: DMPolytopeType ct;
467: PetscInt val;
469: PetscCall(DMPlexGetCellType(dm, p, &ct));
470: switch (ct) {
471: case DM_POLYTOPE_POINT:
472: PetscCall(DMLabelSetValue(trType, p, RT_VERTEX));
473: break;
474: case DM_POLYTOPE_SEGMENT:
475: PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
476: if (val == 1) PetscCall(DMLabelSetValue(trType, p, RT_EDGE_SPLIT));
477: else PetscCall(DMLabelSetValue(trType, p, RT_EDGE));
478: break;
479: case DM_POLYTOPE_TRIANGLE:
480: PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
481: if (val == 2) {
482: const PetscInt *cone;
483: PetscInt vals[3], i;
484: PetscBool prec;
486: PetscCall(DMPlexGetCone(dm, p, &cone));
487: for (i = 0; i < 3; ++i) {
488: PetscCall(DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i]));
489: vals[i] = vals[i] < 0 ? 0 : vals[i];
490: }
491: if (vals[0] && vals[1] && vals[2]) {
492: PetscInt suppSize, k = 0;
494: /* A triangle below a 3D cell must be subdivided compatibly with the bisection cascade of
495: the cells above it, fanning around the midpoint of its first-bisected edge; a maximal
496: (2D) triangle keeps the regular subdivision used since the original 2D implementation */
497: PetscCall(DMPlexGetSupportSize(dm, p, &suppSize));
498: if (suppSize) {
499: for (i = 1; i < 3; ++i) {
500: PetscCall(SBREdgePrecedes_Private(tr, cone[i], cone[k], &prec));
501: if (prec) k = i;
502: }
503: PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_FAN_0 + k));
504: } else PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT));
505: } else if (vals[0] && vals[1]) {
506: PetscCall(SBREdgePrecedes_Private(tr, cone[0], cone[1], &prec));
507: PetscCall(DMLabelSetValue(trType, p, prec ? RT_TRIANGLE_SPLIT_01 : RT_TRIANGLE_SPLIT_10));
508: } else if (vals[1] && vals[2]) {
509: PetscCall(SBREdgePrecedes_Private(tr, cone[1], cone[2], &prec));
510: PetscCall(DMLabelSetValue(trType, p, prec ? RT_TRIANGLE_SPLIT_12 : RT_TRIANGLE_SPLIT_21));
511: } else if (vals[2] && vals[0]) {
512: PetscCall(SBREdgePrecedes_Private(tr, cone[2], cone[0], &prec));
513: PetscCall(DMLabelSetValue(trType, p, prec ? RT_TRIANGLE_SPLIT_20 : RT_TRIANGLE_SPLIT_02));
514: } else if (vals[0]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_0));
515: else if (vals[1]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_1));
516: else if (vals[2]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_2));
517: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT " does not fit any refinement type (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")", p, vals[0], vals[1], vals[2]);
518: } else PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE));
519: break;
520: case DM_POLYTOPE_TETRAHEDRON:
521: PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
522: if (val == 3) {
523: PetscInt edges[6], ord[6], na = 0, code = 0;
525: PetscCall(SBRGetTetEdges_Private(dm, p, edges));
526: for (PetscInt k = 0; k < 6; ++k) {
527: PetscInt eval;
529: PetscCall(DMLabelGetValue(sbr->splitPoints, edges[k], &eval));
530: if (eval == 1) ord[na++] = k;
531: }
532: PetscCheck(na, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT " is marked for subdivision but has no marked edges", p);
533: /* Insertion sort of the marked edges into bisection order */
534: for (PetscInt i = 1; i < na; ++i) {
535: const PetscInt ei = ord[i];
536: PetscInt j;
538: for (j = i - 1; j >= 0; --j) {
539: PetscBool prec;
541: PetscCall(SBREdgePrecedes_Private(tr, edges[ei], edges[ord[j]], &prec));
542: if (!prec) break;
543: ord[j + 1] = ord[j];
544: }
545: ord[j + 1] = ei;
546: }
547: for (PetscInt i = na - 1; i >= 0; --i) code = code * 7 + (ord[i] + 1);
548: PetscCall(DMLabelSetValue(trType, p, RT_TET_SPLIT_BASE + code));
549: } else PetscCall(DMLabelSetValue(trType, p, RT_TET));
550: break;
551: default:
552: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
553: }
554: PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
555: }
556: /* Cleanup */
557: PetscCall(DMPlexPointQueueDestroy(&queue));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: /* Append one cone entry {ft, fn, acp[0..fn-1], r} plus its matching ornt entry o */
562: static PetscErrorCode SBRAppendCone_Private(PetscInt cone[], PetscInt *coff, PetscInt maxCone, PetscInt ornt[], PetscInt *ooff, PetscInt maxOrnt, DMPolytopeType ft, PetscInt fn, const PetscInt acp[], PetscInt r, PetscInt o)
563: {
564: PetscFunctionBeginHot;
565: PetscCheck(*coff + 3 + fn <= maxCone && *ooff + 1 <= maxOrnt, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cone buffer of size (%" PetscInt_FMT ", %" PetscInt_FMT ") is too small", maxCone, maxOrnt);
566: cone[(*coff)++] = (PetscInt)ft;
567: cone[(*coff)++] = fn;
568: for (PetscInt i = 0; i < fn; ++i) cone[(*coff)++] = acp[i];
569: cone[(*coff)++] = r;
570: ornt[(*ooff)++] = o;
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: /*
575: Symbolic subdivision machinery
577: Subdivisions are worked out on a canonical reference cell, in terms of symbols: for a triangle,
578: 0-2 denote its vertices and 3+e the midpoint of its edge e; for a tetrahedron, 0-3 denote its
579: vertices and 4+e the midpoint of its edge e. Everything is purely combinatorial, so the results
580: depend only on the refinement type, never on mesh coordinates.
581: */
583: typedef struct {
584: PetscInt nseg; /* Number of segments added inside the triangle */
585: PetscInt ntri; /* Number of subtriangles */
586: PetscInt seg[3][2]; /* The endpoints of each added segment, in its canonical direction */
587: PetscInt tri[4][3]; /* The vertices of each subtriangle, in canonical cone order */
588: } SBRTriSubdiv;
590: /* Map a triangle symbol authored against arrangement o of the triangle to the actual cone numbering */
591: static PetscInt SBRTriRelabel_Private(PetscInt o, PetscInt s)
592: {
593: const PetscInt *arr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, o);
595: if (s < 3) return triEdgeVert[arr[s * 2]][arr[s * 2 + 1] ? 1 : 0];
596: return 3 + arr[(s - 3) * 2];
597: }
599: /*
600: Compute the symbolic subdivision of a triangle: nsplit = 0 leaves it whole, nsplit = 1 bisects
601: edge 'first', nsplit = 2 bisects edge 'first' and then edge 'second', and nsplit = 3 bisects all
602: three edges starting with edge 'first' (the fan, where the order of the remaining two does not
603: affect the result). The subcells are listed exactly as produced, in replica order, by the identity
604: transform, SBRGetTriangleSplitSingle(), SBRGetTriangleSplitDouble(), and
605: SBRGetTriangleSplitFan_Private(): each case is authored here for the canonical arrangement and
606: relabeled through the same arrangement its template is dispatched with.
607: */
608: static PetscErrorCode SBRTriangleSubdiv_Private(PetscInt nsplit, PetscInt first, PetscInt second, SBRTriSubdiv *sub)
609: {
610: static const SBRTriSubdiv unsplitS = {0, 1, {{0, 0}}, {{0, 1, 2}}};
611: static const SBRTriSubdiv singleS = {
612: 1, 2, {{0, 4}},
613: {{0, 1, 4}, {4, 2, 0}}
614: };
615: static const SBRTriSubdiv doubleS = {
616: 2, 3, {{0, 4}, {4, 5}},
617: {{0, 1, 4}, {4, 2, 5}, {5, 0, 4}}
618: };
619: /* The reflected double split (second = first - 1), authored for first = 2, second = 1: the
620: template emits its subtriangle cones starting from different corners than a plain reflection
621: of the rotated instances, so it gets its own canonical tuples */
622: static const SBRTriSubdiv doubleRS = {
623: 2, 3, {{1, 5}, {5, 4}},
624: {{5, 0, 1}, {4, 2, 5}, {1, 4, 5}}
625: };
626: static const SBRTriSubdiv fanS = {
627: 3, 4, {{0, 4}, {4, 3}, {4, 5}},
628: {{0, 3, 4}, {3, 1, 4}, {5, 4, 2}, {0, 4, 5}}
629: };
630: const SBRTriSubdiv *canon;
631: PetscInt o;
633: PetscFunctionBeginHot;
634: switch (nsplit) {
635: case 0:
636: canon = &unsplitS;
637: o = 0;
638: break;
639: case 1:
640: case 3:
641: canon = nsplit == 1 ? &singleS : &fanS;
642: o = (first + 2) % 3;
643: break;
644: case 2:
645: if (second == (first + 1) % 3) {
646: canon = &doubleS;
647: o = (first + 2) % 3;
648: } else {
649: canon = &doubleRS;
650: o = (first + 1) % 3;
651: }
652: break;
653: default:
654: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of split edges %" PetscInt_FMT, nsplit);
655: }
656: sub->nseg = canon->nseg;
657: sub->ntri = canon->ntri;
658: for (PetscInt s = 0; s < canon->nseg; ++s)
659: for (PetscInt i = 0; i < 2; ++i) sub->seg[s][i] = SBRTriRelabel_Private(o, canon->seg[s][i]);
660: for (PetscInt t = 0; t < canon->ntri; ++t)
661: for (PetscInt i = 0; i < 3; ++i) sub->tri[t][i] = SBRTriRelabel_Private(o, canon->tri[t][i]);
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: #define SBR_TET_MAX_SEGS 16
666: #define SBR_TET_MAX_TRIS 32
667: #define SBR_TET_MAX_LEAVES 16
668: #define SBR_TET_MAX_CONE 2048
669: #define SBR_TET_MAX_ORNT 512
671: typedef struct {
672: PetscInt nseg, ntri, nleaf;
673: PetscInt seg[SBR_TET_MAX_SEGS][2]; /* Interior segments, in their canonical direction */
674: PetscInt tri[SBR_TET_MAX_TRIS][3]; /* Interior triangles, in canonical cone order */
675: PetscInt leaf[SBR_TET_MAX_LEAVES][4]; /* Subtetrahedra, in canonical vertex order */
676: } SBRTetSubdiv;
678: /* The faces of the tetrahedron containing a given symbol, as a bitmask over f0..f3 */
679: static PetscInt SBRTetSymbolFaces_Private(PetscInt s)
680: {
681: return s < 4 ? tetVertFaces[s] : tetEdgeFaces[s - 4];
682: }
684: static PetscInt SBRTupleCompare_Private(PetscInt n, const PetscInt a[], const PetscInt b[])
685: {
686: for (PetscInt i = 0; i < n; ++i)
687: if (a[i] != b[i]) return a[i] < b[i] ? -1 : 1;
688: return 0;
689: }
691: static PetscBool SBRTupleMatch_Private(PetscInt n, const PetscInt a[], const PetscInt b[])
692: {
693: for (PetscInt i = 0; i < n; ++i) {
694: PetscInt j;
696: for (j = 0; j < n; ++j)
697: if (a[i] == b[j]) break;
698: if (j == n) return PETSC_FALSE;
699: }
700: return PETSC_TRUE;
701: }
703: static void SBRTetDecodeOrder_Private(PetscInt code, PetscInt ord[6], PetscInt *na)
704: {
705: *na = 0;
706: while (code) {
707: ord[(*na)++] = code % 7 - 1;
708: code /= 7;
709: }
710: }
712: /* Add the interior segment (P, Q) if it is not already present */
713: static PetscErrorCode SBRTetAddSeg_Private(SBRTetSubdiv *sub, PetscInt P, PetscInt Q)
714: {
715: PetscFunctionBeginHot;
716: for (PetscInt s = 0; s < sub->nseg; ++s)
717: if ((sub->seg[s][0] == P && sub->seg[s][1] == Q) || (sub->seg[s][0] == Q && sub->seg[s][1] == P)) PetscFunctionReturn(PETSC_SUCCESS);
718: PetscCheck(sub->nseg < SBR_TET_MAX_SEGS, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Too many interior segments");
719: sub->seg[sub->nseg][0] = P;
720: sub->seg[sub->nseg][1] = Q;
721: ++sub->nseg;
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: /*
726: Subdivide the reference tetrahedron by bisecting the given edges in the given order. Each
727: bisection splits every current subsimplex containing the edge (Lemma 1.1 of Plaza & Carey (2000)):
728: the leaf subtetrahedra, the interior triangles created by earlier bisections, and implicitly the
729: faces of the tetrahedron, whose subdivision is not tracked here because each face is subdivided by
730: its own transform, consistently, since its refinement type orders its marked edges with the same
731: comparison (see RefinementType). Only the segments and triangles strictly interior to the
732: tetrahedron are recorded as its own subcells; everything on the boundary belongs to a face.
733: */
734: static PetscErrorCode SBRTetSubdivide_Private(const PetscInt ord[], PetscInt na, SBRTetSubdiv *sub)
735: {
736: PetscFunctionBegin;
737: sub->nseg = 0;
738: sub->ntri = 0;
739: sub->nleaf = 1;
740: for (PetscInt i = 0; i < 4; ++i) sub->leaf[0][i] = i;
741: for (PetscInt s = 0; s < na; ++s) {
742: const PetscInt e = ord[s], a = tetEdgeVert[e][0], b = tetEdgeVert[e][1], m = 4 + e;
744: for (PetscInt t = 0; t < sub->ntri; ++t) {
745: PetscInt pa = -1, pb = -1;
747: for (PetscInt i = 0; i < 3; ++i) {
748: if (sub->tri[t][i] == a) pa = i;
749: else if (sub->tri[t][i] == b) pb = i;
750: }
751: if (pa < 0 || pb < 0) continue;
752: PetscCheck(sub->ntri < SBR_TET_MAX_TRIS, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Too many interior triangles");
753: PetscCall(PetscArraymove(&sub->tri[t + 2][0], &sub->tri[t + 1][0], 3 * (sub->ntri - t - 1)));
754: PetscCall(PetscArraycpy(sub->tri[t + 1], sub->tri[t], 3));
755: PetscCall(SBRTetAddSeg_Private(sub, sub->tri[t][3 - pa - pb], m));
756: sub->tri[t][pb] = m;
757: sub->tri[t + 1][pa] = m;
758: ++sub->ntri;
759: ++t;
760: }
761: for (PetscInt l = 0; l < sub->nleaf; ++l) {
762: PetscInt pa = -1, pb = -1, *newtri;
764: for (PetscInt i = 0; i < 4; ++i) {
765: if (sub->leaf[l][i] == a) pa = i;
766: else if (sub->leaf[l][i] == b) pb = i;
767: }
768: if (pa < 0 || pb < 0) continue;
769: PetscCheck(sub->nleaf < SBR_TET_MAX_LEAVES && sub->ntri < SBR_TET_MAX_TRIS, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Too many subcells");
770: PetscCall(PetscArraymove(&sub->leaf[l + 2][0], &sub->leaf[l + 1][0], 4 * (sub->nleaf - l - 1)));
771: PetscCall(PetscArraycpy(sub->leaf[l + 1], sub->leaf[l], 4));
772: sub->leaf[l][pb] = m;
773: sub->leaf[l + 1][pa] = m;
774: ++sub->nleaf;
775: /* Record the interior triangle separating the two children, and its interior median edges */
776: newtri = sub->tri[sub->ntri++];
777: for (PetscInt i = 0; i < 3; ++i) newtri[i] = sub->leaf[l + 1][tetOppFace[pb][i]];
778: for (PetscInt i = 0; i < 3; ++i) {
779: const PetscInt P = newtri[i], Q = newtri[(i + 1) % 3], other = P == m ? Q : P;
781: if (P != m && Q != m) continue;
782: if (SBRTetSymbolFaces_Private(m) & SBRTetSymbolFaces_Private(other)) continue;
783: PetscCall(SBRTetAddSeg_Private(sub, other, m));
784: }
785: ++l;
786: }
787: }
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: /* The subdivision of face f of the tetrahedron, relabeled from face symbols into tetrahedron symbols */
792: static PetscErrorCode SBRTetFaceSubdiv_Private(const PetscInt ord[], PetscInt na, PetscInt f, SBRTriSubdiv *sub)
793: {
794: PetscInt fl[3], nfm = 0;
796: PetscFunctionBeginHot;
797: for (PetscInt s = 0; s < na; ++s)
798: for (PetscInt l = 0; l < 3; ++l)
799: if (tetFaceEdge[f][l] == ord[s]) fl[nfm++] = l;
800: PetscCall(SBRTriangleSubdiv_Private(nfm, nfm > 0 ? fl[0] : -1, nfm > 1 ? fl[1] : -1, sub));
801: for (PetscInt s = 0; s < sub->nseg; ++s)
802: for (PetscInt i = 0; i < 2; ++i) sub->seg[s][i] = sub->seg[s][i] < 3 ? tetFaceVert[f][sub->seg[s][i]] : 4 + tetFaceEdge[f][sub->seg[s][i] - 3];
803: for (PetscInt t = 0; t < sub->ntri; ++t)
804: for (PetscInt i = 0; i < 3; ++i) sub->tri[t][i] = sub->tri[t][i] < 3 ? tetFaceVert[f][sub->tri[t][i]] : 4 + tetFaceEdge[f][sub->tri[t][i] - 3];
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: /* Emit the cone entry for the segment (P, Q): a whole original edge, half of a split edge, a
809: segment added inside a face, or a segment interior to the tetrahedron */
810: static PetscErrorCode SBRTetEmitSegRef_Private(const SBRTetSubdiv *sub, const SBRTriSubdiv fsub[], PetscInt P, PetscInt Q, PetscInt cone[], PetscInt *coff, PetscInt ornt[], PetscInt *ooff)
811: {
812: PetscFunctionBeginHot;
813: if (P < 4 && Q < 4) {
814: const PetscInt e = SBREdgeIndexFromVerts_Private(P, Q);
816: PetscCall(SBRAppendCone_Private(cone, coff, SBR_TET_MAX_CONE, ornt, ooff, SBR_TET_MAX_ORNT, DM_POLYTOPE_SEGMENT, 2, tetEdgeFaceLoc[e], 0, tetEdgeVert[e][0] == P ? 0 : -1));
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
819: if (P < 4 || Q < 4) {
820: const PetscInt v = P < 4 ? P : Q, e = (P < 4 ? Q : P) - 4;
822: if (tetEdgeVert[e][0] == v || tetEdgeVert[e][1] == v) {
823: const PetscInt r = tetEdgeVert[e][0] == v ? 0 : 1;
825: /* Half of split edge e: replica 0 runs from the edge tail to the midpoint, replica 1 from the midpoint to the head */
826: PetscCall(SBRAppendCone_Private(cone, coff, SBR_TET_MAX_CONE, ornt, ooff, SBR_TET_MAX_ORNT, DM_POLYTOPE_SEGMENT, 2, tetEdgeFaceLoc[e], r, (r == 0) == (P == v) ? 0 : -1));
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
829: }
830: {
831: const PetscInt mask = SBRTetSymbolFaces_Private(P) & SBRTetSymbolFaces_Private(Q);
833: if (mask) {
834: PetscInt f = 0;
836: while (!(mask & (1 << f))) ++f;
837: for (PetscInt r = 0; r < fsub[f].nseg; ++r) {
838: if ((fsub[f].seg[r][0] == P && fsub[f].seg[r][1] == Q) || (fsub[f].seg[r][0] == Q && fsub[f].seg[r][1] == P)) {
839: const PetscInt acp[1] = {f};
841: PetscCall(SBRAppendCone_Private(cone, coff, SBR_TET_MAX_CONE, ornt, ooff, SBR_TET_MAX_ORNT, DM_POLYTOPE_SEGMENT, 1, acp, r, fsub[f].seg[r][0] == P ? 0 : -1));
842: PetscFunctionReturn(PETSC_SUCCESS);
843: }
844: }
845: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Segment (%" PetscInt_FMT ", %" PetscInt_FMT ") not found in the subdivision of face %" PetscInt_FMT, P, Q, f);
846: }
847: for (PetscInt r = 0; r < sub->nseg; ++r) {
848: if ((sub->seg[r][0] == P && sub->seg[r][1] == Q) || (sub->seg[r][0] == Q && sub->seg[r][1] == P)) {
849: PetscCall(SBRAppendCone_Private(cone, coff, SBR_TET_MAX_CONE, ornt, ooff, SBR_TET_MAX_ORNT, DM_POLYTOPE_SEGMENT, 0, NULL, r, sub->seg[r][0] == P ? 0 : -1));
850: PetscFunctionReturn(PETSC_SUCCESS);
851: }
852: }
853: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior segment (%" PetscInt_FMT ", %" PetscInt_FMT ") not found", P, Q);
854: }
855: }
857: /* Emit the cone entry for the triangle with vertex cycle C: a subtriangle of a face, or a triangle
858: interior to the tetrahedron */
859: static PetscErrorCode SBRTetEmitTriRef_Private(const SBRTetSubdiv *sub, const SBRTriSubdiv fsub[], const PetscInt C[3], PetscInt cone[], PetscInt *coff, PetscInt ornt[], PetscInt *ooff)
860: {
861: const PetscInt mask = SBRTetSymbolFaces_Private(C[0]) & SBRTetSymbolFaces_Private(C[1]) & SBRTetSymbolFaces_Private(C[2]);
862: PetscInt o;
864: PetscFunctionBeginHot;
865: if (mask) {
866: PetscInt f = 0;
868: while (!(mask & (1 << f))) ++f;
869: for (PetscInt r = 0; r < fsub[f].ntri; ++r) {
870: if (SBRTupleMatch_Private(3, fsub[f].tri[r], C)) {
871: const PetscInt acp[1] = {f};
873: PetscCall(DMPolytopeGetVertexOrientation(DM_POLYTOPE_TRIANGLE, fsub[f].tri[r], C, &o));
874: PetscCall(SBRAppendCone_Private(cone, coff, SBR_TET_MAX_CONE, ornt, ooff, SBR_TET_MAX_ORNT, DM_POLYTOPE_TRIANGLE, 1, acp, r, o));
875: PetscFunctionReturn(PETSC_SUCCESS);
876: }
877: }
878: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Triangle (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") not found in the subdivision of face %" PetscInt_FMT, C[0], C[1], C[2], f);
879: }
880: for (PetscInt r = 0; r < sub->ntri; ++r) {
881: if (SBRTupleMatch_Private(3, sub->tri[r], C)) {
882: PetscCall(DMPolytopeGetVertexOrientation(DM_POLYTOPE_TRIANGLE, sub->tri[r], C, &o));
883: PetscCall(SBRAppendCone_Private(cone, coff, SBR_TET_MAX_CONE, ornt, ooff, SBR_TET_MAX_ORNT, DM_POLYTOPE_TRIANGLE, 0, NULL, r, o));
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
886: }
887: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior triangle (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") not found", C[0], C[1], C[2]);
888: }
890: /* Emit the cone and orientation data for the subdivision of a tetrahedron */
891: static PetscErrorCode SBRTetEmit_Private(const PetscInt ord[], PetscInt na, const SBRTetSubdiv *sub, PetscInt *Nt, DMPolytopeType target[], PetscInt size[], PetscInt cone[], PetscInt *Ncone, PetscInt ornt[], PetscInt *Nornt)
892: {
893: SBRTriSubdiv fsub[4];
894: PetscInt coff = 0, ooff = 0, n = 0;
896: PetscFunctionBegin;
897: for (PetscInt f = 0; f < 4; ++f) PetscCall(SBRTetFaceSubdiv_Private(ord, na, f, &fsub[f]));
898: for (PetscInt s = 0; s < sub->nseg; ++s) {
899: for (PetscInt i = 0; i < 2; ++i) {
900: const PetscInt sym = sub->seg[s][i];
902: if (sym < 4) PetscCall(SBRAppendCone_Private(cone, &coff, SBR_TET_MAX_CONE, ornt, &ooff, SBR_TET_MAX_ORNT, DM_POLYTOPE_POINT, 3, tetVertPath[sym], 0, 0));
903: else PetscCall(SBRAppendCone_Private(cone, &coff, SBR_TET_MAX_CONE, ornt, &ooff, SBR_TET_MAX_ORNT, DM_POLYTOPE_POINT, 2, tetEdgeFaceLoc[sym - 4], 0, 0));
904: }
905: }
906: for (PetscInt t = 0; t < sub->ntri; ++t)
907: for (PetscInt i = 0; i < 3; ++i) PetscCall(SBRTetEmitSegRef_Private(sub, fsub, sub->tri[t][i], sub->tri[t][(i + 1) % 3], cone, &coff, ornt, &ooff));
908: for (PetscInt l = 0; l < sub->nleaf; ++l) {
909: for (PetscInt f = 0; f < 4; ++f) {
910: /* The face f of the child tetrahedron is the face opposite its vertex 3 - f */
911: const PetscInt mp = 3 - f;
912: PetscInt C[3];
914: for (PetscInt i = 0; i < 3; ++i) C[i] = sub->leaf[l][tetOppFace[mp][i]];
915: PetscCall(SBRTetEmitTriRef_Private(sub, fsub, C, cone, &coff, ornt, &ooff));
916: }
917: }
918: if (sub->nseg) {
919: target[n] = DM_POLYTOPE_SEGMENT;
920: size[n] = sub->nseg;
921: ++n;
922: }
923: target[n] = DM_POLYTOPE_TRIANGLE;
924: size[n] = sub->ntri;
925: ++n;
926: target[n] = DM_POLYTOPE_TETRAHEDRON;
927: size[n] = sub->nleaf;
928: ++n;
929: *Nt = n;
930: *Ncone = coff;
931: *Nornt = ooff;
932: PetscFunctionReturn(PETSC_SUCCESS);
933: }
935: /* Get the subdivision data for the encoded bisection order, generating and caching it on first use */
936: static PetscErrorCode SBRGetTetSplit_Private(DMPlexTransform tr, PetscInt code, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
937: {
938: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
939: PetscInt i;
941: PetscFunctionBeginHot;
942: for (i = 0; i < sbr->Ncache; ++i)
943: if (sbr->cacheCode[i] == code) break;
944: if (i == sbr->Ncache) {
945: SBRTetSubdiv sub;
946: DMPolytopeType targetTmp[3];
947: PetscInt sizeTmp[3], coneTmp[SBR_TET_MAX_CONE], orntTmp[SBR_TET_MAX_ORNT];
948: PetscInt ord[6], na, NtTmp, Ncone, Nornt;
950: SBRTetDecodeOrder_Private(code, ord, &na);
951: PetscCall(SBRTetSubdivide_Private(ord, na, &sub));
952: PetscCall(SBRTetEmit_Private(ord, na, &sub, &NtTmp, targetTmp, sizeTmp, coneTmp, &Ncone, orntTmp, &Nornt));
953: if (sbr->Ncache == sbr->maxCache) {
954: const PetscInt newmax = PetscMax(2 * sbr->maxCache, 16);
955: PetscInt *newCode, *newNt, **newSize, **newCone, **newOrnt;
956: DMPolytopeType **newTarget;
958: PetscCall(PetscMalloc6(newmax, &newCode, newmax, &newNt, newmax, &newTarget, newmax, &newSize, newmax, &newCone, newmax, &newOrnt));
959: PetscCall(PetscArraycpy(newCode, sbr->cacheCode, sbr->Ncache));
960: PetscCall(PetscArraycpy(newNt, sbr->cacheNt, sbr->Ncache));
961: PetscCall(PetscArraycpy(newTarget, sbr->cacheTarget, sbr->Ncache));
962: PetscCall(PetscArraycpy(newSize, sbr->cacheSize, sbr->Ncache));
963: PetscCall(PetscArraycpy(newCone, sbr->cacheCone, sbr->Ncache));
964: PetscCall(PetscArraycpy(newOrnt, sbr->cacheOrnt, sbr->Ncache));
965: PetscCall(PetscFree6(sbr->cacheCode, sbr->cacheNt, sbr->cacheTarget, sbr->cacheSize, sbr->cacheCone, sbr->cacheOrnt));
966: sbr->cacheCode = newCode;
967: sbr->cacheNt = newNt;
968: sbr->cacheTarget = newTarget;
969: sbr->cacheSize = newSize;
970: sbr->cacheCone = newCone;
971: sbr->cacheOrnt = newOrnt;
972: sbr->maxCache = newmax;
973: }
974: PetscCall(PetscMalloc1(NtTmp, &sbr->cacheTarget[i]));
975: PetscCall(PetscMalloc1(NtTmp, &sbr->cacheSize[i]));
976: PetscCall(PetscMalloc1(Ncone, &sbr->cacheCone[i]));
977: PetscCall(PetscMalloc1(Nornt, &sbr->cacheOrnt[i]));
978: PetscCall(PetscArraycpy(sbr->cacheTarget[i], targetTmp, NtTmp));
979: PetscCall(PetscArraycpy(sbr->cacheSize[i], sizeTmp, NtTmp));
980: PetscCall(PetscArraycpy(sbr->cacheCone[i], coneTmp, Ncone));
981: PetscCall(PetscArraycpy(sbr->cacheOrnt[i], orntTmp, Nornt));
982: sbr->cacheCode[i] = code;
983: sbr->cacheNt[i] = NtTmp;
984: ++sbr->Ncache;
985: }
986: *Nt = sbr->cacheNt[i];
987: *target = sbr->cacheTarget[i];
988: *size = sbr->cacheSize[i];
989: *cone = sbr->cacheCone[i];
990: *ornt = sbr->cacheOrnt[i];
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: /*
995: Map the subcell (r, o) of a split triangle, authored against arrangement so of the triangle, to
996: the matching subcell of the triangle's own production: both subdivisions are reconstructed
997: symbolically, the authored subcell is relabeled through the arrangement, and matched against the
998: actual subcells
999: */
1000: static PetscErrorCode SBRTriangleOrient_Private(PetscInt rt, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1001: {
1002: const PetscInt *arr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, so);
1003: SBRTriSubdiv subT, subV;
1004: const PetscInt *tupleV, *canonT;
1005: PetscInt inv[3], tuple[3], n, nlist, nsplit, first, second = -1, r2, dO;
1007: PetscFunctionBeginHot;
1008: switch (rt) {
1009: case RT_TRIANGLE_SPLIT_0:
1010: case RT_TRIANGLE_SPLIT_1:
1011: case RT_TRIANGLE_SPLIT_2:
1012: nsplit = 1;
1013: first = rt - RT_TRIANGLE_SPLIT_0;
1014: break;
1015: case RT_TRIANGLE_SPLIT_FAN_0:
1016: case RT_TRIANGLE_SPLIT_FAN_1:
1017: case RT_TRIANGLE_SPLIT_FAN_2:
1018: nsplit = 3;
1019: first = rt - RT_TRIANGLE_SPLIT_FAN_0;
1020: break;
1021: case RT_TRIANGLE_SPLIT_01:
1022: case RT_TRIANGLE_SPLIT_10:
1023: case RT_TRIANGLE_SPLIT_12:
1024: case RT_TRIANGLE_SPLIT_21:
1025: case RT_TRIANGLE_SPLIT_20:
1026: case RT_TRIANGLE_SPLIT_02: {
1027: static const PetscInt pairs[6][2] = {
1028: {0, 1},
1029: {1, 0},
1030: {1, 2},
1031: {2, 1},
1032: {2, 0},
1033: {0, 2}
1034: };
1036: nsplit = 2;
1037: first = pairs[rt - RT_TRIANGLE_SPLIT_01][0];
1038: second = pairs[rt - RT_TRIANGLE_SPLIT_01][1];
1039: } break;
1040: default:
1041: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid refinement type %" PetscInt_FMT, rt);
1042: }
1043: for (PetscInt l = 0; l < 3; ++l) inv[arr[l * 2]] = l;
1044: PetscCall(SBRTriangleSubdiv_Private(nsplit, first, second, &subT));
1045: PetscCall(SBRTriangleSubdiv_Private(nsplit, inv[first], second < 0 ? -1 : inv[second], &subV));
1046: if (tct == DM_POLYTOPE_SEGMENT) {
1047: n = 2;
1048: nlist = subT.nseg;
1049: tupleV = subV.seg[r];
1050: } else {
1051: n = 3;
1052: nlist = subT.ntri;
1053: tupleV = subV.tri[r];
1054: }
1055: for (PetscInt i = 0; i < n; ++i) tuple[i] = SBRTriRelabel_Private(so, tupleV[i]);
1056: for (r2 = 0; r2 < nlist; ++r2) {
1057: canonT = tct == DM_POLYTOPE_SEGMENT ? subT.seg[r2] : subT.tri[r2];
1058: if (SBRTupleMatch_Private(n, canonT, tuple)) break;
1059: }
1060: PetscCheck(r2 < nlist, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No matching %s subcell for replica %" PetscInt_FMT " at orientation %" PetscInt_FMT, DMPolytopeTypes[tct], r, so);
1061: PetscCall(DMPolytopeGetVertexOrientation(tct, canonT, tuple, &dO));
1062: *rnew = r2;
1063: *onew = DMPolytopeTypeComposeOrientation(tct, o, dO);
1064: PetscFunctionReturn(PETSC_SUCCESS);
1065: }
1067: /* The images of the tetrahedron edges under the vertex permutation varr */
1068: static void SBRTetEdgeImage_Private(const PetscInt varr[], PetscInt medge[6])
1069: {
1070: for (PetscInt e = 0; e < 6; ++e) medge[e] = SBREdgeIndexFromVerts_Private(varr[tetEdgeVert[e][0]], varr[tetEdgeVert[e][1]]);
1071: }
1073: /* The tetrahedron analogue of SBRTriangleOrient_Private() */
1074: static PetscErrorCode SBRTetOrient_Private(PetscInt code, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1075: {
1076: const PetscInt *varr = DMPolytopeTypeGetVertexArrangement(DM_POLYTOPE_TETRAHEDRON, so);
1077: SBRTetSubdiv subT, subV;
1078: const PetscInt *tupleV, *canonT = NULL;
1079: PetscInt medge[6], minv[6], ordT[6], ordV[6], na, tuple[4], n, nlist, r2, dO;
1081: PetscFunctionBeginHot;
1082: SBRTetDecodeOrder_Private(code, ordT, &na);
1083: SBRTetEdgeImage_Private(varr, medge);
1084: for (PetscInt e = 0; e < 6; ++e) minv[medge[e]] = e;
1085: for (PetscInt i = 0; i < na; ++i) ordV[i] = minv[ordT[i]];
1086: PetscCall(SBRTetSubdivide_Private(ordT, na, &subT));
1087: PetscCall(SBRTetSubdivide_Private(ordV, na, &subV));
1088: switch (tct) {
1089: case DM_POLYTOPE_SEGMENT:
1090: n = 2;
1091: nlist = subT.nseg;
1092: tupleV = subV.seg[r];
1093: break;
1094: case DM_POLYTOPE_TRIANGLE:
1095: n = 3;
1096: nlist = subT.ntri;
1097: tupleV = subV.tri[r];
1098: break;
1099: case DM_POLYTOPE_TETRAHEDRON:
1100: n = 4;
1101: nlist = subT.nleaf;
1102: tupleV = subV.leaf[r];
1103: break;
1104: default:
1105: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid target type %s", DMPolytopeTypes[tct]);
1106: }
1107: for (PetscInt i = 0; i < n; ++i) tuple[i] = tupleV[i] < 4 ? varr[tupleV[i]] : 4 + medge[tupleV[i] - 4];
1108: for (r2 = 0; r2 < nlist; ++r2) {
1109: canonT = tct == DM_POLYTOPE_SEGMENT ? subT.seg[r2] : (tct == DM_POLYTOPE_TRIANGLE ? subT.tri[r2] : subT.leaf[r2]);
1110: if (SBRTupleMatch_Private(n, canonT, tuple)) break;
1111: }
1112: PetscCheck(r2 < nlist, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No matching %s subcell for replica %" PetscInt_FMT " at orientation %" PetscInt_FMT, DMPolytopeTypes[tct], r, so);
1113: PetscCall(DMPolytopeGetVertexOrientation(tct, canonT, tuple, &dO));
1114: *rnew = r2;
1115: *onew = DMPolytopeTypeComposeOrientation(tct, o, dO);
1116: PetscFunctionReturn(PETSC_SUCCESS);
1117: }
1119: /*
1120: Validate the subdivision generator against the enumeration of Plaza & Carey (2000): over all
1121: possible edge length orders and all conforming sets of marked edges, there are exactly 51
1122: subdivision classes up to rotation, distributed over the number of marked edges na = 1..6 as
1123: 1, 3, 9, 17, 15, 6 (Table 4 of the paper). This runs the generator on every configuration, which
1124: also exercises all internal consistency checks in the emission of the subdivision data.
1125: */
1126: static PetscErrorCode DMPlexTransformSBRValidate_Private(DMPlexTransform tr)
1127: {
1128: const PetscInt expected[6] = {1, 3, 9, 17, 15, 6};
1129: const PetscInt fact[6] = {120, 24, 6, 2, 1, 1};
1130: const PetscInt keyLen = 2 + 4 * SBR_TET_MAX_LEAVES;
1131: PetscBT btCode;
1132: PetscInt *codes, *keys, *nas, counts[6] = {0, 0, 0, 0, 0, 0};
1133: PetscInt Ncodes = 0, maxCodes = 4096;
1135: PetscFunctionBegin;
1136: PetscCall(PetscBTCreate(117649, &btCode)); /* 7^6 possible codes */
1137: PetscCall(PetscMalloc3(maxCodes, &codes, maxCodes * keyLen, &keys, maxCodes, &nas));
1138: for (PetscInt idx = 0; idx < 720; ++idx) {
1139: PetscInt avail[6] = {0, 1, 2, 3, 4, 5}, perm[6], prio[6], k = idx;
1141: for (PetscInt i = 0; i < 6; ++i) {
1142: const PetscInt pos = k / fact[i];
1144: k %= fact[i];
1145: perm[i] = avail[pos];
1146: for (PetscInt j = pos; j < 5 - i; ++j) avail[j] = avail[j + 1];
1147: }
1148: for (PetscInt i = 0; i < 6; ++i) prio[perm[i]] = i;
1149: for (PetscInt mask = 1; mask < 64; ++mask) {
1150: PetscInt ord[6], na = 0, code = 0;
1151: PetscBool stable = (mask & (1 << perm[0])) ? PETSC_TRUE : PETSC_FALSE;
1153: /* The set is conforming when the longest edge overall and of each affected face is marked */
1154: for (PetscInt f = 0; f < 4 && stable; ++f) {
1155: PetscInt fmax = tetFaceEdge[f][0], nm = 0;
1157: for (PetscInt l = 0; l < 3; ++l) {
1158: if (mask & (1 << tetFaceEdge[f][l])) ++nm;
1159: if (prio[tetFaceEdge[f][l]] < prio[fmax]) fmax = tetFaceEdge[f][l];
1160: }
1161: if (nm && !(mask & (1 << fmax))) stable = PETSC_FALSE;
1162: }
1163: if (!stable) continue;
1164: for (PetscInt i = 0; i < 6; ++i)
1165: if (mask & (1 << perm[i])) ord[na++] = perm[i];
1166: for (PetscInt i = na - 1; i >= 0; --i) code = code * 7 + (ord[i] + 1);
1167: if (!PetscBTLookupSet(btCode, code)) {
1168: PetscCheck(Ncodes < maxCodes, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Too many configurations");
1169: codes[Ncodes] = code;
1170: nas[Ncodes] = na;
1171: ++Ncodes;
1172: }
1173: }
1174: }
1175: for (PetscInt c = 0; c < Ncodes; ++c) {
1176: SBRTetSubdiv sub;
1177: DMPolytopeType targetTmp[3];
1178: PetscInt sizeTmp[3], coneTmp[SBR_TET_MAX_CONE], orntTmp[SBR_TET_MAX_ORNT];
1179: PetscInt ord[6], na, NtTmp, Ncone, Nornt;
1180: PetscInt *key = &keys[c * keyLen];
1182: SBRTetDecodeOrder_Private(codes[c], ord, &na);
1183: PetscCall(SBRTetSubdivide_Private(ord, na, &sub));
1184: PetscCall(SBRTetEmit_Private(ord, na, &sub, &NtTmp, targetTmp, sizeTmp, coneTmp, &Ncone, orntTmp, &Nornt));
1185: /* The class key is the smallest relabeling, over all 12 rotations, of the leaf set together
1186: with the first bisected edge: the paper counts length configurations, which distinguish two
1187: marked sets producing the same subdivision when their longest edges do not correspond */
1188: for (PetscInt so = 0; so < 12; ++so) {
1189: const PetscInt *varr = DMPolytopeTypeGetVertexArrangement(DM_POLYTOPE_TETRAHEDRON, so);
1190: PetscInt medge[6], cand[2 + 4 * SBR_TET_MAX_LEAVES];
1192: SBRTetEdgeImage_Private(varr, medge);
1193: for (PetscInt i = 0; i < keyLen; ++i) cand[i] = -1;
1194: cand[0] = sub.nleaf;
1195: cand[1] = medge[ord[0]];
1196: for (PetscInt l = 0; l < sub.nleaf; ++l) {
1197: PetscInt *t = &cand[2 + 4 * l];
1199: for (PetscInt i = 0; i < 4; ++i) t[i] = sub.leaf[l][i] < 4 ? varr[sub.leaf[l][i]] : 4 + medge[sub.leaf[l][i] - 4];
1200: PetscCall(PetscSortInt(4, t));
1201: }
1202: /* Insertion sort of the leaf records, then keep the smallest candidate */
1203: for (PetscInt l = 1; l < sub.nleaf; ++l) {
1204: PetscInt t[4], j;
1206: PetscCall(PetscArraycpy(t, &cand[2 + 4 * l], 4));
1207: for (j = l - 1; j >= 0 && SBRTupleCompare_Private(4, &cand[2 + 4 * j], t) > 0; --j) PetscCall(PetscArraycpy(&cand[2 + 4 * (j + 1)], &cand[2 + 4 * j], 4));
1208: PetscCall(PetscArraycpy(&cand[2 + 4 * (j + 1)], t, 4));
1209: }
1210: if (so == 0 || SBRTupleCompare_Private(keyLen, cand, key) < 0) PetscCall(PetscArraycpy(key, cand, keyLen));
1211: }
1212: }
1213: for (PetscInt c = 0; c < Ncodes; ++c) {
1214: PetscInt d;
1216: for (d = 0; d < c; ++d) {
1217: PetscBool same;
1219: if (nas[d] != nas[c]) continue;
1220: PetscCall(PetscMemcmp(&keys[d * keyLen], &keys[c * keyLen], keyLen * sizeof(PetscInt), &same));
1221: if (same) break;
1222: }
1223: if (d == c) ++counts[nas[c] - 1];
1224: }
1225: for (PetscInt na = 1; na <= 6; ++na)
1226: PetscCheck(counts[na - 1] == expected[na - 1], PetscObjectComm((PetscObject)tr), PETSC_ERR_PLIB, "Found %" PetscInt_FMT " subdivision classes with %" PetscInt_FMT " marked edges, expected %" PetscInt_FMT, counts[na - 1], na, expected[na - 1]);
1227: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tr), "SBR: validated %" PetscInt_FMT " conforming configurations covering 51 subdivision classes (1 3 9 17 15 6 for 1-6 bisected edges)\n", Ncodes));
1228: PetscCall(PetscFree3(codes, keys, nas));
1229: PetscCall(PetscBTDestroy(&btCode));
1230: PetscFunctionReturn(PETSC_SUCCESS);
1231: }
1233: static PetscErrorCode DMPlexTransformGetSubcellOrientation_SBR(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1234: {
1235: PetscInt rt;
1237: PetscFunctionBeginHot;
1238: PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
1239: *rnew = r;
1240: *onew = o;
1241: switch (rt) {
1242: case RT_TRIANGLE_SPLIT_01:
1243: case RT_TRIANGLE_SPLIT_10:
1244: case RT_TRIANGLE_SPLIT_12:
1245: case RT_TRIANGLE_SPLIT_21:
1246: case RT_TRIANGLE_SPLIT_20:
1247: case RT_TRIANGLE_SPLIT_02:
1248: case RT_TRIANGLE_SPLIT_0:
1249: case RT_TRIANGLE_SPLIT_1:
1250: case RT_TRIANGLE_SPLIT_2:
1251: case RT_TRIANGLE_SPLIT_FAN_0:
1252: case RT_TRIANGLE_SPLIT_FAN_1:
1253: case RT_TRIANGLE_SPLIT_FAN_2:
1254: if (so) PetscCall(SBRTriangleOrient_Private(rt, so, tct, r, o, rnew, onew));
1255: break;
1256: case RT_EDGE_SPLIT:
1257: case RT_TRIANGLE_SPLIT:
1258: PetscCall(DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew));
1259: break;
1260: default:
1261: if (rt >= RT_TET_SPLIT_BASE) {
1262: if (so) PetscCall(SBRTetOrient_Private(rt - RT_TET_SPLIT_BASE, so, tct, r, o, rnew, onew));
1263: } else PetscCall(DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew));
1264: }
1265: PetscFunctionReturn(PETSC_SUCCESS);
1266: }
1268: /* Add 1 edge inside this triangle, making 2 new triangles.
1269: 2
1270: |\
1271: | \
1272: | \
1273: | \
1274: | 1
1275: | \
1276: | B \
1277: 2 1
1278: | / \
1279: | ____/ 0
1280: |/ A \
1281: 0-----0-----1
1282: */
1283: static PetscErrorCode SBRGetTriangleSplitSingle(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1284: {
1285: const PetscInt *arr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, o);
1286: static DMPolytopeType triT1[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
1287: static PetscInt triS1[] = {1, 2};
1288: static PetscInt triC1[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1289: DM_POLYTOPE_SEGMENT, 0, 0};
1290: static PetscInt triO1[] = {0, 0, 0, 0, -1, 0, 0, 0};
1292: PetscFunctionBeginHot;
1293: /* To get the other divisions, we reorient the triangle */
1294: triC1[2] = arr[0 * 2];
1295: triC1[7] = arr[1 * 2];
1296: triC1[11] = arr[0 * 2];
1297: triC1[15] = arr[1 * 2];
1298: triC1[22] = arr[1 * 2];
1299: triC1[26] = arr[2 * 2];
1300: *Nt = 2;
1301: *target = triT1;
1302: *size = triS1;
1303: *cone = triC1;
1304: *ornt = triO1;
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: /* Add 2 edges inside this triangle, making 3 new triangles.
1309: RT_TRIANGLE_SPLIT_12
1310: 2
1311: |\
1312: | \
1313: | \
1314: 0 \
1315: | 1
1316: | \
1317: | B \
1318: 2-------1
1319: | C / \
1320: 1 ____/ 0
1321: |/ A \
1322: 0-----0-----1
1323: RT_TRIANGLE_SPLIT_10
1324: 2
1325: |\
1326: | \
1327: | \
1328: 0 \
1329: | 1
1330: | \
1331: | A \
1332: 2 1
1333: | /|\
1334: 1 ____/ / 0
1335: |/ C / B \
1336: 0-----0-----1
1337: RT_TRIANGLE_SPLIT_20
1338: 2
1339: |\
1340: | \
1341: | \
1342: 0 \
1343: | \
1344: | \
1345: | \
1346: 2 A 1
1347: |\ \
1348: 1 ---\ \
1349: |B \_C----\\
1350: 0-----0-----1
1351: RT_TRIANGLE_SPLIT_21
1352: 2
1353: |\
1354: | \
1355: | \
1356: 0 \
1357: | \
1358: | B \
1359: | \
1360: 2-------1
1361: |\ C \
1362: 1 ---\ \
1363: | A ----\\
1364: 0-----0-----1
1365: RT_TRIANGLE_SPLIT_01
1366: 2
1367: |\
1368: |\\
1369: || \
1370: | \ \
1371: | | \
1372: | | \
1373: | | \
1374: 2 \ C 1
1375: | A | / \
1376: | | |B \
1377: | \/ \
1378: 0-----0-----1
1379: RT_TRIANGLE_SPLIT_02
1380: 2
1381: |\
1382: |\\
1383: || \
1384: | \ \
1385: | | \
1386: | | \
1387: | | \
1388: 2 C \ 1
1389: |\ | \
1390: | \__| A \
1391: | B \\ \
1392: 0-----0-----1
1393: */
1394: static PetscErrorCode SBRGetTriangleSplitDouble(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1395: {
1396: PetscInt e0, e1;
1397: const PetscInt *arr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, o);
1398: static DMPolytopeType triT2[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
1399: static PetscInt triS2[] = {2, 3};
1400: static PetscInt triC2[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1};
1401: static PetscInt triO2[] = {0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0};
1403: PetscFunctionBeginHot;
1404: /* To get the other divisions, we reorient the triangle */
1405: triC2[2] = arr[0 * 2];
1406: triC2[3] = arr[0 * 2 + 1] ? 1 : 0;
1407: triC2[7] = arr[1 * 2];
1408: triC2[11] = arr[1 * 2];
1409: triC2[15] = arr[2 * 2];
1410: /* Swap the first two edges if the triangle is reversed */
1411: e0 = o < 0 ? 23 : 19;
1412: e1 = o < 0 ? 19 : 23;
1413: triC2[e0] = arr[0 * 2];
1414: triC2[e0 + 1] = 0;
1415: triC2[e1] = arr[1 * 2];
1416: triC2[e1 + 1] = o < 0 ? 1 : 0;
1417: triO2[6] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]);
1418: /* Swap the first two edges if the triangle is reversed */
1419: e0 = o < 0 ? 34 : 30;
1420: e1 = o < 0 ? 30 : 34;
1421: triC2[e0] = arr[1 * 2];
1422: triC2[e0 + 1] = o < 0 ? 0 : 1;
1423: triC2[e1] = arr[2 * 2];
1424: triC2[e1 + 1] = o < 0 ? 1 : 0;
1425: triO2[9] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]);
1426: /* Swap the last two edges if the triangle is reversed */
1427: triC2[41] = arr[2 * 2];
1428: triC2[42] = o < 0 ? 0 : 1;
1429: triC2[45] = o < 0 ? 1 : 0;
1430: triC2[48] = o < 0 ? 0 : 1;
1431: triO2[11] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[1 * 2 + 1]);
1432: triO2[12] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[2 * 2 + 1]);
1433: *Nt = 2;
1434: *target = triT2;
1435: *size = triS2;
1436: *cone = triC2;
1437: *ornt = triO2;
1438: PetscFunctionReturn(PETSC_SUCCESS);
1439: }
1441: /* Add 3 edges inside this triangle, making 4 new triangles fanned around the midpoint of the
1442: first-bisected edge (the 4T partition of Rivara), e.g. for first = 1:
1443: 2
1444: |\
1445: | \
1446: 2 1
1447: | \ \
1448: | \ \
1449: | \ \
1450: 2 D \ B 1
1451: |\ \ |\
1452: | \ C \ | 0
1453: | 2--__\|A \
1454: 0-----0-----1
1455: This subdivision is used for a fully marked triangle below a 3D cell, where it matches the
1456: subdivision induced on the triangle by bisecting the cells above it (Lemma 1.2 of Plaza & Carey
1457: (2000)); a fully marked maximal (2D) triangle is subdivided regularly instead. */
1458: static PetscErrorCode SBRGetTriangleSplitFan_Private(PetscInt first, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1459: {
1460: static DMPolytopeType fanT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
1461: static PetscInt fanS[] = {3, 4};
1462: static PetscInt fanC[96];
1463: static PetscInt fanO[24];
1464: SBRTriSubdiv sub;
1465: PetscInt coff = 0, ooff = 0;
1467: PetscFunctionBeginHot;
1468: PetscCall(SBRTriangleSubdiv_Private(3, first, -1, &sub));
1469: for (PetscInt s = 0; s < sub.nseg; ++s) {
1470: for (PetscInt i = 0; i < 2; ++i) {
1471: const PetscInt sym = sub.seg[s][i];
1473: if (sym < 3) {
1474: const PetscInt acp[2] = {sym, 0}; /* vertex v is the first vertex of edge v */
1476: PetscCall(SBRAppendCone_Private(fanC, &coff, 96, fanO, &ooff, 24, DM_POLYTOPE_POINT, 2, acp, 0, 0));
1477: } else {
1478: const PetscInt acp[1] = {sym - 3};
1480: PetscCall(SBRAppendCone_Private(fanC, &coff, 96, fanO, &ooff, 24, DM_POLYTOPE_POINT, 1, acp, 0, 0));
1481: }
1482: }
1483: }
1484: for (PetscInt t = 0; t < sub.ntri; ++t) {
1485: for (PetscInt i = 0; i < 3; ++i) {
1486: const PetscInt P = sub.tri[t][i], Q = sub.tri[t][(i + 1) % 3];
1487: const PetscInt v = P < 3 ? P : Q, e = (P < 3 ? Q : P) - 3;
1488: PetscInt r;
1490: if ((P < 3 || Q < 3) && (triEdgeVert[e][0] == v || triEdgeVert[e][1] == v)) {
1491: /* Half of a split edge: replica 0 runs from the edge tail to the midpoint, replica 1 from the midpoint to the head */
1492: const PetscInt acp[1] = {e};
1494: r = triEdgeVert[e][0] == v ? 0 : 1;
1495: PetscCall(SBRAppendCone_Private(fanC, &coff, 96, fanO, &ooff, 24, DM_POLYTOPE_SEGMENT, 1, acp, r, (r == 0) == (P == v) ? 0 : -1));
1496: } else {
1497: for (r = 0; r < sub.nseg; ++r)
1498: if ((sub.seg[r][0] == P && sub.seg[r][1] == Q) || (sub.seg[r][0] == Q && sub.seg[r][1] == P)) break;
1499: PetscCheck(r < sub.nseg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior segment (%" PetscInt_FMT ", %" PetscInt_FMT ") not found", P, Q);
1500: PetscCall(SBRAppendCone_Private(fanC, &coff, 96, fanO, &ooff, 24, DM_POLYTOPE_SEGMENT, 0, NULL, r, sub.seg[r][0] == P ? 0 : -1));
1501: }
1502: }
1503: }
1504: *Nt = 2;
1505: *target = fanT;
1506: *size = fanS;
1507: *cone = fanC;
1508: *ornt = fanO;
1509: PetscFunctionReturn(PETSC_SUCCESS);
1510: }
1512: static PetscErrorCode DMPlexTransformCellTransform_SBR(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1513: {
1514: DMLabel trType = tr->trType;
1515: PetscInt val;
1517: PetscFunctionBeginHot;
1518: PetscCheck(p >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point argument is invalid");
1519: PetscCall(DMLabelGetValue(trType, p, &val));
1520: if (rt) *rt = val;
1521: switch (source) {
1522: case DM_POLYTOPE_POINT:
1523: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1524: case DM_POLYTOPE_QUADRILATERAL:
1525: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1526: case DM_POLYTOPE_HEXAHEDRON:
1527: case DM_POLYTOPE_TRI_PRISM:
1528: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1529: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1530: case DM_POLYTOPE_PYRAMID:
1531: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
1532: break;
1533: case DM_POLYTOPE_SEGMENT:
1534: if (val == RT_EDGE) PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
1535: else PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt));
1536: break;
1537: case DM_POLYTOPE_TRIANGLE:
1538: switch (val) {
1539: case RT_TRIANGLE_SPLIT_0:
1540: PetscCall(SBRGetTriangleSplitSingle(2, Nt, target, size, cone, ornt));
1541: break;
1542: case RT_TRIANGLE_SPLIT_1:
1543: PetscCall(SBRGetTriangleSplitSingle(0, Nt, target, size, cone, ornt));
1544: break;
1545: case RT_TRIANGLE_SPLIT_2:
1546: PetscCall(SBRGetTriangleSplitSingle(1, Nt, target, size, cone, ornt));
1547: break;
1548: case RT_TRIANGLE_SPLIT_21:
1549: PetscCall(SBRGetTriangleSplitDouble(-3, Nt, target, size, cone, ornt));
1550: break;
1551: case RT_TRIANGLE_SPLIT_10:
1552: PetscCall(SBRGetTriangleSplitDouble(-2, Nt, target, size, cone, ornt));
1553: break;
1554: case RT_TRIANGLE_SPLIT_02:
1555: PetscCall(SBRGetTriangleSplitDouble(-1, Nt, target, size, cone, ornt));
1556: break;
1557: case RT_TRIANGLE_SPLIT_12:
1558: PetscCall(SBRGetTriangleSplitDouble(0, Nt, target, size, cone, ornt));
1559: break;
1560: case RT_TRIANGLE_SPLIT_20:
1561: PetscCall(SBRGetTriangleSplitDouble(1, Nt, target, size, cone, ornt));
1562: break;
1563: case RT_TRIANGLE_SPLIT_01:
1564: PetscCall(SBRGetTriangleSplitDouble(2, Nt, target, size, cone, ornt));
1565: break;
1566: case RT_TRIANGLE_SPLIT_FAN_0:
1567: case RT_TRIANGLE_SPLIT_FAN_1:
1568: case RT_TRIANGLE_SPLIT_FAN_2:
1569: PetscCall(SBRGetTriangleSplitFan_Private(val - RT_TRIANGLE_SPLIT_FAN_0, Nt, target, size, cone, ornt));
1570: break;
1571: case RT_TRIANGLE_SPLIT:
1572: PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt));
1573: break;
1574: default:
1575: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
1576: }
1577: break;
1578: case DM_POLYTOPE_TETRAHEDRON:
1579: if (val == RT_TET) PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
1580: else {
1581: PetscCheck(val >= RT_TET_SPLIT_BASE, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid refinement type %" PetscInt_FMT " for tetrahedron %" PetscInt_FMT, val, p);
1582: PetscCall(SBRGetTetSplit_Private(tr, val - RT_TET_SPLIT_BASE, Nt, target, size, cone, ornt));
1583: }
1584: break;
1585: default:
1586: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1587: }
1588: PetscFunctionReturn(PETSC_SUCCESS);
1589: }
1591: static PetscErrorCode DMPlexTransformSetFromOptions_SBR(DMPlexTransform tr, PetscOptionItems PetscOptionsObject)
1592: {
1593: PetscInt cells[256], n = 256, i;
1594: PetscBool flg, validate = PETSC_FALSE;
1596: PetscFunctionBegin;
1597: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Options");
1598: PetscCall(PetscOptionsIntArray("-dm_plex_transform_sbr_ref_cell", "Mark cells for refinement", "", cells, &n, &flg));
1599: if (flg) {
1600: DMLabel active;
1602: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active));
1603: for (i = 0; i < n; ++i) PetscCall(DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE));
1604: PetscCall(DMPlexTransformSetActive(tr, active));
1605: PetscCall(DMLabelDestroy(&active));
1606: }
1607: PetscCall(PetscOptionsBool("-dm_plex_transform_sbr_validate", "Validate the tetrahedron subdivision generator against Plaza & Carey (2000), Table 4", "", validate, &validate, NULL));
1608: PetscOptionsHeadEnd();
1609: if (validate) PetscCall(DMPlexTransformSBRValidate_Private(tr));
1610: PetscFunctionReturn(PETSC_SUCCESS);
1611: }
1613: static PetscErrorCode DMPlexTransformView_SBR(DMPlexTransform tr, PetscViewer viewer)
1614: {
1615: PetscBool isascii;
1617: PetscFunctionBegin;
1620: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1621: if (isascii) {
1622: PetscViewerFormat format;
1623: const char *name;
1625: PetscCall(PetscObjectGetName((PetscObject)tr, &name));
1626: PetscCall(PetscViewerASCIIPrintf(viewer, "SBR refinement %s\n", name ? name : ""));
1627: PetscCall(PetscViewerGetFormat(viewer, &format));
1628: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(DMLabelView(tr->trType, viewer));
1629: } else {
1630: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
1631: }
1632: PetscFunctionReturn(PETSC_SUCCESS);
1633: }
1635: static PetscErrorCode DMPlexTransformDestroy_SBR(DMPlexTransform tr)
1636: {
1637: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
1639: PetscFunctionBegin;
1640: for (PetscInt i = 0; i < sbr->Ncache; ++i) {
1641: PetscCall(PetscFree(sbr->cacheTarget[i]));
1642: PetscCall(PetscFree(sbr->cacheSize[i]));
1643: PetscCall(PetscFree(sbr->cacheCone[i]));
1644: PetscCall(PetscFree(sbr->cacheOrnt[i]));
1645: }
1646: PetscCall(PetscFree6(sbr->cacheCode, sbr->cacheNt, sbr->cacheTarget, sbr->cacheSize, sbr->cacheCone, sbr->cacheOrnt));
1647: PetscCall(PetscFree(sbr->edgeLen));
1648: PetscCall(PetscSectionDestroy(&sbr->secEdgeLen));
1649: PetscCall(DMLabelDestroy(&sbr->splitPoints));
1650: PetscCall(PetscFree(tr->data));
1651: PetscFunctionReturn(PETSC_SUCCESS);
1652: }
1654: static PetscErrorCode DMPlexTransformInitialize_SBR(DMPlexTransform tr)
1655: {
1656: PetscFunctionBegin;
1657: tr->ops->view = DMPlexTransformView_SBR;
1658: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_SBR;
1659: tr->ops->setup = DMPlexTransformSetUp_SBR;
1660: tr->ops->destroy = DMPlexTransformDestroy_SBR;
1661: tr->ops->setdimensions = DMPlexTransformSetDimensions_Internal;
1662: tr->ops->celltransform = DMPlexTransformCellTransform_SBR;
1663: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_SBR;
1664: tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
1665: PetscFunctionReturn(PETSC_SUCCESS);
1666: }
1668: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform tr)
1669: {
1670: DMPlexRefine_SBR *f;
1672: PetscFunctionBegin;
1674: PetscCall(PetscNew(&f));
1675: tr->data = f;
1677: PetscCall(DMPlexTransformInitialize_SBR(tr));
1678: PetscCall(PetscCitationsRegister(SBRCitation, &SBRcite));
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }