1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/ 2 3 /* Regular Refinment of Hybrid Meshes 4 5 We would like to express regular refinement as a small set of rules that can be applied on every point of the Plex 6 to automatically generate a refined Plex. In fact, we would like these rules to be general enough to encompass other 7 transformations, such as changing from one type of cell to another, as simplex to hex. 8 9 To start, we can create a function that takes an original cell type and returns the number of new cells replacing it 10 and the types of the new cells. 11 12 We need the group multiplication table for group actions from the dihedral group for each cell type. 13 14 We need an operator which takes in a cell, and produces a new set of cells with new faces and correct orientations. I think 15 we can just write this operator for faces with identity, and then compose the face orientation actions to get the actual 16 (face, orient) pairs for each subcell. 17 */ 18 19 /*@ 20 DMPlexRefineRegularGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell 21 22 Input Parameters: 23 + cr - The DMPlexCellRefiner object 24 - ct - The cell type 25 26 Output Parameters: 27 + Nf - The number of faces for this cell type 28 . v0 - The translation of the first vertex for each face 29 . J - The Jacobian for each face (map from original cell to subcell) 30 . invJ - The inverse Jacobian for each face 31 - detJ - The determinant of the Jacobian for each face 32 33 Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse. 34 $ v0 + j x_face = x_cell 35 $ invj (x_cell - v0) = x_face 36 37 Level: developer 38 39 .seealso: `DMPlexCellRefinerGetAffineTransforms()`, `Create()` 40 @*/ 41 PetscErrorCode DMPlexRefineRegularGetAffineFaceTransforms(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[]) { 42 /* 43 2 44 |\ 45 | \ 46 | \ 47 | \ 48 | \ 49 | \ 50 | \ 51 2 1 52 | \ 53 | \ 54 | \ 55 0---0-------1 56 v0[Nf][dc]: 3 x 2 57 J[Nf][df][dc]: 3 x 1 x 2 58 invJ[Nf][dc][df]: 3 x 2 x 1 59 detJ[Nf]: 3 60 */ 61 static PetscReal tri_v0[] = {0.0, -1.0, 0.0, 0.0, -1.0, 0.0}; 62 static PetscReal tri_J[] = {1.0, 0.0, -1.0, 1.0, 0.0, -1.0}; 63 static PetscReal tri_invJ[] = {1.0, 0.0, -0.5, 0.5, 0.0, -1.0}; 64 static PetscReal tri_detJ[] = {1.0, 1.414213562373095, 1.0}; 65 /* 66 3---------2---------2 67 | | 68 | | 69 | | 70 3 1 71 | | 72 | | 73 | | 74 0---------0---------1 75 76 v0[Nf][dc]: 4 x 2 77 J[Nf][df][dc]: 4 x 1 x 2 78 invJ[Nf][dc][df]: 4 x 2 x 1 79 detJ[Nf]: 4 80 */ 81 static PetscReal quad_v0[] = {0.0, -1.0, 1.0, 0.0, 0.0, 1.0 - 1.0, 0.0}; 82 static PetscReal quad_J[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0}; 83 static PetscReal quad_invJ[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0}; 84 static PetscReal quad_detJ[] = {1.0, 1.0, 1.0, 1.0}; 85 86 PetscFunctionBegin; 87 switch (ct) { 88 case DM_POLYTOPE_TRIANGLE: 89 if (Nf) *Nf = 3; 90 if (v0) *v0 = tri_v0; 91 if (J) *J = tri_J; 92 if (invJ) *invJ = tri_invJ; 93 if (detJ) *detJ = tri_detJ; 94 break; 95 case DM_POLYTOPE_QUADRILATERAL: 96 if (Nf) *Nf = 4; 97 if (v0) *v0 = quad_v0; 98 if (J) *J = quad_J; 99 if (invJ) *invJ = quad_invJ; 100 if (detJ) *detJ = quad_detJ; 101 break; 102 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]); 103 } 104 PetscFunctionReturn(0); 105 } 106 107 /*@ 108 DMPlexRefineRegularGetAffineTransforms - Gets the affine map from the reference cell to each subcell 109 110 Input Parameters: 111 + cr - The DMPlexCellRefiner object 112 - ct - The cell type 113 114 Output Parameters: 115 + Nc - The number of subcells produced from this cell type 116 . v0 - The translation of the first vertex for each subcell 117 . J - The Jacobian for each subcell (map from reference cell to subcell) 118 - invJ - The inverse Jacobian for each subcell 119 120 Level: developer 121 122 .seealso: `DMPlexRefineRegularGetAffineFaceTransforms()`, `DMPLEXREFINEREGULAR` 123 @*/ 124 PetscErrorCode DMPlexRefineRegularGetAffineTransforms(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[]) { 125 /* 126 2 127 |\ 128 | \ 129 | \ 130 | \ 131 | C \ 132 | \ 133 | \ 134 2---1---1 135 |\ D / \ 136 | 2 0 \ 137 |A \ / B \ 138 0---0-------1 139 */ 140 static PetscReal tri_v0[] = {-1.0, -1.0, 0.0, -1.0, -1.0, 0.0, 0.0, -1.0}; 141 static PetscReal tri_J[] = {0.5, 0.0, 0.0, 0.5, 142 143 0.5, 0.0, 0.0, 0.5, 144 145 0.5, 0.0, 0.0, 0.5, 146 147 0.0, -0.5, 0.5, 0.5}; 148 static PetscReal tri_invJ[] = {2.0, 0.0, 0.0, 2.0, 149 150 2.0, 0.0, 0.0, 2.0, 151 152 2.0, 0.0, 0.0, 2.0, 153 154 2.0, 2.0, -2.0, 0.0}; 155 /* 156 3---------2---------2 157 | | | 158 | D 2 C | 159 | | | 160 3----3----0----1----1 161 | | | 162 | A 0 B | 163 | | | 164 0---------0---------1 165 */ 166 static PetscReal quad_v0[] = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0}; 167 static PetscReal quad_J[] = {0.5, 0.0, 0.0, 0.5, 168 169 0.5, 0.0, 0.0, 0.5, 170 171 0.5, 0.0, 0.0, 0.5, 172 173 0.5, 0.0, 0.0, 0.5}; 174 static PetscReal quad_invJ[] = {2.0, 0.0, 0.0, 2.0, 175 176 2.0, 0.0, 0.0, 2.0, 177 178 2.0, 0.0, 0.0, 2.0, 179 180 2.0, 0.0, 0.0, 2.0}; 181 /* 182 Bottom (viewed from top) Top 183 1---------2---------2 7---------2---------6 184 | | | | | | 185 | B 2 C | | H 2 G | 186 | | | | | | 187 3----3----0----1----1 3----3----0----1----1 188 | | | | | | 189 | A 0 D | | E 0 F | 190 | | | | | | 191 0---------0---------3 4---------0---------5 192 */ 193 static PetscReal hex_v0[] = {-1.0, -1.0, -1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0}; 194 static PetscReal hex_J[] = {0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 195 196 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 197 198 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 199 200 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 201 202 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 203 204 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 205 206 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 207 208 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5}; 209 static PetscReal hex_invJ[] = {2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 210 211 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 212 213 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 214 215 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 216 217 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 218 219 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 220 221 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 222 223 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0}; 224 225 PetscFunctionBegin; 226 switch (ct) { 227 case DM_POLYTOPE_TRIANGLE: 228 if (Nc) *Nc = 4; 229 if (v0) *v0 = tri_v0; 230 if (J) *J = tri_J; 231 if (invJ) *invJ = tri_invJ; 232 break; 233 case DM_POLYTOPE_QUADRILATERAL: 234 if (Nc) *Nc = 4; 235 if (v0) *v0 = quad_v0; 236 if (J) *J = quad_J; 237 if (invJ) *invJ = quad_invJ; 238 break; 239 case DM_POLYTOPE_HEXAHEDRON: 240 if (Nc) *Nc = 8; 241 if (v0) *v0 = hex_v0; 242 if (J) *J = hex_J; 243 if (invJ) *invJ = hex_invJ; 244 break; 245 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]); 246 } 247 PetscFunctionReturn(0); 248 } 249 250 #if 0 251 static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[]) 252 { 253 static PetscReal seg_v[] = {-1.0, 0.0, 1.0}; 254 static PetscReal tri_v[] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0}; 255 static PetscReal quad_v[] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0}; 256 static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0, 257 -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0, 258 -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0}; 259 static PetscReal hex_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0, 260 -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, -1.0, 261 -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, 1.0, -1.0, 262 -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 1.0, -1.0, 0.0, 263 -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 264 -1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 265 -1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, 1.0, 266 -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 267 -1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0}; 268 269 PetscFunctionBegin; 270 switch (ct) { 271 case DM_POLYTOPE_SEGMENT: *Nv = 3; *subcellV = seg_v; break; 272 case DM_POLYTOPE_TRIANGLE: *Nv = 6; *subcellV = tri_v; break; 273 case DM_POLYTOPE_QUADRILATERAL: *Nv = 9; *subcellV = quad_v; break; 274 case DM_POLYTOPE_TETRAHEDRON: *Nv = 10; *subcellV = tet_v; break; 275 case DM_POLYTOPE_HEXAHEDRON: *Nv = 27; *subcellV = hex_v; break; 276 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]); 277 } 278 PetscFunctionReturn(0); 279 } 280 281 static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[]) 282 { 283 static PetscInt seg_v[] = {0, 1, 1, 2}; 284 static PetscInt tri_v[] = {0, 3, 5, 3, 1, 4, 5, 4, 2, 3, 4, 5}; 285 static PetscInt quad_v[] = {0, 4, 8, 7, 4, 1, 5, 8, 8, 5, 2, 6, 7, 8, 6, 3}; 286 static PetscInt tet_v[] = {0, 3, 1, 6, 3, 2, 4, 8, 1, 4, 5, 7, 6, 8, 7, 9, 287 1, 6, 3, 7, 8, 4, 3, 7, 7, 3, 1, 4, 7, 3, 8, 6}; 288 static PetscInt hex_v[] = {0, 3, 4, 1, 9, 10, 13, 12, 3, 6, 7, 4, 12, 13, 16, 15, 4, 7, 8, 5, 13, 14, 17, 16, 1, 4 , 5 , 2, 10, 11, 14, 13, 289 9, 12, 13, 10, 18, 19, 22, 21, 10, 13, 14, 11, 19, 20, 23, 22, 13, 16, 17, 14, 22, 23, 26, 25, 12, 15, 16, 13, 21, 22, 25, 24}; 290 291 PetscFunctionBegin; 292 PetscCheck(ct == rct,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 293 switch (ct) { 294 case DM_POLYTOPE_SEGMENT: *Nv = 2; *subcellV = &seg_v[r*(*Nv)]; break; 295 case DM_POLYTOPE_TRIANGLE: *Nv = 3; *subcellV = &tri_v[r*(*Nv)]; break; 296 case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break; 297 case DM_POLYTOPE_TETRAHEDRON: *Nv = 4; *subcellV = &tet_v[r*(*Nv)]; break; 298 case DM_POLYTOPE_HEXAHEDRON: *Nv = 8; *subcellV = &hex_v[r*(*Nv)]; break; 299 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]); 300 } 301 PetscFunctionReturn(0); 302 } 303 #endif 304 305 static PetscErrorCode DMPlexTransformView_Regular(DMPlexTransform tr, PetscViewer viewer) { 306 PetscBool isascii; 307 308 PetscFunctionBegin; 309 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 310 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 311 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 312 if (isascii) { 313 const char *name; 314 315 PetscCall(PetscObjectGetName((PetscObject)tr, &name)); 316 PetscCall(PetscViewerASCIIPrintf(viewer, "Regular refinement %s\n", name ? name : "")); 317 } else { 318 SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name); 319 } 320 PetscFunctionReturn(0); 321 } 322 323 static PetscErrorCode DMPlexTransformSetUp_Regular(DMPlexTransform tr) { 324 PetscFunctionBegin; 325 PetscFunctionReturn(0); 326 } 327 328 static PetscErrorCode DMPlexTransformDestroy_Regular(DMPlexTransform tr) { 329 DMPlexRefine_Regular *f = (DMPlexRefine_Regular *)tr->data; 330 331 PetscFunctionBegin; 332 PetscCall(PetscFree(f)); 333 PetscFunctionReturn(0); 334 } 335 336 PetscErrorCode DMPlexTransformGetSubcellOrientation_Regular(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) { 337 static PetscInt seg_seg[] = {1, -1, 0, -1, 0, 0, 1, 0}; 338 static PetscInt tri_seg[] = {2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, -1, 2, -1, 1, -1, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0}; 339 static PetscInt tri_tri[] = {1, -3, 0, -3, 2, -3, 3, -2, 0, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 3, -3, 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 0, 1, 3, 1, 2, 2, 0, 2, 1, 2, 3, 2}; 340 static PetscInt quad_seg[] = {1, 0, 0, 0, 3, 0, 2, 0, 0, 0, 3, 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, 1, 0, 3, 0, 0, 0, 1, 0, 2, 0}; 341 static PetscInt quad_quad[] = {2, -4, 1, -4, 0, -4, 3, -4, 1, -3, 0, -3, 3, -3, 2, -3, 0, -2, 3, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 3, 1, 0, 1, 2, 2, 3, 2, 0, 2, 1, 2, 3, 3, 0, 3, 1, 3, 2, 3}; 342 static PetscInt tseg_seg[] = {0, -1, 0, 0, 0, 0, 0, -1}; 343 static PetscInt tseg_tseg[] = {1, -2, 0, -2, 1, -1, 0, -1, 0, 0, 1, 0, 0, 1, 1, 1}; 344 static PetscInt tet_seg[] = {0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0}; 345 static PetscInt tet_tri[] = {2, -1, 3, -1, 1, -3, 0, -2, 6, 1, 7, -3, 5, 2, 4, -3, 3, -1, 1, -1, 2, -3, 0, -1, 4, 0, 5, 0, 6, 0, 7, 0, 1, -1, 2, -1, 3, -3, 0, -3, 4, 0, 5, 0, 6, 0, 7, 0, 3, -2, 2, -3, 0, -1, 1, -1, 7, -3, 6, 1, 4, 2, 5, -3, 346 2, -3, 0, -2, 3, -1, 1, -3, 4, 0, 5, 0, 6, 0, 7, 0, 0, -2, 3, -2, 2, -2, 1, -2, 4, 0, 5, 0, 6, 0, 7, 0, 0, -1, 1, -2, 3, -2, 2, -2, 7, 1, 6, -3, 5, -3, 4, 2, 1, -2, 3, -3, 0, -3, 2, -1, 4, 0, 5, 0, 6, 0, 7, 0, 347 3, -3, 0, -1, 1, -1, 2, -3, 4, 0, 5, 0, 6, 0, 7, 0, 1, -3, 0, -3, 2, -1, 3, -3, 6, -3, 7, 1, 4, -3, 5, 2, 0, -3, 2, -2, 1, -2, 3, -2, 4, 0, 5, 0, 6, 0, 7, 0, 2, -2, 1, -3, 0, -2, 3, -1, 4, 0, 5, 0, 6, 0, 7, 0, 348 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 1, 0, 2, 2, 0, 1, 3, 1, 4, 0, 5, 0, 6, 0, 7, 0, 2, 2, 0, 0, 1, 1, 3, 2, 4, 0, 5, 0, 6, 0, 7, 0, 1, 2, 0, 1, 3, 1, 2, 2, 5, 0, 4, 0, 7, -1, 6, -1, 349 0, 1, 3, 0, 1, 0, 2, 0, 4, 0, 5, 0, 6, 0, 7, 0, 3, 0, 1, 2, 0, 2, 2, 1, 4, 0, 5, 0, 6, 0, 7, 0, 2, 0, 3, 2, 0, 0, 1, 1, 4, -2, 5, -2, 7, 0, 6, 0, 3, 2, 0, 2, 2, 1, 1, 2, 4, 0, 5, 0, 6, 0, 7, 0, 350 0, 2, 2, 0, 3, 0, 1, 0, 4, 0, 5, 0, 6, 0, 7, 0, 3, 1, 2, 1, 1, 2, 0, 2, 5, -2, 4, -2, 6, -1, 7, -1, 2, 1, 1, 1, 3, 2, 0, 0, 4, 0, 5, 0, 6, 0, 7, 0, 1, 1, 3, 1, 2, 2, 0, 1, 4, 0, 5, 0, 6, 0, 7, 0}; 351 static PetscInt tet_tet[] = {2, -12, 3, -12, 1, -12, 0, -12, 6, -9, 7, -9, 5, -12, 4, -12, 3, -11, 1, -11, 2, -11, 0, -11, 4, 0, 5, 0, 6, 0, 7, 0, 1, -10, 2, -10, 3, -10, 0, -10, 4, 0, 5, 0, 6, 0, 7, 0, 3, -9, 2, -9, 0, -9, 1, -9, 352 7, -9, 6, -9, 4, -12, 5, -12, 2, -8, 0, -8, 3, -8, 1, -8, 4, 0, 5, 0, 6, 0, 7, 0, 0, -7, 3, -7, 2, -7, 1, -7, 4, 0, 5, 0, 6, 0, 7, 0, 0, -6, 1, -6, 3, -6, 2, -6, 4, -3, 5, -3, 7, -6, 6, -6, 353 1, -5, 3, -5, 0, -5, 2, -5, 4, 0, 5, 0, 6, 0, 7, 0, 3, -4, 0, -4, 1, -4, 2, -4, 4, 0, 5, 0, 6, 0, 7, 0, 1, -3, 0, -3, 2, -3, 3, -3, 5, -3, 4, -3, 6, -6, 7, -6, 0, -2, 2, -2, 1, -2, 3, -2, 354 4, 0, 5, 0, 6, 0, 7, 0, 2, -1, 1, -1, 0, -1, 3, -1, 4, 0, 5, 0, 6, 0, 7, 0, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 1, 1, 2, 1, 0, 1, 3, 1, 4, 0, 5, 0, 6, 0, 7, 0, 355 2, 2, 0, 2, 1, 2, 3, 2, 4, 0, 5, 0, 6, 0, 7, 0, 1, 3, 0, 3, 3, 3, 2, 3, 5, 0, 4, 0, 7, 0, 6, 0, 0, 4, 3, 4, 1, 4, 2, 4, 4, 0, 5, 0, 6, 0, 7, 0, 3, 5, 1, 5, 0, 5, 2, 5, 356 4, 0, 5, 0, 6, 0, 7, 0, 2, 6, 3, 6, 0, 6, 1, 6, 6, 6, 7, 6, 4, 6, 5, 6, 3, 7, 0, 7, 2, 7, 1, 7, 4, 0, 5, 0, 6, 0, 7, 0, 0, 8, 2, 8, 3, 8, 1, 8, 4, 0, 5, 0, 6, 0, 7, 0, 357 3, 9, 2, 9, 1, 9, 0, 9, 7, 6, 6, 6, 5, 6, 4, 6, 2, 10, 1, 10, 3, 10, 0, 10, 4, 0, 5, 0, 6, 0, 7, 0, 1, 11, 3, 11, 2, 11, 0, 11, 4, 0, 5, 0, 6, 0, 7, 0}; 358 static PetscInt hex_seg[] = {2, 0, 3, 0, 4, 0, 5, 0, 1, 0, 0, 0, 4, 0, 5, 0, 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 1, 0, 0, 0, 3, 0, 2, 0, 3, 0, 2, 0, 4, 0, 5, 0, 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 1, 0, 0, 0, 4, 0, 5, 0, 1, 0, 0, 0, 2, 0, 3, 0, 359 2, 0, 3, 0, 5, 0, 4, 0, 0, 0, 1, 0, 5, 0, 4, 0, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 3, 0, 2, 0, 0, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 4, 0, 5, 0, 360 1, 0, 0, 0, 4, 0, 5, 0, 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 2, 0, 3, 0, 5, 0, 4, 0, 2, 0, 3, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 3, 0, 4, 0, 5, 0, 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 0, 0, 1, 0, 5, 0, 4, 0, 361 1, 0, 0, 0, 3, 0, 2, 0, 5, 0, 4, 0, 2, 0, 3, 0, 1, 0, 0, 0, 5, 0, 4, 0, 0, 0, 1, 0, 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 4, 0, 5, 0, 0, 0, 1, 0, 5, 0, 4, 0, 3, 0, 2, 0, 0, 0, 1, 0, 2, 0, 3, 0, 5, 0, 4, 0, 362 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0, 5, 0, 4, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 0, 0, 1, 0, 4, 0, 5, 0, 3, 0, 2, 0, 2, 0, 3, 0, 1, 0, 0, 0, 4, 0, 5, 0, 1, 0, 0, 0, 3, 0, 2, 0, 4, 0, 5, 0, 363 3, 0, 2, 0, 0, 0, 1, 0, 4, 0, 5, 0, 4, 0, 5, 0, 2, 0, 3, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 3, 0, 5, 0, 4, 0, 5, 0, 4, 0, 2, 0, 3, 0, 0, 0, 1, 0, 1, 0, 0, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 2, 0, 3, 0, 364 2, 0, 3, 0, 0, 0, 1, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 3, 0, 2, 0, 0, 0, 1, 0, 5, 0, 4, 0, 0, 0, 1, 0, 3, 0, 2, 0, 2, 0, 3, 0, 5, 0, 4, 0, 1, 0, 0, 0, 365 4, 0, 5, 0, 1, 0, 0, 0, 3, 0, 2, 0, 3, 0, 2, 0, 5, 0, 4, 0, 0, 0, 1, 0, 3, 0, 2, 0, 4, 0, 5, 0, 1, 0, 0, 0, 5, 0, 4, 0, 1, 0, 0, 0, 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0, 2, 0, 3, 0, 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0}; 366 static PetscInt hex_quad[] = {7, -2, 4, -2, 5, -2, 6, -2, 8, -3, 11, -3, 10, -3, 9, -3, 3, 1, 2, 1, 1, 1, 0, 1, 8, -2, 9, -2, 10, -2, 11, -2, 3, -4, 0, -4, 1, -4, 2, -4, 7, 0, 4, 0, 5, 0, 6, 0, 9, 1, 8, 1, 11, 367 1, 10, 1, 0, 3, 3, 3, 2, 3, 1, 3, 5, 2, 6, 2, 7, 2, 4, 2, 6, 3, 5, 3, 4, 3, 7, 3, 10, -1, 9, -1, 8, -1, 11, -1, 2, -4, 3, -4, 0, -4, 1, -4, 4, 1, 7, 1, 6, 1, 5, 1, 11, 2, 368 8, 2, 9, 2, 10, 2, 1, 3, 0, 3, 3, 3, 2, 3, 10, -4, 11, -4, 8, -4, 9, -4, 2, 1, 1, 1, 0, 1, 3, 1, 6, -1, 5, -1, 4, -1, 7, -1, 5, -4, 6, -4, 7, -4, 4, -4, 9, 0, 10, 0, 11, 0, 8, 369 0, 0, -2, 1, -2, 2, -2, 3, -2, 11, 3, 10, 3, 9, 3, 8, 3, 1, -2, 2, -2, 3, -2, 0, -2, 4, -3, 7, -3, 6, -3, 5, -3, 11, -1, 8, -1, 9, -1, 10, -1, 7, 3, 4, 3, 5, 3, 6, 3, 2, 2, 1, 2, 370 0, 2, 3, 2, 10, 2, 9, 2, 8, 2, 11, 2, 5, 1, 6, 1, 7, 1, 4, 1, 1, -1, 2, -1, 3, -1, 0, -1, 5, 2, 4, 2, 7, 2, 6, 2, 1, 2, 0, 2, 3, 2, 2, 2, 10, -4, 9, -4, 8, -4, 11, -4, 4, 371 -3, 5, -3, 6, -3, 7, -3, 0, -3, 1, -3, 2, -3, 3, -3, 8, -2, 11, -2, 10, -2, 9, -2, 3, 1, 0, 1, 1, 1, 2, 1, 9, -4, 8, -4, 11, -4, 10, -4, 6, 3, 7, 3, 4, 3, 5, 3, 1, 3, 2, 3, 3, 3, 372 0, 3, 10, 1, 11, 1, 8, 1, 9, 1, 5, -4, 4, -4, 7, -4, 6, -4, 8, 0, 11, 0, 10, 0, 9, 0, 4, -4, 7, -4, 6, -4, 5, -4, 0, 0, 3, 0, 2, 0, 1, 0, 0, 0, 1, 0, 2, 0, 3, 0, 5, -1, 4, 373 -1, 7, -1, 6, -1, 9, -3, 8, -3, 11, -3, 10, -3, 9, -3, 10, -3, 11, -3, 8, -3, 6, -2, 5, -2, 4, -2, 7, -2, 3, -3, 0, -3, 1, -3, 2, -3, 7, 0, 6, 0, 5, 0, 4, 0, 2, -1, 3, -1, 0, -1, 1, -1, 374 11, 3, 8, 3, 9, 3, 10, 3, 2, 2, 3, 2, 0, 2, 1, 2, 6, 2, 7, 2, 4, 2, 5, 2, 10, 2, 11, 2, 8, 2, 9, 2, 6, -1, 7, -1, 4, -1, 5, -1, 3, 0, 2, 0, 1, 0, 0, 0, 9, 1, 10, 1, 11, 375 1, 8, 1, 2, -4, 1, -4, 0, -4, 3, -4, 11, -2, 10, -2, 9, -2, 8, -2, 7, -2, 6, -2, 5, -2, 4, -2, 1, -1, 0, -1, 3, -1, 2, -1, 4, 0, 5, 0, 6, 0, 7, 0, 11, -1, 10, -1, 9, -1, 8, -1, 0, -2, 376 3, -2, 2, -2, 1, -2, 8, 3, 9, 3, 10, 3, 11, 3, 4, 1, 5, 1, 6, 1, 7, 1, 3, -3, 2, -3, 1, -3, 0, -3, 7, -3, 6, -3, 5, -3, 4, -3, 8, 0, 9, 0, 10, 0, 11, 0, 0, 0, 1, 0, 2, 0, 3, 377 0, 4, 0, 5, 0, 6, 0, 7, 0, 8, 0, 9, 0, 10, 0, 11, 0, 1, 3, 2, 3, 3, 3, 0, 3, 11, -2, 10, -2, 9, -2, 8, -2, 4, 1, 5, 1, 6, 1, 7, 1, 2, 2, 3, 2, 0, 2, 1, 2, 7, -3, 6, -3, 378 5, -3, 4, -3, 11, -1, 10, -1, 9, -1, 8, -1, 3, 1, 0, 1, 1, 1, 2, 1, 8, 3, 9, 3, 10, 3, 11, 3, 7, -2, 6, -2, 5, -2, 4, -2, 5, 2, 4, 2, 7, 2, 6, 2, 0, -3, 1, -3, 2, -3, 3, -3, 9, 379 1, 10, 1, 11, 1, 8, 1, 1, -1, 0, -1, 3, -1, 2, -1, 5, -1, 4, -1, 7, -1, 6, -1, 10, 2, 11, 2, 8, 2, 9, 2, 4, -3, 5, -3, 6, -3, 7, -3, 1, 2, 0, 2, 3, 2, 2, 2, 11, 3, 8, 3, 9, 3, 380 10, 3, 8, 0, 11, 0, 10, 0, 9, 0, 7, 3, 4, 3, 5, 3, 6, 3, 3, -3, 0, -3, 1, -3, 2, -3, 3, -3, 2, -3, 1, -3, 0, -3, 6, 2, 7, 2, 4, 2, 5, 2, 9, -3, 8, -3, 11, -3, 10, -3, 9, -3, 10, 381 -3, 11, -3, 8, -3, 5, 1, 6, 1, 7, 1, 4, 1, 0, 0, 3, 0, 2, 0, 1, 0, 0, -2, 3, -2, 2, -2, 1, -2, 9, -4, 8, -4, 11, -4, 10, -4, 5, -4, 4, -4, 7, -4, 6, -4, 2, -4, 1, -4, 0, -4, 3, -4, 382 10, 1, 11, 1, 8, 1, 9, 1, 6, 3, 7, 3, 4, 3, 5, 3, 7, 0, 6, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 8, -2, 11, -2, 10, -2, 9, -2, 6, -1, 7, -1, 4, -1, 5, -1, 2, -1, 3, -1, 0, 383 -1, 1, -1, 10, -4, 9, -4, 8, -4, 11, -4, 11, -1, 8, -1, 9, -1, 10, -1, 4, -4, 7, -4, 6, -4, 5, -4, 1, -1, 2, -1, 3, -1, 0, -1, 10, 2, 9, 2, 8, 2, 11, 2, 6, -2, 5, -2, 4, -2, 7, -2, 2, 2, 384 1, 2, 0, 2, 3, 2, 8, -2, 9, -2, 10, -2, 11, -2, 0, 3, 3, 3, 2, 3, 1, 3, 4, -3, 7, -3, 6, -3, 5, -3, 4, 1, 7, 1, 6, 1, 5, 1, 8, -3, 11, -3, 10, -3, 9, -3, 0, -2, 1, -2, 2, -2, 3, 385 -2, 9, 1, 8, 1, 11, 1, 10, 1, 3, -4, 0, -4, 1, -4, 2, -4, 6, -1, 5, -1, 4, -1, 7, -1, 5, -4, 6, -4, 7, -4, 4, -4, 10, -1, 9, -1, 8, -1, 11, -1, 1, 3, 0, 3, 3, 3, 2, 3, 7, -2, 4, -2, 386 5, -2, 6, -2, 11, 2, 8, 2, 9, 2, 10, 2, 2, -4, 3, -4, 0, -4, 1, -4, 10, -4, 11, -4, 8, -4, 9, -4, 1, -2, 2, -2, 3, -2, 0, -2, 5, 2, 6, 2, 7, 2, 4, 2, 11, 3, 10, 3, 9, 3, 8, 3, 2, 387 1, 1, 1, 0, 1, 3, 1, 7, 0, 4, 0, 5, 0, 6, 0, 6, 3, 5, 3, 4, 3, 7, 3, 9, 0, 10, 0, 11, 0, 8, 0, 3, 1, 2, 1, 1, 1, 0, 1}; 388 static PetscInt hex_hex[] = {3, -24, 0, -24, 4, -24, 5, -24, 2, -24, 6, -24, 7, -24, 1, -24, 3, -23, 5, -23, 6, -23, 2, -23, 0, -23, 1, -23, 7, -23, 4, -23, 4, -22, 0, -22, 1, -22, 7, -22, 5, -22, 6, -22, 2, -22, 3, -22, 6, -21, 7, -21, 389 1, -21, 2, -21, 5, -21, 3, -21, 0, -21, 4, -21, 1, -20, 2, -20, 6, -20, 7, -20, 0, -20, 4, -20, 5, -20, 3, -20, 6, -19, 2, -19, 3, -19, 5, -19, 7, -19, 4, -19, 0, -19, 1, -19, 4, -18, 5, -18, 3, -18, 0, -18, 390 7, -18, 1, -18, 2, -18, 6, -18, 1, -17, 7, -17, 4, -17, 0, -17, 2, -17, 3, -17, 5, -17, 6, -17, 2, -16, 3, -16, 5, -16, 6, -16, 1, -16, 7, -16, 4, -16, 0, -16, 7, -15, 4, -15, 0, -15, 1, -15, 6, -15, 2, -15, 391 3, -15, 5, -15, 7, -14, 1, -14, 2, -14, 6, -14, 4, -14, 5, -14, 3, -14, 0, -14, 0, -13, 4, -13, 5, -13, 3, -13, 1, -13, 2, -13, 6, -13, 7, -13, 5, -12, 4, -12, 7, -12, 6, -12, 3, -12, 2, -12, 1, -12, 0, -12, 392 7, -11, 6, -11, 5, -11, 4, -11, 1, -11, 0, -11, 3, -11, 2, -11, 0, -10, 1, -10, 7, -10, 4, -10, 3, -10, 5, -10, 6, -10, 2, -10, 4, -9, 7, -9, 6, -9, 5, -9, 0, -9, 3, -9, 2, -9, 1, -9, 5, -8, 6, -8, 393 2, -8, 3, -8, 4, -8, 0, -8, 1, -8, 7, -8, 2, -7, 6, -7, 7, -7, 1, -7, 3, -7, 0, -7, 4, -7, 5, -7, 6, -6, 5, -6, 4, -6, 7, -6, 2, -6, 1, -6, 0, -6, 3, -6, 5, -5, 3, -5, 0, -5, 4, -5, 394 6, -5, 7, -5, 1, -5, 2, -5, 2, -4, 1, -4, 0, -4, 3, -4, 6, -4, 5, -4, 4, -4, 7, -4, 1, -3, 0, -3, 3, -3, 2, -3, 7, -3, 6, -3, 5, -3, 4, -3, 0, -2, 3, -2, 2, -2, 1, -2, 4, -2, 7, -2, 395 6, -2, 5, -2, 3, -1, 2, -1, 1, -1, 0, -1, 5, -1, 4, -1, 7, -1, 6, -1, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 1, 1, 2, 1, 3, 1, 0, 1, 7, 1, 4, 1, 5, 1, 6, 1, 396 2, 2, 3, 2, 0, 2, 1, 2, 6, 2, 7, 2, 4, 2, 5, 2, 3, 3, 0, 3, 1, 3, 2, 3, 5, 3, 6, 3, 7, 3, 4, 3, 4, 4, 0, 4, 3, 4, 5, 4, 7, 4, 6, 4, 2, 4, 1, 4, 7, 5, 4, 5, 397 5, 5, 6, 5, 1, 5, 2, 5, 3, 5, 0, 5, 1, 6, 7, 6, 6, 6, 2, 6, 0, 6, 3, 6, 5, 6, 4, 6, 3, 7, 2, 7, 6, 7, 5, 7, 0, 7, 4, 7, 7, 7, 1, 7, 5, 8, 6, 8, 7, 8, 4, 8, 398 3, 8, 0, 8, 1, 8, 2, 8, 4, 9, 7, 9, 1, 9, 0, 9, 5, 9, 3, 9, 2, 9, 6, 9, 4, 10, 5, 10, 6, 10, 7, 10, 0, 10, 1, 10, 2, 10, 3, 10, 6, 11, 7, 11, 4, 11, 5, 11, 2, 11, 3, 11, 399 0, 11, 1, 11, 3, 12, 5, 12, 4, 12, 0, 12, 2, 12, 1, 12, 7, 12, 6, 12, 6, 13, 2, 13, 1, 13, 7, 13, 5, 13, 4, 13, 0, 13, 3, 13, 1, 14, 0, 14, 4, 14, 7, 14, 2, 14, 6, 14, 5, 14, 3, 14, 400 6, 15, 5, 15, 3, 15, 2, 15, 7, 15, 1, 15, 0, 15, 4, 15, 0, 16, 4, 16, 7, 16, 1, 16, 3, 16, 2, 16, 6, 16, 5, 16, 0, 17, 3, 17, 5, 17, 4, 17, 1, 17, 7, 17, 6, 17, 2, 17, 5, 18, 3, 18, 401 2, 18, 6, 18, 4, 18, 7, 18, 1, 18, 0, 18, 7, 19, 6, 19, 2, 19, 1, 19, 4, 19, 0, 19, 3, 19, 5, 19, 2, 20, 1, 20, 7, 20, 6, 20, 3, 20, 5, 20, 4, 20, 0, 20, 7, 21, 1, 21, 0, 21, 4, 21, 402 6, 21, 5, 21, 3, 21, 2, 21, 2, 22, 6, 22, 5, 22, 3, 22, 1, 22, 0, 22, 4, 22, 7, 22, 5, 23, 4, 23, 0, 23, 3, 23, 6, 23, 2, 23, 1, 23, 7, 23}; 403 static PetscInt trip_seg[] = {1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 1, 0, 2, 0, 0, -1, 2, -1, 1, -1, 1, -1, 0, -1, 2, -1, 2, -1, 1, -1, 0, -1, 404 0, 0, 1, 0, 2, 0, 2, 0, 0, 0, 1, 0, 1, 0, 2, 0, 0, 0, 2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, -1, 2, -1, 1, -1}; 405 static PetscInt trip_tri[] = {1, 1, 2, 1, 0, 1, 3, 1, 2, 2, 0, 2, 1, 2, 3, 2, 0, 0, 1, 0, 2, 0, 3, 0, 2, -1, 1, -1, 0, -1, 3, -3, 0, -2, 2, -2, 1, -2, 3, -1, 1, -3, 0, -3, 2, -3, 3, -2, 406 0, 0, 1, 0, 2, 0, 3, 0, 2, 2, 0, 2, 1, 2, 3, 2, 1, 1, 2, 1, 0, 1, 3, 1, 1, -3, 0, -3, 2, -3, 3, -2, 0, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 3, -3}; 407 static PetscInt trip_quad[] = {4, -1, 5, -1, 3, -1, 1, -1, 2, -1, 0, -1, 5, -1, 3, -1, 4, -1, 2, -1, 0, -1, 1, -1, 3, -1, 4, -1, 5, -1, 0, -1, 1, -1, 2, -1, 0, -3, 2, -3, 1, -3, 3, -3, 5, -3, 4, -3, 408 1, -3, 0, -3, 2, -3, 4, -3, 3, -3, 5, -3, 2, -3, 1, -3, 0, -3, 5, -3, 4, -3, 3, -3, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 2, 0, 0, 0, 1, 0, 5, 0, 3, 0, 4, 0, 409 1, 0, 2, 0, 0, 0, 4, 0, 5, 0, 3, 0, 5, 2, 4, 2, 3, 2, 2, 2, 1, 2, 0, 2, 4, 2, 3, 2, 5, 2, 1, 2, 0, 2, 2, 2, 3, 2, 5, 2, 4, 2, 0, 2, 2, 2, 1, 2}; 410 static PetscInt trip_trip[] = {5, -6, 6, -6, 4, -6, 7, -6, 1, -6, 2, -6, 0, -6, 3, -6, 6, -5, 4, -5, 5, -5, 7, -5, 2, -5, 0, -5, 1, -5, 3, -5, 4, -4, 5, -4, 6, -4, 7, -4, 0, -4, 1, -4, 2, -4, 3, -4, 411 2, -3, 1, -3, 0, -3, 3, -1, 6, -3, 5, -3, 4, -3, 7, -1, 0, -2, 2, -2, 1, -2, 3, -3, 4, -2, 6, -2, 5, -2, 7, -3, 1, -1, 0, -1, 2, -1, 3, -2, 5, -1, 4, -1, 6, -1, 7, -2, 412 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 2, 1, 0, 1, 1, 1, 3, 1, 6, 1, 4, 1, 5, 1, 7, 1, 1, 2, 2, 2, 0, 2, 3, 2, 5, 2, 6, 2, 4, 2, 7, 2, 413 5, 3, 4, 3, 6, 3, 7, 4, 1, 3, 0, 3, 2, 3, 3, 4, 4, 4, 6, 4, 5, 4, 7, 5, 0, 4, 2, 4, 1, 4, 3, 5, 6, 5, 5, 5, 4, 5, 7, 3, 2, 5, 1, 5, 0, 5, 3, 3}; 414 static PetscInt ttri_tseg[] = {2, -2, 1, -2, 0, -2, 1, -2, 0, -2, 2, -2, 0, -2, 2, -2, 1, -2, 2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, -1, 2, -1, 1, -1, 415 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}; 416 static PetscInt ttri_ttri[] = {1, -6, 0, -6, 2, -6, 3, -5, 0, -5, 2, -5, 1, -5, 3, -4, 2, -4, 1, -4, 0, -4, 3, -6, 1, -3, 0, -3, 2, -3, 3, -2, 0, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 3, -3, 417 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 0, 1, 3, 1, 2, 2, 0, 2, 1, 2, 3, 2, 0, 3, 1, 3, 2, 3, 3, 3, 1, 4, 2, 4, 0, 4, 3, 4, 2, 5, 0, 5, 1, 5, 3, 5}; 418 static PetscInt tquad_tvert[] = {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}; 419 static PetscInt tquad_tseg[] = {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, 0, 2, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 2, 0, 1, 0, 0, 0, 3, 0, 420 0, 0, 1, 0, 2, 0, 3, 0, 1, 0, 2, 0, 3, 0, 0, 0, 2, 0, 3, 0, 0, 0, 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}; 421 static PetscInt tquad_tquad[] = {2, -8, 1, -8, 0, -8, 3, -8, 1, -7, 0, -7, 3, -7, 2, -7, 0, -6, 3, -6, 2, -6, 1, -6, 3, -5, 2, -5, 1, -5, 0, -5, 2, -4, 1, -4, 0, -4, 3, -4, 1, -3, 0, 422 -3, 3, -3, 2, -3, 0, -2, 3, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 3, 1, 0, 1, 2, 2, 3, 2, 0, 2, 423 1, 2, 3, 3, 0, 3, 1, 3, 2, 3, 0, 4, 1, 4, 2, 4, 3, 4, 1, 5, 2, 5, 3, 5, 0, 5, 2, 6, 3, 6, 0, 6, 1, 6, 3, 7, 0, 7, 1, 7, 2, 7}; 424 425 PetscFunctionBeginHot; 426 *rnew = r; 427 *onew = o; 428 if (!so) PetscFunctionReturn(0); 429 switch (sct) { 430 case DM_POLYTOPE_POINT: break; 431 case DM_POLYTOPE_POINT_PRISM_TENSOR: *onew = so < 0 ? -(o + 1) : o; break; 432 case DM_POLYTOPE_SEGMENT: 433 switch (tct) { 434 case DM_POLYTOPE_POINT: break; 435 case DM_POLYTOPE_SEGMENT: 436 *rnew = seg_seg[(so + 1) * 4 + r * 2]; 437 *onew = DMPolytopeTypeComposeOrientation(tct, o, seg_seg[(so + 1) * 4 + r * 2 + 1]); 438 break; 439 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 440 } 441 break; 442 case DM_POLYTOPE_TRIANGLE: 443 switch (tct) { 444 case DM_POLYTOPE_SEGMENT: 445 *rnew = tri_seg[(so + 3) * 6 + r * 2]; 446 *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so + 3) * 6 + r * 2 + 1]); 447 break; 448 case DM_POLYTOPE_TRIANGLE: 449 *rnew = tri_tri[(so + 3) * 8 + r * 2]; 450 *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so + 3) * 8 + r * 2 + 1]); 451 break; 452 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 453 } 454 break; 455 case DM_POLYTOPE_QUADRILATERAL: 456 switch (tct) { 457 case DM_POLYTOPE_POINT: break; 458 case DM_POLYTOPE_SEGMENT: 459 *rnew = quad_seg[(so + 4) * 8 + r * 2]; 460 *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_seg[(so + 4) * 8 + r * 2 + 1]); 461 break; 462 case DM_POLYTOPE_QUADRILATERAL: 463 *rnew = quad_quad[(so + 4) * 8 + r * 2]; 464 *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_quad[(so + 4) * 8 + r * 2 + 1]); 465 break; 466 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 467 } 468 break; 469 case DM_POLYTOPE_SEG_PRISM_TENSOR: 470 switch (tct) { 471 case DM_POLYTOPE_POINT_PRISM_TENSOR: 472 *rnew = tseg_seg[(so + 2) * 2 + r * 2]; 473 *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_seg[(so + 2) * 2 + r * 2 + 1]); 474 break; 475 case DM_POLYTOPE_SEG_PRISM_TENSOR: 476 *rnew = tseg_tseg[(so + 2) * 4 + r * 2]; 477 *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_tseg[(so + 2) * 4 + r * 2 + 1]); 478 break; 479 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 480 } 481 break; 482 case DM_POLYTOPE_TETRAHEDRON: 483 switch (tct) { 484 case DM_POLYTOPE_SEGMENT: 485 *rnew = tet_seg[(so + 12) * 2 + r * 2]; 486 *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so + 12) * 2 + r * 2 + 1]); 487 break; 488 case DM_POLYTOPE_TRIANGLE: 489 *rnew = tet_tri[(so + 12) * 16 + r * 2]; 490 *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so + 12) * 16 + r * 2 + 1]); 491 break; 492 case DM_POLYTOPE_TETRAHEDRON: 493 *rnew = tet_tet[(so + 12) * 16 + r * 2]; 494 *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so + 12) * 16 + r * 2 + 1]); 495 break; 496 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 497 } 498 break; 499 case DM_POLYTOPE_HEXAHEDRON: 500 switch (tct) { 501 case DM_POLYTOPE_POINT: break; 502 case DM_POLYTOPE_SEGMENT: 503 *rnew = hex_seg[(so + 24) * 12 + r * 2]; 504 *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_seg[(so + 24) * 12 + r * 2 + 1]); 505 break; 506 case DM_POLYTOPE_QUADRILATERAL: 507 *rnew = hex_quad[(so + 24) * 24 + r * 2]; 508 *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_quad[(so + 24) * 24 + r * 2 + 1]); 509 break; 510 case DM_POLYTOPE_HEXAHEDRON: 511 *rnew = hex_hex[(so + 24) * 16 + r * 2]; 512 *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_hex[(so + 24) * 16 + r * 2 + 1]); 513 break; 514 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 515 } 516 break; 517 case DM_POLYTOPE_TRI_PRISM: 518 switch (tct) { 519 case DM_POLYTOPE_SEGMENT: 520 *rnew = trip_seg[(so + 6) * 6 + r * 2]; 521 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_seg[(so + 6) * 6 + r * 2 + 1]); 522 break; 523 case DM_POLYTOPE_TRIANGLE: 524 *rnew = trip_tri[(so + 6) * 8 + r * 2]; 525 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_tri[(so + 6) * 8 + r * 2 + 1]); 526 break; 527 case DM_POLYTOPE_QUADRILATERAL: 528 *rnew = trip_quad[(so + 6) * 12 + r * 2]; 529 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_quad[(so + 6) * 12 + r * 2 + 1]); 530 break; 531 case DM_POLYTOPE_TRI_PRISM: 532 *rnew = trip_trip[(so + 6) * 16 + r * 2]; 533 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_trip[(so + 6) * 16 + r * 2 + 1]); 534 break; 535 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 536 } 537 break; 538 case DM_POLYTOPE_TRI_PRISM_TENSOR: 539 switch (tct) { 540 case DM_POLYTOPE_SEG_PRISM_TENSOR: 541 *rnew = ttri_tseg[(so + 6) * 6 + r * 2]; 542 *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_tseg[(so + 6) * 6 + r * 2 + 1]); 543 break; 544 case DM_POLYTOPE_TRI_PRISM_TENSOR: 545 *rnew = ttri_ttri[(so + 6) * 8 + r * 2]; 546 *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_ttri[(so + 6) * 8 + r * 2 + 1]); 547 break; 548 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 549 } 550 break; 551 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 552 switch (tct) { 553 case DM_POLYTOPE_POINT_PRISM_TENSOR: 554 *rnew = tquad_tvert[(so + 8) * 2 + r * 2]; 555 *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tvert[(so + 8) * 2 + r * 2 + 1]); 556 break; 557 case DM_POLYTOPE_SEG_PRISM_TENSOR: 558 *rnew = tquad_tseg[(so + 8) * 8 + r * 2]; 559 *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tseg[(so + 8) * 8 + r * 2 + 1]); 560 break; 561 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 562 *rnew = tquad_tquad[(so + 8) * 8 + r * 2]; 563 *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tquad[(so + 8) * 8 + r * 2 + 1]); 564 break; 565 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 566 } 567 break; 568 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]); 569 } 570 PetscFunctionReturn(0); 571 } 572 573 PetscErrorCode DMPlexTransformCellRefine_Regular(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) { 574 /* All vertices remain in the refined mesh */ 575 static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; 576 static PetscInt vertexS[] = {1}; 577 static PetscInt vertexC[] = {0}; 578 static PetscInt vertexO[] = {0}; 579 /* Split all edges with a new vertex, making two new 2 edges 580 0--0--0--1--1 581 */ 582 static DMPolytopeType segT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}; 583 static PetscInt segS[] = {1, 2}; 584 static PetscInt segC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 585 static PetscInt segO[] = {0, 0, 0, 0}; 586 /* Do not split tensor edges */ 587 static DMPolytopeType tvertT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; 588 static PetscInt tvertS[] = {1}; 589 static PetscInt tvertC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 590 static PetscInt tvertO[] = {0, 0}; 591 /* Add 3 edges inside every triangle, making 4 new triangles. 592 2 593 |\ 594 | \ 595 | \ 596 0 1 597 | C \ 598 | \ 599 | \ 600 2---1---1 601 |\ D / \ 602 1 2 0 0 603 |A \ / B \ 604 0-0-0---1---1 605 */ 606 static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE}; 607 static PetscInt triS[] = {3, 4}; 608 static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 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, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2}; 609 static PetscInt triO[] = {0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0}; 610 /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads. 611 3----1----2----0----2 612 | | | 613 0 D 2 C 1 614 | | | 615 3----3----0----1----1 616 | | | 617 1 A 0 B 0 618 | | | 619 0----0----0----1----1 620 */ 621 static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL}; 622 static PetscInt quadS[] = {1, 4, 4}; 623 static PetscInt quadC[] = {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, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 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, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0}; 624 static PetscInt quadO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, -1, 0, 0}; 625 /* Add 1 edge inside every tensor quad, making 2 new tensor quads 626 2----2----1----3----3 627 | | | 628 | | | 629 | | | 630 4 A 6 B 5 631 | | | 632 | | | 633 | | | 634 0----0----0----1----1 635 */ 636 static DMPolytopeType tsegT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR}; 637 static PetscInt tsegS[] = {1, 2}; 638 static PetscInt tsegC[] = {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, 0, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0}; 639 static PetscInt tsegO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 640 /* Add 1 edge and 8 triangles inside every cell, making 8 new tets 641 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 642 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] 643 called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are 644 [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3] 645 The first four tets just cut off the corners, using the replica notation for new vertices, 646 [v0, (e0, 0), (e2, 0), (e3, 0)] 647 [(e0, 0), v1, (e1, 0), (e4, 0)] 648 [(e2, 0), (e1, 0), v2, (e5, 0)] 649 [(e3, 0), (e4, 0), (e5, 0), v3 ] 650 The next four tets match a vertex to the newly created faces from cutting off those first tets. 651 [(e2, 0), (e3, 0), (e0, 0), (e5, 0)] 652 [(e4, 0), (e1, 0), (e0, 0), (e5, 0)] 653 [(e5, 0), (e0, 0), (e2, 0), (e1, 0)] 654 [(e5, 0), (e0, 0), (e4, 0), (e3, 0)] 655 We can see that a new edge is introduced in the cell [(e0, 0), (e5, 0)] which we call (-1, 0). The first four faces created are 656 [(e2, 0), (e0, 0), (e3, 0)] 657 [(e0, 0), (e1, 0), (e4, 0)] 658 [(e2, 0), (e5, 0), (e1, 0)] 659 [(e3, 0), (e4, 0), (e5, 0)] 660 The next four, from the second group of tets, are 661 [(e2, 0), (e0, 0), (e5, 0)] 662 [(e4, 0), (e0, 0), (e5, 0)] 663 [(e0, 0), (e1, 0), (e5, 0)] 664 [(e5, 0), (e3, 0), (e0, 0)] 665 I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table. 666 */ 667 static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON}; 668 static PetscInt tetS[] = {1, 8, 8}; 669 static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1, DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3}; 670 static PetscInt tetO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, -1, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -1, -1, 1, 0, 0, -1, -1, -3, 2, -1, 0, -1, 1}; 671 /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes 672 The vertices of our reference hex are (-1, -1, -1), (-1, 1, -1), (1, 1, -1), (1, -1, -1), (-1, -1, 1), (1, -1, 1), (1, 1, 1), (-1, 1, 1) which we call [v0, v1, v2, v3, v4, v5, v6, v7]. The fours edges around the bottom [v0, v1], [v1, v2], [v2, v3], [v3, v0] are [e0, e1, e2, e3], and likewise around the top [v4, v5], [v5, v6], [v6, v7], [v7, v4] are [e4, e5, e6, e7]. Finally [v0, v4], [v1, v7], [v2, v6], [v3, v5] are [e9, e10, e11, e8]. The faces of a hex, given in DMPlexGetRawFaces_Internal(), oriented with outward normals, are 673 [v0, v1, v2, v3] f0 bottom 674 [v4, v5, v6, v7] f1 top 675 [v0, v3, v5, v4] f2 front 676 [v2, v1, v7, v6] f3 back 677 [v3, v2, v6, v5] f4 right 678 [v0, v4, v7, v1] f5 left 679 The eight hexes are divided into four on the bottom, and four on the top, 680 [v0, (e0, 0), (f0, 0), (e3, 0), (e9, 0), (f2, 0), (c0, 0), (f5, 0)] 681 [(e0, 0), v1, (e1, 0), (f0, 0), (f5, 0), (c0, 0), (f3, 0), (e10, 0)] 682 [(f0, 0), (e1, 0), v2, (e2, 0), (c0, 0), (f4, 0), (e11, 0), (f3, 0)] 683 [(e3, 0), (f0, 0), (e2, 0), v3, (f2, 0), (e8, 0), (f4, 0), (c0, 0)] 684 [(e9, 0), (f5, 0), (c0, 0), (f2, 0), v4, (e4, 0), (f1, 0), (e7, 0)] 685 [(f2, 0), (c0, 0), (f4, 0), (e8, 0), (e4, 0), v5, (e5, 0), (f1, 0)] 686 [(c0, 0), (f3, 0), (e11, 0), (f4, 0), (f1, 0), (e5, 0), v6, (e6, 0)] 687 [(f5, 0), (e10, 0), (f3, 0), (c0, 0), (e7, 0), (f1, 0), (e6, 0), v7] 688 The 6 internal edges will go from the faces to the central vertex. The 12 internal faces can be divided into groups of 4 by the plane on which they sit. First the faces on the x-y plane are, 689 [(e9, 0), (f2, 0), (c0, 0), (f5, 0)] 690 [(f5, 0), (c0, 0), (f3, 0), (e10, 0)] 691 [(c0, 0), (f4, 0), (e11, 0), (f3, 0)] 692 [(f2, 0), (e8, 0), (f4, 0), (c0, 0)] 693 and on the x-z plane, 694 [(f0, 0), (e0, 0), (f5, 0), (c0, 0)] 695 [(c0, 0), (f5, 0), (e7, 0), (f1, 0)] 696 [(f4, 0), (c0, 0), (f1, 0), (e5, 0)] 697 [(e2, 0), (f0, 0), (c0, 0), (f4, 0)] 698 and on the y-z plane, 699 [(e3, 0), (f2, 0), (c0, 0), (f0, 0)] 700 [(f2, 0), (e4, 0), (f1, 0), (c0, 0)] 701 [(c0, 0), (f1, 0), (e6, 0), (f3, 0)] 702 [(f0, 0), (c0, 0), (f3, 0), (e1, 0)] 703 */ 704 static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 705 static PetscInt hexS[] = {1, 6, 12, 8}; 706 static PetscInt hexC[] = {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_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 11, DM_POLYTOPE_QUADRILATERAL, 1, 5, 3, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, DM_POLYTOPE_QUADRILATERAL, 0, 11, DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 0, 9, DM_POLYTOPE_QUADRILATERAL, 1, 5, 1, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3, DM_POLYTOPE_QUADRILATERAL, 0, 9, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2, DM_POLYTOPE_QUADRILATERAL, 0, 10, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 3, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 0, 10, DM_POLYTOPE_QUADRILATERAL, 1, 5, 2}; 707 static PetscInt hexO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 0, -1, -1, 0, -1, -1, 0, 0, -1, 0, 0, -1, -1, 0, 0, -1, -1, -1, 0, 0, 0, -1, -1, 0, 0, 0, -1, -1, 0, 0, -1, -1, -1, 0, 0, -1, -1, -1, 708 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, -2, 0, 0, 0, -3, 0, -2, 0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, -2, 0, -2, 0, 0, 0, 0, 0, -2, 0, -3, 0, 0, 0, -2, 0, -3, 0, -2, 0}; 709 /* Add 3 quads inside every triangular prism, making 4 new prisms. */ 710 static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM}; 711 static PetscInt tripS[] = {3, 4, 6, 8}; 712 static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 0, 5}; 713 static PetscInt tripO[] = {0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, -1, -1, -1, 0, 0, 0, 0, -1, 0, -1, -1, -1, 0, -1, -1, -1, 0, -1, -1, 0, -1, -1, 0, 0, -1, -1, 0, 0, -1, -1, 714 0, 0, 0, 0, -3, 0, 0, 0, 0, 0, -3, 0, 0, -3, 0, 0, 2, 0, 0, 0, 0, -2, 0, 0, -3, 0, -2, 0, 0, 0, -3, -2, 0, -3, 0, 0, -2, 0, 0, 0, 0}; 715 /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms. 716 2 717 |\ 718 | \ 719 | \ 720 0---1 721 722 2 723 724 0 1 725 726 2 727 |\ 728 | \ 729 | \ 730 0---1 731 */ 732 static DMPolytopeType ttriT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR}; 733 static PetscInt ttriS[] = {3, 4}; 734 static PetscInt ttriC[] = {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, 1, 3, 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, 1, 4, 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, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 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, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2}; 735 static PetscInt ttriO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0}; 736 /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */ 737 static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR}; 738 static PetscInt tquadS[] = {1, 4, 4}; 739 static PetscInt tquadC[] = {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_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 5, 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, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 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, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 1, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0}; 740 static PetscInt tquadO[] = {0, 0, 0, 0, 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, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0}; 741 /* Add 4 edges, 12 triangles, 1 quad, 4 tetrahedra, and 6 pyramids inside every pyramid. */ 742 static DMPolytopeType tpyrT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TETRAHEDRON, DM_POLYTOPE_PYRAMID}; 743 static PetscInt tpyrS[] = {4, 12, 1, 4, 6}; 744 static PetscInt tpyrC[] = { 745 DM_POLYTOPE_POINT, 746 2, 747 1, 748 1, 749 0, 750 DM_POLYTOPE_POINT, 751 1, 752 0, 753 0, 754 DM_POLYTOPE_POINT, 755 2, 756 2, 757 1, 758 0, 759 DM_POLYTOPE_POINT, 760 1, 761 0, 762 0, 763 DM_POLYTOPE_POINT, 764 2, 765 3, 766 1, 767 0, 768 DM_POLYTOPE_POINT, 769 1, 770 0, 771 0, 772 DM_POLYTOPE_POINT, 773 2, 774 4, 775 1, 776 0, 777 DM_POLYTOPE_POINT, 778 1, 779 0, 780 0, 781 /* These four triangle face out of the bottom pyramid */ 782 DM_POLYTOPE_SEGMENT, 783 1, 784 1, 785 1, 786 DM_POLYTOPE_SEGMENT, 787 0, 788 3, 789 DM_POLYTOPE_SEGMENT, 790 0, 791 0, 792 DM_POLYTOPE_SEGMENT, 793 1, 794 2, 795 1, 796 DM_POLYTOPE_SEGMENT, 797 0, 798 0, 799 DM_POLYTOPE_SEGMENT, 800 0, 801 1, 802 DM_POLYTOPE_SEGMENT, 803 1, 804 3, 805 1, 806 DM_POLYTOPE_SEGMENT, 807 0, 808 1, 809 DM_POLYTOPE_SEGMENT, 810 0, 811 2, 812 DM_POLYTOPE_SEGMENT, 813 1, 814 4, 815 1, 816 DM_POLYTOPE_SEGMENT, 817 0, 818 2, 819 DM_POLYTOPE_SEGMENT, 820 0, 821 3, 822 /* These eight triangles face out of the corner pyramids */ 823 DM_POLYTOPE_SEGMENT, 824 1, 825 0, 826 3, 827 DM_POLYTOPE_SEGMENT, 828 0, 829 3, 830 DM_POLYTOPE_SEGMENT, 831 1, 832 1, 833 2, 834 DM_POLYTOPE_SEGMENT, 835 1, 836 0, 837 2, 838 DM_POLYTOPE_SEGMENT, 839 0, 840 0, 841 DM_POLYTOPE_SEGMENT, 842 1, 843 2, 844 2, 845 DM_POLYTOPE_SEGMENT, 846 1, 847 0, 848 1, 849 DM_POLYTOPE_SEGMENT, 850 0, 851 1, 852 DM_POLYTOPE_SEGMENT, 853 1, 854 3, 855 2, 856 DM_POLYTOPE_SEGMENT, 857 1, 858 0, 859 0, 860 DM_POLYTOPE_SEGMENT, 861 0, 862 2, 863 DM_POLYTOPE_SEGMENT, 864 1, 865 4, 866 2, 867 DM_POLYTOPE_SEGMENT, 868 1, 869 0, 870 3, 871 DM_POLYTOPE_SEGMENT, 872 1, 873 1, 874 0, 875 DM_POLYTOPE_SEGMENT, 876 0, 877 0, 878 DM_POLYTOPE_SEGMENT, 879 1, 880 0, 881 2, 882 DM_POLYTOPE_SEGMENT, 883 1, 884 2, 885 0, 886 DM_POLYTOPE_SEGMENT, 887 0, 888 1, 889 DM_POLYTOPE_SEGMENT, 890 1, 891 0, 892 1, 893 DM_POLYTOPE_SEGMENT, 894 1, 895 3, 896 0, 897 DM_POLYTOPE_SEGMENT, 898 0, 899 2, 900 DM_POLYTOPE_SEGMENT, 901 1, 902 0, 903 0, 904 DM_POLYTOPE_SEGMENT, 905 1, 906 4, 907 0, 908 DM_POLYTOPE_SEGMENT, 909 0, 910 3, 911 /* This quad faces out of the bottom pyramid */ 912 DM_POLYTOPE_SEGMENT, 913 1, 914 1, 915 1, 916 DM_POLYTOPE_SEGMENT, 917 1, 918 2, 919 1, 920 DM_POLYTOPE_SEGMENT, 921 1, 922 3, 923 1, 924 DM_POLYTOPE_SEGMENT, 925 1, 926 4, 927 1, 928 /* The bottom face of each tet is on the triangular face */ 929 DM_POLYTOPE_TRIANGLE, 930 1, 931 1, 932 3, 933 DM_POLYTOPE_TRIANGLE, 934 0, 935 8, 936 DM_POLYTOPE_TRIANGLE, 937 0, 938 4, 939 DM_POLYTOPE_TRIANGLE, 940 0, 941 0, 942 DM_POLYTOPE_TRIANGLE, 943 1, 944 2, 945 3, 946 DM_POLYTOPE_TRIANGLE, 947 0, 948 9, 949 DM_POLYTOPE_TRIANGLE, 950 0, 951 5, 952 DM_POLYTOPE_TRIANGLE, 953 0, 954 1, 955 DM_POLYTOPE_TRIANGLE, 956 1, 957 3, 958 3, 959 DM_POLYTOPE_TRIANGLE, 960 0, 961 10, 962 DM_POLYTOPE_TRIANGLE, 963 0, 964 6, 965 DM_POLYTOPE_TRIANGLE, 966 0, 967 2, 968 DM_POLYTOPE_TRIANGLE, 969 1, 970 4, 971 3, 972 DM_POLYTOPE_TRIANGLE, 973 0, 974 11, 975 DM_POLYTOPE_TRIANGLE, 976 0, 977 7, 978 DM_POLYTOPE_TRIANGLE, 979 0, 980 3, 981 /* The front face of all pyramids is toward the front */ 982 DM_POLYTOPE_QUADRILATERAL, 983 1, 984 0, 985 0, 986 DM_POLYTOPE_TRIANGLE, 987 1, 988 1, 989 0, 990 DM_POLYTOPE_TRIANGLE, 991 0, 992 4, 993 DM_POLYTOPE_TRIANGLE, 994 0, 995 11, 996 DM_POLYTOPE_TRIANGLE, 997 1, 998 4, 999 1, 1000 DM_POLYTOPE_QUADRILATERAL, 1001 1, 1002 0, 1003 3, 1004 DM_POLYTOPE_TRIANGLE, 1005 1, 1006 1, 1007 1, 1008 DM_POLYTOPE_TRIANGLE, 1009 1, 1010 2, 1011 0, 1012 DM_POLYTOPE_TRIANGLE, 1013 0, 1014 5, 1015 DM_POLYTOPE_TRIANGLE, 1016 0, 1017 8, 1018 DM_POLYTOPE_QUADRILATERAL, 1019 1, 1020 0, 1021 2, 1022 DM_POLYTOPE_TRIANGLE, 1023 0, 1024 9, 1025 DM_POLYTOPE_TRIANGLE, 1026 1, 1027 2, 1028 1, 1029 DM_POLYTOPE_TRIANGLE, 1030 1, 1031 3, 1032 0, 1033 DM_POLYTOPE_TRIANGLE, 1034 0, 1035 6, 1036 DM_POLYTOPE_QUADRILATERAL, 1037 1, 1038 0, 1039 1, 1040 DM_POLYTOPE_TRIANGLE, 1041 0, 1042 7, 1043 DM_POLYTOPE_TRIANGLE, 1044 0, 1045 10, 1046 DM_POLYTOPE_TRIANGLE, 1047 1, 1048 3, 1049 1, 1050 DM_POLYTOPE_TRIANGLE, 1051 1, 1052 4, 1053 0, 1054 DM_POLYTOPE_QUADRILATERAL, 1055 0, 1056 0, 1057 DM_POLYTOPE_TRIANGLE, 1058 1, 1059 1, 1060 2, 1061 DM_POLYTOPE_TRIANGLE, 1062 1, 1063 2, 1064 2, 1065 DM_POLYTOPE_TRIANGLE, 1066 1, 1067 3, 1068 2, 1069 DM_POLYTOPE_TRIANGLE, 1070 1, 1071 4, 1072 2, 1073 DM_POLYTOPE_QUADRILATERAL, 1074 0, 1075 0, 1076 DM_POLYTOPE_TRIANGLE, 1077 0, 1078 0, 1079 DM_POLYTOPE_TRIANGLE, 1080 0, 1081 3, 1082 DM_POLYTOPE_TRIANGLE, 1083 0, 1084 2, 1085 DM_POLYTOPE_TRIANGLE, 1086 0, 1087 1, 1088 }; 1089 static PetscInt tpyrO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, -1, -1, 1090 -1, 0, -3, -2, -3, 0, -3, -2, -3, 0, -3, -2, -3, 0, -3, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 0, 0}; 1091 1092 PetscFunctionBegin; 1093 if (rt) *rt = 0; 1094 switch (source) { 1095 case DM_POLYTOPE_POINT: 1096 *Nt = 1; 1097 *target = vertexT; 1098 *size = vertexS; 1099 *cone = vertexC; 1100 *ornt = vertexO; 1101 break; 1102 case DM_POLYTOPE_SEGMENT: 1103 *Nt = 2; 1104 *target = segT; 1105 *size = segS; 1106 *cone = segC; 1107 *ornt = segO; 1108 break; 1109 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1110 *Nt = 1; 1111 *target = tvertT; 1112 *size = tvertS; 1113 *cone = tvertC; 1114 *ornt = tvertO; 1115 break; 1116 case DM_POLYTOPE_TRIANGLE: 1117 *Nt = 2; 1118 *target = triT; 1119 *size = triS; 1120 *cone = triC; 1121 *ornt = triO; 1122 break; 1123 case DM_POLYTOPE_QUADRILATERAL: 1124 *Nt = 3; 1125 *target = quadT; 1126 *size = quadS; 1127 *cone = quadC; 1128 *ornt = quadO; 1129 break; 1130 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1131 *Nt = 2; 1132 *target = tsegT; 1133 *size = tsegS; 1134 *cone = tsegC; 1135 *ornt = tsegO; 1136 break; 1137 case DM_POLYTOPE_TETRAHEDRON: 1138 *Nt = 3; 1139 *target = tetT; 1140 *size = tetS; 1141 *cone = tetC; 1142 *ornt = tetO; 1143 break; 1144 case DM_POLYTOPE_HEXAHEDRON: 1145 *Nt = 4; 1146 *target = hexT; 1147 *size = hexS; 1148 *cone = hexC; 1149 *ornt = hexO; 1150 break; 1151 case DM_POLYTOPE_TRI_PRISM: 1152 *Nt = 4; 1153 *target = tripT; 1154 *size = tripS; 1155 *cone = tripC; 1156 *ornt = tripO; 1157 break; 1158 case DM_POLYTOPE_TRI_PRISM_TENSOR: 1159 *Nt = 2; 1160 *target = ttriT; 1161 *size = ttriS; 1162 *cone = ttriC; 1163 *ornt = ttriO; 1164 break; 1165 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1166 *Nt = 3; 1167 *target = tquadT; 1168 *size = tquadS; 1169 *cone = tquadC; 1170 *ornt = tquadO; 1171 break; 1172 case DM_POLYTOPE_PYRAMID: 1173 *Nt = 5; 1174 *target = tpyrT; 1175 *size = tpyrS; 1176 *cone = tpyrC; 1177 *ornt = tpyrO; 1178 break; 1179 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1180 } 1181 PetscFunctionReturn(0); 1182 } 1183 1184 static PetscErrorCode DMPlexTransformInitialize_Regular(DMPlexTransform tr) { 1185 PetscFunctionBegin; 1186 tr->ops->view = DMPlexTransformView_Regular; 1187 tr->ops->setup = DMPlexTransformSetUp_Regular; 1188 tr->ops->destroy = DMPlexTransformDestroy_Regular; 1189 tr->ops->celltransform = DMPlexTransformCellRefine_Regular; 1190 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Regular; 1191 tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal; 1192 PetscFunctionReturn(0); 1193 } 1194 1195 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform tr) { 1196 DMPlexRefine_Regular *f; 1197 1198 PetscFunctionBegin; 1199 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1200 PetscCall(PetscNew(&f)); 1201 tr->data = f; 1202 1203 PetscCall(DMPlexTransformInitialize_Regular(tr)); 1204 PetscFunctionReturn(0); 1205 } 1206