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, 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};
39: 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};
40: 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,
41: 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,
42: 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};
43: 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,
44: 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,
45: 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,
46: 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,
47: 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};
48: 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,
49: 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,
50: 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,
51: 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};
54: *rnew = r;
55: *onew = o;
56: if (!so) return 0;
57: DMPlexTransformGetDM(tr, &dm);
58: DMGetDimension(dm, &dim);
59: if (dim == 2 && sct == DM_POLYTOPE_TRIANGLE) {
60: switch (tct) {
61: case DM_POLYTOPE_POINT:
62: break;
63: case DM_POLYTOPE_SEGMENT:
64: *rnew = tri_seg[(so + 3) * 6 + r * 2];
65: *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so + 3) * 6 + r * 2 + 1]);
66: break;
67: case DM_POLYTOPE_TRIANGLE:
68: *rnew = tri_tri[(so + 3) * 6 + r * 2];
69: *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so + 3) * 6 + r * 2 + 1]);
70: break;
71: default:
72: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
73: }
74: } else if (dim == 3 && sct == DM_POLYTOPE_TETRAHEDRON) {
75: switch (tct) {
76: case DM_POLYTOPE_POINT:
77: break;
78: case DM_POLYTOPE_SEGMENT:
79: *rnew = tet_seg[(so + 12) * 8 + r * 2];
80: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so + 12) * 8 + r * 2 + 1]);
81: break;
82: case DM_POLYTOPE_TRIANGLE:
83: *rnew = tet_tri[(so + 12) * 12 + r * 2];
84: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so + 12) * 12 + r * 2 + 1]);
85: break;
86: case DM_POLYTOPE_TETRAHEDRON:
87: *rnew = tet_tet[(so + 12) * 8 + r * 2];
88: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so + 12) * 8 + r * 2 + 1]);
89: break;
90: default:
91: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
92: }
93: } else {
94: DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
95: }
96: return 0;
97: }
99: static PetscErrorCode DMPlexTransformCellRefine_Alfeld(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
100: {
101: DM dm;
102: PetscInt dim;
103: /* Add 1 vertex, 3 edges inside every triangle, making 3 new triangles.
104: 2
105: |\
106: |\\
107: | |\
108: | \ \
109: | | \
110: | \ \
111: | | \
112: 2 \ \
113: | | 1
114: | 2 \
115: | | \
116: | /\ \
117: | 0 1 |
118: | / \ |
119: |/ \|
120: 0---0----1
121: */
122: static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
123: static PetscInt triS[] = {1, 3, 3};
124: 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};
125: static PetscInt triO[] = {0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, -1};
126: /* Add 6 triangles inside every cell, making 4 new tets
127: 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
128: 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]
129: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
130: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
131: We make a new tet on each face
132: [v0, v1, v2, (c0, 0)]
133: [v0, v3, v1, (c0, 0)]
134: [v0, v2, v3, (c0, 0)]
135: [v2, v1, v3, (c0, 0)]
136: We create a new face for each edge
137: [v0, (c0, 0), v1 ]
138: [v0, v2, (c0, 0)]
139: [v2, v1, (c0, 0)]
140: [v0, (c0, 0), v3 ]
141: [v1, v3, (c0, 0)]
142: [v3, v2, (c0, 0)]
143: */
144: static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
145: static PetscInt tetS[] = {1, 4, 6, 4};
146: 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};
147: 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};
150: if (rt) *rt = 0;
151: DMPlexTransformGetDM(tr, &dm);
152: DMGetDimension(dm, &dim);
153: if (dim == 2 && source == DM_POLYTOPE_TRIANGLE) {
154: *Nt = 3;
155: *target = triT;
156: *size = triS;
157: *cone = triC;
158: *ornt = triO;
159: } else if (dim == 3 && source == DM_POLYTOPE_TETRAHEDRON) {
160: *Nt = 4;
161: *target = tetT;
162: *size = tetS;
163: *cone = tetC;
164: *ornt = tetO;
165: } else {
166: DMPlexTransformCellTransformIdentity(tr, source, p, rt, Nt, target, size, cone, ornt);
167: }
168: return 0;
169: }
171: static PetscErrorCode DMPlexTransformInitialize_Alfeld(DMPlexTransform tr)
172: {
173: tr->ops->view = DMPlexTransformView_Alfeld;
174: tr->ops->setup = DMPlexTransformSetUp_Alfeld;
175: tr->ops->destroy = DMPlexTransformDestroy_Alfeld;
176: tr->ops->celltransform = DMPlexTransformCellRefine_Alfeld;
177: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Alfeld;
178: tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
179: return 0;
180: }
182: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform tr)
183: {
184: DMPlexRefine_Alfeld *f;
187: PetscNew(&f);
188: tr->data = f;
190: DMPlexTransformInitialize_Alfeld(tr);
191: return 0;
192: }