1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/petscfeimpl.h> /* For PetscFEInterpolate_Static() */ 3 #include <petscsf.h> 4 5 const char * const DMPlexCellRefinerTypes[] = {"Regular", "ToBox", "ToSimplex", "DMPlexCellRefinerTypes", "DM_REFINER_", 0}; 6 7 /* 8 Note that j and invj are non-square: 9 v0 + j x_face = x_cell 10 invj (x_cell - v0) = x_face 11 */ 12 static PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[]) 13 { 14 /* 15 2 16 |\ 17 | \ 18 | \ 19 | \ 20 | \ 21 | \ 22 | \ 23 2 1 24 | \ 25 | \ 26 | \ 27 0---0-------1 28 v0[Nf][dc]: 3 x 2 29 J[Nf][df][dc]: 3 x 1 x 2 30 invJ[Nf][dc][df]: 3 x 2 x 1 31 detJ[Nf]: 3 32 */ 33 static PetscReal tri_v0[] = {0.0, -1.0, 0.0, 0.0, -1.0, 0.0}; 34 static PetscReal tri_J[] = {1.0, 0.0, -1.0, 1.0, 0.0, -1.0}; 35 static PetscReal tri_invJ[] = {1.0, 0.0, -0.5, 0.5, 0.0, -1.0}; 36 static PetscReal tri_detJ[] = {1.0, 1.414213562373095, 1.0}; 37 /* 38 3---------2---------2 39 | | 40 | | 41 | | 42 3 1 43 | | 44 | | 45 | | 46 0---------0---------1 47 48 v0[Nf][dc]: 4 x 2 49 J[Nf][df][dc]: 4 x 1 x 2 50 invJ[Nf][dc][df]: 4 x 2 x 1 51 detJ[Nf]: 4 52 */ 53 static PetscReal quad_v0[] = {0.0, -1.0, 1.0, 0.0, 0.0, 1.0 -1.0, 0.0}; 54 static PetscReal quad_J[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0}; 55 static PetscReal quad_invJ[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0}; 56 static PetscReal quad_detJ[] = {1.0, 1.0, 1.0, 1.0}; 57 58 PetscFunctionBegin; 59 switch (ct) { 60 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; 61 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; 62 default: 63 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]); 64 } 65 PetscFunctionReturn(0); 66 } 67 68 /* Gets the affine map from the original cell to each subcell */ 69 static PetscErrorCode DMPlexCellRefinerGetAffineTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[]) 70 { 71 /* 72 2 73 |\ 74 | \ 75 | \ 76 | \ 77 | C \ 78 | \ 79 | \ 80 2---1---1 81 |\ D / \ 82 | 2 0 \ 83 |A \ / B \ 84 0---0-------1 85 */ 86 static PetscReal tri_v0[] = {-1.0, -1.0, 0.0, -1.0, -1.0, 0.0, 0.0, -1.0}; 87 static PetscReal tri_J[] = {0.5, 0.0, 88 0.0, 0.5, 89 90 0.5, 0.0, 91 0.0, 0.5, 92 93 0.5, 0.0, 94 0.0, 0.5, 95 96 0.0, -0.5, 97 0.5, 0.5}; 98 static PetscReal tri_invJ[] = {2.0, 0.0, 99 0.0, 2.0, 100 101 2.0, 0.0, 102 0.0, 2.0, 103 104 2.0, 0.0, 105 0.0, 2.0, 106 107 2.0, 2.0, 108 -2.0, 0.0}; 109 /* 110 3---------2---------2 111 | | | 112 | D 2 C | 113 | | | 114 3----3----0----1----1 115 | | | 116 | A 0 B | 117 | | | 118 0---------0---------1 119 */ 120 static PetscReal quad_v0[] = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0}; 121 static PetscReal quad_J[] = {0.5, 0.0, 122 0.0, 0.5, 123 124 0.5, 0.0, 125 0.0, 0.5, 126 127 0.5, 0.0, 128 0.0, 0.5, 129 130 0.5, 0.0, 131 0.0, 0.5}; 132 static PetscReal quad_invJ[] = {2.0, 0.0, 133 0.0, 2.0, 134 135 2.0, 0.0, 136 0.0, 2.0, 137 138 2.0, 0.0, 139 0.0, 2.0, 140 141 2.0, 0.0, 142 0.0, 2.0}; 143 /* 144 Bottom (viewed from top) Top 145 1---------2---------2 7---------2---------6 146 | | | | | | 147 | B 2 C | | H 2 G | 148 | | | | | | 149 3----3----0----1----1 3----3----0----1----1 150 | | | | | | 151 | A 0 D | | E 0 F | 152 | | | | | | 153 0---------0---------3 4---------0---------5 154 */ 155 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, 156 -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0}; 157 static PetscReal hex_J[] = {0.5, 0.0, 0.0, 158 0.0, 0.5, 0.0, 159 0.0, 0.0, 0.5, 160 161 0.5, 0.0, 0.0, 162 0.0, 0.5, 0.0, 163 0.0, 0.0, 0.5, 164 165 0.5, 0.0, 0.0, 166 0.0, 0.5, 0.0, 167 0.0, 0.0, 0.5, 168 169 0.5, 0.0, 0.0, 170 0.0, 0.5, 0.0, 171 0.0, 0.0, 0.5, 172 173 0.5, 0.0, 0.0, 174 0.0, 0.5, 0.0, 175 0.0, 0.0, 0.5, 176 177 0.5, 0.0, 0.0, 178 0.0, 0.5, 0.0, 179 0.0, 0.0, 0.5, 180 181 0.5, 0.0, 0.0, 182 0.0, 0.5, 0.0, 183 0.0, 0.0, 0.5, 184 185 0.5, 0.0, 0.0, 186 0.0, 0.5, 0.0, 187 0.0, 0.0, 0.5}; 188 static PetscReal hex_invJ[] = {2.0, 0.0, 0.0, 189 0.0, 2.0, 0.0, 190 0.0, 0.0, 2.0, 191 192 2.0, 0.0, 0.0, 193 0.0, 2.0, 0.0, 194 0.0, 0.0, 2.0, 195 196 2.0, 0.0, 0.0, 197 0.0, 2.0, 0.0, 198 0.0, 0.0, 2.0, 199 200 2.0, 0.0, 0.0, 201 0.0, 2.0, 0.0, 202 0.0, 0.0, 2.0, 203 204 2.0, 0.0, 0.0, 205 0.0, 2.0, 0.0, 206 0.0, 0.0, 2.0, 207 208 2.0, 0.0, 0.0, 209 0.0, 2.0, 0.0, 210 0.0, 0.0, 2.0, 211 212 2.0, 0.0, 0.0, 213 0.0, 2.0, 0.0, 214 0.0, 0.0, 2.0, 215 216 2.0, 0.0, 0.0, 217 0.0, 2.0, 0.0, 218 0.0, 0.0, 2.0}; 219 220 PetscFunctionBegin; 221 switch (ct) { 222 case DM_POLYTOPE_TRIANGLE: if (Nc) *Nc = 4; if (v0) *v0 = tri_v0; if (J) *J = tri_J; if (invJ) *invJ = tri_invJ; break; 223 case DM_POLYTOPE_QUADRILATERAL: if (Nc) *Nc = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; break; 224 case DM_POLYTOPE_HEXAHEDRON: if (Nc) *Nc = 8; if (v0) *v0 = hex_v0; if (J) *J = hex_J; if (invJ) *invJ = hex_invJ; break; 225 default: 226 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]); 227 } 228 PetscFunctionReturn(0); 229 } 230 231 /* Should this be here or in the DualSpace somehow? */ 232 PetscErrorCode CellRefinerInCellTest_Internal(DMPolytopeType ct, const PetscReal point[], PetscBool *inside) 233 { 234 PetscReal sum = 0.0; 235 PetscInt d; 236 237 PetscFunctionBegin; 238 *inside = PETSC_TRUE; 239 switch (ct) { 240 case DM_POLYTOPE_TRIANGLE: 241 case DM_POLYTOPE_TETRAHEDRON: 242 for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d) { 243 if (point[d] < -1.0) {*inside = PETSC_FALSE; break;} 244 sum += point[d]; 245 } 246 if (sum > PETSC_SMALL) {*inside = PETSC_FALSE; break;} 247 break; 248 case DM_POLYTOPE_QUADRILATERAL: 249 case DM_POLYTOPE_HEXAHEDRON: 250 for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d) 251 if (PetscAbsReal(point[d]) > 1.+PETSC_SMALL) {*inside = PETSC_FALSE; break;} 252 break; 253 default: 254 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]); 255 } 256 PetscFunctionReturn(0); 257 } 258 259 /* Regular Refinment of Hybrid Meshes 260 261 We would like to express regular refinement as a small set of rules that can be applied on every point of the Plex 262 to automatically generate a refined Plex. In fact, we would like these rules to be general enough to encompass other 263 transformations, such as changing from one type of cell to another, as simplex to hex. 264 265 To start, we can create a function that takes an original cell type and returns the number of new cells replacing it 266 and the types of the new cells. 267 268 We need the group multiplication table for group actions from the dihedral group for each cell type. 269 270 We need an operator which takes in a cell, and produces a new set of cells with new faces and correct orientations. I think 271 we can just write this operator for faces with identity, and then compose the face orientation actions to get the actual 272 (face, orient) pairs for each subcell. 273 */ 274 275 /* 276 Input Parameters: 277 + ct - The type of the input cell 278 . co - The orientation of this cell in it parent 279 - cp - The requested cone point of this cell assuming orientation 0 280 281 Output Parameters: 282 . cpnew - The new cone point, taking inout account the orientation co 283 */ 284 PETSC_STATIC_INLINE PetscErrorCode DMPolytopeMapCell(DMPolytopeType ct, PetscInt co, PetscInt cp, PetscInt *cpnew) 285 { 286 const PetscInt csize = DMPolytopeTypeGetConeSize(ct); 287 288 PetscFunctionBeginHot; 289 if (ct == DM_POLYTOPE_POINT) {*cpnew = cp;} 290 else {*cpnew = (co < 0 ? -(co+1)-cp+csize : co+cp)%csize;} 291 PetscFunctionReturn(0); 292 } 293 294 static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[]) 295 { 296 static PetscReal seg_v[] = {-1.0, 0.0, 1.0}; 297 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}; 298 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}; 299 static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0, 300 -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0, 301 -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0}; 302 static PetscReal hex_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0, 303 -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, -1.0, 304 -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, 1.0, -1.0, 305 -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 1.0, -1.0, 0.0, 306 -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 307 -1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 308 -1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, 1.0, 309 -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 310 -1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0}; 311 312 PetscFunctionBegin; 313 switch (ct) { 314 case DM_POLYTOPE_SEGMENT: *Nv = 3; *subcellV = seg_v; break; 315 case DM_POLYTOPE_TRIANGLE: *Nv = 6; *subcellV = tri_v; break; 316 case DM_POLYTOPE_QUADRILATERAL: *Nv = 9; *subcellV = quad_v; break; 317 case DM_POLYTOPE_TETRAHEDRON: *Nv = 10; *subcellV = tet_v; break; 318 case DM_POLYTOPE_HEXAHEDRON: *Nv = 27; *subcellV = hex_v; break; 319 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]); 320 } 321 PetscFunctionReturn(0); 322 } 323 324 static PetscErrorCode DMPlexCellRefinerGetCellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[]) 325 { 326 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, -1.0/3.0, -1.0/3.0}; 327 static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0, 328 -1.0, 0.0, -1.0, -1.0/3.0, -1.0/3.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0, 329 -1.0, -1.0, 0.0, -1.0/3.0, -1.0, -1.0/3.0, 0.0, -1.0, 0.0, 330 -1.0, -1.0/3.0, -1.0/3.0, -1.0/3.0, -1.0/3.0, -1.0/3.0, -1.0, 0.0, 0.0, 331 -1.0, -1.0, 1.0, -0.5, -0.5, -0.5}; 332 PetscErrorCode ierr; 333 334 PetscFunctionBegin; 335 switch (ct) { 336 case DM_POLYTOPE_SEGMENT: 337 case DM_POLYTOPE_QUADRILATERAL: 338 case DM_POLYTOPE_HEXAHEDRON: 339 ierr = DMPlexCellRefinerGetCellVertices_Regular(cr, ct, Nv, subcellV);CHKERRQ(ierr);break; 340 case DM_POLYTOPE_TRIANGLE: *Nv = 7; *subcellV = tri_v; break; 341 case DM_POLYTOPE_TETRAHEDRON: *Nv = 15; *subcellV = tet_v; break; 342 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]); 343 } 344 PetscFunctionReturn(0); 345 } 346 347 /* 348 DMPlexCellRefinerGetCellVertices - Get the set of refined vertices lying in the closure of a reference cell of given type 349 350 Input Parameters: 351 + cr - The DMPlexCellRefiner object 352 - ct - The cell type 353 354 Output Parameters: 355 + Nv - The number of refined vertices in the closure of the reference cell of given type 356 - subcellV - The coordinates of these vertices in the reference cell 357 358 Level: developer 359 360 .seealso: DMPlexCellRefinerGetSubcellVertices() 361 */ 362 static PetscErrorCode DMPlexCellRefinerGetCellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[]) 363 { 364 PetscErrorCode ierr; 365 366 PetscFunctionBegin; 367 ierr = (*cr->ops->getcellvertices)(cr, ct, Nv, subcellV);CHKERRQ(ierr); 368 PetscFunctionReturn(0); 369 } 370 371 static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[]) 372 { 373 static PetscInt seg_v[] = {0, 1, 1, 2}; 374 static PetscInt tri_v[] = {0, 3, 5, 3, 1, 4, 5, 4, 2, 3, 4, 5}; 375 static PetscInt quad_v[] = {0, 4, 8, 7, 4, 1, 5, 8, 8, 5, 2, 6, 7, 8, 6, 3}; 376 static PetscInt tet_v[] = {0, 3, 1, 6, 3, 2, 4, 8, 1, 4, 5, 7, 6, 8, 7, 9, 377 1, 6, 3, 7, 8, 4, 3, 7, 7, 3, 1, 4, 7, 3, 8, 6}; 378 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, 379 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}; 380 381 PetscFunctionBegin; 382 if (ct != rct) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 383 switch (ct) { 384 case DM_POLYTOPE_SEGMENT: *Nv = 2; *subcellV = &seg_v[r*(*Nv)]; break; 385 case DM_POLYTOPE_TRIANGLE: *Nv = 3; *subcellV = &tri_v[r*(*Nv)]; break; 386 case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break; 387 case DM_POLYTOPE_TETRAHEDRON: *Nv = 4; *subcellV = &tet_v[r*(*Nv)]; break; 388 case DM_POLYTOPE_HEXAHEDRON: *Nv = 8; *subcellV = &hex_v[r*(*Nv)]; break; 389 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]); 390 } 391 PetscFunctionReturn(0); 392 } 393 394 static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[]) 395 { 396 static PetscInt tri_v[] = {0, 3, 6, 5, 3, 1, 4, 6, 5, 6, 4, 2}; 397 static PetscInt tet_v[] = {0, 3, 4, 1, 7, 8, 14, 10, 6, 12, 11, 5, 3, 4, 14, 10, 2, 5, 11, 9, 1, 8, 14, 4, 13, 12 , 10, 7, 9, 8, 14, 11}; 398 PetscErrorCode ierr; 399 400 PetscFunctionBegin; 401 switch (ct) { 402 case DM_POLYTOPE_SEGMENT: 403 case DM_POLYTOPE_QUADRILATERAL: 404 case DM_POLYTOPE_HEXAHEDRON: 405 ierr = DMPlexCellRefinerGetSubcellVertices_Regular(cr, ct, rct, r, Nv, subcellV);CHKERRQ(ierr);break; 406 case DM_POLYTOPE_TRIANGLE: 407 if (rct != DM_POLYTOPE_QUADRILATERAL) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 408 *Nv = 4; *subcellV = &tri_v[r*(*Nv)]; break; 409 case DM_POLYTOPE_TETRAHEDRON: 410 if (rct != DM_POLYTOPE_HEXAHEDRON) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 411 *Nv = 8; *subcellV = &tet_v[r*(*Nv)]; break; 412 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]); 413 } 414 PetscFunctionReturn(0); 415 } 416 417 /* 418 DMPlexCellRefinerGetSubcellVertices - Get the set of refined vertices defining a subcell in the reference cell of given type 419 420 Input Parameters: 421 + cr - The DMPlexCellRefiner object 422 . ct - The cell type 423 . rct - The type of subcell 424 - r - The subcell index 425 426 Output Parameters: 427 + Nv - The number of refined vertices in the subcell 428 - subcellV - The indices of these vertices in the set of vertices returned by DMPlexCellRefinerGetCellVertices() 429 430 Level: developer 431 432 .seealso: DMPlexCellRefinerGetCellVertices() 433 */ 434 static PetscErrorCode DMPlexCellRefinerGetSubcellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[]) 435 { 436 PetscErrorCode ierr; 437 438 PetscFunctionBegin; 439 ierr = (*cr->ops->getsubcellvertices)(cr, ct, rct, r, Nv, subcellV);CHKERRQ(ierr); 440 PetscFunctionReturn(0); 441 } 442 443 /* 444 Input Parameters: 445 + cr - The DMPlexCellRefiner 446 . pct - The cell type of the parent, from whom the new cell is being produced 447 . po - The orientation of the parent cell in its enclosing parent 448 . ct - The type being produced 449 . r - The replica number requested for the produced cell type 450 - o - The relative orientation of the replica 451 452 Output Parameters: 453 + rnew - The replica number, given the orientation of of the parent 454 - onew - The replica orientation, given the orientation of the parent 455 */ 456 static PetscErrorCode DMPlexCellRefinerMapSubcells(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 457 { 458 PetscErrorCode ierr; 459 460 PetscFunctionBeginHot; 461 ierr = (*cr->ops->mapsubcells)(cr, pct, po, ct, r, o, rnew, onew);CHKERRQ(ierr); 462 PetscFunctionReturn(0); 463 } 464 465 static PetscErrorCode DMPlexCellRefinerMapSubcells_Regular(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 466 { 467 /* We shift any input orientation in order to make it non-negative 468 The orientation array o[po][o] gives the orientation the new replica rnew has to have in order to reproduce the face sequence from (r, o) 469 The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po 470 Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po 471 */ 472 PetscInt tri_seg_o[] = {-2, 0, 473 -2, 0, 474 -2, 0, 475 0, -2, 476 0, -2, 477 0, -2}; 478 PetscInt tri_seg_r[] = {1, 0, 2, 479 0, 2, 1, 480 2, 1, 0, 481 0, 1, 2, 482 1, 2, 0, 483 2, 0, 1}; 484 PetscInt tri_tri_o[] = {0, 1, 2, -3, -2, -1, 485 2, 0, 1, -2, -1, -3, 486 1, 2, 0, -1, -3, -2, 487 -3, -2, -1, 0, 1, 2, 488 -1, -3, -2, 1, 2, 0, 489 -2, -1, -3, 2, 0, 1}; 490 /* orientation if the replica is the center triangle */ 491 PetscInt tri_tri_o_c[] = {2, 0, 1, -2, -1, -3, 492 1, 2, 0, -1, -3, -2, 493 0, 1, 2, -3, -2, -1, 494 -3, -2, -1, 0, 1, 2, 495 -1, -3, -2, 1, 2, 0, 496 -2, -1, -3, 2, 0, 1}; 497 PetscInt tri_tri_r[] = {0, 2, 1, 3, 498 2, 1, 0, 3, 499 1, 0, 2, 3, 500 0, 1, 2, 3, 501 1, 2, 0, 3, 502 2, 0, 1, 3}; 503 PetscInt quad_seg_r[] = {3, 2, 1, 0, 504 2, 1, 0, 3, 505 1, 0, 3, 2, 506 0, 3, 2, 1, 507 0, 1, 2, 3, 508 1, 2, 3, 0, 509 2, 3, 0, 1, 510 3, 0, 1, 2}; 511 PetscInt quad_quad_o[] = { 0, 1, 2, 3, -4, -3, -2, -1, 512 4, 0, 1, 2, -3, -2, -1, -4, 513 3, 4, 0, 1, -2, -1, -4, -3, 514 2, 3, 4, 0, -1, -4, -3, -2, 515 -4, -3, -2, -1, 0, 1, 2, 3, 516 -1, -4, -3, -2, 1, 2, 3, 0, 517 -2, -1, -4, -3, 2, 3, 0, 1, 518 -3, -2, -1, -4, 3, 0, 1, 2}; 519 PetscInt quad_quad_r[] = {0, 3, 2, 1, 520 3, 2, 1, 0, 521 2, 1, 0, 3, 522 1, 0, 3, 2, 523 0, 1, 2, 3, 524 1, 2, 3, 0, 525 2, 3, 0, 1, 526 3, 0, 1, 2}; 527 PetscInt tquad_tquad_o[] = { 0, 1, -2, -1, 528 1, 0, -1, -2, 529 -2, -1, 0, 1, 530 -1, -2, 1, 0}; 531 PetscInt tquad_tquad_r[] = {1, 0, 532 1, 0, 533 0, 1, 534 0, 1}; 535 536 PetscFunctionBeginHot; 537 /* The default is no transformation */ 538 *rnew = r; 539 *onew = o; 540 switch (pct) { 541 case DM_POLYTOPE_SEGMENT: 542 if (ct == DM_POLYTOPE_SEGMENT) { 543 if (po == 0 || po == -1) {*rnew = r; *onew = o;} 544 else if (po == 1 || po == -2) {*rnew = (r+1)%2; *onew = (o == 0 || o == -1) ? -2 : 0;} 545 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for segment", po); 546 } 547 break; 548 case DM_POLYTOPE_TRIANGLE: 549 switch (ct) { 550 case DM_POLYTOPE_SEGMENT: 551 if (o == -1) o = 0; 552 if (o == -2) o = 1; 553 *onew = tri_seg_o[(po+3)*2+o]; 554 *rnew = tri_seg_r[(po+3)*3+r]; 555 break; 556 case DM_POLYTOPE_TRIANGLE: 557 *onew = r == 3 && po < 0 ? tri_tri_o_c[((po+3)%3)*6+o+3] : tri_tri_o[(po+3)*6+o+3]; 558 *rnew = tri_tri_r[(po+3)*4+r]; 559 break; 560 default: break; 561 } 562 break; 563 case DM_POLYTOPE_QUADRILATERAL: 564 switch (ct) { 565 case DM_POLYTOPE_SEGMENT: 566 *onew = o; 567 *rnew = quad_seg_r[(po+4)*4+r]; 568 break; 569 case DM_POLYTOPE_QUADRILATERAL: 570 *onew = quad_quad_o[(po+4)*8+o+4]; 571 *rnew = quad_quad_r[(po+4)*4+r]; 572 break; 573 default: break; 574 } 575 break; 576 case DM_POLYTOPE_SEG_PRISM_TENSOR: 577 switch (ct) { 578 /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */ 579 case DM_POLYTOPE_SEG_PRISM_TENSOR: 580 *onew = tquad_tquad_o[(po+2)*4+o+2]; 581 *rnew = tquad_tquad_r[(po+2)*2+r]; 582 break; 583 default: break; 584 } 585 break; 586 default: break; 587 } 588 PetscFunctionReturn(0); 589 } 590 591 static PetscErrorCode DMPlexCellRefinerMapSubcells_ToBox(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 592 { 593 PetscErrorCode ierr; 594 /* We shift any input orientation in order to make it non-negative 595 The orientation array o[po][o] gives the orientation the new replica rnew has to have in order to reproduce the face sequence from (r, o) 596 The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po 597 Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po 598 */ 599 PetscInt tri_seg_o[] = {0, -2, 600 0, -2, 601 0, -2, 602 0, -2, 603 0, -2, 604 0, -2}; 605 PetscInt tri_seg_r[] = {2, 1, 0, 606 1, 0, 2, 607 0, 2, 1, 608 0, 1, 2, 609 1, 2, 0, 610 2, 0, 1}; 611 PetscInt tri_tri_o[] = {0, 1, 2, 3, -4, -3, -2, -1, 612 3, 0, 1, 2, -3, -2, -1, -4, 613 1, 2, 3, 0, -1, -4, -3, -2, 614 -4, -3, -2, -1, 0, 1, 2, 3, 615 -1, -4, -3, -2, 1, 2, 3, 0, 616 -3, -2, -1, -4, 3, 0, 1, 2}; 617 PetscInt tri_tri_r[] = {0, 2, 1, 618 2, 1, 0, 619 1, 0, 2, 620 0, 1, 2, 621 1, 2, 0, 622 2, 0, 1}; 623 PetscInt tquad_tquad_o[] = { 0, 1, -2, -1, 624 1, 0, -1, -2, 625 -2, -1, 0, 1, 626 -1, -2, 1, 0}; 627 PetscInt tquad_tquad_r[] = {1, 0, 628 1, 0, 629 0, 1, 630 0, 1}; 631 PetscInt tquad_quad_o[] = {-2, -1, -4, -3, 2, 3, 0, 1, 632 1, 2, 3, 0, -1, -4, -3, -2, 633 -4, -3, -2, -1, 0, 1, 2, 3, 634 1, 0, 3, 2, -3, -4, -1, -2}; 635 636 PetscFunctionBeginHot; 637 *rnew = r; 638 *onew = o; 639 switch (pct) { 640 case DM_POLYTOPE_TRIANGLE: 641 switch (ct) { 642 case DM_POLYTOPE_SEGMENT: 643 if (o == -1) o = 0; 644 if (o == -2) o = 1; 645 *onew = tri_seg_o[(po+3)*2+o]; 646 *rnew = tri_seg_r[(po+3)*3+r]; 647 break; 648 case DM_POLYTOPE_QUADRILATERAL: 649 *onew = tri_tri_o[(po+3)*8+o+4]; 650 /* TODO See sheet, for fixing po == 1 and 2 */ 651 if (po == 2 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4]; 652 if (po == 2 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)]; 653 if (po == 1 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4]; 654 if (po == 1 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)]; 655 if (po == -1 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4]; 656 if (po == -1 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)]; 657 if (po == -2 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4]; 658 if (po == -2 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)]; 659 *rnew = tri_tri_r[(po+3)*3+r]; 660 break; 661 default: break; 662 } 663 break; 664 case DM_POLYTOPE_SEG_PRISM_TENSOR: 665 switch (ct) { 666 /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */ 667 case DM_POLYTOPE_SEG_PRISM_TENSOR: 668 *onew = tquad_tquad_o[(po+2)*4+o+2]; 669 *rnew = tquad_tquad_r[(po+2)*2+r]; 670 break; 671 case DM_POLYTOPE_QUADRILATERAL: 672 *onew = tquad_quad_o[(po+2)*8+o+4]; 673 *rnew = tquad_tquad_r[(po+2)*2+r]; 674 break; 675 default: break; 676 } 677 break; 678 default: 679 ierr = DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew);CHKERRQ(ierr); 680 } 681 PetscFunctionReturn(0); 682 } 683 684 static PetscErrorCode DMPlexCellRefinerMapSubcells_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 685 { 686 return DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew); 687 } 688 689 /*@ 690 DMPlexCellRefinerRefine - Return a description of the refinement for a given cell type 691 692 Input Parameter: 693 . source - The cell type for a source point 694 695 Output Parameter: 696 + Nt - The number of cell types generated by refinement 697 . target - The cell types generated 698 . size - The number of subcells of each type, ordered by dimension 699 . cone - A list of the faces for each subcell of the same type as source 700 - ornt - A list of the face orientations for each subcell of the same type as source 701 702 Note: The cone array gives the cone of each subcell listed by the first three outputs. For the each cone point, we 703 need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical 704 division (described in these outputs) of the cell in the original mesh. The point identifier is given by 705 $ the number of cones to be taken, or 0 for the current cell 706 $ the cell cone point number at each level from which it is subdivided 707 $ the replica number r of the subdivision. 708 The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is 709 $ Nt = 2 710 $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT} 711 $ size = {1, 2} 712 $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0} 713 $ ornt = { 0, 0, 0, 0} 714 715 Level: developer 716 717 .seealso: DMPlexCellRefinerCreate(), DMPlexRefineUniform() 718 @*/ 719 PetscErrorCode DMPlexCellRefinerRefine(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 720 { 721 PetscErrorCode ierr; 722 723 PetscFunctionBeginHot; 724 ierr = (*cr->ops->refine)(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 725 PetscFunctionReturn(0); 726 } 727 728 static PetscErrorCode DMPlexCellRefinerRefine_Regular(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 729 { 730 /* All vertices remain in the refined mesh */ 731 static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; 732 static PetscInt vertexS[] = {1}; 733 static PetscInt vertexC[] = {0}; 734 static PetscInt vertexO[] = {0}; 735 /* Split all edges with a new vertex, making two new 2 edges 736 0--0--0--1--1 737 */ 738 static DMPolytopeType edgeT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}; 739 static PetscInt edgeS[] = {1, 2}; 740 static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 741 static PetscInt edgeO[] = { 0, 0, 0, 0}; 742 /* Do not split tensor edges */ 743 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; 744 static PetscInt tedgeS[] = {1}; 745 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 746 static PetscInt tedgeO[] = { 0, 0}; 747 /* Add 3 edges inside every triangle, making 4 new triangles. 748 2 749 |\ 750 | \ 751 | \ 752 0 1 753 | C \ 754 | \ 755 | \ 756 2---1---1 757 |\ D / \ 758 1 2 0 0 759 |A \ / B \ 760 0-0-0---1---1 761 */ 762 static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE}; 763 static PetscInt triS[] = {3, 4}; 764 static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 765 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0, 766 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0, 767 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, 768 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, 769 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, 770 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2}; 771 static PetscInt triO[] = {0, 0, 772 0, 0, 773 0, 0, 774 0, -2, 0, 775 0, 0, -2, 776 -2, 0, 0, 777 0, 0, 0}; 778 /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads. 779 3----1----2----0----2 780 | | | 781 0 D 2 C 1 782 | | | 783 3----3----0----1----1 784 | | | 785 1 A 0 B 0 786 | | | 787 0----0----0----1----1 788 */ 789 static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL}; 790 static PetscInt quadS[] = {1, 4, 4}; 791 static PetscInt quadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 792 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 793 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 794 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 795 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, 796 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, 797 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, 798 DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0}; 799 static PetscInt quadO[] = {0, 0, 800 0, 0, 801 0, 0, 802 0, 0, 803 0, 0, -2, 0, 804 0, 0, 0, -2, 805 -2, 0, 0, 0, 806 0, -2, 0, 0}; 807 /* Add 1 edge inside every tensor quad, making 2 new tensor quads 808 2----2----1----3----3 809 | | | 810 | | | 811 | | | 812 4 A 6 B 5 813 | | | 814 | | | 815 | | | 816 0----0----0----1----1 817 */ 818 static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR}; 819 static PetscInt tquadS[] = {1, 2}; 820 static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 821 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, 822 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}; 823 static PetscInt tquadO[] = {0, 0, 824 0, 0, 0, 0, 825 0, 0, 0, 0}; 826 /* Add 1 edge and 8 triangles inside every cell, making 8 new tets 827 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 828 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] 829 called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are 830 [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3] 831 The first four tets just cut off the corners, using the replica notation for new vertices, 832 [v0, (e0, 0), (e2, 0), (e3, 0)] 833 [(e0, 0), v1, (e1, 0), (e4, 0)] 834 [(e2, 0), (e1, 0), v2, (e5, 0)] 835 [(e3, 0), (e4, 0), (e5, 0), v3 ] 836 The next four tets match a vertex to the newly created faces from cutting off those first tets. 837 [(e2, 0), (e3, 0), (e0, 0), (e5, 0)] 838 [(e4, 0), (e1, 0), (e0, 0), (e5, 0)] 839 [(e5, 0), (e0, 0), (e2, 0), (e1, 0)] 840 [(e5, 0), (e0, 0), (e4, 0), (e3, 0)] 841 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 842 [(e2, 0), (e0, 0), (e3, 0)] 843 [(e0, 0), (e1, 0), (e4, 0)] 844 [(e2, 0), (e5, 0), (e1, 0)] 845 [(e3, 0), (e4, 0), (e5, 0)] 846 The next four, from the second group of tets, are 847 [(e2, 0), (e0, 0), (e5, 0)] 848 [(e4, 0), (e0, 0), (e5, 0)] 849 [(e0, 0), (e1, 0), (e5, 0)] 850 [(e5, 0), (e3, 0), (e0, 0)] 851 I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table. 852 */ 853 static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON}; 854 static PetscInt tetS[] = {1, 8, 8}; 855 static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0, 856 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, 857 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, 858 DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1, 859 DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, 860 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 861 DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, 862 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0, 863 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0, 864 DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0, 865 DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1, 866 DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0, 867 DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2, 868 DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7, 869 DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6, 870 DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3, 871 DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3}; 872 static PetscInt tetO[] = {0, 0, 873 0, 0, 0, 874 0, 0, 0, 875 0, 0, 0, 876 0, 0, 0, 877 0, 0, -2, 878 0, 0, -2, 879 0, -2, -2, 880 0, -2, 0, 881 0, 0, 0, 0, 882 0, 0, 0, 0, 883 0, 0, 0, 0, 884 0, 0, 0, 0, 885 -3, 0, 0, -2, 886 -2, 1, 0, 0, 887 -2, -2, -1, 2, 888 -2, 0, -2, 1}; 889 /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes 890 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 891 [v0, v1, v2, v3] f0 bottom 892 [v4, v5, v6, v7] f1 top 893 [v0, v3, v5, v4] f2 front 894 [v2, v1, v7, v6] f3 back 895 [v3, v2, v6, v5] f4 right 896 [v0, v4, v7, v1] f5 left 897 The eight hexes are divided into four on the bottom, and four on the top, 898 [v0, (e0, 0), (f0, 0), (e3, 0), (e9, 0), (f2, 0), (c0, 0), (f5, 0)] 899 [(e0, 0), v1, (e1, 0), (f0, 0), (f5, 0), (c0, 0), (f3, 0), (e10, 0)] 900 [(f0, 0), (e1, 0), v2, (e2, 0), (c0, 0), (f4, 0), (e11, 0), (f3, 0)] 901 [(e3, 0), (f0, 0), (e2, 0), v3, (f2, 0), (e8, 0), (f4, 0), (c0, 0)] 902 [(e9, 0), (f5, 0), (c0, 0), (f2, 0), v4, (e4, 0), (f1, 0), (e7, 0)] 903 [(f2, 0), (c0, 0), (f4, 0), (e8, 0), (e4, 0), v5, (e5, 0), (f1, 0)] 904 [(c0, 0), (f3, 0), (e11, 0), (f4, 0), (f1, 0), (e5, 0), v6, (e6, 0)] 905 [(f5, 0), (e10, 0), (f3, 0), (c0, 0), (e7, 0), (f1, 0), (e6, 0), v7] 906 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, 907 [(e9, 0), (f2, 0), (c0, 0), (f5, 0)] 908 [(f5, 0), (c0, 0), (f3, 0), (e10, 0)] 909 [(c0, 0), (f4, 0), (e11, 0), (f3, 0)] 910 [(f2, 0), (e8, 0), (f4, 0), (c0, 0)] 911 and on the x-z plane, 912 [(f0, 0), (e0, 0), (f5, 0), (c0, 0)] 913 [(c0, 0), (f5, 0), (e7, 0), (f1, 0)] 914 [(f4, 0), (c0, 0), (f1, 0), (e5, 0)] 915 [(e2, 0), (f0, 0), (c0, 0), (f4, 0)] 916 and on the y-z plane, 917 [(e3, 0), (f2, 0), (c0, 0), (f0, 0)] 918 [(f2, 0), (e4, 0), (f1, 0), (c0, 0)] 919 [(c0, 0), (f1, 0), (e6, 0), (f3, 0)] 920 [(f0, 0), (c0, 0), (f3, 0), (e1, 0)] 921 */ 922 static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 923 static PetscInt hexS[] = {1, 6, 12, 8}; 924 static PetscInt hexC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 925 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 926 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 927 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 928 DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0, 929 DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0, 930 DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0, 931 DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2, 932 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, 933 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2, 934 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0, 935 DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1, 936 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, 937 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0, 938 DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3, 939 DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, 940 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, 941 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, 942 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, 943 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, 944 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, 945 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, 946 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, 947 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, 948 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, 949 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}; 950 static PetscInt hexO[] = {0, 0, 951 0, 0, 952 0, 0, 953 0, 0, 954 0, 0, 955 0, 0, 956 0, 0, -2, -2, 957 0, -2, -2, 0, 958 -2, -2, 0, 0, 959 -2, 0, 0, -2, 960 -2, 0, 0, -2, 961 -2, -2, 0, 0, 962 0, -2, -2, 0, 963 0, 0, -2, -2, 964 0, 0, -2, -2, 965 -2, 0, 0, -2, 966 -2, -2, 0, 0, 967 0, -2, -2, 0, 968 0, 0, 0, 0, -4, 0, 969 0, 0, -1, 0, -4, 0, 970 0, 0, -1, 0, 0, 0, 971 0, 0, 0, 0, 0, 0, 972 -4, 0, 0, 0, -4, 0, 973 -4, 0, 0, 0, 0, 0, 974 -4, 0, -1, 0, 0, 0, 975 -4, 0, -1, 0, -4, 0}; 976 /* Add 3 quads inside every triangular prism, making 4 new prisms. */ 977 static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM}; 978 static PetscInt tripS[] = {3, 4, 6, 8}; 979 static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0, 980 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0, 981 DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0, 982 DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1, 983 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0, 984 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, 985 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, 986 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 987 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, 988 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0, 989 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, 990 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, 991 DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, 992 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, 993 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, 994 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, 995 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, 996 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, 997 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, 998 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, 999 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}; 1000 static PetscInt tripO[] = {0, 0, 1001 0, 0, 1002 0, 0, 1003 0, -2, -2, 1004 -2, 0, -2, 1005 -2, -2, 0, 1006 0, 0, 0, 1007 -2, 0, -2, -2, 1008 -2, 0, -2, -2, 1009 -2, 0, -2, -2, 1010 0, -2, -2, 0, 1011 0, -2, -2, 0, 1012 0, -2, -2, 0, 1013 0, 0, 0, -1, 0, 1014 0, 0, 0, 0, -1, 1015 0, 0, -1, 0, 0, 1016 2, 0, 0, 0, 0, 1017 -3, 0, 0, -1, 0, 1018 -3, 0, 0, 0, -1, 1019 -3, 0, -1, 0, 0, 1020 -3, 0, 0, 0, 0}; 1021 /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms. 1022 2 1023 |\ 1024 | \ 1025 | \ 1026 0---1 1027 1028 2 1029 1030 0 1 1031 1032 2 1033 |\ 1034 | \ 1035 | \ 1036 0---1 1037 */ 1038 static DMPolytopeType ttripT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR}; 1039 static PetscInt ttripS[] = {3, 4}; 1040 static PetscInt ttripC[] = {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, 1041 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, 1042 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, 1043 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, 1044 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, 1045 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, 1046 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}; 1047 static PetscInt ttripO[] = {0, 0, 0, 0, 1048 0, 0, 0, 0, 1049 0, 0, 0, 0, 1050 0, 0, 0, -1, 0, 1051 0, 0, 0, 0, -1, 1052 0, 0, -1, 0, 0, 1053 0, 0, 0, 0, 0}; 1054 /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */ 1055 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1056 static PetscInt tquadpS[] = {1, 4, 4}; 1057 static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1058 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, 1059 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, 1060 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, 1061 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, 1062 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, 1063 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, 1064 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, 1065 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}; 1066 static PetscInt tquadpO[] = {0, 0, 1067 0, 0, 0, 0, 1068 0, 0, 0, 0, 1069 0, 0, 0, 0, 1070 0, 0, 0, 0, 1071 0, 0, 0, 0, -1, 0, 1072 0, 0, 0, 0, 0, -1, 1073 0, 0, -1, 0, 0, 0, 1074 0, 0, 0, -1, 0, 0}; 1075 1076 PetscFunctionBegin; 1077 switch (source) { 1078 case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break; 1079 case DM_POLYTOPE_SEGMENT: *Nt = 2; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break; 1080 case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break; 1081 case DM_POLYTOPE_TRIANGLE: *Nt = 2; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 1082 case DM_POLYTOPE_QUADRILATERAL: *Nt = 3; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break; 1083 case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break; 1084 case DM_POLYTOPE_TETRAHEDRON: *Nt = 3; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 1085 case DM_POLYTOPE_HEXAHEDRON: *Nt = 4; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break; 1086 case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; 1087 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 2; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break; 1088 case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break; 1089 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1090 } 1091 PetscFunctionReturn(0); 1092 } 1093 1094 static PetscErrorCode DMPlexCellRefinerRefine_ToBox(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1095 { 1096 PetscErrorCode ierr; 1097 /* Change tensor edges to segments */ 1098 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_SEGMENT}; 1099 static PetscInt tedgeS[] = {1}; 1100 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 1101 static PetscInt tedgeO[] = { 0, 0}; 1102 /* Add 1 vertex, 3 edges inside every triangle, making 3 new quadrilaterals. 1103 2 1104 |\ 1105 | \ 1106 | \ 1107 | \ 1108 0 1 1109 | \ 1110 | \ 1111 2 1 1112 |\ / \ 1113 | 2 1 \ 1114 | \ / \ 1115 1 | 0 1116 | 0 \ 1117 | | \ 1118 | | \ 1119 0-0-0-----1-----1 1120 */ 1121 static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL}; 1122 static PetscInt triS[] = {1, 3, 3}; 1123 static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1124 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 1125 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 1126 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, 1127 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, 1128 DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0}; 1129 static PetscInt triO[] = {0, 0, 1130 0, 0, 1131 0, 0, 1132 0, 0, -2, 0, 1133 0, 0, 0, -2, 1134 0, -2, 0, 0}; 1135 /* Add 1 edge inside every tensor quad, making 2 new quadrilaterals 1136 2----2----1----3----3 1137 | | | 1138 | | | 1139 | | | 1140 4 A 6 B 5 1141 | | | 1142 | | | 1143 | | | 1144 0----0----0----1----1 1145 */ 1146 static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL}; 1147 static PetscInt tquadS[] = {1, 2}; 1148 static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1149 /* TODO Fix these */ 1150 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1151 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0}; 1152 static PetscInt tquadO[] = {0, 0, 1153 0, 0, -2, -2, 1154 0, 0, -2, -2}; 1155 /* Add 6 triangles inside every cell, making 4 new hexs 1156 TODO: Need different SubcellMap(). Need to make a struct with the function pointers in it 1157 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 1158 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] 1159 called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are 1160 [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3] 1161 We make a new hex in each corner 1162 [v0, (e0, 0), (f0, 0), (e2, 0), (e3, 0), (f2, 0), (c0, 0), (f1, 0)] 1163 [v1, (e4, 0), (f3, 0), (e1, 0), (e0, 0), (f0, 0), (c0, 0), (f1, 0)] 1164 [v2, (e1, 0), (f3, 0), (e5, 0), (e2, 0), (f2, 0), (c0, 0), (f0, 0)] 1165 [v3, (e4, 0), (f1, 0), (e3, 0), (e5, 0), (f2, 0), (c0, 0), (f3, 0)] 1166 We create a new face for each edge 1167 [(e3, 0), (f2, 0), (c0, 0), (f1, 0)] 1168 [(f0, 0), (e0, 0), (f1, 0), (c0, 0)] 1169 [(e2, 0), (f0, 0), (c0, 0), (f2, 0)] 1170 [(f3, 0), (e4, 0), (f1, 0), (c0, 0)] 1171 [(e1, 0), (f3, 0), (c0, 0), (f0, 0)] 1172 [(e5, 0), (f3, 0), (c0, 0), (f2, 0)] 1173 I could write a program to generate these from the first hex by acting with the symmetry group to take one subcell into another. 1174 */ 1175 static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 1176 static PetscInt tetS[] = {1, 4, 6, 4}; 1177 static PetscInt tetC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1178 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 1179 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 1180 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 1181 DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, 1182 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, 1183 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1184 DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 3, 1185 DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, 1186 DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, 1187 DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, 1188 DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, 1189 DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, 1190 DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2}; 1191 static PetscInt tetO[] = {0, 0, 1192 0, 0, 1193 0, 0, 1194 0, 0, 1195 0, 0, -2, -2, 1196 -2, 0, 0, -2, 1197 0, 0, -2, -2, 1198 -2, 0, 0, -2, 1199 0, 0, -2, -2, 1200 0, 0, -2, -2, 1201 0, 0, 0, 0, 0, 0, 1202 1, -1, 1, 0, 0, 3, 1203 0, -4, 1, -1, 0, 3, 1204 1, -4, 3, -2, -4, 3}; 1205 /* Add 3 quads inside every triangular prism, making 4 new prisms. */ 1206 static DMPolytopeType tripT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 1207 static PetscInt tripS[] = {1, 5, 9, 6}; 1208 static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1209 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 1210 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 1211 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 1212 DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0, 1213 DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, 1214 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, 1215 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, 1216 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1217 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, 1218 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0, 1219 DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, 1220 DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, 1221 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, 1222 DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, 1223 DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 3, 1224 DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, 1225 DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2, 1226 DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 0, 6, 1227 DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3}; 1228 static PetscInt tripO[] = {0, 0, 1229 0, 0, 1230 0, 0, 1231 0, 0, 1232 0, 0, 1233 0, 0, -2, -2, 1234 -2, 0, 0, -2, 1235 0, -2, -2, 0, 1236 0, 0, -2, -2, 1237 0, 0, -2, -2, 1238 0, 0, -2, -2, 1239 0, -2, -2, 0, 1240 0, -2, -2, 0, 1241 0, -2, -2, 0, 1242 0, 0, 0, -1, 0, 1, 1243 0, 0, 0, 0, 0, -4, 1244 0, 0, 0, 0, -1, 1, 1245 -4, 0, 0, -1, 0, 1, 1246 -4, 0, 0, 0, 0, -4, 1247 -4, 0, 0, 0, -1, 1}; 1248 /* Add 3 tensor quads inside every tensor triangular prism, making 4 new tensor triangular prisms. 1249 2 1250 |\ 1251 | \ 1252 | \ 1253 0---1 1254 1255 2 1256 1257 0 1 1258 1259 2 1260 |\ 1261 | \ 1262 | \ 1263 0---1 1264 */ 1265 static DMPolytopeType ttripT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1266 static PetscInt ttripS[] = {1, 3, 3}; 1267 static PetscInt ttripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1268 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, 1269 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, 1270 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, 1271 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, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, 1272 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, 1273 DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0}; 1274 static PetscInt ttripO[] = {0, 0, 1275 0, 0, 0, 0, 1276 0, 0, 0, 0, 1277 0, 0, 0, 0, 1278 0, 0, 0, 0, -1, 0, 1279 0, 0, 0, 0, 0, -1, 1280 0, 0, 0, -1, 0, 0}; 1281 /* TODO Add 3 quads inside every tensor triangular prism, making 4 new triangular prisms. 1282 2 1283 |\ 1284 | \ 1285 | \ 1286 0---1 1287 1288 2 1289 1290 0 1 1291 1292 2 1293 |\ 1294 | \ 1295 | \ 1296 0---1 1297 */ 1298 static DMPolytopeType ctripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 1299 static PetscInt ctripS[] = {1, 3, 3}; 1300 static PetscInt ctripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1301 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1302 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, 1303 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0, 1304 DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, 1305 DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, 1306 DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0}; 1307 static PetscInt ctripO[] = {0, 0, 1308 0, 0, -2, -2, 1309 0, 0, -2, -2, 1310 0, 0, -2, -2, 1311 -4, 0, 0, -1, 0, 1, 1312 -4, 0, 0, 0, 0, -4, 1313 -4, 0, 0, 0, -1, 1}; 1314 /* Add 1 edge and 4 quads inside every tensor quad prism, making 4 new hexahedra. */ 1315 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1316 static PetscInt tquadpS[] = {1, 4, 4}; 1317 static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1318 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, 1319 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, 1320 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, 1321 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, 1322 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, 1323 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, 1324 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, 1325 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}; 1326 static PetscInt tquadpO[] = {0, 0, 1327 0, 0, 0, 0, 1328 0, 0, 0, 0, 1329 0, 0, 0, 0, 1330 0, 0, 0, 0, 1331 0, 0, 0, 0, -1, 0, 1332 0, 0, 0, 0, 0, -1, 1333 0, 0, -1, 0, 0, 0, 1334 0, 0, 0, -1, 0, 0}; 1335 PetscBool convertTensor = PETSC_TRUE; 1336 1337 PetscFunctionBeginHot; 1338 if (convertTensor) { 1339 switch (source) { 1340 case DM_POLYTOPE_POINT: 1341 case DM_POLYTOPE_SEGMENT: 1342 case DM_POLYTOPE_QUADRILATERAL: 1343 case DM_POLYTOPE_HEXAHEDRON: 1344 ierr = DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 1345 break; 1346 case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break; 1347 case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break; 1348 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ctripT; *size = ctripS; *cone = ctripC; *ornt = ctripO; break; 1349 case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break; 1350 case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 1351 case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 1352 case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; 1353 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1354 } 1355 } else { 1356 switch (source) { 1357 case DM_POLYTOPE_POINT: 1358 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1359 case DM_POLYTOPE_SEGMENT: 1360 case DM_POLYTOPE_QUADRILATERAL: 1361 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1362 case DM_POLYTOPE_HEXAHEDRON: 1363 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1364 ierr = DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 1365 break; 1366 case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 1367 case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 1368 case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; 1369 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break; 1370 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1371 } 1372 } 1373 PetscFunctionReturn(0); 1374 } 1375 1376 static PetscErrorCode DMPlexCellRefinerRefine_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1377 { 1378 PetscErrorCode ierr; 1379 1380 PetscFunctionBeginHot; 1381 switch (source) { 1382 case DM_POLYTOPE_POINT: 1383 case DM_POLYTOPE_SEGMENT: 1384 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1385 case DM_POLYTOPE_TRIANGLE: 1386 case DM_POLYTOPE_TETRAHEDRON: 1387 case DM_POLYTOPE_TRI_PRISM: 1388 case DM_POLYTOPE_TRI_PRISM_TENSOR: 1389 case DM_POLYTOPE_QUADRILATERAL: 1390 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1391 case DM_POLYTOPE_HEXAHEDRON: 1392 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1393 ierr = DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 1394 break; 1395 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1396 } 1397 PetscFunctionReturn(0); 1398 } 1399 1400 static PetscErrorCode CellRefinerCreateOffset_Internal(DMPlexCellRefiner cr, PetscInt ctOrder[], PetscInt ctStart[], PetscInt **offset) 1401 { 1402 PetscInt c, cN, *off; 1403 PetscErrorCode ierr; 1404 1405 PetscFunctionBegin; 1406 ierr = PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);CHKERRQ(ierr); 1407 for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) { 1408 const DMPolytopeType ct = (DMPolytopeType) c; 1409 for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) { 1410 const DMPolytopeType ctNew = (DMPolytopeType) cN; 1411 DMPolytopeType *rct; 1412 PetscInt *rsize, *cone, *ornt; 1413 PetscInt Nct, n, i; 1414 1415 if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[ct*DM_NUM_POLYTOPES+ctNew] = -1; break;} 1416 off[ct*DM_NUM_POLYTOPES+ctNew] = 0; 1417 for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) { 1418 const DMPolytopeType ict = (DMPolytopeType) ctOrder[i]; 1419 const DMPolytopeType ictn = (DMPolytopeType) ctOrder[i+1]; 1420 1421 ierr = DMPlexCellRefinerRefine(cr, ict, &Nct, &rct, &rsize, &cone, &ornt);CHKERRQ(ierr); 1422 if (ict == ct) { 1423 for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break; 1424 if (n == Nct) off[ct*DM_NUM_POLYTOPES+ctNew] = -1; 1425 break; 1426 } 1427 for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct*DM_NUM_POLYTOPES+ctNew] += (ctStart[ictn]-ctStart[ict]) * rsize[n]; 1428 } 1429 } 1430 } 1431 *offset = off; 1432 PetscFunctionReturn(0); 1433 } 1434 1435 static PetscErrorCode DMPlexCellRefinerSetStarts(DMPlexCellRefiner cr, const PetscInt ctStart[], const PetscInt ctStartNew[]) 1436 { 1437 const PetscInt ctSize = DM_NUM_POLYTOPES+1; 1438 PetscErrorCode ierr; 1439 1440 PetscFunctionBegin; 1441 if (cr->setupcalled) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_ARG_WRONGSTATE, "Must call this function before DMPlexCellRefinerSetUp()"); 1442 ierr = PetscCalloc2(ctSize, &cr->ctStart, ctSize, &cr->ctStartNew);CHKERRQ(ierr); 1443 ierr = PetscArraycpy(cr->ctStart, ctStart, ctSize);CHKERRQ(ierr); 1444 ierr = PetscArraycpy(cr->ctStartNew, ctStartNew, ctSize);CHKERRQ(ierr); 1445 PetscFunctionReturn(0); 1446 } 1447 1448 /* Construct cell type order since we must loop over cell types in depth order */ 1449 PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell, PetscInt *ctOrder[], PetscInt *ctOrderInv[]) 1450 { 1451 PetscInt *ctO, *ctOInv; 1452 PetscInt c, d, off = 0; 1453 PetscErrorCode ierr; 1454 1455 PetscFunctionBegin; 1456 ierr = PetscCalloc2(DM_NUM_POLYTOPES+1, &ctO, DM_NUM_POLYTOPES+1, &ctOInv);CHKERRQ(ierr); 1457 for (d = 3; d >= DMPolytopeTypeGetDim(ctCell); --d) { 1458 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1459 if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue; 1460 ctO[off++] = c; 1461 } 1462 } 1463 if (DMPolytopeTypeGetDim(ctCell) != 0) { 1464 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1465 if (DMPolytopeTypeGetDim((DMPolytopeType) c) != 0) continue; 1466 ctO[off++] = c; 1467 } 1468 } 1469 for (d = DMPolytopeTypeGetDim(ctCell)-1; d > 0; --d) { 1470 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1471 if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue; 1472 ctO[off++] = c; 1473 } 1474 } 1475 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1476 if (DMPolytopeTypeGetDim((DMPolytopeType) c) >= 0) continue; 1477 ctO[off++] = c; 1478 } 1479 if (off != DM_NUM_POLYTOPES+1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %D for cell type order", off); 1480 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1481 ctOInv[ctO[c]] = c; 1482 } 1483 *ctOrder = ctO; 1484 *ctOrderInv = ctOInv; 1485 PetscFunctionReturn(0); 1486 } 1487 1488 PetscErrorCode DMPlexCellRefinerSetUp(DMPlexCellRefiner cr) 1489 { 1490 DM dm = cr->dm; 1491 DMPolytopeType ctCell; 1492 PetscInt pStart, pEnd, p, c; 1493 PetscErrorCode ierr; 1494 1495 PetscFunctionBegin; 1496 PetscValidHeader(cr, 1); 1497 if (cr->setupcalled) PetscFunctionReturn(0); 1498 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1499 if (pEnd > pStart) {ierr = DMPlexGetCellType(dm, 0, &ctCell);CHKERRQ(ierr);} 1500 else { 1501 PetscInt dim; 1502 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1503 switch (dim) { 1504 case 0: ctCell = DM_POLYTOPE_POINT;break; 1505 case 1: ctCell = DM_POLYTOPE_SEGMENT;break; 1506 case 2: ctCell = DM_POLYTOPE_TRIANGLE;break; 1507 case 3: ctCell = DM_POLYTOPE_TETRAHEDRON;break; 1508 default: ctCell = DM_POLYTOPE_TETRAHEDRON; 1509 } 1510 } 1511 ierr = DMPlexCreateCellTypeOrder_Internal(ctCell, &cr->ctOrder, &cr->ctOrderInv);CHKERRQ(ierr); 1512 /* Construct sizes and offsets for each cell type */ 1513 if (!cr->ctStart) { 1514 PetscInt *ctS, *ctSN, *ctC, *ctCN; 1515 1516 ierr = PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);CHKERRQ(ierr); 1517 ierr = PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);CHKERRQ(ierr); 1518 for (p = pStart; p < pEnd; ++p) { 1519 DMPolytopeType ct; 1520 DMPolytopeType *rct; 1521 PetscInt *rsize, *cone, *ornt; 1522 PetscInt Nct, n; 1523 1524 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1525 if ((PetscInt) ct < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %D", p); 1526 ++ctC[ct]; 1527 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);CHKERRQ(ierr); 1528 for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n]; 1529 } 1530 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 1531 const PetscInt ct = cr->ctOrder[c]; 1532 const PetscInt ctn = cr->ctOrder[c+1]; 1533 1534 ctS[ctn] = ctS[ct] + ctC[ct]; 1535 ctSN[ctn] = ctSN[ct] + ctCN[ct]; 1536 } 1537 ierr = PetscFree2(ctC, ctCN);CHKERRQ(ierr); 1538 cr->ctStart = ctS; 1539 cr->ctStartNew = ctSN; 1540 } 1541 ierr = CellRefinerCreateOffset_Internal(cr, cr->ctOrder, cr->ctStart, &cr->offset);CHKERRQ(ierr); 1542 cr->setupcalled = PETSC_TRUE; 1543 PetscFunctionReturn(0); 1544 } 1545 1546 static PetscErrorCode DMPlexCellRefinerView_Ascii(DMPlexCellRefiner cr, PetscViewer v) 1547 { 1548 PetscErrorCode ierr; 1549 1550 PetscFunctionBegin; 1551 ierr = PetscViewerASCIIPrintf(v, "Cell Refiner: %s\n", DMPlexCellRefinerTypes[cr->type]);CHKERRQ(ierr); 1552 PetscFunctionReturn(0); 1553 } 1554 1555 /* 1556 DMPlexCellRefinerView - Views a DMPlexCellRefiner object 1557 1558 Collective on cr 1559 1560 Input Parameters: 1561 + cr - The DMPlexCellRefiner object 1562 - viewer - The PetscViewer object 1563 1564 Level: beginner 1565 1566 .seealso: DMPlexCellRefinerCreate() 1567 */ 1568 static PetscErrorCode DMPlexCellRefinerView(DMPlexCellRefiner cr, PetscViewer viewer) 1569 { 1570 PetscBool iascii; 1571 PetscErrorCode ierr; 1572 1573 PetscFunctionBegin; 1574 PetscValidHeader(cr, 1); 1575 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 1576 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) cr), &viewer);CHKERRQ(ierr);} 1577 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 1578 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1579 if (iascii) {ierr = DMPlexCellRefinerView_Ascii(cr, viewer);CHKERRQ(ierr);} 1580 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1581 PetscFunctionReturn(0); 1582 } 1583 1584 PetscErrorCode DMPlexCellRefinerDestroy(DMPlexCellRefiner *cr) 1585 { 1586 PetscInt c; 1587 PetscErrorCode ierr; 1588 1589 PetscFunctionBegin; 1590 if (!*cr) PetscFunctionReturn(0); 1591 PetscValidHeaderSpecific(*cr, DM_CLASSID, 1); 1592 ierr = PetscObjectDereference((PetscObject) (*cr)->dm);CHKERRQ(ierr); 1593 ierr = PetscFree2((*cr)->ctOrder, (*cr)->ctOrderInv);CHKERRQ(ierr); 1594 ierr = PetscFree2((*cr)->ctStart, (*cr)->ctStartNew);CHKERRQ(ierr); 1595 ierr = PetscFree((*cr)->offset);CHKERRQ(ierr); 1596 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 1597 ierr = PetscFEDestroy(&(*cr)->coordFE[c]);CHKERRQ(ierr); 1598 ierr = PetscFEGeomDestroy(&(*cr)->refGeom[c]);CHKERRQ(ierr); 1599 } 1600 ierr = PetscFree2((*cr)->coordFE, (*cr)->refGeom);CHKERRQ(ierr); 1601 ierr = PetscHeaderDestroy(cr);CHKERRQ(ierr); 1602 PetscFunctionReturn(0); 1603 } 1604 1605 PetscErrorCode DMPlexCellRefinerCreate(DM dm, DMPlexCellRefiner *cr) 1606 { 1607 DMPlexCellRefiner tmp; 1608 PetscErrorCode ierr; 1609 1610 PetscFunctionBegin; 1611 PetscValidPointer(cr, 2); 1612 *cr = NULL; 1613 ierr = PetscHeaderCreate(tmp, DM_CLASSID, "DMPlexCellRefiner", "Cell Refiner", "DMPlexCellRefiner", PETSC_COMM_SELF, DMPlexCellRefinerDestroy, DMPlexCellRefinerView);CHKERRQ(ierr); 1614 tmp->setupcalled = PETSC_FALSE; 1615 1616 tmp->dm = dm; 1617 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1618 ierr = DMPlexGetCellRefinerType(dm, &tmp->type);CHKERRQ(ierr); 1619 switch (tmp->type) { 1620 case DM_REFINER_REGULAR: 1621 tmp->ops->refine = DMPlexCellRefinerRefine_Regular; 1622 tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_Regular; 1623 tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_Regular; 1624 tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_Regular; 1625 tmp->ops->getaffinetransforms = DMPlexCellRefinerGetAffineTransforms_Regular; 1626 tmp->ops->getaffinefacetransforms = DMPlexCellRefinerGetAffineFaceTransforms_Regular; 1627 break; 1628 case DM_REFINER_TO_BOX: 1629 tmp->ops->refine = DMPlexCellRefinerRefine_ToBox; 1630 tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToBox; 1631 tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_ToBox; 1632 tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_ToBox; 1633 break; 1634 case DM_REFINER_TO_SIMPLEX: 1635 tmp->ops->refine = DMPlexCellRefinerRefine_ToSimplex; 1636 tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToSimplex; 1637 break; 1638 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid cell refiner type %s", DMPlexCellRefinerTypes[tmp->type]); 1639 } 1640 ierr = PetscCalloc2(DM_NUM_POLYTOPES, &tmp->coordFE, DM_NUM_POLYTOPES, &tmp->refGeom);CHKERRQ(ierr); 1641 *cr = tmp; 1642 PetscFunctionReturn(0); 1643 } 1644 1645 /*@ 1646 DMPlexCellRefinerGetAffineTransforms - Gets the affine map from the reference cell to each subcell 1647 1648 Input Parameters: 1649 + cr - The DMPlexCellRefiner object 1650 - ct - The cell type 1651 1652 Output Parameters: 1653 + Nc - The number of subcells produced from this cell type 1654 . v0 - The translation of the first vertex for each subcell 1655 . J - The Jacobian for each subcell (map from reference cell to subcell) 1656 - invJ - The inverse Jacobian for each subcell 1657 1658 Level: developer 1659 1660 .seealso: DMPlexCellRefinerGetAffineFaceTransforms(), Create() 1661 @*/ 1662 PetscErrorCode DMPlexCellRefinerGetAffineTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[]) 1663 { 1664 PetscErrorCode ierr; 1665 1666 PetscFunctionBegin; 1667 if (!cr->ops->getaffinetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine transforms from this refiner"); 1668 ierr = (*cr->ops->getaffinetransforms)(cr, ct, Nc, v0, J, invJ);CHKERRQ(ierr); 1669 PetscFunctionReturn(0); 1670 } 1671 1672 /*@ 1673 DMPlexCellRefinerGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell 1674 1675 Input Parameters: 1676 + cr - The DMPlexCellRefiner object 1677 - ct - The cell type 1678 1679 Output Parameters: 1680 + Nf - The number of faces for this cell type 1681 . v0 - The translation of the first vertex for each face 1682 . J - The Jacobian for each face (map from original cell to subcell) 1683 . invJ - The inverse Jacobian for each face 1684 - detJ - The determinant of the Jacobian for each face 1685 1686 Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse. 1687 1688 Level: developer 1689 1690 .seealso: DMPlexCellRefinerGetAffineTransforms(), Create() 1691 @*/ 1692 PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[]) 1693 { 1694 PetscErrorCode ierr; 1695 1696 PetscFunctionBegin; 1697 if (!cr->ops->getaffinefacetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine face transforms from this refiner"); 1698 ierr = (*cr->ops->getaffinefacetransforms)(cr, ct, Nf, v0, J, invJ, detJ);CHKERRQ(ierr); 1699 PetscFunctionReturn(0); 1700 } 1701 1702 /* Numbering regularly refined meshes 1703 1704 We want the numbering of the new mesh to respect the same depth stratification as the old mesh. We first compute 1705 the number of new points at each depth. This means that offsets for each depth can be computed, making no assumptions 1706 about the order of different cell types. 1707 1708 However, when we want to order different depth strata, it will be very useful to make assumptions about contiguous 1709 numbering of different cell types, especially if we want to compute new numberings without communication. Therefore, we 1710 will require that cells are numbering contiguously for each cell type, and that those blocks come in the same order as 1711 the cell type enumeration within a given depth stratum. 1712 1713 Thus, at each depth, each cell type will add a certain number of points at that depth. To get the new point number, we 1714 start at the new depth offset, run through all prior cell types incrementing by the total addition from that type, then 1715 offset by the old cell type number and replica number for the insertion. 1716 */ 1717 PetscErrorCode DMPlexCellRefinerGetNewPoint(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew) 1718 { 1719 DMPolytopeType *rct; 1720 PetscInt *rsize, *cone, *ornt; 1721 PetscInt Nct, n; 1722 PetscInt off = cr->offset[ct*DM_NUM_POLYTOPES+ctNew]; 1723 PetscInt ctS = cr->ctStart[ct], ctE = cr->ctStart[cr->ctOrder[cr->ctOrderInv[ct]+1]]; 1724 PetscInt ctSN = cr->ctStartNew[ctNew], ctEN = cr->ctStartNew[cr->ctOrder[cr->ctOrderInv[ctNew]+1]]; 1725 PetscInt newp = ctSN; 1726 PetscErrorCode ierr; 1727 1728 PetscFunctionBeginHot; 1729 if ((p < ctS) || (p >= ctE)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D is not a %s [%D, %D)", p, DMPolytopeTypes[ct], ctS, ctE); 1730 if (off < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s does not produce type %s", DMPolytopeTypes[ct], DMPolytopeTypes[ctNew]); 1731 1732 newp += off; 1733 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);CHKERRQ(ierr); 1734 for (n = 0; n < Nct; ++n) { 1735 if (rct[n] == ctNew) { 1736 if (rsize[n] && r >= rsize[n]) 1737 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %D should be in [0, %D) for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]); 1738 newp += (p - ctS) * rsize[n] + r; 1739 break; 1740 } 1741 } 1742 1743 if ((newp < ctSN) || (newp >= ctEN)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "New point %D is not a %s [%D, %D)", newp, DMPolytopeTypes[ctNew], ctSN, ctEN); 1744 *pNew = newp; 1745 PetscFunctionReturn(0); 1746 } 1747 1748 static PetscErrorCode DMPlexCellRefinerSetConeSizes(DMPlexCellRefiner cr, DM rdm) 1749 { 1750 DM dm = cr->dm; 1751 PetscInt pStart, pEnd, p, pNew; 1752 PetscErrorCode ierr; 1753 1754 PetscFunctionBegin; 1755 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1756 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 1757 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1758 for (p = pStart; p < pEnd; ++p) { 1759 DMPolytopeType ct; 1760 DMPolytopeType *rct; 1761 PetscInt *rsize, *rcone, *rornt; 1762 PetscInt Nct, n, r; 1763 1764 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1765 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1766 for (n = 0; n < Nct; ++n) { 1767 for (r = 0; r < rsize[n]; ++r) { 1768 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 1769 ierr = DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));CHKERRQ(ierr); 1770 ierr = DMPlexSetCellType(rdm, pNew, rct[n]);CHKERRQ(ierr); 1771 } 1772 } 1773 } 1774 { 1775 DMLabel ctLabel; 1776 DM_Plex *plex = (DM_Plex *) rdm->data; 1777 1778 ierr = DMPlexGetCellTypeLabel(rdm, &ctLabel);CHKERRQ(ierr); 1779 ierr = PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);CHKERRQ(ierr); 1780 } 1781 PetscFunctionReturn(0); 1782 } 1783 1784 static PetscErrorCode DMPlexCellRefinerSetCones(DMPlexCellRefiner cr, DM rdm) 1785 { 1786 DM dm = cr->dm; 1787 DMPolytopeType ct; 1788 PetscInt *coneNew, *orntNew; 1789 PetscInt maxConeSize = 0, pStart, pEnd, p, pNew; 1790 PetscErrorCode ierr; 1791 1792 PetscFunctionBegin; 1793 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p)); 1794 ierr = PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew);CHKERRQ(ierr); 1795 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1796 for (p = pStart; p < pEnd; ++p) { 1797 const PetscInt *cone, *ornt; 1798 PetscInt coff, ooff, c; 1799 DMPolytopeType *rct; 1800 PetscInt *rsize, *rcone, *rornt; 1801 PetscInt Nct, n, r; 1802 1803 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1804 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 1805 ierr = DMPlexGetConeOrientation(dm, p, &ornt);CHKERRQ(ierr); 1806 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1807 for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) { 1808 const DMPolytopeType ctNew = rct[n]; 1809 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1810 1811 for (r = 0; r < rsize[n]; ++r) { 1812 /* pNew is a subcell produced by subdividing p */ 1813 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 1814 for (c = 0; c < csizeNew; ++c) { 1815 PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */ 1816 PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */ 1817 PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */ 1818 DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */ 1819 const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */ 1820 PetscInt pr = -1; /* Replica number of pp that produces new cone point */ 1821 const DMPolytopeType ft = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */ 1822 const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1823 PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */ 1824 PetscInt lc; 1825 1826 /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */ 1827 for (lc = 0; lc < fn; ++lc) { 1828 const PetscInt *ppornt; 1829 PetscInt pcp; 1830 1831 ierr = DMPolytopeMapCell(pct, po, rcone[coff++], &pcp);CHKERRQ(ierr); 1832 ppp = pp; 1833 pp = pcone[pcp]; 1834 ierr = DMPlexGetCellType(dm, pp, &pct);CHKERRQ(ierr); 1835 ierr = DMPlexGetCone(dm, pp, &pcone);CHKERRQ(ierr); 1836 ierr = DMPlexGetConeOrientation(dm, ppp, &ppornt);CHKERRQ(ierr); 1837 po = ppornt[pcp]; 1838 } 1839 pr = rcone[coff++]; 1840 /* Orientation po of pp maps (pr, fo) -> (pr', fo') */ 1841 ierr = DMPlexCellRefinerMapSubcells(cr, pct, po, ft, pr, fo, &pr, &fo);CHKERRQ(ierr); 1842 ierr = DMPlexCellRefinerGetNewPoint(cr, pct, ft, pp, pr, &coneNew[c]);CHKERRQ(ierr); 1843 orntNew[c] = fo; 1844 } 1845 ierr = DMPlexSetCone(rdm, pNew, coneNew);CHKERRQ(ierr); 1846 ierr = DMPlexSetConeOrientation(rdm, pNew, orntNew);CHKERRQ(ierr); 1847 } 1848 } 1849 } 1850 ierr = PetscFree2(coneNew, orntNew);CHKERRQ(ierr); 1851 ierr = DMPlexSymmetrize(rdm);CHKERRQ(ierr); 1852 ierr = DMPlexStratify(rdm);CHKERRQ(ierr); 1853 PetscFunctionReturn(0); 1854 } 1855 1856 static PetscErrorCode DMPlexCellRefinerGetCoordinateFE(DMPlexCellRefiner cr, DMPolytopeType ct, PetscFE *fe) 1857 { 1858 PetscErrorCode ierr; 1859 1860 PetscFunctionBegin; 1861 if (!cr->coordFE[ct]) { 1862 PetscInt dim, cdim; 1863 PetscBool isSimplex; 1864 1865 switch (ct) { 1866 case DM_POLYTOPE_SEGMENT: dim = 1; isSimplex = PETSC_TRUE; break; 1867 case DM_POLYTOPE_TRIANGLE: dim = 2; isSimplex = PETSC_TRUE; break; 1868 case DM_POLYTOPE_QUADRILATERAL: dim = 2; isSimplex = PETSC_FALSE; break; 1869 case DM_POLYTOPE_TETRAHEDRON: dim = 3; isSimplex = PETSC_TRUE; break; 1870 case DM_POLYTOPE_HEXAHEDRON: dim = 3; isSimplex = PETSC_FALSE; break; 1871 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No coordinate FE for cell type %s", DMPolytopeTypes[ct]); 1872 } 1873 ierr = DMGetCoordinateDim(cr->dm, &cdim);CHKERRQ(ierr); 1874 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dim, cdim, isSimplex, 1, PETSC_DETERMINE, &cr->coordFE[ct]);CHKERRQ(ierr); 1875 { 1876 PetscDualSpace dsp; 1877 PetscQuadrature quad; 1878 DM K; 1879 PetscFEGeom *cg; 1880 PetscReal *Xq, *xq, *wq; 1881 PetscInt Nq, q; 1882 1883 ierr = DMPlexCellRefinerGetCellVertices(cr, ct, &Nq, &Xq);CHKERRQ(ierr); 1884 ierr = PetscMalloc1(Nq*cdim, &xq);CHKERRQ(ierr); 1885 for (q = 0; q < Nq*cdim; ++q) xq[q] = Xq[q]; 1886 ierr = PetscMalloc1(Nq, &wq);CHKERRQ(ierr); 1887 for (q = 0; q < Nq; ++q) wq[q] = 1.0; 1888 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &quad);CHKERRQ(ierr); 1889 ierr = PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);CHKERRQ(ierr); 1890 ierr = PetscFESetQuadrature(cr->coordFE[ct], quad);CHKERRQ(ierr); 1891 1892 ierr = PetscFEGetDualSpace(cr->coordFE[ct], &dsp);CHKERRQ(ierr); 1893 ierr = PetscDualSpaceGetDM(dsp, &K);CHKERRQ(ierr); 1894 ierr = PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &cr->refGeom[ct]);CHKERRQ(ierr); 1895 cg = cr->refGeom[ct]; 1896 ierr = DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);CHKERRQ(ierr); 1897 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 1898 } 1899 } 1900 *fe = cr->coordFE[ct]; 1901 PetscFunctionReturn(0); 1902 } 1903 1904 /* 1905 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct. 1906 1907 Not collective 1908 1909 Input Parameters: 1910 + cr - The DMPlexCellRefiner 1911 . ct - The type of the parent cell 1912 . rct - The type of the produced cell 1913 . r - The index of the produced cell 1914 - x - The localized coordinates for the parent cell 1915 1916 Output Parameter: 1917 . xr - The localized coordinates for the produced cell 1918 1919 Level: developer 1920 1921 .seealso: DMPlexCellRefinerSetCoordinates() 1922 */ 1923 static PetscErrorCode DMPlexCellRefinerMapLocalizedCoordinates(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[]) 1924 { 1925 PetscFE fe = NULL; 1926 PetscInt cdim, Nv, v, *subcellV; 1927 PetscErrorCode ierr; 1928 1929 PetscFunctionBegin; 1930 ierr = DMPlexCellRefinerGetCoordinateFE(cr, ct, &fe);CHKERRQ(ierr); 1931 ierr = DMPlexCellRefinerGetSubcellVertices(cr, ct, rct, r, &Nv, &subcellV);CHKERRQ(ierr); 1932 ierr = PetscFEGetNumComponents(fe, &cdim);CHKERRQ(ierr); 1933 for (v = 0; v < Nv; ++v) { 1934 ierr = PetscFEInterpolate_Static(fe, x, cr->refGeom[ct], subcellV[v], &xr[v*cdim]);CHKERRQ(ierr); 1935 } 1936 PetscFunctionReturn(0); 1937 } 1938 1939 static PetscErrorCode DMPlexCellRefinerSetCoordinates(DMPlexCellRefiner cr, DM rdm) 1940 { 1941 DM dm = cr->dm, cdm; 1942 PetscSection coordSection, coordSectionNew; 1943 Vec coordsLocal, coordsLocalNew; 1944 const PetscScalar *coords; 1945 PetscScalar *coordsNew; 1946 const DMBoundaryType *bd; 1947 const PetscReal *maxCell, *L; 1948 PetscBool isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE; 1949 PetscInt dE, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd; 1950 PetscErrorCode ierr; 1951 1952 PetscFunctionBegin; 1953 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1954 ierr = DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);CHKERRQ(ierr); 1955 /* Determine if we need to localize coordinates when generating them */ 1956 if (isperiodic) { 1957 localizeVertices = PETSC_TRUE; 1958 if (!maxCell) { 1959 PetscBool localized; 1960 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 1961 if (!localized) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized"); 1962 localizeCells = PETSC_TRUE; 1963 } 1964 } 1965 1966 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1967 ierr = PetscSectionGetFieldComponents(coordSection, 0, &dE);CHKERRQ(ierr); 1968 if (maxCell) { 1969 PetscReal maxCellNew[3]; 1970 1971 for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0; 1972 ierr = DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);CHKERRQ(ierr); 1973 } else { 1974 ierr = DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);CHKERRQ(ierr); 1975 } 1976 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);CHKERRQ(ierr); 1977 ierr = PetscSectionSetNumFields(coordSectionNew, 1);CHKERRQ(ierr); 1978 ierr = PetscSectionSetFieldComponents(coordSectionNew, 0, dE);CHKERRQ(ierr); 1979 ierr = DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);CHKERRQ(ierr); 1980 if (localizeCells) {ierr = PetscSectionSetChart(coordSectionNew, 0, vEndNew);CHKERRQ(ierr);} 1981 else {ierr = PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);CHKERRQ(ierr);} 1982 1983 /* Localization should be inherited */ 1984 /* Stefano calculates parent cells for each new cell for localization */ 1985 /* Localized cells need coordinates of closure */ 1986 for (v = vStartNew; v < vEndNew; ++v) { 1987 ierr = PetscSectionSetDof(coordSectionNew, v, dE);CHKERRQ(ierr); 1988 ierr = PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);CHKERRQ(ierr); 1989 } 1990 if (localizeCells) { 1991 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1992 for (c = cStart; c < cEnd; ++c) { 1993 PetscInt dof; 1994 1995 ierr = PetscSectionGetDof(coordSection, c, &dof); CHKERRQ(ierr); 1996 if (dof) { 1997 DMPolytopeType ct; 1998 DMPolytopeType *rct; 1999 PetscInt *rsize, *rcone, *rornt; 2000 PetscInt dim, cNew, Nct, n, r; 2001 2002 ierr = DMPlexGetCellType(dm, c, &ct);CHKERRQ(ierr); 2003 dim = DMPolytopeTypeGetDim(ct); 2004 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2005 /* This allows for different cell types */ 2006 for (n = 0; n < Nct; ++n) { 2007 if (dim != DMPolytopeTypeGetDim(rct[n])) continue; 2008 for (r = 0; r < rsize[n]; ++r) { 2009 PetscInt *closure = NULL; 2010 PetscInt clSize, cl, Nv = 0; 2011 2012 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], c, r, &cNew);CHKERRQ(ierr); 2013 ierr = DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 2014 for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;} 2015 ierr = DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 2016 ierr = PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);CHKERRQ(ierr); 2017 ierr = PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);CHKERRQ(ierr); 2018 } 2019 } 2020 } 2021 } 2022 } 2023 ierr = PetscSectionSetUp(coordSectionNew);CHKERRQ(ierr); 2024 ierr = DMViewFromOptions(dm, NULL, "-coarse_dm_view");CHKERRQ(ierr); 2025 ierr = DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);CHKERRQ(ierr); 2026 { 2027 VecType vtype; 2028 PetscInt coordSizeNew, bs; 2029 const char *name; 2030 2031 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 2032 ierr = VecCreate(PETSC_COMM_SELF, &coordsLocalNew);CHKERRQ(ierr); 2033 ierr = PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);CHKERRQ(ierr); 2034 ierr = VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);CHKERRQ(ierr); 2035 ierr = PetscObjectGetName((PetscObject) coordsLocal, &name);CHKERRQ(ierr); 2036 ierr = PetscObjectSetName((PetscObject) coordsLocalNew, name);CHKERRQ(ierr); 2037 ierr = VecGetBlockSize(coordsLocal, &bs);CHKERRQ(ierr); 2038 ierr = VecSetBlockSize(coordsLocalNew, bs);CHKERRQ(ierr); 2039 ierr = VecGetType(coordsLocal, &vtype);CHKERRQ(ierr); 2040 ierr = VecSetType(coordsLocalNew, vtype);CHKERRQ(ierr); 2041 } 2042 ierr = VecGetArrayRead(coordsLocal, &coords);CHKERRQ(ierr); 2043 ierr = VecGetArray(coordsLocalNew, &coordsNew);CHKERRQ(ierr); 2044 ierr = PetscSectionGetChart(coordSection, &ocStart, &ocEnd);CHKERRQ(ierr); 2045 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2046 /* First set coordinates for vertices*/ 2047 for (p = pStart; p < pEnd; ++p) { 2048 DMPolytopeType ct; 2049 DMPolytopeType *rct; 2050 PetscInt *rsize, *rcone, *rornt; 2051 PetscInt Nct, n, r; 2052 PetscBool hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE; 2053 2054 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2055 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2056 for (n = 0; n < Nct; ++n) { 2057 if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;} 2058 } 2059 if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) { 2060 PetscInt dof; 2061 ierr = PetscSectionGetDof(coordSection, p, &dof); CHKERRQ(ierr); 2062 if (dof) isLocalized = PETSC_TRUE; 2063 } 2064 if (hasVertex) { 2065 PetscScalar *pcoords = NULL; 2066 PetscScalar vcoords[3] = {0., 0., 0.}; 2067 PetscInt Nc, Nv, v, d; 2068 2069 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);CHKERRQ(ierr); 2070 if (ct == DM_POLYTOPE_POINT) { 2071 for (d = 0; d < dE; ++d) vcoords[d] = pcoords[d]; 2072 } else { 2073 if (localizeVertices) { 2074 PetscScalar anchor[3]; 2075 2076 for (d = 0; d < dE; ++d) anchor[d] = pcoords[d]; 2077 if (!isLocalized) { 2078 Nv = Nc/dE; 2079 for (v = 0; v < Nv; ++v) {ierr = DMLocalizeAddCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], vcoords);CHKERRQ(ierr);} 2080 } else { 2081 Nv = Nc/(2*dE); 2082 for (v = Nv; v < Nv*2; ++v) {ierr = DMLocalizeAddCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], vcoords);CHKERRQ(ierr);} 2083 } 2084 } else { 2085 Nv = Nc/dE; 2086 for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) vcoords[d] += pcoords[v*dE+d]; 2087 } 2088 for (d = 0; d < dE; ++d) vcoords[d] /= Nv; 2089 } 2090 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);CHKERRQ(ierr); 2091 for (n = 0; n < Nct; ++n) { 2092 if (rct[n] != DM_POLYTOPE_POINT) continue; 2093 if (rsize[n] > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Only support creating a single vertex in cell refinement, not %D", rsize[n]); 2094 for (r = 0; r < rsize[n]; ++r) { 2095 PetscInt vNew, off; 2096 2097 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &vNew);CHKERRQ(ierr); 2098 ierr = PetscSectionGetOffset(coordSectionNew, vNew, &off);CHKERRQ(ierr); 2099 ierr = DMPlexSnapToGeomModel(dm, p, vcoords, &coordsNew[off]);CHKERRQ(ierr); 2100 } 2101 } 2102 } 2103 } 2104 /* Then set coordinates for cells by localizing */ 2105 for (p = pStart; p < pEnd; ++p) { 2106 DMPolytopeType ct; 2107 DMPolytopeType *rct; 2108 PetscInt *rsize, *rcone, *rornt; 2109 PetscInt Nct, n, r; 2110 PetscBool isLocalized = PETSC_FALSE; 2111 2112 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2113 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2114 if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) { 2115 PetscInt dof; 2116 ierr = PetscSectionGetDof(coordSection, p, &dof); CHKERRQ(ierr); 2117 if (dof) isLocalized = PETSC_TRUE; 2118 } 2119 if (isLocalized) { 2120 const PetscScalar *pcoords; 2121 2122 ierr = DMPlexPointLocalRead(cdm, p, coords, &pcoords);CHKERRQ(ierr); 2123 for (n = 0; n < Nct; ++n) { 2124 const PetscInt Nr = rsize[n]; 2125 2126 if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue; 2127 for (r = 0; r < Nr; ++r) { 2128 PetscInt pNew, offNew; 2129 2130 /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that 2131 DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger 2132 cell to the ones it produces. */ 2133 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 2134 ierr = PetscSectionGetOffset(coordSectionNew, pNew, &offNew);CHKERRQ(ierr); 2135 ierr = DMPlexCellRefinerMapLocalizedCoordinates(cr, ct, rct[n], r, pcoords, &coordsNew[offNew]);CHKERRQ(ierr); 2136 } 2137 } 2138 } 2139 } 2140 ierr = VecRestoreArrayRead(coordsLocal, &coords);CHKERRQ(ierr); 2141 ierr = VecRestoreArray(coordsLocalNew, &coordsNew);CHKERRQ(ierr); 2142 ierr = DMSetCoordinatesLocal(rdm, coordsLocalNew);CHKERRQ(ierr); 2143 /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */ 2144 ierr = VecDestroy(&coordsLocalNew);CHKERRQ(ierr); 2145 ierr = PetscSectionDestroy(&coordSectionNew);CHKERRQ(ierr); 2146 if (!localizeCells) {ierr = DMLocalizeCoordinates(rdm);CHKERRQ(ierr);} 2147 PetscFunctionReturn(0); 2148 } 2149 2150 /*@ 2151 DMPlexCreateProcessSF - Create an SF which just has process connectivity 2152 2153 Collective on dm 2154 2155 Input Parameters: 2156 + dm - The DM 2157 - sfPoint - The PetscSF which encodes point connectivity 2158 2159 Output Parameters: 2160 + processRanks - A list of process neighbors, or NULL 2161 - sfProcess - An SF encoding the process connectivity, or NULL 2162 2163 Level: developer 2164 2165 .seealso: PetscSFCreate(), DMPlexCreateTwoSidedProcessSF() 2166 @*/ 2167 PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess) 2168 { 2169 PetscInt numRoots, numLeaves, l; 2170 const PetscInt *localPoints; 2171 const PetscSFNode *remotePoints; 2172 PetscInt *localPointsNew; 2173 PetscSFNode *remotePointsNew; 2174 PetscInt *ranks, *ranksNew; 2175 PetscMPIInt size; 2176 PetscErrorCode ierr; 2177 2178 PetscFunctionBegin; 2179 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2180 PetscValidHeaderSpecific(sfPoint, PETSCSF_CLASSID, 2); 2181 if (processRanks) {PetscValidPointer(processRanks, 3);} 2182 if (sfProcess) {PetscValidPointer(sfProcess, 4);} 2183 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2184 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2185 ierr = PetscMalloc1(numLeaves, &ranks);CHKERRQ(ierr); 2186 for (l = 0; l < numLeaves; ++l) { 2187 ranks[l] = remotePoints[l].rank; 2188 } 2189 ierr = PetscSortRemoveDupsInt(&numLeaves, ranks);CHKERRQ(ierr); 2190 ierr = PetscMalloc1(numLeaves, &ranksNew);CHKERRQ(ierr); 2191 ierr = PetscMalloc1(numLeaves, &localPointsNew);CHKERRQ(ierr); 2192 ierr = PetscMalloc1(numLeaves, &remotePointsNew);CHKERRQ(ierr); 2193 for (l = 0; l < numLeaves; ++l) { 2194 ranksNew[l] = ranks[l]; 2195 localPointsNew[l] = l; 2196 remotePointsNew[l].index = 0; 2197 remotePointsNew[l].rank = ranksNew[l]; 2198 } 2199 ierr = PetscFree(ranks);CHKERRQ(ierr); 2200 if (processRanks) {ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);CHKERRQ(ierr);} 2201 else {ierr = PetscFree(ranksNew);CHKERRQ(ierr);} 2202 if (sfProcess) { 2203 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);CHKERRQ(ierr); 2204 ierr = PetscObjectSetName((PetscObject) *sfProcess, "Process SF");CHKERRQ(ierr); 2205 ierr = PetscSFSetFromOptions(*sfProcess);CHKERRQ(ierr); 2206 ierr = PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2207 } 2208 PetscFunctionReturn(0); 2209 } 2210 2211 static PetscErrorCode DMPlexCellRefinerCreateSF(DMPlexCellRefiner cr, DM rdm) 2212 { 2213 DM dm = cr->dm; 2214 DMPlexCellRefiner *crRem; 2215 PetscSF sf, sfNew, sfProcess; 2216 IS processRanks; 2217 MPI_Datatype ctType; 2218 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 2219 const PetscInt *localPoints, *neighbors; 2220 const PetscSFNode *remotePoints; 2221 PetscInt *localPointsNew; 2222 PetscSFNode *remotePointsNew; 2223 PetscInt *ctStartRem, *ctStartNewRem; 2224 PetscInt ctSize = DM_NUM_POLYTOPES+1, numNeighbors, n, pStartNew, pEndNew, pNew, pNewRem; 2225 PetscErrorCode ierr; 2226 2227 PetscFunctionBegin; 2228 ierr = DMPlexGetChart(rdm, &pStartNew, &pEndNew);CHKERRQ(ierr); 2229 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2230 ierr = DMGetPointSF(rdm, &sfNew);CHKERRQ(ierr); 2231 /* Calculate size of new SF */ 2232 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2233 if (numRoots < 0) PetscFunctionReturn(0); 2234 for (l = 0; l < numLeaves; ++l) { 2235 const PetscInt p = localPoints[l]; 2236 DMPolytopeType ct; 2237 DMPolytopeType *rct; 2238 PetscInt *rsize, *rcone, *rornt; 2239 PetscInt Nct, n; 2240 2241 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2242 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2243 for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n]; 2244 } 2245 /* Communicate ctStart and cStartNew for each remote rank */ 2246 ierr = DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);CHKERRQ(ierr); 2247 ierr = ISGetLocalSize(processRanks, &numNeighbors);CHKERRQ(ierr); 2248 ierr = PetscMalloc2(ctSize*numNeighbors, &ctStartRem, ctSize*numNeighbors, &ctStartNewRem);CHKERRQ(ierr); 2249 ierr = MPI_Type_contiguous(ctSize, MPIU_INT, &ctType);CHKERRQ(ierr); 2250 ierr = MPI_Type_commit(&ctType);CHKERRQ(ierr); 2251 ierr = PetscSFBcastBegin(sfProcess, ctType, cr->ctStart, ctStartRem);CHKERRQ(ierr); 2252 ierr = PetscSFBcastEnd(sfProcess, ctType, cr->ctStart, ctStartRem);CHKERRQ(ierr); 2253 ierr = PetscSFBcastBegin(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);CHKERRQ(ierr); 2254 ierr = PetscSFBcastEnd(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);CHKERRQ(ierr); 2255 ierr = MPI_Type_free(&ctType);CHKERRQ(ierr); 2256 ierr = PetscSFDestroy(&sfProcess);CHKERRQ(ierr); 2257 ierr = PetscMalloc1(numNeighbors, &crRem);CHKERRQ(ierr); 2258 for (n = 0; n < numNeighbors; ++n) { 2259 ierr = DMPlexCellRefinerCreate(dm, &crRem[n]);CHKERRQ(ierr); 2260 ierr = DMPlexCellRefinerSetStarts(crRem[n], &ctStartRem[n*ctSize], &ctStartNewRem[n*ctSize]); 2261 ierr = DMPlexCellRefinerSetUp(crRem[n]);CHKERRQ(ierr); 2262 } 2263 ierr = PetscFree2(ctStartRem, ctStartNewRem);CHKERRQ(ierr); 2264 /* Calculate new point SF */ 2265 ierr = PetscMalloc1(numLeavesNew, &localPointsNew);CHKERRQ(ierr); 2266 ierr = PetscMalloc1(numLeavesNew, &remotePointsNew);CHKERRQ(ierr); 2267 ierr = ISGetIndices(processRanks, &neighbors);CHKERRQ(ierr); 2268 for (l = 0, m = 0; l < numLeaves; ++l) { 2269 PetscInt p = localPoints[l]; 2270 PetscInt pRem = remotePoints[l].index; 2271 PetscMPIInt rankRem = remotePoints[l].rank; 2272 DMPolytopeType ct; 2273 DMPolytopeType *rct; 2274 PetscInt *rsize, *rcone, *rornt; 2275 PetscInt neighbor, Nct, n, r; 2276 2277 ierr = PetscFindInt(rankRem, numNeighbors, neighbors, &neighbor);CHKERRQ(ierr); 2278 if (neighbor < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %D", rankRem); 2279 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2280 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2281 for (n = 0; n < Nct; ++n) { 2282 for (r = 0; r < rsize[n]; ++r) { 2283 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 2284 ierr = DMPlexCellRefinerGetNewPoint(crRem[neighbor], ct, rct[n], pRem, r, &pNewRem);CHKERRQ(ierr); 2285 localPointsNew[m] = pNew; 2286 remotePointsNew[m].index = pNewRem; 2287 remotePointsNew[m].rank = rankRem; 2288 ++m; 2289 } 2290 } 2291 } 2292 for (n = 0; n < numNeighbors; ++n) {ierr = DMPlexCellRefinerDestroy(&crRem[n]);CHKERRQ(ierr);} 2293 ierr = PetscFree(crRem);CHKERRQ(ierr); 2294 if (m != numLeavesNew) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of leaf point %D should be %D", m, numLeavesNew); 2295 ierr = ISRestoreIndices(processRanks, &neighbors);CHKERRQ(ierr); 2296 ierr = ISDestroy(&processRanks);CHKERRQ(ierr); 2297 { 2298 PetscSFNode *rp, *rtmp; 2299 PetscInt *lp, *idx, *ltmp, i; 2300 2301 /* SF needs sorted leaves to correct calculate Gather */ 2302 ierr = PetscMalloc1(numLeavesNew, &idx);CHKERRQ(ierr); 2303 ierr = PetscMalloc1(numLeavesNew, &lp);CHKERRQ(ierr); 2304 ierr = PetscMalloc1(numLeavesNew, &rp);CHKERRQ(ierr); 2305 for (i = 0; i < numLeavesNew; ++i) { 2306 if ((localPointsNew[i] < pStartNew) || (localPointsNew[i] >= pEndNew)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local SF point %D (%D) not in [%D, %D)", localPointsNew[i], i, pStartNew, pEndNew); 2307 idx[i] = i; 2308 } 2309 ierr = PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);CHKERRQ(ierr); 2310 for (i = 0; i < numLeavesNew; ++i) { 2311 lp[i] = localPointsNew[idx[i]]; 2312 rp[i] = remotePointsNew[idx[i]]; 2313 } 2314 ltmp = localPointsNew; 2315 localPointsNew = lp; 2316 rtmp = remotePointsNew; 2317 remotePointsNew = rp; 2318 ierr = PetscFree(idx);CHKERRQ(ierr); 2319 ierr = PetscFree(ltmp);CHKERRQ(ierr); 2320 ierr = PetscFree(rtmp);CHKERRQ(ierr); 2321 } 2322 ierr = PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2323 PetscFunctionReturn(0); 2324 } 2325 2326 static PetscErrorCode RefineLabel_Internal(DMPlexCellRefiner cr, DMLabel label, DMLabel labelNew) 2327 { 2328 DM dm = cr->dm; 2329 IS valueIS; 2330 const PetscInt *values; 2331 PetscInt defVal, Nv, val; 2332 PetscErrorCode ierr; 2333 2334 PetscFunctionBegin; 2335 ierr = DMLabelGetDefaultValue(label, &defVal);CHKERRQ(ierr); 2336 ierr = DMLabelSetDefaultValue(labelNew, defVal);CHKERRQ(ierr); 2337 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 2338 ierr = ISGetLocalSize(valueIS, &Nv);CHKERRQ(ierr); 2339 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 2340 for (val = 0; val < Nv; ++val) { 2341 IS pointIS; 2342 const PetscInt *points; 2343 PetscInt numPoints, p; 2344 2345 /* Ensure refined label is created with same number of strata as 2346 * original (even if no entries here). */ 2347 ierr = DMLabelAddStratum(labelNew, values[val]);CHKERRQ(ierr); 2348 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 2349 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 2350 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2351 for (p = 0; p < numPoints; ++p) { 2352 const PetscInt point = points[p]; 2353 DMPolytopeType ct; 2354 DMPolytopeType *rct; 2355 PetscInt *rsize, *rcone, *rornt; 2356 PetscInt Nct, n, r, pNew; 2357 2358 ierr = DMPlexGetCellType(dm, point, &ct);CHKERRQ(ierr); 2359 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2360 for (n = 0; n < Nct; ++n) { 2361 for (r = 0; r < rsize[n]; ++r) { 2362 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], point, r, &pNew);CHKERRQ(ierr); 2363 ierr = DMLabelSetValue(labelNew, pNew, values[val]);CHKERRQ(ierr); 2364 } 2365 } 2366 } 2367 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2368 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2369 } 2370 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 2371 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 2372 PetscFunctionReturn(0); 2373 } 2374 2375 static PetscErrorCode DMPlexCellRefinerCreateLabels(DMPlexCellRefiner cr, DM rdm) 2376 { 2377 DM dm = cr->dm; 2378 PetscInt numLabels, l; 2379 PetscErrorCode ierr; 2380 2381 PetscFunctionBegin; 2382 ierr = DMGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 2383 for (l = 0; l < numLabels; ++l) { 2384 DMLabel label, labelNew; 2385 const char *lname; 2386 PetscBool isDepth, isCellType; 2387 2388 ierr = DMGetLabelName(dm, l, &lname);CHKERRQ(ierr); 2389 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 2390 if (isDepth) continue; 2391 ierr = PetscStrcmp(lname, "celltype", &isCellType);CHKERRQ(ierr); 2392 if (isCellType) continue; 2393 ierr = DMCreateLabel(rdm, lname);CHKERRQ(ierr); 2394 ierr = DMGetLabel(dm, lname, &label);CHKERRQ(ierr); 2395 ierr = DMGetLabel(rdm, lname, &labelNew);CHKERRQ(ierr); 2396 ierr = RefineLabel_Internal(cr, label, labelNew);CHKERRQ(ierr); 2397 } 2398 PetscFunctionReturn(0); 2399 } 2400 2401 /* This will only work for interpolated meshes */ 2402 PetscErrorCode DMPlexRefineUniform(DM dm, DMPlexCellRefiner cr, DM *dmRefined) 2403 { 2404 DM rdm; 2405 PetscInt dim, embedDim, depth; 2406 PetscErrorCode ierr; 2407 2408 PetscFunctionBegin; 2409 PetscValidHeader(cr, 1); 2410 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &rdm);CHKERRQ(ierr); 2411 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2412 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2413 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2414 ierr = DMGetCoordinateDim(dm, &embedDim);CHKERRQ(ierr); 2415 ierr = DMSetCoordinateDim(rdm, embedDim);CHKERRQ(ierr); 2416 /* Calculate number of new points of each depth */ 2417 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2418 if (depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated for regular refinement"); 2419 /* Step 1: Set chart */ 2420 ierr = DMPlexSetChart(rdm, 0, cr->ctStartNew[cr->ctOrder[DM_NUM_POLYTOPES]]);CHKERRQ(ierr); 2421 /* Step 2: Set cone/support sizes (automatically stratifies) */ 2422 ierr = DMPlexCellRefinerSetConeSizes(cr, rdm);CHKERRQ(ierr); 2423 /* Step 3: Setup refined DM */ 2424 ierr = DMSetUp(rdm);CHKERRQ(ierr); 2425 /* Step 4: Set cones and supports (automatically symmetrizes) */ 2426 ierr = DMPlexCellRefinerSetCones(cr, rdm);CHKERRQ(ierr); 2427 /* Step 5: Create pointSF */ 2428 ierr = DMPlexCellRefinerCreateSF(cr, rdm);CHKERRQ(ierr); 2429 /* Step 6: Create labels */ 2430 ierr = DMPlexCellRefinerCreateLabels(cr, rdm);CHKERRQ(ierr); 2431 /* Step 7: Set coordinates */ 2432 ierr = DMPlexCellRefinerSetCoordinates(cr, rdm);CHKERRQ(ierr); 2433 *dmRefined = rdm; 2434 PetscFunctionReturn(0); 2435 } 2436 2437 /*@ 2438 DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data 2439 2440 Input Parameter: 2441 . dm - The coarse DM 2442 2443 Output Parameter: 2444 . fpointIS - The IS of all the fine points which exist in the original coarse mesh 2445 2446 Level: developer 2447 2448 .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS() 2449 @*/ 2450 PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS) 2451 { 2452 DMPlexCellRefiner cr; 2453 PetscInt *fpoints; 2454 PetscInt pStart, pEnd, p, vStart, vEnd, v; 2455 PetscErrorCode ierr; 2456 2457 PetscFunctionBegin; 2458 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2459 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2460 ierr = DMPlexCellRefinerCreate(dm, &cr);CHKERRQ(ierr); 2461 ierr = DMPlexCellRefinerSetUp(cr);CHKERRQ(ierr); 2462 ierr = PetscMalloc1(pEnd-pStart, &fpoints);CHKERRQ(ierr); 2463 for (p = 0; p < pEnd-pStart; ++p) fpoints[p] = -1; 2464 for (v = vStart; v < vEnd; ++v) { 2465 PetscInt vNew = -1; /* silent overzelous may be used uninitialized */ 2466 2467 ierr = DMPlexCellRefinerGetNewPoint(cr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);CHKERRQ(ierr); 2468 fpoints[v-pStart] = vNew; 2469 } 2470 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr); 2471 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);CHKERRQ(ierr); 2472 PetscFunctionReturn(0); 2473 } 2474 2475 /*@ 2476 DMPlexSetRefinementUniform - Set the flag for uniform refinement 2477 2478 Input Parameters: 2479 + dm - The DM 2480 - refinementUniform - The flag for uniform refinement 2481 2482 Level: developer 2483 2484 .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 2485 @*/ 2486 PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform) 2487 { 2488 DM_Plex *mesh = (DM_Plex*) dm->data; 2489 2490 PetscFunctionBegin; 2491 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2492 mesh->refinementUniform = refinementUniform; 2493 PetscFunctionReturn(0); 2494 } 2495 2496 /*@ 2497 DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement 2498 2499 Input Parameter: 2500 . dm - The DM 2501 2502 Output Parameter: 2503 . refinementUniform - The flag for uniform refinement 2504 2505 Level: developer 2506 2507 .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 2508 @*/ 2509 PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform) 2510 { 2511 DM_Plex *mesh = (DM_Plex*) dm->data; 2512 2513 PetscFunctionBegin; 2514 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2515 PetscValidPointer(refinementUniform, 2); 2516 *refinementUniform = mesh->refinementUniform; 2517 PetscFunctionReturn(0); 2518 } 2519 2520 /*@ 2521 DMPlexSetRefinementLimit - Set the maximum cell volume for refinement 2522 2523 Input Parameters: 2524 + dm - The DM 2525 - refinementLimit - The maximum cell volume in the refined mesh 2526 2527 Level: developer 2528 2529 .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform() 2530 @*/ 2531 PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit) 2532 { 2533 DM_Plex *mesh = (DM_Plex*) dm->data; 2534 2535 PetscFunctionBegin; 2536 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2537 mesh->refinementLimit = refinementLimit; 2538 PetscFunctionReturn(0); 2539 } 2540 2541 /*@ 2542 DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement 2543 2544 Input Parameter: 2545 . dm - The DM 2546 2547 Output Parameter: 2548 . refinementLimit - The maximum cell volume in the refined mesh 2549 2550 Level: developer 2551 2552 .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform() 2553 @*/ 2554 PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit) 2555 { 2556 DM_Plex *mesh = (DM_Plex*) dm->data; 2557 2558 PetscFunctionBegin; 2559 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2560 PetscValidPointer(refinementLimit, 2); 2561 /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */ 2562 *refinementLimit = mesh->refinementLimit; 2563 PetscFunctionReturn(0); 2564 } 2565 2566 /*@ 2567 DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement 2568 2569 Input Parameters: 2570 + dm - The DM 2571 - refinementFunc - Function giving the maximum cell volume in the refined mesh 2572 2573 Note: The calling sequence is refinementFunc(coords, limit) 2574 $ coords - Coordinates of the current point, usually a cell centroid 2575 $ limit - The maximum cell volume for a cell containing this point 2576 2577 Level: developer 2578 2579 .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 2580 @*/ 2581 PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *)) 2582 { 2583 DM_Plex *mesh = (DM_Plex*) dm->data; 2584 2585 PetscFunctionBegin; 2586 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2587 mesh->refinementFunc = refinementFunc; 2588 PetscFunctionReturn(0); 2589 } 2590 2591 /*@ 2592 DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement 2593 2594 Input Parameter: 2595 . dm - The DM 2596 2597 Output Parameter: 2598 . refinementFunc - Function giving the maximum cell volume in the refined mesh 2599 2600 Note: The calling sequence is refinementFunc(coords, limit) 2601 $ coords - Coordinates of the current point, usually a cell centroid 2602 $ limit - The maximum cell volume for a cell containing this point 2603 2604 Level: developer 2605 2606 .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 2607 @*/ 2608 PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *)) 2609 { 2610 DM_Plex *mesh = (DM_Plex*) dm->data; 2611 2612 PetscFunctionBegin; 2613 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2614 PetscValidPointer(refinementFunc, 2); 2615 *refinementFunc = mesh->refinementFunc; 2616 PetscFunctionReturn(0); 2617 } 2618 2619 static PetscErrorCode RefineDiscLabels_Internal(DMPlexCellRefiner cr, DM rdm) 2620 { 2621 DM dm = cr->dm; 2622 PetscInt Nf, f, Nds, s; 2623 PetscErrorCode ierr; 2624 2625 PetscFunctionBegin; 2626 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 2627 for (f = 0; f < Nf; ++f) { 2628 DMLabel label, labelNew; 2629 PetscObject obj; 2630 const char *lname; 2631 2632 ierr = DMGetField(rdm, f, &label, &obj);CHKERRQ(ierr); 2633 if (!label) continue; 2634 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr); 2635 ierr = DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);CHKERRQ(ierr); 2636 ierr = RefineLabel_Internal(cr, label, labelNew);CHKERRQ(ierr); 2637 ierr = DMSetField_Internal(rdm, f, labelNew, obj);CHKERRQ(ierr); 2638 ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr); 2639 } 2640 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr); 2641 for (s = 0; s < Nds; ++s) { 2642 DMLabel label, labelNew; 2643 const char *lname; 2644 2645 ierr = DMGetRegionNumDS(rdm, s, &label, NULL, NULL);CHKERRQ(ierr); 2646 if (!label) continue; 2647 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr); 2648 ierr = DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);CHKERRQ(ierr); 2649 ierr = RefineLabel_Internal(cr, label, labelNew);CHKERRQ(ierr); 2650 ierr = DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);CHKERRQ(ierr); 2651 ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr); 2652 } 2653 PetscFunctionReturn(0); 2654 } 2655 2656 PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined) 2657 { 2658 PetscBool isUniform; 2659 DMPlexCellRefiner cr; 2660 PetscErrorCode ierr; 2661 2662 PetscFunctionBegin; 2663 ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr); 2664 ierr = DMViewFromOptions(dm, NULL, "-initref_dm_view");CHKERRQ(ierr); 2665 if (isUniform) { 2666 PetscBool localized; 2667 2668 ierr = DMPlexCellRefinerCreate(dm, &cr);CHKERRQ(ierr); 2669 ierr = DMPlexCellRefinerSetUp(cr);CHKERRQ(ierr); 2670 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 2671 ierr = DMPlexRefineUniform(dm, cr, dmRefined);CHKERRQ(ierr); 2672 ierr = DMPlexSetRegularRefinement(*dmRefined, PETSC_TRUE);CHKERRQ(ierr); 2673 ierr = DMCopyDisc(dm, *dmRefined);CHKERRQ(ierr); 2674 ierr = RefineDiscLabels_Internal(cr, *dmRefined);CHKERRQ(ierr); 2675 ierr = DMCopyBoundary(dm, *dmRefined);CHKERRQ(ierr); 2676 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr); 2677 } else { 2678 ierr = DMPlexRefine_Internal(dm, NULL, dmRefined);CHKERRQ(ierr); 2679 } 2680 PetscFunctionReturn(0); 2681 } 2682 2683 PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]) 2684 { 2685 DM cdm = dm; 2686 PetscInt r; 2687 PetscBool isUniform, localized; 2688 PetscErrorCode ierr; 2689 2690 PetscFunctionBegin; 2691 ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr); 2692 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 2693 if (isUniform) { 2694 for (r = 0; r < nlevels; ++r) { 2695 DMPlexCellRefiner cr; 2696 2697 ierr = DMPlexCellRefinerCreate(cdm, &cr);CHKERRQ(ierr); 2698 ierr = DMPlexCellRefinerSetUp(cr);CHKERRQ(ierr); 2699 ierr = DMPlexRefineUniform(cdm, cr, &dmRefined[r]);CHKERRQ(ierr); 2700 ierr = DMSetCoarsenLevel(dmRefined[r], cdm->leveldown);CHKERRQ(ierr); 2701 ierr = DMSetRefineLevel(dmRefined[r], cdm->levelup+1);CHKERRQ(ierr); 2702 ierr = DMCopyDisc(cdm, dmRefined[r]);CHKERRQ(ierr); 2703 ierr = RefineDiscLabels_Internal(cr, dmRefined[r]);CHKERRQ(ierr); 2704 ierr = DMCopyBoundary(cdm, dmRefined[r]);CHKERRQ(ierr); 2705 ierr = DMSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr); 2706 ierr = DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);CHKERRQ(ierr); 2707 cdm = dmRefined[r]; 2708 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr); 2709 } 2710 } else { 2711 for (r = 0; r < nlevels; ++r) { 2712 ierr = DMRefine(cdm, PetscObjectComm((PetscObject) dm), &dmRefined[r]);CHKERRQ(ierr); 2713 ierr = DMCopyDisc(cdm, dmRefined[r]);CHKERRQ(ierr); 2714 ierr = DMCopyBoundary(cdm, dmRefined[r]);CHKERRQ(ierr); 2715 if (localized) {ierr = DMLocalizeCoordinates(dmRefined[r]);CHKERRQ(ierr);} 2716 ierr = DMSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr); 2717 cdm = dmRefined[r]; 2718 } 2719 } 2720 PetscFunctionReturn(0); 2721 } 2722