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