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