Actual source code: plexrefalfeld.c
1: #include <petsc/private/dmplextransformimpl.h>
3: static PetscErrorCode DMPlexTransformView_Alfeld(DMPlexTransform tr, PetscViewer viewer)
4: {
5: PetscBool isascii;
7: PetscFunctionBegin;
10: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
11: if (isascii) {
12: const char *name;
14: PetscCall(PetscObjectGetName((PetscObject)tr, &name));
15: PetscCall(PetscViewerASCIIPrintf(viewer, "Alfeld refinement %s\n", name ? name : ""));
16: } else {
17: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
18: }
19: PetscFunctionReturn(PETSC_SUCCESS);
20: }
22: static PetscErrorCode DMPlexTransformDestroy_Alfeld(DMPlexTransform tr)
23: {
24: DMPlexRefine_Alfeld *f = (DMPlexRefine_Alfeld *)tr->data;
26: PetscFunctionBegin;
27: PetscCall(PetscFree(f));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Alfeld(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
32: {
33: DM dm;
34: PetscInt dim;
35: static PetscInt tri_seg[] = {1, 0, 0, 0, 2, 0, 0, 0, 2, 0, 1, 0, 2, 0, 1, 0, 0, 0, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0};
36: static PetscInt tri_tri[] = {0, -3, 2, -3, 1, -3, 2, -3, 1, -3, 0, -3, 1, -3, 0, -3, 2, -3, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0};
37: static PetscInt quad_seg[] = {
38: 2, 0, 1, 0, 0, 0, 3, 0, /* o = -4 */
39: 1, 0, 0, 0, 3, 0, 2, 0, /* o = -3 */
40: 0, 0, 3, 0, 2, 0, 1, 0, /* o = -2 */
41: 3, 0, 2, 0, 1, 0, 0, 0, /* o = -1 */
42: 0, 0, 1, 0, 2, 0, 3, 0, /* o = 0 */
43: 1, 0, 2, 0, 3, 0, 0, 0, /* o = 1 */
44: 2, 0, 3, 0, 0, 0, 1, 0, /* o = 2 */
45: 3, 0, 0, 0, 1, 0, 2, 0, /* o = 3 */
46: };
47: static PetscInt quad_tri[] = {
48: 1, -3, 0, -3, 3, -3, 2, -3, /* o = -4 */
49: 0, -3, 3, -3, 2, -3, 1, -3, /* o = -3 */
50: 3, -3, 2, -3, 1, -3, 0, -3, /* o = -2 */
51: 2, -3, 1, -3, 0, -3, 3, -3, /* o = -1 */
52: 0, 0, 1, 0, 2, 0, 3, 0, /* o = 0 */
53: 1, 0, 2, 0, 3, 0, 0, 0, /* o = 1 */
54: 2, 0, 3, 0, 0, 0, 1, 0, /* o = 2 */
55: 3, 0, 0, 0, 1, 0, 2, 0, /* o = 3 */
56: };
57: static PetscInt tet_seg[] = {2, 0, 3, 0, 1, 0, 0, 0, 3, 0, 1, 0, 2, 0, 0, 0, 1, 0, 2, 0, 3, 0, 0, 0, 3, 0, 2, 0, 0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 1, 0, 0, 0, 3, 0, 2, 0, 1, 0, 0, 0, 1, 0, 3, 0, 2, 0, 1, 0, 3, 0, 0, 0, 2, 0,
58: 3, 0, 0, 0, 1, 0, 2, 0, 1, 0, 0, 0, 2, 0, 3, 0, 0, 0, 2, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 3, 0, 0, 0, 1, 0, 2, 0, 3, 0, 1, 0, 2, 0, 0, 0, 3, 0, 2, 0, 0, 0, 1, 0, 3, 0, 1, 0, 0, 0, 3, 0, 2, 0,
59: 0, 0, 3, 0, 1, 0, 2, 0, 3, 0, 1, 0, 0, 0, 2, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 0, 0, 2, 0, 1, 0, 0, 0, 2, 0, 3, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 2, 0, 1, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 0, 0};
60: static PetscInt tet_tri[] = {5, 1, 2, 0, 4, 0, 1, 1, 3, 2, 0, -2, 4, 1, 5, 0, 2, 0, 3, -1, 0, 2, 1, 0, 2, 1, 4, 0, 5, 0, 0, -1, 1, -3, 3, -2, 5, -2, 3, 2, 1, 0, 4, 1, 2, 0, 0, 2, 1, 1, 5, -3, 3, 2, 2, -2, 0, -2,
61: 4, 0, 3, 0, 1, 0, 5, -3, 0, 0, 4, -3, 2, -3, 0, 0, 3, -2, 4, -3, 1, -2, 2, -3, 5, -3, 4, -2, 0, 2, 3, -2, 2, 1, 5, 0, 1, -3, 3, -1, 4, -3, 0, 2, 5, -2, 1, 0, 2, 0, 0, -1, 2, -3, 1, -3, 4, -2,
62: 3, -2, 5, 0, 1, -2, 0, -2, 2, -3, 3, 0, 5, -3, 4, -3, 2, -2, 1, -3, 0, -2, 5, 1, 4, 0, 3, 2, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 2, 1, 0, 2, 1, 0, 4, -2, 5, -3, 3, 2, 1, 1, 2, 0, 0, 2,
63: 5, 1, 3, -2, 4, -3, 0, -1, 4, 0, 3, 2, 2, 1, 1, 0, 5, -3, 3, 0, 0, -2, 4, 0, 1, -2, 5, 0, 2, 0, 4, 1, 3, 2, 0, -2, 5, -2, 2, -3, 1, -3, 5, 1, 1, -3, 3, -2, 2, -2, 4, -3, 0, 2, 3, -1, 5, 0,
64: 1, -3, 4, 1, 0, -2, 2, -3, 1, -2, 3, -2, 5, 0, 0, 0, 2, 0, 4, 0, 5, -2, 4, -3, 2, -3, 3, -1, 1, -3, 0, -2, 2, -2, 5, -3, 4, -3, 1, 1, 0, 2, 3, -2, 4, -2, 2, -3, 5, -3, 0, -1, 3, 2, 1, 0};
65: static PetscInt tet_tet[] = {3, -2, 2, -3, 0, -1, 1, -1, 3, -1, 1, -3, 2, -1, 0, -1, 3, -3, 0, -3, 1, -1, 2, -1, 2, -1, 3, -1, 1, -3, 0, -2, 2, -3, 0, -1, 3, -2, 1, -3, 2, -2, 1, -2, 0, -2, 3, -2,
66: 1, -2, 0, -2, 2, -2, 3, -1, 1, -1, 3, -3, 0, -3, 2, -2, 1, -3, 2, -1, 3, -1, 0, -3, 0, -3, 1, -1, 3, -3, 2, -3, 0, -2, 2, -2, 1, -2, 3, -3, 0, -1, 3, -2, 2, -3, 1, -2,
67: 0, 0, 1, 0, 2, 0, 3, 0, 0, 1, 3, 1, 1, 2, 2, 0, 0, 2, 2, 1, 3, 0, 1, 2, 1, 2, 0, 1, 3, 1, 2, 2, 1, 0, 2, 0, 0, 0, 3, 1, 1, 1, 3, 2, 2, 2, 0, 0,
68: 2, 1, 3, 0, 0, 2, 1, 0, 2, 2, 1, 1, 3, 2, 0, 2, 2, 0, 0, 0, 1, 0, 3, 2, 3, 2, 2, 2, 1, 1, 0, 1, 3, 0, 0, 2, 2, 1, 1, 1, 3, 1, 1, 2, 0, 1, 2, 1};
70: PetscFunctionBeginHot;
71: *rnew = r;
72: *onew = o;
73: if (!so) PetscFunctionReturn(PETSC_SUCCESS);
74: PetscCall(DMPlexTransformGetDM(tr, &dm));
75: PetscCall(DMGetDimension(dm, &dim));
76: if (dim == 2 && sct == DM_POLYTOPE_TRIANGLE) {
77: switch (tct) {
78: case DM_POLYTOPE_POINT:
79: break;
80: case DM_POLYTOPE_SEGMENT:
81: *rnew = tri_seg[(so + 3) * 6 + r * 2];
82: *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so + 3) * 6 + r * 2 + 1]);
83: break;
84: case DM_POLYTOPE_TRIANGLE:
85: *rnew = tri_tri[(so + 3) * 6 + r * 2];
86: *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so + 3) * 6 + r * 2 + 1]);
87: break;
88: default:
89: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
90: }
91: } else if (dim == 2 && sct == DM_POLYTOPE_QUADRILATERAL) {
92: switch (tct) {
93: case DM_POLYTOPE_POINT:
94: break;
95: case DM_POLYTOPE_SEGMENT:
96: *rnew = quad_seg[(so + 4) * 8 + r * 2];
97: *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_seg[(so + 4) * 8 + r * 2 + 1]);
98: break;
99: case DM_POLYTOPE_TRIANGLE:
100: *rnew = quad_tri[(so + 4) * 8 + r * 2];
101: *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_tri[(so + 4) * 8 + r * 2 + 1]);
102: break;
103: default:
104: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
105: }
106: } else if (dim == 3 && sct == DM_POLYTOPE_TETRAHEDRON) {
107: switch (tct) {
108: case DM_POLYTOPE_POINT:
109: break;
110: case DM_POLYTOPE_SEGMENT:
111: *rnew = tet_seg[(so + 12) * 8 + r * 2];
112: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so + 12) * 8 + r * 2 + 1]);
113: break;
114: case DM_POLYTOPE_TRIANGLE:
115: *rnew = tet_tri[(so + 12) * 12 + r * 2];
116: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so + 12) * 12 + r * 2 + 1]);
117: break;
118: case DM_POLYTOPE_TETRAHEDRON:
119: *rnew = tet_tet[(so + 12) * 8 + r * 2];
120: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so + 12) * 8 + r * 2 + 1]);
121: break;
122: default:
123: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
124: }
125: } else {
126: PetscCall(DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew));
127: }
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode DMPlexTransformCellRefine_Alfeld(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
132: {
133: DM dm;
134: PetscInt dim;
135: /* Add 1 vertex, 3 edges inside every triangle, making 3 new triangles.
136: 2
137: |\
138: |\\
139: | |\
140: | \ \
141: | | \
142: | \ \
143: | | \
144: 2 \ \
145: | | 1
146: | 2 \
147: | | \
148: | /\ \
149: | 0 1 |
150: | / \ |
151: |/ \|
152: 0---0----1
153: */
154: static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
155: static PetscInt triS[] = {1, 3, 3};
156: static PetscInt triC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2};
157: static PetscInt triO[] = {0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, -1};
158: /* Add 1 vertex and 4 edges inside every quad, making 4 new triangles.
159: 3---2---2 +-------+
160: | | |\ 2 /|
161: | | | 3 2 |
162: | | | \ / |
163: 3 1 --> | 3 x 1 |
164: | | | / \ |
165: | | | 0 1 |
166: | | |/ 0 \|
167: 0---0---1 +-------+
168: */
169: static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
170: static PetscInt quadS[] = {1, 4, 4};
171: static PetscInt quadC[] = {/* Cone of edge 0 */
172: DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 2, 0, 0, 0,
173: /* Cone of edge 1 */
174: DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 2, 1, 0, 0,
175: /* Cone of edge 2 */
176: DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 2, 2, 0, 0,
177: /* Cone of edge 3 */
178: DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 2, 3, 0, 0,
179: /* Cone of cell 0 */
180: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
181: /* Cone of cell 1 */
182: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1,
183: /* Cone of cell 2 */
184: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2,
185: /* Cone of cell 3 */
186: DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3};
187: static PetscInt quadO[] = {
188: 0, 0, /* Cone of edge 0 */
189: 0, 0, /* Cone of edge 1 */
190: 0, 0, /* Cone of edge 2 */
191: 0, 0, /* Cone of edge 3 */
192: 0, -1, 0, /* Cone of cell 0 */
193: 0, -1, 0, /* Cone of cell 1 */
194: 0, -1, 0, /* Cone of cell 2 */
195: 0, -1, 0, /* Cone of cell 3 */
196: };
197: /* Add 6 triangles inside every cell, making 4 new tets
198: The vertices of our reference tet are [(-1, -1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, 1)], which we call [v0, v1, v2, v3]. The first
199: three edges are [v0, v1], [v1, v2], [v2, v0] called e0, e1, and e2, and then three edges to the top point [v0, v3], [v1, v3], [v2, v3]
200: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
201: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
202: We make a new tet on each face
203: [v0, v1, v2, (c0, 0)]
204: [v0, v3, v1, (c0, 0)]
205: [v0, v2, v3, (c0, 0)]
206: [v2, v1, v3, (c0, 0)]
207: We create a new face for each edge
208: [v0, (c0, 0), v1 ]
209: [v0, v2, (c0, 0)]
210: [v2, v1, (c0, 0)]
211: [v0, (c0, 0), v3 ]
212: [v1, v3, (c0, 0)]
213: [v3, v2, (c0, 0)]
214: */
215: static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
216: static PetscInt tetS[] = {1, 4, 6, 4};
217: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 3, 0, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 3, 0, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 3, 0, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 3, 1, 0, 1, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 2, 0, 0, 0, DM_POLYTOPE_SEGMENT, 2, 0, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 2, 0, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 2, 1, 0, 0, DM_POLYTOPE_SEGMENT, 2, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 4};
218: static PetscInt tetO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, -1, -1, 0, -1, 0, -1, -1, -1, 0, -1, -1, 0, -1, 0, 0, 0, 0, 0, 0, -2, 0, 0, -2, -2, 0, 0, -2, -3, -3};
220: PetscFunctionBeginHot;
221: if (rt) *rt = 0;
222: PetscCall(DMPlexTransformGetDM(tr, &dm));
223: PetscCall(DMGetDimension(dm, &dim));
224: if (dim == 2 && source == DM_POLYTOPE_TRIANGLE) {
225: *Nt = 3;
226: *target = triT;
227: *size = triS;
228: *cone = triC;
229: *ornt = triO;
230: } else if (dim == 2 && source == DM_POLYTOPE_QUADRILATERAL) {
231: *Nt = 3;
232: *target = quadT;
233: *size = quadS;
234: *cone = quadC;
235: *ornt = quadO;
236: } else if (dim == 3 && source == DM_POLYTOPE_TETRAHEDRON) {
237: *Nt = 4;
238: *target = tetT;
239: *size = tetS;
240: *cone = tetC;
241: *ornt = tetO;
242: } else {
243: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, rt, Nt, target, size, cone, ornt));
244: }
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: static PetscErrorCode DMPlexTransformInitialize_Alfeld(DMPlexTransform tr)
249: {
250: PetscFunctionBegin;
251: tr->ops->view = DMPlexTransformView_Alfeld;
252: tr->ops->destroy = DMPlexTransformDestroy_Alfeld;
253: tr->ops->setdimensions = DMPlexTransformSetDimensions_Internal;
254: tr->ops->celltransform = DMPlexTransformCellRefine_Alfeld;
255: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Alfeld;
256: tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform tr)
261: {
262: DMPlexRefine_Alfeld *f;
264: PetscFunctionBegin;
266: PetscCall(PetscNew(&f));
267: tr->redFactor = 1.0;
268: tr->data = f;
270: PetscCall(DMPlexTransformInitialize_Alfeld(tr));
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }