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: SETERRQ1(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: SETERRQ1(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: SETERRQ1(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 if (ct != rct) SETERRQ2(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: SETERRQ1(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 PetscErrorCode ierr; 333 334 PetscFunctionBegin; 335 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 336 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 337 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 338 if (isascii) { 339 const char *name; 340 341 ierr = PetscObjectGetName((PetscObject) tr, &name);CHKERRQ(ierr); 342 ierr = PetscViewerASCIIPrintf(viewer, "Regular refinement %s\n", name ? name : "");CHKERRQ(ierr); 343 } else { 344 SETERRQ1(PetscObjectComm((PetscObject) tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject) viewer)->type_name); 345 } 346 PetscFunctionReturn(0); 347 } 348 349 static PetscErrorCode DMPlexTransformSetUp_Regular(DMPlexTransform tr) 350 { 351 PetscFunctionBegin; 352 PetscFunctionReturn(0); 353 } 354 355 static PetscErrorCode DMPlexTransformDestroy_Regular(DMPlexTransform tr) 356 { 357 DMPlexRefine_Regular *f = (DMPlexRefine_Regular *) tr->data; 358 PetscErrorCode ierr; 359 360 PetscFunctionBegin; 361 ierr = PetscFree(f);CHKERRQ(ierr); 362 PetscFunctionReturn(0); 363 } 364 365 PetscErrorCode DMPlexTransformGetSubcellOrientation_Regular(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 366 { 367 static PetscInt seg_seg[] = {1, -1, 0, -1, 368 0, 0, 1, 0}; 369 static PetscInt tri_seg[] = {2, -1, 1, -1, 0, -1, 370 1, -1, 0, -1, 2, -1, 371 0, -1, 2, -1, 1, -1, 372 0, 0, 1, 0, 2, 0, 373 1, 0, 2, 0, 0, 0, 374 2, 0, 0, 0, 1, 0}; 375 static PetscInt tri_tri[] = {1, -3, 0, -3, 2, -3, 3, -2, 376 0, -2, 2, -2, 1, -2, 3, -1, 377 2, -1, 1, -1, 0, -1, 3, -3, 378 0, 0, 1, 0, 2, 0, 3, 0, 379 1, 1, 2, 1, 0, 1, 3, 1, 380 2, 2, 0, 2, 1, 2, 3, 2}; 381 static PetscInt quad_seg[] = {1, 0, 0, 0, 3, 0, 2, 0, 382 0, 0, 3, 0, 2, 0, 1, 0, 383 3, 0, 2, 0, 1, 0, 0, 0, 384 2, 0, 1, 0, 0, 0, 3, 0, 385 0, 0, 1, 0, 2, 0, 3, 0, 386 1, 0, 2, 0, 3, 0, 0, 0, 387 2, 0, 3, 0, 0, 0, 1, 0, 388 3, 0, 0, 0, 1, 0, 2, 0}; 389 static PetscInt quad_quad[] = {2, -4, 1, -4, 0, -4, 3, -4, 390 1, -3, 0, -3, 3, -3, 2, -3, 391 0, -2, 3, -2, 2, -2, 1, -2, 392 3, -1, 2, -1, 1, -1, 0, -1, 393 0, 0, 1, 0, 2, 0, 3, 0, 394 1, 1, 2, 1, 3, 1, 0, 1, 395 2, 2, 3, 2, 0, 2, 1, 2, 396 3, 3, 0, 3, 1, 3, 2, 3}; 397 static PetscInt tseg_seg[] = {0, -1, 398 0, 0, 399 0, 0, 400 0, -1}; 401 static PetscInt tseg_tseg[] = {1, -2, 0, -2, 402 1, -1, 0, -1, 403 0, 0, 1, 0, 404 0, 1, 1, 1}; 405 static PetscInt tet_seg[] = {0, -1, 406 0, 0, 407 0, 0, 408 0, -1, 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, 0, 422 0, 0, 423 0, -1, 424 0, 0, 425 0, 0, 426 0, -1, 427 0, 0, 428 0, 0}; 429 static PetscInt tet_tri[] = {2, -1, 3, -1, 1, -3, 0, -2, 6, 1, 7, -3, 5, 2, 4, -3, 430 3, -1, 1, -1, 2, -3, 0, -1, 4, 0, 5, 0, 6, 0, 7, 0, 431 1, -1, 2, -1, 3, -3, 0, -3, 4, 0, 5, 0, 6, 0, 7, 0, 432 3, -2, 2, -3, 0, -1, 1, -1, 7, -3, 6, 1, 4, 2, 5, -3, 433 2, -3, 0, -2, 3, -1, 1, -3, 4, 0, 5, 0, 6, 0, 7, 0, 434 0, -2, 3, -2, 2, -2, 1, -2, 4, 0, 5, 0, 6, 0, 7, 0, 435 0, -1, 1, -2, 3, -2, 2, -2, 7, 1, 6, -3, 5, -3, 4, 2, 436 1, -2, 3, -3, 0, -3, 2, -1, 4, 0, 5, 0, 6, 0, 7, 0, 437 3, -3, 0, -1, 1, -1, 2, -3, 4, 0, 5, 0, 6, 0, 7, 0, 438 1, -3, 0, -3, 2, -1, 3, -3, 6, -3, 7, 1, 4, -3, 5, 2, 439 0, -3, 2, -2, 1, -2, 3, -2, 4, 0, 5, 0, 6, 0, 7, 0, 440 2, -2, 1, -3, 0, -2, 3, -1, 4, 0, 5, 0, 6, 0, 7, 0, 441 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 442 1, 0, 2, 2, 0, 1, 3, 1, 4, 0, 5, 0, 6, 0, 7, 0, 443 2, 2, 0, 0, 1, 1, 3, 2, 4, 0, 5, 0, 6, 0, 7, 0, 444 1, 2, 0, 1, 3, 1, 2, 2, 5, 0, 4, 0, 7, -1, 6, -1, 445 0, 1, 3, 0, 1, 0, 2, 0, 4, 0, 5, 0, 6, 0, 7, 0, 446 3, 0, 1, 2, 0, 2, 2, 1, 4, 0, 5, 0, 6, 0, 7, 0, 447 2, 0, 3, 2, 0, 0, 1, 1, 4, -2, 5, -2, 7, 0, 6, 0, 448 3, 2, 0, 2, 2, 1, 1, 2, 4, 0, 5, 0, 6, 0, 7, 0, 449 0, 2, 2, 0, 3, 0, 1, 0, 4, 0, 5, 0, 6, 0, 7, 0, 450 3, 1, 2, 1, 1, 2, 0, 2, 5, -2, 4, -2, 6, -1, 7, -1, 451 2, 1, 1, 1, 3, 2, 0, 0, 4, 0, 5, 0, 6, 0, 7, 0, 452 1, 1, 3, 1, 2, 2, 0, 1, 4, 0, 5, 0, 6, 0, 7, 0}; 453 static PetscInt tet_tet[] = {2, -12, 3, -12, 1, -12, 0, -12, 6, -9, 7, -9, 5, -12, 4, -12, 454 3, -11, 1, -11, 2, -11, 0, -11, 4, 0, 5, 0, 6, 0, 7, 0, 455 1, -10, 2, -10, 3, -10, 0, -10, 4, 0, 5, 0, 6, 0, 7, 0, 456 3, -9, 2, -9, 0, -9, 1, -9, 7, -9, 6, -9, 4, -12, 5, -12, 457 2, -8, 0, -8, 3, -8, 1, -8, 4, 0, 5, 0, 6, 0, 7, 0, 458 0, -7, 3, -7, 2, -7, 1, -7, 4, 0, 5, 0, 6, 0, 7, 0, 459 0, -6, 1, -6, 3, -6, 2, -6, 4, -3, 5, -3, 7, -6, 6, -6, 460 1, -5, 3, -5, 0, -5, 2, -5, 4, 0, 5, 0, 6, 0, 7, 0, 461 3, -4, 0, -4, 1, -4, 2, -4, 4, 0, 5, 0, 6, 0, 7, 0, 462 1, -3, 0, -3, 2, -3, 3, -3, 5, -3, 4, -3, 6, -6, 7, -6, 463 0, -2, 2, -2, 1, -2, 3, -2, 4, 0, 5, 0, 6, 0, 7, 0, 464 2, -1, 1, -1, 0, -1, 3, -1, 4, 0, 5, 0, 6, 0, 7, 0, 465 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 466 1, 1, 2, 1, 0, 1, 3, 1, 4, 0, 5, 0, 6, 0, 7, 0, 467 2, 2, 0, 2, 1, 2, 3, 2, 4, 0, 5, 0, 6, 0, 7, 0, 468 1, 3, 0, 3, 3, 3, 2, 3, 5, 0, 4, 0, 7, 0, 6, 0, 469 0, 4, 3, 4, 1, 4, 2, 4, 4, 0, 5, 0, 6, 0, 7, 0, 470 3, 5, 1, 5, 0, 5, 2, 5, 4, 0, 5, 0, 6, 0, 7, 0, 471 2, 6, 3, 6, 0, 6, 1, 6, 6, 6, 7, 6, 4, 6, 5, 6, 472 3, 7, 0, 7, 2, 7, 1, 7, 4, 0, 5, 0, 6, 0, 7, 0, 473 0, 8, 2, 8, 3, 8, 1, 8, 4, 0, 5, 0, 6, 0, 7, 0, 474 3, 9, 2, 9, 1, 9, 0, 9, 7, 6, 6, 6, 5, 6, 4, 6, 475 2, 10, 1, 10, 3, 10, 0, 10, 4, 0, 5, 0, 6, 0, 7, 0, 476 1, 11, 3, 11, 2, 11, 0, 11, 4, 0, 5, 0, 6, 0, 7, 0}; 477 static PetscInt hex_seg[] = {2, 0, 3, 0, 4, 0, 5, 0, 1, 0, 0, 0, 478 4, 0, 5, 0, 0, 0, 1, 0, 3, 0, 2, 0, 479 5, 0, 4, 0, 1, 0, 0, 0, 3, 0, 2, 0, 480 3, 0, 2, 0, 4, 0, 5, 0, 0, 0, 1, 0, 481 3, 0, 2, 0, 5, 0, 4, 0, 1, 0, 0, 0, 482 4, 0, 5, 0, 1, 0, 0, 0, 2, 0, 3, 0, 483 2, 0, 3, 0, 5, 0, 4, 0, 0, 0, 1, 0, 484 5, 0, 4, 0, 0, 0, 1, 0, 2, 0, 3, 0, 485 4, 0, 5, 0, 3, 0, 2, 0, 1, 0, 0, 0, 486 5, 0, 4, 0, 3, 0, 2, 0, 0, 0, 1, 0, 487 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 488 2, 0, 3, 0, 0, 0, 1, 0, 4, 0, 5, 0, 489 1, 0, 0, 0, 4, 0, 5, 0, 3, 0, 2, 0, 490 1, 0, 0, 0, 5, 0, 4, 0, 2, 0, 3, 0, 491 5, 0, 4, 0, 2, 0, 3, 0, 1, 0, 0, 0, 492 1, 0, 0, 0, 2, 0, 3, 0, 4, 0, 5, 0, 493 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 494 3, 0, 2, 0, 0, 0, 1, 0, 5, 0, 4, 0, 495 1, 0, 0, 0, 3, 0, 2, 0, 5, 0, 4, 0, 496 2, 0, 3, 0, 1, 0, 0, 0, 5, 0, 4, 0, 497 0, 0, 1, 0, 4, 0, 5, 0, 2, 0, 3, 0, 498 0, 0, 1, 0, 3, 0, 2, 0, 4, 0, 5, 0, 499 0, 0, 1, 0, 5, 0, 4, 0, 3, 0, 2, 0, 500 0, 0, 1, 0, 2, 0, 3, 0, 5, 0, 4, 0, 501 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 502 0, 0, 1, 0, 5, 0, 4, 0, 2, 0, 3, 0, 503 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 504 0, 0, 1, 0, 4, 0, 5, 0, 3, 0, 2, 0, 505 2, 0, 3, 0, 1, 0, 0, 0, 4, 0, 5, 0, 506 1, 0, 0, 0, 3, 0, 2, 0, 4, 0, 5, 0, 507 3, 0, 2, 0, 0, 0, 1, 0, 4, 0, 5, 0, 508 4, 0, 5, 0, 2, 0, 3, 0, 1, 0, 0, 0, 509 1, 0, 0, 0, 2, 0, 3, 0, 5, 0, 4, 0, 510 5, 0, 4, 0, 2, 0, 3, 0, 0, 0, 1, 0, 511 1, 0, 0, 0, 5, 0, 4, 0, 3, 0, 2, 0, 512 1, 0, 0, 0, 4, 0, 5, 0, 2, 0, 3, 0, 513 2, 0, 3, 0, 0, 0, 1, 0, 5, 0, 4, 0, 514 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 515 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 516 4, 0, 5, 0, 3, 0, 2, 0, 0, 0, 1, 0, 517 5, 0, 4, 0, 0, 0, 1, 0, 3, 0, 2, 0, 518 2, 0, 3, 0, 5, 0, 4, 0, 1, 0, 0, 0, 519 4, 0, 5, 0, 1, 0, 0, 0, 3, 0, 2, 0, 520 3, 0, 2, 0, 5, 0, 4, 0, 0, 0, 1, 0, 521 3, 0, 2, 0, 4, 0, 5, 0, 1, 0, 0, 0, 522 5, 0, 4, 0, 1, 0, 0, 0, 2, 0, 3, 0, 523 4, 0, 5, 0, 0, 0, 1, 0, 2, 0, 3, 0, 524 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0}; 525 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, 526 8, -2, 9, -2, 10, -2, 11, -2, 3, -4, 0, -4, 1, -4, 2, -4, 7, 0, 4, 0, 5, 0, 6, 0, 527 9, 1, 8, 1, 11, 1, 10, 1, 0, 3, 3, 3, 2, 3, 1, 3, 5, 2, 6, 2, 7, 2, 4, 2, 528 6, 3, 5, 3, 4, 3, 7, 3, 10, -1, 9, -1, 8, -1, 11, -1, 2, -4, 3, -4, 0, -4, 1, -4, 529 4, 1, 7, 1, 6, 1, 5, 1, 11, 2, 8, 2, 9, 2, 10, 2, 1, 3, 0, 3, 3, 3, 2, 3, 530 10, -4, 11, -4, 8, -4, 9, -4, 2, 1, 1, 1, 0, 1, 3, 1, 6, -1, 5, -1, 4, -1, 7, -1, 531 5, -4, 6, -4, 7, -4, 4, -4, 9, 0, 10, 0, 11, 0, 8, 0, 0, -2, 1, -2, 2, -2, 3, -2, 532 11, 3, 10, 3, 9, 3, 8, 3, 1, -2, 2, -2, 3, -2, 0, -2, 4, -3, 7, -3, 6, -3, 5, -3, 533 11, -1, 8, -1, 9, -1, 10, -1, 7, 3, 4, 3, 5, 3, 6, 3, 2, 2, 1, 2, 0, 2, 3, 2, 534 10, 2, 9, 2, 8, 2, 11, 2, 5, 1, 6, 1, 7, 1, 4, 1, 1, -1, 2, -1, 3, -1, 0, -1, 535 5, 2, 4, 2, 7, 2, 6, 2, 1, 2, 0, 2, 3, 2, 2, 2, 10, -4, 9, -4, 8, -4, 11, -4, 536 4, -3, 5, -3, 6, -3, 7, -3, 0, -3, 1, -3, 2, -3, 3, -3, 8, -2, 11, -2, 10, -2, 9, -2, 537 3, 1, 0, 1, 1, 1, 2, 1, 9, -4, 8, -4, 11, -4, 10, -4, 6, 3, 7, 3, 4, 3, 5, 3, 538 1, 3, 2, 3, 3, 3, 0, 3, 10, 1, 11, 1, 8, 1, 9, 1, 5, -4, 4, -4, 7, -4, 6, -4, 539 8, 0, 11, 0, 10, 0, 9, 0, 4, -4, 7, -4, 6, -4, 5, -4, 0, 0, 3, 0, 2, 0, 1, 0, 540 0, 0, 1, 0, 2, 0, 3, 0, 5, -1, 4, -1, 7, -1, 6, -1, 9, -3, 8, -3, 11, -3, 10, -3, 541 9, -3, 10, -3, 11, -3, 8, -3, 6, -2, 5, -2, 4, -2, 7, -2, 3, -3, 0, -3, 1, -3, 2, -3, 542 7, 0, 6, 0, 5, 0, 4, 0, 2, -1, 3, -1, 0, -1, 1, -1, 11, 3, 8, 3, 9, 3, 10, 3, 543 2, 2, 3, 2, 0, 2, 1, 2, 6, 2, 7, 2, 4, 2, 5, 2, 10, 2, 11, 2, 8, 2, 9, 2, 544 6, -1, 7, -1, 4, -1, 5, -1, 3, 0, 2, 0, 1, 0, 0, 0, 9, 1, 10, 1, 11, 1, 8, 1, 545 2, -4, 1, -4, 0, -4, 3, -4, 11, -2, 10, -2, 9, -2, 8, -2, 7, -2, 6, -2, 5, -2, 4, -2, 546 1, -1, 0, -1, 3, -1, 2, -1, 4, 0, 5, 0, 6, 0, 7, 0, 11, -1, 10, -1, 9, -1, 8, -1, 547 0, -2, 3, -2, 2, -2, 1, -2, 8, 3, 9, 3, 10, 3, 11, 3, 4, 1, 5, 1, 6, 1, 7, 1, 548 3, -3, 2, -3, 1, -3, 0, -3, 7, -3, 6, -3, 5, -3, 4, -3, 8, 0, 9, 0, 10, 0, 11, 0, 549 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 8, 0, 9, 0, 10, 0, 11, 0, 550 1, 3, 2, 3, 3, 3, 0, 3, 11, -2, 10, -2, 9, -2, 8, -2, 4, 1, 5, 1, 6, 1, 7, 1, 551 2, 2, 3, 2, 0, 2, 1, 2, 7, -3, 6, -3, 5, -3, 4, -3, 11, -1, 10, -1, 9, -1, 8, -1, 552 3, 1, 0, 1, 1, 1, 2, 1, 8, 3, 9, 3, 10, 3, 11, 3, 7, -2, 6, -2, 5, -2, 4, -2, 553 5, 2, 4, 2, 7, 2, 6, 2, 0, -3, 1, -3, 2, -3, 3, -3, 9, 1, 10, 1, 11, 1, 8, 1, 554 1, -1, 0, -1, 3, -1, 2, -1, 5, -1, 4, -1, 7, -1, 6, -1, 10, 2, 11, 2, 8, 2, 9, 2, 555 4, -3, 5, -3, 6, -3, 7, -3, 1, 2, 0, 2, 3, 2, 2, 2, 11, 3, 8, 3, 9, 3, 10, 3, 556 8, 0, 11, 0, 10, 0, 9, 0, 7, 3, 4, 3, 5, 3, 6, 3, 3, -3, 0, -3, 1, -3, 2, -3, 557 3, -3, 2, -3, 1, -3, 0, -3, 6, 2, 7, 2, 4, 2, 5, 2, 9, -3, 8, -3, 11, -3, 10, -3, 558 9, -3, 10, -3, 11, -3, 8, -3, 5, 1, 6, 1, 7, 1, 4, 1, 0, 0, 3, 0, 2, 0, 1, 0, 559 0, -2, 3, -2, 2, -2, 1, -2, 9, -4, 8, -4, 11, -4, 10, -4, 5, -4, 4, -4, 7, -4, 6, -4, 560 2, -4, 1, -4, 0, -4, 3, -4, 10, 1, 11, 1, 8, 1, 9, 1, 6, 3, 7, 3, 4, 3, 5, 3, 561 7, 0, 6, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 8, -2, 11, -2, 10, -2, 9, -2, 562 6, -1, 7, -1, 4, -1, 5, -1, 2, -1, 3, -1, 0, -1, 1, -1, 10, -4, 9, -4, 8, -4, 11, -4, 563 11, -1, 8, -1, 9, -1, 10, -1, 4, -4, 7, -4, 6, -4, 5, -4, 1, -1, 2, -1, 3, -1, 0, -1, 564 10, 2, 9, 2, 8, 2, 11, 2, 6, -2, 5, -2, 4, -2, 7, -2, 2, 2, 1, 2, 0, 2, 3, 2, 565 8, -2, 9, -2, 10, -2, 11, -2, 0, 3, 3, 3, 2, 3, 1, 3, 4, -3, 7, -3, 6, -3, 5, -3, 566 4, 1, 7, 1, 6, 1, 5, 1, 8, -3, 11, -3, 10, -3, 9, -3, 0, -2, 1, -2, 2, -2, 3, -2, 567 9, 1, 8, 1, 11, 1, 10, 1, 3, -4, 0, -4, 1, -4, 2, -4, 6, -1, 5, -1, 4, -1, 7, -1, 568 5, -4, 6, -4, 7, -4, 4, -4, 10, -1, 9, -1, 8, -1, 11, -1, 1, 3, 0, 3, 3, 3, 2, 3, 569 7, -2, 4, -2, 5, -2, 6, -2, 11, 2, 8, 2, 9, 2, 10, 2, 2, -4, 3, -4, 0, -4, 1, -4, 570 10, -4, 11, -4, 8, -4, 9, -4, 1, -2, 2, -2, 3, -2, 0, -2, 5, 2, 6, 2, 7, 2, 4, 2, 571 11, 3, 10, 3, 9, 3, 8, 3, 2, 1, 1, 1, 0, 1, 3, 1, 7, 0, 4, 0, 5, 0, 6, 0, 572 6, 3, 5, 3, 4, 3, 7, 3, 9, 0, 10, 0, 11, 0, 8, 0, 3, 1, 2, 1, 1, 1, 0, 1}; 573 static PetscInt hex_hex[] = {3, -24, 0, -24, 4, -24, 5, -24, 2, -24, 6, -24, 7, -24, 1, -24, 574 3, -23, 5, -23, 6, -23, 2, -23, 0, -23, 1, -23, 7, -23, 4, -23, 575 4, -22, 0, -22, 1, -22, 7, -22, 5, -22, 6, -22, 2, -22, 3, -22, 576 6, -21, 7, -21, 1, -21, 2, -21, 5, -21, 3, -21, 0, -21, 4, -21, 577 1, -20, 2, -20, 6, -20, 7, -20, 0, -20, 4, -20, 5, -20, 3, -20, 578 6, -19, 2, -19, 3, -19, 5, -19, 7, -19, 4, -19, 0, -19, 1, -19, 579 4, -18, 5, -18, 3, -18, 0, -18, 7, -18, 1, -18, 2, -18, 6, -18, 580 1, -17, 7, -17, 4, -17, 0, -17, 2, -17, 3, -17, 5, -17, 6, -17, 581 2, -16, 3, -16, 5, -16, 6, -16, 1, -16, 7, -16, 4, -16, 0, -16, 582 7, -15, 4, -15, 0, -15, 1, -15, 6, -15, 2, -15, 3, -15, 5, -15, 583 7, -14, 1, -14, 2, -14, 6, -14, 4, -14, 5, -14, 3, -14, 0, -14, 584 0, -13, 4, -13, 5, -13, 3, -13, 1, -13, 2, -13, 6, -13, 7, -13, 585 5, -12, 4, -12, 7, -12, 6, -12, 3, -12, 2, -12, 1, -12, 0, -12, 586 7, -11, 6, -11, 5, -11, 4, -11, 1, -11, 0, -11, 3, -11, 2, -11, 587 0, -10, 1, -10, 7, -10, 4, -10, 3, -10, 5, -10, 6, -10, 2, -10, 588 4, -9, 7, -9, 6, -9, 5, -9, 0, -9, 3, -9, 2, -9, 1, -9, 589 5, -8, 6, -8, 2, -8, 3, -8, 4, -8, 0, -8, 1, -8, 7, -8, 590 2, -7, 6, -7, 7, -7, 1, -7, 3, -7, 0, -7, 4, -7, 5, -7, 591 6, -6, 5, -6, 4, -6, 7, -6, 2, -6, 1, -6, 0, -6, 3, -6, 592 5, -5, 3, -5, 0, -5, 4, -5, 6, -5, 7, -5, 1, -5, 2, -5, 593 2, -4, 1, -4, 0, -4, 3, -4, 6, -4, 5, -4, 4, -4, 7, -4, 594 1, -3, 0, -3, 3, -3, 2, -3, 7, -3, 6, -3, 5, -3, 4, -3, 595 0, -2, 3, -2, 2, -2, 1, -2, 4, -2, 7, -2, 6, -2, 5, -2, 596 3, -1, 2, -1, 1, -1, 0, -1, 5, -1, 4, -1, 7, -1, 6, -1, 597 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 598 1, 1, 2, 1, 3, 1, 0, 1, 7, 1, 4, 1, 5, 1, 6, 1, 599 2, 2, 3, 2, 0, 2, 1, 2, 6, 2, 7, 2, 4, 2, 5, 2, 600 3, 3, 0, 3, 1, 3, 2, 3, 5, 3, 6, 3, 7, 3, 4, 3, 601 4, 4, 0, 4, 3, 4, 5, 4, 7, 4, 6, 4, 2, 4, 1, 4, 602 7, 5, 4, 5, 5, 5, 6, 5, 1, 5, 2, 5, 3, 5, 0, 5, 603 1, 6, 7, 6, 6, 6, 2, 6, 0, 6, 3, 6, 5, 6, 4, 6, 604 3, 7, 2, 7, 6, 7, 5, 7, 0, 7, 4, 7, 7, 7, 1, 7, 605 5, 8, 6, 8, 7, 8, 4, 8, 3, 8, 0, 8, 1, 8, 2, 8, 606 4, 9, 7, 9, 1, 9, 0, 9, 5, 9, 3, 9, 2, 9, 6, 9, 607 4, 10, 5, 10, 6, 10, 7, 10, 0, 10, 1, 10, 2, 10, 3, 10, 608 6, 11, 7, 11, 4, 11, 5, 11, 2, 11, 3, 11, 0, 11, 1, 11, 609 3, 12, 5, 12, 4, 12, 0, 12, 2, 12, 1, 12, 7, 12, 6, 12, 610 6, 13, 2, 13, 1, 13, 7, 13, 5, 13, 4, 13, 0, 13, 3, 13, 611 1, 14, 0, 14, 4, 14, 7, 14, 2, 14, 6, 14, 5, 14, 3, 14, 612 6, 15, 5, 15, 3, 15, 2, 15, 7, 15, 1, 15, 0, 15, 4, 15, 613 0, 16, 4, 16, 7, 16, 1, 16, 3, 16, 2, 16, 6, 16, 5, 16, 614 0, 17, 3, 17, 5, 17, 4, 17, 1, 17, 7, 17, 6, 17, 2, 17, 615 5, 18, 3, 18, 2, 18, 6, 18, 4, 18, 7, 18, 1, 18, 0, 18, 616 7, 19, 6, 19, 2, 19, 1, 19, 4, 19, 0, 19, 3, 19, 5, 19, 617 2, 20, 1, 20, 7, 20, 6, 20, 3, 20, 5, 20, 4, 20, 0, 20, 618 7, 21, 1, 21, 0, 21, 4, 21, 6, 21, 5, 21, 3, 21, 2, 21, 619 2, 22, 6, 22, 5, 22, 3, 22, 1, 22, 0, 22, 4, 22, 7, 22, 620 5, 23, 4, 23, 0, 23, 3, 23, 6, 23, 2, 23, 1, 23, 7, 23}; 621 static PetscInt trip_seg[] = {1, 0, 2, 0, 0, 0, 622 2, 0, 0, 0, 1, 0, 623 0, 0, 1, 0, 2, 0, 624 0, -1, 2, -1, 1, -1, 625 1, -1, 0, -1, 2, -1, 626 2, -1, 1, -1, 0, -1, 627 0, 0, 1, 0, 2, 0, 628 2, 0, 0, 0, 1, 0, 629 1, 0, 2, 0, 0, 0, 630 2, -1, 1, -1, 0, -1, 631 1, -1, 0, -1, 2, -1, 632 0, -1, 2, -1, 1, -1}; 633 static PetscInt trip_tri[] = {1, 1, 2, 1, 0, 1, 3, 1, 634 2, 2, 0, 2, 1, 2, 3, 2, 635 0, 0, 1, 0, 2, 0, 3, 0, 636 2, -1, 1, -1, 0, -1, 3, -3, 637 0, -2, 2, -2, 1, -2, 3, -1, 638 1, -3, 0, -3, 2, -3, 3, -2, 639 0, 0, 1, 0, 2, 0, 3, 0, 640 2, 2, 0, 2, 1, 2, 3, 2, 641 1, 1, 2, 1, 0, 1, 3, 1, 642 1, -3, 0, -3, 2, -3, 3, -2, 643 0, -2, 2, -2, 1, -2, 3, -1, 644 2, -1, 1, -1, 0, -1, 3, -3}; 645 static PetscInt trip_quad[] = {4, -1, 5, -1, 3, -1, 1, -1, 2, -1, 0, -1, 646 5, -1, 3, -1, 4, -1, 2, -1, 0, -1, 1, -1, 647 3, -1, 4, -1, 5, -1, 0, -1, 1, -1, 2, -1, 648 0, -3, 2, -3, 1, -3, 3, -3, 5, -3, 4, -3, 649 1, -3, 0, -3, 2, -3, 4, -3, 3, -3, 5, -3, 650 2, -3, 1, -3, 0, -3, 5, -3, 4, -3, 3, -3, 651 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 652 2, 0, 0, 0, 1, 0, 5, 0, 3, 0, 4, 0, 653 1, 0, 2, 0, 0, 0, 4, 0, 5, 0, 3, 0, 654 5, 2, 4, 2, 3, 2, 2, 2, 1, 2, 0, 2, 655 4, 2, 3, 2, 5, 2, 1, 2, 0, 2, 2, 2, 656 3, 2, 5, 2, 4, 2, 0, 2, 2, 2, 1, 2}; 657 static PetscInt trip_trip[] = {5, -6, 6, -6, 4, -6, 7, -6, 1, -6, 2, -6, 0, -6, 3, -6, 658 6, -5, 4, -5, 5, -5, 7, -5, 2, -5, 0, -5, 1, -5, 3, -5, 659 4, -4, 5, -4, 6, -4, 7, -4, 0, -4, 1, -4, 2, -4, 3, -4, 660 2, -3, 1, -3, 0, -3, 3, -1, 6, -3, 5, -3, 4, -3, 7, -1, 661 0, -2, 2, -2, 1, -2, 3, -3, 4, -2, 6, -2, 5, -2, 7, -3, 662 1, -1, 0, -1, 2, -1, 3, -2, 5, -1, 4, -1, 6, -1, 7, -2, 663 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 664 2, 1, 0, 1, 1, 1, 3, 1, 6, 1, 4, 1, 5, 1, 7, 1, 665 1, 2, 2, 2, 0, 2, 3, 2, 5, 2, 6, 2, 4, 2, 7, 2, 666 5, 3, 4, 3, 6, 3, 7, 4, 1, 3, 0, 3, 2, 3, 3, 4, 667 4, 4, 6, 4, 5, 4, 7, 5, 0, 4, 2, 4, 1, 4, 3, 5, 668 6, 5, 5, 5, 4, 5, 7, 3, 2, 5, 1, 5, 0, 5, 3, 3}; 669 static PetscInt ttri_tseg[] = {2, -2, 1, -2, 0, -2, 670 1, -2, 0, -2, 2, -2, 671 0, -2, 2, -2, 1, -2, 672 2, -1, 1, -1, 0, -1, 673 1, -1, 0, -1, 2, -1, 674 0, -1, 2, -1, 1, -1, 675 0, 0, 1, 0, 2, 0, 676 1, 0, 2, 0, 0, 0, 677 2, 0, 0, 0, 1, 0, 678 0, 1, 1, 1, 2, 1, 679 1, 1, 2, 1, 0, 1, 680 2, 1, 0, 1, 1, 1}; 681 static PetscInt ttri_ttri[] = {1, -6, 0, -6, 2, -6, 3, -5, 682 0, -5, 2, -5, 1, -5, 3, -4, 683 2, -4, 1, -4, 0, -4, 3, -6, 684 1, -3, 0, -3, 2, -3, 3, -2, 685 0, -2, 2, -2, 1, -2, 3, -1, 686 2, -1, 1, -1, 0, -1, 3, -3, 687 0, 0, 1, 0, 2, 0, 3, 0, 688 1, 1, 2, 1, 0, 1, 3, 1, 689 2, 2, 0, 2, 1, 2, 3, 2, 690 0, 3, 1, 3, 2, 3, 3, 3, 691 1, 4, 2, 4, 0, 4, 3, 4, 692 2, 5, 0, 5, 1, 5, 3, 5}; 693 static PetscInt tquad_tvert[] = {0, -1, 694 0, -1, 695 0, -1, 696 0, -1, 697 0, 0, 698 0, 0, 699 0, 0, 700 0, 0, 701 0, 0, 702 0, 0, 703 0, 0, 704 0, 0, 705 0, -1, 706 0, -1, 707 0, -1, 708 0, -1}; 709 static PetscInt tquad_tseg[] = {1, 1, 0, 1, 3, 1, 2, 1, 710 0, 1, 3, 1, 2, 1, 1, 1, 711 3, 1, 2, 1, 1, 1, 0, 1, 712 2, 1, 1, 1, 0, 1, 3, 1, 713 1, 0, 0, 0, 3, 0, 2, 0, 714 0, 0, 3, 0, 2, 0, 1, 0, 715 3, 0, 2, 0, 1, 0, 0, 0, 716 2, 0, 1, 0, 0, 0, 3, 0, 717 0, 0, 1, 0, 2, 0, 3, 0, 718 1, 0, 2, 0, 3, 0, 0, 0, 719 2, 0, 3, 0, 0, 0, 1, 0, 720 3, 0, 0, 0, 1, 0, 2, 0, 721 0, 1, 1, 1, 2, 1, 3, 1, 722 1, 1, 2, 1, 3, 1, 0, 1, 723 2, 1, 3, 1, 0, 1, 1, 1, 724 3, 1, 0, 1, 1, 1, 2, 1}; 725 static PetscInt tquad_tquad[] = {2, -8, 1, -8, 0, -8, 3, -8, 726 1, -7, 0, -7, 3, -7, 2, -7, 727 0, -6, 3, -6, 2, -6, 1, -6, 728 3, -5, 2, -5, 1, -5, 0, -5, 729 2, -4, 1, -4, 0, -4, 3, -4, 730 1, -3, 0, -3, 3, -3, 2, -3, 731 0, -2, 3, -2, 2, -2, 1, -2, 732 3, -1, 2, -1, 1, -1, 0, -1, 733 0, 0, 1, 0, 2, 0, 3, 0, 734 1, 1, 2, 1, 3, 1, 0, 1, 735 2, 2, 3, 2, 0, 2, 1, 2, 736 3, 3, 0, 3, 1, 3, 2, 3, 737 0, 4, 1, 4, 2, 4, 3, 4, 738 1, 5, 2, 5, 3, 5, 0, 5, 739 2, 6, 3, 6, 0, 6, 1, 6, 740 3, 7, 0, 7, 1, 7, 2, 7}; 741 742 PetscFunctionBeginHot; 743 *rnew = r; *onew = o; 744 if (!so) PetscFunctionReturn(0); 745 switch (sct) { 746 case DM_POLYTOPE_POINT: break; 747 case DM_POLYTOPE_POINT_PRISM_TENSOR: 748 *onew = so < 0 ? -(o+1) : o; 749 break; 750 case DM_POLYTOPE_SEGMENT: 751 switch (tct) { 752 case DM_POLYTOPE_POINT: break; 753 case DM_POLYTOPE_SEGMENT: 754 *rnew = seg_seg[(so+1)*4 + r*2]; 755 *onew = DMPolytopeTypeComposeOrientation(tct, o, seg_seg[(so+1)*4 + r*2 + 1]); 756 break; 757 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 758 } 759 break; 760 case DM_POLYTOPE_TRIANGLE: 761 switch (tct) { 762 case DM_POLYTOPE_SEGMENT: 763 *rnew = tri_seg[(so+3)*6 + r*2]; 764 *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so+3)*6 + r*2 + 1]); 765 break; 766 case DM_POLYTOPE_TRIANGLE: 767 *rnew = tri_tri[(so+3)*8 + r*2]; 768 *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so+3)*8 + r*2 + 1]); 769 break; 770 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 771 } 772 break; 773 case DM_POLYTOPE_QUADRILATERAL: 774 switch (tct) { 775 case DM_POLYTOPE_POINT: break; 776 case DM_POLYTOPE_SEGMENT: 777 *rnew = quad_seg[(so+4)*8 + r*2]; 778 *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_seg[(so+4)*8 + r*2 + 1]); 779 break; 780 case DM_POLYTOPE_QUADRILATERAL: 781 *rnew = quad_quad[(so+4)*8 + r*2]; 782 *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_quad[(so+4)*8 + r*2 + 1]); 783 break; 784 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 785 } 786 break; 787 case DM_POLYTOPE_SEG_PRISM_TENSOR: 788 switch (tct) { 789 case DM_POLYTOPE_POINT_PRISM_TENSOR: 790 *rnew = tseg_seg[(so+2)*2 + r*2]; 791 *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_seg[(so+2)*2 + r*2 + 1]); 792 break; 793 case DM_POLYTOPE_SEG_PRISM_TENSOR: 794 *rnew = tseg_tseg[(so+2)*4 + r*2]; 795 *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_tseg[(so+2)*4 + r*2 + 1]); 796 break; 797 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 798 } 799 break; 800 case DM_POLYTOPE_TETRAHEDRON: 801 switch (tct) { 802 case DM_POLYTOPE_SEGMENT: 803 *rnew = tet_seg[(so+12)*2 + r*2]; 804 *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so+12)*2 + r*2 + 1]); 805 break; 806 case DM_POLYTOPE_TRIANGLE: 807 *rnew = tet_tri[(so+12)*16 + r*2]; 808 *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so+12)*16 + r*2 + 1]); 809 break; 810 case DM_POLYTOPE_TETRAHEDRON: 811 *rnew = tet_tet[(so+12)*16 + r*2]; 812 *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so+12)*16 + r*2 + 1]); 813 break; 814 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 815 } 816 break; 817 case DM_POLYTOPE_HEXAHEDRON: 818 switch (tct) { 819 case DM_POLYTOPE_POINT: break; 820 case DM_POLYTOPE_SEGMENT: 821 *rnew = hex_seg[(so+24)*12 + r*2]; 822 *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_seg[(so+24)*12 + r*2 + 1]); 823 break; 824 case DM_POLYTOPE_QUADRILATERAL: 825 *rnew = hex_quad[(so+24)*24 + r*2]; 826 *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_quad[(so+24)*24 + r*2 + 1]); 827 break; 828 case DM_POLYTOPE_HEXAHEDRON: 829 *rnew = hex_hex[(so+24)*16 + r*2]; 830 *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_hex[(so+24)*16 + r*2 + 1]); 831 break; 832 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 833 } 834 break; 835 case DM_POLYTOPE_TRI_PRISM: 836 switch (tct) { 837 case DM_POLYTOPE_SEGMENT: 838 *rnew = trip_seg[(so+6)*6 + r*2]; 839 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_seg[(so+6)*6 + r*2 + 1]); 840 break; 841 case DM_POLYTOPE_TRIANGLE: 842 *rnew = trip_tri[(so+6)*8 + r*2]; 843 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_tri[(so+6)*8 + r*2 + 1]); 844 break; 845 case DM_POLYTOPE_QUADRILATERAL: 846 *rnew = trip_quad[(so+6)*12 + r*2]; 847 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_quad[(so+6)*12 + r*2 + 1]); 848 break; 849 case DM_POLYTOPE_TRI_PRISM: 850 *rnew = trip_trip[(so+6)*16 + r*2]; 851 *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_trip[(so+6)*16 + r*2 + 1]); 852 break; 853 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 854 } 855 break; 856 case DM_POLYTOPE_TRI_PRISM_TENSOR: 857 switch (tct) { 858 case DM_POLYTOPE_SEG_PRISM_TENSOR: 859 *rnew = ttri_tseg[(so+6)*6 + r*2]; 860 *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_tseg[(so+6)*6 + r*2 + 1]); 861 break; 862 case DM_POLYTOPE_TRI_PRISM_TENSOR: 863 *rnew = ttri_ttri[(so+6)*8 + r*2]; 864 *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_ttri[(so+6)*8 + r*2 + 1]); 865 break; 866 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 867 } 868 break; 869 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 870 switch (tct) { 871 case DM_POLYTOPE_POINT_PRISM_TENSOR: 872 *rnew = tquad_tvert[(so+8)*2 + r*2]; 873 *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tvert[(so+8)*2 + r*2 + 1]); 874 break; 875 case DM_POLYTOPE_SEG_PRISM_TENSOR: 876 *rnew = tquad_tseg[(so+8)*8 + r*2]; 877 *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tseg[(so+8)*8 + r*2 + 1]); 878 break; 879 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 880 *rnew = tquad_tquad[(so+8)*8 + r*2]; 881 *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tquad[(so+8)*8 + r*2 + 1]); 882 break; 883 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 884 } 885 break; 886 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]); 887 } 888 PetscFunctionReturn(0); 889 } 890 891 PetscErrorCode DMPlexTransformCellRefine_Regular(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 892 { 893 /* All vertices remain in the refined mesh */ 894 static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; 895 static PetscInt vertexS[] = {1}; 896 static PetscInt vertexC[] = {0}; 897 static PetscInt vertexO[] = {0}; 898 /* Split all edges with a new vertex, making two new 2 edges 899 0--0--0--1--1 900 */ 901 static DMPolytopeType segT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}; 902 static PetscInt segS[] = {1, 2}; 903 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}; 904 static PetscInt segO[] = { 0, 0, 0, 0}; 905 /* Do not split tensor edges */ 906 static DMPolytopeType tvertT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; 907 static PetscInt tvertS[] = {1}; 908 static PetscInt tvertC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 909 static PetscInt tvertO[] = { 0, 0}; 910 /* Add 3 edges inside every triangle, making 4 new triangles. 911 2 912 |\ 913 | \ 914 | \ 915 0 1 916 | C \ 917 | \ 918 | \ 919 2---1---1 920 |\ D / \ 921 1 2 0 0 922 |A \ / B \ 923 0-0-0---1---1 924 */ 925 static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE}; 926 static PetscInt triS[] = {3, 4}; 927 static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 928 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0, 929 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0, 930 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, 931 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, 932 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, 933 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2}; 934 static PetscInt triO[] = {0, 0, 935 0, 0, 936 0, 0, 937 0, -1, 0, 938 0, 0, -1, 939 -1, 0, 0, 940 0, 0, 0}; 941 /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads. 942 3----1----2----0----2 943 | | | 944 0 D 2 C 1 945 | | | 946 3----3----0----1----1 947 | | | 948 1 A 0 B 0 949 | | | 950 0----0----0----1----1 951 */ 952 static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL}; 953 static PetscInt quadS[] = {1, 4, 4}; 954 static PetscInt quadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 955 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 956 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 957 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 958 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, 959 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, 960 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, 961 DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0}; 962 static PetscInt quadO[] = {0, 0, 963 0, 0, 964 0, 0, 965 0, 0, 966 0, 0, -1, 0, 967 0, 0, 0, -1, 968 -1, 0, 0, 0, 969 0, -1, 0, 0}; 970 /* Add 1 edge inside every tensor quad, making 2 new tensor quads 971 2----2----1----3----3 972 | | | 973 | | | 974 | | | 975 4 A 6 B 5 976 | | | 977 | | | 978 | | | 979 0----0----0----1----1 980 */ 981 static DMPolytopeType tsegT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR}; 982 static PetscInt tsegS[] = {1, 2}; 983 static PetscInt tsegC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 984 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, 985 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}; 986 static PetscInt tsegO[] = {0, 0, 987 0, 0, 0, 0, 988 0, 0, 0, 0}; 989 /* Add 1 edge and 8 triangles inside every cell, making 8 new tets 990 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 991 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] 992 called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are 993 [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3] 994 The first four tets just cut off the corners, using the replica notation for new vertices, 995 [v0, (e0, 0), (e2, 0), (e3, 0)] 996 [(e0, 0), v1, (e1, 0), (e4, 0)] 997 [(e2, 0), (e1, 0), v2, (e5, 0)] 998 [(e3, 0), (e4, 0), (e5, 0), v3 ] 999 The next four tets match a vertex to the newly created faces from cutting off those first tets. 1000 [(e2, 0), (e3, 0), (e0, 0), (e5, 0)] 1001 [(e4, 0), (e1, 0), (e0, 0), (e5, 0)] 1002 [(e5, 0), (e0, 0), (e2, 0), (e1, 0)] 1003 [(e5, 0), (e0, 0), (e4, 0), (e3, 0)] 1004 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 1005 [(e2, 0), (e0, 0), (e3, 0)] 1006 [(e0, 0), (e1, 0), (e4, 0)] 1007 [(e2, 0), (e5, 0), (e1, 0)] 1008 [(e3, 0), (e4, 0), (e5, 0)] 1009 The next four, from the second group of tets, are 1010 [(e2, 0), (e0, 0), (e5, 0)] 1011 [(e4, 0), (e0, 0), (e5, 0)] 1012 [(e0, 0), (e1, 0), (e5, 0)] 1013 [(e5, 0), (e3, 0), (e0, 0)] 1014 I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table. 1015 */ 1016 static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON}; 1017 static PetscInt tetS[] = {1, 8, 8}; 1018 static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0, 1019 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, 1020 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, 1021 DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1, 1022 DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, 1023 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1024 DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, 1025 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0, 1026 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0, 1027 DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0, 1028 DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1, 1029 DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0, 1030 DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2, 1031 DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7, 1032 DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6, 1033 DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3, 1034 DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3}; 1035 static PetscInt tetO[] = {0, 0, 1036 0, 0, 0, 1037 0, 0, 0, 1038 0, 0, 0, 1039 0, 0, 0, 1040 0, 0, -1, 1041 0, 0, -1, 1042 0, -1, -1, 1043 0, -1, 0, 1044 0, 0, 0, 0, 1045 0, 0, 0, 0, 1046 0, 0, 0, 0, 1047 0, 0, 0, 0, 1048 -2, 0, 0, -1, 1049 -1, 1, 0, 0, 1050 -1, -1, -3, 2, 1051 -1, 0, -1, 1}; 1052 /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes 1053 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 1054 [v0, v1, v2, v3] f0 bottom 1055 [v4, v5, v6, v7] f1 top 1056 [v0, v3, v5, v4] f2 front 1057 [v2, v1, v7, v6] f3 back 1058 [v3, v2, v6, v5] f4 right 1059 [v0, v4, v7, v1] f5 left 1060 The eight hexes are divided into four on the bottom, and four on the top, 1061 [v0, (e0, 0), (f0, 0), (e3, 0), (e9, 0), (f2, 0), (c0, 0), (f5, 0)] 1062 [(e0, 0), v1, (e1, 0), (f0, 0), (f5, 0), (c0, 0), (f3, 0), (e10, 0)] 1063 [(f0, 0), (e1, 0), v2, (e2, 0), (c0, 0), (f4, 0), (e11, 0), (f3, 0)] 1064 [(e3, 0), (f0, 0), (e2, 0), v3, (f2, 0), (e8, 0), (f4, 0), (c0, 0)] 1065 [(e9, 0), (f5, 0), (c0, 0), (f2, 0), v4, (e4, 0), (f1, 0), (e7, 0)] 1066 [(f2, 0), (c0, 0), (f4, 0), (e8, 0), (e4, 0), v5, (e5, 0), (f1, 0)] 1067 [(c0, 0), (f3, 0), (e11, 0), (f4, 0), (f1, 0), (e5, 0), v6, (e6, 0)] 1068 [(f5, 0), (e10, 0), (f3, 0), (c0, 0), (e7, 0), (f1, 0), (e6, 0), v7] 1069 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, 1070 [(e9, 0), (f2, 0), (c0, 0), (f5, 0)] 1071 [(f5, 0), (c0, 0), (f3, 0), (e10, 0)] 1072 [(c0, 0), (f4, 0), (e11, 0), (f3, 0)] 1073 [(f2, 0), (e8, 0), (f4, 0), (c0, 0)] 1074 and on the x-z plane, 1075 [(f0, 0), (e0, 0), (f5, 0), (c0, 0)] 1076 [(c0, 0), (f5, 0), (e7, 0), (f1, 0)] 1077 [(f4, 0), (c0, 0), (f1, 0), (e5, 0)] 1078 [(e2, 0), (f0, 0), (c0, 0), (f4, 0)] 1079 and on the y-z plane, 1080 [(e3, 0), (f2, 0), (c0, 0), (f0, 0)] 1081 [(f2, 0), (e4, 0), (f1, 0), (c0, 0)] 1082 [(c0, 0), (f1, 0), (e6, 0), (f3, 0)] 1083 [(f0, 0), (c0, 0), (f3, 0), (e1, 0)] 1084 */ 1085 static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 1086 static PetscInt hexS[] = {1, 6, 12, 8}; 1087 static PetscInt hexC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1088 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 1089 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 1090 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 1091 DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0, 1092 DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0, 1093 DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0, 1094 DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2, 1095 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, 1096 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2, 1097 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0, 1098 DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1, 1099 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, 1100 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0, 1101 DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3, 1102 DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, 1103 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, 1104 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, 1105 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, 1106 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, 1107 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, 1108 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, 1109 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, 1110 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, 1111 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, 1112 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}; 1113 static PetscInt hexO[] = {0, 0, 1114 0, 0, 1115 0, 0, 1116 0, 0, 1117 0, 0, 1118 0, 0, 1119 0, 0, -1, -1, 1120 0, -1, -1, 0, 1121 -1, -1, 0, 0, 1122 -1, 0, 0, -1, 1123 -1, 0, 0, -1, 1124 -1, -1, 0, 0, 1125 0, -1, -1, 0, 1126 0, 0, -1, -1, 1127 0, 0, -1, -1, 1128 -1, 0, 0, -1, 1129 -1, -1, 0, 0, 1130 0, -1, -1, 0, 1131 0, 0, 0, 0, -2, 0, 1132 0, 0, -3, 0, -2, 0, 1133 0, 0, -3, 0, 0, 0, 1134 0, 0, 0, 0, 0, 0, 1135 -2, 0, 0, 0, -2, 0, 1136 -2, 0, 0, 0, 0, 0, 1137 -2, 0, -3, 0, 0, 0, 1138 -2, 0, -3, 0, -2, 0}; 1139 /* Add 3 quads inside every triangular prism, making 4 new prisms. */ 1140 static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM}; 1141 static PetscInt tripS[] = {3, 4, 6, 8}; 1142 static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0, 1143 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0, 1144 DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0, 1145 DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1, 1146 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0, 1147 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, 1148 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, 1149 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1150 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, 1151 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0, 1152 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, 1153 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, 1154 DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, 1155 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, 1156 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, 1157 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, 1158 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, 1159 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, 1160 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, 1161 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, 1162 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}; 1163 static PetscInt tripO[] = {0, 0, 1164 0, 0, 1165 0, 0, 1166 0, -1, -1, 1167 -1, 0, -1, 1168 -1, -1, 0, 1169 0, 0, 0, 1170 -1, 0, -1, -1, 1171 -1, 0, -1, -1, 1172 -1, 0, -1, -1, 1173 0, -1, -1, 0, 1174 0, -1, -1, 0, 1175 0, -1, -1, 0, 1176 0, 0, 0, -3, 0, 1177 0, 0, 0, 0, -3, 1178 0, 0, -3, 0, 0, 1179 2, 0, 0, 0, 0, 1180 -2, 0, 0, -3, 0, 1181 -2, 0, 0, 0, -3, 1182 -2, 0, -3, 0, 0, 1183 -2, 0, 0, 0, 0}; 1184 /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms. 1185 2 1186 |\ 1187 | \ 1188 | \ 1189 0---1 1190 1191 2 1192 1193 0 1 1194 1195 2 1196 |\ 1197 | \ 1198 | \ 1199 0---1 1200 */ 1201 static DMPolytopeType ttriT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR}; 1202 static PetscInt ttriS[] = {3, 4}; 1203 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, 1204 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, 1205 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, 1206 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, 1207 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, 1208 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, 1209 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}; 1210 static PetscInt ttriO[] = {0, 0, 0, 0, 1211 0, 0, 0, 0, 1212 0, 0, 0, 0, 1213 0, 0, 0, -1, 0, 1214 0, 0, 0, 0, -1, 1215 0, 0, -1, 0, 0, 1216 0, 0, 0, 0, 0}; 1217 /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */ 1218 static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1219 static PetscInt tquadS[] = {1, 4, 4}; 1220 static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1221 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, 1222 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, 1223 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, 1224 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, 1225 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, 1226 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, 1227 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, 1228 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}; 1229 static PetscInt tquadO[] = {0, 0, 1230 0, 0, 0, 0, 1231 0, 0, 0, 0, 1232 0, 0, 0, 0, 1233 0, 0, 0, 0, 1234 0, 0, 0, 0, -1, 0, 1235 0, 0, 0, 0, 0, -1, 1236 0, 0, -1, 0, 0, 0, 1237 0, 0, 0, -1, 0, 0}; 1238 /* Add 4 edges, 12 triangles, 1 quad, 4 tetrahedra, and 6 pyramids inside every pyramid. */ 1239 static DMPolytopeType tpyrT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TETRAHEDRON, DM_POLYTOPE_PYRAMID}; 1240 static PetscInt tpyrS[] = {4, 12, 1, 4, 6}; 1241 static PetscInt tpyrC[] = {DM_POLYTOPE_POINT, 2, 1, 1, 0, DM_POLYTOPE_POINT, 1, 0, 0, 1242 DM_POLYTOPE_POINT, 2, 2, 1, 0, DM_POLYTOPE_POINT, 1, 0, 0, 1243 DM_POLYTOPE_POINT, 2, 3, 1, 0, DM_POLYTOPE_POINT, 1, 0, 0, 1244 DM_POLYTOPE_POINT, 2, 4, 1, 0, DM_POLYTOPE_POINT, 1, 0, 0, 1245 /* These four triangle face out of the bottom pyramid */ 1246 DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, 1247 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, 1248 DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, 1249 DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3, 1250 /* These eight triangles face out of the corner pyramids */ 1251 DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 1, 2, 1252 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, 1253 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, 1254 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, 1255 DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, 1256 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 1, 1257 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 2, 1258 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 3, 1259 /* This quad faces out of the bottom pyramid */ 1260 DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 1, 1261 /* The bottom face of each tet is on the triangular face */ 1262 DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_TRIANGLE, 0, 8, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 0, 1263 DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 9, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 1, 1264 DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 10, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, 1265 DM_POLYTOPE_TRIANGLE, 1, 4, 3, DM_POLYTOPE_TRIANGLE, 0, 11, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, 1266 /* The front face of all pyramids is toward the front */ 1267 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, 1268 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, 1269 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, 1270 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, 1271 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, 1272 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, 1273 }; 1274 static PetscInt tpyrO[] = {0, 0, 1275 0, 0, 1276 0, 0, 1277 0, 0, 1278 0, 0, -1, 1279 0, 0, -1, 1280 0, 0, -1, 1281 0, 0, -1, 1282 0, -1, 0, 1283 0, -1, 0, 1284 0, -1, 0, 1285 0, -1, 0, 1286 -1, 0, 0, 1287 -1, 0, 0, 1288 -1, 0, 0, 1289 -1, 0, 0, 1290 -1, -1, -1, -1, 1291 0, -3, -2, -3, 1292 0, -3, -2, -3, 1293 0, -3, -2, -3, 1294 0, -3, -2, -3, 1295 0, 0, 0, 0, 0, 1296 0, 0, 0, 0, 0, 1297 0, 0, 0, 0, 0, 1298 0, 0, 0, 0, 0, 1299 -2, 0, 0, 0, 0, 1300 1, 0, 0, 0, 0}; 1301 1302 PetscFunctionBegin; 1303 if (rt) *rt = 0; 1304 switch (source) { 1305 case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break; 1306 case DM_POLYTOPE_SEGMENT: *Nt = 2; *target = segT; *size = segS; *cone = segC; *ornt = segO; break; 1307 case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tvertT; *size = tvertS; *cone = tvertC; *ornt = tvertO; break; 1308 case DM_POLYTOPE_TRIANGLE: *Nt = 2; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 1309 case DM_POLYTOPE_QUADRILATERAL: *Nt = 3; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break; 1310 case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tsegT; *size = tsegS; *cone = tsegC; *ornt = tsegO; break; 1311 case DM_POLYTOPE_TETRAHEDRON: *Nt = 3; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 1312 case DM_POLYTOPE_HEXAHEDRON: *Nt = 4; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break; 1313 case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; 1314 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 2; *target = ttriT; *size = ttriS; *cone = ttriC; *ornt = ttriO; break; 1315 case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break; 1316 case DM_POLYTOPE_PYRAMID: *Nt = 5; *target = tpyrT; *size = tpyrS; *cone = tpyrC; *ornt = tpyrO; break; 1317 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1318 } 1319 PetscFunctionReturn(0); 1320 } 1321 1322 static PetscErrorCode DMPlexTransformInitialize_Regular(DMPlexTransform tr) 1323 { 1324 PetscFunctionBegin; 1325 tr->ops->view = DMPlexTransformView_Regular; 1326 tr->ops->setup = DMPlexTransformSetUp_Regular; 1327 tr->ops->destroy = DMPlexTransformDestroy_Regular; 1328 tr->ops->celltransform = DMPlexTransformCellRefine_Regular; 1329 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Regular; 1330 tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal; 1331 PetscFunctionReturn(0); 1332 } 1333 1334 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform tr) 1335 { 1336 DMPlexRefine_Regular *f; 1337 PetscErrorCode ierr; 1338 1339 PetscFunctionBegin; 1340 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1341 ierr = PetscNewLog(tr, &f);CHKERRQ(ierr); 1342 tr->data = f; 1343 1344 ierr = DMPlexTransformInitialize_Regular(tr);CHKERRQ(ierr); 1345 PetscFunctionReturn(0); 1346 } 1347