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