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+0, coneNew);CHKERRQ(ierr); 2102 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 2103 #if 1 2104 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); 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+1, coneNew);CHKERRQ(ierr); 2123 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 2124 #if 1 2125 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); 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+2, coneNew);CHKERRQ(ierr); 2144 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 2145 #if 1 2146 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); 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+3, coneNew);CHKERRQ(ierr); 2165 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 2166 #if 1 2167 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); 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 /* 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], coneSize, c, supportSize, s; 2183 2184 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 2185 /* TODO: Redo using orientation information */ 2186 coneNew[0] = eStartNew + (cone[(r+3)%4] - eStart)*2 + 1; 2187 coneNew[1] = eStartNew + (cone[r] - eStart)*2 + 0; 2188 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + r; 2189 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + (r+3)%4; 2190 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2191 #if 1 2192 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2193 for (p = 0; p < 4; ++p) { 2194 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); 2195 } 2196 #endif 2197 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 2198 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 2199 for (s = 0; s < supportSize; ++s) { 2200 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 2201 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2202 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 2203 for (c = 0; c < coneSize; ++c) { 2204 if (cone[c] == f) break; 2205 } 2206 /* TODO: Redo using orientation information */ 2207 supportRef[s] = cStartNew + (support[s] - cStart)*8 + newCells[c*4+r]; 2208 } 2209 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2210 #if 1 2211 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2212 for (p = 0; p < supportSize; ++p) { 2213 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); 2214 } 2215 #endif 2216 } 2217 } 2218 /* Interior faces have 4 edges and 2 cells */ 2219 for (c = cStart; c < cEnd; ++c) { 2220 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}; 2221 const PetscInt *cone; 2222 PetscInt coneNew[4], supportNew[2]; 2223 2224 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2225 for (r = 0; r < 12; ++r) { 2226 const PetscInt newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + r; 2227 2228 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[r] - fStart); 2229 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart) + (c - cStart); 2230 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[r] - fStart); 2231 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (cone[r] - fStart); 2232 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2233 #if 1 2234 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2235 for (p = 0; p < 4; ++p) { 2236 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); 2237 } 2238 #endif 2239 supportNew[0] = cStartNew + (c - cStart)*8 + newCells[r*2+0]; 2240 supportNew[1] = cStartNew + (c - cStart)*8 + newCells[r*2+1]; 2241 ierr = DMPlexSetSupport(rdm, newp, supportNew);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 < 2; ++p) { 2245 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); 2246 } 2247 #endif 2248 } 2249 } 2250 /* Split edges have 2 vertices and the same faces as the parent */ 2251 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 2252 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 2253 for (e = eStart; e < eEnd; ++e) { 2254 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 2255 2256 for (r = 0; r < 2; ++r) { 2257 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 2258 const PetscInt *cone, *ornt, *support; 2259 PetscInt coneNew[2], coneSize, c, supportSize, s; 2260 2261 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 2262 coneNew[0] = vStartNew + (cone[0] - vStart); 2263 coneNew[1] = vStartNew + (cone[1] - vStart); 2264 coneNew[(r+1)%2] = newv; 2265 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2266 #if 1 2267 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2268 for (p = 0; p < 2; ++p) { 2269 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); 2270 } 2271 #endif 2272 ierr = DMPlexGetSupportSize(dm, e, &supportSize);CHKERRQ(ierr); 2273 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 2274 for (s = 0; s < supportSize; ++s) { 2275 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 2276 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2277 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 2278 for (c = 0; c < coneSize; ++c) { 2279 if (cone[c] == e) break; 2280 } 2281 supportRef[s] = fStartNew + (support[s] - fStart)*4 + (ornt[c] < 0 ? (c+1-r)%4 : (c+r)%4); 2282 } 2283 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2284 #if 1 2285 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2286 for (p = 0; p < supportSize; ++p) { 2287 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); 2288 } 2289 #endif 2290 } 2291 } 2292 /* Face edges have 2 vertices and 2+cells faces */ 2293 for (f = fStart; f < fEnd; ++f) { 2294 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}; 2295 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 2296 const PetscInt *cone, *coneCell, *support; 2297 PetscInt coneNew[2], coneSize, c, supportSize, s; 2298 2299 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 2300 for (r = 0; r < 4; ++r) { 2301 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + r; 2302 2303 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - eStart); 2304 coneNew[1] = newv; 2305 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2306 #if 1 2307 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2308 for (p = 0; p < 2; ++p) { 2309 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); 2310 } 2311 #endif 2312 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 2313 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 2314 supportRef[0] = fStartNew + (f - fStart)*4 + r; 2315 supportRef[1] = fStartNew + (f - fStart)*4 + (r+1)%4; 2316 for (s = 0; s < supportSize; ++s) { 2317 ierr = DMPlexGetConeSize(dm, f, &coneSize);CHKERRQ(ierr); 2318 ierr = DMPlexGetCone(dm, f, &coneCell);CHKERRQ(ierr); 2319 for (c = 0; c < coneSize; ++c) if (coneCell[c] == f) break; 2320 supportRef[2+s] = fStartNew + (f - fStart)*4 + newFaces[c*4+r]; 2321 } 2322 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2323 #if 1 2324 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2325 for (p = 0; p < 2+supportSize; ++p) { 2326 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); 2327 } 2328 #endif 2329 } 2330 } 2331 /* Cell edges have 2 vertices and 4 faces */ 2332 for (c = cStart; c < cEnd; ++c) { 2333 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}; 2334 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (c - cStart); 2335 const PetscInt *cone; 2336 PetscInt coneNew[2], supportNew[4]; 2337 2338 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2339 for (r = 0; r < 6; ++r) { 2340 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + r; 2341 2342 coneNew[0] = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (cone[r] - fStart); 2343 coneNew[1] = newv; 2344 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2345 #if 1 2346 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2347 for (p = 0; p < 2; ++p) { 2348 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); 2349 } 2350 #endif 2351 for (f = 0; f < 4; ++f) supportNew[f] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + newFaces[r*4+f]; 2352 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 2353 #if 1 2354 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2355 for (p = 0; p < 4; ++p) { 2356 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); 2357 } 2358 #endif 2359 } 2360 } 2361 /* Old vertices have identical supports */ 2362 for (v = vStart; v < vEnd; ++v) { 2363 const PetscInt newp = vStartNew + (v - vStart); 2364 const PetscInt *support, *cone; 2365 PetscInt size, s; 2366 2367 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 2368 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 2369 for (s = 0; s < size; ++s) { 2370 PetscInt r = 0; 2371 2372 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2373 if (cone[1] == v) r = 1; 2374 supportRef[s] = eStartNew + (support[s] - eStart)*2 + r; 2375 } 2376 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2377 #if 1 2378 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2379 for (p = 0; p < size; ++p) { 2380 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); 2381 } 2382 #endif 2383 } 2384 /* Edge vertices have 2 + faces supports */ 2385 for (e = eStart; e < eEnd; ++e) { 2386 const PetscInt newp = vStartNew + (vEnd - vStart) + (e - eStart); 2387 const PetscInt *cone, *support; 2388 PetscInt size, s; 2389 2390 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 2391 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 2392 supportRef[0] = eStartNew + (e - eStart)*2 + 0; 2393 supportRef[1] = eStartNew + (e - eStart)*2 + 1; 2394 for (s = 0; s < size; ++s) { 2395 PetscInt r; 2396 2397 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2398 for (r = 0; r < 4; ++r) if (cone[r] == e) break; 2399 supportRef[2+s] = eStartNew + (eEnd - eStart)*2 + (support[s] - fStart)*4 + r; 2400 } 2401 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2402 #if 1 2403 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2404 for (p = 0; p < 2+size; ++p) { 2405 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); 2406 } 2407 #endif 2408 } 2409 /* Face vertices have 4 + cells supports */ 2410 for (f = fStart; f < fEnd; ++f) { 2411 const PetscInt newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 2412 const PetscInt *cone, *support; 2413 PetscInt size, s; 2414 2415 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 2416 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 2417 for (r = 0; r < 4; ++r) supportRef[r] = eStartNew + (e - eStart)*2 + (f - fStart)*4 + r; 2418 for (s = 0; s < size; ++s) { 2419 PetscInt r; 2420 2421 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2422 for (r = 0; r < 6; ++r) if (cone[r] == f) break; 2423 supportRef[4+s] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (support[s] - cStart)*6 + r; 2424 } 2425 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2426 #if 1 2427 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2428 for (p = 0; p < 4+size; ++p) { 2429 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); 2430 } 2431 #endif 2432 } 2433 /* Cell vertices have 6 supports */ 2434 for (c = cStart; c < cEnd; ++c) { 2435 const PetscInt newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (c - cStart); 2436 PetscInt supportNew[6]; 2437 2438 for (r = 0; r < 6; ++r) { 2439 supportNew[r] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + r; 2440 } 2441 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 2442 } 2443 ierr = PetscFree(supportRef);CHKERRQ(ierr); 2444 break; 2445 default: 2446 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2447 } 2448 PetscFunctionReturn(0); 2449 } 2450 2451 #undef __FUNCT__ 2452 #define __FUNCT__ "CellRefinerSetCoordinates" 2453 PetscErrorCode CellRefinerSetCoordinates(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 2454 { 2455 PetscSection coordSection, coordSectionNew; 2456 Vec coordinates, coordinatesNew; 2457 PetscScalar *coords, *coordsNew; 2458 PetscInt dim, depth, coordSizeNew, cStart, cEnd, c, vStart, vStartNew, vEnd, v, eStart, eEnd, eMax, e, fStart, fEnd, fMax, f; 2459 PetscErrorCode ierr; 2460 2461 PetscFunctionBegin; 2462 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2463 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2464 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2465 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 2466 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2467 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2468 ierr = DMPlexGetHybridBounds(dm, NULL, &fMax, &eMax, NULL);CHKERRQ(ierr); 2469 ierr = GetDepthStart_Private(depth, depthSize, NULL, NULL, NULL, &vStartNew);CHKERRQ(ierr); 2470 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2471 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &coordSectionNew);CHKERRQ(ierr); 2472 ierr = PetscSectionSetNumFields(coordSectionNew, 1);CHKERRQ(ierr); 2473 ierr = PetscSectionSetFieldComponents(coordSectionNew, 0, dim);CHKERRQ(ierr); 2474 ierr = PetscSectionSetChart(coordSectionNew, vStartNew, vStartNew+depthSize[0]);CHKERRQ(ierr); 2475 if (fMax < 0) fMax = fEnd; 2476 if (eMax < 0) eMax = eEnd; 2477 /* All vertices have the dim coordinates */ 2478 for (v = vStartNew; v < vStartNew+depthSize[0]; ++v) { 2479 ierr = PetscSectionSetDof(coordSectionNew, v, dim);CHKERRQ(ierr); 2480 ierr = PetscSectionSetFieldDof(coordSectionNew, v, 0, dim);CHKERRQ(ierr); 2481 } 2482 ierr = PetscSectionSetUp(coordSectionNew);CHKERRQ(ierr); 2483 ierr = DMPlexSetCoordinateSection(rdm, coordSectionNew);CHKERRQ(ierr); 2484 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2485 ierr = PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);CHKERRQ(ierr); 2486 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinatesNew);CHKERRQ(ierr); 2487 ierr = PetscObjectSetName((PetscObject) coordinatesNew, "coordinates");CHKERRQ(ierr); 2488 ierr = VecSetSizes(coordinatesNew, coordSizeNew, PETSC_DETERMINE);CHKERRQ(ierr); 2489 ierr = VecSetFromOptions(coordinatesNew);CHKERRQ(ierr); 2490 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2491 ierr = VecGetArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 2492 switch (refiner) { 2493 case 6: /* Hex 3D */ 2494 /* Face vertices have the average of corner coordinates */ 2495 for (f = fStart; f < fEnd; ++f) { 2496 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 2497 PetscInt *cone = NULL; 2498 PetscInt closureSize, coneSize = 0, off[8], offnew, p, d; 2499 2500 ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2501 for (p = 0; p < closureSize*2; p += 2) { 2502 const PetscInt point = cone[p]; 2503 if ((point >= vStart) && (point < vEnd)) cone[coneSize++] = point; 2504 } 2505 for (v = 0; v < coneSize; ++v) { 2506 ierr = PetscSectionGetOffset(coordSection, cone[v], &off[v]);CHKERRQ(ierr); 2507 } 2508 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2509 for (d = 0; d < dim; ++d) { 2510 coordsNew[offnew+d] = 0.0; 2511 for (v = 0; v < coneSize; ++v) coordsNew[offnew+d] += coords[off[v]+d]; 2512 coordsNew[offnew+d] /= coneSize; 2513 } 2514 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2515 } 2516 case 2: /* Hex 2D */ 2517 /* Cell vertices have the average of corner coordinates */ 2518 for (c = cStart; c < cEnd; ++c) { 2519 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (c - cStart) + (dim > 2 ? (fEnd - fStart) : 0); 2520 PetscInt *cone = NULL; 2521 PetscInt closureSize, coneSize = 0, off[8], offnew, p, d; 2522 2523 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2524 for (p = 0; p < closureSize*2; p += 2) { 2525 const PetscInt point = cone[p]; 2526 if ((point >= vStart) && (point < vEnd)) cone[coneSize++] = point; 2527 } 2528 for (v = 0; v < coneSize; ++v) { 2529 ierr = PetscSectionGetOffset(coordSection, cone[v], &off[v]);CHKERRQ(ierr); 2530 } 2531 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2532 for (d = 0; d < dim; ++d) { 2533 coordsNew[offnew+d] = 0.0; 2534 for (v = 0; v < coneSize; ++v) coordsNew[offnew+d] += coords[off[v]+d]; 2535 coordsNew[offnew+d] /= coneSize; 2536 } 2537 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2538 } 2539 case 1: /* Simplicial 2D */ 2540 case 3: /* Hybrid Simplicial 2D */ 2541 case 5: /* Simplicial 3D */ 2542 /* Edge vertices have the average of endpoint coordinates */ 2543 for (e = eStart; e < eMax; ++e) { 2544 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 2545 const PetscInt *cone; 2546 PetscInt coneSize, offA, offB, offnew, d; 2547 2548 ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr); 2549 if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %d cone should have two vertices, not %d", e, coneSize); 2550 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 2551 ierr = PetscSectionGetOffset(coordSection, cone[0], &offA);CHKERRQ(ierr); 2552 ierr = PetscSectionGetOffset(coordSection, cone[1], &offB);CHKERRQ(ierr); 2553 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2554 for (d = 0; d < dim; ++d) { 2555 coordsNew[offnew+d] = 0.5*(coords[offA+d] + coords[offB+d]); 2556 } 2557 } 2558 /* Old vertices have the same coordinates */ 2559 for (v = vStart; v < vEnd; ++v) { 2560 const PetscInt newv = vStartNew + (v - vStart); 2561 PetscInt off, offnew, d; 2562 2563 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 2564 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2565 for (d = 0; d < dim; ++d) { 2566 coordsNew[offnew+d] = coords[off+d]; 2567 } 2568 } 2569 break; 2570 default: 2571 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2572 } 2573 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2574 ierr = VecRestoreArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 2575 ierr = DMSetCoordinatesLocal(rdm, coordinatesNew);CHKERRQ(ierr); 2576 ierr = VecDestroy(&coordinatesNew);CHKERRQ(ierr); 2577 ierr = PetscSectionDestroy(&coordSectionNew);CHKERRQ(ierr); 2578 PetscFunctionReturn(0); 2579 } 2580 2581 #undef __FUNCT__ 2582 #define __FUNCT__ "DMPlexCreateProcessSF" 2583 PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess) 2584 { 2585 PetscInt numRoots, numLeaves, l; 2586 const PetscInt *localPoints; 2587 const PetscSFNode *remotePoints; 2588 PetscInt *localPointsNew; 2589 PetscSFNode *remotePointsNew; 2590 PetscInt *ranks, *ranksNew; 2591 PetscErrorCode ierr; 2592 2593 PetscFunctionBegin; 2594 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2595 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranks);CHKERRQ(ierr); 2596 for (l = 0; l < numLeaves; ++l) { 2597 ranks[l] = remotePoints[l].rank; 2598 } 2599 ierr = PetscSortRemoveDupsInt(&numLeaves, ranks);CHKERRQ(ierr); 2600 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranksNew);CHKERRQ(ierr); 2601 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 2602 ierr = PetscMalloc(numLeaves * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 2603 for (l = 0; l < numLeaves; ++l) { 2604 ranksNew[l] = ranks[l]; 2605 localPointsNew[l] = l; 2606 remotePointsNew[l].index = 0; 2607 remotePointsNew[l].rank = ranksNew[l]; 2608 } 2609 ierr = PetscFree(ranks);CHKERRQ(ierr); 2610 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);CHKERRQ(ierr); 2611 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);CHKERRQ(ierr); 2612 ierr = PetscSFSetFromOptions(*sfProcess);CHKERRQ(ierr); 2613 ierr = PetscSFSetGraph(*sfProcess, 1, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2614 PetscFunctionReturn(0); 2615 } 2616 2617 #undef __FUNCT__ 2618 #define __FUNCT__ "CellRefinerCreateSF" 2619 PetscErrorCode CellRefinerCreateSF(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 2620 { 2621 PetscSF sf, sfNew, sfProcess; 2622 IS processRanks; 2623 MPI_Datatype depthType; 2624 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 2625 const PetscInt *localPoints, *neighbors; 2626 const PetscSFNode *remotePoints; 2627 PetscInt *localPointsNew; 2628 PetscSFNode *remotePointsNew; 2629 PetscInt *depthSizeOld, *rdepthSize, *rdepthSizeOld, *rdepthMaxOld, *rvStart, *rvStartNew, *reStart, *reStartNew, *rfStart, *rfStartNew, *rcStart, *rcStartNew; 2630 PetscInt depth, numNeighbors, pStartNew, pEndNew, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eStartNew, eEnd, eMax, r, n; 2631 PetscErrorCode ierr; 2632 2633 PetscFunctionBegin; 2634 ierr = DMPlexGetChart(rdm, &pStartNew, &pEndNew);CHKERRQ(ierr); 2635 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2636 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2637 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 2638 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2639 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2640 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 2641 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 2642 switch (refiner) { 2643 case 3: 2644 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 2645 cMax = PetscMin(cEnd, cMax); 2646 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 2647 fMax = PetscMin(fEnd, fMax); 2648 } 2649 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2650 ierr = DMGetPointSF(rdm, &sfNew);CHKERRQ(ierr); 2651 /* Caculate size of new SF */ 2652 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2653 if (numRoots < 0) PetscFunctionReturn(0); 2654 for (l = 0; l < numLeaves; ++l) { 2655 const PetscInt p = localPoints[l]; 2656 2657 switch (refiner) { 2658 case 1: 2659 /* Simplicial 2D */ 2660 if ((p >= vStart) && (p < vEnd)) { 2661 /* Old vertices stay the same */ 2662 ++numLeavesNew; 2663 } else if ((p >= fStart) && (p < fEnd)) { 2664 /* Old faces add new faces and vertex */ 2665 numLeavesNew += 2 + 1; 2666 } else if ((p >= cStart) && (p < cEnd)) { 2667 /* Old cells add new cells and interior faces */ 2668 numLeavesNew += 4 + 3; 2669 } 2670 break; 2671 case 2: 2672 /* Hex 2D */ 2673 if ((p >= vStart) && (p < vEnd)) { 2674 /* Old vertices stay the same */ 2675 ++numLeavesNew; 2676 } else if ((p >= fStart) && (p < fEnd)) { 2677 /* Old faces add new faces and vertex */ 2678 numLeavesNew += 2 + 1; 2679 } else if ((p >= cStart) && (p < cEnd)) { 2680 /* Old cells add new cells, interior faces, and vertex */ 2681 numLeavesNew += 4 + 4 + 1; 2682 } 2683 break; 2684 case 5: 2685 /* Simplicial 3D */ 2686 if ((p >= vStart) && (p < vEnd)) { 2687 /* Old vertices stay the same */ 2688 ++numLeavesNew; 2689 } else if ((p >= eStart) && (p < eEnd)) { 2690 /* Old edges add new edges and vertex */ 2691 numLeavesNew += 2 + 1; 2692 } else if ((p >= fStart) && (p < fEnd)) { 2693 /* Old faces add new faces and face edges */ 2694 numLeavesNew += 4 + 3; 2695 } else if ((p >= cStart) && (p < cEnd)) { 2696 /* Old cells add new cells and interior faces and edges */ 2697 numLeavesNew += 8 + 8 + 1; 2698 } 2699 break; 2700 case 6: 2701 /* Hex 3D */ 2702 if ((p >= vStart) && (p < vEnd)) { 2703 /* Old vertices stay the same */ 2704 ++numLeavesNew; 2705 } else if ((p >= eStart) && (p < eEnd)) { 2706 /* Old edges add new edges, and vertex */ 2707 numLeavesNew += 2 + 1; 2708 } else if ((p >= fStart) && (p < fEnd)) { 2709 /* Old faces add new faces, edges, and vertex */ 2710 numLeavesNew += 4 + 4 + 1; 2711 } else if ((p >= cStart) && (p < cEnd)) { 2712 /* Old cells add new cells, faces, edges, and vertex */ 2713 numLeavesNew += 8 + 12 + 6 + 1; 2714 } 2715 break; 2716 default: 2717 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2718 } 2719 } 2720 /* Communicate depthSizes for each remote rank */ 2721 ierr = DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);CHKERRQ(ierr); 2722 ierr = ISGetLocalSize(processRanks, &numNeighbors);CHKERRQ(ierr); 2723 ierr = PetscMalloc5((depth+1)*numNeighbors,PetscInt,&rdepthSize,numNeighbors,PetscInt,&rvStartNew,numNeighbors,PetscInt,&reStartNew,numNeighbors,PetscInt,&rfStartNew,numNeighbors,PetscInt,&rcStartNew);CHKERRQ(ierr); 2724 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); 2725 ierr = MPI_Type_contiguous(depth+1, MPIU_INT, &depthType);CHKERRQ(ierr); 2726 ierr = MPI_Type_commit(&depthType);CHKERRQ(ierr); 2727 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 2728 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 2729 for (n = 0; n < numNeighbors; ++n) { 2730 ierr = GetDepthStart_Private(depth, &rdepthSize[n*(depth+1)], &rcStartNew[n], &rfStartNew[n], &reStartNew[n], &rvStartNew[n]);CHKERRQ(ierr); 2731 } 2732 depthSizeOld[depth] = cMax; 2733 depthSizeOld[0] = vMax; 2734 depthSizeOld[depth-1] = fMax; 2735 depthSizeOld[1] = eMax; 2736 2737 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 2738 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 2739 2740 depthSizeOld[depth] = cEnd - cStart; 2741 depthSizeOld[0] = vEnd - vStart; 2742 depthSizeOld[depth-1] = fEnd - fStart; 2743 depthSizeOld[1] = eEnd - eStart; 2744 2745 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 2746 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 2747 for (n = 0; n < numNeighbors; ++n) { 2748 ierr = GetDepthStart_Private(depth, &rdepthSizeOld[n*(depth+1)], &rcStart[n], &rfStart[n], &reStart[n], &rvStart[n]);CHKERRQ(ierr); 2749 } 2750 ierr = MPI_Type_free(&depthType);CHKERRQ(ierr); 2751 ierr = PetscSFDestroy(&sfProcess);CHKERRQ(ierr); 2752 /* Calculate new point SF */ 2753 ierr = PetscMalloc(numLeavesNew * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 2754 ierr = PetscMalloc(numLeavesNew * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 2755 ierr = ISGetIndices(processRanks, &neighbors);CHKERRQ(ierr); 2756 for (l = 0, m = 0; l < numLeaves; ++l) { 2757 PetscInt p = localPoints[l]; 2758 PetscInt rp = remotePoints[l].index, n; 2759 PetscMPIInt rrank = remotePoints[l].rank; 2760 2761 ierr = PetscFindInt(rrank, numNeighbors, neighbors, &n);CHKERRQ(ierr); 2762 if (n < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %d", rrank); 2763 switch (refiner) { 2764 case 1: 2765 /* Simplicial 2D */ 2766 if ((p >= vStart) && (p < vEnd)) { 2767 /* Old vertices stay the same */ 2768 localPointsNew[m] = vStartNew + (p - vStart); 2769 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2770 remotePointsNew[m].rank = rrank; 2771 ++m; 2772 } else if ((p >= fStart) && (p < fEnd)) { 2773 /* Old faces add new faces and vertex */ 2774 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 2775 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 2776 remotePointsNew[m].rank = rrank; 2777 ++m; 2778 for (r = 0; r < 2; ++r, ++m) { 2779 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 2780 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 2781 remotePointsNew[m].rank = rrank; 2782 } 2783 } else if ((p >= cStart) && (p < cEnd)) { 2784 /* Old cells add new cells and interior faces */ 2785 for (r = 0; r < 4; ++r, ++m) { 2786 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2787 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2788 remotePointsNew[m].rank = rrank; 2789 } 2790 for (r = 0; r < 3; ++r, ++m) { 2791 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 2792 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*3 + r; 2793 remotePointsNew[m].rank = rrank; 2794 } 2795 } 2796 break; 2797 case 2: 2798 /* Hex 2D */ 2799 if ((p >= vStart) && (p < vEnd)) { 2800 /* Old vertices stay the same */ 2801 localPointsNew[m] = vStartNew + (p - vStart); 2802 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2803 remotePointsNew[m].rank = rrank; 2804 ++m; 2805 } else if ((p >= fStart) && (p < fEnd)) { 2806 /* Old faces add new faces and vertex */ 2807 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 2808 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 2809 remotePointsNew[m].rank = rrank; 2810 ++m; 2811 for (r = 0; r < 2; ++r, ++m) { 2812 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 2813 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 2814 remotePointsNew[m].rank = rrank; 2815 } 2816 } else if ((p >= cStart) && (p < cEnd)) { 2817 /* Old cells add new cells, interior faces, and vertex */ 2818 for (r = 0; r < 4; ++r, ++m) { 2819 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2820 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2821 remotePointsNew[m].rank = rrank; 2822 } 2823 for (r = 0; r < 4; ++r, ++m) { 2824 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 2825 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*4 + r; 2826 remotePointsNew[m].rank = rrank; 2827 } 2828 for (r = 0; r < 1; ++r, ++m) { 2829 localPointsNew[m] = vStartNew + (fEnd - fStart) + (p - cStart) + r; 2830 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1] + (rp - rcStart[n]) + r; 2831 remotePointsNew[m].rank = rrank; 2832 } 2833 } 2834 break; 2835 case 3: 2836 /* Hybrid simplicial 2D */ 2837 if ((p >= vStart) && (p < vEnd)) { 2838 /* Old vertices stay the same */ 2839 localPointsNew[m] = vStartNew + (p - vStart); 2840 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2841 remotePointsNew[m].rank = rrank; 2842 ++m; 2843 } else if ((p >= fStart) && (p < fMax)) { 2844 /* Old interior faces add new faces and vertex */ 2845 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 2846 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 2847 remotePointsNew[m].rank = rrank; 2848 ++m; 2849 for (r = 0; r < 2; ++r, ++m) { 2850 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 2851 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 2852 remotePointsNew[m].rank = rrank; 2853 } 2854 } else if ((p >= fMax) && (p < fEnd)) { 2855 /* Old hybrid faces stay the same */ 2856 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - fMax); 2857 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rdepthMaxOld[n*(depth+1)+depth-1]); 2858 remotePointsNew[m].rank = rrank; 2859 ++m; 2860 } else if ((p >= cStart) && (p < cMax)) { 2861 /* Old interior cells add new cells and interior faces */ 2862 for (r = 0; r < 4; ++r, ++m) { 2863 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2864 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2865 remotePointsNew[m].rank = rrank; 2866 } 2867 for (r = 0; r < 3; ++r, ++m) { 2868 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - cStart)*3 + r; 2869 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rcStart[n])*3 + r; 2870 remotePointsNew[m].rank = rrank; 2871 } 2872 } else if ((p >= cStart) && (p < cMax)) { 2873 /* Old hybrid cells add new cells and hybrid face */ 2874 for (r = 0; r < 2; ++r, ++m) { 2875 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2876 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2877 remotePointsNew[m].rank = rrank; 2878 } 2879 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 2880 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]); 2881 remotePointsNew[m].rank = rrank; 2882 ++m; 2883 } 2884 break; 2885 case 5: 2886 /* Simplicial 3D */ 2887 if ((p >= vStart) && (p < vEnd)) { 2888 /* Old vertices stay the same */ 2889 localPointsNew[m] = vStartNew + (p - vStart); 2890 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2891 remotePointsNew[m].rank = rrank; 2892 ++m; 2893 } else if ((p >= eStart) && (p < eEnd)) { 2894 /* Old edges add new edges and vertex */ 2895 for (r = 0; r < 2; ++r, ++m) { 2896 localPointsNew[m] = eStartNew + (p - eStart)*2 + r; 2897 remotePointsNew[m].index = reStartNew[n] + (rp - reStart[n])*2 + r; 2898 remotePointsNew[m].rank = rrank; 2899 } 2900 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - eStart); 2901 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - reStart[n]); 2902 remotePointsNew[m].rank = rrank; 2903 ++m; 2904 } else if ((p >= fStart) && (p < fEnd)) { 2905 /* Old faces add new faces and face edges */ 2906 for (r = 0; r < 4; ++r, ++m) { 2907 localPointsNew[m] = fStartNew + (p - fStart)*4 + r; 2908 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*4 + r; 2909 remotePointsNew[m].rank = rrank; 2910 } 2911 for (r = 0; r < 3; ++r, ++m) { 2912 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (p - fStart)*3 + r; 2913 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + (rp - rfStart[n])*3 + r; 2914 remotePointsNew[m].rank = rrank; 2915 } 2916 } else if ((p >= cStart) && (p < cEnd)) { 2917 /* Old cells add new cells and interior faces and edges */ 2918 for (r = 0; r < 8; ++r, ++m) { 2919 localPointsNew[m] = cStartNew + (p - cStart)*8 + r; 2920 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*8 + r; 2921 remotePointsNew[m].rank = rrank; 2922 } 2923 for (r = 0; r < 8; ++r, ++m) { 2924 localPointsNew[m] = fStartNew + (fEnd - fStart)*4 + (p - cStart)*8 + r; 2925 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*8 + r; 2926 remotePointsNew[m].rank = rrank; 2927 } 2928 for (r = 0; r < 1; ++r, ++m) { 2929 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (p - cStart)*1 + r; 2930 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + rdepthSizeOld[n*(depth+1)+depth-1]*3 + (rp - rcStart[n])*1 + r; 2931 remotePointsNew[m].rank = rrank; 2932 } 2933 } 2934 break; 2935 case 6: 2936 /* Hex 3D */ 2937 if ((p >= vStart) && (p < vEnd)) { 2938 /* Old vertices stay the same */ 2939 localPointsNew[m] = vStartNew + (p - vStart); 2940 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2941 remotePointsNew[m].rank = rrank; 2942 ++m; 2943 } else if ((p >= eStart) && (p < eEnd)) { 2944 /* Old edges add new edges and vertex */ 2945 for (r = 0; r < 2; ++r, ++m) { 2946 localPointsNew[m] = eStartNew + (p - eStart)*2 + r; 2947 remotePointsNew[m].index = reStartNew[n] + (rp - reStart[n])*2 + r; 2948 remotePointsNew[m].rank = rrank; 2949 } 2950 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - eStart); 2951 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - reStart[n]); 2952 remotePointsNew[m].rank = rrank; 2953 ++m; 2954 } else if ((p >= fStart) && (p < fEnd)) { 2955 /* Old faces add new faces, edges, and vertex */ 2956 for (r = 0; r < 4; ++r, ++m) { 2957 localPointsNew[m] = fStartNew + (p - fStart)*4 + r; 2958 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*4 + r; 2959 remotePointsNew[m].rank = rrank; 2960 } 2961 for (r = 0; r < 4; ++r, ++m) { 2962 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (p - fStart)*4 + r; 2963 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + (rp - rfStart[n])*4 + r; 2964 remotePointsNew[m].rank = rrank; 2965 } 2966 localPointsNew[m] = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (p - fStart); 2967 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + rdepthSizeOld[n*(depth+1)+1] + (rp - rfStart[n]); 2968 remotePointsNew[m].rank = rrank; 2969 ++m; 2970 } else if ((p >= cStart) && (p < cEnd)) { 2971 /* Old cells add new cells, faces, edges, and vertex */ 2972 for (r = 0; r < 8; ++r, ++m) { 2973 localPointsNew[m] = cStartNew + (p - cStart)*8 + r; 2974 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*8 + r; 2975 remotePointsNew[m].rank = rrank; 2976 } 2977 for (r = 0; r < 12; ++r, ++m) { 2978 localPointsNew[m] = fStartNew + (fEnd - fStart)*4 + (p - cStart)*12 + r; 2979 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*12 + r; 2980 remotePointsNew[m].rank = rrank; 2981 } 2982 for (r = 0; r < 6; ++r, ++m) { 2983 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (p - cStart)*6 + r; 2984 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*6 + r; 2985 remotePointsNew[m].rank = rrank; 2986 } 2987 for (r = 0; r < 1; ++r, ++m) { 2988 localPointsNew[m] = vStartNew + (eEnd - eStart) + (fEnd - fStart) + (p - cStart) + r; 2989 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+1] + rdepthSizeOld[n*(depth+1)+depth-1] + (rp - rcStart[n]) + r; 2990 remotePointsNew[m].rank = rrank; 2991 } 2992 } 2993 break; 2994 default: 2995 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2996 } 2997 } 2998 ierr = ISRestoreIndices(processRanks, &neighbors);CHKERRQ(ierr); 2999 ierr = ISDestroy(&processRanks);CHKERRQ(ierr); 3000 ierr = PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 3001 ierr = PetscFree5(rdepthSize,rvStartNew,reStartNew,rfStartNew,rcStartNew);CHKERRQ(ierr); 3002 ierr = PetscFree7(depthSizeOld,rdepthSizeOld,rdepthMaxOld,rvStart,reStart,rfStart,rcStart);CHKERRQ(ierr); 3003 PetscFunctionReturn(0); 3004 } 3005 3006 #undef __FUNCT__ 3007 #define __FUNCT__ "CellRefinerCreateLabels" 3008 PetscErrorCode CellRefinerCreateLabels(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 3009 { 3010 PetscInt numLabels, l; 3011 PetscInt depth, newp, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eStartNew, eEnd, eMax, r; 3012 PetscErrorCode ierr; 3013 3014 PetscFunctionBegin; 3015 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3016 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 3017 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 3018 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 3019 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3020 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 3021 ierr = DMPlexGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 3022 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 3023 switch (refiner) { 3024 case 3: 3025 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 3026 cMax = PetscMin(cEnd, cMax); 3027 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 3028 fMax = PetscMin(fEnd, fMax); 3029 } 3030 for (l = 0; l < numLabels; ++l) { 3031 DMLabel label, labelNew; 3032 const char *lname; 3033 PetscBool isDepth; 3034 IS valueIS; 3035 const PetscInt *values; 3036 PetscInt numValues, val; 3037 3038 ierr = DMPlexGetLabelName(dm, l, &lname);CHKERRQ(ierr); 3039 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 3040 if (isDepth) continue; 3041 ierr = DMPlexCreateLabel(rdm, lname);CHKERRQ(ierr); 3042 ierr = DMPlexGetLabel(dm, lname, &label);CHKERRQ(ierr); 3043 ierr = DMPlexGetLabel(rdm, lname, &labelNew);CHKERRQ(ierr); 3044 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 3045 ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr); 3046 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 3047 for (val = 0; val < numValues; ++val) { 3048 IS pointIS; 3049 const PetscInt *points; 3050 PetscInt numPoints, n; 3051 3052 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 3053 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 3054 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 3055 for (n = 0; n < numPoints; ++n) { 3056 const PetscInt p = points[n]; 3057 switch (refiner) { 3058 case 1: 3059 /* Simplicial 2D */ 3060 if ((p >= vStart) && (p < vEnd)) { 3061 /* Old vertices stay the same */ 3062 newp = vStartNew + (p - vStart); 3063 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3064 } else if ((p >= fStart) && (p < fEnd)) { 3065 /* Old faces add new faces and vertex */ 3066 newp = vStartNew + (vEnd - vStart) + (p - fStart); 3067 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3068 for (r = 0; r < 2; ++r) { 3069 newp = fStartNew + (p - fStart)*2 + r; 3070 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3071 } 3072 } else if ((p >= cStart) && (p < cEnd)) { 3073 /* Old cells add new cells and interior faces */ 3074 for (r = 0; r < 4; ++r) { 3075 newp = cStartNew + (p - cStart)*4 + r; 3076 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3077 } 3078 for (r = 0; r < 3; ++r) { 3079 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 3080 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3081 } 3082 } 3083 break; 3084 case 2: 3085 /* Hex 2D */ 3086 if ((p >= vStart) && (p < vEnd)) { 3087 /* Old vertices stay the same */ 3088 newp = vStartNew + (p - vStart); 3089 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3090 } else if ((p >= fStart) && (p < fEnd)) { 3091 /* Old faces add new faces and vertex */ 3092 newp = vStartNew + (vEnd - vStart) + (p - fStart); 3093 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3094 for (r = 0; r < 2; ++r) { 3095 newp = fStartNew + (p - fStart)*2 + r; 3096 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3097 } 3098 } else if ((p >= cStart) && (p < cEnd)) { 3099 /* Old cells add new cells and interior faces and vertex */ 3100 for (r = 0; r < 4; ++r) { 3101 newp = cStartNew + (p - cStart)*4 + r; 3102 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3103 } 3104 for (r = 0; r < 4; ++r) { 3105 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 3106 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3107 } 3108 newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (p - cStart); 3109 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3110 } 3111 break; 3112 case 3: 3113 /* Hybrid simplicial 2D */ 3114 if ((p >= vStart) && (p < vEnd)) { 3115 /* Old vertices stay the same */ 3116 newp = vStartNew + (p - vStart); 3117 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3118 } else if ((p >= fStart) && (p < fMax)) { 3119 /* Old interior faces add new faces and vertex */ 3120 newp = vStartNew + (vEnd - vStart) + (p - fStart); 3121 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3122 for (r = 0; r < 2; ++r) { 3123 newp = fStartNew + (p - fStart)*2 + r; 3124 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3125 } 3126 } else if ((p >= fMax) && (p < fEnd)) { 3127 /* Old hybrid faces stay the same */ 3128 newp = fStartNew + (fMax - fStart)*2 + (p - fMax); 3129 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3130 } else if ((p >= cStart) && (p < cMax)) { 3131 /* Old interior cells add new cells and interior faces */ 3132 for (r = 0; r < 4; ++r) { 3133 newp = cStartNew + (p - cStart)*4 + r; 3134 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3135 } 3136 for (r = 0; r < 3; ++r) { 3137 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 3138 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3139 } 3140 } else if ((p >= cMax) && (p < cEnd)) { 3141 /* Old hybrid cells add new cells and hybrid face */ 3142 for (r = 0; r < 2; ++r) { 3143 newp = cStartNew + (cMax - cStart)*4 + (p - cMax)*2 + r; 3144 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3145 } 3146 newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 3147 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3148 } 3149 break; 3150 case 5: 3151 /* Simplicial 3D */ 3152 if ((p >= vStart) && (p < vEnd)) { 3153 /* Old vertices stay the same */ 3154 newp = vStartNew + (p - vStart); 3155 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3156 } else if ((p >= eStart) && (p < eEnd)) { 3157 /* Old edges add new edges and vertex */ 3158 for (r = 0; r < 2; ++r) { 3159 newp = eStartNew + (p - eStart)*2 + r; 3160 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3161 } 3162 newp = vStartNew + (vEnd - vStart) + (p - eStart); 3163 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3164 } else if ((p >= fStart) && (p < fEnd)) { 3165 /* Old faces add new faces and edges */ 3166 for (r = 0; r < 4; ++r) { 3167 newp = fStartNew + (p - fStart)*4 + r; 3168 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3169 } 3170 for (r = 0; r < 3; ++r) { 3171 newp = eStartNew + (eEnd - eStart)*2 + (p - fStart)*3 + r; 3172 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3173 } 3174 } else if ((p >= cStart) && (p < cEnd)) { 3175 /* Old cells add new cells and interior faces and edges */ 3176 for (r = 0; r < 8; ++r) { 3177 newp = cStartNew + (p - cStart)*8 + r; 3178 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3179 } 3180 for (r = 0; r < 8; ++r) { 3181 newp = fStartNew + (fEnd - fStart)*4 + (p - cStart)*8 + r; 3182 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3183 } 3184 for (r = 0; r < 1; ++r) { 3185 newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (p - cStart)*1 + r; 3186 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3187 } 3188 } 3189 break; 3190 case 6: 3191 /* Hex 3D */ 3192 if ((p >= vStart) && (p < vEnd)) { 3193 /* Old vertices stay the same */ 3194 newp = vStartNew + (p - vStart); 3195 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3196 } else if ((p >= fStart) && (p < fEnd)) { 3197 /* Old edges add new edges and vertex */ 3198 for (r = 0; r < 2; ++r) { 3199 newp = eStartNew + (p - eStart)*2 + r; 3200 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3201 } 3202 newp = vStartNew + (vEnd - vStart) + (p - eStart); 3203 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3204 } else if ((p >= fStart) && (p < fEnd)) { 3205 /* Old faces add new faces, edges, and vertex */ 3206 for (r = 0; r < 4; ++r) { 3207 newp = fStartNew + (p - fStart)*4 + r; 3208 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3209 } 3210 for (r = 0; r < 4; ++r) { 3211 newp = eStartNew + (eEnd - eStart)*2 + (p - fStart)*4 + r; 3212 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3213 } 3214 newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (p - fStart); 3215 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3216 } else if ((p >= cStart) && (p < cEnd)) { 3217 /* Old cells add new cells, faces, edges, and vertex */ 3218 for (r = 0; r < 8; ++r) { 3219 newp = cStartNew + (p - cStart)*8 + r; 3220 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3221 } 3222 for (r = 0; r < 12; ++r) { 3223 newp = fStartNew + (fEnd - fStart)*4 + (p - cStart)*12 + r; 3224 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3225 } 3226 for (r = 0; r < 6; ++r) { 3227 newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (p - cStart)*6 + r; 3228 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3229 } 3230 newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (p - cStart); 3231 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3232 } 3233 break; 3234 default: 3235 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 3236 } 3237 } 3238 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 3239 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 3240 } 3241 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 3242 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 3243 if (0) { 3244 ierr = PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);CHKERRQ(ierr); 3245 ierr = DMLabelView(labelNew, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 3246 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 3247 } 3248 } 3249 PetscFunctionReturn(0); 3250 } 3251 3252 #undef __FUNCT__ 3253 #define __FUNCT__ "DMPlexRefineUniform_Internal" 3254 /* This will only work for interpolated meshes */ 3255 PetscErrorCode DMPlexRefineUniform_Internal(DM dm, CellRefiner cellRefiner, DM *dmRefined) 3256 { 3257 DM rdm; 3258 PetscInt *depthSize; 3259 PetscInt dim, depth = 0, d, pStart = 0, pEnd = 0; 3260 PetscErrorCode ierr; 3261 3262 PetscFunctionBegin; 3263 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &rdm);CHKERRQ(ierr); 3264 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3265 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 3266 ierr = DMPlexSetDimension(rdm, dim);CHKERRQ(ierr); 3267 /* Calculate number of new points of each depth */ 3268 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3269 ierr = PetscMalloc((depth+1) * sizeof(PetscInt), &depthSize);CHKERRQ(ierr); 3270 ierr = PetscMemzero(depthSize, (depth+1) * sizeof(PetscInt));CHKERRQ(ierr); 3271 ierr = CellRefinerGetSizes(cellRefiner, dm, depthSize);CHKERRQ(ierr); 3272 /* Step 1: Set chart */ 3273 for (d = 0; d <= depth; ++d) pEnd += depthSize[d]; 3274 ierr = DMPlexSetChart(rdm, pStart, pEnd);CHKERRQ(ierr); 3275 /* Step 2: Set cone/support sizes */ 3276 ierr = CellRefinerSetConeSizes(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3277 /* Step 3: Setup refined DM */ 3278 ierr = DMSetUp(rdm);CHKERRQ(ierr); 3279 /* Step 4: Set cones and supports */ 3280 ierr = CellRefinerSetCones(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3281 /* Step 5: Stratify */ 3282 ierr = DMPlexStratify(rdm);CHKERRQ(ierr); 3283 /* Step 6: Set coordinates for vertices */ 3284 ierr = CellRefinerSetCoordinates(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3285 /* Step 7: Create pointSF */ 3286 ierr = CellRefinerCreateSF(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3287 /* Step 8: Create labels */ 3288 ierr = CellRefinerCreateLabels(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3289 ierr = PetscFree(depthSize);CHKERRQ(ierr); 3290 3291 *dmRefined = rdm; 3292 PetscFunctionReturn(0); 3293 } 3294 3295 #undef __FUNCT__ 3296 #define __FUNCT__ "DMPlexSetRefinementUniform" 3297 PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform) 3298 { 3299 DM_Plex *mesh = (DM_Plex*) dm->data; 3300 3301 PetscFunctionBegin; 3302 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3303 mesh->refinementUniform = refinementUniform; 3304 PetscFunctionReturn(0); 3305 } 3306 3307 #undef __FUNCT__ 3308 #define __FUNCT__ "DMPlexGetRefinementUniform" 3309 PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform) 3310 { 3311 DM_Plex *mesh = (DM_Plex*) dm->data; 3312 3313 PetscFunctionBegin; 3314 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3315 PetscValidPointer(refinementUniform, 2); 3316 *refinementUniform = mesh->refinementUniform; 3317 PetscFunctionReturn(0); 3318 } 3319 3320 #undef __FUNCT__ 3321 #define __FUNCT__ "DMPlexSetRefinementLimit" 3322 PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit) 3323 { 3324 DM_Plex *mesh = (DM_Plex*) dm->data; 3325 3326 PetscFunctionBegin; 3327 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3328 mesh->refinementLimit = refinementLimit; 3329 PetscFunctionReturn(0); 3330 } 3331 3332 #undef __FUNCT__ 3333 #define __FUNCT__ "DMPlexGetRefinementLimit" 3334 PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit) 3335 { 3336 DM_Plex *mesh = (DM_Plex*) dm->data; 3337 3338 PetscFunctionBegin; 3339 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3340 PetscValidPointer(refinementLimit, 2); 3341 /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */ 3342 *refinementLimit = mesh->refinementLimit; 3343 PetscFunctionReturn(0); 3344 } 3345 3346 #undef __FUNCT__ 3347 #define __FUNCT__ "DMPlexGetCellRefiner_Internal" 3348 PetscErrorCode DMPlexGetCellRefiner_Internal(DM dm, CellRefiner *cellRefiner) 3349 { 3350 PetscInt dim, cStart, coneSize, cMax; 3351 PetscErrorCode ierr; 3352 3353 PetscFunctionBegin; 3354 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 3355 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 3356 ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr); 3357 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 3358 switch (dim) { 3359 case 2: 3360 switch (coneSize) { 3361 case 3: 3362 if (cMax >= 0) *cellRefiner = 3; /* Hybrid */ 3363 else *cellRefiner = 1; /* Triangular */ 3364 break; 3365 case 4: 3366 if (cMax >= 0) *cellRefiner = 4; /* Hybrid */ 3367 else *cellRefiner = 2; /* Quadrilateral */ 3368 break; 3369 default: 3370 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown coneSize %d in dimension %d for cell refiner", coneSize, dim); 3371 } 3372 break; 3373 case 3: 3374 switch (coneSize) { 3375 case 4: 3376 if (cMax >= 0) *cellRefiner = 7; /* Hybrid */ 3377 else *cellRefiner = 5; /* Tetrahedral */ 3378 break; 3379 case 6: 3380 if (cMax >= 0) *cellRefiner = 8; /* Hybrid */ 3381 else *cellRefiner = 6; /* hexahedral */ 3382 break; 3383 default: 3384 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown coneSize %d in dimension %d for cell refiner", coneSize, dim); 3385 } 3386 break; 3387 default: 3388 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown dimension %d for cell refiner", dim); 3389 } 3390 PetscFunctionReturn(0); 3391 } 3392