1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "GetDepthStart_Private" 6 PETSC_STATIC_INLINE PetscErrorCode GetDepthStart_Private(PetscInt depth, PetscInt depthSize[], PetscInt *cStart, PetscInt *fStart, PetscInt *eStart, PetscInt *vStart) 7 { 8 PetscFunctionBegin; 9 if (cStart) *cStart = 0; 10 if (vStart) *vStart = depthSize[depth]; 11 if (fStart) *fStart = depthSize[depth] + depthSize[0]; 12 if (eStart) *eStart = depthSize[depth] + depthSize[0] + depthSize[depth-1]; 13 PetscFunctionReturn(0); 14 } 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "GetDepthEnd_Private" 18 PETSC_STATIC_INLINE PetscErrorCode GetDepthEnd_Private(PetscInt depth, PetscInt depthSize[], PetscInt *cEnd, PetscInt *fEnd, PetscInt *eEnd, PetscInt *vEnd) 19 { 20 PetscFunctionBegin; 21 if (cEnd) *cEnd = depthSize[depth]; 22 if (vEnd) *vEnd = depthSize[depth] + depthSize[0]; 23 if (fEnd) *fEnd = depthSize[depth] + depthSize[0] + depthSize[depth-1]; 24 if (eEnd) *eEnd = depthSize[depth] + depthSize[0] + depthSize[depth-1] + depthSize[1]; 25 PetscFunctionReturn(0); 26 } 27 28 #undef __FUNCT__ 29 #define __FUNCT__ "DMPlexGetNumHybridFaces_Internal" 30 /* This is a stopgap since we do not currently keep track of faces for hybrid cells */ 31 static PetscErrorCode DMPlexGetNumHybridFaces_Internal(DM dm, PetscInt *numHybridFaces) 32 { 33 PetscInt eStart, eEnd, eMax, cEnd, cMax, c, hEdges = 0; 34 PetscErrorCode ierr; 35 36 PetscFunctionBegin; 37 *numHybridFaces = 0; 38 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 39 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 40 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, &eMax, NULL);CHKERRQ(ierr); 41 if (cMax < 0) PetscFunctionReturn(0); 42 /* Count interior edges in hybrid cells */ 43 for (c = cMax; c < cEnd; ++c) { 44 PetscInt *closure = NULL, closureSize, cl; 45 46 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 47 for (cl = 0; cl < closureSize*2; cl += 2) { 48 const PetscInt p = closure[cl]; 49 50 if ((p >= eStart) && (p < eMax)) ++hEdges; 51 } 52 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 53 } 54 if (hEdges%2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of interior edges in hybrid cells cannot be odd: %d", hEdges); 55 *numHybridFaces = hEdges/2; 56 PetscFunctionReturn(0); 57 } 58 59 #undef __FUNCT__ 60 #define __FUNCT__ "CellRefinerGetSizes" 61 PetscErrorCode CellRefinerGetSizes(CellRefiner refiner, DM dm, PetscInt depthSize[]) 62 { 63 PetscInt numHybridFaces, cStart, cEnd, cMax, vStart, vEnd, vMax, fStart, fEnd, fMax, eStart, eEnd, eMax; 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 68 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 69 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 70 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 71 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 72 switch (refiner) { 73 case 1: 74 /* Simplicial 2D */ 75 depthSize[0] = vEnd - vStart + fEnd - fStart; /* Add a vertex on every face */ 76 depthSize[1] = 2*(fEnd - fStart) + 3*(cEnd - cStart); /* Every face is split into 2 faces and 3 faces are added for each cell */ 77 depthSize[2] = 4*(cEnd - cStart); /* Every cell split into 4 cells */ 78 break; 79 case 3: 80 /* Hybrid Simplicial 2D */ 81 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 82 cMax = PetscMin(cEnd, cMax); 83 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 84 fMax = PetscMin(fEnd, fMax); 85 depthSize[0] = vEnd - vStart + fMax - fStart; /* Add a vertex on every face, but not hybrid faces */ 86 depthSize[1] = 2*(fMax - fStart) + 3*(cMax - cStart) + (fEnd - fMax) + (cEnd - cMax); /* Every interior face is split into 2 faces, 3 faces are added for each interior cell, and one in each hybrid cell */ 87 depthSize[2] = 4*(cMax - cStart) + 2*(cEnd - cMax); /* Interior cells split into 4 cells, Hybrid cells split into 2 cells */ 88 break; 89 case 2: 90 /* Hex 2D */ 91 depthSize[0] = vEnd - vStart + cEnd - cStart + fEnd - fStart; /* Add a vertex on every face and cell */ 92 depthSize[1] = 2*(fEnd - fStart) + 4*(cEnd - cStart); /* Every face is split into 2 faces and 4 faces are added for each cell */ 93 depthSize[2] = 4*(cEnd - cStart); /* Every cell split into 4 cells */ 94 break; 95 case 5: 96 /* Simplicial 3D */ 97 depthSize[0] = vEnd - vStart + eEnd - eStart; /* Add a vertex on every edge */ 98 depthSize[1] = 2*(eEnd - eStart) + 3*(fEnd - fStart) + (cEnd - cStart); /* Every edge is split into 2 edges, 3 edges are added for each face, and 1 edge for each cell */ 99 depthSize[2] = 4*(fEnd - fStart) + 8*(cEnd - cStart); /* Every face split into 4 faces and 8 faces are added for each cell */ 100 depthSize[3] = 8*(cEnd - cStart); /* Every cell split into 8 cells */ 101 break; 102 case 6: 103 /* Hex 3D */ 104 depthSize[0] = vEnd - vStart + eEnd - eStart + fEnd - fStart + cEnd - cStart; /* Add a vertex on every edge, face and cell */ 105 depthSize[1] = 2*(eEnd - eStart) + 4*(fEnd - fStart) + 6*(cEnd - cStart); /* Every edge is split into 2 edge, 4 edges are added for each face, and 6 edges for each cell */ 106 depthSize[2] = 4*(fEnd - fStart) + 12*(cEnd - cStart); /* Every face is split into 4 faces, and 12 faces are added for each cell */ 107 depthSize[3] = 8*(cEnd - cStart); /* Every cell split into 8 cells */ 108 break; 109 case 7: 110 /* Hybrid Simplicial 3D */ 111 ierr = DMPlexGetNumHybridFaces_Internal(dm, &numHybridFaces);CHKERRQ(ierr); 112 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 113 cMax = PetscMin(cEnd, cMax); 114 if (eMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No edge maximum specified in hybrid mesh"); 115 eMax = PetscMin(eEnd, eMax); 116 /* Tetrahedra */ 117 depthSize[0] = vEnd - vStart + eMax - eStart; /* Add a vertex on every interior edge */ 118 depthSize[1] = 2*(eMax - eStart) + 3*(fMax - fStart) + (cMax - cStart); /* Every interior edge split into 2 edges, 3 edges added for each interior face, 1 edge for each interior cell */ 119 depthSize[2] = 4*(fMax - fStart) + 8*(cMax - cStart); /* Every interior face split into 4 faces, 8 faces added for each interior cell */ 120 depthSize[3] = 8*(cMax - cStart); /* Every interior cell split into 8 cells */ 121 /* Triangular Prisms */ 122 depthSize[0] += 0; /* No hybrid vertices */ 123 depthSize[1] += (eEnd - eMax) + (fEnd - fMax); /* Every hybrid edge remains, 1 edge for every hybrid face */ 124 depthSize[2] += 4*(fEnd - fMax) + 3*(cEnd - cMax); /* Every hybrid face split into 2 faces and 3 faces are added for each hybrid cell */ 125 depthSize[3] += 4*(cEnd - cMax); /* Every hybrid cell split into 4 cells */ 126 break; 127 default: 128 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 129 } 130 PetscFunctionReturn(0); 131 } 132 133 /* Return triangle edge for orientation o, if it is r for o == 0 */ 134 PETSC_STATIC_INLINE PetscInt GetTriEdge_Static(PetscInt o, PetscInt r) { 135 return (o < 0 ? 2-(o+r) : o+r)%3; 136 } 137 138 /* Return triangle subface for orientation o, if it is r for o == 0 */ 139 PETSC_STATIC_INLINE PetscInt GetTriSubface_Static(PetscInt o, PetscInt r) { 140 return (o < 0 ? 0-(o+r) : o+r)%3; 141 } 142 143 /* Return quad edge for orientation o, if it is r for o == 0 */ 144 PETSC_STATIC_INLINE PetscInt GetQuadEdge_Static(PetscInt o, PetscInt r) { 145 return (o < 0 ? 3-(o+r) : o+r)%4; 146 } 147 148 /* Return quad subface for orientation o, if it is r for o == 0 */ 149 PETSC_STATIC_INLINE PetscInt GetQuadSubface_Static(PetscInt o, PetscInt r) { 150 return (o < 0 ? 0-(o+r) : o+r)%4; 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "CellRefinerSetConeSizes" 155 PetscErrorCode CellRefinerSetConeSizes(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 156 { 157 PetscInt depth, cStart, cStartNew, cEnd, cMax, c, vStart, vStartNew, vEnd, vMax, v, fStart, fStartNew, fEnd, fMax, f, eStart, eStartNew, eEnd, eMax, e, r; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 162 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 163 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 164 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 165 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 166 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 167 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 168 switch (refiner) { 169 case 1: 170 /* Simplicial 2D */ 171 /* All cells have 3 faces */ 172 for (c = cStart; c < cEnd; ++c) { 173 for (r = 0; r < 4; ++r) { 174 const PetscInt newp = (c - cStart)*4 + r; 175 176 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 177 } 178 } 179 /* Split faces have 2 vertices and the same cells as the parent */ 180 for (f = fStart; f < fEnd; ++f) { 181 for (r = 0; r < 2; ++r) { 182 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 183 PetscInt size; 184 185 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 186 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 187 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 188 } 189 } 190 /* Interior faces have 2 vertices and 2 cells */ 191 for (c = cStart; c < cEnd; ++c) { 192 for (r = 0; r < 3; ++r) { 193 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + r; 194 195 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 196 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 197 } 198 } 199 /* Old vertices have identical supports */ 200 for (v = vStart; v < vEnd; ++v) { 201 const PetscInt newp = vStartNew + (v - vStart); 202 PetscInt size; 203 204 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 205 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 206 } 207 /* Face vertices have 2 + cells*2 supports */ 208 for (f = fStart; f < fEnd; ++f) { 209 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 210 PetscInt size; 211 212 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 213 ierr = DMPlexSetSupportSize(rdm, newp, 2 + size*2);CHKERRQ(ierr); 214 } 215 break; 216 case 2: 217 /* Hex 2D */ 218 /* All cells have 4 faces */ 219 for (c = cStart; c < cEnd; ++c) { 220 for (r = 0; r < 4; ++r) { 221 const PetscInt newp = (c - cStart)*4 + r; 222 223 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 224 } 225 } 226 /* Split faces have 2 vertices and the same cells as the parent */ 227 for (f = fStart; f < fEnd; ++f) { 228 for (r = 0; r < 2; ++r) { 229 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 230 PetscInt size; 231 232 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 233 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 234 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 235 } 236 } 237 /* Interior faces have 2 vertices and 2 cells */ 238 for (c = cStart; c < cEnd; ++c) { 239 for (r = 0; r < 4; ++r) { 240 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + r; 241 242 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 243 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 244 } 245 } 246 /* Old vertices have identical supports */ 247 for (v = vStart; v < vEnd; ++v) { 248 const PetscInt newp = vStartNew + (v - vStart); 249 PetscInt size; 250 251 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 252 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 253 } 254 /* Face vertices have 2 + cells supports */ 255 for (f = fStart; f < fEnd; ++f) { 256 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 257 PetscInt size; 258 259 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 260 ierr = DMPlexSetSupportSize(rdm, newp, 2 + size);CHKERRQ(ierr); 261 } 262 /* Cell vertices have 4 supports */ 263 for (c = cStart; c < cEnd; ++c) { 264 const PetscInt newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 265 266 ierr = DMPlexSetSupportSize(rdm, newp, 4);CHKERRQ(ierr); 267 } 268 break; 269 case 3: 270 /* Hybrid Simplicial 2D */ 271 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 272 cMax = PetscMin(cEnd, cMax); 273 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 274 fMax = PetscMin(fEnd, fMax); 275 ierr = DMPlexSetHybridBounds(rdm, cStartNew + (cMax - cStart)*4, fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 276 /* Interior cells have 3 faces */ 277 for (c = cStart; c < cMax; ++c) { 278 for (r = 0; r < 4; ++r) { 279 const PetscInt newp = cStartNew + (c - cStart)*4 + r; 280 281 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 282 } 283 } 284 /* Hybrid cells have 4 faces */ 285 for (c = cMax; c < cEnd; ++c) { 286 for (r = 0; r < 2; ++r) { 287 const PetscInt newp = cStartNew + (cMax - cStart)*4 + (c - cMax)*2 + r; 288 289 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 290 } 291 } 292 /* Interior split faces have 2 vertices and the same cells as the parent */ 293 for (f = fStart; f < fMax; ++f) { 294 for (r = 0; r < 2; ++r) { 295 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 296 PetscInt size; 297 298 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 299 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 300 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 301 } 302 } 303 /* Interior cell faces have 2 vertices and 2 cells */ 304 for (c = cStart; c < cMax; ++c) { 305 for (r = 0; r < 3; ++r) { 306 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + r; 307 308 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 309 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 310 } 311 } 312 /* Hybrid faces have 2 vertices and the same cells */ 313 for (f = fMax; f < fEnd; ++f) { 314 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (f - fMax); 315 PetscInt size; 316 317 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 318 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 319 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 320 } 321 /* Hybrid cell faces have 2 vertices and 2 cells */ 322 for (c = cMax; c < cEnd; ++c) { 323 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 324 325 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 326 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 327 } 328 /* Old vertices have identical supports */ 329 for (v = vStart; v < vEnd; ++v) { 330 const PetscInt newp = vStartNew + (v - vStart); 331 PetscInt size; 332 333 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 334 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 335 } 336 /* Face vertices have 2 + (2 interior, 1 hybrid) supports */ 337 for (f = fStart; f < fMax; ++f) { 338 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 339 const PetscInt *support; 340 PetscInt size, newSize = 2, s; 341 342 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 343 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 344 for (s = 0; s < size; ++s) { 345 if (support[s] >= cMax) newSize += 1; 346 else newSize += 2; 347 } 348 ierr = DMPlexSetSupportSize(rdm, newp, newSize);CHKERRQ(ierr); 349 } 350 break; 351 case 5: 352 /* Simplicial 3D */ 353 /* All cells have 4 faces */ 354 for (c = cStart; c < cEnd; ++c) { 355 for (r = 0; r < 8; ++r) { 356 const PetscInt newp = cStartNew + (c - cStart)*8 + r; 357 358 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 359 } 360 } 361 /* Split faces have 3 edges and the same cells as the parent */ 362 for (f = fStart; f < fEnd; ++f) { 363 for (r = 0; r < 4; ++r) { 364 const PetscInt newp = fStartNew + (f - fStart)*4 + r; 365 PetscInt size; 366 367 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 368 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 369 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 370 } 371 } 372 /* Interior faces have 3 edges and 2 cells */ 373 for (c = cStart; c < cEnd; ++c) { 374 for (r = 0; r < 8; ++r) { 375 const PetscInt newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + r; 376 377 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 378 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 379 } 380 } 381 /* Split edges have 2 vertices and the same faces */ 382 for (e = eStart; e < eEnd; ++e) { 383 for (r = 0; r < 2; ++r) { 384 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 385 PetscInt size; 386 387 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 388 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 389 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 390 } 391 } 392 /* Face edges have 2 vertices and 2+cells*(1/2) faces */ 393 for (f = fStart; f < fEnd; ++f) { 394 for (r = 0; r < 3; ++r) { 395 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + r; 396 const PetscInt *cone, *ornt, *support, eint[4] = {1, 0, 2, 0}; 397 PetscInt coneSize, c, supportSize, s, er, intFaces = 0; 398 399 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 400 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 401 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 402 for (s = 0; s < supportSize; ++s) { 403 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 404 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 405 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 406 for (c = 0; c < coneSize; ++c) {if (cone[c] == f) break;} 407 /* Here we want to determine whether edge newp contains a vertex which is part of the cross-tet edge */ 408 er = ornt[c] < 0 ? (-(ornt[c]+1) + 2-r)%3 : (ornt[c] + r)%3; 409 if (er == eint[c]) { 410 intFaces += 1; 411 } else { 412 intFaces += 2; 413 } 414 } 415 ierr = DMPlexSetSupportSize(rdm, newp, 2+intFaces);CHKERRQ(ierr); 416 } 417 } 418 /* Interior edges have 2 vertices and 4 faces */ 419 for (c = cStart; c < cEnd; ++c) { 420 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 421 422 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 423 ierr = DMPlexSetSupportSize(rdm, newp, 4);CHKERRQ(ierr); 424 } 425 /* Old vertices have identical supports */ 426 for (v = vStart; v < vEnd; ++v) { 427 const PetscInt newp = vStartNew + (v - vStart); 428 PetscInt size; 429 430 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 431 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 432 } 433 /* Edge vertices have 2 + faces*2 + cells*0/1 supports */ 434 for (e = eStart; e < eEnd; ++e) { 435 const PetscInt newp = vStartNew + (vEnd - vStart) + (e - eStart); 436 PetscInt size, *star = NULL, starSize, s, cellSize = 0; 437 438 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 439 ierr = DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 440 for (s = 0; s < starSize*2; s += 2) { 441 const PetscInt *cone, *ornt; 442 PetscInt e01, e23; 443 444 if ((star[s] >= cStart) && (star[s] < cEnd)) { 445 /* Check edge 0-1 */ 446 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 447 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 448 ierr = DMPlexGetCone(dm, cone[0], &cone);CHKERRQ(ierr); 449 e01 = cone[GetTriEdge_Static(ornt[0], 0)]; 450 /* Check edge 2-3 */ 451 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 452 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 453 ierr = DMPlexGetCone(dm, cone[2], &cone);CHKERRQ(ierr); 454 e23 = cone[GetTriEdge_Static(ornt[2], 1)]; 455 if ((e01 == e) || (e23 == e)) ++cellSize; 456 } 457 } 458 ierr = DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 459 ierr = DMPlexSetSupportSize(rdm, newp, 2 + size*2 + cellSize);CHKERRQ(ierr); 460 } 461 break; 462 case 7: 463 /* Hybrid Simplicial 2D */ 464 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 465 cMax = PetscMin(cEnd, cMax); 466 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 467 fMax = PetscMin(fEnd, fMax); 468 ierr = DMPlexSetHybridBounds(rdm, cStartNew + (cMax-cStart)*8, PETSC_DETERMINE, eStartNew + (eMax-eStart)*2 + (cMax-cStart)*13, PETSC_DETERMINE);CHKERRQ(ierr); 469 /* Interior cells have 4 faces */ 470 for (c = cStart; c < cMax; ++c) { 471 for (r = 0; r < 8; ++r) { 472 const PetscInt newp = cStartNew + (c - cStart)*8 + r; 473 474 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 475 } 476 } 477 /* Hybrid cells have 5 faces */ 478 for (c = cMax; c < cEnd; ++c) { 479 for (r = 0; r < 4; ++r) { 480 const PetscInt newp = cStartNew + (cMax - cStart)*8 + (c - cMax)*4 + r; 481 482 ierr = DMPlexSetConeSize(rdm, newp, 5);CHKERRQ(ierr); 483 } 484 } 485 /* Interior/Hybrid split faces have 3 edges and the same cells as the parent */ 486 for (f = fStart; f < fMax; ++f) { 487 for (r = 0; r < 4; ++r) { 488 const PetscInt newp = fStartNew + (f - fStart)*4 + r; 489 PetscInt size; 490 491 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 492 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 493 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 494 } 495 } 496 /* Interior cell faces have 3 edges and 2 cells */ 497 for (c = cStart; c < cMax; ++c) { 498 for (r = 0; r < 8; ++r) { 499 const PetscInt newp = fStartNew + (fMax - fStart)*4 + (c - cStart)*8 + r; 500 501 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 502 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 503 } 504 } 505 /* Split edges have 2 vertices and the same faces */ 506 for (e = eStart; e < eMax; ++e) { 507 for (r = 0; r < 2; ++r) { 508 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 509 PetscInt size; 510 511 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 512 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 513 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 514 } 515 } 516 /* Face edges have 2 vertices and 2+cells*(1/2) faces */ 517 for (f = fStart; f < fMax; ++f) { 518 for (r = 0; r < 3; ++r) { 519 const PetscInt newp = eStartNew + (eMax - eStart)*2 + (f - fStart)*3 + r; 520 const PetscInt *cone, *ornt, *support, eint[4] = {1, 0, 2, 0}; 521 PetscInt coneSize, c, supportSize, s, er, intFaces = 0; 522 523 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 524 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 525 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 526 for (s = 0; s < supportSize; ++s) { 527 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 528 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 529 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 530 for (c = 0; c < coneSize; ++c) {if (cone[c] == f) break;} 531 /* Here we want to determine whether edge newp contains a vertex which is part of the cross-tet edge */ 532 er = ornt[c] < 0 ? (-(ornt[c]+1) + 2-r)%3 : (ornt[c] + r)%3; 533 if (er == eint[c]) { 534 intFaces += 1; 535 } else { 536 intFaces += 2; 537 } 538 } 539 ierr = DMPlexSetSupportSize(rdm, newp, 2+intFaces);CHKERRQ(ierr); 540 } 541 } 542 /* Interior edges have 2 vertices and 4 faces */ 543 for (c = cStart; c < cMax; ++c) { 544 const PetscInt newp = eStartNew + (eMax - eStart)*2 + (fMax - fStart)*3 + (c - cStart); 545 546 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 547 ierr = DMPlexSetSupportSize(rdm, newp, 4);CHKERRQ(ierr); 548 } 549 /* Hybrid edges have 2 vertices and the same cells */ 550 for (e = eMax; e < eEnd; ++e) { 551 const PetscInt newp = eStartNew + (eMax - eStart)*2 + (fMax - fStart)*3 + (cMax - cStart) + (e - eEnd); 552 PetscInt size; 553 554 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 555 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 556 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 557 } 558 /* Hybrid split edges have 2 vertices and 2 faces */ 559 for (f = fMax; f < fEnd; ++f) { 560 const PetscInt newp = eStartNew + (eMax - eStart)*2 + (fMax - fStart)*3 + (cMax - cStart) + (eEnd - eMax) + (f - fMax); 561 562 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 563 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 564 } 565 /* Old vertices have identical supports */ 566 for (v = vStart; v < vEnd; ++v) { 567 const PetscInt newp = vStartNew + (v - vStart); 568 PetscInt size; 569 570 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 571 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 572 } 573 /* Face vertices have 2 + (2 interior, 1 hybrid) supports */ 574 for (f = fStart; f < fMax; ++f) { 575 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 576 const PetscInt *support; 577 PetscInt size, newSize = 2, s; 578 579 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 580 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 581 for (s = 0; s < size; ++s) { 582 if (support[s] >= cMax) newSize += 1; 583 else newSize += 2; 584 } 585 ierr = DMPlexSetSupportSize(rdm, newp, newSize);CHKERRQ(ierr); 586 } 587 break; 588 case 6: 589 /* Hex 3D */ 590 /* All cells have 6 faces */ 591 for (c = cStart; c < cEnd; ++c) { 592 for (r = 0; r < 8; ++r) { 593 const PetscInt newp = (c - cStart)*8 + r; 594 595 ierr = DMPlexSetConeSize(rdm, newp, 6);CHKERRQ(ierr); 596 } 597 } 598 /* Split faces have 4 edges and the same cells as the parent */ 599 for (f = fStart; f < fEnd; ++f) { 600 for (r = 0; r < 4; ++r) { 601 const PetscInt newp = fStartNew + (f - fStart)*4 + r; 602 PetscInt size; 603 604 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 605 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 606 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 607 } 608 } 609 /* Interior faces have 4 edges and 2 cells */ 610 for (c = cStart; c < cEnd; ++c) { 611 for (r = 0; r < 12; ++r) { 612 const PetscInt newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + r; 613 614 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 615 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 616 } 617 } 618 /* Split edges have 2 vertices and the same faces as the parent */ 619 for (e = eStart; e < eEnd; ++e) { 620 for (r = 0; r < 2; ++r) { 621 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 622 PetscInt size; 623 624 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 625 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 626 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 627 } 628 } 629 /* Face edges have 2 vertices and 2+cells faces */ 630 for (f = fStart; f < fEnd; ++f) { 631 for (r = 0; r < 4; ++r) { 632 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + r; 633 PetscInt size; 634 635 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 636 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 637 ierr = DMPlexSetSupportSize(rdm, newp, 2+size);CHKERRQ(ierr); 638 } 639 } 640 /* Cell edges have 2 vertices and 4 faces */ 641 for (c = cStart; c < cEnd; ++c) { 642 for (r = 0; r < 6; ++r) { 643 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + r; 644 645 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 646 ierr = DMPlexSetSupportSize(rdm, newp, 4);CHKERRQ(ierr); 647 } 648 } 649 /* Old vertices have identical supports */ 650 for (v = vStart; v < vEnd; ++v) { 651 const PetscInt newp = vStartNew + (v - vStart); 652 PetscInt size; 653 654 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 655 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 656 } 657 /* Edge vertices have 2 + faces supports */ 658 for (e = eStart; e < eEnd; ++e) { 659 const PetscInt newp = vStartNew + (vEnd - vStart) + (e - eStart); 660 PetscInt size; 661 662 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 663 ierr = DMPlexSetSupportSize(rdm, newp, 2 + size);CHKERRQ(ierr); 664 } 665 /* Face vertices have 4 + cells supports */ 666 for (f = fStart; f < fEnd; ++f) { 667 const PetscInt newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 668 PetscInt size; 669 670 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 671 ierr = DMPlexSetSupportSize(rdm, newp, 4 + size);CHKERRQ(ierr); 672 } 673 /* Cell vertices have 6 supports */ 674 for (c = cStart; c < cEnd; ++c) { 675 const PetscInt newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (c - cStart); 676 677 ierr = DMPlexSetSupportSize(rdm, newp, 6);CHKERRQ(ierr); 678 } 679 break; 680 default: 681 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 682 } 683 PetscFunctionReturn(0); 684 } 685 686 #undef __FUNCT__ 687 #define __FUNCT__ "CellRefinerSetCones" 688 PetscErrorCode CellRefinerSetCones(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 689 { 690 const PetscInt *faces, cellInd[4] = {0, 1, 2, 3}; 691 PetscInt depth, cStart, cEnd, cMax, cStartNew, cEndNew, c, vStart, vEnd, vMax, vStartNew, vEndNew, v, fStart, fEnd, fMax, fStartNew, fEndNew, f, eStart, eEnd, eMax, eStartNew, eEndNew, e, r, p; 692 PetscInt maxSupportSize, *supportRef; 693 PetscErrorCode ierr; 694 695 PetscFunctionBegin; 696 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 697 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 698 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 699 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 700 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 701 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 702 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 703 ierr = GetDepthEnd_Private(depth, depthSize, &cEndNew, &fEndNew, &eEndNew, &vEndNew);CHKERRQ(ierr); 704 switch (refiner) { 705 case 1: 706 /* Simplicial 2D */ 707 /* 708 2 709 |\ 710 | \ 711 | \ 712 | \ 713 | C \ 714 | \ 715 | \ 716 2---1---1 717 |\ D / \ 718 | 2 0 \ 719 |A \ / B \ 720 0---0-------1 721 */ 722 /* All cells have 3 faces */ 723 for (c = cStart; c < cEnd; ++c) { 724 const PetscInt newp = cStartNew + (c - cStart)*4; 725 const PetscInt *cone, *ornt; 726 PetscInt coneNew[3], orntNew[3]; 727 728 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 729 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 730 /* A triangle */ 731 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 732 orntNew[0] = ornt[0]; 733 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 2; 734 orntNew[1] = -2; 735 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 736 orntNew[2] = ornt[2]; 737 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 738 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 739 #if 1 740 if ((newp+0 < cStartNew) || (newp+0 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+0, cStartNew, cEndNew); 741 for (p = 0; p < 3; ++p) { 742 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 743 } 744 #endif 745 /* B triangle */ 746 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 747 orntNew[0] = ornt[0]; 748 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 749 orntNew[1] = ornt[1]; 750 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 0; 751 orntNew[2] = -2; 752 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 753 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 754 #if 1 755 if ((newp+1 < cStartNew) || (newp+1 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+1, cStartNew, cEndNew); 756 for (p = 0; p < 3; ++p) { 757 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 758 } 759 #endif 760 /* C triangle */ 761 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 1; 762 orntNew[0] = -2; 763 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 764 orntNew[1] = ornt[1]; 765 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 766 orntNew[2] = ornt[2]; 767 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 768 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 769 #if 1 770 if ((newp+2 < cStartNew) || (newp+2 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+2, cStartNew, cEndNew); 771 for (p = 0; p < 3; ++p) { 772 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 773 } 774 #endif 775 /* D triangle */ 776 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 0; 777 orntNew[0] = 0; 778 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 1; 779 orntNew[1] = 0; 780 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 2; 781 orntNew[2] = 0; 782 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 783 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 784 #if 1 785 if ((newp+3 < cStartNew) || (newp+3 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+3, cStartNew, cEndNew); 786 for (p = 0; p < 3; ++p) { 787 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 788 } 789 #endif 790 } 791 /* Split faces have 2 vertices and the same cells as the parent */ 792 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 793 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 794 for (f = fStart; f < fEnd; ++f) { 795 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 796 797 for (r = 0; r < 2; ++r) { 798 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 799 const PetscInt *cone, *ornt, *support; 800 PetscInt coneNew[2], coneSize, c, supportSize, s; 801 802 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 803 coneNew[0] = vStartNew + (cone[0] - vStart); 804 coneNew[1] = vStartNew + (cone[1] - vStart); 805 coneNew[(r+1)%2] = newv; 806 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 807 #if 1 808 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 809 for (p = 0; p < 2; ++p) { 810 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 811 } 812 #endif 813 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 814 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 815 for (s = 0; s < supportSize; ++s) { 816 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 817 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 818 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 819 for (c = 0; c < coneSize; ++c) { 820 if (cone[c] == f) break; 821 } 822 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (ornt[c] < 0 ? (c+1-r)%3 : (c+r)%3); 823 } 824 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 825 #if 1 826 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 827 for (p = 0; p < supportSize; ++p) { 828 if ((supportRef[p] < cStartNew) || (supportRef[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportRef[p], cStartNew, cEndNew); 829 } 830 #endif 831 } 832 } 833 /* Interior faces have 2 vertices and 2 cells */ 834 for (c = cStart; c < cEnd; ++c) { 835 const PetscInt *cone; 836 837 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 838 for (r = 0; r < 3; ++r) { 839 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + r; 840 PetscInt coneNew[2]; 841 PetscInt supportNew[2]; 842 843 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 844 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - fStart); 845 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 846 #if 1 847 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 848 for (p = 0; p < 2; ++p) { 849 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 850 } 851 #endif 852 supportNew[0] = (c - cStart)*4 + (r+1)%3; 853 supportNew[1] = (c - cStart)*4 + 3; 854 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 855 #if 1 856 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 857 for (p = 0; p < 2; ++p) { 858 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 859 } 860 #endif 861 } 862 } 863 /* Old vertices have identical supports */ 864 for (v = vStart; v < vEnd; ++v) { 865 const PetscInt newp = vStartNew + (v - vStart); 866 const PetscInt *support, *cone; 867 PetscInt size, s; 868 869 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 870 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 871 for (s = 0; s < size; ++s) { 872 PetscInt r = 0; 873 874 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 875 if (cone[1] == v) r = 1; 876 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 877 } 878 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 879 #if 1 880 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 881 for (p = 0; p < size; ++p) { 882 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 883 } 884 #endif 885 } 886 /* Face vertices have 2 + cells*2 supports */ 887 for (f = fStart; f < fEnd; ++f) { 888 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 889 const PetscInt *cone, *support; 890 PetscInt size, s; 891 892 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 893 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 894 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 895 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 896 for (s = 0; s < size; ++s) { 897 PetscInt r = 0; 898 899 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 900 if (cone[1] == f) r = 1; 901 else if (cone[2] == f) r = 2; 902 supportRef[2+s*2+0] = fStartNew + (fEnd - fStart)*2 + (support[s] - cStart)*3 + (r+2)%3; 903 supportRef[2+s*2+1] = fStartNew + (fEnd - fStart)*2 + (support[s] - cStart)*3 + r; 904 } 905 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 906 #if 1 907 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 908 for (p = 0; p < 2+size*2; ++p) { 909 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 910 } 911 #endif 912 } 913 ierr = PetscFree(supportRef);CHKERRQ(ierr); 914 break; 915 case 2: 916 /* Hex 2D */ 917 /* 918 3---------2---------2 919 | | | 920 | D 2 C | 921 | | | 922 3----3----0----1----1 923 | | | 924 | A 0 B | 925 | | | 926 0---------0---------1 927 */ 928 /* All cells have 4 faces */ 929 for (c = cStart; c < cEnd; ++c) { 930 const PetscInt newp = (c - cStart)*4; 931 const PetscInt *cone, *ornt; 932 PetscInt coneNew[4], orntNew[4]; 933 934 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 935 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 936 /* A quad */ 937 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 938 orntNew[0] = ornt[0]; 939 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 0; 940 orntNew[1] = 0; 941 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 3; 942 orntNew[2] = -2; 943 coneNew[3] = fStartNew + (cone[3] - fStart)*2 + (ornt[3] < 0 ? 0 : 1); 944 orntNew[3] = ornt[3]; 945 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 946 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 947 #if 1 948 if ((newp+0 < cStartNew) || (newp+0 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+0, cStartNew, cEndNew); 949 for (p = 0; p < 4; ++p) { 950 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 951 } 952 #endif 953 /* B quad */ 954 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 955 orntNew[0] = ornt[0]; 956 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 957 orntNew[1] = ornt[1]; 958 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 1; 959 orntNew[2] = 0; 960 coneNew[3] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 0; 961 orntNew[3] = -2; 962 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 963 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 964 #if 1 965 if ((newp+1 < cStartNew) || (newp+1 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+1, cStartNew, cEndNew); 966 for (p = 0; p < 4; ++p) { 967 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 968 } 969 #endif 970 /* C quad */ 971 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 1; 972 orntNew[0] = -2; 973 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 974 orntNew[1] = ornt[1]; 975 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 976 orntNew[2] = ornt[2]; 977 coneNew[3] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 2; 978 orntNew[3] = 0; 979 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 980 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 981 #if 1 982 if ((newp+2 < cStartNew) || (newp+2 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+2, cStartNew, cEndNew); 983 for (p = 0; p < 4; ++p) { 984 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 985 } 986 #endif 987 /* D quad */ 988 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 3; 989 orntNew[0] = 0; 990 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 2; 991 orntNew[1] = -2; 992 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 993 orntNew[2] = ornt[2]; 994 coneNew[3] = fStartNew + (cone[3] - fStart)*2 + (ornt[3] < 0 ? 1 : 0); 995 orntNew[3] = ornt[3]; 996 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 997 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 998 #if 1 999 if ((newp+3 < cStartNew) || (newp+3 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+3, cStartNew, cEndNew); 1000 for (p = 0; p < 4; ++p) { 1001 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1002 } 1003 #endif 1004 } 1005 /* Split faces have 2 vertices and the same cells as the parent */ 1006 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 1007 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 1008 for (f = fStart; f < fEnd; ++f) { 1009 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 1010 1011 for (r = 0; r < 2; ++r) { 1012 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 1013 const PetscInt *cone, *ornt, *support; 1014 PetscInt coneNew[2], coneSize, c, supportSize, s; 1015 1016 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1017 coneNew[0] = vStartNew + (cone[0] - vStart); 1018 coneNew[1] = vStartNew + (cone[1] - vStart); 1019 coneNew[(r+1)%2] = newv; 1020 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1021 #if 1 1022 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1023 for (p = 0; p < 2; ++p) { 1024 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 1025 } 1026 #endif 1027 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1028 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1029 for (s = 0; s < supportSize; ++s) { 1030 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1031 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1032 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1033 for (c = 0; c < coneSize; ++c) { 1034 if (cone[c] == f) break; 1035 } 1036 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (ornt[c] < 0 ? (c+1-r)%4 : (c+r)%4); 1037 } 1038 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1039 #if 1 1040 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1041 for (p = 0; p < supportSize; ++p) { 1042 if ((supportRef[p] < cStartNew) || (supportRef[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportRef[p], cStartNew, cEndNew); 1043 } 1044 #endif 1045 } 1046 } 1047 /* Interior faces have 2 vertices and 2 cells */ 1048 for (c = cStart; c < cEnd; ++c) { 1049 const PetscInt *cone; 1050 PetscInt coneNew[2], supportNew[2]; 1051 1052 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1053 for (r = 0; r < 4; ++r) { 1054 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + r; 1055 1056 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 1057 coneNew[1] = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 1058 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1059 #if 1 1060 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1061 for (p = 0; p < 2; ++p) { 1062 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 1063 } 1064 #endif 1065 supportNew[0] = (c - cStart)*4 + r; 1066 supportNew[1] = (c - cStart)*4 + (r+1)%4; 1067 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1068 #if 1 1069 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1070 for (p = 0; p < 2; ++p) { 1071 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1072 } 1073 #endif 1074 } 1075 } 1076 /* Old vertices have identical supports */ 1077 for (v = vStart; v < vEnd; ++v) { 1078 const PetscInt newp = vStartNew + (v - vStart); 1079 const PetscInt *support, *cone; 1080 PetscInt size, s; 1081 1082 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 1083 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 1084 for (s = 0; s < size; ++s) { 1085 PetscInt r = 0; 1086 1087 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1088 if (cone[1] == v) r = 1; 1089 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 1090 } 1091 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1092 #if 1 1093 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1094 for (p = 0; p < size; ++p) { 1095 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 1096 } 1097 #endif 1098 } 1099 /* Face vertices have 2 + cells supports */ 1100 for (f = fStart; f < fEnd; ++f) { 1101 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 1102 const PetscInt *cone, *support; 1103 PetscInt size, s; 1104 1105 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 1106 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1107 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 1108 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 1109 for (s = 0; s < size; ++s) { 1110 PetscInt r = 0; 1111 1112 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1113 if (cone[1] == f) r = 1; 1114 else if (cone[2] == f) r = 2; 1115 else if (cone[3] == f) r = 3; 1116 supportRef[2+s] = fStartNew + (fEnd - fStart)*2 + (support[s] - cStart)*4 + r; 1117 } 1118 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1119 #if 1 1120 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1121 for (p = 0; p < 2+size; ++p) { 1122 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 1123 } 1124 #endif 1125 } 1126 /* Cell vertices have 4 supports */ 1127 for (c = cStart; c < cEnd; ++c) { 1128 const PetscInt newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 1129 PetscInt supportNew[4]; 1130 1131 for (r = 0; r < 4; ++r) { 1132 supportNew[r] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + r; 1133 } 1134 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1135 } 1136 ierr = PetscFree(supportRef);CHKERRQ(ierr); 1137 break; 1138 case 3: 1139 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 1140 cMax = PetscMin(cEnd, cMax); 1141 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 1142 fMax = PetscMin(fEnd, fMax); 1143 /* Interior cells have 3 faces */ 1144 for (c = cStart; c < cMax; ++c) { 1145 const PetscInt newp = cStartNew + (c - cStart)*4; 1146 const PetscInt *cone, *ornt; 1147 PetscInt coneNew[3], orntNew[3]; 1148 1149 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1150 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1151 /* A triangle */ 1152 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 1153 orntNew[0] = ornt[0]; 1154 coneNew[1] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 2; 1155 orntNew[1] = -2; 1156 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 1157 orntNew[2] = ornt[2]; 1158 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1159 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1160 #if 1 1161 if ((newp+0 < cStartNew) || (newp+0 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+0, cStartNew, cEndNew); 1162 for (p = 0; p < 3; ++p) { 1163 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1164 } 1165 #endif 1166 /* B triangle */ 1167 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 1168 orntNew[0] = ornt[0]; 1169 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 1170 orntNew[1] = ornt[1]; 1171 coneNew[2] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 0; 1172 orntNew[2] = -2; 1173 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1174 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1175 #if 1 1176 if ((newp+1 < cStartNew) || (newp+1 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+1, cStartNew, cEndNew); 1177 for (p = 0; p < 3; ++p) { 1178 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1179 } 1180 #endif 1181 /* C triangle */ 1182 coneNew[0] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 1; 1183 orntNew[0] = -2; 1184 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 1185 orntNew[1] = ornt[1]; 1186 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 1187 orntNew[2] = ornt[2]; 1188 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 1189 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 1190 #if 1 1191 if ((newp+2 < cStartNew) || (newp+2 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+2, cStartNew, cEndNew); 1192 for (p = 0; p < 3; ++p) { 1193 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1194 } 1195 #endif 1196 /* D triangle */ 1197 coneNew[0] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 0; 1198 orntNew[0] = 0; 1199 coneNew[1] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 1; 1200 orntNew[1] = 0; 1201 coneNew[2] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 2; 1202 orntNew[2] = 0; 1203 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 1204 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 1205 #if 1 1206 if ((newp+3 < cStartNew) || (newp+3 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+3, cStartNew, cEndNew); 1207 for (p = 0; p < 3; ++p) { 1208 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1209 } 1210 #endif 1211 } 1212 /* 1213 2----3----3 1214 | | 1215 | B | 1216 | | 1217 0----4--- 1 1218 | | 1219 | A | 1220 | | 1221 0----2----1 1222 */ 1223 /* Hybrid cells have 4 faces */ 1224 for (c = cMax; c < cEnd; ++c) { 1225 const PetscInt newp = cStartNew + (cMax - cStart)*4 + (c - cMax)*2; 1226 const PetscInt *cone, *ornt; 1227 PetscInt coneNew[4], orntNew[4]; 1228 1229 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1230 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1231 /* A quad */ 1232 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 1233 orntNew[0] = ornt[0]; 1234 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 1235 orntNew[1] = ornt[1]; 1236 coneNew[2] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (cone[2] - fMax); 1237 orntNew[2] = 0; 1238 coneNew[3] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 1239 orntNew[3] = 0; 1240 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1241 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1242 #if 1 1243 if ((newp+0 < cStartNew) || (newp+0 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+0, cStartNew, cEndNew); 1244 for (p = 0; p < 4; ++p) { 1245 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1246 } 1247 #endif 1248 /* B quad */ 1249 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 1250 orntNew[0] = ornt[0]; 1251 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 1252 orntNew[1] = ornt[1]; 1253 coneNew[2] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 1254 orntNew[2] = 0; 1255 coneNew[3] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (cone[3] - fMax); 1256 orntNew[3] = 0; 1257 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1258 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1259 #if 1 1260 if ((newp+1 < cStartNew) || (newp+1 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+1, cStartNew, cEndNew); 1261 for (p = 0; p < 4; ++p) { 1262 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1263 } 1264 #endif 1265 } 1266 /* Interior split faces have 2 vertices and the same cells as the parent */ 1267 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 1268 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 1269 for (f = fStart; f < fMax; ++f) { 1270 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 1271 1272 for (r = 0; r < 2; ++r) { 1273 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 1274 const PetscInt *cone, *ornt, *support; 1275 PetscInt coneNew[2], coneSize, c, supportSize, s; 1276 1277 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1278 coneNew[0] = vStartNew + (cone[0] - vStart); 1279 coneNew[1] = vStartNew + (cone[1] - vStart); 1280 coneNew[(r+1)%2] = newv; 1281 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1282 #if 1 1283 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1284 for (p = 0; p < 2; ++p) { 1285 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 1286 } 1287 #endif 1288 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1289 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1290 for (s = 0; s < supportSize; ++s) { 1291 if (support[s] >= cMax) { 1292 supportRef[s] = cStartNew + (cMax - cStart)*4 + (support[s] - cMax)*2 + r; 1293 } else { 1294 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1295 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1296 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1297 for (c = 0; c < coneSize; ++c) { 1298 if (cone[c] == f) break; 1299 } 1300 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (ornt[c] < 0 ? (c+1-r)%3 : (c+r)%3); 1301 } 1302 } 1303 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1304 #if 1 1305 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1306 for (p = 0; p < supportSize; ++p) { 1307 if ((supportRef[p] < cStartNew) || (supportRef[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportRef[p], cStartNew, cEndNew); 1308 } 1309 #endif 1310 } 1311 } 1312 /* Interior cell faces have 2 vertices and 2 cells */ 1313 for (c = cStart; c < cMax; ++c) { 1314 const PetscInt *cone; 1315 1316 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1317 for (r = 0; r < 3; ++r) { 1318 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + r; 1319 PetscInt coneNew[2]; 1320 PetscInt supportNew[2]; 1321 1322 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 1323 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - fStart); 1324 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1325 #if 1 1326 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1327 for (p = 0; p < 2; ++p) { 1328 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 1329 } 1330 #endif 1331 supportNew[0] = (c - cStart)*4 + (r+1)%3; 1332 supportNew[1] = (c - cStart)*4 + 3; 1333 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1334 #if 1 1335 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1336 for (p = 0; p < 2; ++p) { 1337 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1338 } 1339 #endif 1340 } 1341 } 1342 /* Interior hybrid faces have 2 vertices and the same cells */ 1343 for (f = fMax; f < fEnd; ++f) { 1344 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (f - fMax); 1345 const PetscInt *cone; 1346 const PetscInt *support; 1347 PetscInt coneNew[2]; 1348 PetscInt supportNew[2]; 1349 PetscInt size, s, r; 1350 1351 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1352 coneNew[0] = vStartNew + (cone[0] - vStart); 1353 coneNew[1] = vStartNew + (cone[1] - vStart); 1354 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1355 #if 1 1356 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1357 for (p = 0; p < 2; ++p) { 1358 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 1359 } 1360 #endif 1361 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 1362 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1363 for (s = 0; s < size; ++s) { 1364 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1365 for (r = 0; r < 2; ++r) { 1366 if (cone[r+2] == f) break; 1367 } 1368 supportNew[s] = (cMax - cStart)*4 + (support[s] - cMax)*2 + r; 1369 } 1370 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1371 #if 1 1372 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1373 for (p = 0; p < size; ++p) { 1374 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1375 } 1376 #endif 1377 } 1378 /* Cell hybrid faces have 2 vertices and 2 cells */ 1379 for (c = cMax; c < cEnd; ++c) { 1380 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 1381 const PetscInt *cone; 1382 PetscInt coneNew[2]; 1383 PetscInt supportNew[2]; 1384 1385 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1386 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[0] - fStart); 1387 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[1] - fStart); 1388 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1389 #if 1 1390 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1391 for (p = 0; p < 2; ++p) { 1392 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 1393 } 1394 #endif 1395 supportNew[0] = (cMax - cStart)*4 + (c - cMax)*2 + 0; 1396 supportNew[1] = (cMax - cStart)*4 + (c - cMax)*2 + 1; 1397 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1398 #if 1 1399 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1400 for (p = 0; p < 2; ++p) { 1401 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1402 } 1403 #endif 1404 } 1405 /* Old vertices have identical supports */ 1406 for (v = vStart; v < vEnd; ++v) { 1407 const PetscInt newp = vStartNew + (v - vStart); 1408 const PetscInt *support, *cone; 1409 PetscInt size, s; 1410 1411 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 1412 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 1413 for (s = 0; s < size; ++s) { 1414 if (support[s] >= fMax) { 1415 supportRef[s] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (support[s] - fMax); 1416 } else { 1417 PetscInt r = 0; 1418 1419 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1420 if (cone[1] == v) r = 1; 1421 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 1422 } 1423 } 1424 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1425 #if 1 1426 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1427 for (p = 0; p < size; ++p) { 1428 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 1429 } 1430 #endif 1431 } 1432 /* Face vertices have 2 + (2 interior, 1 hybrid) supports */ 1433 for (f = fStart; f < fMax; ++f) { 1434 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 1435 const PetscInt *cone, *support; 1436 PetscInt size, newSize = 2, s; 1437 1438 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 1439 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1440 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 1441 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 1442 for (s = 0; s < size; ++s) { 1443 PetscInt r = 0; 1444 1445 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1446 if (support[s] >= cMax) { 1447 supportRef[newSize+0] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (support[s] - cMax); 1448 1449 newSize += 1; 1450 } else { 1451 if (cone[1] == f) r = 1; 1452 else if (cone[2] == f) r = 2; 1453 supportRef[newSize+0] = fStartNew + (fMax - fStart)*2 + (support[s] - cStart)*3 + (r+2)%3; 1454 supportRef[newSize+1] = fStartNew + (fMax - fStart)*2 + (support[s] - cStart)*3 + r; 1455 1456 newSize += 2; 1457 } 1458 } 1459 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1460 #if 1 1461 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1462 for (p = 0; p < newSize; ++p) { 1463 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 1464 } 1465 #endif 1466 } 1467 ierr = PetscFree(supportRef);CHKERRQ(ierr); 1468 break; 1469 case 5: 1470 /* Simplicial 3D */ 1471 /* All cells have 4 faces: Tet face order is prescribed in DMPlexGetFaces_Internal() */ 1472 ierr = DMPlexGetRawFaces_Internal(dm, 3, 4, cellInd, NULL, NULL, &faces);CHKERRQ(ierr); 1473 for (c = cStart; c < cEnd; ++c) { 1474 const PetscInt newp = cStartNew + (c - cStart)*8; 1475 const PetscInt *cone, *ornt; 1476 PetscInt coneNew[4], orntNew[4]; 1477 1478 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1479 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1480 /* A tetrahedron: {0, a, c, d} */ 1481 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetTriSubface_Static(ornt[0], 0); /* A */ 1482 orntNew[0] = ornt[0]; 1483 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetTriSubface_Static(ornt[1], 0); /* A */ 1484 orntNew[1] = ornt[1]; 1485 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetTriSubface_Static(ornt[2], 0); /* A */ 1486 orntNew[2] = ornt[2]; 1487 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 0; 1488 orntNew[3] = 0; 1489 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1490 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1491 #if 1 1492 if ((newp+0 < cStartNew) || (newp+0 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+0, cStartNew, cEndNew); 1493 for (p = 0; p < 4; ++p) { 1494 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1495 } 1496 #endif 1497 /* B tetrahedron: {a, 1, b, e} */ 1498 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetTriSubface_Static(ornt[0], 1); /* B */ 1499 orntNew[0] = ornt[0]; 1500 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetTriSubface_Static(ornt[1], 2); /* C */ 1501 orntNew[1] = ornt[1]; 1502 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 1; 1503 orntNew[2] = 0; 1504 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetTriSubface_Static(ornt[3], 1); /* B */ 1505 orntNew[3] = ornt[3]; 1506 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1507 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1508 #if 1 1509 if ((newp+1 < cStartNew) || (newp+1 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+1, cStartNew, cEndNew); 1510 for (p = 0; p < 4; ++p) { 1511 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1512 } 1513 #endif 1514 /* C tetrahedron: {c, b, 2, f} */ 1515 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetTriSubface_Static(ornt[0], 2); /* C */ 1516 orntNew[0] = ornt[0]; 1517 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 2; 1518 orntNew[1] = 0; 1519 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetTriSubface_Static(ornt[2], 1); /* B */ 1520 orntNew[2] = ornt[2]; 1521 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetTriSubface_Static(ornt[3], 0); /* A */ 1522 orntNew[3] = ornt[3]; 1523 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 1524 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 1525 #if 1 1526 if ((newp+2 < cStartNew) || (newp+2 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+2, cStartNew, cEndNew); 1527 for (p = 0; p < 4; ++p) { 1528 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1529 } 1530 #endif 1531 /* D tetrahedron: {d, e, f, 3} */ 1532 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 3; 1533 orntNew[0] = 0; 1534 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetTriSubface_Static(ornt[1], 1); /* B */ 1535 orntNew[1] = ornt[1]; 1536 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetTriSubface_Static(ornt[2], 2); /* C */ 1537 orntNew[2] = ornt[2]; 1538 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetTriSubface_Static(ornt[3], 2); /* C */ 1539 orntNew[3] = ornt[3]; 1540 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 1541 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 1542 #if 1 1543 if ((newp+3 < cStartNew) || (newp+3 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+3, cStartNew, cEndNew); 1544 for (p = 0; p < 4; ++p) { 1545 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1546 } 1547 #endif 1548 /* A' tetrahedron: {d, a, c, f} */ 1549 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 0; 1550 orntNew[0] = -3; 1551 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 4; 1552 orntNew[1] = 0; 1553 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + 3; 1554 orntNew[2] = ornt[2] < 0 ? -((-(ornt[2]+1)+2)%3+1) : (ornt[2]+2)%3; 1555 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 5; 1556 orntNew[3] = 0; 1557 ierr = DMPlexSetCone(rdm, newp+4, coneNew);CHKERRQ(ierr); 1558 ierr = DMPlexSetConeOrientation(rdm, newp+4, orntNew);CHKERRQ(ierr); 1559 #if 1 1560 if ((newp+4 < cStartNew) || (newp+4 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+4, cStartNew, cEndNew); 1561 for (p = 0; p < 4; ++p) { 1562 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1563 } 1564 #endif 1565 /* B' tetrahedron: {e, b, a, f} */ 1566 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 1; 1567 orntNew[0] = -3; 1568 coneNew[1] = fStartNew + (cone[3] - fStart)*4 + 3; 1569 orntNew[1] = ornt[3] < 0 ? -((-(ornt[3]+1)+1)%3+1) : (ornt[3]+1)%3; 1570 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 6; 1571 orntNew[2] = 0; 1572 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 7; 1573 orntNew[3] = 0; 1574 ierr = DMPlexSetCone(rdm, newp+5, coneNew);CHKERRQ(ierr); 1575 ierr = DMPlexSetConeOrientation(rdm, newp+5, orntNew);CHKERRQ(ierr); 1576 #if 1 1577 if ((newp+5 < cStartNew) || (newp+5 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+5, cStartNew, cEndNew); 1578 for (p = 0; p < 4; ++p) { 1579 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1580 } 1581 #endif 1582 /* C' tetrahedron: {b, f, c, a} */ 1583 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 2; 1584 orntNew[0] = -3; 1585 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 7; 1586 orntNew[1] = -2; 1587 coneNew[2] = fStartNew + (cone[0] - fStart)*4 + 3; 1588 orntNew[2] = ornt[0] < 0 ? (-(ornt[0]+1)+1)%3 : -((ornt[0]+1)%3+1); 1589 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 5; 1590 orntNew[3] = -1; 1591 ierr = DMPlexSetCone(rdm, newp+6, coneNew);CHKERRQ(ierr); 1592 ierr = DMPlexSetConeOrientation(rdm, newp+6, orntNew);CHKERRQ(ierr); 1593 #if 1 1594 if ((newp+6 < cStartNew) || (newp+6 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+6, cStartNew, cEndNew); 1595 for (p = 0; p < 4; ++p) { 1596 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1597 } 1598 #endif 1599 /* D' tetrahedron: {f, e, d, a} */ 1600 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 3; 1601 orntNew[0] = -3; 1602 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 6; 1603 orntNew[1] = -3; 1604 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 4; 1605 orntNew[2] = -2; 1606 coneNew[3] = fStartNew + (cone[1] - fStart)*4 + 3; 1607 orntNew[3] = ornt[2]; 1608 ierr = DMPlexSetCone(rdm, newp+7, coneNew);CHKERRQ(ierr); 1609 ierr = DMPlexSetConeOrientation(rdm, newp+7, orntNew);CHKERRQ(ierr); 1610 #if 1 1611 if ((newp+7 < cStartNew) || (newp+7 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+7, cStartNew, cEndNew); 1612 for (p = 0; p < 4; ++p) { 1613 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 1614 } 1615 #endif 1616 } 1617 /* Split faces have 3 edges and the same cells as the parent */ 1618 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 1619 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 1620 for (f = fStart; f < fEnd; ++f) { 1621 const PetscInt newp = fStartNew + (f - fStart)*4; 1622 const PetscInt *cone, *ornt, *support; 1623 PetscInt coneNew[3], orntNew[3], coneSize, supportSize, s; 1624 1625 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1626 ierr = DMPlexGetConeOrientation(dm, f, &ornt);CHKERRQ(ierr); 1627 /* A triangle */ 1628 coneNew[0] = eStartNew + (cone[0] - eStart)*2 + (ornt[0] < 0 ? 1 : 0); 1629 orntNew[0] = ornt[0]; 1630 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 2; 1631 orntNew[1] = -2; 1632 coneNew[2] = eStartNew + (cone[2] - eStart)*2 + (ornt[2] < 0 ? 0 : 1); 1633 orntNew[2] = ornt[2]; 1634 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1635 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1636 #if 1 1637 if ((newp+0 < fStartNew) || (newp+0 >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp+0, fStartNew, fEndNew); 1638 for (p = 0; p < 3; ++p) { 1639 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1640 } 1641 #endif 1642 /* B triangle */ 1643 coneNew[0] = eStartNew + (cone[0] - eStart)*2 + (ornt[0] < 0 ? 0 : 1); 1644 orntNew[0] = ornt[0]; 1645 coneNew[1] = eStartNew + (cone[1] - eStart)*2 + (ornt[1] < 0 ? 1 : 0); 1646 orntNew[1] = ornt[1]; 1647 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 0; 1648 orntNew[2] = -2; 1649 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1650 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1651 #if 1 1652 if ((newp+1 < fStartNew) || (newp+1 >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp+1, fStartNew, fEndNew); 1653 for (p = 0; p < 3; ++p) { 1654 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1655 } 1656 #endif 1657 /* C triangle */ 1658 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 1; 1659 orntNew[0] = -2; 1660 coneNew[1] = eStartNew + (cone[1] - eStart)*2 + (ornt[1] < 0 ? 0 : 1); 1661 orntNew[1] = ornt[1]; 1662 coneNew[2] = eStartNew + (cone[2] - eStart)*2 + (ornt[2] < 0 ? 1 : 0); 1663 orntNew[2] = ornt[2]; 1664 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 1665 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 1666 #if 1 1667 if ((newp+2 < fStartNew) || (newp+2 >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp+2, fStartNew, fEndNew); 1668 for (p = 0; p < 3; ++p) { 1669 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1670 } 1671 #endif 1672 /* D triangle */ 1673 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 0; 1674 orntNew[0] = 0; 1675 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 1; 1676 orntNew[1] = 0; 1677 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 2; 1678 orntNew[2] = 0; 1679 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 1680 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 1681 #if 1 1682 if ((newp+3 < fStartNew) || (newp+3 >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp+3, fStartNew, fEndNew); 1683 for (p = 0; p < 3; ++p) { 1684 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1685 } 1686 #endif 1687 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1688 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1689 for (r = 0; r < 4; ++r) { 1690 for (s = 0; s < supportSize; ++s) { 1691 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1692 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1693 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1694 for (c = 0; c < coneSize; ++c) { 1695 if (cone[c] == f) break; 1696 } 1697 supportRef[s] = cStartNew + (support[s] - cStart)*8 + (r==3 ? (c+2)%4 + 4 : (ornt[c] < 0 ? faces[c*3+(-(ornt[c]+1)+1+3-r)%3] : faces[c*3+r])); 1698 } 1699 ierr = DMPlexSetSupport(rdm, newp+r, supportRef);CHKERRQ(ierr); 1700 #if 1 1701 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1702 for (p = 0; p < supportSize; ++p) { 1703 if ((supportRef[p] < cStartNew) || (supportRef[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportRef[p], cStartNew, cEndNew); 1704 } 1705 #endif 1706 } 1707 } 1708 /* Interior faces have 3 edges and 2 cells */ 1709 for (c = cStart; c < cEnd; ++c) { 1710 PetscInt newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8; 1711 const PetscInt *cone, *ornt; 1712 PetscInt coneNew[3], orntNew[3]; 1713 PetscInt supportNew[2]; 1714 1715 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1716 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1717 /* Face A: {c, a, d} */ 1718 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+0)%3 : (ornt[0]+2)%3); 1719 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1720 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+0)%3 : (ornt[1]+2)%3); 1721 orntNew[1] = ornt[1] < 0 ? -2 : 0; 1722 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+0)%3 : (ornt[2]+2)%3); 1723 orntNew[2] = ornt[2] < 0 ? -2 : 0; 1724 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1725 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1726 #if 1 1727 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1728 for (p = 0; p < 3; ++p) { 1729 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1730 } 1731 #endif 1732 supportNew[0] = (c - cStart)*8 + 0; 1733 supportNew[1] = (c - cStart)*8 + 0+4; 1734 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1735 #if 1 1736 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1737 for (p = 0; p < 2; ++p) { 1738 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1739 } 1740 #endif 1741 ++newp; 1742 /* Face B: {a, b, e} */ 1743 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+2)%3 : (ornt[0]+0)%3); 1744 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1745 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+2)%3 : (ornt[3]+0)%3); 1746 orntNew[1] = ornt[3] < 0 ? -2 : 0; 1747 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+1)%3 : (ornt[1]+1)%3); 1748 orntNew[2] = ornt[1] < 0 ? -2 : 0; 1749 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1750 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1751 #if 1 1752 if ((newp+1 < fStartNew) || (newp+1 >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp+1, fStartNew, fEndNew); 1753 for (p = 0; p < 3; ++p) { 1754 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1755 } 1756 #endif 1757 supportNew[0] = (c - cStart)*8 + 1; 1758 supportNew[1] = (c - cStart)*8 + 1+4; 1759 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1760 #if 1 1761 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1762 for (p = 0; p < 2; ++p) { 1763 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1764 } 1765 #endif 1766 ++newp; 1767 /* Face C: {c, f, b} */ 1768 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+2)%3 : (ornt[2]+0)%3); 1769 orntNew[0] = ornt[2] < 0 ? -2 : 0; 1770 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+0)%3 : (ornt[3]+2)%3); 1771 orntNew[1] = ornt[3] < 0 ? -2 : 0; 1772 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+1)%3 : (ornt[0]+1)%3); 1773 orntNew[2] = ornt[0] < 0 ? -2 : 0; 1774 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1775 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1776 #if 1 1777 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1778 for (p = 0; p < 3; ++p) { 1779 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1780 } 1781 #endif 1782 supportNew[0] = (c - cStart)*8 + 2; 1783 supportNew[1] = (c - cStart)*8 + 2+4; 1784 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1785 #if 1 1786 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1787 for (p = 0; p < 2; ++p) { 1788 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1789 } 1790 #endif 1791 ++newp; 1792 /* Face D: {d, e, f} */ 1793 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+2)%3 : (ornt[1]+0)%3); 1794 orntNew[0] = ornt[1] < 0 ? -2 : 0; 1795 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+1)%3 : (ornt[3]+1)%3); 1796 orntNew[1] = ornt[3] < 0 ? -2 : 0; 1797 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+1)%3 : (ornt[2]+1)%3); 1798 orntNew[2] = ornt[2] < 0 ? -2 : 0; 1799 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1800 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1801 #if 1 1802 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1803 for (p = 0; p < 3; ++p) { 1804 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1805 } 1806 #endif 1807 supportNew[0] = (c - cStart)*8 + 3; 1808 supportNew[1] = (c - cStart)*8 + 3+4; 1809 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1810 #if 1 1811 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1812 for (p = 0; p < 2; ++p) { 1813 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1814 } 1815 #endif 1816 ++newp; 1817 /* Face E: {d, f, a} */ 1818 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+1)%3 : (ornt[2]+1)%3); 1819 orntNew[0] = ornt[2] < 0 ? 0 : -2; 1820 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1821 orntNew[1] = 0; 1822 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+0)%3 : (ornt[1]+2)%3); 1823 orntNew[2] = ornt[1] < 0 ? -2 : 0; 1824 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1825 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1826 #if 1 1827 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1828 for (p = 0; p < 3; ++p) { 1829 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1830 } 1831 #endif 1832 supportNew[0] = (c - cStart)*8 + 0+4; 1833 supportNew[1] = (c - cStart)*8 + 3+4; 1834 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1835 #if 1 1836 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1837 for (p = 0; p < 2; ++p) { 1838 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1839 } 1840 #endif 1841 ++newp; 1842 /* Face F: {c, a, f} */ 1843 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+0)%3 : (ornt[0]+2)%3); 1844 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1845 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1846 orntNew[1] = -2; 1847 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+2)%3 : (ornt[2]+0)%3); 1848 orntNew[2] = ornt[1] < 0 ? 0 : -2; 1849 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1850 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1851 #if 1 1852 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1853 for (p = 0; p < 3; ++p) { 1854 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1855 } 1856 #endif 1857 supportNew[0] = (c - cStart)*8 + 0+4; 1858 supportNew[1] = (c - cStart)*8 + 2+4; 1859 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1860 #if 1 1861 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1862 for (p = 0; p < 2; ++p) { 1863 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1864 } 1865 #endif 1866 ++newp; 1867 /* Face G: {e, a, f} */ 1868 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+1)%3 : (ornt[1]+1)%3); 1869 orntNew[0] = ornt[1] < 0 ? -2 : 0; 1870 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1871 orntNew[1] = -2; 1872 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+1)%3 : (ornt[3]+1)%3); 1873 orntNew[2] = ornt[3] < 0 ? 0 : -2; 1874 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1875 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1876 #if 1 1877 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1878 for (p = 0; p < 3; ++p) { 1879 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1880 } 1881 #endif 1882 supportNew[0] = (c - cStart)*8 + 1+4; 1883 supportNew[1] = (c - cStart)*8 + 3+4; 1884 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1885 #if 1 1886 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1887 for (p = 0; p < 2; ++p) { 1888 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1889 } 1890 #endif 1891 ++newp; 1892 /* Face H: {a, b, f} */ 1893 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+2)%3 : (ornt[0]+0)%3); 1894 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1895 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+0)%3 : (ornt[3]+2)%3); 1896 orntNew[1] = ornt[3] < 0 ? 0 : -2; 1897 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1898 orntNew[2] = 0; 1899 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1900 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1901 #if 1 1902 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1903 for (p = 0; p < 3; ++p) { 1904 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1905 } 1906 #endif 1907 supportNew[0] = (c - cStart)*8 + 1+4; 1908 supportNew[1] = (c - cStart)*8 + 2+4; 1909 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1910 #if 1 1911 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1912 for (p = 0; p < 2; ++p) { 1913 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 1914 } 1915 #endif 1916 ++newp; 1917 } 1918 /* Split Edges have 2 vertices and the same faces as the parent */ 1919 for (e = eStart; e < eEnd; ++e) { 1920 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 1921 1922 for (r = 0; r < 2; ++r) { 1923 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 1924 const PetscInt *cone, *ornt, *support; 1925 PetscInt coneNew[2], coneSize, c, supportSize, s; 1926 1927 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 1928 coneNew[0] = vStartNew + (cone[0] - vStart); 1929 coneNew[1] = vStartNew + (cone[1] - vStart); 1930 coneNew[(r+1)%2] = newv; 1931 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1932 #if 1 1933 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1934 for (p = 0; p < 2; ++p) { 1935 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 1936 } 1937 #endif 1938 ierr = DMPlexGetSupportSize(dm, e, &supportSize);CHKERRQ(ierr); 1939 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 1940 for (s = 0; s < supportSize; ++s) { 1941 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1942 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1943 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1944 for (c = 0; c < coneSize; ++c) { 1945 if (cone[c] == e) break; 1946 } 1947 supportRef[s] = fStartNew + (support[s] - fStart)*4 + (c + (ornt[c] < 0 ? 1-r : r))%3; 1948 } 1949 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1950 #if 1 1951 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1952 for (p = 0; p < supportSize; ++p) { 1953 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 1954 } 1955 #endif 1956 } 1957 } 1958 /* Face edges have 2 vertices and 2+cells*(1/2) faces */ 1959 for (f = fStart; f < fEnd; ++f) { 1960 const PetscInt *cone, *ornt, *support; 1961 PetscInt coneSize, supportSize, s; 1962 1963 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1964 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1965 for (r = 0; r < 3; ++r) { 1966 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + r; 1967 PetscInt coneNew[2], intFaces = 0, er, eint[4] = {1, 0, 2, 0}; 1968 PetscInt fint[24] = { 1, 7, -1, -1, 0, 5, 1969 -1, -1, 1, 6, 0, 4, 1970 2, 5, 3, 4, -1, -1, 1971 -1, -1, 3, 6, 2, 7}; 1972 1973 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1974 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[(r+0)%3] - eStart); 1975 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - eStart); 1976 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1977 #if 1 1978 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1979 for (p = 0; p < 2; ++p) { 1980 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 1981 } 1982 #endif 1983 supportRef[0] = fStartNew + (f - fStart)*4 + (r+1)%3; 1984 supportRef[1] = fStartNew + (f - fStart)*4 + 3; 1985 for (s = 0; s < supportSize; ++s) { 1986 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1987 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1988 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1989 for (c = 0; c < coneSize; ++c) {if (cone[c] == f) break;} 1990 /* Here we want to determine whether edge newp contains a vertex which is part of the cross-tet edge */ 1991 er = ornt[c] < 0 ? (-(ornt[c]+1) + 2-r)%3 : (ornt[c] + r)%3; 1992 if (er == eint[c]) { 1993 supportRef[2+intFaces++] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*8 + (c + 2)%4; 1994 } else { 1995 supportRef[2+intFaces++] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*8 + fint[(c*3 + er)*2 + 0]; 1996 supportRef[2+intFaces++] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*8 + fint[(c*3 + er)*2 + 1]; 1997 } 1998 } 1999 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2000 #if 1 2001 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2002 for (p = 0; p < intFaces; ++p) { 2003 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 2004 } 2005 #endif 2006 } 2007 } 2008 /* Interior edges have 2 vertices and 4 faces */ 2009 for (c = cStart; c < cEnd; ++c) { 2010 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 2011 const PetscInt *cone, *ornt, *fcone; 2012 PetscInt coneNew[2], supportNew[4], find; 2013 2014 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2015 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 2016 ierr = DMPlexGetCone(dm, cone[0], &fcone);CHKERRQ(ierr); 2017 find = GetTriEdge_Static(ornt[0], 0); 2018 coneNew[0] = vStartNew + (vEnd - vStart) + (fcone[find] - eStart); 2019 ierr = DMPlexGetCone(dm, cone[2], &fcone);CHKERRQ(ierr); 2020 find = GetTriEdge_Static(ornt[2], 1); 2021 coneNew[1] = vStartNew + (vEnd - vStart) + (fcone[find] - eStart); 2022 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2023 #if 1 2024 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2025 for (p = 0; p < 2; ++p) { 2026 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 2027 } 2028 #endif 2029 supportNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 4; 2030 supportNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 5; 2031 supportNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 6; 2032 supportNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 7; 2033 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 2034 #if 1 2035 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2036 for (p = 0; p < 4; ++p) { 2037 if ((supportNew[p] < fStartNew) || (supportNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportNew[p], fStartNew, fEndNew); 2038 } 2039 #endif 2040 } 2041 /* Old vertices have identical supports */ 2042 for (v = vStart; v < vEnd; ++v) { 2043 const PetscInt newp = vStartNew + (v - vStart); 2044 const PetscInt *support, *cone; 2045 PetscInt size, s; 2046 2047 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 2048 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 2049 for (s = 0; s < size; ++s) { 2050 PetscInt r = 0; 2051 2052 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2053 if (cone[1] == v) r = 1; 2054 supportRef[s] = eStartNew + (support[s] - eStart)*2 + r; 2055 } 2056 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2057 #if 1 2058 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2059 for (p = 0; p < size; ++p) { 2060 if ((supportRef[p] < eStartNew) || (supportRef[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", supportRef[p], eStartNew, eEndNew); 2061 } 2062 #endif 2063 } 2064 /* Edge vertices have 2 + face*2 + 0/1 supports */ 2065 for (e = eStart; e < eEnd; ++e) { 2066 const PetscInt newp = vStartNew + (vEnd - vStart) + (e - eStart); 2067 const PetscInt *cone, *support; 2068 PetscInt *star = NULL, starSize, cellSize = 0, coneSize, size, s; 2069 2070 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 2071 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 2072 supportRef[0] = eStartNew + (e - eStart)*2 + 0; 2073 supportRef[1] = eStartNew + (e - eStart)*2 + 1; 2074 for (s = 0; s < size; ++s) { 2075 PetscInt r = 0; 2076 2077 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 2078 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2079 for (r = 0; r < coneSize; ++r) {if (cone[r] == e) break;} 2080 supportRef[2+s*2+0] = eStartNew + (eEnd - eStart)*2 + (support[s] - fStart)*3 + (r+0)%3; 2081 supportRef[2+s*2+1] = eStartNew + (eEnd - eStart)*2 + (support[s] - fStart)*3 + (r+2)%3; 2082 } 2083 ierr = DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2084 for (s = 0; s < starSize*2; s += 2) { 2085 const PetscInt *cone, *ornt; 2086 PetscInt e01, e23; 2087 2088 if ((star[s] >= cStart) && (star[s] < cEnd)) { 2089 /* Check edge 0-1 */ 2090 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 2091 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 2092 ierr = DMPlexGetCone(dm, cone[0], &cone);CHKERRQ(ierr); 2093 e01 = cone[GetTriEdge_Static(ornt[0], 0)]; 2094 /* Check edge 2-3 */ 2095 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 2096 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 2097 ierr = DMPlexGetCone(dm, cone[2], &cone);CHKERRQ(ierr); 2098 e23 = cone[GetTriEdge_Static(ornt[2], 1)]; 2099 if ((e01 == e) || (e23 == e)) {supportRef[2+size*2+cellSize++] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (star[s] - cStart);} 2100 } 2101 } 2102 ierr = DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2103 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2104 #if 1 2105 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2106 for (p = 0; p < 2+size*2+cellSize; ++p) { 2107 if ((supportRef[p] < eStartNew) || (supportRef[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", supportRef[p], eStartNew, eEndNew); 2108 } 2109 #endif 2110 } 2111 ierr = PetscFree(supportRef);CHKERRQ(ierr); 2112 ierr = DMPlexRestoreFaces_Internal(dm, 3, cStart, NULL, NULL, &faces);CHKERRQ(ierr); 2113 break; 2114 case 6: 2115 /* Hex 3D */ 2116 /* 2117 Bottom (viewed from top) Top 2118 1---------2---------2 7---------2---------6 2119 | | | | | | 2120 | B 2 C | | H 2 G | 2121 | | | | | | 2122 3----3----0----1----1 3----3----0----1----1 2123 | | | | | | 2124 | A 0 D | | E 0 F | 2125 | | | | | | 2126 0---------0---------3 4---------0---------5 2127 */ 2128 /* All cells have 6 faces: Bottom, Top, Front, Back, Right, Left */ 2129 for (c = cStart; c < cEnd; ++c) { 2130 const PetscInt newp = (c - cStart)*8; 2131 const PetscInt *cone, *ornt; 2132 PetscInt coneNew[6], orntNew[6]; 2133 2134 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2135 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 2136 /* A hex */ 2137 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetQuadSubface_Static(ornt[0], 0); 2138 orntNew[0] = ornt[0]; 2139 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 8; /* AE */ 2140 orntNew[1] = 0; 2141 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetQuadSubface_Static(ornt[2], 0); 2142 orntNew[2] = ornt[2]; 2143 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 3; /* AB */ 2144 orntNew[3] = 0; 2145 coneNew[4] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 0; /* AD */ 2146 orntNew[4] = 0; 2147 coneNew[5] = fStartNew + (cone[5] - fStart)*4 + GetQuadSubface_Static(ornt[5], 0); 2148 orntNew[5] = ornt[5]; 2149 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 2150 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 2151 #if 1 2152 if ((newp+0 < cStartNew) || (newp+0 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+0, cStartNew, cEndNew); 2153 for (p = 0; p < 6; ++p) { 2154 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 2155 } 2156 #endif 2157 /* B hex */ 2158 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetQuadSubface_Static(ornt[0], 1); 2159 orntNew[0] = ornt[0]; 2160 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 11; /* BH */ 2161 orntNew[1] = 0; 2162 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 3; /* AB */ 2163 orntNew[2] = -3; 2164 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetQuadSubface_Static(ornt[3], 1); 2165 orntNew[3] = ornt[3]; 2166 coneNew[4] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 2; /* BC */ 2167 orntNew[4] = 0; 2168 coneNew[5] = fStartNew + (cone[5] - fStart)*4 + GetQuadSubface_Static(ornt[5], 3); 2169 orntNew[5] = ornt[5]; 2170 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 2171 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 2172 #if 1 2173 if ((newp+1 < cStartNew) || (newp+1 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+1, cStartNew, cEndNew); 2174 for (p = 0; p < 6; ++p) { 2175 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 2176 } 2177 #endif 2178 /* C hex */ 2179 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetQuadSubface_Static(ornt[0], 2); 2180 orntNew[0] = ornt[0]; 2181 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 10; /* CG */ 2182 orntNew[1] = 0; 2183 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 1; /* CD */ 2184 orntNew[2] = 0; 2185 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetQuadSubface_Static(ornt[3], 0); 2186 orntNew[3] = ornt[3]; 2187 coneNew[4] = fStartNew + (cone[4] - fStart)*4 + GetQuadSubface_Static(ornt[4], 1); 2188 orntNew[4] = ornt[4]; 2189 coneNew[5] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 2; /* BC */ 2190 orntNew[5] = -3; 2191 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 2192 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 2193 #if 1 2194 if ((newp+2 < cStartNew) || (newp+2 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+2, cStartNew, cEndNew); 2195 for (p = 0; p < 6; ++p) { 2196 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 2197 } 2198 #endif 2199 /* D hex */ 2200 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetQuadSubface_Static(ornt[0], 3); 2201 orntNew[0] = ornt[0]; 2202 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 9; /* DF */ 2203 orntNew[1] = 0; 2204 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetQuadSubface_Static(ornt[2], 1); 2205 orntNew[2] = ornt[2]; 2206 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 1; /* CD */ 2207 orntNew[3] = -3; 2208 coneNew[4] = fStartNew + (cone[4] - fStart)*4 + GetQuadSubface_Static(ornt[4], 0); 2209 orntNew[4] = ornt[4]; 2210 coneNew[5] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 0; /* AD */ 2211 orntNew[5] = -3; 2212 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 2213 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 2214 #if 1 2215 if ((newp+3 < cStartNew) || (newp+3 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+3, cStartNew, cEndNew); 2216 for (p = 0; p < 6; ++p) { 2217 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 2218 } 2219 #endif 2220 /* E hex */ 2221 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 8; /* AE */ 2222 orntNew[0] = -3; 2223 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetQuadSubface_Static(ornt[1], 0); 2224 orntNew[1] = ornt[1]; 2225 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetQuadSubface_Static(ornt[2], 3); 2226 orntNew[2] = ornt[2]; 2227 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 7; /* EH */ 2228 orntNew[3] = 0; 2229 coneNew[4] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 4; /* EF */ 2230 orntNew[4] = 0; 2231 coneNew[5] = fStartNew + (cone[5] - fStart)*4 + GetQuadSubface_Static(ornt[5], 1); 2232 orntNew[5] = ornt[5]; 2233 ierr = DMPlexSetCone(rdm, newp+4, coneNew);CHKERRQ(ierr); 2234 ierr = DMPlexSetConeOrientation(rdm, newp+4, orntNew);CHKERRQ(ierr); 2235 #if 1 2236 if ((newp+4 < cStartNew) || (newp+4 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+4, cStartNew, cEndNew); 2237 for (p = 0; p < 6; ++p) { 2238 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 2239 } 2240 #endif 2241 /* F hex */ 2242 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 9; /* DF */ 2243 orntNew[0] = -3; 2244 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetQuadSubface_Static(ornt[1], 1); 2245 orntNew[1] = ornt[1]; 2246 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetQuadSubface_Static(ornt[2], 2); 2247 orntNew[2] = ornt[2]; 2248 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 5; /* FG */ 2249 orntNew[3] = 0; 2250 coneNew[4] = fStartNew + (cone[4] - fStart)*4 + GetQuadSubface_Static(ornt[4], 3); 2251 orntNew[4] = ornt[4]; 2252 coneNew[5] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 4; /* EF */ 2253 orntNew[5] = -3; 2254 ierr = DMPlexSetCone(rdm, newp+5, coneNew);CHKERRQ(ierr); 2255 ierr = DMPlexSetConeOrientation(rdm, newp+5, orntNew);CHKERRQ(ierr); 2256 #if 1 2257 if ((newp+5 < cStartNew) || (newp+5 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+5, cStartNew, cEndNew); 2258 for (p = 0; p < 6; ++p) { 2259 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 2260 } 2261 #endif 2262 /* G hex */ 2263 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 10; /* CG */ 2264 orntNew[0] = -3; 2265 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetQuadSubface_Static(ornt[1], 2); 2266 orntNew[1] = ornt[1]; 2267 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 5; /* FG */ 2268 orntNew[2] = -3; 2269 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetQuadSubface_Static(ornt[3], 3); 2270 orntNew[3] = ornt[3]; 2271 coneNew[4] = fStartNew + (cone[4] - fStart)*4 + GetQuadSubface_Static(ornt[4], 2); 2272 orntNew[4] = ornt[4]; 2273 coneNew[5] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 6; /* GH */ 2274 orntNew[5] = 0; 2275 ierr = DMPlexSetCone(rdm, newp+6, coneNew);CHKERRQ(ierr); 2276 ierr = DMPlexSetConeOrientation(rdm, newp+6, orntNew);CHKERRQ(ierr); 2277 #if 1 2278 if ((newp+6 < cStartNew) || (newp+6 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+6, cStartNew, cEndNew); 2279 for (p = 0; p < 6; ++p) { 2280 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 2281 } 2282 #endif 2283 /* H hex */ 2284 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 11; /* BH */ 2285 orntNew[0] = -3; 2286 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetQuadSubface_Static(ornt[1], 3); 2287 orntNew[1] = ornt[1]; 2288 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 7; /* EH */ 2289 orntNew[2] = -3; 2290 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetQuadSubface_Static(ornt[3], 2); 2291 orntNew[3] = ornt[3]; 2292 coneNew[4] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 6; /* GH */ 2293 orntNew[4] = -3; 2294 coneNew[5] = fStartNew + (cone[5] - fStart)*4 + GetQuadSubface_Static(ornt[5], 2); 2295 orntNew[5] = ornt[5]; 2296 ierr = DMPlexSetCone(rdm, newp+7, coneNew);CHKERRQ(ierr); 2297 ierr = DMPlexSetConeOrientation(rdm, newp+7, orntNew);CHKERRQ(ierr); 2298 #if 1 2299 if ((newp+7 < cStartNew) || (newp+7 >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", newp+7, cStartNew, cEndNew); 2300 for (p = 0; p < 6; ++p) { 2301 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 2302 } 2303 #endif 2304 } 2305 /* Split faces have 4 edges and the same cells as the parent */ 2306 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 2307 ierr = PetscMalloc((4 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 2308 for (f = fStart; f < fEnd; ++f) { 2309 for (r = 0; r < 4; ++r) { 2310 /* TODO: This can come from GetFaces_Internal() */ 2311 const PetscInt newCells[24] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 3, 5, 4, 2, 1, 7, 6, 3, 2, 6, 5, 0, 4, 7, 1}; 2312 const PetscInt newp = fStartNew + (f - fStart)*4 + r; 2313 const PetscInt *cone, *ornt, *support; 2314 PetscInt coneNew[4], orntNew[4], coneSize, c, supportSize, s; 2315 2316 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 2317 ierr = DMPlexGetConeOrientation(dm, f, &ornt);CHKERRQ(ierr); 2318 coneNew[0] = eStartNew + (cone[(r+3)%4] - eStart)*2 + (ornt[(r+3)%4] < 0 ? 0 : 1); 2319 orntNew[0] = ornt[(r+3)%4]; 2320 coneNew[1] = eStartNew + (cone[r] - eStart)*2 + (ornt[r] < 0 ? 1 : 0); 2321 orntNew[1] = ornt[r]; 2322 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + r; 2323 orntNew[2] = 0; 2324 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + (r+3)%4; 2325 orntNew[3] = -2; 2326 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2327 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2328 #if 1 2329 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2330 for (p = 0; p < 4; ++p) { 2331 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2332 } 2333 #endif 2334 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 2335 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 2336 for (s = 0; s < supportSize; ++s) { 2337 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 2338 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2339 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 2340 for (c = 0; c < coneSize; ++c) { 2341 if (cone[c] == f) break; 2342 } 2343 supportRef[s] = cStartNew + (support[s] - cStart)*8 + newCells[c*4+GetQuadSubface_Static(ornt[c], r)]; 2344 } 2345 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2346 #if 1 2347 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2348 for (p = 0; p < supportSize; ++p) { 2349 if ((supportRef[p] < cStartNew) || (supportRef[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportRef[p], cStartNew, cEndNew); 2350 } 2351 #endif 2352 } 2353 } 2354 /* Interior faces have 4 edges and 2 cells */ 2355 for (c = cStart; c < cEnd; ++c) { 2356 const PetscInt newCells[24] = {0, 3, 2, 3, 1, 2, 0, 1, 4, 5, 5, 6, 6, 7, 4, 7, 0, 4, 3, 5, 2, 6, 1, 7}; 2357 const PetscInt *cone, *ornt; 2358 PetscInt newp, coneNew[4], orntNew[4], supportNew[2]; 2359 2360 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2361 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 2362 /* A-D face */ 2363 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 0; 2364 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*4 + GetQuadEdge_Static(ornt[2], 0); 2365 orntNew[0] = -2; 2366 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*4 + GetQuadEdge_Static(ornt[0], 3); 2367 orntNew[1] = 0; 2368 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 0; 2369 orntNew[2] = 0; 2370 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 2; 2371 orntNew[3] = -2; 2372 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2373 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2374 #if 1 2375 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2376 for (p = 0; p < 4; ++p) { 2377 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2378 } 2379 #endif 2380 /* C-D face */ 2381 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 1; 2382 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[4] - fStart)*4 + GetQuadEdge_Static(ornt[4], 0); 2383 orntNew[0] = -2; 2384 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*4 + GetQuadEdge_Static(ornt[0], 2); 2385 orntNew[1] = 0; 2386 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 0; 2387 orntNew[2] = 0; 2388 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 4; 2389 orntNew[3] = -2; 2390 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2391 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2392 #if 1 2393 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2394 for (p = 0; p < 4; ++p) { 2395 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2396 } 2397 #endif 2398 /* B-C face */ 2399 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 2; 2400 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*4 + GetQuadEdge_Static(ornt[0], 1); 2401 orntNew[0] = -2; 2402 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*4 + GetQuadEdge_Static(ornt[3], 0); 2403 orntNew[1] = 0; 2404 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 3; 2405 orntNew[2] = 0; 2406 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 0; 2407 orntNew[3] = -2; 2408 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2409 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2410 #if 1 2411 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2412 for (p = 0; p < 4; ++p) { 2413 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2414 } 2415 #endif 2416 /* A-B face */ 2417 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 3; 2418 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*4 + GetQuadEdge_Static(ornt[0], 0); 2419 orntNew[0] = -2; 2420 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[5] - fStart)*4 + GetQuadEdge_Static(ornt[5], 3); 2421 orntNew[1] = 0; 2422 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 5; 2423 orntNew[2] = 0; 2424 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 0; 2425 orntNew[3] = -2; 2426 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2427 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2428 #if 1 2429 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2430 for (p = 0; p < 4; ++p) { 2431 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2432 } 2433 #endif 2434 /* E-F face */ 2435 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 4; 2436 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*4 + GetQuadEdge_Static(ornt[2], 2); 2437 orntNew[0] = -2; 2438 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*4 + GetQuadEdge_Static(ornt[1], 0); 2439 orntNew[1] = 0; 2440 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 1; 2441 orntNew[2] = 0; 2442 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 2; 2443 orntNew[3] = -2; 2444 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2445 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2446 #if 1 2447 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2448 for (p = 0; p < 4; ++p) { 2449 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2450 } 2451 #endif 2452 /* F-G face */ 2453 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 5; 2454 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[4] - fStart)*4 + GetQuadEdge_Static(ornt[4], 2); 2455 orntNew[0] = -2; 2456 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*4 + GetQuadEdge_Static(ornt[1], 1); 2457 orntNew[1] = 0; 2458 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 1; 2459 orntNew[2] = 0; 2460 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 4; 2461 orntNew[3] = -2; 2462 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2463 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2464 #if 1 2465 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2466 for (p = 0; p < 4; ++p) { 2467 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2468 } 2469 #endif 2470 /* G-H face */ 2471 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 6; 2472 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*4 + GetQuadEdge_Static(ornt[3], 2); 2473 orntNew[0] = -2; 2474 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*4 + GetQuadEdge_Static(ornt[1], 2); 2475 orntNew[1] = 0; 2476 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 1; 2477 orntNew[2] = 0; 2478 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 3; 2479 orntNew[3] = -2; 2480 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2481 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2482 #if 1 2483 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2484 for (p = 0; p < 4; ++p) { 2485 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2486 } 2487 #endif 2488 /* E-H face */ 2489 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 7; 2490 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[5] - fStart)*4 + GetQuadEdge_Static(ornt[5], 1); 2491 orntNew[0] = -2; 2492 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*4 + GetQuadEdge_Static(ornt[1], 3); 2493 orntNew[1] = 0; 2494 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 1; 2495 orntNew[2] = 0; 2496 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 5; 2497 orntNew[3] = -2; 2498 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2499 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2500 #if 1 2501 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2502 for (p = 0; p < 4; ++p) { 2503 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2504 } 2505 #endif 2506 /* A-E face */ 2507 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 8; 2508 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[5] - fStart)*4 + GetQuadEdge_Static(ornt[5], 0); 2509 orntNew[0] = -2; 2510 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*4 + GetQuadEdge_Static(ornt[2], 3); 2511 orntNew[1] = 0; 2512 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 2; 2513 orntNew[2] = 0; 2514 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 5; 2515 orntNew[3] = -2; 2516 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2517 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2518 #if 1 2519 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2520 for (p = 0; p < 4; ++p) { 2521 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2522 } 2523 #endif 2524 /* D-F face */ 2525 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 9; 2526 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*4 + GetQuadEdge_Static(ornt[2], 1); 2527 orntNew[0] = -2; 2528 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[4] - fStart)*4 + GetQuadEdge_Static(ornt[4], 3); 2529 orntNew[1] = 0; 2530 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 4; 2531 orntNew[2] = 0; 2532 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 2; 2533 orntNew[3] = -2; 2534 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2535 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2536 #if 1 2537 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2538 for (p = 0; p < 4; ++p) { 2539 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2540 } 2541 #endif 2542 /* C-G face */ 2543 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 10; 2544 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[4] - fStart)*4 + GetQuadEdge_Static(ornt[4], 1); 2545 orntNew[0] = -2; 2546 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*4 + GetQuadEdge_Static(ornt[3], 3); 2547 orntNew[1] = 0; 2548 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 3; 2549 orntNew[2] = 0; 2550 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 4; 2551 orntNew[3] = -2; 2552 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2553 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2554 #if 1 2555 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2556 for (p = 0; p < 4; ++p) { 2557 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2558 } 2559 #endif 2560 /* B-H face */ 2561 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 11; 2562 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*4 + GetQuadEdge_Static(ornt[3], 1); 2563 orntNew[0] = -2; 2564 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[5] - fStart)*4 + GetQuadEdge_Static(ornt[5], 2); 2565 orntNew[1] = 0; 2566 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 5; 2567 orntNew[2] = 0; 2568 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + 3; 2569 orntNew[3] = -2; 2570 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2571 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 2572 #if 1 2573 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2574 for (p = 0; p < 4; ++p) { 2575 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 2576 } 2577 #endif 2578 for (r = 0; r < 12; ++r) { 2579 newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + r; 2580 supportNew[0] = cStartNew + (c - cStart)*8 + newCells[r*2+0]; 2581 supportNew[1] = cStartNew + (c - cStart)*8 + newCells[r*2+1]; 2582 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 2583 #if 1 2584 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2585 for (p = 0; p < 2; ++p) { 2586 if ((supportNew[p] < cStartNew) || (supportNew[p] >= cEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a cell [%d, %d)", supportNew[p], cStartNew, cEndNew); 2587 } 2588 #endif 2589 } 2590 } 2591 /* Split edges have 2 vertices and the same faces as the parent */ 2592 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 2593 for (e = eStart; e < eEnd; ++e) { 2594 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 2595 2596 for (r = 0; r < 2; ++r) { 2597 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 2598 const PetscInt *cone, *ornt, *support; 2599 PetscInt coneNew[2], coneSize, c, supportSize, s; 2600 2601 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 2602 coneNew[0] = vStartNew + (cone[0] - vStart); 2603 coneNew[1] = vStartNew + (cone[1] - vStart); 2604 coneNew[(r+1)%2] = newv; 2605 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2606 #if 1 2607 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2608 for (p = 0; p < 2; ++p) { 2609 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 2610 } 2611 #endif 2612 ierr = DMPlexGetSupportSize(dm, e, &supportSize);CHKERRQ(ierr); 2613 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 2614 for (s = 0; s < supportSize; ++s) { 2615 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 2616 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2617 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 2618 for (c = 0; c < coneSize; ++c) { 2619 if (cone[c] == e) break; 2620 } 2621 supportRef[s] = fStartNew + (support[s] - fStart)*4 + (ornt[c] < 0 ? (c+1-r)%4 : (c+r)%4); 2622 } 2623 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2624 #if 1 2625 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2626 for (p = 0; p < supportSize; ++p) { 2627 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 2628 } 2629 #endif 2630 } 2631 } 2632 /* Face edges have 2 vertices and 2+cells faces */ 2633 for (f = fStart; f < fEnd; ++f) { 2634 const PetscInt newFaces[24] = {3, 2, 1, 0, 4, 5, 6, 7, 0, 9, 4, 8, 2, 11, 6, 10, 1, 10, 5, 9, 8, 7, 11, 3}; 2635 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 2636 const PetscInt *cone, *coneCell, *orntCell, *support; 2637 PetscInt coneNew[2], coneSize, c, supportSize, s; 2638 2639 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 2640 for (r = 0; r < 4; ++r) { 2641 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + r; 2642 2643 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - eStart); 2644 coneNew[1] = newv; 2645 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2646 #if 1 2647 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2648 for (p = 0; p < 2; ++p) { 2649 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 2650 } 2651 #endif 2652 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 2653 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 2654 supportRef[0] = fStartNew + (f - fStart)*4 + r; 2655 supportRef[1] = fStartNew + (f - fStart)*4 + (r+1)%4; 2656 for (s = 0; s < supportSize; ++s) { 2657 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 2658 ierr = DMPlexGetCone(dm, support[s], &coneCell);CHKERRQ(ierr); 2659 ierr = DMPlexGetConeOrientation(dm, support[s], &orntCell);CHKERRQ(ierr); 2660 for (c = 0; c < coneSize; ++c) if (coneCell[c] == f) break; 2661 supportRef[2+s] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*12 + newFaces[c*4 + GetQuadEdge_Static(orntCell[c], r)]; 2662 } 2663 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2664 #if 1 2665 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2666 for (p = 0; p < 2+supportSize; ++p) { 2667 if ((supportRef[p] < fStartNew) || (supportRef[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportRef[p], fStartNew, fEndNew); 2668 } 2669 #endif 2670 } 2671 } 2672 /* Cell edges have 2 vertices and 4 faces */ 2673 for (c = cStart; c < cEnd; ++c) { 2674 const PetscInt newFaces[24] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 9, 4, 8, 2, 11, 6, 10, 1, 10, 5, 9, 3, 8, 7, 11}; 2675 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (c - cStart); 2676 const PetscInt *cone; 2677 PetscInt coneNew[2], supportNew[4]; 2678 2679 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2680 for (r = 0; r < 6; ++r) { 2681 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + r; 2682 2683 coneNew[0] = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (cone[r] - fStart); 2684 coneNew[1] = newv; 2685 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2686 #if 1 2687 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2688 for (p = 0; p < 2; ++p) { 2689 if ((coneNew[p] < vStartNew) || (coneNew[p] >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", coneNew[p], vStartNew, vEndNew); 2690 } 2691 #endif 2692 for (f = 0; f < 4; ++f) supportNew[f] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + newFaces[r*4+f]; 2693 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 2694 #if 1 2695 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2696 for (p = 0; p < 4; ++p) { 2697 if ((supportNew[p] < fStartNew) || (supportNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", supportNew[p], fStartNew, fEndNew); 2698 } 2699 #endif 2700 } 2701 } 2702 /* Old vertices have identical supports */ 2703 for (v = vStart; v < vEnd; ++v) { 2704 const PetscInt newp = vStartNew + (v - vStart); 2705 const PetscInt *support, *cone; 2706 PetscInt size, s; 2707 2708 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 2709 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 2710 for (s = 0; s < size; ++s) { 2711 PetscInt r = 0; 2712 2713 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2714 if (cone[1] == v) r = 1; 2715 supportRef[s] = eStartNew + (support[s] - eStart)*2 + r; 2716 } 2717 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2718 #if 1 2719 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2720 for (p = 0; p < size; ++p) { 2721 if ((supportRef[p] < eStartNew) || (supportRef[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", supportRef[p], eStartNew, eEndNew); 2722 } 2723 #endif 2724 } 2725 /* Edge vertices have 2 + faces supports */ 2726 for (e = eStart; e < eEnd; ++e) { 2727 const PetscInt newp = vStartNew + (vEnd - vStart) + (e - eStart); 2728 const PetscInt *cone, *support; 2729 PetscInt size, s; 2730 2731 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 2732 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 2733 supportRef[0] = eStartNew + (e - eStart)*2 + 0; 2734 supportRef[1] = eStartNew + (e - eStart)*2 + 1; 2735 for (s = 0; s < size; ++s) { 2736 PetscInt r; 2737 2738 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2739 for (r = 0; r < 4; ++r) if (cone[r] == e) break; 2740 supportRef[2+s] = eStartNew + (eEnd - eStart)*2 + (support[s] - fStart)*4 + r; 2741 } 2742 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2743 #if 1 2744 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2745 for (p = 0; p < 2+size; ++p) { 2746 if ((supportRef[p] < eStartNew) || (supportRef[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", supportRef[p], eStartNew, eEndNew); 2747 } 2748 #endif 2749 } 2750 /* Face vertices have 4 + cells supports */ 2751 for (f = fStart; f < fEnd; ++f) { 2752 const PetscInt newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 2753 const PetscInt *cone, *support; 2754 PetscInt size, s; 2755 2756 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 2757 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 2758 for (r = 0; r < 4; ++r) supportRef[r] = eStartNew + (e - eStart)*2 + (f - fStart)*4 + r; 2759 for (s = 0; s < size; ++s) { 2760 PetscInt r; 2761 2762 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2763 for (r = 0; r < 6; ++r) if (cone[r] == f) break; 2764 supportRef[4+s] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (support[s] - cStart)*6 + r; 2765 } 2766 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2767 #if 1 2768 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2769 for (p = 0; p < 4+size; ++p) { 2770 if ((supportRef[p] < eStartNew) || (supportRef[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", supportRef[p], eStartNew, eEndNew); 2771 } 2772 #endif 2773 } 2774 /* Cell vertices have 6 supports */ 2775 for (c = cStart; c < cEnd; ++c) { 2776 const PetscInt newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (c - cStart); 2777 PetscInt supportNew[6]; 2778 2779 for (r = 0; r < 6; ++r) { 2780 supportNew[r] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + r; 2781 } 2782 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 2783 } 2784 ierr = PetscFree(supportRef);CHKERRQ(ierr); 2785 break; 2786 default: 2787 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2788 } 2789 PetscFunctionReturn(0); 2790 } 2791 2792 #undef __FUNCT__ 2793 #define __FUNCT__ "CellRefinerSetCoordinates" 2794 PetscErrorCode CellRefinerSetCoordinates(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 2795 { 2796 PetscSection coordSection, coordSectionNew; 2797 Vec coordinates, coordinatesNew; 2798 PetscScalar *coords, *coordsNew; 2799 PetscInt dim, depth, coordSizeNew, cStart, cEnd, c, vStart, vStartNew, vEnd, v, eStart, eEnd, eMax, e, fStart, fEnd, fMax, f; 2800 PetscErrorCode ierr; 2801 2802 PetscFunctionBegin; 2803 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2804 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2805 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2806 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 2807 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2808 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2809 ierr = DMPlexGetHybridBounds(dm, NULL, &fMax, &eMax, NULL);CHKERRQ(ierr); 2810 ierr = GetDepthStart_Private(depth, depthSize, NULL, NULL, NULL, &vStartNew);CHKERRQ(ierr); 2811 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2812 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &coordSectionNew);CHKERRQ(ierr); 2813 ierr = PetscSectionSetNumFields(coordSectionNew, 1);CHKERRQ(ierr); 2814 ierr = PetscSectionSetFieldComponents(coordSectionNew, 0, dim);CHKERRQ(ierr); 2815 ierr = PetscSectionSetChart(coordSectionNew, vStartNew, vStartNew+depthSize[0]);CHKERRQ(ierr); 2816 if (fMax < 0) fMax = fEnd; 2817 if (eMax < 0) eMax = eEnd; 2818 /* All vertices have the dim coordinates */ 2819 for (v = vStartNew; v < vStartNew+depthSize[0]; ++v) { 2820 ierr = PetscSectionSetDof(coordSectionNew, v, dim);CHKERRQ(ierr); 2821 ierr = PetscSectionSetFieldDof(coordSectionNew, v, 0, dim);CHKERRQ(ierr); 2822 } 2823 ierr = PetscSectionSetUp(coordSectionNew);CHKERRQ(ierr); 2824 ierr = DMPlexSetCoordinateSection(rdm, coordSectionNew);CHKERRQ(ierr); 2825 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2826 ierr = PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);CHKERRQ(ierr); 2827 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinatesNew);CHKERRQ(ierr); 2828 ierr = PetscObjectSetName((PetscObject) coordinatesNew, "coordinates");CHKERRQ(ierr); 2829 ierr = VecSetSizes(coordinatesNew, coordSizeNew, PETSC_DETERMINE);CHKERRQ(ierr); 2830 ierr = VecSetFromOptions(coordinatesNew);CHKERRQ(ierr); 2831 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2832 ierr = VecGetArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 2833 switch (refiner) { 2834 case 6: /* Hex 3D */ 2835 /* Face vertices have the average of corner coordinates */ 2836 for (f = fStart; f < fEnd; ++f) { 2837 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 2838 PetscInt *cone = NULL; 2839 PetscInt closureSize, coneSize = 0, off[8], offnew, p, d; 2840 2841 ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2842 for (p = 0; p < closureSize*2; p += 2) { 2843 const PetscInt point = cone[p]; 2844 if ((point >= vStart) && (point < vEnd)) cone[coneSize++] = point; 2845 } 2846 for (v = 0; v < coneSize; ++v) { 2847 ierr = PetscSectionGetOffset(coordSection, cone[v], &off[v]);CHKERRQ(ierr); 2848 } 2849 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2850 for (d = 0; d < dim; ++d) { 2851 coordsNew[offnew+d] = 0.0; 2852 for (v = 0; v < coneSize; ++v) coordsNew[offnew+d] += coords[off[v]+d]; 2853 coordsNew[offnew+d] /= coneSize; 2854 } 2855 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2856 } 2857 case 2: /* Hex 2D */ 2858 /* Cell vertices have the average of corner coordinates */ 2859 for (c = cStart; c < cEnd; ++c) { 2860 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (c - cStart) + (dim > 2 ? (fEnd - fStart) : 0); 2861 PetscInt *cone = NULL; 2862 PetscInt closureSize, coneSize = 0, off[8], offnew, p, d; 2863 2864 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2865 for (p = 0; p < closureSize*2; p += 2) { 2866 const PetscInt point = cone[p]; 2867 if ((point >= vStart) && (point < vEnd)) cone[coneSize++] = point; 2868 } 2869 for (v = 0; v < coneSize; ++v) { 2870 ierr = PetscSectionGetOffset(coordSection, cone[v], &off[v]);CHKERRQ(ierr); 2871 } 2872 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2873 for (d = 0; d < dim; ++d) { 2874 coordsNew[offnew+d] = 0.0; 2875 for (v = 0; v < coneSize; ++v) coordsNew[offnew+d] += coords[off[v]+d]; 2876 coordsNew[offnew+d] /= coneSize; 2877 } 2878 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2879 } 2880 case 1: /* Simplicial 2D */ 2881 case 3: /* Hybrid Simplicial 2D */ 2882 case 5: /* Simplicial 3D */ 2883 /* Edge vertices have the average of endpoint coordinates */ 2884 for (e = eStart; e < eMax; ++e) { 2885 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 2886 const PetscInt *cone; 2887 PetscInt coneSize, offA, offB, offnew, d; 2888 2889 ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr); 2890 if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %d cone should have two vertices, not %d", e, coneSize); 2891 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 2892 ierr = PetscSectionGetOffset(coordSection, cone[0], &offA);CHKERRQ(ierr); 2893 ierr = PetscSectionGetOffset(coordSection, cone[1], &offB);CHKERRQ(ierr); 2894 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2895 for (d = 0; d < dim; ++d) { 2896 coordsNew[offnew+d] = 0.5*(coords[offA+d] + coords[offB+d]); 2897 } 2898 } 2899 /* Old vertices have the same coordinates */ 2900 for (v = vStart; v < vEnd; ++v) { 2901 const PetscInt newv = vStartNew + (v - vStart); 2902 PetscInt off, offnew, d; 2903 2904 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 2905 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2906 for (d = 0; d < dim; ++d) { 2907 coordsNew[offnew+d] = coords[off+d]; 2908 } 2909 } 2910 break; 2911 default: 2912 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2913 } 2914 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2915 ierr = VecRestoreArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 2916 ierr = DMSetCoordinatesLocal(rdm, coordinatesNew);CHKERRQ(ierr); 2917 ierr = VecDestroy(&coordinatesNew);CHKERRQ(ierr); 2918 ierr = PetscSectionDestroy(&coordSectionNew);CHKERRQ(ierr); 2919 PetscFunctionReturn(0); 2920 } 2921 2922 #undef __FUNCT__ 2923 #define __FUNCT__ "DMPlexCreateProcessSF" 2924 PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess) 2925 { 2926 PetscInt numRoots, numLeaves, l; 2927 const PetscInt *localPoints; 2928 const PetscSFNode *remotePoints; 2929 PetscInt *localPointsNew; 2930 PetscSFNode *remotePointsNew; 2931 PetscInt *ranks, *ranksNew; 2932 PetscErrorCode ierr; 2933 2934 PetscFunctionBegin; 2935 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2936 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranks);CHKERRQ(ierr); 2937 for (l = 0; l < numLeaves; ++l) { 2938 ranks[l] = remotePoints[l].rank; 2939 } 2940 ierr = PetscSortRemoveDupsInt(&numLeaves, ranks);CHKERRQ(ierr); 2941 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranksNew);CHKERRQ(ierr); 2942 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 2943 ierr = PetscMalloc(numLeaves * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 2944 for (l = 0; l < numLeaves; ++l) { 2945 ranksNew[l] = ranks[l]; 2946 localPointsNew[l] = l; 2947 remotePointsNew[l].index = 0; 2948 remotePointsNew[l].rank = ranksNew[l]; 2949 } 2950 ierr = PetscFree(ranks);CHKERRQ(ierr); 2951 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);CHKERRQ(ierr); 2952 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);CHKERRQ(ierr); 2953 ierr = PetscSFSetFromOptions(*sfProcess);CHKERRQ(ierr); 2954 ierr = PetscSFSetGraph(*sfProcess, 1, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2955 PetscFunctionReturn(0); 2956 } 2957 2958 #undef __FUNCT__ 2959 #define __FUNCT__ "CellRefinerCreateSF" 2960 PetscErrorCode CellRefinerCreateSF(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 2961 { 2962 PetscSF sf, sfNew, sfProcess; 2963 IS processRanks; 2964 MPI_Datatype depthType; 2965 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 2966 const PetscInt *localPoints, *neighbors; 2967 const PetscSFNode *remotePoints; 2968 PetscInt *localPointsNew; 2969 PetscSFNode *remotePointsNew; 2970 PetscInt *depthSizeOld, *rdepthSize, *rdepthSizeOld, *rdepthMaxOld, *rvStart, *rvStartNew, *reStart, *reStartNew, *rfStart, *rfStartNew, *rcStart, *rcStartNew; 2971 PetscInt depth, numNeighbors, pStartNew, pEndNew, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eStartNew, eEnd, eMax, r, n; 2972 PetscErrorCode ierr; 2973 2974 PetscFunctionBegin; 2975 ierr = DMPlexGetChart(rdm, &pStartNew, &pEndNew);CHKERRQ(ierr); 2976 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2977 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2978 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 2979 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2980 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2981 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 2982 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 2983 switch (refiner) { 2984 case 3: 2985 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 2986 cMax = PetscMin(cEnd, cMax); 2987 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 2988 fMax = PetscMin(fEnd, fMax); 2989 } 2990 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2991 ierr = DMGetPointSF(rdm, &sfNew);CHKERRQ(ierr); 2992 /* Caculate size of new SF */ 2993 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2994 if (numRoots < 0) PetscFunctionReturn(0); 2995 for (l = 0; l < numLeaves; ++l) { 2996 const PetscInt p = localPoints[l]; 2997 2998 switch (refiner) { 2999 case 1: 3000 /* Simplicial 2D */ 3001 if ((p >= vStart) && (p < vEnd)) { 3002 /* Old vertices stay the same */ 3003 ++numLeavesNew; 3004 } else if ((p >= fStart) && (p < fEnd)) { 3005 /* Old faces add new faces and vertex */ 3006 numLeavesNew += 2 + 1; 3007 } else if ((p >= cStart) && (p < cEnd)) { 3008 /* Old cells add new cells and interior faces */ 3009 numLeavesNew += 4 + 3; 3010 } 3011 break; 3012 case 2: 3013 /* Hex 2D */ 3014 if ((p >= vStart) && (p < vEnd)) { 3015 /* Old vertices stay the same */ 3016 ++numLeavesNew; 3017 } else if ((p >= fStart) && (p < fEnd)) { 3018 /* Old faces add new faces and vertex */ 3019 numLeavesNew += 2 + 1; 3020 } else if ((p >= cStart) && (p < cEnd)) { 3021 /* Old cells add new cells, interior faces, and vertex */ 3022 numLeavesNew += 4 + 4 + 1; 3023 } 3024 break; 3025 case 5: 3026 /* Simplicial 3D */ 3027 if ((p >= vStart) && (p < vEnd)) { 3028 /* Old vertices stay the same */ 3029 ++numLeavesNew; 3030 } else if ((p >= eStart) && (p < eEnd)) { 3031 /* Old edges add new edges and vertex */ 3032 numLeavesNew += 2 + 1; 3033 } else if ((p >= fStart) && (p < fEnd)) { 3034 /* Old faces add new faces and face edges */ 3035 numLeavesNew += 4 + 3; 3036 } else if ((p >= cStart) && (p < cEnd)) { 3037 /* Old cells add new cells and interior faces and edges */ 3038 numLeavesNew += 8 + 8 + 1; 3039 } 3040 break; 3041 case 6: 3042 /* Hex 3D */ 3043 if ((p >= vStart) && (p < vEnd)) { 3044 /* Old vertices stay the same */ 3045 ++numLeavesNew; 3046 } else if ((p >= eStart) && (p < eEnd)) { 3047 /* Old edges add new edges, and vertex */ 3048 numLeavesNew += 2 + 1; 3049 } else if ((p >= fStart) && (p < fEnd)) { 3050 /* Old faces add new faces, edges, and vertex */ 3051 numLeavesNew += 4 + 4 + 1; 3052 } else if ((p >= cStart) && (p < cEnd)) { 3053 /* Old cells add new cells, faces, edges, and vertex */ 3054 numLeavesNew += 8 + 12 + 6 + 1; 3055 } 3056 break; 3057 default: 3058 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 3059 } 3060 } 3061 /* Communicate depthSizes for each remote rank */ 3062 ierr = DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);CHKERRQ(ierr); 3063 ierr = ISGetLocalSize(processRanks, &numNeighbors);CHKERRQ(ierr); 3064 ierr = PetscMalloc5((depth+1)*numNeighbors,PetscInt,&rdepthSize,numNeighbors,PetscInt,&rvStartNew,numNeighbors,PetscInt,&reStartNew,numNeighbors,PetscInt,&rfStartNew,numNeighbors,PetscInt,&rcStartNew);CHKERRQ(ierr); 3065 ierr = PetscMalloc7(depth+1,PetscInt,&depthSizeOld,(depth+1)*numNeighbors,PetscInt,&rdepthSizeOld,(depth+1)*numNeighbors,PetscInt,&rdepthMaxOld,numNeighbors,PetscInt,&rvStart,numNeighbors,PetscInt,&reStart,numNeighbors,PetscInt,&rfStart,numNeighbors,PetscInt,&rcStart);CHKERRQ(ierr); 3066 ierr = MPI_Type_contiguous(depth+1, MPIU_INT, &depthType);CHKERRQ(ierr); 3067 ierr = MPI_Type_commit(&depthType);CHKERRQ(ierr); 3068 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 3069 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 3070 for (n = 0; n < numNeighbors; ++n) { 3071 ierr = GetDepthStart_Private(depth, &rdepthSize[n*(depth+1)], &rcStartNew[n], &rfStartNew[n], &reStartNew[n], &rvStartNew[n]);CHKERRQ(ierr); 3072 } 3073 depthSizeOld[depth] = cMax; 3074 depthSizeOld[0] = vMax; 3075 depthSizeOld[depth-1] = fMax; 3076 depthSizeOld[1] = eMax; 3077 3078 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 3079 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 3080 3081 depthSizeOld[depth] = cEnd - cStart; 3082 depthSizeOld[0] = vEnd - vStart; 3083 depthSizeOld[depth-1] = fEnd - fStart; 3084 depthSizeOld[1] = eEnd - eStart; 3085 3086 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 3087 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 3088 for (n = 0; n < numNeighbors; ++n) { 3089 ierr = GetDepthStart_Private(depth, &rdepthSizeOld[n*(depth+1)], &rcStart[n], &rfStart[n], &reStart[n], &rvStart[n]);CHKERRQ(ierr); 3090 } 3091 ierr = MPI_Type_free(&depthType);CHKERRQ(ierr); 3092 ierr = PetscSFDestroy(&sfProcess);CHKERRQ(ierr); 3093 /* Calculate new point SF */ 3094 ierr = PetscMalloc(numLeavesNew * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 3095 ierr = PetscMalloc(numLeavesNew * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 3096 ierr = ISGetIndices(processRanks, &neighbors);CHKERRQ(ierr); 3097 for (l = 0, m = 0; l < numLeaves; ++l) { 3098 PetscInt p = localPoints[l]; 3099 PetscInt rp = remotePoints[l].index, n; 3100 PetscMPIInt rrank = remotePoints[l].rank; 3101 3102 ierr = PetscFindInt(rrank, numNeighbors, neighbors, &n);CHKERRQ(ierr); 3103 if (n < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %d", rrank); 3104 switch (refiner) { 3105 case 1: 3106 /* Simplicial 2D */ 3107 if ((p >= vStart) && (p < vEnd)) { 3108 /* Old vertices stay the same */ 3109 localPointsNew[m] = vStartNew + (p - vStart); 3110 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 3111 remotePointsNew[m].rank = rrank; 3112 ++m; 3113 } else if ((p >= fStart) && (p < fEnd)) { 3114 /* Old faces add new faces and vertex */ 3115 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 3116 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 3117 remotePointsNew[m].rank = rrank; 3118 ++m; 3119 for (r = 0; r < 2; ++r, ++m) { 3120 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 3121 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 3122 remotePointsNew[m].rank = rrank; 3123 } 3124 } else if ((p >= cStart) && (p < cEnd)) { 3125 /* Old cells add new cells and interior faces */ 3126 for (r = 0; r < 4; ++r, ++m) { 3127 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 3128 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 3129 remotePointsNew[m].rank = rrank; 3130 } 3131 for (r = 0; r < 3; ++r, ++m) { 3132 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 3133 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*3 + r; 3134 remotePointsNew[m].rank = rrank; 3135 } 3136 } 3137 break; 3138 case 2: 3139 /* Hex 2D */ 3140 if ((p >= vStart) && (p < vEnd)) { 3141 /* Old vertices stay the same */ 3142 localPointsNew[m] = vStartNew + (p - vStart); 3143 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 3144 remotePointsNew[m].rank = rrank; 3145 ++m; 3146 } else if ((p >= fStart) && (p < fEnd)) { 3147 /* Old faces add new faces and vertex */ 3148 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 3149 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 3150 remotePointsNew[m].rank = rrank; 3151 ++m; 3152 for (r = 0; r < 2; ++r, ++m) { 3153 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 3154 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 3155 remotePointsNew[m].rank = rrank; 3156 } 3157 } else if ((p >= cStart) && (p < cEnd)) { 3158 /* Old cells add new cells, interior faces, and vertex */ 3159 for (r = 0; r < 4; ++r, ++m) { 3160 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 3161 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 3162 remotePointsNew[m].rank = rrank; 3163 } 3164 for (r = 0; r < 4; ++r, ++m) { 3165 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 3166 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*4 + r; 3167 remotePointsNew[m].rank = rrank; 3168 } 3169 for (r = 0; r < 1; ++r, ++m) { 3170 localPointsNew[m] = vStartNew + (fEnd - fStart) + (p - cStart) + r; 3171 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1] + (rp - rcStart[n]) + r; 3172 remotePointsNew[m].rank = rrank; 3173 } 3174 } 3175 break; 3176 case 3: 3177 /* Hybrid simplicial 2D */ 3178 if ((p >= vStart) && (p < vEnd)) { 3179 /* Old vertices stay the same */ 3180 localPointsNew[m] = vStartNew + (p - vStart); 3181 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 3182 remotePointsNew[m].rank = rrank; 3183 ++m; 3184 } else if ((p >= fStart) && (p < fMax)) { 3185 /* Old interior faces add new faces and vertex */ 3186 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 3187 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 3188 remotePointsNew[m].rank = rrank; 3189 ++m; 3190 for (r = 0; r < 2; ++r, ++m) { 3191 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 3192 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 3193 remotePointsNew[m].rank = rrank; 3194 } 3195 } else if ((p >= fMax) && (p < fEnd)) { 3196 /* Old hybrid faces stay the same */ 3197 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - fMax); 3198 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rdepthMaxOld[n*(depth+1)+depth-1]); 3199 remotePointsNew[m].rank = rrank; 3200 ++m; 3201 } else if ((p >= cStart) && (p < cMax)) { 3202 /* Old interior cells add new cells and interior faces */ 3203 for (r = 0; r < 4; ++r, ++m) { 3204 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 3205 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 3206 remotePointsNew[m].rank = rrank; 3207 } 3208 for (r = 0; r < 3; ++r, ++m) { 3209 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - cStart)*3 + r; 3210 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rcStart[n])*3 + r; 3211 remotePointsNew[m].rank = rrank; 3212 } 3213 } else if ((p >= cStart) && (p < cMax)) { 3214 /* Old hybrid cells add new cells and hybrid face */ 3215 for (r = 0; r < 2; ++r, ++m) { 3216 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 3217 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 3218 remotePointsNew[m].rank = rrank; 3219 } 3220 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 3221 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rdepthMaxOld[n*(depth+1)+depth] - rcStart[n])*3 + (rp - rdepthMaxOld[n*(depth+1)+depth]); 3222 remotePointsNew[m].rank = rrank; 3223 ++m; 3224 } 3225 break; 3226 case 5: 3227 /* Simplicial 3D */ 3228 if ((p >= vStart) && (p < vEnd)) { 3229 /* Old vertices stay the same */ 3230 localPointsNew[m] = vStartNew + (p - vStart); 3231 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 3232 remotePointsNew[m].rank = rrank; 3233 ++m; 3234 } else if ((p >= eStart) && (p < eEnd)) { 3235 /* Old edges add new edges and vertex */ 3236 for (r = 0; r < 2; ++r, ++m) { 3237 localPointsNew[m] = eStartNew + (p - eStart)*2 + r; 3238 remotePointsNew[m].index = reStartNew[n] + (rp - reStart[n])*2 + r; 3239 remotePointsNew[m].rank = rrank; 3240 } 3241 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - eStart); 3242 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - reStart[n]); 3243 remotePointsNew[m].rank = rrank; 3244 ++m; 3245 } else if ((p >= fStart) && (p < fEnd)) { 3246 /* Old faces add new faces and face edges */ 3247 for (r = 0; r < 4; ++r, ++m) { 3248 localPointsNew[m] = fStartNew + (p - fStart)*4 + r; 3249 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*4 + r; 3250 remotePointsNew[m].rank = rrank; 3251 } 3252 for (r = 0; r < 3; ++r, ++m) { 3253 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (p - fStart)*3 + r; 3254 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + (rp - rfStart[n])*3 + r; 3255 remotePointsNew[m].rank = rrank; 3256 } 3257 } else if ((p >= cStart) && (p < cEnd)) { 3258 /* Old cells add new cells and interior faces and edges */ 3259 for (r = 0; r < 8; ++r, ++m) { 3260 localPointsNew[m] = cStartNew + (p - cStart)*8 + r; 3261 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*8 + r; 3262 remotePointsNew[m].rank = rrank; 3263 } 3264 for (r = 0; r < 8; ++r, ++m) { 3265 localPointsNew[m] = fStartNew + (fEnd - fStart)*4 + (p - cStart)*8 + r; 3266 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*8 + r; 3267 remotePointsNew[m].rank = rrank; 3268 } 3269 for (r = 0; r < 1; ++r, ++m) { 3270 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (p - cStart)*1 + r; 3271 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + rdepthSizeOld[n*(depth+1)+depth-1]*3 + (rp - rcStart[n])*1 + r; 3272 remotePointsNew[m].rank = rrank; 3273 } 3274 } 3275 break; 3276 case 6: 3277 /* Hex 3D */ 3278 if ((p >= vStart) && (p < vEnd)) { 3279 /* Old vertices stay the same */ 3280 localPointsNew[m] = vStartNew + (p - vStart); 3281 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 3282 remotePointsNew[m].rank = rrank; 3283 ++m; 3284 } else if ((p >= eStart) && (p < eEnd)) { 3285 /* Old edges add new edges and vertex */ 3286 for (r = 0; r < 2; ++r, ++m) { 3287 localPointsNew[m] = eStartNew + (p - eStart)*2 + r; 3288 remotePointsNew[m].index = reStartNew[n] + (rp - reStart[n])*2 + r; 3289 remotePointsNew[m].rank = rrank; 3290 } 3291 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - eStart); 3292 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - reStart[n]); 3293 remotePointsNew[m].rank = rrank; 3294 ++m; 3295 } else if ((p >= fStart) && (p < fEnd)) { 3296 /* Old faces add new faces, edges, and vertex */ 3297 for (r = 0; r < 4; ++r, ++m) { 3298 localPointsNew[m] = fStartNew + (p - fStart)*4 + r; 3299 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*4 + r; 3300 remotePointsNew[m].rank = rrank; 3301 } 3302 for (r = 0; r < 4; ++r, ++m) { 3303 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (p - fStart)*4 + r; 3304 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + (rp - rfStart[n])*4 + r; 3305 remotePointsNew[m].rank = rrank; 3306 } 3307 localPointsNew[m] = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (p - fStart); 3308 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + rdepthSizeOld[n*(depth+1)+1] + (rp - rfStart[n]); 3309 remotePointsNew[m].rank = rrank; 3310 ++m; 3311 } else if ((p >= cStart) && (p < cEnd)) { 3312 /* Old cells add new cells, faces, edges, and vertex */ 3313 for (r = 0; r < 8; ++r, ++m) { 3314 localPointsNew[m] = cStartNew + (p - cStart)*8 + r; 3315 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*8 + r; 3316 remotePointsNew[m].rank = rrank; 3317 } 3318 for (r = 0; r < 12; ++r, ++m) { 3319 localPointsNew[m] = fStartNew + (fEnd - fStart)*4 + (p - cStart)*12 + r; 3320 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*12 + r; 3321 remotePointsNew[m].rank = rrank; 3322 } 3323 for (r = 0; r < 6; ++r, ++m) { 3324 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (p - cStart)*6 + r; 3325 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*6 + r; 3326 remotePointsNew[m].rank = rrank; 3327 } 3328 for (r = 0; r < 1; ++r, ++m) { 3329 localPointsNew[m] = vStartNew + (eEnd - eStart) + (fEnd - fStart) + (p - cStart) + r; 3330 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+1] + rdepthSizeOld[n*(depth+1)+depth-1] + (rp - rcStart[n]) + r; 3331 remotePointsNew[m].rank = rrank; 3332 } 3333 } 3334 break; 3335 default: 3336 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 3337 } 3338 } 3339 ierr = ISRestoreIndices(processRanks, &neighbors);CHKERRQ(ierr); 3340 ierr = ISDestroy(&processRanks);CHKERRQ(ierr); 3341 ierr = PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 3342 ierr = PetscFree5(rdepthSize,rvStartNew,reStartNew,rfStartNew,rcStartNew);CHKERRQ(ierr); 3343 ierr = PetscFree7(depthSizeOld,rdepthSizeOld,rdepthMaxOld,rvStart,reStart,rfStart,rcStart);CHKERRQ(ierr); 3344 PetscFunctionReturn(0); 3345 } 3346 3347 #undef __FUNCT__ 3348 #define __FUNCT__ "CellRefinerCreateLabels" 3349 PetscErrorCode CellRefinerCreateLabels(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 3350 { 3351 PetscInt numLabels, l; 3352 PetscInt depth, newp, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eStartNew, eEnd, eMax, r; 3353 PetscErrorCode ierr; 3354 3355 PetscFunctionBegin; 3356 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3357 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 3358 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 3359 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 3360 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3361 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 3362 ierr = DMPlexGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 3363 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 3364 switch (refiner) { 3365 case 3: 3366 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 3367 cMax = PetscMin(cEnd, cMax); 3368 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 3369 fMax = PetscMin(fEnd, fMax); 3370 } 3371 for (l = 0; l < numLabels; ++l) { 3372 DMLabel label, labelNew; 3373 const char *lname; 3374 PetscBool isDepth; 3375 IS valueIS; 3376 const PetscInt *values; 3377 PetscInt numValues, val; 3378 3379 ierr = DMPlexGetLabelName(dm, l, &lname);CHKERRQ(ierr); 3380 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 3381 if (isDepth) continue; 3382 ierr = DMPlexCreateLabel(rdm, lname);CHKERRQ(ierr); 3383 ierr = DMPlexGetLabel(dm, lname, &label);CHKERRQ(ierr); 3384 ierr = DMPlexGetLabel(rdm, lname, &labelNew);CHKERRQ(ierr); 3385 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 3386 ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr); 3387 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 3388 for (val = 0; val < numValues; ++val) { 3389 IS pointIS; 3390 const PetscInt *points; 3391 PetscInt numPoints, n; 3392 3393 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 3394 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 3395 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 3396 for (n = 0; n < numPoints; ++n) { 3397 const PetscInt p = points[n]; 3398 switch (refiner) { 3399 case 1: 3400 /* Simplicial 2D */ 3401 if ((p >= vStart) && (p < vEnd)) { 3402 /* Old vertices stay the same */ 3403 newp = vStartNew + (p - vStart); 3404 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3405 } else if ((p >= fStart) && (p < fEnd)) { 3406 /* Old faces add new faces and vertex */ 3407 newp = vStartNew + (vEnd - vStart) + (p - fStart); 3408 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3409 for (r = 0; r < 2; ++r) { 3410 newp = fStartNew + (p - fStart)*2 + r; 3411 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3412 } 3413 } else if ((p >= cStart) && (p < cEnd)) { 3414 /* Old cells add new cells and interior faces */ 3415 for (r = 0; r < 4; ++r) { 3416 newp = cStartNew + (p - cStart)*4 + r; 3417 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3418 } 3419 for (r = 0; r < 3; ++r) { 3420 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 3421 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3422 } 3423 } 3424 break; 3425 case 2: 3426 /* Hex 2D */ 3427 if ((p >= vStart) && (p < vEnd)) { 3428 /* Old vertices stay the same */ 3429 newp = vStartNew + (p - vStart); 3430 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3431 } else if ((p >= fStart) && (p < fEnd)) { 3432 /* Old faces add new faces and vertex */ 3433 newp = vStartNew + (vEnd - vStart) + (p - fStart); 3434 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3435 for (r = 0; r < 2; ++r) { 3436 newp = fStartNew + (p - fStart)*2 + r; 3437 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3438 } 3439 } else if ((p >= cStart) && (p < cEnd)) { 3440 /* Old cells add new cells and interior faces and vertex */ 3441 for (r = 0; r < 4; ++r) { 3442 newp = cStartNew + (p - cStart)*4 + r; 3443 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3444 } 3445 for (r = 0; r < 4; ++r) { 3446 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 3447 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3448 } 3449 newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (p - cStart); 3450 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3451 } 3452 break; 3453 case 3: 3454 /* Hybrid simplicial 2D */ 3455 if ((p >= vStart) && (p < vEnd)) { 3456 /* Old vertices stay the same */ 3457 newp = vStartNew + (p - vStart); 3458 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3459 } else if ((p >= fStart) && (p < fMax)) { 3460 /* Old interior faces add new faces and vertex */ 3461 newp = vStartNew + (vEnd - vStart) + (p - fStart); 3462 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3463 for (r = 0; r < 2; ++r) { 3464 newp = fStartNew + (p - fStart)*2 + r; 3465 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3466 } 3467 } else if ((p >= fMax) && (p < fEnd)) { 3468 /* Old hybrid faces stay the same */ 3469 newp = fStartNew + (fMax - fStart)*2 + (p - fMax); 3470 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3471 } else if ((p >= cStart) && (p < cMax)) { 3472 /* Old interior cells add new cells and interior faces */ 3473 for (r = 0; r < 4; ++r) { 3474 newp = cStartNew + (p - cStart)*4 + r; 3475 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3476 } 3477 for (r = 0; r < 3; ++r) { 3478 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 3479 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3480 } 3481 } else if ((p >= cMax) && (p < cEnd)) { 3482 /* Old hybrid cells add new cells and hybrid face */ 3483 for (r = 0; r < 2; ++r) { 3484 newp = cStartNew + (cMax - cStart)*4 + (p - cMax)*2 + r; 3485 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3486 } 3487 newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 3488 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3489 } 3490 break; 3491 case 5: 3492 /* Simplicial 3D */ 3493 if ((p >= vStart) && (p < vEnd)) { 3494 /* Old vertices stay the same */ 3495 newp = vStartNew + (p - vStart); 3496 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3497 } else if ((p >= eStart) && (p < eEnd)) { 3498 /* Old edges add new edges and vertex */ 3499 for (r = 0; r < 2; ++r) { 3500 newp = eStartNew + (p - eStart)*2 + r; 3501 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3502 } 3503 newp = vStartNew + (vEnd - vStart) + (p - eStart); 3504 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3505 } else if ((p >= fStart) && (p < fEnd)) { 3506 /* Old faces add new faces and edges */ 3507 for (r = 0; r < 4; ++r) { 3508 newp = fStartNew + (p - fStart)*4 + r; 3509 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3510 } 3511 for (r = 0; r < 3; ++r) { 3512 newp = eStartNew + (eEnd - eStart)*2 + (p - fStart)*3 + r; 3513 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3514 } 3515 } else if ((p >= cStart) && (p < cEnd)) { 3516 /* Old cells add new cells and interior faces and edges */ 3517 for (r = 0; r < 8; ++r) { 3518 newp = cStartNew + (p - cStart)*8 + r; 3519 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3520 } 3521 for (r = 0; r < 8; ++r) { 3522 newp = fStartNew + (fEnd - fStart)*4 + (p - cStart)*8 + r; 3523 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3524 } 3525 for (r = 0; r < 1; ++r) { 3526 newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (p - cStart)*1 + r; 3527 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3528 } 3529 } 3530 break; 3531 case 6: 3532 /* Hex 3D */ 3533 if ((p >= vStart) && (p < vEnd)) { 3534 /* Old vertices stay the same */ 3535 newp = vStartNew + (p - vStart); 3536 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3537 } else if ((p >= fStart) && (p < fEnd)) { 3538 /* Old edges add new edges and vertex */ 3539 for (r = 0; r < 2; ++r) { 3540 newp = eStartNew + (p - eStart)*2 + r; 3541 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3542 } 3543 newp = vStartNew + (vEnd - vStart) + (p - eStart); 3544 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3545 } else if ((p >= fStart) && (p < fEnd)) { 3546 /* Old faces add new faces, edges, and vertex */ 3547 for (r = 0; r < 4; ++r) { 3548 newp = fStartNew + (p - fStart)*4 + r; 3549 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3550 } 3551 for (r = 0; r < 4; ++r) { 3552 newp = eStartNew + (eEnd - eStart)*2 + (p - fStart)*4 + r; 3553 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3554 } 3555 newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (p - fStart); 3556 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3557 } else if ((p >= cStart) && (p < cEnd)) { 3558 /* Old cells add new cells, faces, edges, and vertex */ 3559 for (r = 0; r < 8; ++r) { 3560 newp = cStartNew + (p - cStart)*8 + r; 3561 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3562 } 3563 for (r = 0; r < 12; ++r) { 3564 newp = fStartNew + (fEnd - fStart)*4 + (p - cStart)*12 + r; 3565 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3566 } 3567 for (r = 0; r < 6; ++r) { 3568 newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (p - cStart)*6 + r; 3569 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3570 } 3571 newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (p - cStart); 3572 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3573 } 3574 break; 3575 default: 3576 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 3577 } 3578 } 3579 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 3580 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 3581 } 3582 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 3583 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 3584 if (0) { 3585 ierr = PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);CHKERRQ(ierr); 3586 ierr = DMLabelView(labelNew, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 3587 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 3588 } 3589 } 3590 PetscFunctionReturn(0); 3591 } 3592 3593 #undef __FUNCT__ 3594 #define __FUNCT__ "DMPlexRefineUniform_Internal" 3595 /* This will only work for interpolated meshes */ 3596 PetscErrorCode DMPlexRefineUniform_Internal(DM dm, CellRefiner cellRefiner, DM *dmRefined) 3597 { 3598 DM rdm; 3599 PetscInt *depthSize; 3600 PetscInt dim, depth = 0, d, pStart = 0, pEnd = 0; 3601 PetscErrorCode ierr; 3602 3603 PetscFunctionBegin; 3604 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &rdm);CHKERRQ(ierr); 3605 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3606 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 3607 ierr = DMPlexSetDimension(rdm, dim);CHKERRQ(ierr); 3608 /* Calculate number of new points of each depth */ 3609 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3610 ierr = PetscMalloc((depth+1) * sizeof(PetscInt), &depthSize);CHKERRQ(ierr); 3611 ierr = PetscMemzero(depthSize, (depth+1) * sizeof(PetscInt));CHKERRQ(ierr); 3612 ierr = CellRefinerGetSizes(cellRefiner, dm, depthSize);CHKERRQ(ierr); 3613 /* Step 1: Set chart */ 3614 for (d = 0; d <= depth; ++d) pEnd += depthSize[d]; 3615 ierr = DMPlexSetChart(rdm, pStart, pEnd);CHKERRQ(ierr); 3616 /* Step 2: Set cone/support sizes */ 3617 ierr = CellRefinerSetConeSizes(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3618 /* Step 3: Setup refined DM */ 3619 ierr = DMSetUp(rdm);CHKERRQ(ierr); 3620 /* Step 4: Set cones and supports */ 3621 ierr = CellRefinerSetCones(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3622 /* Step 5: Stratify */ 3623 ierr = DMPlexStratify(rdm);CHKERRQ(ierr); 3624 /* Step 6: Set coordinates for vertices */ 3625 ierr = CellRefinerSetCoordinates(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3626 /* Step 7: Create pointSF */ 3627 ierr = CellRefinerCreateSF(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3628 /* Step 8: Create labels */ 3629 ierr = CellRefinerCreateLabels(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3630 ierr = PetscFree(depthSize);CHKERRQ(ierr); 3631 3632 *dmRefined = rdm; 3633 PetscFunctionReturn(0); 3634 } 3635 3636 #undef __FUNCT__ 3637 #define __FUNCT__ "DMPlexSetRefinementUniform" 3638 PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform) 3639 { 3640 DM_Plex *mesh = (DM_Plex*) dm->data; 3641 3642 PetscFunctionBegin; 3643 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3644 mesh->refinementUniform = refinementUniform; 3645 PetscFunctionReturn(0); 3646 } 3647 3648 #undef __FUNCT__ 3649 #define __FUNCT__ "DMPlexGetRefinementUniform" 3650 PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform) 3651 { 3652 DM_Plex *mesh = (DM_Plex*) dm->data; 3653 3654 PetscFunctionBegin; 3655 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3656 PetscValidPointer(refinementUniform, 2); 3657 *refinementUniform = mesh->refinementUniform; 3658 PetscFunctionReturn(0); 3659 } 3660 3661 #undef __FUNCT__ 3662 #define __FUNCT__ "DMPlexSetRefinementLimit" 3663 PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit) 3664 { 3665 DM_Plex *mesh = (DM_Plex*) dm->data; 3666 3667 PetscFunctionBegin; 3668 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3669 mesh->refinementLimit = refinementLimit; 3670 PetscFunctionReturn(0); 3671 } 3672 3673 #undef __FUNCT__ 3674 #define __FUNCT__ "DMPlexGetRefinementLimit" 3675 PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit) 3676 { 3677 DM_Plex *mesh = (DM_Plex*) dm->data; 3678 3679 PetscFunctionBegin; 3680 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3681 PetscValidPointer(refinementLimit, 2); 3682 /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */ 3683 *refinementLimit = mesh->refinementLimit; 3684 PetscFunctionReturn(0); 3685 } 3686 3687 #undef __FUNCT__ 3688 #define __FUNCT__ "DMPlexGetCellRefiner_Internal" 3689 PetscErrorCode DMPlexGetCellRefiner_Internal(DM dm, CellRefiner *cellRefiner) 3690 { 3691 PetscInt dim, cStart, coneSize, cMax; 3692 PetscErrorCode ierr; 3693 3694 PetscFunctionBegin; 3695 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 3696 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 3697 ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr); 3698 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 3699 switch (dim) { 3700 case 2: 3701 switch (coneSize) { 3702 case 3: 3703 if (cMax >= 0) *cellRefiner = 3; /* Hybrid */ 3704 else *cellRefiner = 1; /* Triangular */ 3705 break; 3706 case 4: 3707 if (cMax >= 0) *cellRefiner = 4; /* Hybrid */ 3708 else *cellRefiner = 2; /* Quadrilateral */ 3709 break; 3710 default: 3711 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown coneSize %d in dimension %d for cell refiner", coneSize, dim); 3712 } 3713 break; 3714 case 3: 3715 switch (coneSize) { 3716 case 4: 3717 if (cMax >= 0) *cellRefiner = 7; /* Hybrid */ 3718 else *cellRefiner = 5; /* Tetrahedral */ 3719 break; 3720 case 6: 3721 if (cMax >= 0) *cellRefiner = 8; /* Hybrid */ 3722 else *cellRefiner = 6; /* hexahedral */ 3723 break; 3724 default: 3725 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown coneSize %d in dimension %d for cell refiner", coneSize, dim); 3726 } 3727 break; 3728 default: 3729 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown dimension %d for cell refiner", dim); 3730 } 3731 PetscFunctionReturn(0); 3732 } 3733