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