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