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