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