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