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