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