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 DMPlexTransformDestroy_Alfeld(DMPlexTransform tr)
23: {
24: DMPlexRefine_Alfeld *f = (DMPlexRefine_Alfeld *)tr->data;
26: PetscFunctionBegin;
27: PetscCall(PetscFree(f));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
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};
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: }
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};
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: }
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: }
180: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform tr)
181: {
182: DMPlexRefine_Alfeld *f;
184: PetscFunctionBegin;
186: PetscCall(PetscNew(&f));
187: tr->data = f;
189: PetscCall(DMPlexTransformInitialize_Alfeld(tr));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }