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 7: 103 /* Hybrid Simplicial 3D */ 104 ierr = DMPlexGetNumHybridFaces_Internal(dm, &numHybridFaces);CHKERRQ(ierr); 105 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 106 cMax = PetscMin(cEnd, cMax); 107 if (eMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No edge maximum specified in hybrid mesh"); 108 eMax = PetscMin(eEnd, eMax); 109 depthSize[0] = vEnd - vStart + eMax - eStart; /* Add a vertex on every edge, but not hybrid edges */ 110 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 */ 111 depthSize[2] = 4*(fEnd - fStart) + 8*(cEnd - cStart); /* Every face split into 4 faces and 8 faces are added for each cell */ 112 depthSize[3] = 8*(cMax - cStart) + 4*(cEnd - cMax); /* Every interior cell split into 8 cells, and every hybrid cell split into 4 cells */ 113 break; 114 default: 115 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 116 } 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "CellRefinerSetConeSizes" 122 PetscErrorCode CellRefinerSetConeSizes(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 123 { 124 PetscInt depth, cStart, cStartNew, cEnd, cMax, c, vStart, vStartNew, vEnd, vMax, v, fStart, fStartNew, fEnd, fMax, f, eStart, eStartNew, eEnd, eMax, e, r; 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 129 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 130 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 131 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 132 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 133 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 134 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 135 switch (refiner) { 136 case 1: 137 /* Simplicial 2D */ 138 /* All cells have 3 faces */ 139 for (c = cStart; c < cEnd; ++c) { 140 for (r = 0; r < 4; ++r) { 141 const PetscInt newp = (c - cStart)*4 + r; 142 143 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 144 } 145 } 146 /* Split faces have 2 vertices and the same cells as the parent */ 147 for (f = fStart; f < fEnd; ++f) { 148 for (r = 0; r < 2; ++r) { 149 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 150 PetscInt size; 151 152 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 153 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 154 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 155 } 156 } 157 /* Interior faces have 2 vertices and 2 cells */ 158 for (c = cStart; c < cEnd; ++c) { 159 for (r = 0; r < 3; ++r) { 160 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + r; 161 162 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 163 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 164 } 165 } 166 /* Old vertices have identical supports */ 167 for (v = vStart; v < vEnd; ++v) { 168 const PetscInt newp = vStartNew + (v - vStart); 169 PetscInt size; 170 171 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 172 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 173 } 174 /* Face vertices have 2 + cells*2 supports */ 175 for (f = fStart; f < fEnd; ++f) { 176 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 177 PetscInt size; 178 179 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 180 ierr = DMPlexSetSupportSize(rdm, newp, 2 + size*2);CHKERRQ(ierr); 181 } 182 break; 183 case 2: 184 /* Hex 2D */ 185 /* All cells have 4 faces */ 186 for (c = cStart; c < cEnd; ++c) { 187 for (r = 0; r < 4; ++r) { 188 const PetscInt newp = (c - cStart)*4 + r; 189 190 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 191 } 192 } 193 /* Split faces have 2 vertices and the same cells as the parent */ 194 for (f = fStart; f < fEnd; ++f) { 195 for (r = 0; r < 2; ++r) { 196 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 197 PetscInt size; 198 199 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 200 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 201 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 202 } 203 } 204 /* Interior faces have 2 vertices and 2 cells */ 205 for (c = cStart; c < cEnd; ++c) { 206 for (r = 0; r < 4; ++r) { 207 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + r; 208 209 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 210 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 211 } 212 } 213 /* Old vertices have identical supports */ 214 for (v = vStart; v < vEnd; ++v) { 215 const PetscInt newp = vStartNew + (v - vStart); 216 PetscInt size; 217 218 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 219 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 220 } 221 /* Face vertices have 2 + cells supports */ 222 for (f = fStart; f < fEnd; ++f) { 223 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 224 PetscInt size; 225 226 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 227 ierr = DMPlexSetSupportSize(rdm, newp, 2 + size);CHKERRQ(ierr); 228 } 229 /* Cell vertices have 4 supports */ 230 for (c = cStart; c < cEnd; ++c) { 231 const PetscInt newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 232 233 ierr = DMPlexSetSupportSize(rdm, newp, 4);CHKERRQ(ierr); 234 } 235 break; 236 case 3: 237 /* Hybrid Simplicial 2D */ 238 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 239 cMax = PetscMin(cEnd, cMax); 240 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 241 fMax = PetscMin(fEnd, fMax); 242 ierr = DMPlexSetHybridBounds(rdm, cStartNew + (cMax - cStart)*4, fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 243 /* Interior cells have 3 faces */ 244 for (c = cStart; c < cMax; ++c) { 245 for (r = 0; r < 4; ++r) { 246 const PetscInt newp = cStartNew + (c - cStart)*4 + r; 247 248 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 249 } 250 } 251 /* Hybrid cells have 4 faces */ 252 for (c = cMax; c < cEnd; ++c) { 253 for (r = 0; r < 2; ++r) { 254 const PetscInt newp = cStartNew + (cMax - cStart)*4 + (c - cMax)*2 + r; 255 256 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 257 } 258 } 259 /* Interior split faces have 2 vertices and the same cells as the parent */ 260 for (f = fStart; f < fMax; ++f) { 261 for (r = 0; r < 2; ++r) { 262 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 263 PetscInt size; 264 265 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 266 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 267 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 268 } 269 } 270 /* Interior cell faces have 2 vertices and 2 cells */ 271 for (c = cStart; c < cMax; ++c) { 272 for (r = 0; r < 3; ++r) { 273 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + r; 274 275 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 276 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 277 } 278 } 279 /* Hybrid faces have 2 vertices and the same cells */ 280 for (f = fMax; f < fEnd; ++f) { 281 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (f - fMax); 282 PetscInt size; 283 284 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 285 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 286 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 287 } 288 /* Hybrid cell faces have 2 vertices and 2 cells */ 289 for (c = cMax; c < cEnd; ++c) { 290 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 291 292 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 293 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 294 } 295 /* Old vertices have identical supports */ 296 for (v = vStart; v < vEnd; ++v) { 297 const PetscInt newp = vStartNew + (v - vStart); 298 PetscInt size; 299 300 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 301 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 302 } 303 /* Face vertices have 2 + (2 interior, 1 hybrid) supports */ 304 for (f = fStart; f < fMax; ++f) { 305 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 306 const PetscInt *support; 307 PetscInt size, newSize = 2, s; 308 309 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 310 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 311 for (s = 0; s < size; ++s) { 312 if (support[s] >= cMax) newSize += 1; 313 else newSize += 2; 314 } 315 ierr = DMPlexSetSupportSize(rdm, newp, newSize);CHKERRQ(ierr); 316 } 317 break; 318 case 5: 319 /* Simplicial 3D */ 320 /* All cells have 4 faces */ 321 for (c = cStart; c < cEnd; ++c) { 322 for (r = 0; r < 8; ++r) { 323 const PetscInt newp = (c - cStart)*8 + r; 324 325 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 326 } 327 } 328 /* Split faces have 3 edges and the same cells as the parent */ 329 for (f = fStart; f < fEnd; ++f) { 330 for (r = 0; r < 4; ++r) { 331 const PetscInt newp = fStartNew + (f - fStart)*4 + r; 332 PetscInt size; 333 334 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 335 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 336 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 337 } 338 } 339 /* Interior faces have 3 edges and 2 cells */ 340 for (c = cStart; c < cEnd; ++c) { 341 for (r = 0; r < 8; ++r) { 342 const PetscInt newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + r; 343 344 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 345 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 346 } 347 } 348 /* Split edges have 2 vertices and the same faces */ 349 for (e = eStart; e < eEnd; ++e) { 350 for (r = 0; r < 2; ++r) { 351 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 352 PetscInt size; 353 354 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 355 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 356 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 357 } 358 } 359 /* Face edges have 2 vertices and 2+cells*(1/2) faces */ 360 for (f = fStart; f < fEnd; ++f) { 361 for (r = 0; r < 3; ++r) { 362 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + r; 363 const PetscInt *cone, *ornt, *support, eint[4] = {1, 0, 2, 0}; 364 PetscInt coneSize, c, supportSize, s, er, intFaces = 0; 365 366 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 367 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 368 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 369 for (s = 0; s < supportSize; ++s) { 370 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 371 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 372 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 373 for (c = 0; c < coneSize; ++c) {if (cone[c] == f) break;} 374 /* Here we want to determine whether edge newp contains a vertex which is part of the cross-tet edge */ 375 er = ornt[c] < 0 ? (-(ornt[c]+1) + 2-r)%3 : (ornt[c] + r)%3; 376 if (er == eint[c]) { 377 intFaces += 1; 378 } else { 379 intFaces += 2; 380 } 381 } 382 ierr = DMPlexSetSupportSize(rdm, newp, 2+intFaces);CHKERRQ(ierr); 383 } 384 } 385 /* Interior edges have 2 vertices and 4 faces */ 386 for (c = cStart; c < cEnd; ++c) { 387 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 388 389 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 390 ierr = DMPlexSetSupportSize(rdm, newp, 4);CHKERRQ(ierr); 391 } 392 /* Old vertices have identical supports */ 393 for (v = vStart; v < vEnd; ++v) { 394 const PetscInt newp = vStartNew + (v - vStart); 395 PetscInt size; 396 397 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 398 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 399 } 400 /* Edge vertices have 2 + faces*2 + cells*0/1 supports */ 401 for (e = eStart; e < eEnd; ++e) { 402 const PetscInt newp = vStartNew + (vEnd - vStart) + (e - eStart); 403 PetscInt size, *star = NULL, starSize, s, cellSize = 0; 404 405 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 406 ierr = DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 407 for (s = 0; s < starSize*2; s += 2) { 408 const PetscInt *cone, *ornt; 409 PetscInt e01, e23; 410 411 if ((star[s] >= cStart) && (star[s] < cEnd)) { 412 /* Check edge 0-1 */ 413 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 414 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 415 ierr = DMPlexGetCone(dm, cone[0], &cone);CHKERRQ(ierr); 416 e01 = cone[ornt[0] < 0 ? (-(ornt[0]+1) + 0)%3 : ornt[0]]; 417 /* Check edge 2-3 */ 418 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 419 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 420 ierr = DMPlexGetCone(dm, cone[3], &cone);CHKERRQ(ierr); 421 e23 = cone[ornt[3] < 0 ? (-(ornt[3]+1) + 2)%3 : (ornt[3] + 2)%3]; 422 if ((e01 == e) || (e23 == e)) ++cellSize; 423 } 424 } 425 ierr = DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 426 ierr = DMPlexSetSupportSize(rdm, newp, 2 + size*2 + cellSize);CHKERRQ(ierr); 427 } 428 break; 429 default: 430 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 431 } 432 PetscFunctionReturn(0); 433 } 434 435 #undef __FUNCT__ 436 #define __FUNCT__ "CellRefinerSetCones" 437 PetscErrorCode CellRefinerSetCones(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 438 { 439 const PetscInt *faces, cellInd[4] = {0, 1, 2, 3}; 440 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; 441 PetscInt maxSupportSize, *supportRef; 442 PetscErrorCode ierr; 443 444 PetscFunctionBegin; 445 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 446 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 447 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 448 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 449 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 450 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 451 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 452 ierr = GetDepthEnd_Private(depth, depthSize, &cEndNew, &fEndNew, &eEndNew, &vEndNew);CHKERRQ(ierr); 453 switch (refiner) { 454 case 1: 455 /* Simplicial 2D */ 456 /* 457 2 458 |\ 459 | \ 460 | \ 461 | \ 462 | C \ 463 | \ 464 | \ 465 2---1---1 466 |\ D / \ 467 | 2 0 \ 468 |A \ / B \ 469 0---0-------1 470 */ 471 /* All cells have 3 faces */ 472 for (c = cStart; c < cEnd; ++c) { 473 const PetscInt newp = cStartNew + (c - cStart)*4; 474 const PetscInt *cone, *ornt; 475 PetscInt coneNew[3], orntNew[3]; 476 477 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 478 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 479 /* A triangle */ 480 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 481 orntNew[0] = ornt[0]; 482 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 2; 483 orntNew[1] = -2; 484 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 485 orntNew[2] = ornt[2]; 486 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 487 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 488 #if 1 489 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); 490 for (p = 0; p < 3; ++p) { 491 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); 492 } 493 #endif 494 /* B triangle */ 495 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 496 orntNew[0] = ornt[0]; 497 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 498 orntNew[1] = ornt[1]; 499 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 0; 500 orntNew[2] = -2; 501 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 502 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 503 #if 1 504 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); 505 for (p = 0; p < 3; ++p) { 506 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); 507 } 508 #endif 509 /* C triangle */ 510 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 1; 511 orntNew[0] = -2; 512 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 513 orntNew[1] = ornt[1]; 514 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 515 orntNew[2] = ornt[2]; 516 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 517 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 518 #if 1 519 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); 520 for (p = 0; p < 3; ++p) { 521 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); 522 } 523 #endif 524 /* D triangle */ 525 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 0; 526 orntNew[0] = 0; 527 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 1; 528 orntNew[1] = 0; 529 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 2; 530 orntNew[2] = 0; 531 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 532 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 533 #if 1 534 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); 535 for (p = 0; p < 3; ++p) { 536 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); 537 } 538 #endif 539 } 540 /* Split faces have 2 vertices and the same cells as the parent */ 541 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 542 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 543 for (f = fStart; f < fEnd; ++f) { 544 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 545 546 for (r = 0; r < 2; ++r) { 547 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 548 const PetscInt *cone, *ornt, *support; 549 PetscInt coneNew[2], coneSize, c, supportSize, s; 550 551 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 552 coneNew[0] = vStartNew + (cone[0] - vStart); 553 coneNew[1] = vStartNew + (cone[1] - vStart); 554 coneNew[(r+1)%2] = newv; 555 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 556 #if 1 557 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 558 for (p = 0; p < 2; ++p) { 559 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); 560 } 561 #endif 562 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 563 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 564 for (s = 0; s < supportSize; ++s) { 565 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 566 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 567 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 568 for (c = 0; c < coneSize; ++c) { 569 if (cone[c] == f) break; 570 } 571 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (ornt[c] < 0 ? (c+1-r)%3 : (c+r)%3); 572 } 573 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 574 #if 1 575 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 576 for (p = 0; p < supportSize; ++p) { 577 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); 578 } 579 #endif 580 } 581 } 582 /* Interior faces have 2 vertices and 2 cells */ 583 for (c = cStart; c < cEnd; ++c) { 584 const PetscInt *cone; 585 586 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 587 for (r = 0; r < 3; ++r) { 588 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + r; 589 PetscInt coneNew[2]; 590 PetscInt supportNew[2]; 591 592 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 593 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - fStart); 594 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 595 #if 1 596 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 597 for (p = 0; p < 2; ++p) { 598 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); 599 } 600 #endif 601 supportNew[0] = (c - cStart)*4 + (r+1)%3; 602 supportNew[1] = (c - cStart)*4 + 3; 603 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 604 #if 1 605 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 606 for (p = 0; p < 2; ++p) { 607 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); 608 } 609 #endif 610 } 611 } 612 /* Old vertices have identical supports */ 613 for (v = vStart; v < vEnd; ++v) { 614 const PetscInt newp = vStartNew + (v - vStart); 615 const PetscInt *support, *cone; 616 PetscInt size, s; 617 618 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 619 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 620 for (s = 0; s < size; ++s) { 621 PetscInt r = 0; 622 623 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 624 if (cone[1] == v) r = 1; 625 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 626 } 627 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 628 #if 1 629 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 630 for (p = 0; p < size; ++p) { 631 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); 632 } 633 #endif 634 } 635 /* Face vertices have 2 + cells*2 supports */ 636 for (f = fStart; f < fEnd; ++f) { 637 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 638 const PetscInt *cone, *support; 639 PetscInt size, s; 640 641 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 642 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 643 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 644 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 645 for (s = 0; s < size; ++s) { 646 PetscInt r = 0; 647 648 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 649 if (cone[1] == f) r = 1; 650 else if (cone[2] == f) r = 2; 651 supportRef[2+s*2+0] = fStartNew + (fEnd - fStart)*2 + (support[s] - cStart)*3 + (r+2)%3; 652 supportRef[2+s*2+1] = fStartNew + (fEnd - fStart)*2 + (support[s] - cStart)*3 + r; 653 } 654 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 655 #if 1 656 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 657 for (p = 0; p < 2+size*2; ++p) { 658 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); 659 } 660 #endif 661 } 662 ierr = PetscFree(supportRef);CHKERRQ(ierr); 663 break; 664 case 2: 665 /* Hex 2D */ 666 /* 667 3---------2---------2 668 | | | 669 | D 2 C | 670 | | | 671 3----3----0----1----1 672 | | | 673 | A 0 B | 674 | | | 675 0---------0---------1 676 */ 677 /* All cells have 4 faces */ 678 for (c = cStart; c < cEnd; ++c) { 679 const PetscInt newp = (c - cStart)*4; 680 const PetscInt *cone, *ornt; 681 PetscInt coneNew[4], orntNew[4]; 682 683 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 684 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 685 /* A quad */ 686 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 687 orntNew[0] = ornt[0]; 688 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 0; 689 orntNew[1] = 0; 690 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 3; 691 orntNew[2] = -2; 692 coneNew[3] = fStartNew + (cone[3] - fStart)*2 + (ornt[3] < 0 ? 0 : 1); 693 orntNew[3] = ornt[3]; 694 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 695 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 696 #if 1 697 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); 698 for (p = 0; p < 4; ++p) { 699 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); 700 } 701 #endif 702 /* B quad */ 703 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 704 orntNew[0] = ornt[0]; 705 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 706 orntNew[1] = ornt[1]; 707 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 1; 708 orntNew[2] = 0; 709 coneNew[3] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 0; 710 orntNew[3] = -2; 711 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 712 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 713 #if 1 714 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); 715 for (p = 0; p < 4; ++p) { 716 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); 717 } 718 #endif 719 /* C quad */ 720 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 1; 721 orntNew[0] = -2; 722 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 723 orntNew[1] = ornt[1]; 724 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 725 orntNew[2] = ornt[2]; 726 coneNew[3] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 2; 727 orntNew[3] = 0; 728 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 729 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 730 #if 1 731 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); 732 for (p = 0; p < 4; ++p) { 733 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); 734 } 735 #endif 736 /* D quad */ 737 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 3; 738 orntNew[0] = 0; 739 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 2; 740 orntNew[1] = -2; 741 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 742 orntNew[2] = ornt[2]; 743 coneNew[3] = fStartNew + (cone[3] - fStart)*2 + (ornt[3] < 0 ? 1 : 0); 744 orntNew[3] = ornt[3]; 745 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 746 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 747 #if 1 748 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); 749 for (p = 0; p < 4; ++p) { 750 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); 751 } 752 #endif 753 } 754 /* Split faces have 2 vertices and the same cells as the parent */ 755 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 756 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 757 for (f = fStart; f < fEnd; ++f) { 758 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 759 760 for (r = 0; r < 2; ++r) { 761 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 762 const PetscInt *cone, *ornt, *support; 763 PetscInt coneNew[2], coneSize, c, supportSize, s; 764 765 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 766 coneNew[0] = vStartNew + (cone[0] - vStart); 767 coneNew[1] = vStartNew + (cone[1] - vStart); 768 coneNew[(r+1)%2] = newv; 769 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 770 #if 1 771 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 772 for (p = 0; p < 2; ++p) { 773 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); 774 } 775 #endif 776 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 777 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 778 for (s = 0; s < supportSize; ++s) { 779 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 780 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 781 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 782 for (c = 0; c < coneSize; ++c) { 783 if (cone[c] == f) break; 784 } 785 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (ornt[c] < 0 ? (c+1-r)%4 : (c+r)%4); 786 } 787 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 788 #if 1 789 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 790 for (p = 0; p < supportSize; ++p) { 791 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); 792 } 793 #endif 794 } 795 } 796 /* Interior faces have 2 vertices and 2 cells */ 797 for (c = cStart; c < cEnd; ++c) { 798 const PetscInt *cone; 799 PetscInt coneNew[2], supportNew[2]; 800 801 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 802 for (r = 0; r < 4; ++r) { 803 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + r; 804 805 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 806 coneNew[1] = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 807 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 808 #if 1 809 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 810 for (p = 0; p < 2; ++p) { 811 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); 812 } 813 #endif 814 supportNew[0] = (c - cStart)*4 + r; 815 supportNew[1] = (c - cStart)*4 + (r+1)%4; 816 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 817 #if 1 818 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 819 for (p = 0; p < 2; ++p) { 820 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); 821 } 822 #endif 823 } 824 } 825 /* Old vertices have identical supports */ 826 for (v = vStart; v < vEnd; ++v) { 827 const PetscInt newp = vStartNew + (v - vStart); 828 const PetscInt *support, *cone; 829 PetscInt size, s; 830 831 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 832 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 833 for (s = 0; s < size; ++s) { 834 PetscInt r = 0; 835 836 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 837 if (cone[1] == v) r = 1; 838 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 839 } 840 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 841 #if 1 842 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 843 for (p = 0; p < size; ++p) { 844 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); 845 } 846 #endif 847 } 848 /* Face vertices have 2 + cells supports */ 849 for (f = fStart; f < fEnd; ++f) { 850 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 851 const PetscInt *cone, *support; 852 PetscInt size, s; 853 854 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 855 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 856 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 857 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 858 for (s = 0; s < size; ++s) { 859 PetscInt r = 0; 860 861 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 862 if (cone[1] == f) r = 1; 863 else if (cone[2] == f) r = 2; 864 else if (cone[3] == f) r = 3; 865 supportRef[2+s] = fStartNew + (fEnd - fStart)*2 + (support[s] - cStart)*4 + r; 866 } 867 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 868 #if 1 869 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 870 for (p = 0; p < 2+size; ++p) { 871 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); 872 } 873 #endif 874 } 875 /* Cell vertices have 4 supports */ 876 for (c = cStart; c < cEnd; ++c) { 877 const PetscInt newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 878 PetscInt supportNew[4]; 879 880 for (r = 0; r < 4; ++r) { 881 supportNew[r] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + r; 882 } 883 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 884 } 885 break; 886 case 3: 887 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 888 cMax = PetscMin(cEnd, cMax); 889 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 890 fMax = PetscMin(fEnd, fMax); 891 /* Interior cells have 3 faces */ 892 for (c = cStart; c < cMax; ++c) { 893 const PetscInt newp = cStartNew + (c - cStart)*4; 894 const PetscInt *cone, *ornt; 895 PetscInt coneNew[3], orntNew[3]; 896 897 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 898 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 899 /* A triangle */ 900 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 901 orntNew[0] = ornt[0]; 902 coneNew[1] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 2; 903 orntNew[1] = -2; 904 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 905 orntNew[2] = ornt[2]; 906 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 907 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 908 #if 1 909 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); 910 for (p = 0; p < 3; ++p) { 911 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); 912 } 913 #endif 914 /* B triangle */ 915 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 916 orntNew[0] = ornt[0]; 917 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 918 orntNew[1] = ornt[1]; 919 coneNew[2] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 0; 920 orntNew[2] = -2; 921 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 922 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 923 #if 1 924 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); 925 for (p = 0; p < 3; ++p) { 926 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); 927 } 928 #endif 929 /* C triangle */ 930 coneNew[0] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 1; 931 orntNew[0] = -2; 932 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 933 orntNew[1] = ornt[1]; 934 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 935 orntNew[2] = ornt[2]; 936 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 937 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 938 #if 1 939 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); 940 for (p = 0; p < 3; ++p) { 941 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); 942 } 943 #endif 944 /* D triangle */ 945 coneNew[0] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 0; 946 orntNew[0] = 0; 947 coneNew[1] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 1; 948 orntNew[1] = 0; 949 coneNew[2] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 2; 950 orntNew[2] = 0; 951 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 952 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 953 #if 1 954 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); 955 for (p = 0; p < 3; ++p) { 956 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); 957 } 958 #endif 959 } 960 /* 961 2----3----3 962 | | 963 | B | 964 | | 965 0----4--- 1 966 | | 967 | A | 968 | | 969 0----2----1 970 */ 971 /* Hybrid cells have 4 faces */ 972 for (c = cMax; c < cEnd; ++c) { 973 const PetscInt newp = cStartNew + (cMax - cStart)*4 + (c - cMax)*2; 974 const PetscInt *cone, *ornt; 975 PetscInt coneNew[4], orntNew[4]; 976 977 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 978 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 979 /* A quad */ 980 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 981 orntNew[0] = ornt[0]; 982 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 983 orntNew[1] = ornt[1]; 984 coneNew[2] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (cone[2] - fMax); 985 orntNew[2] = 0; 986 coneNew[3] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 987 orntNew[3] = 0; 988 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 989 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 990 #if 1 991 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); 992 for (p = 0; p < 4; ++p) { 993 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); 994 } 995 #endif 996 /* B quad */ 997 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 998 orntNew[0] = ornt[0]; 999 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 1000 orntNew[1] = ornt[1]; 1001 coneNew[2] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 1002 orntNew[2] = 0; 1003 coneNew[3] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (cone[3] - fMax); 1004 orntNew[3] = 0; 1005 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1006 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1007 #if 1 1008 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); 1009 for (p = 0; p < 4; ++p) { 1010 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); 1011 } 1012 #endif 1013 } 1014 /* Interior split faces have 2 vertices and the same cells as the parent */ 1015 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 1016 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 1017 for (f = fStart; f < fMax; ++f) { 1018 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 1019 1020 for (r = 0; r < 2; ++r) { 1021 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 1022 const PetscInt *cone, *ornt, *support; 1023 PetscInt coneNew[2], coneSize, c, supportSize, s; 1024 1025 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1026 coneNew[0] = vStartNew + (cone[0] - vStart); 1027 coneNew[1] = vStartNew + (cone[1] - vStart); 1028 coneNew[(r+1)%2] = newv; 1029 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1030 #if 1 1031 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1032 for (p = 0; p < 2; ++p) { 1033 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); 1034 } 1035 #endif 1036 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1037 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1038 for (s = 0; s < supportSize; ++s) { 1039 if (support[s] >= cMax) { 1040 supportRef[s] = cStartNew + (cMax - cStart)*4 + (support[s] - cMax)*2 + r; 1041 } else { 1042 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1043 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1044 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1045 for (c = 0; c < coneSize; ++c) { 1046 if (cone[c] == f) break; 1047 } 1048 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (ornt[c] < 0 ? (c+1-r)%3 : (c+r)%3); 1049 } 1050 } 1051 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1052 #if 1 1053 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1054 for (p = 0; p < supportSize; ++p) { 1055 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); 1056 } 1057 #endif 1058 } 1059 } 1060 /* Interior cell faces have 2 vertices and 2 cells */ 1061 for (c = cStart; c < cMax; ++c) { 1062 const PetscInt *cone; 1063 1064 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1065 for (r = 0; r < 3; ++r) { 1066 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + r; 1067 PetscInt coneNew[2]; 1068 PetscInt supportNew[2]; 1069 1070 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 1071 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - fStart); 1072 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1073 #if 1 1074 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1075 for (p = 0; p < 2; ++p) { 1076 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); 1077 } 1078 #endif 1079 supportNew[0] = (c - cStart)*4 + (r+1)%3; 1080 supportNew[1] = (c - cStart)*4 + 3; 1081 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1082 #if 1 1083 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1084 for (p = 0; p < 2; ++p) { 1085 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); 1086 } 1087 #endif 1088 } 1089 } 1090 /* Interior hybrid faces have 2 vertices and the same cells */ 1091 for (f = fMax; f < fEnd; ++f) { 1092 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (f - fMax); 1093 const PetscInt *cone; 1094 const PetscInt *support; 1095 PetscInt coneNew[2]; 1096 PetscInt supportNew[2]; 1097 PetscInt size, s, r; 1098 1099 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1100 coneNew[0] = vStartNew + (cone[0] - vStart); 1101 coneNew[1] = vStartNew + (cone[1] - vStart); 1102 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1103 #if 1 1104 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1105 for (p = 0; p < 2; ++p) { 1106 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); 1107 } 1108 #endif 1109 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 1110 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1111 for (s = 0; s < size; ++s) { 1112 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1113 for (r = 0; r < 2; ++r) { 1114 if (cone[r+2] == f) break; 1115 } 1116 supportNew[s] = (cMax - cStart)*4 + (support[s] - cMax)*2 + r; 1117 } 1118 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1119 #if 1 1120 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1121 for (p = 0; p < size; ++p) { 1122 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); 1123 } 1124 #endif 1125 } 1126 /* Cell hybrid faces have 2 vertices and 2 cells */ 1127 for (c = cMax; c < cEnd; ++c) { 1128 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 1129 const PetscInt *cone; 1130 PetscInt coneNew[2]; 1131 PetscInt supportNew[2]; 1132 1133 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1134 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[0] - fStart); 1135 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[1] - fStart); 1136 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1137 #if 1 1138 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1139 for (p = 0; p < 2; ++p) { 1140 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); 1141 } 1142 #endif 1143 supportNew[0] = (cMax - cStart)*4 + (c - cMax)*2 + 0; 1144 supportNew[1] = (cMax - cStart)*4 + (c - cMax)*2 + 1; 1145 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1146 #if 1 1147 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1148 for (p = 0; p < 2; ++p) { 1149 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); 1150 } 1151 #endif 1152 } 1153 /* Old vertices have identical supports */ 1154 for (v = vStart; v < vEnd; ++v) { 1155 const PetscInt newp = vStartNew + (v - vStart); 1156 const PetscInt *support, *cone; 1157 PetscInt size, s; 1158 1159 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 1160 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 1161 for (s = 0; s < size; ++s) { 1162 if (support[s] >= fMax) { 1163 supportRef[s] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (support[s] - fMax); 1164 } else { 1165 PetscInt r = 0; 1166 1167 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1168 if (cone[1] == v) r = 1; 1169 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 1170 } 1171 } 1172 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1173 #if 1 1174 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1175 for (p = 0; p < size; ++p) { 1176 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); 1177 } 1178 #endif 1179 } 1180 /* Face vertices have 2 + (2 interior, 1 hybrid) supports */ 1181 for (f = fStart; f < fMax; ++f) { 1182 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 1183 const PetscInt *cone, *support; 1184 PetscInt size, newSize = 2, s; 1185 1186 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 1187 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1188 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 1189 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 1190 for (s = 0; s < size; ++s) { 1191 PetscInt r = 0; 1192 1193 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1194 if (support[s] >= cMax) { 1195 supportRef[newSize+0] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (support[s] - cMax); 1196 1197 newSize += 1; 1198 } else { 1199 if (cone[1] == f) r = 1; 1200 else if (cone[2] == f) r = 2; 1201 supportRef[newSize+0] = fStartNew + (fMax - fStart)*2 + (support[s] - cStart)*3 + (r+2)%3; 1202 supportRef[newSize+1] = fStartNew + (fMax - fStart)*2 + (support[s] - cStart)*3 + r; 1203 1204 newSize += 2; 1205 } 1206 } 1207 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1208 #if 1 1209 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1210 for (p = 0; p < newSize; ++p) { 1211 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); 1212 } 1213 #endif 1214 } 1215 ierr = PetscFree(supportRef);CHKERRQ(ierr); 1216 break; 1217 case 5: 1218 /* Simplicial 3D */ 1219 /* All cells have 4 faces: Tet face order is prescribed in DMPlexGetFaces_Internal() */ 1220 ierr = DMPlexGetRawFaces_Internal(dm, 3, 4, cellInd, NULL, NULL, &faces);CHKERRQ(ierr); 1221 for (c = cStart; c < cEnd; ++c) { 1222 const PetscInt newp = cStartNew + (c - cStart)*8; 1223 const PetscInt *cone, *ornt; 1224 PetscInt coneNew[4], orntNew[4]; 1225 1226 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1227 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1228 /* A tetrahedron: {0, a, c, d} */ 1229 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + (ornt[0] < 0 ? (-(ornt[0]+1)+1)%3 : ornt[0]); /* A */ 1230 orntNew[0] = ornt[0]; 1231 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + (ornt[1] < 0 ? (-(ornt[1]+1)+1)%3 : ornt[1]); /* A */ 1232 orntNew[1] = ornt[1]; 1233 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + (ornt[2] < 0 ? (-(ornt[2]+1)+1)%3 : ornt[2]); /* A */ 1234 orntNew[2] = ornt[2]; 1235 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 0; 1236 orntNew[3] = 0; 1237 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1238 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1239 #if 1 1240 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); 1241 for (p = 0; p < 4; ++p) { 1242 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); 1243 } 1244 #endif 1245 /* B tetrahedron: {a, 1, b, e} */ 1246 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + (ornt[0] < 0 ? (-(ornt[0]+1)+0)%3 : (ornt[0]+1)%3); /* B */ 1247 orntNew[0] = ornt[0]; 1248 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + (ornt[1] < 0 ? (-(ornt[1]+1)+2)%3 : (ornt[1]+2)%3); /* C */ 1249 orntNew[1] = ornt[1]; 1250 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 1; 1251 orntNew[2] = 0; 1252 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + (ornt[3] < 0 ? (-(ornt[3]+1)+0)%3 : (ornt[3]+1)%3); /* B */ 1253 orntNew[3] = ornt[3]; 1254 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1255 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1256 #if 1 1257 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); 1258 for (p = 0; p < 4; ++p) { 1259 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); 1260 } 1261 #endif 1262 /* C tetrahedron: {c, b, 2, f} */ 1263 coneNew[0] = fStartNew + (cone[0] - fStart)*4 + (ornt[0] < 0 ? (-(ornt[0]+1)+2)%3 : (ornt[0]+2)%3); /* C */ 1264 orntNew[0] = ornt[0]; 1265 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 2; 1266 orntNew[1] = 0; 1267 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + (ornt[2] < 0 ? (-(ornt[2]+1)+0)%3 : (ornt[2]+1)%3); /* B */ 1268 orntNew[2] = ornt[2]; 1269 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + (ornt[3] < 0 ? (-(ornt[3]+1)+1)%3 : (ornt[3]+0)%3); /* A */ 1270 orntNew[3] = ornt[3]; 1271 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 1272 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 1273 #if 1 1274 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); 1275 for (p = 0; p < 4; ++p) { 1276 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); 1277 } 1278 #endif 1279 /* D tetrahedron: {d, e, f, 3} */ 1280 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 3; 1281 orntNew[0] = 0; 1282 coneNew[1] = fStartNew + (cone[1] - fStart)*4 + (ornt[1] < 0 ? (-(ornt[1]+1)+0)%3 : (ornt[1]+1)%3); /* B */ 1283 orntNew[1] = ornt[1]; 1284 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + (ornt[2] < 0 ? (-(ornt[2]+1)+2)%3 : (ornt[2]+2)%3); /* C */ 1285 orntNew[2] = ornt[2]; 1286 coneNew[3] = fStartNew + (cone[3] - fStart)*4 + (ornt[3] < 0 ? (-(ornt[3]+1)+2)%3 : (ornt[3]+2)%3); /* C */ 1287 orntNew[3] = ornt[3]; 1288 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 1289 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 1290 #if 1 1291 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); 1292 for (p = 0; p < 4; ++p) { 1293 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); 1294 } 1295 #endif 1296 /* A' tetrahedron: {d, a, c, f} */ 1297 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 0; 1298 orntNew[0] = -3; 1299 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 4; 1300 orntNew[1] = 0; 1301 coneNew[2] = fStartNew + (cone[2] - fStart)*4 + 3; 1302 orntNew[2] = ornt[2] < 0 ? -((-(ornt[2]+1)+2)%3+1) : (ornt[2]+2)%3; 1303 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 5; 1304 orntNew[3] = 0; 1305 ierr = DMPlexSetCone(rdm, newp+4, coneNew);CHKERRQ(ierr); 1306 ierr = DMPlexSetConeOrientation(rdm, newp+4, orntNew);CHKERRQ(ierr); 1307 #if 1 1308 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); 1309 for (p = 0; p < 4; ++p) { 1310 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); 1311 } 1312 #endif 1313 /* B' tetrahedron: {e, b, a, f} */ 1314 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 1; 1315 orntNew[0] = -3; 1316 coneNew[1] = fStartNew + (cone[3] - fStart)*4 + 3; 1317 orntNew[1] = ornt[3] < 0 ? -((-(ornt[3]+1)+1)%3+1) : (ornt[3]+1)%3; 1318 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 6; 1319 orntNew[2] = 0; 1320 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 7; 1321 orntNew[3] = 0; 1322 ierr = DMPlexSetCone(rdm, newp+5, coneNew);CHKERRQ(ierr); 1323 ierr = DMPlexSetConeOrientation(rdm, newp+5, orntNew);CHKERRQ(ierr); 1324 #if 1 1325 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); 1326 for (p = 0; p < 4; ++p) { 1327 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); 1328 } 1329 #endif 1330 /* C' tetrahedron: {b, f, c, a} */ 1331 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 2; 1332 orntNew[0] = -3; 1333 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 7; 1334 orntNew[1] = -2; 1335 coneNew[2] = fStartNew + (cone[0] - fStart)*4 + 3; 1336 orntNew[2] = ornt[0] < 0 ? (-(ornt[0]+1)+1)%3 : -((ornt[0]+1)%3+1); 1337 coneNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 5; 1338 orntNew[3] = -1; 1339 ierr = DMPlexSetCone(rdm, newp+6, coneNew);CHKERRQ(ierr); 1340 ierr = DMPlexSetConeOrientation(rdm, newp+6, orntNew);CHKERRQ(ierr); 1341 #if 1 1342 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); 1343 for (p = 0; p < 4; ++p) { 1344 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); 1345 } 1346 #endif 1347 /* D' tetrahedron: {f, e, d, a} */ 1348 coneNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 3; 1349 orntNew[0] = -3; 1350 coneNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 6; 1351 orntNew[1] = -3; 1352 coneNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 4; 1353 orntNew[2] = -2; 1354 coneNew[3] = fStartNew + (cone[1] - fStart)*4 + 3; 1355 orntNew[3] = ornt[2]; 1356 ierr = DMPlexSetCone(rdm, newp+7, coneNew);CHKERRQ(ierr); 1357 ierr = DMPlexSetConeOrientation(rdm, newp+7, orntNew);CHKERRQ(ierr); 1358 #if 1 1359 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); 1360 for (p = 0; p < 4; ++p) { 1361 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); 1362 } 1363 #endif 1364 } 1365 /* Split faces have 3 edges and the same cells as the parent */ 1366 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 1367 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 1368 for (f = fStart; f < fEnd; ++f) { 1369 const PetscInt newp = fStartNew + (f - fStart)*4; 1370 const PetscInt *cone, *ornt, *support; 1371 PetscInt coneNew[3], orntNew[3], coneSize, supportSize, s; 1372 1373 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1374 ierr = DMPlexGetConeOrientation(dm, f, &ornt);CHKERRQ(ierr); 1375 /* A triangle */ 1376 coneNew[0] = eStartNew + (cone[0] - eStart)*2 + (ornt[0] < 0 ? 1 : 0); 1377 orntNew[0] = ornt[0]; 1378 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 2; 1379 orntNew[1] = -2; 1380 coneNew[2] = eStartNew + (cone[2] - eStart)*2 + (ornt[2] < 0 ? 0 : 1); 1381 orntNew[2] = ornt[2]; 1382 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 1383 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 1384 #if 1 1385 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); 1386 for (p = 0; p < 3; ++p) { 1387 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); 1388 } 1389 #endif 1390 /* B triangle */ 1391 coneNew[0] = eStartNew + (cone[0] - eStart)*2 + (ornt[0] < 0 ? 0 : 1); 1392 orntNew[0] = ornt[0]; 1393 coneNew[1] = eStartNew + (cone[1] - eStart)*2 + (ornt[1] < 0 ? 1 : 0); 1394 orntNew[1] = ornt[1]; 1395 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 0; 1396 orntNew[2] = -2; 1397 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 1398 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 1399 #if 1 1400 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); 1401 for (p = 0; p < 3; ++p) { 1402 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); 1403 } 1404 #endif 1405 /* C triangle */ 1406 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 1; 1407 orntNew[0] = -2; 1408 coneNew[1] = eStartNew + (cone[1] - eStart)*2 + (ornt[1] < 0 ? 0 : 1); 1409 orntNew[1] = ornt[1]; 1410 coneNew[2] = eStartNew + (cone[2] - eStart)*2 + (ornt[2] < 0 ? 1 : 0); 1411 orntNew[2] = ornt[2]; 1412 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 1413 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 1414 #if 1 1415 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); 1416 for (p = 0; p < 3; ++p) { 1417 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); 1418 } 1419 #endif 1420 /* D triangle */ 1421 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 0; 1422 orntNew[0] = 0; 1423 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 1; 1424 orntNew[1] = 0; 1425 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + 2; 1426 orntNew[2] = 0; 1427 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 1428 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 1429 #if 1 1430 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); 1431 for (p = 0; p < 3; ++p) { 1432 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); 1433 } 1434 #endif 1435 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1436 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1437 for (r = 0; r < 4; ++r) { 1438 for (s = 0; s < supportSize; ++s) { 1439 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1440 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1441 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1442 for (c = 0; c < coneSize; ++c) { 1443 if (cone[c] == f) break; 1444 } 1445 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])); 1446 } 1447 ierr = DMPlexSetSupport(rdm, newp+r, supportRef);CHKERRQ(ierr); 1448 #if 1 1449 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1450 for (p = 0; p < supportSize; ++p) { 1451 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); 1452 } 1453 #endif 1454 } 1455 } 1456 /* Interior faces have 3 edges and 2 cells */ 1457 for (c = cStart; c < cEnd; ++c) { 1458 PetscInt newp = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8; 1459 const PetscInt *cone, *ornt; 1460 PetscInt coneNew[3], orntNew[3]; 1461 PetscInt supportNew[2]; 1462 1463 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1464 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1465 /* Face A: {c, a, d} */ 1466 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+0)%3 : (ornt[0]+2)%3); 1467 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1468 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+0)%3 : (ornt[1]+2)%3); 1469 orntNew[1] = ornt[1] < 0 ? -2 : 0; 1470 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+0)%3 : (ornt[2]+2)%3); 1471 orntNew[2] = ornt[2] < 0 ? -2 : 0; 1472 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1473 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1474 #if 1 1475 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1476 for (p = 0; p < 3; ++p) { 1477 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); 1478 } 1479 #endif 1480 supportNew[0] = (c - cStart)*8 + 0; 1481 supportNew[1] = (c - cStart)*8 + 0+4; 1482 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1483 #if 1 1484 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1485 for (p = 0; p < 2; ++p) { 1486 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); 1487 } 1488 #endif 1489 ++newp; 1490 /* Face B: {a, b, e} */ 1491 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+2)%3 : (ornt[0]+0)%3); 1492 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1493 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+2)%3 : (ornt[3]+0)%3); 1494 orntNew[1] = ornt[3] < 0 ? -2 : 0; 1495 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+1)%3 : (ornt[1]+1)%3); 1496 orntNew[2] = ornt[1] < 0 ? -2 : 0; 1497 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1498 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1499 #if 1 1500 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); 1501 for (p = 0; p < 3; ++p) { 1502 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); 1503 } 1504 #endif 1505 supportNew[0] = (c - cStart)*8 + 1; 1506 supportNew[1] = (c - cStart)*8 + 1+4; 1507 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1508 #if 1 1509 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1510 for (p = 0; p < 2; ++p) { 1511 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); 1512 } 1513 #endif 1514 ++newp; 1515 /* Face C: {c, f, b} */ 1516 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+2)%3 : (ornt[2]+0)%3); 1517 orntNew[0] = ornt[2] < 0 ? -2 : 0; 1518 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+0)%3 : (ornt[3]+2)%3); 1519 orntNew[1] = ornt[3] < 0 ? -2 : 0; 1520 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+1)%3 : (ornt[0]+1)%3); 1521 orntNew[2] = ornt[0] < 0 ? -2 : 0; 1522 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1523 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1524 #if 1 1525 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1526 for (p = 0; p < 3; ++p) { 1527 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); 1528 } 1529 #endif 1530 supportNew[0] = (c - cStart)*8 + 2; 1531 supportNew[1] = (c - cStart)*8 + 2+4; 1532 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1533 #if 1 1534 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1535 for (p = 0; p < 2; ++p) { 1536 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); 1537 } 1538 #endif 1539 ++newp; 1540 /* Face D: {d, e, f} */ 1541 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+2)%3 : (ornt[1]+0)%3); 1542 orntNew[0] = ornt[1] < 0 ? -2 : 0; 1543 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+1)%3 : (ornt[3]+1)%3); 1544 orntNew[1] = ornt[3] < 0 ? -2 : 0; 1545 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+1)%3 : (ornt[2]+1)%3); 1546 orntNew[2] = ornt[2] < 0 ? -2 : 0; 1547 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1548 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1549 #if 1 1550 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1551 for (p = 0; p < 3; ++p) { 1552 if ((coneNew[p] < eStartNew) || (coneNew[p] >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", coneNew[p], eStartNew, eEndNew); 1553 } 1554 #endif 1555 supportNew[0] = (c - cStart)*8 + 3; 1556 supportNew[1] = (c - cStart)*8 + 3+4; 1557 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1558 #if 1 1559 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1560 for (p = 0; p < 2; ++p) { 1561 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); 1562 } 1563 #endif 1564 ++newp; 1565 /* Face E: {d, f, a} */ 1566 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+1)%3 : (ornt[2]+1)%3); 1567 orntNew[0] = ornt[2] < 0 ? 0 : -2; 1568 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1569 orntNew[1] = 0; 1570 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+0)%3 : (ornt[1]+2)%3); 1571 orntNew[2] = ornt[1] < 0 ? -2 : 0; 1572 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1573 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1574 #if 1 1575 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1576 for (p = 0; p < 3; ++p) { 1577 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); 1578 } 1579 #endif 1580 supportNew[0] = (c - cStart)*8 + 0+4; 1581 supportNew[1] = (c - cStart)*8 + 3+4; 1582 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1583 #if 1 1584 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1585 for (p = 0; p < 2; ++p) { 1586 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); 1587 } 1588 #endif 1589 ++newp; 1590 /* Face F: {c, a, f} */ 1591 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+0)%3 : (ornt[0]+2)%3); 1592 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1593 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1594 orntNew[1] = -2; 1595 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[2] - fStart)*3 + (ornt[2] < 0 ? (-(ornt[2]+1)+2)%3 : (ornt[2]+0)%3); 1596 orntNew[2] = ornt[1] < 0 ? 0 : -2; 1597 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1598 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1599 #if 1 1600 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1601 for (p = 0; p < 3; ++p) { 1602 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); 1603 } 1604 #endif 1605 supportNew[0] = (c - cStart)*8 + 0+4; 1606 supportNew[1] = (c - cStart)*8 + 2+4; 1607 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1608 #if 1 1609 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1610 for (p = 0; p < 2; ++p) { 1611 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); 1612 } 1613 #endif 1614 ++newp; 1615 /* Face G: {e, a, f} */ 1616 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[1] - fStart)*3 + (ornt[1] < 0 ? (-(ornt[1]+1)+1)%3 : (ornt[1]+1)%3); 1617 orntNew[0] = ornt[1] < 0 ? -2 : 0; 1618 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1619 orntNew[1] = -2; 1620 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+1)%3 : (ornt[3]+1)%3); 1621 orntNew[2] = ornt[3] < 0 ? 0 : -2; 1622 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1623 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1624 #if 1 1625 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1626 for (p = 0; p < 3; ++p) { 1627 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); 1628 } 1629 #endif 1630 supportNew[0] = (c - cStart)*8 + 1+4; 1631 supportNew[1] = (c - cStart)*8 + 3+4; 1632 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1633 #if 1 1634 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1635 for (p = 0; p < 2; ++p) { 1636 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); 1637 } 1638 #endif 1639 ++newp; 1640 /* Face H: {a, b, f} */ 1641 coneNew[0] = eStartNew + (eEnd - eStart)*2 + (cone[0] - fStart)*3 + (ornt[0] < 0 ? (-(ornt[0]+1)+2)%3 : (ornt[0]+0)%3); 1642 orntNew[0] = ornt[0] < 0 ? -2 : 0; 1643 coneNew[1] = eStartNew + (eEnd - eStart)*2 + (cone[3] - fStart)*3 + (ornt[3] < 0 ? (-(ornt[3]+1)+0)%3 : (ornt[3]+2)%3); 1644 orntNew[1] = ornt[3] < 0 ? 0 : -2; 1645 coneNew[2] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1646 orntNew[2] = 0; 1647 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1648 ierr = DMPlexSetConeOrientation(rdm, newp, orntNew);CHKERRQ(ierr); 1649 #if 1 1650 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1651 for (p = 0; p < 3; ++p) { 1652 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); 1653 } 1654 #endif 1655 supportNew[0] = (c - cStart)*8 + 1+4; 1656 supportNew[1] = (c - cStart)*8 + 2+4; 1657 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1658 #if 1 1659 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 1660 for (p = 0; p < 2; ++p) { 1661 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); 1662 } 1663 #endif 1664 ++newp; 1665 } 1666 /* Split Edges have 2 vertices and the same faces as the parent */ 1667 for (e = eStart; e < eEnd; ++e) { 1668 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 1669 1670 for (r = 0; r < 2; ++r) { 1671 const PetscInt newp = eStartNew + (e - eStart)*2 + r; 1672 const PetscInt *cone, *ornt, *support; 1673 PetscInt coneNew[2], coneSize, c, supportSize, s; 1674 1675 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 1676 coneNew[0] = vStartNew + (cone[0] - vStart); 1677 coneNew[1] = vStartNew + (cone[1] - vStart); 1678 coneNew[(r+1)%2] = newv; 1679 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1680 #if 1 1681 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1682 for (p = 0; p < 2; ++p) { 1683 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); 1684 } 1685 #endif 1686 ierr = DMPlexGetSupportSize(dm, e, &supportSize);CHKERRQ(ierr); 1687 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 1688 for (s = 0; s < supportSize; ++s) { 1689 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1690 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1691 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1692 for (c = 0; c < coneSize; ++c) { 1693 if (cone[c] == e) break; 1694 } 1695 supportRef[s] = fStartNew + (support[s] - fStart)*4 + (c + (ornt[c] < 0 ? 1-r : r))%3; 1696 } 1697 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1698 #if 1 1699 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1700 for (p = 0; p < supportSize; ++p) { 1701 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); 1702 } 1703 #endif 1704 } 1705 } 1706 /* Face edges have 2 vertices and 2+cells*(1/2) faces */ 1707 for (f = fStart; f < fEnd; ++f) { 1708 const PetscInt *cone, *ornt, *support; 1709 PetscInt coneSize, supportSize, s; 1710 1711 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 1712 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1713 for (r = 0; r < 3; ++r) { 1714 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (f - fStart)*3 + r; 1715 PetscInt coneNew[2], intFaces = 0, er, eint[4] = {1, 0, 2, 0}; 1716 PetscInt fint[24] = { 1, 7, -1, -1, 0, 5, 1717 -1, -1, 1, 6, 0, 4, 1718 2, 5, 3, 4, -1, -1, 1719 -1, -1, 3, 6, 2, 7}; 1720 1721 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1722 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[(r+0)%3] - eStart); 1723 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - eStart); 1724 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1725 #if 1 1726 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1727 for (p = 0; p < 2; ++p) { 1728 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); 1729 } 1730 #endif 1731 supportRef[0] = fStartNew + (f - fStart)*4 + (r+1)%3; 1732 supportRef[1] = fStartNew + (f - fStart)*4 + 3; 1733 for (s = 0; s < supportSize; ++s) { 1734 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1735 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1736 ierr = DMPlexGetConeOrientation(dm, support[s], &ornt);CHKERRQ(ierr); 1737 for (c = 0; c < coneSize; ++c) {if (cone[c] == f) break;} 1738 /* Here we want to determine whether edge newp contains a vertex which is part of the cross-tet edge */ 1739 er = ornt[c] < 0 ? (-(ornt[c]+1) + 2-r)%3 : (ornt[c] + r)%3; 1740 if (er == eint[c]) { 1741 supportRef[2+intFaces++] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*8 + (c + 2)%4; 1742 } else { 1743 supportRef[2+intFaces++] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*8 + fint[(c*3 + er)*2 + 0]; 1744 supportRef[2+intFaces++] = fStartNew + (fEnd - fStart)*4 + (support[s] - cStart)*8 + fint[(c*3 + er)*2 + 1]; 1745 } 1746 } 1747 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1748 #if 1 1749 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1750 for (p = 0; p < intFaces; ++p) { 1751 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); 1752 } 1753 #endif 1754 } 1755 } 1756 /* Interior edges have 2 vertices and 4 faces */ 1757 for (c = cStart; c < cEnd; ++c) { 1758 const PetscInt newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (c - cStart); 1759 const PetscInt *cone, *ornt, *fcone; 1760 PetscInt coneNew[2], supportNew[4], find; 1761 1762 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 1763 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 1764 ierr = DMPlexGetCone(dm, cone[0], &fcone);CHKERRQ(ierr); 1765 find = ornt[0] < 0 ? (-(ornt[0]+1) + 0)%3 : ornt[0]; 1766 coneNew[0] = vStartNew + (vEnd - vStart) + (fcone[find] - eStart); 1767 ierr = DMPlexGetCone(dm, cone[2], &fcone);CHKERRQ(ierr); 1768 find = ornt[2] < 0 ? (-(ornt[2]+1) + 1)%3 : (ornt[2]+1)%3; 1769 coneNew[1] = vStartNew + (vEnd - vStart) + (fcone[find] - eStart); 1770 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 1771 #if 1 1772 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1773 for (p = 0; p < 2; ++p) { 1774 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); 1775 } 1776 #endif 1777 supportNew[0] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 4; 1778 supportNew[1] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 5; 1779 supportNew[2] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 6; 1780 supportNew[3] = fStartNew + (fEnd - fStart)*4 + (c - cStart)*8 + 7; 1781 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 1782 #if 1 1783 if ((newp < eStartNew) || (newp >= eEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not an edge [%d, %d)", newp, eStartNew, eEndNew); 1784 for (p = 0; p < 4; ++p) { 1785 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); 1786 } 1787 #endif 1788 } 1789 /* Old vertices have identical supports */ 1790 for (v = vStart; v < vEnd; ++v) { 1791 const PetscInt newp = vStartNew + (v - vStart); 1792 const PetscInt *support, *cone; 1793 PetscInt size, s; 1794 1795 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 1796 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 1797 for (s = 0; s < size; ++s) { 1798 PetscInt r = 0; 1799 1800 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1801 if (cone[1] == v) r = 1; 1802 supportRef[s] = eStartNew + (support[s] - eStart)*2 + r; 1803 } 1804 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1805 #if 1 1806 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1807 for (p = 0; p < size; ++p) { 1808 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); 1809 } 1810 #endif 1811 } 1812 /* Edge vertices have 2 + face*2 + 0/1 supports */ 1813 for (e = eStart; e < eEnd; ++e) { 1814 const PetscInt newp = vStartNew + (vEnd - vStart) + (e - eStart); 1815 const PetscInt *cone, *support; 1816 PetscInt *star = NULL, starSize, cellSize = 0, coneSize, size, s; 1817 1818 ierr = DMPlexGetSupportSize(dm, e, &size);CHKERRQ(ierr); 1819 ierr = DMPlexGetSupport(dm, e, &support);CHKERRQ(ierr); 1820 supportRef[0] = eStartNew + (e - eStart)*2 + 0; 1821 supportRef[1] = eStartNew + (e - eStart)*2 + 1; 1822 for (s = 0; s < size; ++s) { 1823 PetscInt r = 0; 1824 1825 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1826 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1827 for (r = 0; r < coneSize; ++r) {if (cone[r] == e) break;} 1828 supportRef[2+s*2+0] = eStartNew + (eEnd - eStart)*2 + (support[s] - fStart)*3 + (r+0)%3; 1829 supportRef[2+s*2+1] = eStartNew + (eEnd - eStart)*2 + (support[s] - fStart)*3 + (r+2)%3; 1830 } 1831 ierr = DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1832 for (s = 0; s < starSize*2; s += 2) { 1833 const PetscInt *cone, *ornt; 1834 PetscInt e01, e23; 1835 1836 if ((star[s] >= cStart) && (star[s] < cEnd)) { 1837 /* Check edge 0-1 */ 1838 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 1839 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 1840 ierr = DMPlexGetCone(dm, cone[0], &cone);CHKERRQ(ierr); 1841 e01 = cone[ornt[0] < 0 ? (-(ornt[0]+1) + 0)%3 : ornt[0]]; 1842 /* Check edge 2-3 */ 1843 ierr = DMPlexGetCone(dm, star[s], &cone);CHKERRQ(ierr); 1844 ierr = DMPlexGetConeOrientation(dm, star[s], &ornt);CHKERRQ(ierr); 1845 ierr = DMPlexGetCone(dm, cone[3], &cone);CHKERRQ(ierr); 1846 e23 = cone[ornt[3] < 0 ? (-(ornt[3]+1) + 2)%3 : (ornt[3] + 2)%3]; 1847 if ((e01 == e) || (e23 == e)) {supportRef[2+size*2+cellSize++] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (star[s] - cStart);} 1848 } 1849 } 1850 ierr = DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1851 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1852 #if 1 1853 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1854 for (p = 0; p < 2+size*2+cellSize; ++p) { 1855 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); 1856 } 1857 #endif 1858 } 1859 ierr = PetscFree(supportRef);CHKERRQ(ierr); 1860 ierr = DMPlexRestoreFaces_Internal(dm, 3, cStart, NULL, NULL, &faces);CHKERRQ(ierr); 1861 break; 1862 default: 1863 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 1864 } 1865 PetscFunctionReturn(0); 1866 } 1867 1868 #undef __FUNCT__ 1869 #define __FUNCT__ "CellRefinerSetCoordinates" 1870 PetscErrorCode CellRefinerSetCoordinates(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 1871 { 1872 PetscSection coordSection, coordSectionNew; 1873 Vec coordinates, coordinatesNew; 1874 PetscScalar *coords, *coordsNew; 1875 PetscInt dim, depth, coordSizeNew, cStart, cEnd, c, vStart, vStartNew, vEnd, v, eStart, eEnd, eMax, e, fStart, fEnd, fMax, f; 1876 PetscErrorCode ierr; 1877 1878 PetscFunctionBegin; 1879 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1880 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1881 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1882 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 1883 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1884 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1885 ierr = DMPlexGetHybridBounds(dm, NULL, &fMax, &eMax, NULL);CHKERRQ(ierr); 1886 ierr = GetDepthStart_Private(depth, depthSize, NULL, NULL, NULL, &vStartNew);CHKERRQ(ierr); 1887 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1888 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &coordSectionNew);CHKERRQ(ierr); 1889 ierr = PetscSectionSetNumFields(coordSectionNew, 1);CHKERRQ(ierr); 1890 ierr = PetscSectionSetFieldComponents(coordSectionNew, 0, dim);CHKERRQ(ierr); 1891 ierr = PetscSectionSetChart(coordSectionNew, vStartNew, vStartNew+depthSize[0]);CHKERRQ(ierr); 1892 if (fMax < 0) fMax = fEnd; 1893 if (eMax < 0) eMax = eEnd; 1894 /* All vertices have the dim coordinates */ 1895 for (v = vStartNew; v < vStartNew+depthSize[0]; ++v) { 1896 ierr = PetscSectionSetDof(coordSectionNew, v, dim);CHKERRQ(ierr); 1897 ierr = PetscSectionSetFieldDof(coordSectionNew, v, 0, dim);CHKERRQ(ierr); 1898 } 1899 ierr = PetscSectionSetUp(coordSectionNew);CHKERRQ(ierr); 1900 ierr = DMPlexSetCoordinateSection(rdm, coordSectionNew);CHKERRQ(ierr); 1901 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1902 ierr = PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);CHKERRQ(ierr); 1903 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinatesNew);CHKERRQ(ierr); 1904 ierr = PetscObjectSetName((PetscObject) coordinatesNew, "coordinates");CHKERRQ(ierr); 1905 ierr = VecSetSizes(coordinatesNew, coordSizeNew, PETSC_DETERMINE);CHKERRQ(ierr); 1906 ierr = VecSetFromOptions(coordinatesNew);CHKERRQ(ierr); 1907 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1908 ierr = VecGetArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 1909 switch (refiner) { 1910 case 6: /* Hex 3D */ 1911 /* Face vertices have the average of corner coordinates */ 1912 for (f = fStart; f < fEnd; ++f) { 1913 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (cEnd - cStart) + (f - fStart); 1914 PetscInt *cone = NULL; 1915 PetscInt closureSize, coneSize = 0, off[8], offnew, p, d; 1916 1917 ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 1918 for (p = 0; p < closureSize*2; p += 2) { 1919 const PetscInt point = cone[p]; 1920 if ((point >= vStart) && (point < vEnd)) cone[coneSize++] = point; 1921 } 1922 for (v = 0; v < coneSize; ++v) { 1923 ierr = PetscSectionGetOffset(coordSection, cone[v], &off[v]);CHKERRQ(ierr); 1924 } 1925 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 1926 for (d = 0; d < dim; ++d) { 1927 coordsNew[offnew+d] = 0.0; 1928 for (v = 0; v < coneSize; ++v) coordsNew[offnew+d] += coords[off[v]+d]; 1929 coordsNew[offnew+d] /= coneSize; 1930 } 1931 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 1932 } 1933 case 2: /* Hex 2D */ 1934 /* Cell vertices have the average of corner coordinates */ 1935 for (c = cStart; c < cEnd; ++c) { 1936 const PetscInt newv = vStartNew + (vEnd - vStart) + (eEnd - eStart) + (c - cStart); 1937 PetscInt *cone = NULL; 1938 PetscInt closureSize, coneSize = 0, off[8], offnew, p, d; 1939 1940 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 1941 for (p = 0; p < closureSize*2; p += 2) { 1942 const PetscInt point = cone[p]; 1943 if ((point >= vStart) && (point < vEnd)) cone[coneSize++] = point; 1944 } 1945 for (v = 0; v < coneSize; ++v) { 1946 ierr = PetscSectionGetOffset(coordSection, cone[v], &off[v]);CHKERRQ(ierr); 1947 } 1948 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 1949 for (d = 0; d < dim; ++d) { 1950 coordsNew[offnew+d] = 0.0; 1951 for (v = 0; v < coneSize; ++v) coordsNew[offnew+d] += coords[off[v]+d]; 1952 coordsNew[offnew+d] /= coneSize; 1953 } 1954 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 1955 } 1956 case 1: /* Simplicial 2D */ 1957 case 3: /* Hybrid Simplicial 2D */ 1958 case 5: /* Simplicial 3D */ 1959 /* Edge vertices have the average of endpoint coordinates */ 1960 for (e = eStart; e < eMax; ++e) { 1961 const PetscInt newv = vStartNew + (vEnd - vStart) + (e - eStart); 1962 const PetscInt *cone; 1963 PetscInt coneSize, offA, offB, offnew, d; 1964 1965 ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr); 1966 if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %d cone should have two vertices, not %d", e, coneSize); 1967 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 1968 ierr = PetscSectionGetOffset(coordSection, cone[0], &offA);CHKERRQ(ierr); 1969 ierr = PetscSectionGetOffset(coordSection, cone[1], &offB);CHKERRQ(ierr); 1970 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 1971 for (d = 0; d < dim; ++d) { 1972 coordsNew[offnew+d] = 0.5*(coords[offA+d] + coords[offB+d]); 1973 } 1974 } 1975 /* Old vertices have the same coordinates */ 1976 for (v = vStart; v < vEnd; ++v) { 1977 const PetscInt newv = vStartNew + (v - vStart); 1978 PetscInt off, offnew, d; 1979 1980 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1981 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 1982 for (d = 0; d < dim; ++d) { 1983 coordsNew[offnew+d] = coords[off+d]; 1984 } 1985 } 1986 break; 1987 default: 1988 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 1989 } 1990 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1991 ierr = VecRestoreArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 1992 ierr = DMSetCoordinatesLocal(rdm, coordinatesNew);CHKERRQ(ierr); 1993 ierr = VecDestroy(&coordinatesNew);CHKERRQ(ierr); 1994 ierr = PetscSectionDestroy(&coordSectionNew);CHKERRQ(ierr); 1995 PetscFunctionReturn(0); 1996 } 1997 1998 #undef __FUNCT__ 1999 #define __FUNCT__ "DMPlexCreateProcessSF" 2000 PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess) 2001 { 2002 PetscInt numRoots, numLeaves, l; 2003 const PetscInt *localPoints; 2004 const PetscSFNode *remotePoints; 2005 PetscInt *localPointsNew; 2006 PetscSFNode *remotePointsNew; 2007 PetscInt *ranks, *ranksNew; 2008 PetscErrorCode ierr; 2009 2010 PetscFunctionBegin; 2011 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2012 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranks);CHKERRQ(ierr); 2013 for (l = 0; l < numLeaves; ++l) { 2014 ranks[l] = remotePoints[l].rank; 2015 } 2016 ierr = PetscSortRemoveDupsInt(&numLeaves, ranks);CHKERRQ(ierr); 2017 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranksNew);CHKERRQ(ierr); 2018 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 2019 ierr = PetscMalloc(numLeaves * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 2020 for (l = 0; l < numLeaves; ++l) { 2021 ranksNew[l] = ranks[l]; 2022 localPointsNew[l] = l; 2023 remotePointsNew[l].index = 0; 2024 remotePointsNew[l].rank = ranksNew[l]; 2025 } 2026 ierr = PetscFree(ranks);CHKERRQ(ierr); 2027 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);CHKERRQ(ierr); 2028 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);CHKERRQ(ierr); 2029 ierr = PetscSFSetFromOptions(*sfProcess);CHKERRQ(ierr); 2030 ierr = PetscSFSetGraph(*sfProcess, 1, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2031 PetscFunctionReturn(0); 2032 } 2033 2034 #undef __FUNCT__ 2035 #define __FUNCT__ "CellRefinerCreateSF" 2036 PetscErrorCode CellRefinerCreateSF(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 2037 { 2038 PetscSF sf, sfNew, sfProcess; 2039 IS processRanks; 2040 MPI_Datatype depthType; 2041 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 2042 const PetscInt *localPoints, *neighbors; 2043 const PetscSFNode *remotePoints; 2044 PetscInt *localPointsNew; 2045 PetscSFNode *remotePointsNew; 2046 PetscInt *depthSizeOld, *rdepthSize, *rdepthSizeOld, *rdepthMaxOld, *rvStart, *rvStartNew, *reStart, *reStartNew, *rfStart, *rfStartNew, *rcStart, *rcStartNew; 2047 PetscInt depth, numNeighbors, pStartNew, pEndNew, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eStartNew, eEnd, eMax, r, n; 2048 PetscErrorCode ierr; 2049 2050 PetscFunctionBegin; 2051 ierr = DMPlexGetChart(rdm, &pStartNew, &pEndNew);CHKERRQ(ierr); 2052 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2053 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2054 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 2055 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2056 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2057 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 2058 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 2059 switch (refiner) { 2060 case 3: 2061 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 2062 cMax = PetscMin(cEnd, cMax); 2063 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 2064 fMax = PetscMin(fEnd, fMax); 2065 } 2066 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2067 ierr = DMGetPointSF(rdm, &sfNew);CHKERRQ(ierr); 2068 /* Caculate size of new SF */ 2069 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2070 if (numRoots < 0) PetscFunctionReturn(0); 2071 for (l = 0; l < numLeaves; ++l) { 2072 const PetscInt p = localPoints[l]; 2073 2074 switch (refiner) { 2075 case 1: 2076 /* Simplicial 2D */ 2077 if ((p >= vStart) && (p < vEnd)) { 2078 /* Old vertices stay the same */ 2079 ++numLeavesNew; 2080 } else if ((p >= fStart) && (p < fEnd)) { 2081 /* Old faces add new faces and vertex */ 2082 numLeavesNew += 2 + 1; 2083 } else if ((p >= cStart) && (p < cEnd)) { 2084 /* Old cells add new cells and interior faces */ 2085 numLeavesNew += 4 + 3; 2086 } 2087 break; 2088 case 2: 2089 /* Hex 2D */ 2090 if ((p >= vStart) && (p < vEnd)) { 2091 /* Old vertices stay the same */ 2092 ++numLeavesNew; 2093 } else if ((p >= fStart) && (p < fEnd)) { 2094 /* Old faces add new faces and vertex */ 2095 numLeavesNew += 2 + 1; 2096 } else if ((p >= cStart) && (p < cEnd)) { 2097 /* Old cells add new cells, interior faces, and vertex */ 2098 numLeavesNew += 4 + 4 + 1; 2099 } 2100 break; 2101 case 5: 2102 /* Simplicial 3D */ 2103 if ((p >= vStart) && (p < vEnd)) { 2104 /* Old vertices stay the same */ 2105 ++numLeavesNew; 2106 } else if ((p >= eStart) && (p < eEnd)) { 2107 /* Old edges add new edges and vertex */ 2108 numLeavesNew += 2 + 1; 2109 } else if ((p >= fStart) && (p < fEnd)) { 2110 /* Old faces add new faces and face edges */ 2111 numLeavesNew += 4 + 3; 2112 } else if ((p >= cStart) && (p < cEnd)) { 2113 /* Old cells add new cells and interior faces and edges */ 2114 numLeavesNew += 8 + 8 + 1; 2115 } 2116 break; 2117 default: 2118 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2119 } 2120 } 2121 /* Communicate depthSizes for each remote rank */ 2122 ierr = DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);CHKERRQ(ierr); 2123 ierr = ISGetLocalSize(processRanks, &numNeighbors);CHKERRQ(ierr); 2124 ierr = PetscMalloc5((depth+1)*numNeighbors,PetscInt,&rdepthSize,numNeighbors,PetscInt,&rvStartNew,numNeighbors,PetscInt,&reStartNew,numNeighbors,PetscInt,&rfStartNew,numNeighbors,PetscInt,&rcStartNew);CHKERRQ(ierr); 2125 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); 2126 ierr = MPI_Type_contiguous(depth+1, MPIU_INT, &depthType);CHKERRQ(ierr); 2127 ierr = MPI_Type_commit(&depthType);CHKERRQ(ierr); 2128 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 2129 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 2130 for (n = 0; n < numNeighbors; ++n) { 2131 ierr = GetDepthStart_Private(depth, &rdepthSize[n*(depth+1)], &rcStartNew[n], &rfStartNew[n], &reStartNew[n], &rvStartNew[n]);CHKERRQ(ierr); 2132 } 2133 depthSizeOld[depth] = cMax; 2134 depthSizeOld[0] = vMax; 2135 depthSizeOld[depth-1] = fMax; 2136 depthSizeOld[1] = eMax; 2137 2138 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 2139 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 2140 2141 depthSizeOld[depth] = cEnd - cStart; 2142 depthSizeOld[0] = vEnd - vStart; 2143 depthSizeOld[depth-1] = fEnd - fStart; 2144 depthSizeOld[1] = eEnd - eStart; 2145 2146 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 2147 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 2148 for (n = 0; n < numNeighbors; ++n) { 2149 ierr = GetDepthStart_Private(depth, &rdepthSizeOld[n*(depth+1)], &rcStart[n], &rfStart[n], &reStart[n], &rvStart[n]);CHKERRQ(ierr); 2150 } 2151 ierr = MPI_Type_free(&depthType);CHKERRQ(ierr); 2152 ierr = PetscSFDestroy(&sfProcess);CHKERRQ(ierr); 2153 /* Calculate new point SF */ 2154 ierr = PetscMalloc(numLeavesNew * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 2155 ierr = PetscMalloc(numLeavesNew * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 2156 ierr = ISGetIndices(processRanks, &neighbors);CHKERRQ(ierr); 2157 for (l = 0, m = 0; l < numLeaves; ++l) { 2158 PetscInt p = localPoints[l]; 2159 PetscInt rp = remotePoints[l].index, n; 2160 PetscMPIInt rrank = remotePoints[l].rank; 2161 2162 ierr = PetscFindInt(rrank, numNeighbors, neighbors, &n);CHKERRQ(ierr); 2163 if (n < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %d", rrank); 2164 switch (refiner) { 2165 case 1: 2166 /* Simplicial 2D */ 2167 if ((p >= vStart) && (p < vEnd)) { 2168 /* Old vertices stay the same */ 2169 localPointsNew[m] = vStartNew + (p - vStart); 2170 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2171 remotePointsNew[m].rank = rrank; 2172 ++m; 2173 } else if ((p >= fStart) && (p < fEnd)) { 2174 /* Old faces add new faces and vertex */ 2175 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 2176 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 2177 remotePointsNew[m].rank = rrank; 2178 ++m; 2179 for (r = 0; r < 2; ++r, ++m) { 2180 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 2181 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 2182 remotePointsNew[m].rank = rrank; 2183 } 2184 } else if ((p >= cStart) && (p < cEnd)) { 2185 /* Old cells add new cells and interior faces */ 2186 for (r = 0; r < 4; ++r, ++m) { 2187 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2188 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2189 remotePointsNew[m].rank = rrank; 2190 } 2191 for (r = 0; r < 3; ++r, ++m) { 2192 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 2193 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*3 + r; 2194 remotePointsNew[m].rank = rrank; 2195 } 2196 } 2197 break; 2198 case 2: 2199 /* Hex 2D */ 2200 if ((p >= vStart) && (p < vEnd)) { 2201 /* Old vertices stay the same */ 2202 localPointsNew[m] = vStartNew + (p - vStart); 2203 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2204 remotePointsNew[m].rank = rrank; 2205 ++m; 2206 } else if ((p >= fStart) && (p < fEnd)) { 2207 /* Old faces add new faces and vertex */ 2208 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 2209 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 2210 remotePointsNew[m].rank = rrank; 2211 ++m; 2212 for (r = 0; r < 2; ++r, ++m) { 2213 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 2214 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 2215 remotePointsNew[m].rank = rrank; 2216 } 2217 } else if ((p >= cStart) && (p < cEnd)) { 2218 /* Old cells add new cells, interior faces, and vertex */ 2219 for (r = 0; r < 4; ++r, ++m) { 2220 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2221 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2222 remotePointsNew[m].rank = rrank; 2223 } 2224 for (r = 0; r < 4; ++r, ++m) { 2225 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 2226 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*4 + r; 2227 remotePointsNew[m].rank = rrank; 2228 } 2229 for (r = 0; r < 1; ++r, ++m) { 2230 localPointsNew[m] = vStartNew + (fEnd - fStart) + (p - cStart) + r; 2231 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1] + (rp - rcStart[n]) + r; 2232 remotePointsNew[m].rank = rrank; 2233 } 2234 } 2235 break; 2236 case 3: 2237 /* Hybrid simplicial 2D */ 2238 if ((p >= vStart) && (p < vEnd)) { 2239 /* Old vertices stay the same */ 2240 localPointsNew[m] = vStartNew + (p - vStart); 2241 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2242 remotePointsNew[m].rank = rrank; 2243 ++m; 2244 } else if ((p >= fStart) && (p < fMax)) { 2245 /* Old interior faces add new faces and vertex */ 2246 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 2247 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 2248 remotePointsNew[m].rank = rrank; 2249 ++m; 2250 for (r = 0; r < 2; ++r, ++m) { 2251 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 2252 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 2253 remotePointsNew[m].rank = rrank; 2254 } 2255 } else if ((p >= fMax) && (p < fEnd)) { 2256 /* Old hybrid faces stay the same */ 2257 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - fMax); 2258 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rdepthMaxOld[n*(depth+1)+depth-1]); 2259 remotePointsNew[m].rank = rrank; 2260 ++m; 2261 } else if ((p >= cStart) && (p < cMax)) { 2262 /* Old interior cells add new cells and interior faces */ 2263 for (r = 0; r < 4; ++r, ++m) { 2264 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2265 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2266 remotePointsNew[m].rank = rrank; 2267 } 2268 for (r = 0; r < 3; ++r, ++m) { 2269 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - cStart)*3 + r; 2270 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rcStart[n])*3 + r; 2271 remotePointsNew[m].rank = rrank; 2272 } 2273 } else if ((p >= cStart) && (p < cMax)) { 2274 /* Old hybrid cells add new cells and hybrid face */ 2275 for (r = 0; r < 2; ++r, ++m) { 2276 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 2277 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 2278 remotePointsNew[m].rank = rrank; 2279 } 2280 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 2281 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]); 2282 remotePointsNew[m].rank = rrank; 2283 ++m; 2284 } 2285 break; 2286 case 5: 2287 /* Simplicial 3D */ 2288 if ((p >= vStart) && (p < vEnd)) { 2289 /* Old vertices stay the same */ 2290 localPointsNew[m] = vStartNew + (p - vStart); 2291 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 2292 remotePointsNew[m].rank = rrank; 2293 ++m; 2294 } else if ((p >= eStart) && (p < eEnd)) { 2295 /* Old edges add new edges and vertex */ 2296 for (r = 0; r < 2; ++r, ++m) { 2297 localPointsNew[m] = eStartNew + (p - eStart)*2 + r; 2298 remotePointsNew[m].index = reStartNew[n] + (rp - reStart[n])*2 + r; 2299 remotePointsNew[m].rank = rrank; 2300 } 2301 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - eStart); 2302 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - reStart[n]); 2303 remotePointsNew[m].rank = rrank; 2304 ++m; 2305 } else if ((p >= fStart) && (p < fEnd)) { 2306 /* Old faces add new faces and face edges */ 2307 for (r = 0; r < 4; ++r, ++m) { 2308 localPointsNew[m] = fStartNew + (p - fStart)*4 + r; 2309 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*4 + r; 2310 remotePointsNew[m].rank = rrank; 2311 } 2312 for (r = 0; r < 3; ++r, ++m) { 2313 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (p - fStart)*3 + r; 2314 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + (rp - rfStart[n])*3 + r; 2315 remotePointsNew[m].rank = rrank; 2316 } 2317 } else if ((p >= cStart) && (p < cEnd)) { 2318 /* Old cells add new cells and interior faces and edges */ 2319 for (r = 0; r < 8; ++r, ++m) { 2320 localPointsNew[m] = cStartNew + (p - cStart)*8 + r; 2321 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*8 + r; 2322 remotePointsNew[m].rank = rrank; 2323 } 2324 for (r = 0; r < 8; ++r, ++m) { 2325 localPointsNew[m] = fStartNew + (fEnd - fStart)*4 + (p - cStart)*8 + r; 2326 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*4 + (rp - rcStart[n])*8 + r; 2327 remotePointsNew[m].rank = rrank; 2328 } 2329 for (r = 0; r < 1; ++r, ++m) { 2330 localPointsNew[m] = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (p - cStart)*1 + r; 2331 remotePointsNew[m].index = reStartNew[n] + rdepthSizeOld[n*(depth+1)+1]*2 + rdepthSizeOld[n*(depth+1)+depth-1]*3 + (rp - rcStart[n])*1 + r; 2332 remotePointsNew[m].rank = rrank; 2333 } 2334 } 2335 break; 2336 default: 2337 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2338 } 2339 } 2340 ierr = ISRestoreIndices(processRanks, &neighbors);CHKERRQ(ierr); 2341 ierr = ISDestroy(&processRanks);CHKERRQ(ierr); 2342 ierr = PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 2343 ierr = PetscFree5(rdepthSize,rvStartNew,reStartNew,rfStartNew,rcStartNew);CHKERRQ(ierr); 2344 ierr = PetscFree6(depthSizeOld,rdepthSizeOld,rvStart,reStart,rfStart,rcStart);CHKERRQ(ierr); 2345 PetscFunctionReturn(0); 2346 } 2347 2348 #undef __FUNCT__ 2349 #define __FUNCT__ "CellRefinerCreateLabels" 2350 PetscErrorCode CellRefinerCreateLabels(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 2351 { 2352 PetscInt numLabels, l; 2353 PetscInt depth, newp, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eStartNew, eEnd, eMax, r; 2354 PetscErrorCode ierr; 2355 2356 PetscFunctionBegin; 2357 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2358 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 2359 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2360 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2361 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2362 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 2363 ierr = DMPlexGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 2364 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 2365 switch (refiner) { 2366 case 3: 2367 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 2368 cMax = PetscMin(cEnd, cMax); 2369 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 2370 fMax = PetscMin(fEnd, fMax); 2371 } 2372 for (l = 0; l < numLabels; ++l) { 2373 DMLabel label, labelNew; 2374 const char *lname; 2375 PetscBool isDepth; 2376 IS valueIS; 2377 const PetscInt *values; 2378 PetscInt numValues, val; 2379 2380 ierr = DMPlexGetLabelName(dm, l, &lname);CHKERRQ(ierr); 2381 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 2382 if (isDepth) continue; 2383 ierr = DMPlexCreateLabel(rdm, lname);CHKERRQ(ierr); 2384 ierr = DMPlexGetLabel(dm, lname, &label);CHKERRQ(ierr); 2385 ierr = DMPlexGetLabel(rdm, lname, &labelNew);CHKERRQ(ierr); 2386 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 2387 ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr); 2388 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 2389 for (val = 0; val < numValues; ++val) { 2390 IS pointIS; 2391 const PetscInt *points; 2392 PetscInt numPoints, n; 2393 2394 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 2395 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 2396 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2397 for (n = 0; n < numPoints; ++n) { 2398 const PetscInt p = points[n]; 2399 switch (refiner) { 2400 case 1: 2401 /* Simplicial 2D */ 2402 if ((p >= vStart) && (p < vEnd)) { 2403 /* Old vertices stay the same */ 2404 newp = vStartNew + (p - vStart); 2405 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2406 } else if ((p >= fStart) && (p < fEnd)) { 2407 /* Old faces add new faces and vertex */ 2408 newp = vStartNew + (vEnd - vStart) + (p - fStart); 2409 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2410 for (r = 0; r < 2; ++r) { 2411 newp = fStartNew + (p - fStart)*2 + r; 2412 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2413 } 2414 } else if ((p >= cStart) && (p < cEnd)) { 2415 /* Old cells add new cells and interior faces */ 2416 for (r = 0; r < 4; ++r) { 2417 newp = cStartNew + (p - cStart)*4 + r; 2418 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2419 } 2420 for (r = 0; r < 3; ++r) { 2421 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 2422 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2423 } 2424 } 2425 break; 2426 case 2: 2427 /* Hex 2D */ 2428 if ((p >= vStart) && (p < vEnd)) { 2429 /* Old vertices stay the same */ 2430 newp = vStartNew + (p - vStart); 2431 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2432 } else if ((p >= fStart) && (p < fEnd)) { 2433 /* Old faces add new faces and vertex */ 2434 newp = vStartNew + (vEnd - vStart) + (p - fStart); 2435 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2436 for (r = 0; r < 2; ++r) { 2437 newp = fStartNew + (p - fStart)*2 + r; 2438 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2439 } 2440 } else if ((p >= cStart) && (p < cEnd)) { 2441 /* Old cells add new cells and interior faces and vertex */ 2442 for (r = 0; r < 4; ++r) { 2443 newp = cStartNew + (p - cStart)*4 + r; 2444 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2445 } 2446 for (r = 0; r < 4; ++r) { 2447 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 2448 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2449 } 2450 newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (p - cStart); 2451 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2452 } 2453 break; 2454 case 3: 2455 /* Hybrid simplicial 2D */ 2456 if ((p >= vStart) && (p < vEnd)) { 2457 /* Old vertices stay the same */ 2458 newp = vStartNew + (p - vStart); 2459 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2460 } else if ((p >= fStart) && (p < fMax)) { 2461 /* Old interior faces add new faces and vertex */ 2462 newp = vStartNew + (vEnd - vStart) + (p - fStart); 2463 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2464 for (r = 0; r < 2; ++r) { 2465 newp = fStartNew + (p - fStart)*2 + r; 2466 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2467 } 2468 } else if ((p >= fMax) && (p < fEnd)) { 2469 /* Old hybrid faces stay the same */ 2470 newp = fStartNew + (fMax - fStart)*2 + (p - fMax); 2471 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2472 } else if ((p >= cStart) && (p < cMax)) { 2473 /* Old interior cells add new cells and interior faces */ 2474 for (r = 0; r < 4; ++r) { 2475 newp = cStartNew + (p - cStart)*4 + r; 2476 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2477 } 2478 for (r = 0; r < 3; ++r) { 2479 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 2480 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2481 } 2482 } else if ((p >= cMax) && (p < cEnd)) { 2483 /* Old hybrid cells add new cells and hybrid face */ 2484 for (r = 0; r < 2; ++r) { 2485 newp = cStartNew + (cMax - cStart)*4 + (p - cMax)*2 + r; 2486 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2487 } 2488 newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 2489 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2490 } 2491 break; 2492 case 5: 2493 /* Simplicial 3D */ 2494 if ((p >= vStart) && (p < vEnd)) { 2495 /* Old vertices stay the same */ 2496 newp = vStartNew + (p - vStart); 2497 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2498 } else if ((p >= eStart) && (p < eEnd)) { 2499 /* Old edges add new edges and vertex */ 2500 for (r = 0; r < 2; ++r) { 2501 newp = eStartNew + (p - eStart)*2 + r; 2502 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2503 } 2504 newp = vStartNew + (vEnd - vStart) + (p - eStart); 2505 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2506 } else if ((p >= fStart) && (p < fEnd)) { 2507 /* Old faces add new faces and edges */ 2508 for (r = 0; r < 4; ++r) { 2509 newp = fStartNew + (p - fStart)*4 + r; 2510 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2511 } 2512 for (r = 0; r < 3; ++r) { 2513 newp = eStartNew + (eEnd - eStart)*2 + (p - fStart)*3 + r; 2514 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2515 } 2516 } else if ((p >= cStart) && (p < cEnd)) { 2517 /* Old cells add new cells and interior faces and edges */ 2518 for (r = 0; r < 8; ++r) { 2519 newp = cStartNew + (p - cStart)*8 + r; 2520 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2521 } 2522 for (r = 0; r < 8; ++r) { 2523 newp = fStartNew + (fEnd - fStart)*4 + (p - cStart)*8 + r; 2524 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2525 } 2526 for (r = 0; r < 1; ++r) { 2527 newp = eStartNew + (eEnd - eStart)*2 + (fEnd - fStart)*3 + (p - cStart)*1 + r; 2528 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 2529 } 2530 } 2531 break; 2532 default: 2533 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 2534 } 2535 } 2536 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2537 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2538 } 2539 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 2540 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 2541 if (0) { 2542 ierr = PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);CHKERRQ(ierr); 2543 ierr = DMLabelView(labelNew, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 2544 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 2545 } 2546 } 2547 PetscFunctionReturn(0); 2548 } 2549 2550 #undef __FUNCT__ 2551 #define __FUNCT__ "DMPlexRefineUniform_Internal" 2552 /* This will only work for interpolated meshes */ 2553 PetscErrorCode DMPlexRefineUniform_Internal(DM dm, CellRefiner cellRefiner, DM *dmRefined) 2554 { 2555 DM rdm; 2556 PetscInt *depthSize; 2557 PetscInt dim, depth = 0, d, pStart = 0, pEnd = 0; 2558 PetscErrorCode ierr; 2559 2560 PetscFunctionBegin; 2561 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &rdm);CHKERRQ(ierr); 2562 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2563 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2564 ierr = DMPlexSetDimension(rdm, dim);CHKERRQ(ierr); 2565 /* Calculate number of new points of each depth */ 2566 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2567 ierr = PetscMalloc((depth+1) * sizeof(PetscInt), &depthSize);CHKERRQ(ierr); 2568 ierr = PetscMemzero(depthSize, (depth+1) * sizeof(PetscInt));CHKERRQ(ierr); 2569 ierr = CellRefinerGetSizes(cellRefiner, dm, depthSize);CHKERRQ(ierr); 2570 /* Step 1: Set chart */ 2571 for (d = 0; d <= depth; ++d) pEnd += depthSize[d]; 2572 ierr = DMPlexSetChart(rdm, pStart, pEnd);CHKERRQ(ierr); 2573 /* Step 2: Set cone/support sizes */ 2574 ierr = CellRefinerSetConeSizes(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 2575 /* Step 3: Setup refined DM */ 2576 ierr = DMSetUp(rdm);CHKERRQ(ierr); 2577 /* Step 4: Set cones and supports */ 2578 ierr = CellRefinerSetCones(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 2579 /* Step 5: Stratify */ 2580 ierr = DMPlexStratify(rdm);CHKERRQ(ierr); 2581 /* Step 6: Set coordinates for vertices */ 2582 ierr = CellRefinerSetCoordinates(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 2583 /* Step 7: Create pointSF */ 2584 ierr = CellRefinerCreateSF(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 2585 /* Step 8: Create labels */ 2586 ierr = CellRefinerCreateLabels(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 2587 ierr = PetscFree(depthSize);CHKERRQ(ierr); 2588 2589 *dmRefined = rdm; 2590 PetscFunctionReturn(0); 2591 } 2592 2593 #undef __FUNCT__ 2594 #define __FUNCT__ "DMPlexSetRefinementUniform" 2595 PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform) 2596 { 2597 DM_Plex *mesh = (DM_Plex*) dm->data; 2598 2599 PetscFunctionBegin; 2600 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2601 mesh->refinementUniform = refinementUniform; 2602 PetscFunctionReturn(0); 2603 } 2604 2605 #undef __FUNCT__ 2606 #define __FUNCT__ "DMPlexGetRefinementUniform" 2607 PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform) 2608 { 2609 DM_Plex *mesh = (DM_Plex*) dm->data; 2610 2611 PetscFunctionBegin; 2612 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2613 PetscValidPointer(refinementUniform, 2); 2614 *refinementUniform = mesh->refinementUniform; 2615 PetscFunctionReturn(0); 2616 } 2617 2618 #undef __FUNCT__ 2619 #define __FUNCT__ "DMPlexSetRefinementLimit" 2620 PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit) 2621 { 2622 DM_Plex *mesh = (DM_Plex*) dm->data; 2623 2624 PetscFunctionBegin; 2625 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2626 mesh->refinementLimit = refinementLimit; 2627 PetscFunctionReturn(0); 2628 } 2629 2630 #undef __FUNCT__ 2631 #define __FUNCT__ "DMPlexGetRefinementLimit" 2632 PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit) 2633 { 2634 DM_Plex *mesh = (DM_Plex*) dm->data; 2635 2636 PetscFunctionBegin; 2637 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2638 PetscValidPointer(refinementLimit, 2); 2639 /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */ 2640 *refinementLimit = mesh->refinementLimit; 2641 PetscFunctionReturn(0); 2642 } 2643 2644 #undef __FUNCT__ 2645 #define __FUNCT__ "DMPlexGetCellRefiner_Internal" 2646 PetscErrorCode DMPlexGetCellRefiner_Internal(DM dm, CellRefiner *cellRefiner) 2647 { 2648 PetscInt dim, cStart, coneSize, cMax; 2649 PetscErrorCode ierr; 2650 2651 PetscFunctionBegin; 2652 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 2653 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 2654 ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr); 2655 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2656 switch (dim) { 2657 case 2: 2658 switch (coneSize) { 2659 case 3: 2660 if (cMax >= 0) *cellRefiner = 3; /* Hybrid */ 2661 else *cellRefiner = 1; /* Triangular */ 2662 break; 2663 case 4: 2664 if (cMax >= 0) *cellRefiner = 4; /* Hybrid */ 2665 else *cellRefiner = 2; /* Quadrilateral */ 2666 break; 2667 default: 2668 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown coneSize %d in dimension %d for cell refiner", coneSize, dim); 2669 } 2670 break; 2671 case 3: 2672 switch (coneSize) { 2673 case 4: 2674 if (cMax >= 0) *cellRefiner = 7; /* Hybrid */ 2675 else *cellRefiner = 5; /* Tetrahedral */ 2676 break; 2677 default: 2678 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown coneSize %d in dimension %d for cell refiner", coneSize, dim); 2679 } 2680 break; 2681 default: 2682 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown dimension %d for cell refiner", dim); 2683 } 2684 PetscFunctionReturn(0); 2685 } 2686