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