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