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 ierr = PetscFree(supportRef);CHKERRQ(ierr); 999 break; 1000 case 3: 1001 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 1002 cMax = PetscMin(cEnd, cMax); 1003 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 1004 fMax = PetscMin(fEnd, fMax); 1005 /* Interior cells have 3 faces */ 1006 for (c = cStart; c < cMax; ++c) { 1007 const PetscInt newp = cStartNew + (c - cStart)*4; 1008 const PetscInt *cone, *ornt; 1009 PetscInt coneNew[3], orntNew[3]; 1010 1011 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1012 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1013 /* A triangle */ 1014 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 1015 orntNew[0] = ornt[0]; 1016 coneNew[1] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 2; 1017 orntNew[1] = -2; 1018 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 1019 orntNew[2] = ornt[2]; 1020 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1021 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1022 #if 1 1023 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); 1024 for (p = 0; p < 3; ++p) { 1025 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); 1026 } 1027 #endif 1028 /* B triangle */ 1029 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 1030 orntNew[0] = ornt[0]; 1031 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 1032 orntNew[1] = ornt[1]; 1033 coneNew[2] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 0; 1034 orntNew[2] = -2; 1035 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1036 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1037 #if 1 1038 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); 1039 for (p = 0; p < 3; ++p) { 1040 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); 1041 } 1042 #endif 1043 /* C triangle */ 1044 coneNew[0] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 1; 1045 orntNew[0] = -2; 1046 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 1047 orntNew[1] = ornt[1]; 1048 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 1049 orntNew[2] = ornt[2]; 1050 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 1051 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 1052 #if 1 1053 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); 1054 for (p = 0; p < 3; ++p) { 1055 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); 1056 } 1057 #endif 1058 /* D triangle */ 1059 coneNew[0] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 0; 1060 orntNew[0] = 0; 1061 coneNew[1] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 1; 1062 orntNew[1] = 0; 1063 coneNew[2] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 2; 1064 orntNew[2] = 0; 1065 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 1066 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 1067 #if 1 1068 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); 1069 for (p = 0; p < 3; ++p) { 1070 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); 1071 } 1072 #endif 1073 } 1074 /* 1075 2----3----3 1076 | | 1077 | B | 1078 | | 1079 0----4--- 1 1080 | | 1081 | A | 1082 | | 1083 0----2----1 1084 */ 1085 /* Hybrid cells have 4 faces */ 1086 for (c = cMax; c < cEnd; ++c) { 1087 const PetscInt newp = cStartNew + (cMax - cStart)*4 + (c - cMax)*2; 1088 const PetscInt *cone, *ornt; 1089 PetscInt coneNew[4], orntNew[4]; 1090 1091 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1092 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1093 /* A quad */ 1094 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 1095 orntNew[0] = ornt[0]; 1096 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 1097 orntNew[1] = ornt[1]; 1098 coneNew[2] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (cone[2] - fMax); 1099 orntNew[2] = 0; 1100 coneNew[3] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 1101 orntNew[3] = 0; 1102 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1103 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1104 #if 1 1105 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); 1106 for (p = 0; p < 4; ++p) { 1107 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); 1108 } 1109 #endif 1110 /* B quad */ 1111 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 1112 orntNew[0] = ornt[0]; 1113 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 1114 orntNew[1] = ornt[1]; 1115 coneNew[2] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 1116 orntNew[2] = 0; 1117 coneNew[3] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (cone[3] - fMax); 1118 orntNew[3] = 0; 1119 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1120 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1121 #if 1 1122 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); 1123 for (p = 0; p < 4; ++p) { 1124 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); 1125 } 1126 #endif 1127 } 1128 /* Interior split faces have 2 vertices and the same cells as the parent */ 1129 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 1130 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 1131 for (f = fStart; f < fMax; ++f) { 1132 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 1133 1134 for (r = 0; r < 2; ++r) { 1135 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 1136 const PetscInt *cone, *ornt, *support; 1137 PetscInt coneNew[2], coneSize, c, supportSize, s; 1138 1139 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1140 coneNew[0] = vStartNew + (cone[0] - vStart); 1141 coneNew[1] = vStartNew + (cone[1] - vStart); 1142 coneNew[(r+1)%2] = newv; 1143 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1144 #if 1 1145 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1146 for (p = 0; p < 2; ++p) { 1147 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); 1148 } 1149 #endif 1150 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1151 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1152 for (s = 0; s < supportSize; ++s) { 1153 if (support[s] >= cMax) { 1154 supportRef[s] = cStartNew + (cMax - cStart)*4 + (support[s] - cMax)*2 + r; 1155 } else { 1156 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1157 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1158 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1159 for (c = 0; c < coneSize; ++c) { 1160 if (cone[c] == f) break; 1161 } 1162 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (ornt[c] < 0 ? (c+1-r)%3 : (c+r)%3); 1163 } 1164 } 1165 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1166 #if 1 1167 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1168 for (p = 0; p < supportSize; ++p) { 1169 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); 1170 } 1171 #endif 1172 } 1173 } 1174 /* Interior cell faces have 2 vertices and 2 cells */ 1175 for (c = cStart; c < cMax; ++c) { 1176 const PetscInt *cone; 1177 1178 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1179 for (r = 0; r < 3; ++r) { 1180 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + r; 1181 PetscInt coneNew[2]; 1182 PetscInt supportNew[2]; 1183 1184 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 1185 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - fStart); 1186 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1187 #if 1 1188 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1189 for (p = 0; p < 2; ++p) { 1190 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); 1191 } 1192 #endif 1193 supportNew[0] = (c - cStart)*4 + (r+1)%3; 1194 supportNew[1] = (c - cStart)*4 + 3; 1195 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1196 #if 1 1197 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1198 for (p = 0; p < 2; ++p) { 1199 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); 1200 } 1201 #endif 1202 } 1203 } 1204 /* Interior hybrid faces have 2 vertices and the same cells */ 1205 for (f = fMax; f < fEnd; ++f) { 1206 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (f - fMax); 1207 const PetscInt *cone; 1208 const PetscInt *support; 1209 PetscInt coneNew[2]; 1210 PetscInt supportNew[2]; 1211 PetscInt size, s, r; 1212 1213 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1214 coneNew[0] = vStartNew + (cone[0] - vStart); 1215 coneNew[1] = vStartNew + (cone[1] - vStart); 1216 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1217 #if 1 1218 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1219 for (p = 0; p < 2; ++p) { 1220 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); 1221 } 1222 #endif 1223 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 1224 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1225 for (s = 0; s < size; ++s) { 1226 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1227 for (r = 0; r < 2; ++r) { 1228 if (cone[r+2] == f) break; 1229 } 1230 supportNew[s] = (cMax - cStart)*4 + (support[s] - cMax)*2 + r; 1231 } 1232 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1233 #if 1 1234 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1235 for (p = 0; p < size; ++p) { 1236 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); 1237 } 1238 #endif 1239 } 1240 /* Cell hybrid faces have 2 vertices and 2 cells */ 1241 for (c = cMax; c < cEnd; ++c) { 1242 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 1243 const PetscInt *cone; 1244 PetscInt coneNew[2]; 1245 PetscInt supportNew[2]; 1246 1247 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1248 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[0] - fStart); 1249 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[1] - fStart); 1250 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1251 #if 1 1252 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1253 for (p = 0; p < 2; ++p) { 1254 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); 1255 } 1256 #endif 1257 supportNew[0] = (cMax - cStart)*4 + (c - cMax)*2 + 0; 1258 supportNew[1] = (cMax - cStart)*4 + (c - cMax)*2 + 1; 1259 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1260 #if 1 1261 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1262 for (p = 0; p < 2; ++p) { 1263 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); 1264 } 1265 #endif 1266 } 1267 /* Old vertices have identical supports */ 1268 for (v = vStart; v < vEnd; ++v) { 1269 const PetscInt newp = vStartNew + (v - vStart); 1270 const PetscInt *support, *cone; 1271 PetscInt size, s; 1272 1273 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 1274 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 1275 for (s = 0; s < size; ++s) { 1276 if (support[s] >= fMax) { 1277 supportRef[s] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (support[s] - fMax); 1278 } else { 1279 PetscInt r = 0; 1280 1281 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1282 if (cone[1] == v) r = 1; 1283 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 1284 } 1285 } 1286 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1287 #if 1 1288 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1289 for (p = 0; p < size; ++p) { 1290 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); 1291 } 1292 #endif 1293 } 1294 /* Face vertices have 2 + (2 interior, 1 hybrid) supports */ 1295 for (f = fStart; f < fMax; ++f) { 1296 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 1297 const PetscInt *cone, *support; 1298 PetscInt size, newSize = 2, s; 1299 1300 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 1301 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1302 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 1303 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 1304 for (s = 0; s < size; ++s) { 1305 PetscInt r = 0; 1306 1307 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1308 if (support[s] >= cMax) { 1309 supportRef[newSize+0] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (support[s] - cMax); 1310 1311 newSize += 1; 1312 } else { 1313 if (cone[1] == f) r = 1; 1314 else if (cone[2] == f) r = 2; 1315 supportRef[newSize+0] = fStartNew + (fMax - fStart)*2 + (support[s] - cStart)*3 + (r+2)%3; 1316 supportRef[newSize+1] = fStartNew + (fMax - fStart)*2 + (support[s] - cStart)*3 + r; 1317 1318 newSize += 2; 1319 } 1320 } 1321 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1322 #if 1 1323 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1324 for (p = 0; p < newSize; ++p) { 1325 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); 1326 } 1327 #endif 1328 } 1329 ierr = PetscFree(supportRef);CHKERRQ(ierr); 1330 break; 1331 case 5: 1332 /* Simplicial 3D */ 1333 /* All cells have 4 faces: Tet face order is prescribed in DMPlexGetFaces_Internal() */ 1334 ierr = DMPlexGetRawFaces_Internal(dm, 3, 4, cellInd, NULL, NULL, &faces);CHKERRQ(ierr); 1335 for (c = cStart; c < cEnd; ++c) { 1336 const PetscInt newp = cStartNew + (c - cStart)*8; 1337 const PetscInt *cone, *ornt; 1338 PetscInt coneNew[4], orntNew[4]; 1339 1340 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1341 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1342 /* A tetrahedron: {0, a, c, d} */ 1343 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetTriSubface_Static(ornt[0], 0); /* A */ 1344 orntNew[0] = ornt[0]; 1345 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetTriSubface_Static(ornt[1], 0); /* A */ 1346 orntNew[1] = ornt[1]; 1347 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetTriSubface_Static(ornt[2], 0); /* A */ 1348 orntNew[2] = ornt[2]; 1349 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 0; 1350 orntNew[3] = 0; 1351 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1352 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1353 #if 1 1354 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); 1355 for (p = 0; p < 4; ++p) { 1356 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); 1357 } 1358 #endif 1359 /* B tetrahedron: {a, 1, b, e} */ 1360 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetTriSubface_Static(ornt[0], 1); /* B */ 1361 orntNew[0] = ornt[0]; 1362 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetTriSubface_Static(ornt[1], 2); /* C */ 1363 orntNew[1] = ornt[1]; 1364 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 1; 1365 orntNew[2] = 0; 1366 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetTriSubface_Static(ornt[3], 1); /* B */ 1367 orntNew[3] = ornt[3]; 1368 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1369 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1370 #if 1 1371 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); 1372 for (p = 0; p < 4; ++p) { 1373 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); 1374 } 1375 #endif 1376 /* C tetrahedron: {c, b, 2, f} */ 1377 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetTriSubface_Static(ornt[0], 2); /* C */ 1378 orntNew[0] = ornt[0]; 1379 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 2; 1380 orntNew[1] = 0; 1381 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetTriSubface_Static(ornt[2], 1); /* B */ 1382 orntNew[2] = ornt[2]; 1383 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetTriSubface_Static(ornt[3], 0); /* A */ 1384 orntNew[3] = ornt[3]; 1385 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 1386 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 1387 #if 1 1388 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); 1389 for (p = 0; p < 4; ++p) { 1390 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); 1391 } 1392 #endif 1393 /* D tetrahedron: {d, e, f, 3} */ 1394 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 3; 1395 orntNew[0] = 0; 1396 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetTriSubface_Static(ornt[1], 1); /* B */ 1397 orntNew[1] = ornt[1]; 1398 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetTriSubface_Static(ornt[2], 2); /* C */ 1399 orntNew[2] = ornt[2]; 1400 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetTriSubface_Static(ornt[3], 2); /* C */ 1401 orntNew[3] = ornt[3]; 1402 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 1403 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 1404 #if 1 1405 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); 1406 for (p = 0; p < 4; ++p) { 1407 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); 1408 } 1409 #endif 1410 /* A' tetrahedron: {d, a, c, f} */ 1411 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 0; 1412 orntNew[0] = -3; 1413 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 4; 1414 orntNew[1] = 0; 1415 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + 3; 1416 orntNew[2] = ornt[2] < 0 ? -((-(ornt[2]+1)+2)%3+1) : (ornt[2]+2)%3; 1417 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 5; 1418 orntNew[3] = 0; 1419 ierr = DMPlexSetCone(rdm, newp+4, coneNew);CHKERRQ(ierr); 1420 ierr = DMPlexSetConeOrientation(rdm, newp+4, orntNew);CHKERRQ(ierr); 1421 #if 1 1422 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); 1423 for (p = 0; p < 4; ++p) { 1424 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); 1425 } 1426 #endif 1427 /* B' tetrahedron: {e, b, a, f} */ 1428 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 1; 1429 orntNew[0] = -3; 1430 coneNew[1] = fStartNew + (cone[3] - fStart)*4 + 3; 1431 orntNew[1] = ornt[3] < 0 ? -((-(ornt[3]+1)+1)%3+1) : (ornt[3]+1)%3; 1432 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 6; 1433 orntNew[2] = 0; 1434 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 7; 1435 orntNew[3] = 0; 1436 ierr = DMPlexSetCone(rdm, newp+5, coneNew);CHKERRQ(ierr); 1437 ierr = DMPlexSetConeOrientation(rdm, newp+5, orntNew);CHKERRQ(ierr); 1438 #if 1 1439 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); 1440 for (p = 0; p < 4; ++p) { 1441 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); 1442 } 1443 #endif 1444 /* C' tetrahedron: {b, f, c, a} */ 1445 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 2; 1446 orntNew[0] = -3; 1447 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 7; 1448 orntNew[1] = -2; 1449 coneNew[2] = fStartNew + (cone[0] - fStart)*4 + 3; 1450 orntNew[2] = ornt[0] < 0 ? (-(ornt[0]+1)+1)%3 : -((ornt[0]+1)%3+1); 1451 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 5; 1452 orntNew[3] = -1; 1453 ierr = DMPlexSetCone(rdm, newp+6, coneNew);CHKERRQ(ierr); 1454 ierr = DMPlexSetConeOrientation(rdm, newp+6, orntNew);CHKERRQ(ierr); 1455 #if 1 1456 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); 1457 for (p = 0; p < 4; ++p) { 1458 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); 1459 } 1460 #endif 1461 /* D' tetrahedron: {f, e, d, a} */ 1462 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 3; 1463 orntNew[0] = -3; 1464 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 6; 1465 orntNew[1] = -3; 1466 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 4; 1467 orntNew[2] = -2; 1468 coneNew[3] = fStartNew + (cone[1] - fStart)*4 + 3; 1469 orntNew[3] = ornt[2]; 1470 ierr = DMPlexSetCone(rdm, newp+7, coneNew);CHKERRQ(ierr); 1471 ierr = DMPlexSetConeOrientation(rdm, newp+7, orntNew);CHKERRQ(ierr); 1472 #if 1 1473 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); 1474 for (p = 0; p < 4; ++p) { 1475 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); 1476 } 1477 #endif 1478 } 1479 /* Split faces have 3 edges and the same cells as the parent */ 1480 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 1481 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 1482 for (f = fStart; f < fEnd; ++f) { 1483 const PetscInt newp = fStartNew + (f - fStart)*4; 1484 const PetscInt *cone, *ornt, *support; 1485 PetscInt coneNew[3], orntNew[3], coneSize, supportSize, s; 1486 1487 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1488 ierr = DMPlexGetConeOrientation(dm, f, &ornt);CHKERRQ(ierr); 1489 /* A triangle */ 1490 coneNew[0] = eStartNew + (cone[0] - eStart)*2 + (ornt[0] < 0 ? 1 : 0); 1491 orntNew[0] = ornt[0]; 1492 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 2; 1493 orntNew[1] = -2; 1494 coneNew[2] = eStartNew + (cone[2] - eStart)*2 + (ornt[2] < 0 ? 0 : 1); 1495 orntNew[2] = ornt[2]; 1496 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1497 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1498 #if 1 1499 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); 1500 for (p = 0; p < 3; ++p) { 1501 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); 1502 } 1503 #endif 1504 /* B triangle */ 1505 coneNew[0] = eStartNew + (cone[0] - eStart)*2 + (ornt[0] < 0 ? 0 : 1); 1506 orntNew[0] = ornt[0]; 1507 coneNew[1] = eStartNew + (cone[1] - eStart)*2 + (ornt[1] < 0 ? 1 : 0); 1508 orntNew[1] = ornt[1]; 1509 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 0; 1510 orntNew[2] = -2; 1511 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1512 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1513 #if 1 1514 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); 1515 for (p = 0; p < 3; ++p) { 1516 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); 1517 } 1518 #endif 1519 /* C triangle */ 1520 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 1; 1521 orntNew[0] = -2; 1522 coneNew[1] = eStartNew + (cone[1] - eStart)*2 + (ornt[1] < 0 ? 0 : 1); 1523 orntNew[1] = ornt[1]; 1524 coneNew[2] = eStartNew + (cone[2] - eStart)*2 + (ornt[2] < 0 ? 1 : 0); 1525 orntNew[2] = ornt[2]; 1526 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 1527 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 1528 #if 1 1529 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); 1530 for (p = 0; p < 3; ++p) { 1531 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); 1532 } 1533 #endif 1534 /* D triangle */ 1535 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 0; 1536 orntNew[0] = 0; 1537 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 1; 1538 orntNew[1] = 0; 1539 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 2; 1540 orntNew[2] = 0; 1541 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 1542 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 1543 #if 1 1544 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); 1545 for (p = 0; p < 3; ++p) { 1546 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); 1547 } 1548 #endif 1549 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1550 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1551 for (r = 0; r < 4; ++r) { 1552 for (s = 0; s < supportSize; ++s) { 1553 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1554 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1555 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1556 for (c = 0; c < coneSize; ++c) { 1557 if (cone[c] == f) break; 1558 } 1559 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])); 1560 } 1561 ierr = DMPlexSetSupport(rdm, newp+r, supportRef);CHKERRQ(ierr); 1562 #if 1 1563 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1564 for (p = 0; p < supportSize; ++p) { 1565 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); 1566 } 1567 #endif 1568 } 1569 } 1570 /* Interior faces have 3 edges and 2 cells */ 1571 for (c = cStart; c < cEnd; ++c) { 1572 PetscInt newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8; 1573 const PetscInt *cone, *ornt; 1574 PetscInt coneNew[3], orntNew[3]; 1575 PetscInt supportNew[2]; 1576 1577 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1578 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1579 /* Face A: {c, a, d} */ 1580 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+0)%3 : (ornt[0]+2)%3); 1581 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1582 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+0)%3 : (ornt[1]+2)%3); 1583 orntNew[1] = ornt[1] < 0 ? -2 : 0; 1584 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+0)%3 : (ornt[2]+2)%3); 1585 orntNew[2] = ornt[2] < 0 ? -2 : 0; 1586 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1587 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1588 #if 1 1589 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1590 for (p = 0; p < 3; ++p) { 1591 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); 1592 } 1593 #endif 1594 supportNew[0] = (c - cStart)*8 + 0; 1595 supportNew[1] = (c - cStart)*8 + 0+4; 1596 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1597 #if 1 1598 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1599 for (p = 0; p < 2; ++p) { 1600 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); 1601 } 1602 #endif 1603 ++newp; 1604 /* Face B: {a, b, e} */ 1605 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+2)%3 : (ornt[0]+0)%3); 1606 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1607 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+2)%3 : (ornt[3]+0)%3); 1608 orntNew[1] = ornt[3] < 0 ? -2 : 0; 1609 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+1)%3 : (ornt[1]+1)%3); 1610 orntNew[2] = ornt[1] < 0 ? -2 : 0; 1611 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1612 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1613 #if 1 1614 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); 1615 for (p = 0; p < 3; ++p) { 1616 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); 1617 } 1618 #endif 1619 supportNew[0] = (c - cStart)*8 + 1; 1620 supportNew[1] = (c - cStart)*8 + 1+4; 1621 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1622 #if 1 1623 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1624 for (p = 0; p < 2; ++p) { 1625 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); 1626 } 1627 #endif 1628 ++newp; 1629 /* Face C: {c, f, b} */ 1630 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+2)%3 : (ornt[2]+0)%3); 1631 orntNew[0] = ornt[2] < 0 ? -2 : 0; 1632 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+0)%3 : (ornt[3]+2)%3); 1633 orntNew[1] = ornt[3] < 0 ? -2 : 0; 1634 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+1)%3 : (ornt[0]+1)%3); 1635 orntNew[2] = ornt[0] < 0 ? -2 : 0; 1636 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1637 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1638 #if 1 1639 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1640 for (p = 0; p < 3; ++p) { 1641 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); 1642 } 1643 #endif 1644 supportNew[0] = (c - cStart)*8 + 2; 1645 supportNew[1] = (c - cStart)*8 + 2+4; 1646 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1647 #if 1 1648 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1649 for (p = 0; p < 2; ++p) { 1650 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); 1651 } 1652 #endif 1653 ++newp; 1654 /* Face D: {d, e, f} */ 1655 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+2)%3 : (ornt[1]+0)%3); 1656 orntNew[0] = ornt[1] < 0 ? -2 : 0; 1657 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+1)%3 : (ornt[3]+1)%3); 1658 orntNew[1] = ornt[3] < 0 ? -2 : 0; 1659 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+1)%3 : (ornt[2]+1)%3); 1660 orntNew[2] = ornt[2] < 0 ? -2 : 0; 1661 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1662 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1663 #if 1 1664 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1665 for (p = 0; p < 3; ++p) { 1666 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); 1667 } 1668 #endif 1669 supportNew[0] = (c - cStart)*8 + 3; 1670 supportNew[1] = (c - cStart)*8 + 3+4; 1671 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1672 #if 1 1673 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1674 for (p = 0; p < 2; ++p) { 1675 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); 1676 } 1677 #endif 1678 ++newp; 1679 /* Face E: {d, f, a} */ 1680 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+1)%3 : (ornt[2]+1)%3); 1681 orntNew[0] = ornt[2] < 0 ? 0 : -2; 1682 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1683 orntNew[1] = 0; 1684 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+0)%3 : (ornt[1]+2)%3); 1685 orntNew[2] = ornt[1] < 0 ? -2 : 0; 1686 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1687 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1688 #if 1 1689 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1690 for (p = 0; p < 3; ++p) { 1691 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); 1692 } 1693 #endif 1694 supportNew[0] = (c - cStart)*8 + 0+4; 1695 supportNew[1] = (c - cStart)*8 + 3+4; 1696 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1697 #if 1 1698 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1699 for (p = 0; p < 2; ++p) { 1700 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); 1701 } 1702 #endif 1703 ++newp; 1704 /* Face F: {c, a, f} */ 1705 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+0)%3 : (ornt[0]+2)%3); 1706 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1707 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1708 orntNew[1] = -2; 1709 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+2)%3 : (ornt[2]+0)%3); 1710 orntNew[2] = ornt[1] < 0 ? 0 : -2; 1711 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1712 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1713 #if 1 1714 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1715 for (p = 0; p < 3; ++p) { 1716 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); 1717 } 1718 #endif 1719 supportNew[0] = (c - cStart)*8 + 0+4; 1720 supportNew[1] = (c - cStart)*8 + 2+4; 1721 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1722 #if 1 1723 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1724 for (p = 0; p < 2; ++p) { 1725 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); 1726 } 1727 #endif 1728 ++newp; 1729 /* Face G: {e, a, f} */ 1730 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+1)%3 : (ornt[1]+1)%3); 1731 orntNew[0] = ornt[1] < 0 ? -2 : 0; 1732 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1733 orntNew[1] = -2; 1734 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+1)%3 : (ornt[3]+1)%3); 1735 orntNew[2] = ornt[3] < 0 ? 0 : -2; 1736 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1737 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1738 #if 1 1739 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1740 for (p = 0; p < 3; ++p) { 1741 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); 1742 } 1743 #endif 1744 supportNew[0] = (c - cStart)*8 + 1+4; 1745 supportNew[1] = (c - cStart)*8 + 3+4; 1746 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1747 #if 1 1748 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1749 for (p = 0; p < 2; ++p) { 1750 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); 1751 } 1752 #endif 1753 ++newp; 1754 /* Face H: {a, b, f} */ 1755 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+2)%3 : (ornt[0]+0)%3); 1756 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1757 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+0)%3 : (ornt[3]+2)%3); 1758 orntNew[1] = ornt[3] < 0 ? 0 : -2; 1759 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1760 orntNew[2] = 0; 1761 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1762 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1763 #if 1 1764 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1765 for (p = 0; p < 3; ++p) { 1766 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); 1767 } 1768 #endif 1769 supportNew[0] = (c - cStart)*8 + 1+4; 1770 supportNew[1] = (c - cStart)*8 + 2+4; 1771 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1772 #if 1 1773 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1774 for (p = 0; p < 2; ++p) { 1775 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); 1776 } 1777 #endif 1778 ++newp; 1779 } 1780 /* Split Edges have 2 vertices and the same faces as the parent */ 1781 for (e = eStart; e < eEnd; ++e) { 1782 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 1783 1784 for (r = 0; r < 2; ++r) { 1785 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 1786 const PetscInt *cone, *ornt, *support; 1787 PetscInt coneNew[2], coneSize, c, supportSize, s; 1788 1789 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 1790 coneNew[0] = vStartNew + (cone[0] - vStart); 1791 coneNew[1] = vStartNew + (cone[1] - vStart); 1792 coneNew[(r+1)%2] = newv; 1793 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1794 #if 1 1795 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1796 for (p = 0; p < 2; ++p) { 1797 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); 1798 } 1799 #endif 1800 ierr = DMPlexGetSupportSize(dm, e, &supportSize);CHKERRQ(ierr); 1801 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 1802 for (s = 0; s < supportSize; ++s) { 1803 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1804 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1805 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1806 for (c = 0; c < coneSize; ++c) { 1807 if (cone[c] == e) break; 1808 } 1809 supportRef[s] = fStartNew + (support[s] - fStart)*4 + (c + (ornt[c] < 0 ? 1-r : r))%3; 1810 } 1811 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1812 #if 1 1813 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1814 for (p = 0; p < supportSize; ++p) { 1815 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); 1816 } 1817 #endif 1818 } 1819 } 1820 /* Face edges have 2 vertices and 2+cells*(1/2) faces */ 1821 for (f = fStart; f < fEnd; ++f) { 1822 const PetscInt *cone, *ornt, *support; 1823 PetscInt coneSize, supportSize, s; 1824 1825 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1826 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1827 for (r = 0; r < 3; ++r) { 1828 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + r; 1829 PetscInt coneNew[2], intFaces = 0, er, eint[4] = {1, 0, 2, 0}; 1830 PetscInt fint[24] = { 1, 7, -1, -1, 0, 5, 1831 -1, -1, 1, 6, 0, 4, 1832 2, 5, 3, 4, -1, -1, 1833 -1, -1, 3, 6, 2, 7}; 1834 1835 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1836 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[(r+0)%3] - eStart); 1837 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - eStart); 1838 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1839 #if 1 1840 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1841 for (p = 0; p < 2; ++p) { 1842 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); 1843 } 1844 #endif 1845 supportRef[0] = fStartNew + (f - fStart)*4 + (r+1)%3; 1846 supportRef[1] = fStartNew + (f - fStart)*4 + 3; 1847 for (s = 0; s < supportSize; ++s) { 1848 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1849 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1850 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1851 for (c = 0; c < coneSize; ++c) {if (cone[c] == f) break;} 1852 /* Here we want to determine whether edge newp contains a vertex which is part of the cross-tet edge */ 1853 er = ornt[c] < 0 ? (-(ornt[c]+1) + 2-r)%3 : (ornt[c] + r)%3; 1854 if (er == eint[c]) { 1855 supportRef[2+intFaces++] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*8 + (c + 2)%4; 1856 } else { 1857 supportRef[2+intFaces++] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*8 + fint[(c*3 + er)*2 + 0]; 1858 supportRef[2+intFaces++] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*8 + fint[(c*3 + er)*2 + 1]; 1859 } 1860 } 1861 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1862 #if 1 1863 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1864 for (p = 0; p < intFaces; ++p) { 1865 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); 1866 } 1867 #endif 1868 } 1869 } 1870 /* Interior edges have 2 vertices and 4 faces */ 1871 for (c = cStart; c < cEnd; ++c) { 1872 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1873 const PetscInt *cone, *ornt, *fcone; 1874 PetscInt coneNew[2], supportNew[4], find; 1875 1876 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1877 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1878 ierr = DMPlexGetCone(dm, cone[0], &fcone);CHKERRQ(ierr); 1879 find = GetTriEdge_Static(ornt[0], 0); 1880 coneNew[0] = vStartNew + (vEnd - vStart) + (fcone[find] - eStart); 1881 ierr = DMPlexGetCone(dm, cone[2], &fcone);CHKERRQ(ierr); 1882 find = GetTriEdge_Static(ornt[2], 1); 1883 coneNew[1] = vStartNew + (vEnd - vStart) + (fcone[find] - eStart); 1884 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1885 #if 1 1886 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1887 for (p = 0; p < 2; ++p) { 1888 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); 1889 } 1890 #endif 1891 supportNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 4; 1892 supportNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 5; 1893 supportNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 6; 1894 supportNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 7; 1895 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1896 #if 1 1897 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1898 for (p = 0; p < 4; ++p) { 1899 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); 1900 } 1901 #endif 1902 } 1903 /* Old vertices have identical supports */ 1904 for (v = vStart; v < vEnd; ++v) { 1905 const PetscInt newp = vStartNew + (v - vStart); 1906 const PetscInt *support, *cone; 1907 PetscInt size, s; 1908 1909 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 1910 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 1911 for (s = 0; s < size; ++s) { 1912 PetscInt r = 0; 1913 1914 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1915 if (cone[1] == v) r = 1; 1916 supportRef[s] = eStartNew + (support[s] - eStart)*2 + r; 1917 } 1918 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1919 #if 1 1920 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1921 for (p = 0; p < size; ++p) { 1922 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); 1923 } 1924 #endif 1925 } 1926 /* Edge vertices have 2 + face*2 + 0/1 supports */ 1927 for (e = eStart; e < eEnd; ++e) { 1928 const PetscInt newp = vStartNew + (vEnd - vStart) + (e - eStart); 1929 const PetscInt *cone, *support; 1930 PetscInt *star = NULL, starSize, cellSize = 0, coneSize, size, s; 1931 1932 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 1933 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 1934 supportRef[0] = eStartNew + (e - eStart)*2 + 0; 1935 supportRef[1] = eStartNew + (e - eStart)*2 + 1; 1936 for (s = 0; s < size; ++s) { 1937 PetscInt r = 0; 1938 1939 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1940 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1941 for (r = 0; r < coneSize; ++r) {if (cone[r] == e) break;} 1942 supportRef[2+s*2+0] = eStartNew + (eEnd - eStart)*2 + (support[s] - fStart)*3 + (r+0)%3; 1943 supportRef[2+s*2+1] = eStartNew + (eEnd - eStart)*2 + (support[s] - fStart)*3 + (r+2)%3; 1944 } 1945 ierr = DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1946 for (s = 0; s < starSize*2; s += 2) { 1947 const PetscInt *cone, *ornt; 1948 PetscInt e01, e23; 1949 1950 if ((star[s] >= cStart) && (star[s] < cEnd)) { 1951 /* Check edge 0-1 */ 1952 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 1953 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 1954 ierr = DMPlexGetCone(dm, cone[0], &cone);CHKERRQ(ierr); 1955 e01 = cone[GetTriEdge_Static(ornt[0], 0)]; 1956 /* Check edge 2-3 */ 1957 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 1958 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 1959 ierr = DMPlexGetCone(dm, cone[2], &cone);CHKERRQ(ierr); 1960 e23 = cone[GetTriEdge_Static(ornt[2], 1)]; 1961 if ((e01 == e) || (e23 == e)) {supportRef[2+size*2+cellSize++] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (star[s] - cStart);} 1962 } 1963 } 1964 ierr = DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1965 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1966 #if 1 1967 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1968 for (p = 0; p < 2+size*2+cellSize; ++p) { 1969 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); 1970 } 1971 #endif 1972 } 1973 ierr = PetscFree(supportRef);CHKERRQ(ierr); 1974 ierr = DMPlexRestoreFaces_Internal(dm, 3, cStart, NULL, NULL, &faces);CHKERRQ(ierr); 1975 break; 1976 case 6: 1977 /* Hex 3D */ 1978 /* 1979 Bottom (viewed from top) Top 1980 1---------2---------2 7---------2---------6 1981 | | | | | | 1982 | B 2 C | | H 2 G | 1983 | | | | | | 1984 3----3----0----1----1 3----3----0----1----1 1985 | | | | | | 1986 | A 0 D | | E 0 F | 1987 | | | | | | 1988 0---------0---------3 4---------0---------5 1989 */ 1990 /* All cells have 6 faces: Bottom, Top, Front, Back, Right, Left */ 1991 for (c = cStart; c < cEnd; ++c) { 1992 const PetscInt newp = (c - cStart)*8; 1993 const PetscInt *cone, *ornt; 1994 PetscInt coneNew[6], orntNew[6]; 1995 1996 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1997 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1998 /* A hex */ 1999 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetRefHexFace_Static(ornt[0], 0); 2000 orntNew[0] = ornt[0]; 2001 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 8; /* AE */ 2002 orntNew[1] = 0; 2003 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetRefHexFace_Static(ornt[2], 0); 2004 orntNew[2] = ornt[2]; 2005 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 3; /* AB */ 2006 orntNew[3] = 0; 2007 coneNew[4] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 0; /* AD */ 2008 orntNew[4] = 0; 2009 coneNew[5] = fStartNew + (cone[5] - fStart)*4 + GetRefHexFace_Static(ornt[5], 0); 2010 orntNew[5] = ornt[5]; 2011 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 2012 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 2013 #if 1 2014 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); 2015 for (p = 0; p < 6; ++p) { 2016 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); 2017 } 2018 #endif 2019 /* B hex */ 2020 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetRefHexFace_Static(ornt[0], 1); 2021 orntNew[0] = ornt[0]; 2022 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 11; /* BH */ 2023 orntNew[1] = 0; 2024 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 3; /* AB */ 2025 orntNew[2] = -3; 2026 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetRefHexFace_Static(ornt[3], 1); 2027 orntNew[3] = ornt[3]; 2028 coneNew[4] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 2; /* BC */ 2029 orntNew[4] = 0; 2030 coneNew[5] = fStartNew + (cone[5] - fStart)*4 + GetRefHexFace_Static(ornt[5], 3); 2031 orntNew[5] = ornt[5]; 2032 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 2033 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 2034 #if 1 2035 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); 2036 for (p = 0; p < 6; ++p) { 2037 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); 2038 } 2039 #endif 2040 /* C hex */ 2041 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetRefHexFace_Static(ornt[0], 2); 2042 orntNew[0] = ornt[0]; 2043 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 10; /* CG */ 2044 orntNew[1] = 0; 2045 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 1; /* CD */ 2046 orntNew[2] = 0; 2047 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetRefHexFace_Static(ornt[3], 1); 2048 orntNew[3] = ornt[3]; 2049 coneNew[4] = fStartNew + (cone[4] - fStart)*4 + GetRefHexFace_Static(ornt[4], 1); 2050 orntNew[4] = ornt[4]; 2051 coneNew[5] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 2; /* BC */ 2052 orntNew[5] = -3; 2053 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 2054 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 2055 #if 1 2056 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); 2057 for (p = 0; p < 6; ++p) { 2058 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); 2059 } 2060 #endif 2061 /* D hex */ 2062 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + GetRefHexFace_Static(ornt[0], 3); 2063 orntNew[0] = ornt[0]; 2064 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 9; /* DF */ 2065 orntNew[1] = 0; 2066 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetRefHexFace_Static(ornt[2], 1); 2067 orntNew[2] = ornt[2]; 2068 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 1; /* CD */ 2069 orntNew[3] = -3; 2070 coneNew[4] = fStartNew + (cone[4] - fStart)*4 + GetRefHexFace_Static(ornt[4], 0); 2071 orntNew[4] = ornt[4]; 2072 coneNew[5] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 0; /* AD */ 2073 orntNew[5] = -3; 2074 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 2075 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 2076 #if 1 2077 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); 2078 for (p = 0; p < 6; ++p) { 2079 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); 2080 } 2081 #endif 2082 /* E hex */ 2083 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 8; /* AE */ 2084 orntNew[0] = -3; 2085 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetRefHexFace_Static(ornt[1], 0); 2086 orntNew[1] = ornt[1]; 2087 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetRefHexFace_Static(ornt[2], 3); 2088 orntNew[2] = ornt[2]; 2089 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 7; /* EH */ 2090 orntNew[3] = 0; 2091 coneNew[4] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 4; /* EF */ 2092 orntNew[4] = 0; 2093 coneNew[5] = fStartNew + (cone[5] - fStart)*4 + GetRefHexFace_Static(ornt[5], 1); 2094 orntNew[5] = ornt[5]; 2095 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 2096 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 2097 #if 1 2098 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); 2099 for (p = 0; p < 6; ++p) { 2100 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); 2101 } 2102 #endif 2103 /* F hex */ 2104 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 9; /* DF */ 2105 orntNew[0] = -3; 2106 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetRefHexFace_Static(ornt[1], 1); 2107 orntNew[1] = ornt[1]; 2108 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + GetRefHexFace_Static(ornt[2], 2); 2109 orntNew[2] = ornt[2]; 2110 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 5; /* FG */ 2111 orntNew[3] = 0; 2112 coneNew[4] = fStartNew + (cone[4] - fStart)*4 + GetRefHexFace_Static(ornt[4], 3); 2113 orntNew[4] = ornt[4]; 2114 coneNew[5] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 4; /* EF */ 2115 orntNew[5] = -3; 2116 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 2117 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 2118 #if 1 2119 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); 2120 for (p = 0; p < 6; ++p) { 2121 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); 2122 } 2123 #endif 2124 /* G hex */ 2125 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 10; /* CG */ 2126 orntNew[0] = -3; 2127 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetRefHexFace_Static(ornt[1], 2); 2128 orntNew[1] = ornt[1]; 2129 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 5; /* FG */ 2130 orntNew[2] = -3; 2131 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetRefHexFace_Static(ornt[3], 3); 2132 orntNew[3] = ornt[3]; 2133 coneNew[4] = fStartNew + (cone[4] - fStart)*4 + GetRefHexFace_Static(ornt[4], 2); 2134 orntNew[4] = ornt[4]; 2135 coneNew[5] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 6; /* GH */ 2136 orntNew[5] = 0; 2137 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 2138 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 2139 #if 1 2140 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); 2141 for (p = 0; p < 6; ++p) { 2142 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); 2143 } 2144 #endif 2145 /* H hex */ 2146 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 11; /* BH */ 2147 orntNew[0] = -3; 2148 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + GetRefHexFace_Static(ornt[1], 3); 2149 orntNew[1] = ornt[1]; 2150 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 7; /* EH */ 2151 orntNew[2] = -3; 2152 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + GetRefHexFace_Static(ornt[3], 2); 2153 orntNew[3] = ornt[3]; 2154 coneNew[4] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + 6; /* GH */ 2155 orntNew[4] = -3; 2156 coneNew[5] = fStartNew + (cone[5] - fStart)*4 + GetRefHexFace_Static(ornt[5], 2); 2157 orntNew[5] = ornt[5]; 2158 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 2159 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 2160 #if 1 2161 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); 2162 for (p = 0; p < 6; ++p) { 2163 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); 2164 } 2165 #endif 2166 } 2167 /* Split faces have 4 edges and the same cells as the parent */ 2168 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 2169 ierr = PetscMalloc((4 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 2170 for (f = fStart; f < fEnd; ++f) { 2171 for (r = 0; r < 4; ++r) { 2172 /* This can come from GetFaces_Internal() */ 2173 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}; 2174 const PetscInt newp = fStartNew + (f - fStart)*4 + r; 2175 const PetscInt *cone, *ornt, *support; 2176 PetscInt coneNew[4], coneSize, c, supportSize, s; 2177 2178 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 2179 /* TODO: Redo using orientation information */ 2180 coneNew[0] = eStartNew + (cone[(r+3)%4] - eStart)*2 + 1; 2181 coneNew[1] = eStartNew + (cone[r] - eStart)*2 + 0; 2182 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + r; 2183 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + (r+3)%4; 2184 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2185 #if 1 2186 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2187 for (p = 0; p < 4; ++p) { 2188 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); 2189 } 2190 #endif 2191 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 2192 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 2193 for (s = 0; s < supportSize; ++s) { 2194 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 2195 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2196 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 2197 for (c = 0; c < coneSize; ++c) { 2198 if (cone[c] == f) break; 2199 } 2200 /* TODO: Redo using orientation information */ 2201 supportRef[s] = cStartNew + (support[s] - cStart)*8 + newCells[c*4+r]; 2202 } 2203 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2204 #if 1 2205 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2206 for (p = 0; p < supportSize; ++p) { 2207 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); 2208 } 2209 #endif 2210 } 2211 } 2212 /* Interior faces have 4 edges and 2 cells */ 2213 for (c = cStart; c < cEnd; ++c) { 2214 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}; 2215 const PetscInt *cone; 2216 PetscInt coneNew[4], supportNew[2]; 2217 2218 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2219 for (r = 0; r < 12; ++r) { 2220 const PetscInt newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + r; 2221 2222 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[r] - fStart); 2223 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart) + (c - cStart); 2224 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[r] - fStart); 2225 coneNew[3] = eStartNew + (eEnd - eStart)*2 + (cone[r] - fStart); 2226 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2227 #if 1 2228 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2229 for (p = 0; p < 4; ++p) { 2230 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); 2231 } 2232 #endif 2233 supportNew[0] = cStartNew + (c - cStart)*8 + newCells[r*2+0]; 2234 supportNew[1] = cStartNew + (c - cStart)*8 + newCells[r*2+1]; 2235 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 2236 #if 1 2237 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 2238 for (p = 0; p < 2; ++p) { 2239 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); 2240 } 2241 #endif 2242 } 2243 } 2244 /* Split edges have 2 vertices and the same faces as the parent */ 2245 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 2246 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 2247 for (e = eStart; e < eEnd; ++e) { 2248 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 2249 2250 for (r = 0; r < 2; ++r) { 2251 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 2252 const PetscInt *cone, *ornt, *support; 2253 PetscInt coneNew[2], coneSize, c, supportSize, s; 2254 2255 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 2256 coneNew[0] = vStartNew + (cone[0] - vStart); 2257 coneNew[1] = vStartNew + (cone[1] - vStart); 2258 coneNew[(r+1)%2] = newv; 2259 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2260 #if 1 2261 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2262 for (p = 0; p < 2; ++p) { 2263 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); 2264 } 2265 #endif 2266 ierr = DMPlexGetSupportSize(dm, e, &supportSize);CHKERRQ(ierr); 2267 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 2268 for (s = 0; s < supportSize; ++s) { 2269 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 2270 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2271 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 2272 for (c = 0; c < coneSize; ++c) { 2273 if (cone[c] == e) break; 2274 } 2275 supportRef[s] = fStartNew + (support[s] - fStart)*4 + (ornt[c] < 0 ? (c+1-r)%4 : (c+r)%4); 2276 } 2277 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2278 #if 1 2279 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2280 for (p = 0; p < supportSize; ++p) { 2281 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); 2282 } 2283 #endif 2284 } 2285 } 2286 /* Face edges have 2 vertices and 2+cells faces */ 2287 for (f = fStart; f < fEnd; ++f) { 2288 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}; 2289 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 2290 const PetscInt *cone, *coneCell, *support; 2291 PetscInt coneNew[2], coneSize, c, supportSize, s; 2292 2293 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 2294 for (r = 0; r < 4; ++r) { 2295 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (f - fStart)*4 + r; 2296 2297 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - eStart); 2298 coneNew[1] = newv; 2299 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2300 #if 1 2301 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2302 for (p = 0; p < 2; ++p) { 2303 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); 2304 } 2305 #endif 2306 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 2307 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 2308 supportRef[0] = fStartNew + (f - fStart)*4 + r; 2309 supportRef[1] = fStartNew + (f - fStart)*4 + (r+1)%4; 2310 for (s = 0; s < supportSize; ++s) { 2311 ierr = DMPlexGetConeSize(dm, f, &coneSize);CHKERRQ(ierr); 2312 ierr = DMPlexGetCone(dm, f, &coneCell);CHKERRQ(ierr); 2313 for (c = 0; c < coneSize; ++c) if (coneCell[c] == f) break; 2314 supportRef[2+s] = fStartNew + (f - fStart)*4 + newFaces[c*4+r]; 2315 } 2316 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2317 #if 1 2318 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2319 for (p = 0; p < 2+supportSize; ++p) { 2320 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); 2321 } 2322 #endif 2323 } 2324 } 2325 /* Cell edges have 2 vertices and 4 faces */ 2326 for (c = cStart; c < cEnd; ++c) { 2327 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}; 2328 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (c - cStart); 2329 const PetscInt *cone; 2330 PetscInt coneNew[2], supportNew[4]; 2331 2332 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2333 for (r = 0; r < 6; ++r) { 2334 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + r; 2335 2336 coneNew[0] = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (cone[r] - fStart); 2337 coneNew[1] = newv; 2338 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 2339 #if 1 2340 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2341 for (p = 0; p < 2; ++p) { 2342 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); 2343 } 2344 #endif 2345 for (f = 0; f < 4; ++f) supportNew[f] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*12 + newFaces[r*4+f]; 2346 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 2347 #if 1 2348 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 2349 for (p = 0; p < 4; ++p) { 2350 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); 2351 } 2352 #endif 2353 } 2354 } 2355 /* Old vertices have identical supports */ 2356 for (v = vStart; v < vEnd; ++v) { 2357 const PetscInt newp = vStartNew + (v - vStart); 2358 const PetscInt *support, *cone; 2359 PetscInt size, s; 2360 2361 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 2362 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 2363 for (s = 0; s < size; ++s) { 2364 PetscInt r = 0; 2365 2366 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2367 if (cone[1] == v) r = 1; 2368 supportRef[s] = eStartNew + (support[s] - eStart)*2 + r; 2369 } 2370 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2371 #if 1 2372 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2373 for (p = 0; p < size; ++p) { 2374 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); 2375 } 2376 #endif 2377 } 2378 /* Edge vertices have 2 + faces supports */ 2379 for (e = eStart; e < eEnd; ++e) { 2380 const PetscInt newp = vStartNew + (vEnd - vStart) + (e - eStart); 2381 const PetscInt *cone, *support; 2382 PetscInt size, s; 2383 2384 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 2385 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 2386 supportRef[0] = eStartNew + (e - eStart)*2 + 0; 2387 supportRef[1] = eStartNew + (e - eStart)*2 + 1; 2388 for (s = 0; s < size; ++s) { 2389 PetscInt r; 2390 2391 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2392 for (r = 0; r < 4; ++r) if (cone[r] == f) break; 2393 supportRef[2+s] = eStartNew + (eEnd - eStart)*2 + (support[s] - fStart)*4 + r; 2394 } 2395 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2396 #if 1 2397 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2398 for (p = 0; p < 2+size; ++p) { 2399 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); 2400 } 2401 #endif 2402 } 2403 /* Face vertices have 4 + cells supports */ 2404 for (f = fStart; f < fEnd; ++f) { 2405 const PetscInt newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 2406 const PetscInt *cone, *support; 2407 PetscInt size, s; 2408 2409 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 2410 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 2411 for (r = 0; r < 4; ++r) supportRef[r] = eStartNew + (e - eStart)*2 + (f - fStart)*4 + r; 2412 for (s = 0; s < size; ++s) { 2413 PetscInt r; 2414 2415 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 2416 for (r = 0; r < 6; ++r) if (cone[r] == f) break; 2417 supportRef[4+s] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (support[s] - cStart)*6 + r; 2418 } 2419 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 2420 #if 1 2421 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 2422 for (p = 0; p < 4+size; ++p) { 2423 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); 2424 } 2425 #endif 2426 } 2427 /* Cell vertices have 6 supports */ 2428 for (c = cStart; c < cEnd; ++c) { 2429 const PetscInt newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (c - cStart); 2430 PetscInt supportNew[6]; 2431 2432 for (r = 0; r < 6; ++r) { 2433 supportNew[r] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (c - cStart)*6 + r; 2434 } 2435 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 2436 } 2437 ierr = PetscFree(supportRef);CHKERRQ(ierr); 2438 break; 2439 default: 2440 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2441 } 2442 PetscFunctionReturn(0); 2443 } 2444 2445 #undef __FUNCT__ 2446 #define __FUNCT__ "CellRefinerSetCoordinates" 2447 PetscErrorCode CellRefinerSetCoordinates(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 2448 { 2449 PetscSection coordSection, coordSectionNew; 2450 Vec coordinates, coordinatesNew; 2451 PetscScalar *coords, *coordsNew; 2452 PetscInt dim, depth, coordSizeNew, cStart, cEnd, c, vStart, vStartNew, vEnd, v, eStart, eEnd, eMax, e, fStart, fEnd, fMax, f; 2453 PetscErrorCode ierr; 2454 2455 PetscFunctionBegin; 2456 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2457 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2458 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2459 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 2460 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2461 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2462 ierr = DMPlexGetHybridBounds(dm, NULL, &fMax, &eMax, NULL);CHKERRQ(ierr); 2463 ierr = GetDepthStart_Private(depth, depthSize, NULL, NULL, NULL, &vStartNew);CHKERRQ(ierr); 2464 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2465 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &coordSectionNew);CHKERRQ(ierr); 2466 ierr = PetscSectionSetNumFields(coordSectionNew, 1);CHKERRQ(ierr); 2467 ierr = PetscSectionSetFieldComponents(coordSectionNew, 0, dim);CHKERRQ(ierr); 2468 ierr = PetscSectionSetChart(coordSectionNew, vStartNew, vStartNew+depthSize[0]);CHKERRQ(ierr); 2469 if (fMax < 0) fMax = fEnd; 2470 if (eMax < 0) eMax = eEnd; 2471 /* All vertices have the dim coordinates */ 2472 for (v = vStartNew; v < vStartNew+depthSize[0]; ++v) { 2473 ierr = PetscSectionSetDof(coordSectionNew, v, dim);CHKERRQ(ierr); 2474 ierr = PetscSectionSetFieldDof(coordSectionNew, v, 0, dim);CHKERRQ(ierr); 2475 } 2476 ierr = PetscSectionSetUp(coordSectionNew);CHKERRQ(ierr); 2477 ierr = DMPlexSetCoordinateSection(rdm, coordSectionNew);CHKERRQ(ierr); 2478 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2479 ierr = PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);CHKERRQ(ierr); 2480 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinatesNew);CHKERRQ(ierr); 2481 ierr = PetscObjectSetName((PetscObject) coordinatesNew, "coordinates");CHKERRQ(ierr); 2482 ierr = VecSetSizes(coordinatesNew, coordSizeNew, PETSC_DETERMINE);CHKERRQ(ierr); 2483 ierr = VecSetFromOptions(coordinatesNew);CHKERRQ(ierr); 2484 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2485 ierr = VecGetArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 2486 switch (refiner) { 2487 case 6: /* Hex 3D */ 2488 /* Face vertices have the average of corner coordinates */ 2489 for (f = fStart; f < fEnd; ++f) { 2490 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (f - fStart); 2491 PetscInt *cone = NULL; 2492 PetscInt closureSize, coneSize = 0, off[8], offnew, p, d; 2493 2494 ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2495 for (p = 0; p < closureSize*2; p += 2) { 2496 const PetscInt point = cone[p]; 2497 if ((point >= vStart) && (point < vEnd)) cone[coneSize++] = point; 2498 } 2499 for (v = 0; v < coneSize; ++v) { 2500 ierr = PetscSectionGetOffset(coordSection, cone[v], &off[v]);CHKERRQ(ierr); 2501 } 2502 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2503 for (d = 0; d < dim; ++d) { 2504 coordsNew[offnew+d] = 0.0; 2505 for (v = 0; v < coneSize; ++v) coordsNew[offnew+d] += coords[off[v]+d]; 2506 coordsNew[offnew+d] /= coneSize; 2507 } 2508 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2509 } 2510 case 2: /* Hex 2D */ 2511 /* Cell vertices have the average of corner coordinates */ 2512 for (c = cStart; c < cEnd; ++c) { 2513 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (c - cStart) + (dim > 2 ? (fEnd - fStart) : 0); 2514 PetscInt *cone = NULL; 2515 PetscInt closureSize, coneSize = 0, off[8], offnew, p, d; 2516 2517 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2518 for (p = 0; p < closureSize*2; p += 2) { 2519 const PetscInt point = cone[p]; 2520 if ((point >= vStart) && (point < vEnd)) cone[coneSize++] = point; 2521 } 2522 for (v = 0; v < coneSize; ++v) { 2523 ierr = PetscSectionGetOffset(coordSection, cone[v], &off[v]);CHKERRQ(ierr); 2524 } 2525 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2526 for (d = 0; d < dim; ++d) { 2527 coordsNew[offnew+d] = 0.0; 2528 for (v = 0; v < coneSize; ++v) coordsNew[offnew+d] += coords[off[v]+d]; 2529 coordsNew[offnew+d] /= coneSize; 2530 } 2531 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 2532 } 2533 case 1: /* Simplicial 2D */ 2534 case 3: /* Hybrid Simplicial 2D */ 2535 case 5: /* Simplicial 3D */ 2536 /* Edge vertices have the average of endpoint coordinates */ 2537 for (e = eStart; e < eMax; ++e) { 2538 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 2539 const PetscInt *cone; 2540 PetscInt coneSize, offA, offB, offnew, d; 2541 2542 ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr); 2543 if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %d cone should have two vertices, not %d", e, coneSize); 2544 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 2545 ierr = PetscSectionGetOffset(coordSection, cone[0], &offA);CHKERRQ(ierr); 2546 ierr = PetscSectionGetOffset(coordSection, cone[1], &offB);CHKERRQ(ierr); 2547 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2548 for (d = 0; d < dim; ++d) { 2549 coordsNew[offnew+d] = 0.5*(coords[offA+d] + coords[offB+d]); 2550 } 2551 } 2552 /* Old vertices have the same coordinates */ 2553 for (v = vStart; v < vEnd; ++v) { 2554 const PetscInt newv = vStartNew + (v - vStart); 2555 PetscInt off, offnew, d; 2556 2557 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 2558 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 2559 for (d = 0; d < dim; ++d) { 2560 coordsNew[offnew+d] = coords[off+d]; 2561 } 2562 } 2563 break; 2564 default: 2565 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2566 } 2567 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2568 ierr = VecRestoreArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 2569 ierr = DMSetCoordinatesLocal(rdm, coordinatesNew);CHKERRQ(ierr); 2570 ierr = VecDestroy(&coordinatesNew);CHKERRQ(ierr); 2571 ierr = PetscSectionDestroy(&coordSectionNew);CHKERRQ(ierr); 2572 PetscFunctionReturn(0); 2573 } 2574 2575 #undef __FUNCT__ 2576 #define __FUNCT__ "DMPlexCreateProcessSF" 2577 PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess) 2578 { 2579 PetscInt numRoots, numLeaves, l; 2580 const PetscInt *localPoints; 2581 const PetscSFNode *remotePoints; 2582 PetscInt *localPointsNew; 2583 PetscSFNode *remotePointsNew; 2584 PetscInt *ranks, *ranksNew; 2585 PetscErrorCode ierr; 2586 2587 PetscFunctionBegin; 2588 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2589 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranks);CHKERRQ(ierr); 2590 for (l = 0; l < numLeaves; ++l) { 2591 ranks[l] = remotePoints[l].rank; 2592 } 2593 ierr = PetscSortRemoveDupsInt(&numLeaves, ranks);CHKERRQ(ierr); 2594 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranksNew);CHKERRQ(ierr); 2595 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 2596 ierr = PetscMalloc(numLeaves * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 2597 for (l = 0; l < numLeaves; ++l) { 2598 ranksNew[l] = ranks[l]; 2599 localPointsNew[l] = l; 2600 remotePointsNew[l].index = 0; 2601 remotePointsNew[l].rank = ranksNew[l]; 2602 } 2603 ierr = PetscFree(ranks);CHKERRQ(ierr); 2604 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);CHKERRQ(ierr); 2605 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);CHKERRQ(ierr); 2606 ierr = PetscSFSetFromOptions(*sfProcess);CHKERRQ(ierr); 2607 ierr = PetscSFSetGraph(*sfProcess, 1, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2608 PetscFunctionReturn(0); 2609 } 2610 2611 #undef __FUNCT__ 2612 #define __FUNCT__ "CellRefinerCreateSF" 2613 PetscErrorCode CellRefinerCreateSF(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 2614 { 2615 PetscSF sf, sfNew, sfProcess; 2616 IS processRanks; 2617 MPI_Datatype depthType; 2618 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 2619 const PetscInt *localPoints, *neighbors; 2620 const PetscSFNode *remotePoints; 2621 PetscInt *localPointsNew; 2622 PetscSFNode *remotePointsNew; 2623 PetscInt *depthSizeOld, *rdepthSize, *rdepthSizeOld, *rdepthMaxOld, *rvStart, *rvStartNew, *reStart, *reStartNew, *rfStart, *rfStartNew, *rcStart, *rcStartNew; 2624 PetscInt depth, numNeighbors, pStartNew, pEndNew, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eStartNew, eEnd, eMax, r, n; 2625 PetscErrorCode ierr; 2626 2627 PetscFunctionBegin; 2628 ierr = DMPlexGetChart(rdm, &pStartNew, &pEndNew);CHKERRQ(ierr); 2629 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2630 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2631 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 2632 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2633 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2634 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 2635 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 2636 switch (refiner) { 2637 case 3: 2638 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 2639 cMax = PetscMin(cEnd, cMax); 2640 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 2641 fMax = PetscMin(fEnd, fMax); 2642 } 2643 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2644 ierr = DMGetPointSF(rdm, &sfNew);CHKERRQ(ierr); 2645 /* Caculate size of new SF */ 2646 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2647 if (numRoots < 0) PetscFunctionReturn(0); 2648 for (l = 0; l < numLeaves; ++l) { 2649 const PetscInt p = localPoints[l]; 2650 2651 switch (refiner) { 2652 case 1: 2653 /* Simplicial 2D */ 2654 if ((p >= vStart) && (p < vEnd)) { 2655 /* Old vertices stay the same */ 2656 ++numLeavesNew; 2657 } else if ((p >= fStart) && (p < fEnd)) { 2658 /* Old faces add new faces and vertex */ 2659 numLeavesNew += 2 + 1; 2660 } else if ((p >= cStart) && (p < cEnd)) { 2661 /* Old cells add new cells and interior faces */ 2662 numLeavesNew += 4 + 3; 2663 } 2664 break; 2665 case 2: 2666 /* Hex 2D */ 2667 if ((p >= vStart) && (p < vEnd)) { 2668 /* Old vertices stay the same */ 2669 ++numLeavesNew; 2670 } else if ((p >= fStart) && (p < fEnd)) { 2671 /* Old faces add new faces and vertex */ 2672 numLeavesNew += 2 + 1; 2673 } else if ((p >= cStart) && (p < cEnd)) { 2674 /* Old cells add new cells, interior faces, and vertex */ 2675 numLeavesNew += 4 + 4 + 1; 2676 } 2677 break; 2678 case 5: 2679 /* Simplicial 3D */ 2680 if ((p >= vStart) && (p < vEnd)) { 2681 /* Old vertices stay the same */ 2682 ++numLeavesNew; 2683 } else if ((p >= eStart) && (p < eEnd)) { 2684 /* Old edges add new edges and vertex */ 2685 numLeavesNew += 2 + 1; 2686 } else if ((p >= fStart) && (p < fEnd)) { 2687 /* Old faces add new faces and face edges */ 2688 numLeavesNew += 4 + 3; 2689 } else if ((p >= cStart) && (p < cEnd)) { 2690 /* Old cells add new cells and interior faces and edges */ 2691 numLeavesNew += 8 + 8 + 1; 2692 } 2693 break; 2694 case 6: 2695 /* Hex 3D */ 2696 if ((p >= vStart) && (p < vEnd)) { 2697 /* Old vertices stay the same */ 2698 ++numLeavesNew; 2699 } else if ((p >= eStart) && (p < eEnd)) { 2700 /* Old edges add new edges, and vertex */ 2701 numLeavesNew += 2 + 1; 2702 } else if ((p >= fStart) && (p < fEnd)) { 2703 /* Old faces add new faces, edges, and vertex */ 2704 numLeavesNew += 4 + 4 + 1; 2705 } else if ((p >= cStart) && (p < cEnd)) { 2706 /* Old cells add new cells, faces, edges, and vertex */ 2707 numLeavesNew += 8 + 12 + 6 + 1; 2708 } 2709 break; 2710 default: 2711 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2712 } 2713 } 2714 /* Communicate depthSizes for each remote rank */ 2715 ierr = DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);CHKERRQ(ierr); 2716 ierr = ISGetLocalSize(processRanks, &numNeighbors);CHKERRQ(ierr); 2717 ierr = PetscMalloc5((depth+1)*numNeighbors,PetscInt,&rdepthSize,numNeighbors,PetscInt,&rvStartNew,numNeighbors,PetscInt,&reStartNew,numNeighbors,PetscInt,&rfStartNew,numNeighbors,PetscInt,&rcStartNew);CHKERRQ(ierr); 2718 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); 2719 ierr = MPI_Type_contiguous(depth+1, MPIU_INT, &depthType);CHKERRQ(ierr); 2720 ierr = MPI_Type_commit(&depthType);CHKERRQ(ierr); 2721 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 2722 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 2723 for (n = 0; n < numNeighbors; ++n) { 2724 ierr = GetDepthStart_Private(depth, &rdepthSize[n*(depth+1)], &rcStartNew[n], &rfStartNew[n], &reStartNew[n], &rvStartNew[n]);CHKERRQ(ierr); 2725 } 2726 depthSizeOld[depth] = cMax; 2727 depthSizeOld[0] = vMax; 2728 depthSizeOld[depth-1] = fMax; 2729 depthSizeOld[1] = eMax; 2730 2731 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 2732 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 2733 2734 depthSizeOld[depth] = cEnd - cStart; 2735 depthSizeOld[0] = vEnd - vStart; 2736 depthSizeOld[depth-1] = fEnd - fStart; 2737 depthSizeOld[1] = eEnd - eStart; 2738 2739 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 2740 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 2741 for (n = 0; n < numNeighbors; ++n) { 2742 ierr = GetDepthStart_Private(depth, &rdepthSizeOld[n*(depth+1)], &rcStart[n], &rfStart[n], &reStart[n], &rvStart[n]);CHKERRQ(ierr); 2743 } 2744 ierr = MPI_Type_free(&depthType);CHKERRQ(ierr); 2745 ierr = PetscSFDestroy(&sfProcess);CHKERRQ(ierr); 2746 /* Calculate new point SF */ 2747 ierr = PetscMalloc(numLeavesNew * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 2748 ierr = PetscMalloc(numLeavesNew * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 2749 ierr = ISGetIndices(processRanks, &neighbors);CHKERRQ(ierr); 2750 for (l = 0, m = 0; l < numLeaves; ++l) { 2751 PetscInt p = localPoints[l]; 2752 PetscInt rp = remotePoints[l].index, n; 2753 PetscMPIInt rrank = remotePoints[l].rank; 2754 2755 ierr = PetscFindInt(rrank, numNeighbors, neighbors, &n);CHKERRQ(ierr); 2756 if (n < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %d", rrank); 2757 switch (refiner) { 2758 case 1: 2759 /* Simplicial 2D */ 2760 if ((p >= vStart) && (p < vEnd)) { 2761 /* Old vertices stay the same */ 2762 localPointsNew[m] = vStartNew + (p - vStart); 2763 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2764 remotePointsNew[m].rank = rrank; 2765 ++m; 2766 } else if ((p >= fStart) && (p < fEnd)) { 2767 /* Old faces add new faces and vertex */ 2768 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 2769 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 2770 remotePointsNew[m].rank = rrank; 2771 ++m; 2772 for (r = 0; r < 2; ++r, ++m) { 2773 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 2774 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 2775 remotePointsNew[m].rank = rrank; 2776 } 2777 } else if ((p >= cStart) && (p < cEnd)) { 2778 /* Old cells add new cells and interior faces */ 2779 for (r = 0; r < 4; ++r, ++m) { 2780 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2781 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2782 remotePointsNew[m].rank = rrank; 2783 } 2784 for (r = 0; r < 3; ++r, ++m) { 2785 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 2786 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*3 + r; 2787 remotePointsNew[m].rank = rrank; 2788 } 2789 } 2790 break; 2791 case 2: 2792 /* Hex 2D */ 2793 if ((p >= vStart) && (p < vEnd)) { 2794 /* Old vertices stay the same */ 2795 localPointsNew[m] = vStartNew + (p - vStart); 2796 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2797 remotePointsNew[m].rank = rrank; 2798 ++m; 2799 } else if ((p >= fStart) && (p < fEnd)) { 2800 /* Old faces add new faces and vertex */ 2801 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 2802 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 2803 remotePointsNew[m].rank = rrank; 2804 ++m; 2805 for (r = 0; r < 2; ++r, ++m) { 2806 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 2807 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 2808 remotePointsNew[m].rank = rrank; 2809 } 2810 } else if ((p >= cStart) && (p < cEnd)) { 2811 /* Old cells add new cells, interior faces, and vertex */ 2812 for (r = 0; r < 4; ++r, ++m) { 2813 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2814 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2815 remotePointsNew[m].rank = rrank; 2816 } 2817 for (r = 0; r < 4; ++r, ++m) { 2818 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 2819 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*4 + r; 2820 remotePointsNew[m].rank = rrank; 2821 } 2822 for (r = 0; r < 1; ++r, ++m) { 2823 localPointsNew[m] = vStartNew + (fEnd - fStart) + (p - cStart) + r; 2824 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1] + (rp - rcStart[n]) + r; 2825 remotePointsNew[m].rank = rrank; 2826 } 2827 } 2828 break; 2829 case 3: 2830 /* Hybrid simplicial 2D */ 2831 if ((p >= vStart) && (p < vEnd)) { 2832 /* Old vertices stay the same */ 2833 localPointsNew[m] = vStartNew + (p - vStart); 2834 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2835 remotePointsNew[m].rank = rrank; 2836 ++m; 2837 } else if ((p >= fStart) && (p < fMax)) { 2838 /* Old interior faces add new faces and vertex */ 2839 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 2840 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 2841 remotePointsNew[m].rank = rrank; 2842 ++m; 2843 for (r = 0; r < 2; ++r, ++m) { 2844 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 2845 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 2846 remotePointsNew[m].rank = rrank; 2847 } 2848 } else if ((p >= fMax) && (p < fEnd)) { 2849 /* Old hybrid faces stay the same */ 2850 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - fMax); 2851 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rdepthMaxOld[n*(depth+1)+depth-1]); 2852 remotePointsNew[m].rank = rrank; 2853 ++m; 2854 } else if ((p >= cStart) && (p < cMax)) { 2855 /* Old interior cells add new cells and interior faces */ 2856 for (r = 0; r < 4; ++r, ++m) { 2857 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2858 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2859 remotePointsNew[m].rank = rrank; 2860 } 2861 for (r = 0; r < 3; ++r, ++m) { 2862 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - cStart)*3 + r; 2863 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rcStart[n])*3 + r; 2864 remotePointsNew[m].rank = rrank; 2865 } 2866 } else if ((p >= cStart) && (p < cMax)) { 2867 /* Old hybrid cells add new cells and hybrid face */ 2868 for (r = 0; r < 2; ++r, ++m) { 2869 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2870 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2871 remotePointsNew[m].rank = rrank; 2872 } 2873 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 2874 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]); 2875 remotePointsNew[m].rank = rrank; 2876 ++m; 2877 } 2878 break; 2879 case 5: 2880 /* Simplicial 3D */ 2881 if ((p >= vStart) && (p < vEnd)) { 2882 /* Old vertices stay the same */ 2883 localPointsNew[m] = vStartNew + (p - vStart); 2884 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2885 remotePointsNew[m].rank = rrank; 2886 ++m; 2887 } else if ((p >= eStart) && (p < eEnd)) { 2888 /* Old edges add new edges and vertex */ 2889 for (r = 0; r < 2; ++r, ++m) { 2890 localPointsNew[m] = eStartNew + (p - eStart)*2 + r; 2891 remotePointsNew[m].index = reStartNew[n] + (rp - reStart[n])*2 + r; 2892 remotePointsNew[m].rank = rrank; 2893 } 2894 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - eStart); 2895 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - reStart[n]); 2896 remotePointsNew[m].rank = rrank; 2897 ++m; 2898 } else if ((p >= fStart) && (p < fEnd)) { 2899 /* Old faces add new faces and face edges */ 2900 for (r = 0; r < 4; ++r, ++m) { 2901 localPointsNew[m] = fStartNew + (p - fStart)*4 + r; 2902 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*4 + r; 2903 remotePointsNew[m].rank = rrank; 2904 } 2905 for (r = 0; r < 3; ++r, ++m) { 2906 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (p - fStart)*3 + r; 2907 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + (rp - rfStart[n])*3 + r; 2908 remotePointsNew[m].rank = rrank; 2909 } 2910 } else if ((p >= cStart) && (p < cEnd)) { 2911 /* Old cells add new cells and interior faces and edges */ 2912 for (r = 0; r < 8; ++r, ++m) { 2913 localPointsNew[m] = cStartNew + (p - cStart)*8 + r; 2914 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*8 + r; 2915 remotePointsNew[m].rank = rrank; 2916 } 2917 for (r = 0; r < 8; ++r, ++m) { 2918 localPointsNew[m] = fStartNew + (fEnd - fStart)*4 + (p - cStart)*8 + r; 2919 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*8 + r; 2920 remotePointsNew[m].rank = rrank; 2921 } 2922 for (r = 0; r < 1; ++r, ++m) { 2923 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (p - cStart)*1 + r; 2924 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + rdepthSizeOld[n*(depth+1)+depth-1]*3 + (rp - rcStart[n])*1 + r; 2925 remotePointsNew[m].rank = rrank; 2926 } 2927 } 2928 break; 2929 case 6: 2930 /* Hex 3D */ 2931 if ((p >= vStart) && (p < vEnd)) { 2932 /* Old vertices stay the same */ 2933 localPointsNew[m] = vStartNew + (p - vStart); 2934 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2935 remotePointsNew[m].rank = rrank; 2936 ++m; 2937 } else if ((p >= eStart) && (p < eEnd)) { 2938 /* Old edges add new edges and vertex */ 2939 for (r = 0; r < 2; ++r, ++m) { 2940 localPointsNew[m] = eStartNew + (p - eStart)*2 + r; 2941 remotePointsNew[m].index = reStartNew[n] + (rp - reStart[n])*2 + r; 2942 remotePointsNew[m].rank = rrank; 2943 } 2944 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - eStart); 2945 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - reStart[n]); 2946 remotePointsNew[m].rank = rrank; 2947 ++m; 2948 } else if ((p >= fStart) && (p < fEnd)) { 2949 /* Old faces add new faces, edges, and vertex */ 2950 for (r = 0; r < 4; ++r, ++m) { 2951 localPointsNew[m] = fStartNew + (p - fStart)*4 + r; 2952 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*4 + r; 2953 remotePointsNew[m].rank = rrank; 2954 } 2955 for (r = 0; r < 4; ++r, ++m) { 2956 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (p - fStart)*4 + r; 2957 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + (rp - rfStart[n])*4 + r; 2958 remotePointsNew[m].rank = rrank; 2959 } 2960 localPointsNew[m] = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (p - fStart); 2961 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + rdepthSizeOld[n*(depth+1)+1] + (rp - rfStart[n]); 2962 remotePointsNew[m].rank = rrank; 2963 ++m; 2964 } else if ((p >= cStart) && (p < cEnd)) { 2965 /* Old cells add new cells, faces, edges, and vertex */ 2966 for (r = 0; r < 8; ++r, ++m) { 2967 localPointsNew[m] = cStartNew + (p - cStart)*8 + r; 2968 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*8 + r; 2969 remotePointsNew[m].rank = rrank; 2970 } 2971 for (r = 0; r < 12; ++r, ++m) { 2972 localPointsNew[m] = fStartNew + (fEnd - fStart)*4 + (p - cStart)*12 + r; 2973 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*12 + r; 2974 remotePointsNew[m].rank = rrank; 2975 } 2976 for (r = 0; r < 6; ++r, ++m) { 2977 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (p - cStart)*6 + r; 2978 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*6 + r; 2979 remotePointsNew[m].rank = rrank; 2980 } 2981 for (r = 0; r < 1; ++r, ++m) { 2982 localPointsNew[m] = vStartNew + (eEnd - eStart) + (fEnd - fStart) + (p - cStart) + r; 2983 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+1] + rdepthSizeOld[n*(depth+1)+depth-1] + (rp - rcStart[n]) + r; 2984 remotePointsNew[m].rank = rrank; 2985 } 2986 } 2987 break; 2988 default: 2989 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2990 } 2991 } 2992 ierr = ISRestoreIndices(processRanks, &neighbors);CHKERRQ(ierr); 2993 ierr = ISDestroy(&processRanks);CHKERRQ(ierr); 2994 ierr = PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2995 ierr = PetscFree5(rdepthSize,rvStartNew,reStartNew,rfStartNew,rcStartNew);CHKERRQ(ierr); 2996 ierr = PetscFree7(depthSizeOld,rdepthSizeOld,rdepthMaxOld,rvStart,reStart,rfStart,rcStart);CHKERRQ(ierr); 2997 PetscFunctionReturn(0); 2998 } 2999 3000 #undef __FUNCT__ 3001 #define __FUNCT__ "CellRefinerCreateLabels" 3002 PetscErrorCode CellRefinerCreateLabels(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 3003 { 3004 PetscInt numLabels, l; 3005 PetscInt depth, newp, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eStartNew, eEnd, eMax, r; 3006 PetscErrorCode ierr; 3007 3008 PetscFunctionBegin; 3009 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3010 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 3011 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 3012 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 3013 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3014 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 3015 ierr = DMPlexGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 3016 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 3017 switch (refiner) { 3018 case 3: 3019 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 3020 cMax = PetscMin(cEnd, cMax); 3021 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 3022 fMax = PetscMin(fEnd, fMax); 3023 } 3024 for (l = 0; l < numLabels; ++l) { 3025 DMLabel label, labelNew; 3026 const char *lname; 3027 PetscBool isDepth; 3028 IS valueIS; 3029 const PetscInt *values; 3030 PetscInt numValues, val; 3031 3032 ierr = DMPlexGetLabelName(dm, l, &lname);CHKERRQ(ierr); 3033 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 3034 if (isDepth) continue; 3035 ierr = DMPlexCreateLabel(rdm, lname);CHKERRQ(ierr); 3036 ierr = DMPlexGetLabel(dm, lname, &label);CHKERRQ(ierr); 3037 ierr = DMPlexGetLabel(rdm, lname, &labelNew);CHKERRQ(ierr); 3038 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 3039 ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr); 3040 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 3041 for (val = 0; val < numValues; ++val) { 3042 IS pointIS; 3043 const PetscInt *points; 3044 PetscInt numPoints, n; 3045 3046 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 3047 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 3048 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 3049 for (n = 0; n < numPoints; ++n) { 3050 const PetscInt p = points[n]; 3051 switch (refiner) { 3052 case 1: 3053 /* Simplicial 2D */ 3054 if ((p >= vStart) && (p < vEnd)) { 3055 /* Old vertices stay the same */ 3056 newp = vStartNew + (p - vStart); 3057 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3058 } else if ((p >= fStart) && (p < fEnd)) { 3059 /* Old faces add new faces and vertex */ 3060 newp = vStartNew + (vEnd - vStart) + (p - fStart); 3061 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3062 for (r = 0; r < 2; ++r) { 3063 newp = fStartNew + (p - fStart)*2 + r; 3064 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3065 } 3066 } else if ((p >= cStart) && (p < cEnd)) { 3067 /* Old cells add new cells and interior faces */ 3068 for (r = 0; r < 4; ++r) { 3069 newp = cStartNew + (p - cStart)*4 + r; 3070 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3071 } 3072 for (r = 0; r < 3; ++r) { 3073 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 3074 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3075 } 3076 } 3077 break; 3078 case 2: 3079 /* Hex 2D */ 3080 if ((p >= vStart) && (p < vEnd)) { 3081 /* Old vertices stay the same */ 3082 newp = vStartNew + (p - vStart); 3083 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3084 } else if ((p >= fStart) && (p < fEnd)) { 3085 /* Old faces add new faces and vertex */ 3086 newp = vStartNew + (vEnd - vStart) + (p - fStart); 3087 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3088 for (r = 0; r < 2; ++r) { 3089 newp = fStartNew + (p - fStart)*2 + r; 3090 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3091 } 3092 } else if ((p >= cStart) && (p < cEnd)) { 3093 /* Old cells add new cells and interior faces and vertex */ 3094 for (r = 0; r < 4; ++r) { 3095 newp = cStartNew + (p - cStart)*4 + r; 3096 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3097 } 3098 for (r = 0; r < 4; ++r) { 3099 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 3100 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3101 } 3102 newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (p - cStart); 3103 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3104 } 3105 break; 3106 case 3: 3107 /* Hybrid simplicial 2D */ 3108 if ((p >= vStart) && (p < vEnd)) { 3109 /* Old vertices stay the same */ 3110 newp = vStartNew + (p - vStart); 3111 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3112 } else if ((p >= fStart) && (p < fMax)) { 3113 /* Old interior faces add new faces and vertex */ 3114 newp = vStartNew + (vEnd - vStart) + (p - fStart); 3115 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3116 for (r = 0; r < 2; ++r) { 3117 newp = fStartNew + (p - fStart)*2 + r; 3118 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3119 } 3120 } else if ((p >= fMax) && (p < fEnd)) { 3121 /* Old hybrid faces stay the same */ 3122 newp = fStartNew + (fMax - fStart)*2 + (p - fMax); 3123 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3124 } else if ((p >= cStart) && (p < cMax)) { 3125 /* Old interior cells add new cells and interior faces */ 3126 for (r = 0; r < 4; ++r) { 3127 newp = cStartNew + (p - cStart)*4 + r; 3128 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3129 } 3130 for (r = 0; r < 3; ++r) { 3131 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 3132 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3133 } 3134 } else if ((p >= cMax) && (p < cEnd)) { 3135 /* Old hybrid cells add new cells and hybrid face */ 3136 for (r = 0; r < 2; ++r) { 3137 newp = cStartNew + (cMax - cStart)*4 + (p - cMax)*2 + r; 3138 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3139 } 3140 newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 3141 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3142 } 3143 break; 3144 case 5: 3145 /* Simplicial 3D */ 3146 if ((p >= vStart) && (p < vEnd)) { 3147 /* Old vertices stay the same */ 3148 newp = vStartNew + (p - vStart); 3149 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3150 } else if ((p >= eStart) && (p < eEnd)) { 3151 /* Old edges add new edges and vertex */ 3152 for (r = 0; r < 2; ++r) { 3153 newp = eStartNew + (p - eStart)*2 + r; 3154 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3155 } 3156 newp = vStartNew + (vEnd - vStart) + (p - eStart); 3157 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3158 } else if ((p >= fStart) && (p < fEnd)) { 3159 /* Old faces add new faces and edges */ 3160 for (r = 0; r < 4; ++r) { 3161 newp = fStartNew + (p - fStart)*4 + r; 3162 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3163 } 3164 for (r = 0; r < 3; ++r) { 3165 newp = eStartNew + (eEnd - eStart)*2 + (p - fStart)*3 + r; 3166 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3167 } 3168 } else if ((p >= cStart) && (p < cEnd)) { 3169 /* Old cells add new cells and interior faces and edges */ 3170 for (r = 0; r < 8; ++r) { 3171 newp = cStartNew + (p - cStart)*8 + r; 3172 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3173 } 3174 for (r = 0; r < 8; ++r) { 3175 newp = fStartNew + (fEnd - fStart)*4 + (p - cStart)*8 + r; 3176 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3177 } 3178 for (r = 0; r < 1; ++r) { 3179 newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (p - cStart)*1 + r; 3180 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3181 } 3182 } 3183 break; 3184 case 6: 3185 /* Hex 3D */ 3186 if ((p >= vStart) && (p < vEnd)) { 3187 /* Old vertices stay the same */ 3188 newp = vStartNew + (p - vStart); 3189 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3190 } else if ((p >= fStart) && (p < fEnd)) { 3191 /* Old edges add new edges and vertex */ 3192 for (r = 0; r < 2; ++r) { 3193 newp = eStartNew + (p - eStart)*2 + r; 3194 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3195 } 3196 newp = vStartNew + (vEnd - vStart) + (p - eStart); 3197 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3198 } else if ((p >= fStart) && (p < fEnd)) { 3199 /* Old faces add new faces, edges, and vertex */ 3200 for (r = 0; r < 4; ++r) { 3201 newp = fStartNew + (p - fStart)*4 + r; 3202 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3203 } 3204 for (r = 0; r < 4; ++r) { 3205 newp = eStartNew + (eEnd - eStart)*2 + (p - fStart)*4 + r; 3206 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3207 } 3208 newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (p - fStart); 3209 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3210 } else if ((p >= cStart) && (p < cEnd)) { 3211 /* Old cells add new cells, faces, edges, and vertex */ 3212 for (r = 0; r < 8; ++r) { 3213 newp = cStartNew + (p - cStart)*8 + r; 3214 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3215 } 3216 for (r = 0; r < 12; ++r) { 3217 newp = fStartNew + (fEnd - fStart)*4 + (p - cStart)*12 + r; 3218 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3219 } 3220 for (r = 0; r < 6; ++r) { 3221 newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*4 + (p - cStart)*6 + r; 3222 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3223 } 3224 newp = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (fEnd - fStart) + (p - cStart); 3225 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 3226 } 3227 break; 3228 default: 3229 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 3230 } 3231 } 3232 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 3233 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 3234 } 3235 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 3236 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 3237 if (0) { 3238 ierr = PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);CHKERRQ(ierr); 3239 ierr = DMLabelView(labelNew, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 3240 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 3241 } 3242 } 3243 PetscFunctionReturn(0); 3244 } 3245 3246 #undef __FUNCT__ 3247 #define __FUNCT__ "DMPlexRefineUniform_Internal" 3248 /* This will only work for interpolated meshes */ 3249 PetscErrorCode DMPlexRefineUniform_Internal(DM dm, CellRefiner cellRefiner, DM *dmRefined) 3250 { 3251 DM rdm; 3252 PetscInt *depthSize; 3253 PetscInt dim, depth = 0, d, pStart = 0, pEnd = 0; 3254 PetscErrorCode ierr; 3255 3256 PetscFunctionBegin; 3257 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &rdm);CHKERRQ(ierr); 3258 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3259 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 3260 ierr = DMPlexSetDimension(rdm, dim);CHKERRQ(ierr); 3261 /* Calculate number of new points of each depth */ 3262 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3263 ierr = PetscMalloc((depth+1) * sizeof(PetscInt), &depthSize);CHKERRQ(ierr); 3264 ierr = PetscMemzero(depthSize, (depth+1) * sizeof(PetscInt));CHKERRQ(ierr); 3265 ierr = CellRefinerGetSizes(cellRefiner, dm, depthSize);CHKERRQ(ierr); 3266 /* Step 1: Set chart */ 3267 for (d = 0; d <= depth; ++d) pEnd += depthSize[d]; 3268 ierr = DMPlexSetChart(rdm, pStart, pEnd);CHKERRQ(ierr); 3269 /* Step 2: Set cone/support sizes */ 3270 ierr = CellRefinerSetConeSizes(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3271 /* Step 3: Setup refined DM */ 3272 ierr = DMSetUp(rdm);CHKERRQ(ierr); 3273 /* Step 4: Set cones and supports */ 3274 ierr = CellRefinerSetCones(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3275 /* Step 5: Stratify */ 3276 ierr = DMPlexStratify(rdm);CHKERRQ(ierr); 3277 /* Step 6: Set coordinates for vertices */ 3278 ierr = CellRefinerSetCoordinates(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3279 /* Step 7: Create pointSF */ 3280 ierr = CellRefinerCreateSF(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3281 /* Step 8: Create labels */ 3282 ierr = CellRefinerCreateLabels(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 3283 ierr = PetscFree(depthSize);CHKERRQ(ierr); 3284 3285 *dmRefined = rdm; 3286 PetscFunctionReturn(0); 3287 } 3288 3289 #undef __FUNCT__ 3290 #define __FUNCT__ "DMPlexSetRefinementUniform" 3291 PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform) 3292 { 3293 DM_Plex *mesh = (DM_Plex*) dm->data; 3294 3295 PetscFunctionBegin; 3296 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3297 mesh->refinementUniform = refinementUniform; 3298 PetscFunctionReturn(0); 3299 } 3300 3301 #undef __FUNCT__ 3302 #define __FUNCT__ "DMPlexGetRefinementUniform" 3303 PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform) 3304 { 3305 DM_Plex *mesh = (DM_Plex*) dm->data; 3306 3307 PetscFunctionBegin; 3308 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3309 PetscValidPointer(refinementUniform, 2); 3310 *refinementUniform = mesh->refinementUniform; 3311 PetscFunctionReturn(0); 3312 } 3313 3314 #undef __FUNCT__ 3315 #define __FUNCT__ "DMPlexSetRefinementLimit" 3316 PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit) 3317 { 3318 DM_Plex *mesh = (DM_Plex*) dm->data; 3319 3320 PetscFunctionBegin; 3321 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3322 mesh->refinementLimit = refinementLimit; 3323 PetscFunctionReturn(0); 3324 } 3325 3326 #undef __FUNCT__ 3327 #define __FUNCT__ "DMPlexGetRefinementLimit" 3328 PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit) 3329 { 3330 DM_Plex *mesh = (DM_Plex*) dm->data; 3331 3332 PetscFunctionBegin; 3333 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3334 PetscValidPointer(refinementLimit, 2); 3335 /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */ 3336 *refinementLimit = mesh->refinementLimit; 3337 PetscFunctionReturn(0); 3338 } 3339 3340 #undef __FUNCT__ 3341 #define __FUNCT__ "DMPlexGetCellRefiner_Internal" 3342 PetscErrorCode DMPlexGetCellRefiner_Internal(DM dm, CellRefiner *cellRefiner) 3343 { 3344 PetscInt dim, cStart, coneSize, cMax; 3345 PetscErrorCode ierr; 3346 3347 PetscFunctionBegin; 3348 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 3349 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 3350 ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr); 3351 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 3352 switch (dim) { 3353 case 2: 3354 switch (coneSize) { 3355 case 3: 3356 if (cMax >= 0) *cellRefiner = 3; /* Hybrid */ 3357 else *cellRefiner = 1; /* Triangular */ 3358 break; 3359 case 4: 3360 if (cMax >= 0) *cellRefiner = 4; /* Hybrid */ 3361 else *cellRefiner = 2; /* Quadrilateral */ 3362 break; 3363 default: 3364 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown coneSize %d in dimension %d for cell refiner", coneSize, dim); 3365 } 3366 break; 3367 case 3: 3368 switch (coneSize) { 3369 case 4: 3370 if (cMax >= 0) *cellRefiner = 7; /* Hybrid */ 3371 else *cellRefiner = 5; /* Tetrahedral */ 3372 break; 3373 case 6: 3374 if (cMax >= 0) *cellRefiner = 8; /* Hybrid */ 3375 else *cellRefiner = 6; /* hexahedral */ 3376 break; 3377 default: 3378 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown coneSize %d in dimension %d for cell refiner", coneSize, dim); 3379 } 3380 break; 3381 default: 3382 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown dimension %d for cell refiner", dim); 3383 } 3384 PetscFunctionReturn(0); 3385 } 3386