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