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