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", "Alfeld2D", "Alfeld3D", "PowellSabin", "BoundaryLayer", "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 /* 466 This is the group multiplication table for the dihedral group of the cell. 467 */ 468 static PetscErrorCode ComposeOrientation_Private(PetscInt n, PetscInt o1, PetscInt o2, PetscInt *o) 469 { 470 PetscFunctionBeginHot; 471 if (!n) {*o = 0;} 472 else if (o1 >= 0 && o2 >= 0) {*o = ( o1 + o2) % n;} 473 else if (o1 < 0 && o2 < 0) {*o = (-o1 - o2) % n;} 474 else if (o1 < 0) {*o = -((-(o1+1) + o2) % n + 1);} 475 else if (o2 < 0) {*o = -((-(o2+1) + o1) % n + 1);} 476 PetscFunctionReturn(0); 477 } 478 479 static PetscErrorCode DMPlexCellRefinerMapSubcells_None(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 480 { 481 PetscErrorCode ierr; 482 483 PetscFunctionBeginHot; 484 *rnew = r; 485 ierr = ComposeOrientation_Private(DMPolytopeTypeGetConeSize(ct), po, o, onew);CHKERRQ(ierr); 486 PetscFunctionReturn(0); 487 } 488 489 static PetscErrorCode DMPlexCellRefinerMapSubcells_Regular(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 490 { 491 /* We shift any input orientation in order to make it non-negative 492 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) 493 The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po 494 Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po 495 */ 496 PetscInt tri_seg_o[] = {-2, 0, 497 -2, 0, 498 -2, 0, 499 0, -2, 500 0, -2, 501 0, -2}; 502 PetscInt tri_seg_r[] = {1, 0, 2, 503 0, 2, 1, 504 2, 1, 0, 505 0, 1, 2, 506 1, 2, 0, 507 2, 0, 1}; 508 PetscInt tri_tri_o[] = {0, 1, 2, -3, -2, -1, 509 2, 0, 1, -2, -1, -3, 510 1, 2, 0, -1, -3, -2, 511 -3, -2, -1, 0, 1, 2, 512 -1, -3, -2, 1, 2, 0, 513 -2, -1, -3, 2, 0, 1}; 514 /* orientation if the replica is the center triangle */ 515 PetscInt tri_tri_o_c[] = {2, 0, 1, -2, -1, -3, 516 1, 2, 0, -1, -3, -2, 517 0, 1, 2, -3, -2, -1, 518 -3, -2, -1, 0, 1, 2, 519 -1, -3, -2, 1, 2, 0, 520 -2, -1, -3, 2, 0, 1}; 521 PetscInt tri_tri_r[] = {0, 2, 1, 3, 522 2, 1, 0, 3, 523 1, 0, 2, 3, 524 0, 1, 2, 3, 525 1, 2, 0, 3, 526 2, 0, 1, 3}; 527 PetscInt quad_seg_r[] = {3, 2, 1, 0, 528 2, 1, 0, 3, 529 1, 0, 3, 2, 530 0, 3, 2, 1, 531 0, 1, 2, 3, 532 1, 2, 3, 0, 533 2, 3, 0, 1, 534 3, 0, 1, 2}; 535 PetscInt quad_quad_o[] = { 0, 1, 2, 3, -4, -3, -2, -1, 536 4, 0, 1, 2, -3, -2, -1, -4, 537 3, 4, 0, 1, -2, -1, -4, -3, 538 2, 3, 4, 0, -1, -4, -3, -2, 539 -4, -3, -2, -1, 0, 1, 2, 3, 540 -1, -4, -3, -2, 1, 2, 3, 0, 541 -2, -1, -4, -3, 2, 3, 0, 1, 542 -3, -2, -1, -4, 3, 0, 1, 2}; 543 PetscInt quad_quad_r[] = {0, 3, 2, 1, 544 3, 2, 1, 0, 545 2, 1, 0, 3, 546 1, 0, 3, 2, 547 0, 1, 2, 3, 548 1, 2, 3, 0, 549 2, 3, 0, 1, 550 3, 0, 1, 2}; 551 PetscInt tquad_tquad_o[] = { 0, 1, -2, -1, 552 1, 0, -1, -2, 553 -2, -1, 0, 1, 554 -1, -2, 1, 0}; 555 PetscInt tquad_tquad_r[] = {1, 0, 556 1, 0, 557 0, 1, 558 0, 1}; 559 560 PetscFunctionBeginHot; 561 /* The default is no transformation */ 562 *rnew = r; 563 *onew = o; 564 switch (pct) { 565 case DM_POLYTOPE_SEGMENT: 566 if (ct == DM_POLYTOPE_SEGMENT) { 567 if (po == 0 || po == -1) {*rnew = r; *onew = o;} 568 else if (po == 1 || po == -2) {*rnew = (r+1)%2; *onew = (o == 0 || o == -1) ? -2 : 0;} 569 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for segment", po); 570 } 571 break; 572 case DM_POLYTOPE_TRIANGLE: 573 switch (ct) { 574 case DM_POLYTOPE_SEGMENT: 575 if (o == -1) o = 0; 576 if (o == -2) o = 1; 577 *onew = tri_seg_o[(po+3)*2+o]; 578 *rnew = tri_seg_r[(po+3)*3+r]; 579 break; 580 case DM_POLYTOPE_TRIANGLE: 581 *onew = r == 3 && po < 0 ? tri_tri_o_c[((po+3)%3)*6+o+3] : tri_tri_o[(po+3)*6+o+3]; 582 *rnew = tri_tri_r[(po+3)*4+r]; 583 break; 584 default: break; 585 } 586 break; 587 case DM_POLYTOPE_QUADRILATERAL: 588 switch (ct) { 589 case DM_POLYTOPE_SEGMENT: 590 *onew = o; 591 *rnew = quad_seg_r[(po+4)*4+r]; 592 break; 593 case DM_POLYTOPE_QUADRILATERAL: 594 *onew = quad_quad_o[(po+4)*8+o+4]; 595 *rnew = quad_quad_r[(po+4)*4+r]; 596 break; 597 default: break; 598 } 599 break; 600 case DM_POLYTOPE_SEG_PRISM_TENSOR: 601 switch (ct) { 602 /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */ 603 case DM_POLYTOPE_SEG_PRISM_TENSOR: 604 *onew = tquad_tquad_o[(po+2)*4+o+2]; 605 *rnew = tquad_tquad_r[(po+2)*2+r]; 606 break; 607 default: break; 608 } 609 break; 610 default: break; 611 } 612 PetscFunctionReturn(0); 613 } 614 615 static PetscErrorCode DMPlexCellRefinerMapSubcells_ToBox(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 616 { 617 PetscErrorCode ierr; 618 /* We shift any input orientation in order to make it non-negative 619 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) 620 The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po 621 Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po 622 */ 623 PetscInt tri_seg_o[] = {0, -2, 624 0, -2, 625 0, -2, 626 0, -2, 627 0, -2, 628 0, -2}; 629 PetscInt tri_seg_r[] = {2, 1, 0, 630 1, 0, 2, 631 0, 2, 1, 632 0, 1, 2, 633 1, 2, 0, 634 2, 0, 1}; 635 PetscInt tri_tri_o[] = {0, 1, 2, 3, -4, -3, -2, -1, 636 3, 0, 1, 2, -3, -2, -1, -4, 637 1, 2, 3, 0, -1, -4, -3, -2, 638 -4, -3, -2, -1, 0, 1, 2, 3, 639 -1, -4, -3, -2, 1, 2, 3, 0, 640 -3, -2, -1, -4, 3, 0, 1, 2}; 641 PetscInt tri_tri_r[] = {0, 2, 1, 642 2, 1, 0, 643 1, 0, 2, 644 0, 1, 2, 645 1, 2, 0, 646 2, 0, 1}; 647 PetscInt tquad_tquad_o[] = { 0, 1, -2, -1, 648 1, 0, -1, -2, 649 -2, -1, 0, 1, 650 -1, -2, 1, 0}; 651 PetscInt tquad_tquad_r[] = {1, 0, 652 1, 0, 653 0, 1, 654 0, 1}; 655 PetscInt tquad_quad_o[] = {-2, -1, -4, -3, 2, 3, 0, 1, 656 1, 2, 3, 0, -1, -4, -3, -2, 657 -4, -3, -2, -1, 0, 1, 2, 3, 658 1, 0, 3, 2, -3, -4, -1, -2}; 659 660 PetscFunctionBeginHot; 661 *rnew = r; 662 *onew = o; 663 switch (pct) { 664 case DM_POLYTOPE_TRIANGLE: 665 switch (ct) { 666 case DM_POLYTOPE_SEGMENT: 667 if (o == -1) o = 0; 668 if (o == -2) o = 1; 669 *onew = tri_seg_o[(po+3)*2+o]; 670 *rnew = tri_seg_r[(po+3)*3+r]; 671 break; 672 case DM_POLYTOPE_QUADRILATERAL: 673 *onew = tri_tri_o[(po+3)*8+o+4]; 674 /* TODO See sheet, for fixing po == 1 and 2 */ 675 if (po == 2 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4]; 676 if (po == 2 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)]; 677 if (po == 1 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4]; 678 if (po == 1 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)]; 679 if (po == -1 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4]; 680 if (po == -1 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)]; 681 if (po == -2 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4]; 682 if (po == -2 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)]; 683 *rnew = tri_tri_r[(po+3)*3+r]; 684 break; 685 default: break; 686 } 687 break; 688 case DM_POLYTOPE_SEG_PRISM_TENSOR: 689 switch (ct) { 690 /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */ 691 case DM_POLYTOPE_SEG_PRISM_TENSOR: 692 *onew = tquad_tquad_o[(po+2)*4+o+2]; 693 *rnew = tquad_tquad_r[(po+2)*2+r]; 694 break; 695 case DM_POLYTOPE_QUADRILATERAL: 696 *onew = tquad_quad_o[(po+2)*8+o+4]; 697 *rnew = tquad_tquad_r[(po+2)*2+r]; 698 break; 699 default: break; 700 } 701 break; 702 default: 703 ierr = DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew);CHKERRQ(ierr); 704 } 705 PetscFunctionReturn(0); 706 } 707 708 static PetscErrorCode DMPlexCellRefinerMapSubcells_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 709 { 710 return DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew); 711 } 712 713 static PetscErrorCode DMPlexCellRefinerMapSubcells_BL(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 714 { 715 /* We shift any input orientation in order to make it non-negative 716 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) 717 The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po 718 Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po 719 */ 720 PetscInt tquad_seg_o[] = { 0, 1, -2, -1, 721 0, 1, -2, -1, 722 -2, -1, 0, 1, 723 -2, -1, 0, 1}; 724 PetscInt tquad_tquad_o[] = { 0, 1, -2, -1, 725 1, 0, -1, -2, 726 -2, -1, 0, 1, 727 -1, -2, 1, 0}; 728 PetscInt tquad_tquad_r[] = {1, 0, 729 0, 1, 730 0, 1, 731 1, 0}; 732 PetscErrorCode ierr; 733 734 PetscFunctionBeginHot; 735 *rnew = r; 736 *onew = o; 737 switch (pct) { 738 case DM_POLYTOPE_POINT_PRISM_TENSOR: 739 if (ct == DM_POLYTOPE_POINT_PRISM_TENSOR) { 740 if (po == 0 || po == -1) {*rnew = r; *onew = o;} 741 else if (po == 1 || po == -2) {*rnew = (r+1)%2; *onew = (o == 0 || o == -1) ? -2 : 0;} 742 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for tensor segment", po); 743 } 744 break; 745 case DM_POLYTOPE_SEG_PRISM_TENSOR: 746 switch (ct) { 747 case DM_POLYTOPE_SEGMENT: 748 *onew = tquad_seg_o[(po+2)*4+o+2]; 749 break; 750 case DM_POLYTOPE_SEG_PRISM_TENSOR: 751 *onew = tquad_tquad_o[(po+2)*4+o+2]; 752 *rnew = tquad_tquad_r[(po+2)*2+r]; 753 break; 754 default: break; 755 } 756 break; 757 default: 758 ierr = DMPlexCellRefinerMapSubcells_None(cr, pct, po, ct, r, o, rnew, onew);CHKERRQ(ierr); 759 } 760 PetscFunctionReturn(0); 761 } 762 763 /*@ 764 DMPlexCellRefinerRefine - Return a description of the refinement for a given cell type 765 766 Input Parameter: 767 . source - The cell type for a source point 768 769 Output Parameter: 770 + Nt - The number of cell types generated by refinement 771 . target - The cell types generated 772 . size - The number of subcells of each type, ordered by dimension 773 . cone - A list of the faces for each subcell of the same type as source 774 - ornt - A list of the face orientations for each subcell of the same type as source 775 776 Note: The cone array gives the cone of each subcell listed by the first three outputs. For the each cone point, we 777 need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical 778 division (described in these outputs) of the cell in the original mesh. The point identifier is given by 779 $ the number of cones to be taken, or 0 for the current cell 780 $ the cell cone point number at each level from which it is subdivided 781 $ the replica number r of the subdivision. 782 The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is 783 $ Nt = 2 784 $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT} 785 $ size = {1, 2} 786 $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0} 787 $ ornt = { 0, 0, 0, 0} 788 789 Level: developer 790 791 .seealso: DMPlexCellRefinerCreate(), DMPlexRefineUniform() 792 @*/ 793 PetscErrorCode DMPlexCellRefinerRefine(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 794 { 795 PetscErrorCode ierr; 796 797 PetscFunctionBeginHot; 798 ierr = (*cr->ops->refine)(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 799 PetscFunctionReturn(0); 800 } 801 802 static PetscErrorCode DMPlexCellRefinerRefine_None(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 803 { 804 static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; 805 static PetscInt vertexS[] = {1}; 806 static PetscInt vertexC[] = {0}; 807 static PetscInt vertexO[] = {0}; 808 static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT}; 809 static PetscInt edgeS[] = {1}; 810 static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 811 static PetscInt edgeO[] = {0, 0}; 812 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; 813 static PetscInt tedgeS[] = {1}; 814 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 815 static PetscInt tedgeO[] = {0, 0}; 816 static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE}; 817 static PetscInt triS[] = {1}; 818 static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0}; 819 static PetscInt triO[] = {0, 0, 0}; 820 static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL}; 821 static PetscInt quadS[] = {1}; 822 static PetscInt quadC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0}; 823 static PetscInt quadO[] = {0, 0, 0, 0}; 824 static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR}; 825 static PetscInt tquadS[] = {1}; 826 static PetscInt tquadC[] = {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}; 827 static PetscInt tquadO[] = {0, 0, 0, 0}; 828 static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON}; 829 static PetscInt tetS[] = {1}; 830 static PetscInt tetC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0}; 831 static PetscInt tetO[] = {0, 0, 0, 0}; 832 static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON}; 833 static PetscInt hexS[] = {1}; 834 static PetscInt hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, 835 DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0}; 836 static PetscInt hexO[] = {0, 0, 0, 0, 0, 0}; 837 static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM}; 838 static PetscInt tripS[] = {1}; 839 static PetscInt tripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, 840 DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0}; 841 static PetscInt tripO[] = {0, 0, 0, 0, 0}; 842 static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR}; 843 static PetscInt ttripS[] = {1}; 844 static PetscInt ttripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, 845 DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0}; 846 static PetscInt ttripO[] = {0, 0, 0, 0, 0}; 847 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR}; 848 static PetscInt tquadpS[] = {1}; 849 static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, 850 DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0}; 851 static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0}; 852 853 PetscFunctionBegin; 854 switch (source) { 855 case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break; 856 case DM_POLYTOPE_SEGMENT: *Nt = 1; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break; 857 case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break; 858 case DM_POLYTOPE_TRIANGLE: *Nt = 1; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 859 case DM_POLYTOPE_QUADRILATERAL: *Nt = 1; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break; 860 case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 1; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break; 861 case DM_POLYTOPE_TETRAHEDRON: *Nt = 1; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 862 case DM_POLYTOPE_HEXAHEDRON: *Nt = 1; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break; 863 case DM_POLYTOPE_TRI_PRISM: *Nt = 1; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; 864 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 1; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break; 865 case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 1; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break; 866 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 867 } 868 PetscFunctionReturn(0); 869 } 870 871 static PetscErrorCode DMPlexCellRefinerRefine_Regular(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 872 { 873 /* All vertices remain in the refined mesh */ 874 static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; 875 static PetscInt vertexS[] = {1}; 876 static PetscInt vertexC[] = {0}; 877 static PetscInt vertexO[] = {0}; 878 /* Split all edges with a new vertex, making two new 2 edges 879 0--0--0--1--1 880 */ 881 static DMPolytopeType edgeT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}; 882 static PetscInt edgeS[] = {1, 2}; 883 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}; 884 static PetscInt edgeO[] = { 0, 0, 0, 0}; 885 /* Do not split tensor edges */ 886 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; 887 static PetscInt tedgeS[] = {1}; 888 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 889 static PetscInt tedgeO[] = { 0, 0}; 890 /* Add 3 edges inside every triangle, making 4 new triangles. 891 2 892 |\ 893 | \ 894 | \ 895 0 1 896 | C \ 897 | \ 898 | \ 899 2---1---1 900 |\ D / \ 901 1 2 0 0 902 |A \ / B \ 903 0-0-0---1---1 904 */ 905 static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE}; 906 static PetscInt triS[] = {3, 4}; 907 static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 908 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0, 909 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0, 910 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, 911 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, 912 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, 913 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2}; 914 static PetscInt triO[] = {0, 0, 915 0, 0, 916 0, 0, 917 0, -2, 0, 918 0, 0, -2, 919 -2, 0, 0, 920 0, 0, 0}; 921 /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads. 922 3----1----2----0----2 923 | | | 924 0 D 2 C 1 925 | | | 926 3----3----0----1----1 927 | | | 928 1 A 0 B 0 929 | | | 930 0----0----0----1----1 931 */ 932 static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL}; 933 static PetscInt quadS[] = {1, 4, 4}; 934 static PetscInt quadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 935 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 936 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 937 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 938 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, 939 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, 940 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, 941 DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0}; 942 static PetscInt quadO[] = {0, 0, 943 0, 0, 944 0, 0, 945 0, 0, 946 0, 0, -2, 0, 947 0, 0, 0, -2, 948 -2, 0, 0, 0, 949 0, -2, 0, 0}; 950 /* Add 1 edge inside every tensor quad, making 2 new tensor quads 951 2----2----1----3----3 952 | | | 953 | | | 954 | | | 955 4 A 6 B 5 956 | | | 957 | | | 958 | | | 959 0----0----0----1----1 960 */ 961 static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR}; 962 static PetscInt tquadS[] = {1, 2}; 963 static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 964 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, 965 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}; 966 static PetscInt tquadO[] = {0, 0, 967 0, 0, 0, 0, 968 0, 0, 0, 0}; 969 /* Add 1 edge and 8 triangles inside every cell, making 8 new tets 970 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 971 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] 972 called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are 973 [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3] 974 The first four tets just cut off the corners, using the replica notation for new vertices, 975 [v0, (e0, 0), (e2, 0), (e3, 0)] 976 [(e0, 0), v1, (e1, 0), (e4, 0)] 977 [(e2, 0), (e1, 0), v2, (e5, 0)] 978 [(e3, 0), (e4, 0), (e5, 0), v3 ] 979 The next four tets match a vertex to the newly created faces from cutting off those first tets. 980 [(e2, 0), (e3, 0), (e0, 0), (e5, 0)] 981 [(e4, 0), (e1, 0), (e0, 0), (e5, 0)] 982 [(e5, 0), (e0, 0), (e2, 0), (e1, 0)] 983 [(e5, 0), (e0, 0), (e4, 0), (e3, 0)] 984 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 985 [(e2, 0), (e0, 0), (e3, 0)] 986 [(e0, 0), (e1, 0), (e4, 0)] 987 [(e2, 0), (e5, 0), (e1, 0)] 988 [(e3, 0), (e4, 0), (e5, 0)] 989 The next four, from the second group of tets, are 990 [(e2, 0), (e0, 0), (e5, 0)] 991 [(e4, 0), (e0, 0), (e5, 0)] 992 [(e0, 0), (e1, 0), (e5, 0)] 993 [(e5, 0), (e3, 0), (e0, 0)] 994 I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table. 995 */ 996 static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON}; 997 static PetscInt tetS[] = {1, 8, 8}; 998 static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0, 999 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, 1000 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, 1001 DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1, 1002 DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, 1003 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1004 DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, 1005 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0, 1006 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0, 1007 DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0, 1008 DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1, 1009 DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0, 1010 DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2, 1011 DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7, 1012 DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6, 1013 DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3, 1014 DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3}; 1015 static PetscInt tetO[] = {0, 0, 1016 0, 0, 0, 1017 0, 0, 0, 1018 0, 0, 0, 1019 0, 0, 0, 1020 0, 0, -2, 1021 0, 0, -2, 1022 0, -2, -2, 1023 0, -2, 0, 1024 0, 0, 0, 0, 1025 0, 0, 0, 0, 1026 0, 0, 0, 0, 1027 0, 0, 0, 0, 1028 -3, 0, 0, -2, 1029 -2, 1, 0, 0, 1030 -2, -2, -1, 2, 1031 -2, 0, -2, 1}; 1032 /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes 1033 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 1034 [v0, v1, v2, v3] f0 bottom 1035 [v4, v5, v6, v7] f1 top 1036 [v0, v3, v5, v4] f2 front 1037 [v2, v1, v7, v6] f3 back 1038 [v3, v2, v6, v5] f4 right 1039 [v0, v4, v7, v1] f5 left 1040 The eight hexes are divided into four on the bottom, and four on the top, 1041 [v0, (e0, 0), (f0, 0), (e3, 0), (e9, 0), (f2, 0), (c0, 0), (f5, 0)] 1042 [(e0, 0), v1, (e1, 0), (f0, 0), (f5, 0), (c0, 0), (f3, 0), (e10, 0)] 1043 [(f0, 0), (e1, 0), v2, (e2, 0), (c0, 0), (f4, 0), (e11, 0), (f3, 0)] 1044 [(e3, 0), (f0, 0), (e2, 0), v3, (f2, 0), (e8, 0), (f4, 0), (c0, 0)] 1045 [(e9, 0), (f5, 0), (c0, 0), (f2, 0), v4, (e4, 0), (f1, 0), (e7, 0)] 1046 [(f2, 0), (c0, 0), (f4, 0), (e8, 0), (e4, 0), v5, (e5, 0), (f1, 0)] 1047 [(c0, 0), (f3, 0), (e11, 0), (f4, 0), (f1, 0), (e5, 0), v6, (e6, 0)] 1048 [(f5, 0), (e10, 0), (f3, 0), (c0, 0), (e7, 0), (f1, 0), (e6, 0), v7] 1049 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, 1050 [(e9, 0), (f2, 0), (c0, 0), (f5, 0)] 1051 [(f5, 0), (c0, 0), (f3, 0), (e10, 0)] 1052 [(c0, 0), (f4, 0), (e11, 0), (f3, 0)] 1053 [(f2, 0), (e8, 0), (f4, 0), (c0, 0)] 1054 and on the x-z plane, 1055 [(f0, 0), (e0, 0), (f5, 0), (c0, 0)] 1056 [(c0, 0), (f5, 0), (e7, 0), (f1, 0)] 1057 [(f4, 0), (c0, 0), (f1, 0), (e5, 0)] 1058 [(e2, 0), (f0, 0), (c0, 0), (f4, 0)] 1059 and on the y-z plane, 1060 [(e3, 0), (f2, 0), (c0, 0), (f0, 0)] 1061 [(f2, 0), (e4, 0), (f1, 0), (c0, 0)] 1062 [(c0, 0), (f1, 0), (e6, 0), (f3, 0)] 1063 [(f0, 0), (c0, 0), (f3, 0), (e1, 0)] 1064 */ 1065 static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 1066 static PetscInt hexS[] = {1, 6, 12, 8}; 1067 static PetscInt hexC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1068 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 1069 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 1070 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 1071 DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0, 1072 DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0, 1073 DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0, 1074 DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2, 1075 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, 1076 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2, 1077 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0, 1078 DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1, 1079 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, 1080 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0, 1081 DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3, 1082 DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, 1083 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, 1084 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, 1085 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, 1086 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, 1087 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, 1088 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, 1089 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, 1090 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, 1091 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, 1092 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}; 1093 static PetscInt hexO[] = {0, 0, 1094 0, 0, 1095 0, 0, 1096 0, 0, 1097 0, 0, 1098 0, 0, 1099 0, 0, -2, -2, 1100 0, -2, -2, 0, 1101 -2, -2, 0, 0, 1102 -2, 0, 0, -2, 1103 -2, 0, 0, -2, 1104 -2, -2, 0, 0, 1105 0, -2, -2, 0, 1106 0, 0, -2, -2, 1107 0, 0, -2, -2, 1108 -2, 0, 0, -2, 1109 -2, -2, 0, 0, 1110 0, -2, -2, 0, 1111 0, 0, 0, 0, -4, 0, 1112 0, 0, -1, 0, -4, 0, 1113 0, 0, -1, 0, 0, 0, 1114 0, 0, 0, 0, 0, 0, 1115 -4, 0, 0, 0, -4, 0, 1116 -4, 0, 0, 0, 0, 0, 1117 -4, 0, -1, 0, 0, 0, 1118 -4, 0, -1, 0, -4, 0}; 1119 /* Add 3 quads inside every triangular prism, making 4 new prisms. */ 1120 static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM}; 1121 static PetscInt tripS[] = {3, 4, 6, 8}; 1122 static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0, 1123 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0, 1124 DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0, 1125 DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1, 1126 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0, 1127 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, 1128 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, 1129 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1130 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, 1131 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0, 1132 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, 1133 DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, 1134 DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, 1135 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, 1136 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, 1137 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, 1138 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, 1139 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, 1140 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, 1141 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, 1142 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}; 1143 static PetscInt tripO[] = {0, 0, 1144 0, 0, 1145 0, 0, 1146 0, -2, -2, 1147 -2, 0, -2, 1148 -2, -2, 0, 1149 0, 0, 0, 1150 -2, 0, -2, -2, 1151 -2, 0, -2, -2, 1152 -2, 0, -2, -2, 1153 0, -2, -2, 0, 1154 0, -2, -2, 0, 1155 0, -2, -2, 0, 1156 0, 0, 0, -1, 0, 1157 0, 0, 0, 0, -1, 1158 0, 0, -1, 0, 0, 1159 2, 0, 0, 0, 0, 1160 -3, 0, 0, -1, 0, 1161 -3, 0, 0, 0, -1, 1162 -3, 0, -1, 0, 0, 1163 -3, 0, 0, 0, 0}; 1164 /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms. 1165 2 1166 |\ 1167 | \ 1168 | \ 1169 0---1 1170 1171 2 1172 1173 0 1 1174 1175 2 1176 |\ 1177 | \ 1178 | \ 1179 0---1 1180 */ 1181 static DMPolytopeType ttripT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR}; 1182 static PetscInt ttripS[] = {3, 4}; 1183 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, 1184 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, 1185 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, 1186 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, 1187 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, 1188 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, 1189 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}; 1190 static PetscInt ttripO[] = {0, 0, 0, 0, 1191 0, 0, 0, 0, 1192 0, 0, 0, 0, 1193 0, 0, 0, -1, 0, 1194 0, 0, 0, 0, -1, 1195 0, 0, -1, 0, 0, 1196 0, 0, 0, 0, 0}; 1197 /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */ 1198 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1199 static PetscInt tquadpS[] = {1, 4, 4}; 1200 static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1201 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, 1202 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, 1203 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, 1204 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, 1205 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, 1206 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, 1207 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, 1208 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}; 1209 static PetscInt tquadpO[] = {0, 0, 1210 0, 0, 0, 0, 1211 0, 0, 0, 0, 1212 0, 0, 0, 0, 1213 0, 0, 0, 0, 1214 0, 0, 0, 0, -1, 0, 1215 0, 0, 0, 0, 0, -1, 1216 0, 0, -1, 0, 0, 0, 1217 0, 0, 0, -1, 0, 0}; 1218 1219 PetscFunctionBegin; 1220 switch (source) { 1221 case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break; 1222 case DM_POLYTOPE_SEGMENT: *Nt = 2; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break; 1223 case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break; 1224 case DM_POLYTOPE_TRIANGLE: *Nt = 2; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 1225 case DM_POLYTOPE_QUADRILATERAL: *Nt = 3; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break; 1226 case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break; 1227 case DM_POLYTOPE_TETRAHEDRON: *Nt = 3; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 1228 case DM_POLYTOPE_HEXAHEDRON: *Nt = 4; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break; 1229 case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; 1230 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 2; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break; 1231 case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break; 1232 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1233 } 1234 PetscFunctionReturn(0); 1235 } 1236 1237 static PetscErrorCode DMPlexCellRefinerRefine_ToBox(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1238 { 1239 PetscErrorCode ierr; 1240 /* Change tensor edges to segments */ 1241 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_SEGMENT}; 1242 static PetscInt tedgeS[] = {1}; 1243 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 1244 static PetscInt tedgeO[] = { 0, 0}; 1245 /* Add 1 vertex, 3 edges inside every triangle, making 3 new quadrilaterals. 1246 2 1247 |\ 1248 | \ 1249 | \ 1250 | \ 1251 0 1 1252 | \ 1253 | \ 1254 2 1 1255 |\ / \ 1256 | 2 1 \ 1257 | \ / \ 1258 1 | 0 1259 | 0 \ 1260 | | \ 1261 | | \ 1262 0-0-0-----1-----1 1263 */ 1264 static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL}; 1265 static PetscInt triS[] = {1, 3, 3}; 1266 static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1267 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 1268 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 1269 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, 1270 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, 1271 DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0}; 1272 static PetscInt triO[] = {0, 0, 1273 0, 0, 1274 0, 0, 1275 0, 0, -2, 0, 1276 0, 0, 0, -2, 1277 0, -2, 0, 0}; 1278 /* Add 1 edge inside every tensor quad, making 2 new quadrilaterals 1279 2----2----1----3----3 1280 | | | 1281 | | | 1282 | | | 1283 4 A 6 B 5 1284 | | | 1285 | | | 1286 | | | 1287 0----0----0----1----1 1288 */ 1289 static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL}; 1290 static PetscInt tquadS[] = {1, 2}; 1291 static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1292 /* TODO Fix these */ 1293 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1294 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0}; 1295 static PetscInt tquadO[] = {0, 0, 1296 0, 0, -2, -2, 1297 0, 0, -2, -2}; 1298 /* Add 6 triangles inside every cell, making 4 new hexs 1299 TODO: Need different SubcellMap(). Need to make a struct with the function pointers in it 1300 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 1301 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] 1302 called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are 1303 [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3] 1304 We make a new hex in each corner 1305 [v0, (e0, 0), (f0, 0), (e2, 0), (e3, 0), (f2, 0), (c0, 0), (f1, 0)] 1306 [v1, (e4, 0), (f3, 0), (e1, 0), (e0, 0), (f0, 0), (c0, 0), (f1, 0)] 1307 [v2, (e1, 0), (f3, 0), (e5, 0), (e2, 0), (f2, 0), (c0, 0), (f0, 0)] 1308 [v3, (e4, 0), (f1, 0), (e3, 0), (e5, 0), (f2, 0), (c0, 0), (f3, 0)] 1309 We create a new face for each edge 1310 [(e3, 0), (f2, 0), (c0, 0), (f1, 0)] 1311 [(f0, 0), (e0, 0), (f1, 0), (c0, 0)] 1312 [(e2, 0), (f0, 0), (c0, 0), (f2, 0)] 1313 [(f3, 0), (e4, 0), (f1, 0), (c0, 0)] 1314 [(e1, 0), (f3, 0), (c0, 0), (f0, 0)] 1315 [(e5, 0), (f3, 0), (c0, 0), (f2, 0)] 1316 I could write a program to generate these from the first hex by acting with the symmetry group to take one subcell into another. 1317 */ 1318 static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 1319 static PetscInt tetS[] = {1, 4, 6, 4}; 1320 static PetscInt tetC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1321 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 1322 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 1323 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 1324 DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, 1325 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, 1326 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1327 DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 3, 1328 DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, 1329 DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, 1330 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, 1331 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, 1332 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, 1333 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}; 1334 static PetscInt tetO[] = {0, 0, 1335 0, 0, 1336 0, 0, 1337 0, 0, 1338 0, 0, -2, -2, 1339 -2, 0, 0, -2, 1340 0, 0, -2, -2, 1341 -2, 0, 0, -2, 1342 0, 0, -2, -2, 1343 0, 0, -2, -2, 1344 0, 0, 0, 0, 0, 0, 1345 1, -1, 1, 0, 0, 3, 1346 0, -4, 1, -1, 0, 3, 1347 1, -4, 3, -2, -4, 3}; 1348 /* Add 3 quads inside every triangular prism, making 4 new prisms. */ 1349 static DMPolytopeType tripT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 1350 static PetscInt tripS[] = {1, 5, 9, 6}; 1351 static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1352 DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, 1353 DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, 1354 DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, 1355 DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0, 1356 DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, 1357 DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, 1358 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, 1359 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1360 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, 1361 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0, 1362 DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, 1363 DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, 1364 DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, 1365 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, 1366 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, 1367 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, 1368 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, 1369 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, 1370 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}; 1371 static PetscInt tripO[] = {0, 0, 1372 0, 0, 1373 0, 0, 1374 0, 0, 1375 0, 0, 1376 0, 0, -2, -2, 1377 -2, 0, 0, -2, 1378 0, -2, -2, 0, 1379 0, 0, -2, -2, 1380 0, 0, -2, -2, 1381 0, 0, -2, -2, 1382 0, -2, -2, 0, 1383 0, -2, -2, 0, 1384 0, -2, -2, 0, 1385 0, 0, 0, -1, 0, 1, 1386 0, 0, 0, 0, 0, -4, 1387 0, 0, 0, 0, -1, 1, 1388 -4, 0, 0, -1, 0, 1, 1389 -4, 0, 0, 0, 0, -4, 1390 -4, 0, 0, 0, -1, 1}; 1391 /* Add 3 tensor quads inside every tensor triangular prism, making 4 new tensor triangular prisms. 1392 2 1393 |\ 1394 | \ 1395 | \ 1396 0---1 1397 1398 2 1399 1400 0 1 1401 1402 2 1403 |\ 1404 | \ 1405 | \ 1406 0---1 1407 */ 1408 static DMPolytopeType ttripT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1409 static PetscInt ttripS[] = {1, 3, 3}; 1410 static PetscInt ttripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1411 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, 1412 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, 1413 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, 1414 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, 1415 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, 1416 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}; 1417 static PetscInt ttripO[] = {0, 0, 1418 0, 0, 0, 0, 1419 0, 0, 0, 0, 1420 0, 0, 0, 0, 1421 0, 0, 0, 0, -1, 0, 1422 0, 0, 0, 0, 0, -1, 1423 0, 0, 0, -1, 0, 0}; 1424 /* TODO Add 3 quads inside every tensor triangular prism, making 4 new triangular prisms. 1425 2 1426 |\ 1427 | \ 1428 | \ 1429 0---1 1430 1431 2 1432 1433 0 1 1434 1435 2 1436 |\ 1437 | \ 1438 | \ 1439 0---1 1440 */ 1441 static DMPolytopeType ctripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON}; 1442 static PetscInt ctripS[] = {1, 3, 3}; 1443 static PetscInt ctripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1444 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, 1445 DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, 1446 DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0, 1447 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, 1448 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, 1449 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}; 1450 static PetscInt ctripO[] = {0, 0, 1451 0, 0, -2, -2, 1452 0, 0, -2, -2, 1453 0, 0, -2, -2, 1454 -4, 0, 0, -1, 0, 1, 1455 -4, 0, 0, 0, 0, -4, 1456 -4, 0, 0, 0, -1, 1}; 1457 /* Add 1 edge and 4 quads inside every tensor quad prism, making 4 new hexahedra. */ 1458 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1459 static PetscInt tquadpS[] = {1, 4, 4}; 1460 static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, 1461 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, 1462 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, 1463 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, 1464 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, 1465 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, 1466 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, 1467 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, 1468 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}; 1469 static PetscInt tquadpO[] = {0, 0, 1470 0, 0, 0, 0, 1471 0, 0, 0, 0, 1472 0, 0, 0, 0, 1473 0, 0, 0, 0, 1474 0, 0, 0, 0, -1, 0, 1475 0, 0, 0, 0, 0, -1, 1476 0, 0, -1, 0, 0, 0, 1477 0, 0, 0, -1, 0, 0}; 1478 PetscBool convertTensor = PETSC_TRUE; 1479 1480 PetscFunctionBeginHot; 1481 if (convertTensor) { 1482 switch (source) { 1483 case DM_POLYTOPE_POINT: 1484 case DM_POLYTOPE_SEGMENT: 1485 case DM_POLYTOPE_QUADRILATERAL: 1486 case DM_POLYTOPE_HEXAHEDRON: 1487 ierr = DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 1488 break; 1489 case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break; 1490 case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break; 1491 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ctripT; *size = ctripS; *cone = ctripC; *ornt = ctripO; break; 1492 case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break; 1493 case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 1494 case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 1495 case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; 1496 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1497 } 1498 } else { 1499 switch (source) { 1500 case DM_POLYTOPE_POINT: 1501 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1502 case DM_POLYTOPE_SEGMENT: 1503 case DM_POLYTOPE_QUADRILATERAL: 1504 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1505 case DM_POLYTOPE_HEXAHEDRON: 1506 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1507 ierr = DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 1508 break; 1509 case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 1510 case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 1511 case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; 1512 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break; 1513 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1514 } 1515 } 1516 PetscFunctionReturn(0); 1517 } 1518 1519 static PetscErrorCode DMPlexCellRefinerRefine_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1520 { 1521 PetscErrorCode ierr; 1522 1523 PetscFunctionBeginHot; 1524 switch (source) { 1525 case DM_POLYTOPE_POINT: 1526 case DM_POLYTOPE_SEGMENT: 1527 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1528 case DM_POLYTOPE_TRIANGLE: 1529 case DM_POLYTOPE_TETRAHEDRON: 1530 case DM_POLYTOPE_TRI_PRISM: 1531 case DM_POLYTOPE_TRI_PRISM_TENSOR: 1532 case DM_POLYTOPE_QUADRILATERAL: 1533 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1534 case DM_POLYTOPE_HEXAHEDRON: 1535 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1536 ierr = DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 1537 break; 1538 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1539 } 1540 PetscFunctionReturn(0); 1541 } 1542 1543 static PetscErrorCode DMPlexCellRefinerRefine_Alfeld2D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1544 { 1545 PetscErrorCode ierr; 1546 /* Add 1 vertex, 3 edges inside every triangle, making 3 new triangles. 1547 2 1548 |\ 1549 |\\ 1550 | |\ 1551 | \ \ 1552 | | \ 1553 | \ \ 1554 | | \ 1555 2 \ \ 1556 | | 1 1557 | 2 \ 1558 | | \ 1559 | /\ \ 1560 | 0 1 | 1561 | / \ | 1562 |/ \| 1563 0---0----1 1564 */ 1565 static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE}; 1566 static PetscInt triS[] = {1, 3, 3}; 1567 static PetscInt triC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1568 DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1569 DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1570 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, 1571 DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, 1572 DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2}; 1573 static PetscInt triO[] = {0, 0, 1574 0, 0, 1575 0, 0, 1576 0, 0, -2, 1577 0, 0, -2, 1578 0, 0, -2}; 1579 1580 PetscFunctionBeginHot; 1581 switch (source) { 1582 case DM_POLYTOPE_POINT: 1583 case DM_POLYTOPE_SEGMENT: 1584 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1585 case DM_POLYTOPE_QUADRILATERAL: 1586 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1587 case DM_POLYTOPE_TETRAHEDRON: 1588 case DM_POLYTOPE_HEXAHEDRON: 1589 case DM_POLYTOPE_TRI_PRISM: 1590 case DM_POLYTOPE_TRI_PRISM_TENSOR: 1591 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1592 ierr = DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 1593 break; 1594 case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 1595 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1596 } 1597 PetscFunctionReturn(0); 1598 } 1599 1600 static PetscErrorCode DMPlexCellRefinerRefine_Alfeld3D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1601 { 1602 PetscErrorCode ierr; 1603 /* Add 6 triangles inside every cell, making 4 new tets 1604 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 1605 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] 1606 called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are 1607 [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3] 1608 We make a new tet on each face 1609 [v0, v1, v2, (c0, 0)] 1610 [v0, v3, v1, (c0, 0)] 1611 [v0, v2, v3, (c0, 0)] 1612 [v2, v1, v3, (c0, 0)] 1613 We create a new face for each edge 1614 [v0, (c0, 0), v1 ] 1615 [v0, v2, (c0, 0)] 1616 [v2, v1, (c0, 0)] 1617 [v0, (c0, 0), v3 ] 1618 [v1, v3, (c0, 0)] 1619 [v3, v2, (c0, 0)] 1620 */ 1621 static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON}; 1622 static PetscInt tetS[] = {1, 4, 6, 4}; 1623 static PetscInt tetC[] = {DM_POLYTOPE_POINT, 3, 0, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1624 DM_POLYTOPE_POINT, 3, 0, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1625 DM_POLYTOPE_POINT, 3, 0, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0, 1626 DM_POLYTOPE_POINT, 3, 1, 0, 1, 0, DM_POLYTOPE_POINT, 0, 0, 1627 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 2, 0, 0, 0, 1628 DM_POLYTOPE_SEGMENT, 2, 0, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, 1629 DM_POLYTOPE_SEGMENT, 2, 0, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, 1630 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 2, 1, 0, 0, 1631 DM_POLYTOPE_SEGMENT, 2, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1, 1632 DM_POLYTOPE_SEGMENT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3, 1633 DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2, 1634 DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 4, 1635 DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5, 1636 DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 4}; 1637 static PetscInt tetO[] = {0, 0, 1638 0, 0, 1639 0, 0, 1640 0, 0, 1641 0, -2, -2, 1642 -2, 0, -2, 1643 -2, 0, -2, 1644 0, -2, -2, 1645 -2, 0, -2, 1646 -2, 0, -2, 1647 0, 0, 0, 0, 1648 0, 0, -3, 0, 1649 0, -3, -3, 0, 1650 0, -3, -1, -1}; 1651 1652 PetscFunctionBeginHot; 1653 switch (source) { 1654 case DM_POLYTOPE_POINT: 1655 case DM_POLYTOPE_SEGMENT: 1656 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1657 case DM_POLYTOPE_TRIANGLE: 1658 case DM_POLYTOPE_QUADRILATERAL: 1659 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1660 case DM_POLYTOPE_HEXAHEDRON: 1661 case DM_POLYTOPE_TRI_PRISM: 1662 case DM_POLYTOPE_TRI_PRISM_TENSOR: 1663 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1664 ierr = DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 1665 break; 1666 case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 1667 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1668 } 1669 PetscFunctionReturn(0); 1670 } 1671 1672 static PetscErrorCode DMPlexCellRefinerRefine_BL(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1673 { 1674 PetscErrorCode ierr; 1675 /* Split all tensor edges with a new vertex, making two new 2 tensor edges 1676 0--0--0--1--1 1677 */ 1678 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_POINT_PRISM_TENSOR}; 1679 static PetscInt tedgeS[] = {1, 2}; 1680 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 1681 static PetscInt tedgeO[] = { 0, 0, 0, 0}; 1682 /* Add 1 segment inside every tensor quad, making 2 new tensor quads 1683 2---------1---------3 1684 | | 1685 | B | 1686 | | 1687 2---------2---------3 1688 | | 1689 | A | 1690 | | 1691 0---------0---------1 1692 */ 1693 static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_SEG_PRISM_TENSOR}; 1694 static PetscInt tquadS[] = {1, 2}; 1695 static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0, 1696 DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, 1697 DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 1}; 1698 static PetscInt tquadO[] = {0, 0, 1699 0, 0, 0, 0, 1700 0, 0, 0, 0}; 1701 /* Add 1 triangle and make 2 new tensor triangular prisms */ 1702 static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TRI_PRISM_TENSOR}; 1703 static PetscInt ttripS[] = {1, 2}; 1704 static PetscInt ttripC[] = {DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, 1705 DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, 1706 DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1}; 1707 static PetscInt ttripO[] = {0, 0, 0, 1708 0, 0, 0, 0, 0, 1709 0, 0, 0, 0, 0, 1710 0, 0, 0, 0, 0, 1711 0, 0, 0, 0, 0}; 1712 /* Add 1 quad and make 2 new tensor quad prisms */ 1713 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1714 static PetscInt tquadpS[] = {1, 2}; 1715 static PetscInt tquadpC[] = {DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 1, 5, 0, 1716 DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0, 1717 DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 1}; 1718 static PetscInt tquadpO[] = {0, 0, 0, 0, 1719 0, 0, 0, 0, 0, 0, 1720 0, 0, 0, 0, 0, 0, 1721 0, 0, 0, 0, 0, 0, 1722 0, 0, 0, 0, 0, 0}; 1723 1724 PetscFunctionBeginHot; 1725 switch (source) { 1726 case DM_POLYTOPE_POINT: 1727 case DM_POLYTOPE_SEGMENT: 1728 case DM_POLYTOPE_TRIANGLE: 1729 case DM_POLYTOPE_QUADRILATERAL: 1730 case DM_POLYTOPE_TETRAHEDRON: 1731 case DM_POLYTOPE_HEXAHEDRON: 1732 case DM_POLYTOPE_TRI_PRISM: 1733 ierr = DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr); 1734 break; 1735 case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 2; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break; 1736 case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break; 1737 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 2; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break; 1738 case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 2; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break; 1739 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1740 } 1741 PetscFunctionReturn(0); 1742 } 1743 1744 static PetscErrorCode CellRefinerCreateOffset_Internal(DMPlexCellRefiner cr, PetscInt ctOrder[], PetscInt ctStart[], PetscInt **offset) 1745 { 1746 PetscInt c, cN, *off; 1747 PetscErrorCode ierr; 1748 1749 PetscFunctionBegin; 1750 ierr = PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);CHKERRQ(ierr); 1751 for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) { 1752 const DMPolytopeType ct = (DMPolytopeType) c; 1753 for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) { 1754 const DMPolytopeType ctNew = (DMPolytopeType) cN; 1755 DMPolytopeType *rct; 1756 PetscInt *rsize, *cone, *ornt; 1757 PetscInt Nct, n, i; 1758 1759 if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[ct*DM_NUM_POLYTOPES+ctNew] = -1; break;} 1760 off[ct*DM_NUM_POLYTOPES+ctNew] = 0; 1761 for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) { 1762 const DMPolytopeType ict = (DMPolytopeType) ctOrder[i]; 1763 const DMPolytopeType ictn = (DMPolytopeType) ctOrder[i+1]; 1764 1765 ierr = DMPlexCellRefinerRefine(cr, ict, &Nct, &rct, &rsize, &cone, &ornt);CHKERRQ(ierr); 1766 if (ict == ct) { 1767 for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break; 1768 if (n == Nct) off[ct*DM_NUM_POLYTOPES+ctNew] = -1; 1769 break; 1770 } 1771 for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct*DM_NUM_POLYTOPES+ctNew] += (ctStart[ictn]-ctStart[ict]) * rsize[n]; 1772 } 1773 } 1774 } 1775 *offset = off; 1776 PetscFunctionReturn(0); 1777 } 1778 1779 static PetscErrorCode DMPlexCellRefinerSetStarts(DMPlexCellRefiner cr, const PetscInt ctStart[], const PetscInt ctStartNew[]) 1780 { 1781 const PetscInt ctSize = DM_NUM_POLYTOPES+1; 1782 PetscErrorCode ierr; 1783 1784 PetscFunctionBegin; 1785 if (cr->setupcalled) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_ARG_WRONGSTATE, "Must call this function before DMPlexCellRefinerSetUp()"); 1786 ierr = PetscCalloc2(ctSize, &cr->ctStart, ctSize, &cr->ctStartNew);CHKERRQ(ierr); 1787 ierr = PetscArraycpy(cr->ctStart, ctStart, ctSize);CHKERRQ(ierr); 1788 ierr = PetscArraycpy(cr->ctStartNew, ctStartNew, ctSize);CHKERRQ(ierr); 1789 PetscFunctionReturn(0); 1790 } 1791 1792 /* Construct cell type order since we must loop over cell types in depth order */ 1793 PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell, PetscInt *ctOrder[], PetscInt *ctOrderInv[]) 1794 { 1795 PetscInt *ctO, *ctOInv; 1796 PetscInt c, d, off = 0; 1797 PetscErrorCode ierr; 1798 1799 PetscFunctionBegin; 1800 ierr = PetscCalloc2(DM_NUM_POLYTOPES+1, &ctO, DM_NUM_POLYTOPES+1, &ctOInv);CHKERRQ(ierr); 1801 for (d = 3; d >= DMPolytopeTypeGetDim(ctCell); --d) { 1802 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1803 if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue; 1804 ctO[off++] = c; 1805 } 1806 } 1807 if (DMPolytopeTypeGetDim(ctCell) != 0) { 1808 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1809 if (DMPolytopeTypeGetDim((DMPolytopeType) c) != 0) continue; 1810 ctO[off++] = c; 1811 } 1812 } 1813 for (d = DMPolytopeTypeGetDim(ctCell)-1; d > 0; --d) { 1814 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1815 if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue; 1816 ctO[off++] = c; 1817 } 1818 } 1819 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1820 if (DMPolytopeTypeGetDim((DMPolytopeType) c) >= 0) continue; 1821 ctO[off++] = c; 1822 } 1823 if (off != DM_NUM_POLYTOPES+1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %D for cell type order", off); 1824 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 1825 ctOInv[ctO[c]] = c; 1826 } 1827 *ctOrder = ctO; 1828 *ctOrderInv = ctOInv; 1829 PetscFunctionReturn(0); 1830 } 1831 1832 PetscErrorCode DMPlexCellRefinerSetUp(DMPlexCellRefiner cr) 1833 { 1834 DM dm = cr->dm; 1835 DMPolytopeType ctCell; 1836 PetscInt pStart, pEnd, p, c; 1837 PetscErrorCode ierr; 1838 1839 PetscFunctionBegin; 1840 PetscValidHeader(cr, 1); 1841 if (cr->setupcalled) PetscFunctionReturn(0); 1842 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1843 if (pEnd > pStart) {ierr = DMPlexGetCellType(dm, 0, &ctCell);CHKERRQ(ierr);} 1844 else { 1845 PetscInt dim; 1846 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1847 switch (dim) { 1848 case 0: ctCell = DM_POLYTOPE_POINT;break; 1849 case 1: ctCell = DM_POLYTOPE_SEGMENT;break; 1850 case 2: ctCell = DM_POLYTOPE_TRIANGLE;break; 1851 case 3: ctCell = DM_POLYTOPE_TETRAHEDRON;break; 1852 default: ctCell = DM_POLYTOPE_TETRAHEDRON; 1853 } 1854 } 1855 ierr = DMPlexCreateCellTypeOrder_Internal(ctCell, &cr->ctOrder, &cr->ctOrderInv);CHKERRQ(ierr); 1856 /* Construct sizes and offsets for each cell type */ 1857 if (!cr->ctStart) { 1858 PetscInt *ctS, *ctSN, *ctC, *ctCN; 1859 1860 ierr = PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);CHKERRQ(ierr); 1861 ierr = PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);CHKERRQ(ierr); 1862 for (p = pStart; p < pEnd; ++p) { 1863 DMPolytopeType ct; 1864 DMPolytopeType *rct; 1865 PetscInt *rsize, *cone, *ornt; 1866 PetscInt Nct, n; 1867 1868 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1869 if ((PetscInt) ct < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %D", p); 1870 ++ctC[ct]; 1871 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);CHKERRQ(ierr); 1872 for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n]; 1873 } 1874 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 1875 const PetscInt ct = cr->ctOrder[c]; 1876 const PetscInt ctn = cr->ctOrder[c+1]; 1877 1878 ctS[ctn] = ctS[ct] + ctC[ct]; 1879 ctSN[ctn] = ctSN[ct] + ctCN[ct]; 1880 } 1881 ierr = PetscFree2(ctC, ctCN);CHKERRQ(ierr); 1882 cr->ctStart = ctS; 1883 cr->ctStartNew = ctSN; 1884 } 1885 ierr = CellRefinerCreateOffset_Internal(cr, cr->ctOrder, cr->ctStart, &cr->offset);CHKERRQ(ierr); 1886 cr->setupcalled = PETSC_TRUE; 1887 PetscFunctionReturn(0); 1888 } 1889 1890 static PetscErrorCode DMPlexCellRefinerView_Ascii(DMPlexCellRefiner cr, PetscViewer v) 1891 { 1892 PetscErrorCode ierr; 1893 1894 PetscFunctionBegin; 1895 ierr = PetscViewerASCIIPrintf(v, "Cell Refiner: %s\n", DMPlexCellRefinerTypes[cr->type]);CHKERRQ(ierr); 1896 PetscFunctionReturn(0); 1897 } 1898 1899 /* 1900 DMPlexCellRefinerView - Views a DMPlexCellRefiner object 1901 1902 Collective on cr 1903 1904 Input Parameters: 1905 + cr - The DMPlexCellRefiner object 1906 - viewer - The PetscViewer object 1907 1908 Level: beginner 1909 1910 .seealso: DMPlexCellRefinerCreate() 1911 */ 1912 static PetscErrorCode DMPlexCellRefinerView(DMPlexCellRefiner cr, PetscViewer viewer) 1913 { 1914 PetscBool iascii; 1915 PetscErrorCode ierr; 1916 1917 PetscFunctionBegin; 1918 PetscValidHeader(cr, 1); 1919 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 1920 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) cr), &viewer);CHKERRQ(ierr);} 1921 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 1922 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1923 if (iascii) {ierr = DMPlexCellRefinerView_Ascii(cr, viewer);CHKERRQ(ierr);} 1924 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1925 PetscFunctionReturn(0); 1926 } 1927 1928 PetscErrorCode DMPlexCellRefinerDestroy(DMPlexCellRefiner *cr) 1929 { 1930 PetscInt c; 1931 PetscErrorCode ierr; 1932 1933 PetscFunctionBegin; 1934 if (!*cr) PetscFunctionReturn(0); 1935 PetscValidHeaderSpecific(*cr, DM_CLASSID, 1); 1936 ierr = PetscObjectDereference((PetscObject) (*cr)->dm);CHKERRQ(ierr); 1937 ierr = PetscFree2((*cr)->ctOrder, (*cr)->ctOrderInv);CHKERRQ(ierr); 1938 ierr = PetscFree2((*cr)->ctStart, (*cr)->ctStartNew);CHKERRQ(ierr); 1939 ierr = PetscFree((*cr)->offset);CHKERRQ(ierr); 1940 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 1941 ierr = PetscFEDestroy(&(*cr)->coordFE[c]);CHKERRQ(ierr); 1942 ierr = PetscFEGeomDestroy(&(*cr)->refGeom[c]);CHKERRQ(ierr); 1943 } 1944 ierr = PetscFree2((*cr)->coordFE, (*cr)->refGeom);CHKERRQ(ierr); 1945 ierr = PetscHeaderDestroy(cr);CHKERRQ(ierr); 1946 PetscFunctionReturn(0); 1947 } 1948 1949 PetscErrorCode DMPlexCellRefinerCreate(DM dm, DMPlexCellRefiner *cr) 1950 { 1951 DMPlexCellRefiner tmp; 1952 PetscErrorCode ierr; 1953 1954 PetscFunctionBegin; 1955 PetscValidPointer(cr, 2); 1956 *cr = NULL; 1957 ierr = PetscHeaderCreate(tmp, DM_CLASSID, "DMPlexCellRefiner", "Cell Refiner", "DMPlexCellRefiner", PETSC_COMM_SELF, DMPlexCellRefinerDestroy, DMPlexCellRefinerView);CHKERRQ(ierr); 1958 tmp->setupcalled = PETSC_FALSE; 1959 1960 tmp->dm = dm; 1961 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1962 ierr = DMPlexGetCellRefinerType(dm, &tmp->type);CHKERRQ(ierr); 1963 switch (tmp->type) { 1964 case DM_REFINER_REGULAR: 1965 tmp->ops->refine = DMPlexCellRefinerRefine_Regular; 1966 tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_Regular; 1967 tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_Regular; 1968 tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_Regular; 1969 tmp->ops->getaffinetransforms = DMPlexCellRefinerGetAffineTransforms_Regular; 1970 tmp->ops->getaffinefacetransforms = DMPlexCellRefinerGetAffineFaceTransforms_Regular; 1971 break; 1972 case DM_REFINER_TO_BOX: 1973 tmp->ops->refine = DMPlexCellRefinerRefine_ToBox; 1974 tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToBox; 1975 tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_ToBox; 1976 tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_ToBox; 1977 break; 1978 case DM_REFINER_TO_SIMPLEX: 1979 tmp->ops->refine = DMPlexCellRefinerRefine_ToSimplex; 1980 tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToSimplex; 1981 break; 1982 case DM_REFINER_ALFELD2D: 1983 tmp->ops->refine = DMPlexCellRefinerRefine_Alfeld2D; 1984 tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_None; 1985 break; 1986 case DM_REFINER_ALFELD3D: 1987 tmp->ops->refine = DMPlexCellRefinerRefine_Alfeld3D; 1988 tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_None; 1989 break; 1990 case DM_REFINER_BL: 1991 tmp->ops->refine = DMPlexCellRefinerRefine_BL; 1992 tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_BL; 1993 break; 1994 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid cell refiner type %s", DMPlexCellRefinerTypes[tmp->type]); 1995 } 1996 ierr = PetscCalloc2(DM_NUM_POLYTOPES, &tmp->coordFE, DM_NUM_POLYTOPES, &tmp->refGeom);CHKERRQ(ierr); 1997 *cr = tmp; 1998 PetscFunctionReturn(0); 1999 } 2000 2001 /*@ 2002 DMPlexCellRefinerGetAffineTransforms - Gets the affine map from the reference cell to each subcell 2003 2004 Input Parameters: 2005 + cr - The DMPlexCellRefiner object 2006 - ct - The cell type 2007 2008 Output Parameters: 2009 + Nc - The number of subcells produced from this cell type 2010 . v0 - The translation of the first vertex for each subcell 2011 . J - The Jacobian for each subcell (map from reference cell to subcell) 2012 - invJ - The inverse Jacobian for each subcell 2013 2014 Level: developer 2015 2016 .seealso: DMPlexCellRefinerGetAffineFaceTransforms(), Create() 2017 @*/ 2018 PetscErrorCode DMPlexCellRefinerGetAffineTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[]) 2019 { 2020 PetscErrorCode ierr; 2021 2022 PetscFunctionBegin; 2023 if (!cr->ops->getaffinetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine transforms from this refiner"); 2024 ierr = (*cr->ops->getaffinetransforms)(cr, ct, Nc, v0, J, invJ);CHKERRQ(ierr); 2025 PetscFunctionReturn(0); 2026 } 2027 2028 /*@ 2029 DMPlexCellRefinerGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell 2030 2031 Input Parameters: 2032 + cr - The DMPlexCellRefiner object 2033 - ct - The cell type 2034 2035 Output Parameters: 2036 + Nf - The number of faces for this cell type 2037 . v0 - The translation of the first vertex for each face 2038 . J - The Jacobian for each face (map from original cell to subcell) 2039 . invJ - The inverse Jacobian for each face 2040 - detJ - The determinant of the Jacobian for each face 2041 2042 Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse. 2043 2044 Level: developer 2045 2046 .seealso: DMPlexCellRefinerGetAffineTransforms(), Create() 2047 @*/ 2048 PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[]) 2049 { 2050 PetscErrorCode ierr; 2051 2052 PetscFunctionBegin; 2053 if (!cr->ops->getaffinefacetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine face transforms from this refiner"); 2054 ierr = (*cr->ops->getaffinefacetransforms)(cr, ct, Nf, v0, J, invJ, detJ);CHKERRQ(ierr); 2055 PetscFunctionReturn(0); 2056 } 2057 2058 /* Numbering regularly refined meshes 2059 2060 We want the numbering of the new mesh to respect the same depth stratification as the old mesh. We first compute 2061 the number of new points at each depth. This means that offsets for each depth can be computed, making no assumptions 2062 about the order of different cell types. 2063 2064 However, when we want to order different depth strata, it will be very useful to make assumptions about contiguous 2065 numbering of different cell types, especially if we want to compute new numberings without communication. Therefore, we 2066 will require that cells are numbering contiguously for each cell type, and that those blocks come in the same order as 2067 the cell type enumeration within a given depth stratum. 2068 2069 Thus, at each depth, each cell type will add a certain number of points at that depth. To get the new point number, we 2070 start at the new depth offset, run through all prior cell types incrementing by the total addition from that type, then 2071 offset by the old cell type number and replica number for the insertion. 2072 */ 2073 PetscErrorCode DMPlexCellRefinerGetNewPoint(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew) 2074 { 2075 DMPolytopeType *rct; 2076 PetscInt *rsize, *cone, *ornt; 2077 PetscInt Nct, n; 2078 PetscInt off = cr->offset[ct*DM_NUM_POLYTOPES+ctNew]; 2079 PetscInt ctS = cr->ctStart[ct], ctE = cr->ctStart[cr->ctOrder[cr->ctOrderInv[ct]+1]]; 2080 PetscInt ctSN = cr->ctStartNew[ctNew], ctEN = cr->ctStartNew[cr->ctOrder[cr->ctOrderInv[ctNew]+1]]; 2081 PetscInt newp = ctSN; 2082 PetscErrorCode ierr; 2083 2084 PetscFunctionBeginHot; 2085 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); 2086 if (off < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s does not produce type %s", DMPolytopeTypes[ct], DMPolytopeTypes[ctNew]); 2087 2088 newp += off; 2089 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);CHKERRQ(ierr); 2090 for (n = 0; n < Nct; ++n) { 2091 if (rct[n] == ctNew) { 2092 if (rsize[n] && r >= rsize[n]) 2093 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]); 2094 newp += (p - ctS) * rsize[n] + r; 2095 break; 2096 } 2097 } 2098 2099 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); 2100 *pNew = newp; 2101 PetscFunctionReturn(0); 2102 } 2103 2104 static PetscErrorCode DMPlexCellRefinerSetConeSizes(DMPlexCellRefiner cr, DM rdm) 2105 { 2106 DM dm = cr->dm; 2107 PetscInt pStart, pEnd, p, pNew; 2108 PetscErrorCode ierr; 2109 2110 PetscFunctionBegin; 2111 /* Must create the celltype label here so that we do not automatically try to compute the types */ 2112 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 2113 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2114 for (p = pStart; p < pEnd; ++p) { 2115 DMPolytopeType ct; 2116 DMPolytopeType *rct; 2117 PetscInt *rsize, *rcone, *rornt; 2118 PetscInt Nct, n, r; 2119 2120 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2121 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2122 for (n = 0; n < Nct; ++n) { 2123 for (r = 0; r < rsize[n]; ++r) { 2124 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 2125 ierr = DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));CHKERRQ(ierr); 2126 ierr = DMPlexSetCellType(rdm, pNew, rct[n]);CHKERRQ(ierr); 2127 } 2128 } 2129 } 2130 { 2131 DMLabel ctLabel; 2132 DM_Plex *plex = (DM_Plex *) rdm->data; 2133 2134 ierr = DMPlexGetCellTypeLabel(rdm, &ctLabel);CHKERRQ(ierr); 2135 ierr = PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);CHKERRQ(ierr); 2136 } 2137 PetscFunctionReturn(0); 2138 } 2139 2140 static PetscErrorCode DMPlexCellRefinerSetCones(DMPlexCellRefiner cr, DM rdm) 2141 { 2142 DM dm = cr->dm; 2143 DMPolytopeType ct; 2144 PetscInt *coneNew, *orntNew; 2145 PetscInt maxConeSize = 0, pStart, pEnd, p, pNew; 2146 PetscErrorCode ierr; 2147 2148 PetscFunctionBegin; 2149 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p)); 2150 ierr = PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew);CHKERRQ(ierr); 2151 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2152 for (p = pStart; p < pEnd; ++p) { 2153 const PetscInt *cone, *ornt; 2154 PetscInt coff, ooff, c; 2155 DMPolytopeType *rct; 2156 PetscInt *rsize, *rcone, *rornt; 2157 PetscInt Nct, n, r; 2158 2159 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2160 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 2161 ierr = DMPlexGetConeOrientation(dm, p, &ornt);CHKERRQ(ierr); 2162 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2163 for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) { 2164 const DMPolytopeType ctNew = rct[n]; 2165 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 2166 2167 for (r = 0; r < rsize[n]; ++r) { 2168 /* pNew is a subcell produced by subdividing p */ 2169 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 2170 for (c = 0; c < csizeNew; ++c) { 2171 PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */ 2172 PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */ 2173 PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */ 2174 DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */ 2175 const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */ 2176 PetscInt pr = -1; /* Replica number of pp that produces new cone point */ 2177 const DMPolytopeType ft = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */ 2178 const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 2179 PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */ 2180 PetscInt lc; 2181 2182 /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */ 2183 for (lc = 0; lc < fn; ++lc) { 2184 const PetscInt *ppornt; 2185 PetscInt pcp; 2186 2187 ierr = DMPolytopeMapCell(pct, po, rcone[coff++], &pcp);CHKERRQ(ierr); 2188 ppp = pp; 2189 pp = pcone[pcp]; 2190 ierr = DMPlexGetCellType(dm, pp, &pct);CHKERRQ(ierr); 2191 ierr = DMPlexGetCone(dm, pp, &pcone);CHKERRQ(ierr); 2192 ierr = DMPlexGetConeOrientation(dm, ppp, &ppornt);CHKERRQ(ierr); 2193 if (po < 0 && pct != DM_POLYTOPE_POINT) { 2194 const PetscInt pornt = ppornt[pcp]; 2195 const PetscInt pcsize = DMPolytopeTypeGetConeSize(pct); 2196 const PetscInt pcstart = pornt < 0 ? -(pornt+1) : pornt; 2197 const PetscInt rcstart = (pcstart+pcsize-1)%pcsize; 2198 po = pornt < 0 ? -(rcstart+1) : rcstart; 2199 } else { 2200 po = ppornt[pcp]; 2201 } 2202 } 2203 pr = rcone[coff++]; 2204 /* Orientation po of pp maps (pr, fo) -> (pr', fo') */ 2205 ierr = DMPlexCellRefinerMapSubcells(cr, pct, po, ft, pr, fo, &pr, &fo);CHKERRQ(ierr); 2206 ierr = DMPlexCellRefinerGetNewPoint(cr, pct, ft, pp, pr, &coneNew[c]);CHKERRQ(ierr); 2207 orntNew[c] = fo; 2208 } 2209 ierr = DMPlexSetCone(rdm, pNew, coneNew);CHKERRQ(ierr); 2210 ierr = DMPlexSetConeOrientation(rdm, pNew, orntNew);CHKERRQ(ierr); 2211 } 2212 } 2213 } 2214 ierr = PetscFree2(coneNew, orntNew);CHKERRQ(ierr); 2215 ierr = DMPlexSymmetrize(rdm);CHKERRQ(ierr); 2216 ierr = DMPlexStratify(rdm);CHKERRQ(ierr); 2217 PetscFunctionReturn(0); 2218 } 2219 2220 static PetscErrorCode DMPlexCellRefinerGetCoordinateFE(DMPlexCellRefiner cr, DMPolytopeType ct, PetscFE *fe) 2221 { 2222 PetscErrorCode ierr; 2223 2224 PetscFunctionBegin; 2225 if (!cr->coordFE[ct]) { 2226 PetscInt dim, cdim; 2227 PetscBool isSimplex; 2228 2229 switch (ct) { 2230 case DM_POLYTOPE_SEGMENT: dim = 1; isSimplex = PETSC_TRUE; break; 2231 case DM_POLYTOPE_TRIANGLE: dim = 2; isSimplex = PETSC_TRUE; break; 2232 case DM_POLYTOPE_QUADRILATERAL: dim = 2; isSimplex = PETSC_FALSE; break; 2233 case DM_POLYTOPE_TETRAHEDRON: dim = 3; isSimplex = PETSC_TRUE; break; 2234 case DM_POLYTOPE_HEXAHEDRON: dim = 3; isSimplex = PETSC_FALSE; break; 2235 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No coordinate FE for cell type %s", DMPolytopeTypes[ct]); 2236 } 2237 ierr = DMGetCoordinateDim(cr->dm, &cdim);CHKERRQ(ierr); 2238 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dim, cdim, isSimplex, 1, PETSC_DETERMINE, &cr->coordFE[ct]);CHKERRQ(ierr); 2239 { 2240 PetscDualSpace dsp; 2241 PetscQuadrature quad; 2242 DM K; 2243 PetscFEGeom *cg; 2244 PetscReal *Xq, *xq, *wq; 2245 PetscInt Nq, q; 2246 2247 ierr = DMPlexCellRefinerGetCellVertices(cr, ct, &Nq, &Xq);CHKERRQ(ierr); 2248 ierr = PetscMalloc1(Nq*cdim, &xq);CHKERRQ(ierr); 2249 for (q = 0; q < Nq*cdim; ++q) xq[q] = Xq[q]; 2250 ierr = PetscMalloc1(Nq, &wq);CHKERRQ(ierr); 2251 for (q = 0; q < Nq; ++q) wq[q] = 1.0; 2252 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &quad);CHKERRQ(ierr); 2253 ierr = PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);CHKERRQ(ierr); 2254 ierr = PetscFESetQuadrature(cr->coordFE[ct], quad);CHKERRQ(ierr); 2255 2256 ierr = PetscFEGetDualSpace(cr->coordFE[ct], &dsp);CHKERRQ(ierr); 2257 ierr = PetscDualSpaceGetDM(dsp, &K);CHKERRQ(ierr); 2258 ierr = PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &cr->refGeom[ct]);CHKERRQ(ierr); 2259 cg = cr->refGeom[ct]; 2260 ierr = DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);CHKERRQ(ierr); 2261 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 2262 } 2263 } 2264 *fe = cr->coordFE[ct]; 2265 PetscFunctionReturn(0); 2266 } 2267 2268 /* 2269 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct. 2270 2271 Not collective 2272 2273 Input Parameters: 2274 + cr - The DMPlexCellRefiner 2275 . ct - The type of the parent cell 2276 . rct - The type of the produced cell 2277 . r - The index of the produced cell 2278 - x - The localized coordinates for the parent cell 2279 2280 Output Parameter: 2281 . xr - The localized coordinates for the produced cell 2282 2283 Level: developer 2284 2285 .seealso: DMPlexCellRefinerSetCoordinates() 2286 */ 2287 static PetscErrorCode DMPlexCellRefinerMapLocalizedCoordinates(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[]) 2288 { 2289 PetscFE fe = NULL; 2290 PetscInt cdim, Nv, v, *subcellV; 2291 PetscErrorCode ierr; 2292 2293 PetscFunctionBegin; 2294 ierr = DMPlexCellRefinerGetCoordinateFE(cr, ct, &fe);CHKERRQ(ierr); 2295 ierr = DMPlexCellRefinerGetSubcellVertices(cr, ct, rct, r, &Nv, &subcellV);CHKERRQ(ierr); 2296 ierr = PetscFEGetNumComponents(fe, &cdim);CHKERRQ(ierr); 2297 for (v = 0; v < Nv; ++v) { 2298 ierr = PetscFEInterpolate_Static(fe, x, cr->refGeom[ct], subcellV[v], &xr[v*cdim]);CHKERRQ(ierr); 2299 } 2300 PetscFunctionReturn(0); 2301 } 2302 2303 static PetscErrorCode DMPlexCellRefinerSetCoordinates(DMPlexCellRefiner cr, DM rdm) 2304 { 2305 DM dm = cr->dm, cdm; 2306 PetscSection coordSection, coordSectionNew; 2307 Vec coordsLocal, coordsLocalNew; 2308 const PetscScalar *coords; 2309 PetscScalar *coordsNew; 2310 const DMBoundaryType *bd; 2311 const PetscReal *maxCell, *L; 2312 PetscBool isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE; 2313 PetscInt dE, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd; 2314 PetscErrorCode ierr; 2315 2316 PetscFunctionBegin; 2317 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2318 ierr = DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);CHKERRQ(ierr); 2319 /* Determine if we need to localize coordinates when generating them */ 2320 if (isperiodic) { 2321 localizeVertices = PETSC_TRUE; 2322 if (!maxCell) { 2323 PetscBool localized; 2324 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 2325 if (!localized) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized"); 2326 localizeCells = PETSC_TRUE; 2327 } 2328 } 2329 2330 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2331 ierr = PetscSectionGetFieldComponents(coordSection, 0, &dE);CHKERRQ(ierr); 2332 if (maxCell) { 2333 PetscReal maxCellNew[3]; 2334 2335 for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0; 2336 ierr = DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);CHKERRQ(ierr); 2337 } else { 2338 ierr = DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);CHKERRQ(ierr); 2339 } 2340 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);CHKERRQ(ierr); 2341 ierr = PetscSectionSetNumFields(coordSectionNew, 1);CHKERRQ(ierr); 2342 ierr = PetscSectionSetFieldComponents(coordSectionNew, 0, dE);CHKERRQ(ierr); 2343 ierr = DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);CHKERRQ(ierr); 2344 if (localizeCells) {ierr = PetscSectionSetChart(coordSectionNew, 0, vEndNew);CHKERRQ(ierr);} 2345 else {ierr = PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);CHKERRQ(ierr);} 2346 2347 /* Localization should be inherited */ 2348 /* Stefano calculates parent cells for each new cell for localization */ 2349 /* Localized cells need coordinates of closure */ 2350 for (v = vStartNew; v < vEndNew; ++v) { 2351 ierr = PetscSectionSetDof(coordSectionNew, v, dE);CHKERRQ(ierr); 2352 ierr = PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);CHKERRQ(ierr); 2353 } 2354 if (localizeCells) { 2355 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2356 for (c = cStart; c < cEnd; ++c) { 2357 PetscInt dof; 2358 2359 ierr = PetscSectionGetDof(coordSection, c, &dof); CHKERRQ(ierr); 2360 if (dof) { 2361 DMPolytopeType ct; 2362 DMPolytopeType *rct; 2363 PetscInt *rsize, *rcone, *rornt; 2364 PetscInt dim, cNew, Nct, n, r; 2365 2366 ierr = DMPlexGetCellType(dm, c, &ct);CHKERRQ(ierr); 2367 dim = DMPolytopeTypeGetDim(ct); 2368 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2369 /* This allows for different cell types */ 2370 for (n = 0; n < Nct; ++n) { 2371 if (dim != DMPolytopeTypeGetDim(rct[n])) continue; 2372 for (r = 0; r < rsize[n]; ++r) { 2373 PetscInt *closure = NULL; 2374 PetscInt clSize, cl, Nv = 0; 2375 2376 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], c, r, &cNew);CHKERRQ(ierr); 2377 ierr = DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 2378 for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;} 2379 ierr = DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 2380 ierr = PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);CHKERRQ(ierr); 2381 ierr = PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);CHKERRQ(ierr); 2382 } 2383 } 2384 } 2385 } 2386 } 2387 ierr = PetscSectionSetUp(coordSectionNew);CHKERRQ(ierr); 2388 ierr = DMViewFromOptions(dm, NULL, "-coarse_dm_view");CHKERRQ(ierr); 2389 ierr = DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);CHKERRQ(ierr); 2390 { 2391 VecType vtype; 2392 PetscInt coordSizeNew, bs; 2393 const char *name; 2394 2395 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 2396 ierr = VecCreate(PETSC_COMM_SELF, &coordsLocalNew);CHKERRQ(ierr); 2397 ierr = PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);CHKERRQ(ierr); 2398 ierr = VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);CHKERRQ(ierr); 2399 ierr = PetscObjectGetName((PetscObject) coordsLocal, &name);CHKERRQ(ierr); 2400 ierr = PetscObjectSetName((PetscObject) coordsLocalNew, name);CHKERRQ(ierr); 2401 ierr = VecGetBlockSize(coordsLocal, &bs);CHKERRQ(ierr); 2402 ierr = VecSetBlockSize(coordsLocalNew, bs);CHKERRQ(ierr); 2403 ierr = VecGetType(coordsLocal, &vtype);CHKERRQ(ierr); 2404 ierr = VecSetType(coordsLocalNew, vtype);CHKERRQ(ierr); 2405 } 2406 ierr = VecGetArrayRead(coordsLocal, &coords);CHKERRQ(ierr); 2407 ierr = VecGetArray(coordsLocalNew, &coordsNew);CHKERRQ(ierr); 2408 ierr = PetscSectionGetChart(coordSection, &ocStart, &ocEnd);CHKERRQ(ierr); 2409 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2410 /* First set coordinates for vertices*/ 2411 for (p = pStart; p < pEnd; ++p) { 2412 DMPolytopeType ct; 2413 DMPolytopeType *rct; 2414 PetscInt *rsize, *rcone, *rornt; 2415 PetscInt Nct, n, r; 2416 PetscBool hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE; 2417 2418 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2419 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2420 for (n = 0; n < Nct; ++n) { 2421 if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;} 2422 } 2423 if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) { 2424 PetscInt dof; 2425 ierr = PetscSectionGetDof(coordSection, p, &dof); CHKERRQ(ierr); 2426 if (dof) isLocalized = PETSC_TRUE; 2427 } 2428 if (hasVertex) { 2429 PetscScalar *pcoords = NULL; 2430 PetscScalar vcoords[3] = {0., 0., 0.}; 2431 PetscInt Nc, Nv, v, d; 2432 2433 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);CHKERRQ(ierr); 2434 if (ct == DM_POLYTOPE_POINT) { 2435 for (d = 0; d < dE; ++d) vcoords[d] = pcoords[d]; 2436 } else { 2437 if (localizeVertices) { 2438 PetscScalar anchor[3]; 2439 2440 for (d = 0; d < dE; ++d) anchor[d] = pcoords[d]; 2441 if (!isLocalized) { 2442 Nv = Nc/dE; 2443 for (v = 0; v < Nv; ++v) {ierr = DMLocalizeAddCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], vcoords);CHKERRQ(ierr);} 2444 } else { 2445 Nv = Nc/(2*dE); 2446 for (v = Nv; v < Nv*2; ++v) {ierr = DMLocalizeAddCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], vcoords);CHKERRQ(ierr);} 2447 } 2448 } else { 2449 Nv = Nc/dE; 2450 for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) vcoords[d] += pcoords[v*dE+d]; 2451 } 2452 for (d = 0; d < dE; ++d) vcoords[d] /= Nv; 2453 } 2454 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);CHKERRQ(ierr); 2455 for (n = 0; n < Nct; ++n) { 2456 if (rct[n] != DM_POLYTOPE_POINT) continue; 2457 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]); 2458 for (r = 0; r < rsize[n]; ++r) { 2459 PetscInt vNew, off; 2460 2461 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &vNew);CHKERRQ(ierr); 2462 ierr = PetscSectionGetOffset(coordSectionNew, vNew, &off);CHKERRQ(ierr); 2463 ierr = DMPlexSnapToGeomModel(dm, p, vcoords, &coordsNew[off]);CHKERRQ(ierr); 2464 } 2465 } 2466 } 2467 } 2468 /* Then set coordinates for cells by localizing */ 2469 for (p = pStart; p < pEnd; ++p) { 2470 DMPolytopeType ct; 2471 DMPolytopeType *rct; 2472 PetscInt *rsize, *rcone, *rornt; 2473 PetscInt Nct, n, r; 2474 PetscBool isLocalized = PETSC_FALSE; 2475 2476 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2477 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2478 if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) { 2479 PetscInt dof; 2480 ierr = PetscSectionGetDof(coordSection, p, &dof); CHKERRQ(ierr); 2481 if (dof) isLocalized = PETSC_TRUE; 2482 } 2483 if (isLocalized) { 2484 const PetscScalar *pcoords; 2485 2486 ierr = DMPlexPointLocalRead(cdm, p, coords, &pcoords);CHKERRQ(ierr); 2487 for (n = 0; n < Nct; ++n) { 2488 const PetscInt Nr = rsize[n]; 2489 2490 if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue; 2491 for (r = 0; r < Nr; ++r) { 2492 PetscInt pNew, offNew; 2493 2494 /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that 2495 DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger 2496 cell to the ones it produces. */ 2497 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 2498 ierr = PetscSectionGetOffset(coordSectionNew, pNew, &offNew);CHKERRQ(ierr); 2499 ierr = DMPlexCellRefinerMapLocalizedCoordinates(cr, ct, rct[n], r, pcoords, &coordsNew[offNew]);CHKERRQ(ierr); 2500 } 2501 } 2502 } 2503 } 2504 ierr = VecRestoreArrayRead(coordsLocal, &coords);CHKERRQ(ierr); 2505 ierr = VecRestoreArray(coordsLocalNew, &coordsNew);CHKERRQ(ierr); 2506 ierr = DMSetCoordinatesLocal(rdm, coordsLocalNew);CHKERRQ(ierr); 2507 /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */ 2508 ierr = VecDestroy(&coordsLocalNew);CHKERRQ(ierr); 2509 ierr = PetscSectionDestroy(&coordSectionNew);CHKERRQ(ierr); 2510 if (!localizeCells) {ierr = DMLocalizeCoordinates(rdm);CHKERRQ(ierr);} 2511 PetscFunctionReturn(0); 2512 } 2513 2514 /*@ 2515 DMPlexCreateProcessSF - Create an SF which just has process connectivity 2516 2517 Collective on dm 2518 2519 Input Parameters: 2520 + dm - The DM 2521 - sfPoint - The PetscSF which encodes point connectivity 2522 2523 Output Parameters: 2524 + processRanks - A list of process neighbors, or NULL 2525 - sfProcess - An SF encoding the process connectivity, or NULL 2526 2527 Level: developer 2528 2529 .seealso: PetscSFCreate(), DMPlexCreateTwoSidedProcessSF() 2530 @*/ 2531 PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess) 2532 { 2533 PetscInt numRoots, numLeaves, l; 2534 const PetscInt *localPoints; 2535 const PetscSFNode *remotePoints; 2536 PetscInt *localPointsNew; 2537 PetscSFNode *remotePointsNew; 2538 PetscInt *ranks, *ranksNew; 2539 PetscMPIInt size; 2540 PetscErrorCode ierr; 2541 2542 PetscFunctionBegin; 2543 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2544 PetscValidHeaderSpecific(sfPoint, PETSCSF_CLASSID, 2); 2545 if (processRanks) {PetscValidPointer(processRanks, 3);} 2546 if (sfProcess) {PetscValidPointer(sfProcess, 4);} 2547 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2548 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2549 ierr = PetscMalloc1(numLeaves, &ranks);CHKERRQ(ierr); 2550 for (l = 0; l < numLeaves; ++l) { 2551 ranks[l] = remotePoints[l].rank; 2552 } 2553 ierr = PetscSortRemoveDupsInt(&numLeaves, ranks);CHKERRQ(ierr); 2554 ierr = PetscMalloc1(numLeaves, &ranksNew);CHKERRQ(ierr); 2555 ierr = PetscMalloc1(numLeaves, &localPointsNew);CHKERRQ(ierr); 2556 ierr = PetscMalloc1(numLeaves, &remotePointsNew);CHKERRQ(ierr); 2557 for (l = 0; l < numLeaves; ++l) { 2558 ranksNew[l] = ranks[l]; 2559 localPointsNew[l] = l; 2560 remotePointsNew[l].index = 0; 2561 remotePointsNew[l].rank = ranksNew[l]; 2562 } 2563 ierr = PetscFree(ranks);CHKERRQ(ierr); 2564 if (processRanks) {ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);CHKERRQ(ierr);} 2565 else {ierr = PetscFree(ranksNew);CHKERRQ(ierr);} 2566 if (sfProcess) { 2567 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);CHKERRQ(ierr); 2568 ierr = PetscObjectSetName((PetscObject) *sfProcess, "Process SF");CHKERRQ(ierr); 2569 ierr = PetscSFSetFromOptions(*sfProcess);CHKERRQ(ierr); 2570 ierr = PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2571 } 2572 PetscFunctionReturn(0); 2573 } 2574 2575 static PetscErrorCode DMPlexCellRefinerCreateSF(DMPlexCellRefiner cr, DM rdm) 2576 { 2577 DM dm = cr->dm; 2578 DMPlexCellRefiner *crRem; 2579 PetscSF sf, sfNew, sfProcess; 2580 IS processRanks; 2581 MPI_Datatype ctType; 2582 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 2583 const PetscInt *localPoints, *neighbors; 2584 const PetscSFNode *remotePoints; 2585 PetscInt *localPointsNew; 2586 PetscSFNode *remotePointsNew; 2587 PetscInt *ctStartRem, *ctStartNewRem; 2588 PetscInt ctSize = DM_NUM_POLYTOPES+1, numNeighbors, n, pStartNew, pEndNew, pNew, pNewRem; 2589 PetscErrorCode ierr; 2590 2591 PetscFunctionBegin; 2592 ierr = DMPlexGetChart(rdm, &pStartNew, &pEndNew);CHKERRQ(ierr); 2593 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2594 ierr = DMGetPointSF(rdm, &sfNew);CHKERRQ(ierr); 2595 /* Calculate size of new SF */ 2596 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2597 if (numRoots < 0) PetscFunctionReturn(0); 2598 for (l = 0; l < numLeaves; ++l) { 2599 const PetscInt p = localPoints[l]; 2600 DMPolytopeType ct; 2601 DMPolytopeType *rct; 2602 PetscInt *rsize, *rcone, *rornt; 2603 PetscInt Nct, n; 2604 2605 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2606 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2607 for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n]; 2608 } 2609 /* Communicate ctStart and cStartNew for each remote rank */ 2610 ierr = DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);CHKERRQ(ierr); 2611 ierr = ISGetLocalSize(processRanks, &numNeighbors);CHKERRQ(ierr); 2612 ierr = PetscMalloc2(ctSize*numNeighbors, &ctStartRem, ctSize*numNeighbors, &ctStartNewRem);CHKERRQ(ierr); 2613 ierr = MPI_Type_contiguous(ctSize, MPIU_INT, &ctType);CHKERRQ(ierr); 2614 ierr = MPI_Type_commit(&ctType);CHKERRQ(ierr); 2615 ierr = PetscSFBcastBegin(sfProcess, ctType, cr->ctStart, ctStartRem);CHKERRQ(ierr); 2616 ierr = PetscSFBcastEnd(sfProcess, ctType, cr->ctStart, ctStartRem);CHKERRQ(ierr); 2617 ierr = PetscSFBcastBegin(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);CHKERRQ(ierr); 2618 ierr = PetscSFBcastEnd(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);CHKERRQ(ierr); 2619 ierr = MPI_Type_free(&ctType);CHKERRQ(ierr); 2620 ierr = PetscSFDestroy(&sfProcess);CHKERRQ(ierr); 2621 ierr = PetscMalloc1(numNeighbors, &crRem);CHKERRQ(ierr); 2622 for (n = 0; n < numNeighbors; ++n) { 2623 ierr = DMPlexCellRefinerCreate(dm, &crRem[n]);CHKERRQ(ierr); 2624 ierr = DMPlexCellRefinerSetStarts(crRem[n], &ctStartRem[n*ctSize], &ctStartNewRem[n*ctSize]); 2625 ierr = DMPlexCellRefinerSetUp(crRem[n]);CHKERRQ(ierr); 2626 } 2627 ierr = PetscFree2(ctStartRem, ctStartNewRem);CHKERRQ(ierr); 2628 /* Calculate new point SF */ 2629 ierr = PetscMalloc1(numLeavesNew, &localPointsNew);CHKERRQ(ierr); 2630 ierr = PetscMalloc1(numLeavesNew, &remotePointsNew);CHKERRQ(ierr); 2631 ierr = ISGetIndices(processRanks, &neighbors);CHKERRQ(ierr); 2632 for (l = 0, m = 0; l < numLeaves; ++l) { 2633 PetscInt p = localPoints[l]; 2634 PetscInt pRem = remotePoints[l].index; 2635 PetscMPIInt rankRem = remotePoints[l].rank; 2636 DMPolytopeType ct; 2637 DMPolytopeType *rct; 2638 PetscInt *rsize, *rcone, *rornt; 2639 PetscInt neighbor, Nct, n, r; 2640 2641 ierr = PetscFindInt(rankRem, numNeighbors, neighbors, &neighbor);CHKERRQ(ierr); 2642 if (neighbor < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %D", rankRem); 2643 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 2644 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2645 for (n = 0; n < Nct; ++n) { 2646 for (r = 0; r < rsize[n]; ++r) { 2647 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 2648 ierr = DMPlexCellRefinerGetNewPoint(crRem[neighbor], ct, rct[n], pRem, r, &pNewRem);CHKERRQ(ierr); 2649 localPointsNew[m] = pNew; 2650 remotePointsNew[m].index = pNewRem; 2651 remotePointsNew[m].rank = rankRem; 2652 ++m; 2653 } 2654 } 2655 } 2656 for (n = 0; n < numNeighbors; ++n) {ierr = DMPlexCellRefinerDestroy(&crRem[n]);CHKERRQ(ierr);} 2657 ierr = PetscFree(crRem);CHKERRQ(ierr); 2658 if (m != numLeavesNew) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of leaf point %D should be %D", m, numLeavesNew); 2659 ierr = ISRestoreIndices(processRanks, &neighbors);CHKERRQ(ierr); 2660 ierr = ISDestroy(&processRanks);CHKERRQ(ierr); 2661 { 2662 PetscSFNode *rp, *rtmp; 2663 PetscInt *lp, *idx, *ltmp, i; 2664 2665 /* SF needs sorted leaves to correct calculate Gather */ 2666 ierr = PetscMalloc1(numLeavesNew, &idx);CHKERRQ(ierr); 2667 ierr = PetscMalloc1(numLeavesNew, &lp);CHKERRQ(ierr); 2668 ierr = PetscMalloc1(numLeavesNew, &rp);CHKERRQ(ierr); 2669 for (i = 0; i < numLeavesNew; ++i) { 2670 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); 2671 idx[i] = i; 2672 } 2673 ierr = PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);CHKERRQ(ierr); 2674 for (i = 0; i < numLeavesNew; ++i) { 2675 lp[i] = localPointsNew[idx[i]]; 2676 rp[i] = remotePointsNew[idx[i]]; 2677 } 2678 ltmp = localPointsNew; 2679 localPointsNew = lp; 2680 rtmp = remotePointsNew; 2681 remotePointsNew = rp; 2682 ierr = PetscFree(idx);CHKERRQ(ierr); 2683 ierr = PetscFree(ltmp);CHKERRQ(ierr); 2684 ierr = PetscFree(rtmp);CHKERRQ(ierr); 2685 } 2686 ierr = PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2687 PetscFunctionReturn(0); 2688 } 2689 2690 static PetscErrorCode RefineLabel_Internal(DMPlexCellRefiner cr, DMLabel label, DMLabel labelNew) 2691 { 2692 DM dm = cr->dm; 2693 IS valueIS; 2694 const PetscInt *values; 2695 PetscInt defVal, Nv, val; 2696 PetscErrorCode ierr; 2697 2698 PetscFunctionBegin; 2699 ierr = DMLabelGetDefaultValue(label, &defVal);CHKERRQ(ierr); 2700 ierr = DMLabelSetDefaultValue(labelNew, defVal);CHKERRQ(ierr); 2701 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 2702 ierr = ISGetLocalSize(valueIS, &Nv);CHKERRQ(ierr); 2703 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 2704 for (val = 0; val < Nv; ++val) { 2705 IS pointIS; 2706 const PetscInt *points; 2707 PetscInt numPoints, p; 2708 2709 /* Ensure refined label is created with same number of strata as 2710 * original (even if no entries here). */ 2711 ierr = DMLabelAddStratum(labelNew, values[val]);CHKERRQ(ierr); 2712 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 2713 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 2714 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2715 for (p = 0; p < numPoints; ++p) { 2716 const PetscInt point = points[p]; 2717 DMPolytopeType ct; 2718 DMPolytopeType *rct; 2719 PetscInt *rsize, *rcone, *rornt; 2720 PetscInt Nct, n, r, pNew; 2721 2722 ierr = DMPlexGetCellType(dm, point, &ct);CHKERRQ(ierr); 2723 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 2724 for (n = 0; n < Nct; ++n) { 2725 for (r = 0; r < rsize[n]; ++r) { 2726 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], point, r, &pNew);CHKERRQ(ierr); 2727 ierr = DMLabelSetValue(labelNew, pNew, values[val]);CHKERRQ(ierr); 2728 } 2729 } 2730 } 2731 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2732 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2733 } 2734 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 2735 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 2736 PetscFunctionReturn(0); 2737 } 2738 2739 static PetscErrorCode DMPlexCellRefinerCreateLabels(DMPlexCellRefiner cr, DM rdm) 2740 { 2741 DM dm = cr->dm; 2742 PetscInt numLabels, l; 2743 PetscErrorCode ierr; 2744 2745 PetscFunctionBegin; 2746 ierr = DMGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 2747 for (l = 0; l < numLabels; ++l) { 2748 DMLabel label, labelNew; 2749 const char *lname; 2750 PetscBool isDepth, isCellType; 2751 2752 ierr = DMGetLabelName(dm, l, &lname);CHKERRQ(ierr); 2753 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 2754 if (isDepth) continue; 2755 ierr = PetscStrcmp(lname, "celltype", &isCellType);CHKERRQ(ierr); 2756 if (isCellType) continue; 2757 ierr = DMCreateLabel(rdm, lname);CHKERRQ(ierr); 2758 ierr = DMGetLabel(dm, lname, &label);CHKERRQ(ierr); 2759 ierr = DMGetLabel(rdm, lname, &labelNew);CHKERRQ(ierr); 2760 ierr = RefineLabel_Internal(cr, label, labelNew);CHKERRQ(ierr); 2761 } 2762 PetscFunctionReturn(0); 2763 } 2764 2765 /* This will only work for interpolated meshes */ 2766 PetscErrorCode DMPlexRefineUniform(DM dm, DMPlexCellRefiner cr, DM *dmRefined) 2767 { 2768 DM rdm; 2769 PetscInt dim, embedDim, depth; 2770 PetscErrorCode ierr; 2771 2772 PetscFunctionBegin; 2773 PetscValidHeader(cr, 1); 2774 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &rdm);CHKERRQ(ierr); 2775 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2776 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2777 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2778 ierr = DMGetCoordinateDim(dm, &embedDim);CHKERRQ(ierr); 2779 ierr = DMSetCoordinateDim(rdm, embedDim);CHKERRQ(ierr); 2780 /* Calculate number of new points of each depth */ 2781 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2782 if (depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated for regular refinement"); 2783 /* Step 1: Set chart */ 2784 ierr = DMPlexSetChart(rdm, 0, cr->ctStartNew[cr->ctOrder[DM_NUM_POLYTOPES]]);CHKERRQ(ierr); 2785 /* Step 2: Set cone/support sizes (automatically stratifies) */ 2786 ierr = DMPlexCellRefinerSetConeSizes(cr, rdm);CHKERRQ(ierr); 2787 /* Step 3: Setup refined DM */ 2788 ierr = DMSetUp(rdm);CHKERRQ(ierr); 2789 /* Step 4: Set cones and supports (automatically symmetrizes) */ 2790 ierr = DMPlexCellRefinerSetCones(cr, rdm);CHKERRQ(ierr); 2791 /* Step 5: Create pointSF */ 2792 ierr = DMPlexCellRefinerCreateSF(cr, rdm);CHKERRQ(ierr); 2793 /* Step 6: Create labels */ 2794 ierr = DMPlexCellRefinerCreateLabels(cr, rdm);CHKERRQ(ierr); 2795 /* Step 7: Set coordinates */ 2796 ierr = DMPlexCellRefinerSetCoordinates(cr, rdm);CHKERRQ(ierr); 2797 *dmRefined = rdm; 2798 PetscFunctionReturn(0); 2799 } 2800 2801 /*@ 2802 DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data 2803 2804 Input Parameter: 2805 . dm - The coarse DM 2806 2807 Output Parameter: 2808 . fpointIS - The IS of all the fine points which exist in the original coarse mesh 2809 2810 Level: developer 2811 2812 .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS() 2813 @*/ 2814 PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS) 2815 { 2816 DMPlexCellRefiner cr; 2817 PetscInt *fpoints; 2818 PetscInt pStart, pEnd, p, vStart, vEnd, v; 2819 PetscErrorCode ierr; 2820 2821 PetscFunctionBegin; 2822 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2823 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2824 ierr = DMPlexCellRefinerCreate(dm, &cr);CHKERRQ(ierr); 2825 ierr = DMPlexCellRefinerSetUp(cr);CHKERRQ(ierr); 2826 ierr = PetscMalloc1(pEnd-pStart, &fpoints);CHKERRQ(ierr); 2827 for (p = 0; p < pEnd-pStart; ++p) fpoints[p] = -1; 2828 for (v = vStart; v < vEnd; ++v) { 2829 PetscInt vNew = -1; /* silent overzelous may be used uninitialized */ 2830 2831 ierr = DMPlexCellRefinerGetNewPoint(cr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);CHKERRQ(ierr); 2832 fpoints[v-pStart] = vNew; 2833 } 2834 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr); 2835 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);CHKERRQ(ierr); 2836 PetscFunctionReturn(0); 2837 } 2838 2839 /*@ 2840 DMPlexSetRefinementUniform - Set the flag for uniform refinement 2841 2842 Input Parameters: 2843 + dm - The DM 2844 - refinementUniform - The flag for uniform refinement 2845 2846 Level: developer 2847 2848 .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 2849 @*/ 2850 PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform) 2851 { 2852 DM_Plex *mesh = (DM_Plex*) dm->data; 2853 2854 PetscFunctionBegin; 2855 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2856 mesh->refinementUniform = refinementUniform; 2857 PetscFunctionReturn(0); 2858 } 2859 2860 /*@ 2861 DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement 2862 2863 Input Parameter: 2864 . dm - The DM 2865 2866 Output Parameter: 2867 . refinementUniform - The flag for uniform refinement 2868 2869 Level: developer 2870 2871 .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 2872 @*/ 2873 PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform) 2874 { 2875 DM_Plex *mesh = (DM_Plex*) dm->data; 2876 2877 PetscFunctionBegin; 2878 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2879 PetscValidPointer(refinementUniform, 2); 2880 *refinementUniform = mesh->refinementUniform; 2881 PetscFunctionReturn(0); 2882 } 2883 2884 /*@ 2885 DMPlexSetRefinementLimit - Set the maximum cell volume for refinement 2886 2887 Input Parameters: 2888 + dm - The DM 2889 - refinementLimit - The maximum cell volume in the refined mesh 2890 2891 Level: developer 2892 2893 .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform() 2894 @*/ 2895 PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit) 2896 { 2897 DM_Plex *mesh = (DM_Plex*) dm->data; 2898 2899 PetscFunctionBegin; 2900 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2901 mesh->refinementLimit = refinementLimit; 2902 PetscFunctionReturn(0); 2903 } 2904 2905 /*@ 2906 DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement 2907 2908 Input Parameter: 2909 . dm - The DM 2910 2911 Output Parameter: 2912 . refinementLimit - The maximum cell volume in the refined mesh 2913 2914 Level: developer 2915 2916 .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform() 2917 @*/ 2918 PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit) 2919 { 2920 DM_Plex *mesh = (DM_Plex*) dm->data; 2921 2922 PetscFunctionBegin; 2923 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2924 PetscValidPointer(refinementLimit, 2); 2925 /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */ 2926 *refinementLimit = mesh->refinementLimit; 2927 PetscFunctionReturn(0); 2928 } 2929 2930 /*@ 2931 DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement 2932 2933 Input Parameters: 2934 + dm - The DM 2935 - refinementFunc - Function giving the maximum cell volume in the refined mesh 2936 2937 Note: The calling sequence is refinementFunc(coords, limit) 2938 $ coords - Coordinates of the current point, usually a cell centroid 2939 $ limit - The maximum cell volume for a cell containing this point 2940 2941 Level: developer 2942 2943 .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 2944 @*/ 2945 PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *)) 2946 { 2947 DM_Plex *mesh = (DM_Plex*) dm->data; 2948 2949 PetscFunctionBegin; 2950 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2951 mesh->refinementFunc = refinementFunc; 2952 PetscFunctionReturn(0); 2953 } 2954 2955 /*@ 2956 DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement 2957 2958 Input Parameter: 2959 . dm - The DM 2960 2961 Output Parameter: 2962 . refinementFunc - Function giving the maximum cell volume in the refined mesh 2963 2964 Note: The calling sequence is refinementFunc(coords, limit) 2965 $ coords - Coordinates of the current point, usually a cell centroid 2966 $ limit - The maximum cell volume for a cell containing this point 2967 2968 Level: developer 2969 2970 .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 2971 @*/ 2972 PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *)) 2973 { 2974 DM_Plex *mesh = (DM_Plex*) dm->data; 2975 2976 PetscFunctionBegin; 2977 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2978 PetscValidPointer(refinementFunc, 2); 2979 *refinementFunc = mesh->refinementFunc; 2980 PetscFunctionReturn(0); 2981 } 2982 2983 static PetscErrorCode RefineDiscLabels_Internal(DMPlexCellRefiner cr, DM rdm) 2984 { 2985 DM dm = cr->dm; 2986 PetscInt Nf, f, Nds, s; 2987 PetscErrorCode ierr; 2988 2989 PetscFunctionBegin; 2990 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 2991 for (f = 0; f < Nf; ++f) { 2992 DMLabel label, labelNew; 2993 PetscObject obj; 2994 const char *lname; 2995 2996 ierr = DMGetField(rdm, f, &label, &obj);CHKERRQ(ierr); 2997 if (!label) continue; 2998 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr); 2999 ierr = DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);CHKERRQ(ierr); 3000 ierr = RefineLabel_Internal(cr, label, labelNew);CHKERRQ(ierr); 3001 ierr = DMSetField_Internal(rdm, f, labelNew, obj);CHKERRQ(ierr); 3002 ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr); 3003 } 3004 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr); 3005 for (s = 0; s < Nds; ++s) { 3006 DMLabel label, labelNew; 3007 const char *lname; 3008 3009 ierr = DMGetRegionNumDS(rdm, s, &label, NULL, NULL);CHKERRQ(ierr); 3010 if (!label) continue; 3011 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr); 3012 ierr = DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);CHKERRQ(ierr); 3013 ierr = RefineLabel_Internal(cr, label, labelNew);CHKERRQ(ierr); 3014 ierr = DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);CHKERRQ(ierr); 3015 ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr); 3016 } 3017 PetscFunctionReturn(0); 3018 } 3019 3020 PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined) 3021 { 3022 PetscBool isUniform; 3023 DMPlexCellRefiner cr; 3024 PetscErrorCode ierr; 3025 3026 PetscFunctionBegin; 3027 ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr); 3028 ierr = DMViewFromOptions(dm, NULL, "-initref_dm_view");CHKERRQ(ierr); 3029 if (isUniform) { 3030 PetscBool localized; 3031 3032 ierr = DMPlexCellRefinerCreate(dm, &cr);CHKERRQ(ierr); 3033 ierr = DMPlexCellRefinerSetUp(cr);CHKERRQ(ierr); 3034 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 3035 ierr = DMPlexRefineUniform(dm, cr, dmRefined);CHKERRQ(ierr); 3036 ierr = DMPlexSetRegularRefinement(*dmRefined, PETSC_TRUE);CHKERRQ(ierr); 3037 ierr = DMCopyDisc(dm, *dmRefined);CHKERRQ(ierr); 3038 ierr = RefineDiscLabels_Internal(cr, *dmRefined);CHKERRQ(ierr); 3039 ierr = DMCopyBoundary(dm, *dmRefined);CHKERRQ(ierr); 3040 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr); 3041 } else { 3042 ierr = DMPlexRefine_Internal(dm, NULL, dmRefined);CHKERRQ(ierr); 3043 } 3044 PetscFunctionReturn(0); 3045 } 3046 3047 PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]) 3048 { 3049 DM cdm = dm; 3050 PetscInt r; 3051 PetscBool isUniform, localized; 3052 PetscErrorCode ierr; 3053 3054 PetscFunctionBegin; 3055 ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr); 3056 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 3057 if (isUniform) { 3058 for (r = 0; r < nlevels; ++r) { 3059 DMPlexCellRefiner cr; 3060 3061 ierr = DMPlexCellRefinerCreate(cdm, &cr);CHKERRQ(ierr); 3062 ierr = DMPlexCellRefinerSetUp(cr);CHKERRQ(ierr); 3063 ierr = DMPlexRefineUniform(cdm, cr, &dmRefined[r]);CHKERRQ(ierr); 3064 ierr = DMSetCoarsenLevel(dmRefined[r], cdm->leveldown);CHKERRQ(ierr); 3065 ierr = DMSetRefineLevel(dmRefined[r], cdm->levelup+1);CHKERRQ(ierr); 3066 ierr = DMCopyDisc(cdm, dmRefined[r]);CHKERRQ(ierr); 3067 ierr = RefineDiscLabels_Internal(cr, dmRefined[r]);CHKERRQ(ierr); 3068 ierr = DMCopyBoundary(cdm, dmRefined[r]);CHKERRQ(ierr); 3069 ierr = DMSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr); 3070 ierr = DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);CHKERRQ(ierr); 3071 cdm = dmRefined[r]; 3072 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr); 3073 } 3074 } else { 3075 for (r = 0; r < nlevels; ++r) { 3076 ierr = DMRefine(cdm, PetscObjectComm((PetscObject) dm), &dmRefined[r]);CHKERRQ(ierr); 3077 ierr = DMCopyDisc(cdm, dmRefined[r]);CHKERRQ(ierr); 3078 ierr = DMCopyBoundary(cdm, dmRefined[r]);CHKERRQ(ierr); 3079 if (localized) {ierr = DMLocalizeCoordinates(dmRefined[r]);CHKERRQ(ierr);} 3080 ierr = DMSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr); 3081 cdm = dmRefined[r]; 3082 } 3083 } 3084 PetscFunctionReturn(0); 3085 } 3086