1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/
2
DMPlexTransformView_ToBox(DMPlexTransform tr,PetscViewer viewer)3 static PetscErrorCode DMPlexTransformView_ToBox(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, "Transformation to box cells %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_ToBox(DMPlexTransform tr)22 static PetscErrorCode DMPlexTransformDestroy_ToBox(DMPlexTransform tr)
23 {
24 DMPlexRefine_ToBox *f = (DMPlexRefine_ToBox *)tr->data;
25
26 PetscFunctionBegin;
27 PetscCall(PetscFree(f));
28 PetscFunctionReturn(PETSC_SUCCESS);
29 }
30
DMPlexTransformGetSubcellOrientation_ToBox(DMPlexTransform tr,DMPolytopeType sct,PetscInt sp,PetscInt so,DMPolytopeType tct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)31 static PetscErrorCode DMPlexTransformGetSubcellOrientation_ToBox(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
32 {
33 PetscBool convertTensor = PETSC_TRUE;
34 static PetscInt tri_seg[] = {0, 0, 2, 0, 1, 0, 2, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0};
35 static PetscInt tri_quad[] = {1, -3, 0, -3, 2, -4, 0, -2, 2, -2, 1, -2, 2, -1, 1, -4, 0, -1, 0, 0, 1, 0, 2, 0, 1, 1, 2, 2, 0, 1, 2, 3, 0, 3, 1, 2};
36 static PetscInt tseg_seg[] = {0, -1, 0, 0, 0, 0, 0, -1};
37 static PetscInt tseg_quad[] = {1, 2, 0, 2, 1, -3, 0, -3, 0, 0, 1, 0, 0, -1, 1, -1};
38 static PetscInt tet_seg[] = {3, 0, 2, 0, 0, 0, 1, 0, 3, 0, 1, 0, 2, 0, 0, 0, 3, 0, 0, 0, 1, 0, 2, 0, 2, 0, 3, 0, 1, 0, 0, 0, 2, 0, 0, 0, 3, 0, 1, 0, 2, 0, 1, 0, 0, 0, 3, 0, 1, 0, 0, 0, 2, 0, 3, 0, 1, 0, 3, 0, 0, 0, 2, 0,
39 1, 0, 2, 0, 3, 0, 0, 0, 0, 0, 1, 0, 3, 0, 2, 0, 0, 0, 2, 0, 1, 0, 3, 0, 0, 0, 3, 0, 2, 0, 1, 0, 0, 0, 1, 0, 2, 0, 3, 0, 0, 0, 3, 0, 1, 0, 2, 0, 0, 0, 2, 0, 3, 0, 1, 0, 1, 0, 0, 0, 3, 0, 2, 0,
40 1, 0, 2, 0, 0, 0, 3, 0, 1, 0, 3, 0, 2, 0, 0, 0, 2, 0, 3, 0, 0, 0, 1, 0, 2, 0, 1, 0, 3, 0, 0, 0, 2, 0, 0, 0, 1, 0, 3, 0, 3, 0, 2, 0, 1, 0, 0, 0, 3, 0, 0, 0, 2, 0, 1, 0, 3, 0, 1, 0, 0, 0, 2, 0};
41 static PetscInt tet_quad[] = {2, 0, 5, -3, 4, 0, 0, 3, 3, 1, 1, 1, 0, 0, 3, 0, 5, 0, 1, 0, 4, -2, 2, 0, 1, 1, 4, -3, 3, -3, 2, 3, 5, -2, 0, 0, 3, 1, 5, 3, 0, 0, 4, 3, 2, 0, 1, -3, 4, 0, 2, 3, 5, -2, 1, -4, 0, -2,
42 3, 1, 1, -3, 0, -3, 2, -2, 3, 0, 5, 0, 4, 0, 2, -2, 1, -4, 0, -2, 4, -3, 3, -3, 5, 0, 4, -2, 3, -4, 1, 1, 5, 3, 0, 0, 2, -2, 5, 0, 0, 3, 3, 1, 2, -3, 1, -3, 4, -2, 3, -3, 1, 0, 4, -2, 0, -3,
43 2, -2, 5, -2, 0, -2, 2, -3, 1, -3, 5, -3, 4, 0, 3, -3, 5, -2, 4, 3, 2, 0, 3, -4, 1, 1, 0, -2, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 3, 1, 4, 3, 1, -3, 5, 3, 2, -2, 0, 0, 5, 0, 2, -3, 4, -2,
44 0, 3, 1, 1, 3, 1, 4, 0, 1, -4, 3, 1, 2, 3, 0, 0, 5, -2, 2, 0, 0, 3, 1, 1, 5, -3, 3, -3, 4, 0, 5, -2, 3, -4, 0, -2, 4, 3, 1, -3, 2, 0, 4, -2, 5, 3, 2, -2, 3, -4, 0, -2, 1, 1, 3, -3, 0, -3,
45 5, -2, 1, 0, 2, 0, 4, -2, 1, 1, 2, 3, 0, 0, 4, -3, 5, 0, 3, -3, 0, -2, 5, -3, 3, -3, 2, -3, 4, -2, 1, -3, 2, -2, 4, -3, 5, 0, 1, -4, 3, 1, 0, -2, 1, -3, 3, 0, 4, 0, 0, -3, 5, -2, 2, -2};
46 static PetscInt tet_hex[] = {2, -2, 3, -2, 1, -10, 0, -13, 3, -10, 1, -13, 2, -10, 0, -10, 1, -2, 2, -13, 3, -13, 0, -2, 3, -13, 2, -10, 0, -2, 1, -2, 2, -13, 0, -10, 3, -2, 1, -13, 0, -13, 3, -10, 2, -2, 1, -10,
47 0, -10, 1, -2, 3, -10, 2, -10, 1, -10, 3, -13, 0, -13, 2, -2, 3, -2, 0, -2, 1, -13, 2, -13, 1, -13, 0, -13, 2, -13, 3, -2, 0, -2, 2, -2, 1, -2, 3, -13, 2, -10, 1, -10, 0, -10, 3, -10,
48 0, 0, 1, 0, 2, 0, 3, 0, 1, 17, 2, 17, 0, 17, 3, 16, 2, 16, 0, 16, 1, 16, 3, 17, 1, 16, 0, 17, 3, 17, 2, 16, 0, 16, 3, 16, 1, 0, 2, 17, 3, 0, 1, 17, 0, 0, 2, 0,
49 2, 17, 3, 0, 0, 16, 1, 0, 3, 17, 0, 0, 2, 16, 1, 16, 0, 17, 2, 0, 3, 16, 1, 17, 3, 16, 2, 16, 1, 17, 0, 17, 2, 0, 1, 16, 3, 0, 0, 0, 1, 0, 3, 17, 2, 17, 0, 16};
50 static PetscInt trip_seg[] = {1, 0, 0, 0, 3, 0, 4, 0, 2, 0, 1, 0, 0, 0, 4, 0, 2, 0, 3, 0, 1, 0, 0, 0, 2, 0, 3, 0, 4, 0, 0, 0, 1, 0, 3, 0, 2, 0, 4, 0, 0, 0, 1, 0, 4, 0, 3, 0, 2, 0, 0, 0, 1, 0, 2, 0, 4, 0, 3, 0,
51 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 0, 0, 1, 0, 4, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 4, 0, 2, 0, 1, 0, 0, 0, 2, 0, 4, 0, 3, 0, 1, 0, 0, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 3, 0, 2, 0, 4, 0};
52 static PetscInt trip_quad[] = {1, 1, 2, 2, 0, 1, 7, -1, 8, -1, 6, -1, 4, -1, 5, -1, 3, -1, 2, 3, 0, 3, 1, 2, 8, -1, 6, -1, 7, -1, 5, -1, 3, -1, 4, -1, 0, 0, 1, 0, 2, 0, 6, -1, 7, -1, 8, -1, 3, -1, 4, -1, 5, -1,
53 2, -1, 1, -4, 0, -1, 4, 0, 3, 0, 5, 0, 7, 0, 6, 0, 8, 0, 0, -2, 2, -2, 1, -2, 5, 0, 4, 0, 3, 0, 8, 0, 7, 0, 6, 0, 1, -3, 0, -3, 2, -4, 3, 0, 5, 0, 4, 0, 6, 0, 8, 0, 7, 0,
54 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 8, 0, 2, 3, 0, 3, 1, 2, 5, 0, 3, 0, 4, 0, 8, 0, 6, 0, 7, 0, 1, 1, 2, 2, 0, 1, 4, 0, 5, 0, 3, 0, 7, 0, 8, 0, 6, 0,
55 1, -3, 0, -3, 2, -4, 6, -1, 8, -1, 7, -1, 3, -1, 5, -1, 4, -1, 0, -2, 2, -2, 1, -2, 8, -1, 7, -1, 6, -1, 5, -1, 4, -1, 3, -1, 2, -1, 1, -4, 0, -1, 7, -1, 6, -1, 8, -1, 4, -1, 3, -1, 5, -1};
56 static PetscInt trip_hex[] = {4, -12, 5, -6, 3, -12, 1, -12, 2, -6, 0, -12, 5, -11, 3, -11, 4, -6, 2, -11, 0, -11, 1, -6, 3, -9, 4, -9, 5, -9, 0, -9, 1, -9, 2, -9, 2, -3, 1, -4, 0, -3, 5, -3, 4, -4, 3, -3,
57 0, -2, 2, -2, 1, -2, 3, -2, 5, -2, 4, -2, 1, -1, 0, -1, 2, -4, 4, -1, 3, -1, 5, -4, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 2, 1, 0, 1, 1, 2, 5, 1, 3, 1, 4, 2,
58 1, 3, 2, 2, 0, 3, 4, 3, 5, 2, 3, 3, 4, 8, 3, 8, 5, 11, 1, 8, 0, 8, 2, 11, 3, 10, 5, 10, 4, 10, 0, 10, 2, 10, 1, 10, 5, 5, 4, 11, 3, 5, 2, 5, 1, 11, 0, 5};
59 static PetscInt ctrip_seg[] = {0, -1, 0, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, -1};
60 static PetscInt ctrip_quad[] = {0, -1, 2, -1, 1, -1, 2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, 0, 2, 0, 1, 0, 2, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0,
61 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0, 0, -1, 1, -1, 2, -1, 1, -1, 2, -1, 0, -1, 2, -1, 0, -1, 1, -1};
62 static PetscInt ctrip_hex[] = {1, 8, 0, 8, 2, 11, 0, 10, 2, 10, 1, 10, 2, 5, 1, 11, 0, 5, 1, -1, 0, -1, 2, -4, 0, -2, 2, -2, 1, -2, 2, -3, 1, -4, 0, -3,
63 0, 0, 1, 0, 2, 0, 1, 3, 2, 2, 0, 3, 2, 1, 0, 1, 1, 2, 0, -9, 1, -9, 2, -9, 1, -12, 2, -6, 0, -12, 2, -11, 0, -11, 1, -6};
64 static PetscInt tquadp_seg[] = {0, -1, 0, -1, 0, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, -1, 0, -1};
65 static PetscInt tquadp_quad[] = {1, -1, 0, -1, 3, -1, 2, -1, 0, -1, 3, -1, 2, -1, 1, -1, 3, -1, 2, -1, 1, -1, 0, -1, 2, -1, 1, -1, 0, -1, 3, -1, 1, 0, 0, 0, 3, 0, 2, 0, 0, 0, 3,
66 0, 2, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 2, 0, 1, 0, 0, 0, 3, 0, 0, 0, 1, 0, 2, 0, 3, 0, 1, 0, 2, 0, 3, 0, 0, 0, 2, 0, 3, 0, 0, 0,
67 1, 0, 3, 0, 0, 0, 1, 0, 2, 0, 0, -1, 1, -1, 2, -1, 3, -1, 1, -1, 2, -1, 3, -1, 0, -1, 2, -1, 3, -1, 0, -1, 1, -1, 3, -1, 0, -1, 1, -1, 2, -1};
68 static PetscInt tquadp_hex[] = {2, 11, 1, 11, 0, 11, 3, 11, 1, 8, 0, 8, 3, 8, 2, 8, 0, 10, 3, 10, 2, 10, 1, 10, 3, 5, 2, 5, 1, 5, 0, 5, 2, -4, 1, -4, 0, -4, 3, -4, 1, -1, 0,
69 -1, 3, -1, 2, -1, 0, -2, 3, -2, 2, -2, 1, -2, 3, -3, 2, -3, 1, -3, 0, -3, 0, 0, 1, 0, 2, 0, 3, 0, 1, 3, 2, 3, 3, 3, 0, 3, 2, 2, 3, 2, 0, 2,
70 1, 2, 3, 1, 0, 1, 1, 1, 2, 1, 0, -9, 1, -9, 2, -9, 3, -9, 1, -12, 2, -12, 3, -12, 0, -12, 2, -6, 3, -6, 0, -6, 1, -6, 3, -11, 0, -11, 1, -11, 2, -11};
71
72 PetscFunctionBeginHot;
73 *rnew = r;
74 *onew = o;
75 if (!so) PetscFunctionReturn(PETSC_SUCCESS);
76 if (convertTensor) {
77 switch (sct) {
78 case DM_POLYTOPE_POINT:
79 case DM_POLYTOPE_SEGMENT:
80 case DM_POLYTOPE_POINT_PRISM_TENSOR:
81 case DM_POLYTOPE_QUADRILATERAL:
82 case DM_POLYTOPE_HEXAHEDRON:
83 PetscCall(DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew));
84 break;
85 case DM_POLYTOPE_TRIANGLE:
86 switch (tct) {
87 case DM_POLYTOPE_POINT:
88 break;
89 case DM_POLYTOPE_SEGMENT:
90 *rnew = tri_seg[(so + 3) * 6 + r * 2];
91 *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so + 3) * 6 + r * 2 + 1]);
92 break;
93 case DM_POLYTOPE_QUADRILATERAL:
94 *rnew = tri_quad[(so + 3) * 6 + r * 2];
95 *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_quad[(so + 3) * 6 + r * 2 + 1]);
96 break;
97 default:
98 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
99 }
100 break;
101 case DM_POLYTOPE_SEG_PRISM_TENSOR:
102 switch (tct) {
103 case DM_POLYTOPE_SEGMENT:
104 case DM_POLYTOPE_POINT_PRISM_TENSOR:
105 *rnew = tseg_seg[(so + 2) * 2 + r * 2];
106 *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_seg[(so + 2) * 2 + r * 2 + 1]);
107 break;
108 case DM_POLYTOPE_QUADRILATERAL:
109 *rnew = tseg_quad[(so + 2) * 4 + r * 2];
110 *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_quad[(so + 2) * 4 + r * 2 + 1]);
111 break;
112 default:
113 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
114 }
115 break;
116 case DM_POLYTOPE_TETRAHEDRON:
117 switch (tct) {
118 case DM_POLYTOPE_POINT:
119 break;
120 case DM_POLYTOPE_SEGMENT:
121 *rnew = tet_seg[(so + 12) * 8 + r * 2];
122 *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so + 12) * 8 + r * 2 + 1]);
123 break;
124 case DM_POLYTOPE_QUADRILATERAL:
125 *rnew = tet_quad[(so + 12) * 12 + r * 2];
126 *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_quad[(so + 12) * 12 + r * 2 + 1]);
127 break;
128 case DM_POLYTOPE_HEXAHEDRON:
129 *rnew = tet_hex[(so + 12) * 8 + r * 2];
130 *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_hex[(so + 12) * 8 + r * 2 + 1]);
131 break;
132 default:
133 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
134 }
135 break;
136 case DM_POLYTOPE_TRI_PRISM:
137 switch (tct) {
138 case DM_POLYTOPE_POINT:
139 break;
140 case DM_POLYTOPE_SEGMENT:
141 *rnew = trip_seg[(so + 6) * 10 + r * 2];
142 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_seg[(so + 6) * 10 + r * 2 + 1]);
143 break;
144 case DM_POLYTOPE_QUADRILATERAL:
145 *rnew = trip_quad[(so + 6) * 18 + r * 2];
146 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_quad[(so + 6) * 18 + r * 2 + 1]);
147 break;
148 case DM_POLYTOPE_HEXAHEDRON:
149 *rnew = trip_hex[(so + 6) * 12 + r * 2];
150 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_hex[(so + 6) * 12 + r * 2 + 1]);
151 break;
152 default:
153 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
154 }
155 break;
156 case DM_POLYTOPE_TRI_PRISM_TENSOR:
157 switch (tct) {
158 case DM_POLYTOPE_SEGMENT:
159 *rnew = ctrip_seg[(so + 6) * 2 + r * 2];
160 *onew = DMPolytopeTypeComposeOrientation(tct, o, ctrip_seg[(so + 6) * 2 + r * 2 + 1]);
161 break;
162 case DM_POLYTOPE_QUADRILATERAL:
163 *rnew = ctrip_quad[(so + 6) * 6 + r * 2];
164 *onew = DMPolytopeTypeComposeOrientation(tct, o, ctrip_quad[(so + 6) * 6 + r * 2 + 1]);
165 break;
166 case DM_POLYTOPE_HEXAHEDRON:
167 *rnew = ctrip_hex[(so + 6) * 6 + r * 2];
168 *onew = DMPolytopeTypeComposeOrientation(tct, o, ctrip_hex[(so + 6) * 6 + r * 2 + 1]);
169 break;
170 default:
171 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
172 }
173 break;
174 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
175 switch (tct) {
176 case DM_POLYTOPE_SEGMENT:
177 *rnew = tquadp_seg[(so + 8) * 2 + r * 2];
178 *onew = DMPolytopeTypeComposeOrientation(tct, o, tquadp_seg[(so + 8) * 2 + r * 2 + 1]);
179 break;
180 case DM_POLYTOPE_QUADRILATERAL:
181 *rnew = tquadp_quad[(so + 8) * 8 + r * 2];
182 *onew = DMPolytopeTypeComposeOrientation(tct, o, tquadp_quad[(so + 8) * 8 + r * 2 + 1]);
183 break;
184 case DM_POLYTOPE_HEXAHEDRON:
185 *rnew = tquadp_hex[(so + 8) * 8 + r * 2];
186 *onew = DMPolytopeTypeComposeOrientation(tct, o, tquadp_hex[(so + 8) * 8 + r * 2 + 1]);
187 break;
188 default:
189 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
190 }
191 break;
192 default:
193 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
194 }
195 } else {
196 switch (sct) {
197 case DM_POLYTOPE_POINT:
198 case DM_POLYTOPE_SEGMENT:
199 case DM_POLYTOPE_POINT_PRISM_TENSOR:
200 case DM_POLYTOPE_QUADRILATERAL:
201 case DM_POLYTOPE_SEG_PRISM_TENSOR:
202 case DM_POLYTOPE_HEXAHEDRON:
203 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
204 PetscCall(DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew));
205 break;
206 default:
207 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
208 }
209 }
210 PetscFunctionReturn(PETSC_SUCCESS);
211 }
212
DMPlexTransformCellRefine_ToBox(DMPlexTransform tr,DMPolytopeType source,PetscInt p,PetscInt * rt,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])213 static PetscErrorCode DMPlexTransformCellRefine_ToBox(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
214 {
215 PetscBool convertTensor = PETSC_TRUE;
216 /* Change tensor edges to segments */
217 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_SEGMENT};
218 static PetscInt tedgeS[] = {1};
219 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
220 static PetscInt tedgeO[] = {0, 0};
221 /* Add 1 vertex, 3 edges inside every triangle, making 3 new quadrilaterals.
222 2
223 |\
224 | \
225 | \
226 | \
227 0 1
228 | \
229 | \
230 2 1
231 |\ / \
232 | 2 1 \
233 | \ / \
234 1 | 0
235 | 0 \
236 | | \
237 | | \
238 0-0-0-----1-----1
239 */
240 static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
241 static PetscInt triS[] = {1, 3, 3};
242 static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0};
243 static PetscInt triO[] = {0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, -1, 0, 0};
244 /* Add 1 edge inside every tensor quad, making 2 new quadrilaterals
245 2----2----1----3----3
246 | | |
247 | | |
248 | | |
249 4 A 6 B 5
250 | | |
251 | | |
252 | | |
253 0----0----0----1----1
254 */
255 static DMPolytopeType tsegT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
256 static PetscInt tsegS[] = {1, 2};
257 static PetscInt tsegC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
258 /* TODO Fix these */
259 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0};
260 static PetscInt tsegO[] = {0, 0, 0, 0, -1, -1, 0, 0, -1, -1};
261 /* Add 6 triangles inside every cell, making 4 new hexs
262 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
263 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]
264 called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
265 [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
266 We make a new hex in each corner
267 [v0, (e0, 0), (f0, 0), (e2, 0), (e3, 0), (f2, 0), (c0, 0), (f1, 0)]
268 [v1, (e4, 0), (f3, 0), (e1, 0), (e0, 0), (f0, 0), (c0, 0), (f1, 0)]
269 [v2, (e1, 0), (f3, 0), (e5, 0), (e2, 0), (f2, 0), (c0, 0), (f0, 0)]
270 [v3, (e4, 0), (f1, 0), (e3, 0), (e5, 0), (f2, 0), (c0, 0), (f3, 0)]
271 We create a new face for each edge
272 [(e3, 0), (f2, 0), (c0, 0), (f1, 0)]
273 [(f0, 0), (e0, 0), (f1, 0), (c0, 0)]
274 [(e2, 0), (f0, 0), (c0, 0), (f2, 0)]
275 [(f3, 0), (e4, 0), (f1, 0), (c0, 0)]
276 [(e1, 0), (f3, 0), (c0, 0), (f0, 0)]
277 [(e5, 0), (f3, 0), (c0, 0), (f2, 0)]
278 I could write a program to generate these from the first hex by acting with the symmetry group to take one subcell into another.
279 */
280 static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
281 static PetscInt tetS[] = {1, 4, 6, 4};
282 static PetscInt tetC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2};
283 static PetscInt tetO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, 0, -1, 0, 0, -1, -1, -1, 0, 0, -1, 0, 0, -1, -1, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 1, -3, 1, 0, 0, 3, 0, -2, 1, -3, 0, 3, 1, -2, 3, -4, -2, 3};
284 /* Add 3 quads inside every triangular prism, making 4 new prisms. */
285 static DMPolytopeType tripT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
286 static PetscInt tripS[] = {1, 5, 9, 6};
287 static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3};
288 static PetscInt tripO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, 0, -1, 0, -1, -1, 0, 0, 0, -1, -1, 0, 0, -1, -1, 0, 0, -1, -1, 0, -1, -1, 0, 0, -1, -1,
289 0, 0, -1, -1, 0, 0, 0, 0, -3, 0, 1, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, -3, 1, -2, 0, 0, -3, 0, 1, -2, 0, 0, 0, 0, -2, -2, 0, 0, 0, -3, 1};
290 /* Add 3 tensor quads inside every tensor triangular prism, making 4 new tensor triangular prisms.
291 2
292 |\
293 | \
294 | \
295 0---1
296
297 2
298
299 0 1
300
301 2
302 |\
303 | \
304 | \
305 0---1
306 */
307 static DMPolytopeType ttriT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
308 static PetscInt ttriS[] = {1, 3, 3};
309 static PetscInt ttriC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
310 static PetscInt ttriO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0};
311 /* TODO Add 3 quads inside every tensor triangular prism, making 4 new triangular prisms.
312 2
313 |\
314 | \
315 | \
316 0---1
317
318 2
319
320 0 1
321
322 2
323 |\
324 | \
325 | \
326 0---1
327 */
328 static DMPolytopeType ctripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
329 static PetscInt ctripS[] = {1, 3, 3};
330 static PetscInt ctripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
331 static PetscInt ctripO[] = {0, 0, 0, 0, -1, -1, 0, 0, -1, -1, 0, 0, -1, -1, -2, 0, 0, -3, 0, 1, -2, 0, 0, 0, 0, -2, -2, 0, 0, 0, -3, 1};
332 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
333 static PetscInt tquadpS[] = {1, 4, 4};
334 static PetscInt tquadpC[] = {
335 DM_POLYTOPE_POINT,
336 1,
337 0,
338 0,
339 DM_POLYTOPE_POINT,
340 1,
341 1,
342 0,
343 DM_POLYTOPE_SEGMENT,
344 1,
345 0,
346 0,
347 DM_POLYTOPE_SEGMENT,
348 0,
349 0,
350 DM_POLYTOPE_SEGMENT,
351 1,
352 1,
353 0,
354 DM_POLYTOPE_SEGMENT,
355 1,
356 2,
357 0,
358 DM_POLYTOPE_SEGMENT,
359 1,
360 0,
361 1,
362 DM_POLYTOPE_SEGMENT,
363 0,
364 0,
365 DM_POLYTOPE_SEGMENT,
366 1,
367 1,
368 1,
369 DM_POLYTOPE_SEGMENT,
370 1,
371 3,
372 0,
373 DM_POLYTOPE_SEGMENT,
374 1,
375 0,
376 2,
377 DM_POLYTOPE_SEGMENT,
378 0,
379 0,
380 DM_POLYTOPE_SEGMENT,
381 1,
382 1,
383 2,
384 DM_POLYTOPE_SEGMENT,
385 1,
386 4,
387 0,
388 DM_POLYTOPE_SEGMENT,
389 1,
390 0,
391 3,
392 DM_POLYTOPE_SEGMENT,
393 0,
394 0,
395 DM_POLYTOPE_SEGMENT,
396 1,
397 1,
398 3,
399 DM_POLYTOPE_SEGMENT,
400 1,
401 5,
402 0,
403 DM_POLYTOPE_QUADRILATERAL,
404 1,
405 0,
406 0,
407 DM_POLYTOPE_QUADRILATERAL,
408 1,
409 1,
410 0,
411 DM_POLYTOPE_QUADRILATERAL,
412 1,
413 2,
414 0,
415 DM_POLYTOPE_QUADRILATERAL,
416 0,
417 3,
418 DM_POLYTOPE_QUADRILATERAL,
419 0,
420 0,
421 DM_POLYTOPE_QUADRILATERAL,
422 1,
423 5,
424 1,
425 DM_POLYTOPE_QUADRILATERAL,
426 1,
427 0,
428 1,
429 DM_POLYTOPE_QUADRILATERAL,
430 1,
431 1,
432 1,
433 DM_POLYTOPE_QUADRILATERAL,
434 1,
435 2,
436 1,
437 DM_POLYTOPE_QUADRILATERAL,
438 0,
439 1,
440 DM_POLYTOPE_QUADRILATERAL,
441 1,
442 3,
443 0,
444 DM_POLYTOPE_QUADRILATERAL,
445 0,
446 0,
447 DM_POLYTOPE_QUADRILATERAL,
448 1,
449 0,
450 2,
451 DM_POLYTOPE_QUADRILATERAL,
452 1,
453 1,
454 2,
455 DM_POLYTOPE_QUADRILATERAL,
456 0,
457 1,
458 DM_POLYTOPE_QUADRILATERAL,
459 1,
460 4,
461 0,
462 DM_POLYTOPE_QUADRILATERAL,
463 1,
464 3,
465 1,
466 DM_POLYTOPE_QUADRILATERAL,
467 0,
468 2,
469 DM_POLYTOPE_QUADRILATERAL,
470 1,
471 0,
472 3,
473 DM_POLYTOPE_QUADRILATERAL,
474 1,
475 1,
476 3,
477 DM_POLYTOPE_QUADRILATERAL,
478 0,
479 3,
480 DM_POLYTOPE_QUADRILATERAL,
481 1,
482 4,
483 1,
484 DM_POLYTOPE_QUADRILATERAL,
485 0,
486 2,
487 DM_POLYTOPE_QUADRILATERAL,
488 1,
489 5,
490 0,
491 };
492 static PetscInt tquadpO[] = {0, 0, 0, 0, -1, -1, 0, 0, -1, -1, 0, 0, -1, -1, 0, 0, -1, -1, -2, 0, 0, -3, 0, 1, -2, 0, 0, 0, 0, -2, -2, 0, -3, 0, 0, 1, -2, 0, 0, 0, -3, 1};
493
494 PetscFunctionBeginHot;
495 if (rt) *rt = 0;
496 if (convertTensor) {
497 switch (source) {
498 case DM_POLYTOPE_POINT:
499 case DM_POLYTOPE_SEGMENT:
500 case DM_POLYTOPE_QUADRILATERAL:
501 case DM_POLYTOPE_HEXAHEDRON:
502 PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, rt, Nt, target, size, cone, ornt));
503 break;
504 case DM_POLYTOPE_POINT_PRISM_TENSOR:
505 *Nt = 1;
506 *target = tedgeT;
507 *size = tedgeS;
508 *cone = tedgeC;
509 *ornt = tedgeO;
510 break;
511 case DM_POLYTOPE_SEG_PRISM_TENSOR:
512 *Nt = 2;
513 *target = tsegT;
514 *size = tsegS;
515 *cone = tsegC;
516 *ornt = tsegO;
517 break;
518 case DM_POLYTOPE_TRI_PRISM_TENSOR:
519 *Nt = 3;
520 *target = ctripT;
521 *size = ctripS;
522 *cone = ctripC;
523 *ornt = ctripO;
524 break;
525 case DM_POLYTOPE_TRIANGLE:
526 *Nt = 3;
527 *target = triT;
528 *size = triS;
529 *cone = triC;
530 *ornt = triO;
531 break;
532 case DM_POLYTOPE_TETRAHEDRON:
533 *Nt = 4;
534 *target = tetT;
535 *size = tetS;
536 *cone = tetC;
537 *ornt = tetO;
538 break;
539 case DM_POLYTOPE_TRI_PRISM:
540 *Nt = 4;
541 *target = tripT;
542 *size = tripS;
543 *cone = tripC;
544 *ornt = tripO;
545 break;
546 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
547 *Nt = 3;
548 *target = tquadpT;
549 *size = tquadpS;
550 *cone = tquadpC;
551 *ornt = tquadpO;
552 break;
553 case DM_POLYTOPE_PYRAMID:
554 *Nt = 0;
555 *target = NULL;
556 *size = NULL;
557 *cone = NULL;
558 *ornt = NULL;
559 break;
560 default:
561 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
562 }
563 } else {
564 switch (source) {
565 case DM_POLYTOPE_POINT:
566 case DM_POLYTOPE_POINT_PRISM_TENSOR:
567 case DM_POLYTOPE_SEGMENT:
568 case DM_POLYTOPE_QUADRILATERAL:
569 case DM_POLYTOPE_SEG_PRISM_TENSOR:
570 case DM_POLYTOPE_HEXAHEDRON:
571 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
572 PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, rt, Nt, target, size, cone, ornt));
573 break;
574 case DM_POLYTOPE_TRIANGLE:
575 *Nt = 3;
576 *target = triT;
577 *size = triS;
578 *cone = triC;
579 *ornt = triO;
580 break;
581 case DM_POLYTOPE_TETRAHEDRON:
582 *Nt = 4;
583 *target = tetT;
584 *size = tetS;
585 *cone = tetC;
586 *ornt = tetO;
587 break;
588 case DM_POLYTOPE_TRI_PRISM:
589 *Nt = 4;
590 *target = tripT;
591 *size = tripS;
592 *cone = tripC;
593 *ornt = tripO;
594 break;
595 case DM_POLYTOPE_TRI_PRISM_TENSOR:
596 *Nt = 3;
597 *target = ttriT;
598 *size = ttriS;
599 *cone = ttriC;
600 *ornt = ttriO;
601 break;
602 case DM_POLYTOPE_PYRAMID:
603 *Nt = 0;
604 *target = NULL;
605 *size = NULL;
606 *cone = NULL;
607 *ornt = NULL;
608 break;
609 default:
610 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
611 }
612 }
613 PetscFunctionReturn(PETSC_SUCCESS);
614 }
615
DMPlexTransformInitialize_ToBox(DMPlexTransform tr)616 static PetscErrorCode DMPlexTransformInitialize_ToBox(DMPlexTransform tr)
617 {
618 PetscFunctionBegin;
619 tr->ops->view = DMPlexTransformView_ToBox;
620 tr->ops->destroy = DMPlexTransformDestroy_ToBox;
621 tr->ops->setdimensions = DMPlexTransformSetDimensions_Internal;
622 tr->ops->celltransform = DMPlexTransformCellRefine_ToBox;
623 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_ToBox;
624 tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
625 PetscFunctionReturn(PETSC_SUCCESS);
626 }
627
DMPlexTransformCreate_ToBox(DMPlexTransform tr)628 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform tr)
629 {
630 DMPlexRefine_ToBox *f;
631
632 PetscFunctionBegin;
633 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
634 PetscCall(PetscNew(&f));
635 tr->data = f;
636
637 PetscCall(DMPlexTransformInitialize_ToBox(tr));
638 PetscFunctionReturn(PETSC_SUCCESS);
639 }
640