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