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