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