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__ "CellRefinerGetSizes" 30 PetscErrorCode CellRefinerGetSizes(CellRefiner refiner, DM dm, PetscInt depthSize[]) 31 { 32 PetscInt cStart, cEnd, cMax, vStart, vEnd, vMax, fStart, fEnd, fMax, eStart, eEnd, eMax; 33 PetscErrorCode ierr; 34 35 PetscFunctionBegin; 36 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 37 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 38 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 39 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 40 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 41 switch (refiner) { 42 case 1: 43 /* Simplicial 2D */ 44 depthSize[0] = vEnd - vStart + fEnd - fStart; /* Add a vertex on every face */ 45 depthSize[1] = 2*(fEnd - fStart) + 3*(cEnd - cStart); /* Every face is split into 2 faces and 3 faces are added for each cell */ 46 depthSize[2] = 4*(cEnd - cStart); /* Every cell split into 4 cells */ 47 break; 48 case 3: 49 /* Hybrid 2D */ 50 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 51 cMax = PetscMin(cEnd, cMax); 52 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 53 fMax = PetscMin(fEnd, fMax); 54 depthSize[0] = vEnd - vStart + fMax - fStart; /* Add a vertex on every face, but not hybrid faces */ 55 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 */ 56 depthSize[2] = 4*(cMax - cStart) + 2*(cEnd - cMax); /* Interior cells split into 4 cells, Hybrid cells split into 2 cells */ 57 break; 58 case 2: 59 /* Hex 2D */ 60 depthSize[0] = vEnd - vStart + cEnd - cStart + fEnd - fStart; /* Add a vertex on every face and cell */ 61 depthSize[1] = 2*(fEnd - fStart) + 4*(cEnd - cStart); /* Every face is split into 2 faces and 4 faces are added for each cell */ 62 depthSize[2] = 4*(cEnd - cStart); /* Every cell split into 4 cells */ 63 break; 64 default: 65 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 66 } 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "CellRefinerSetConeSizes" 72 PetscErrorCode CellRefinerSetConeSizes(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 73 { 74 PetscInt depth, cStart, cStartNew, cEnd, cMax, c, vStart, vStartNew, vEnd, vMax, v, fStart, fStartNew, fEnd, fMax, f, eStart, eStartNew, eEnd, eMax, r; 75 PetscErrorCode ierr; 76 77 PetscFunctionBegin; 78 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 79 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 80 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 81 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 82 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 83 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 84 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 85 switch (refiner) { 86 case 1: 87 /* Simplicial 2D */ 88 /* All cells have 3 faces */ 89 for (c = cStart; c < cEnd; ++c) { 90 for (r = 0; r < 4; ++r) { 91 const PetscInt newp = (c - cStart)*4 + r; 92 93 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 94 } 95 } 96 /* Split faces have 2 vertices and the same cells as the parent */ 97 for (f = fStart; f < fEnd; ++f) { 98 for (r = 0; r < 2; ++r) { 99 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 100 PetscInt size; 101 102 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 103 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 104 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 105 } 106 } 107 /* Interior faces have 2 vertices and 2 cells */ 108 for (c = cStart; c < cEnd; ++c) { 109 for (r = 0; r < 3; ++r) { 110 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + r; 111 112 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 113 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 114 } 115 } 116 /* Old vertices have identical supports */ 117 for (v = vStart; v < vEnd; ++v) { 118 const PetscInt newp = vStartNew + (v - vStart); 119 PetscInt size; 120 121 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 122 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 123 } 124 /* Face vertices have 2 + cells*2 supports */ 125 for (f = fStart; f < fEnd; ++f) { 126 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 127 PetscInt size; 128 129 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 130 ierr = DMPlexSetSupportSize(rdm, newp, 2 + size*2);CHKERRQ(ierr); 131 } 132 break; 133 case 2: 134 /* Hex 2D */ 135 /* All cells have 4 faces */ 136 for (c = cStart; c < cEnd; ++c) { 137 for (r = 0; r < 4; ++r) { 138 const PetscInt newp = (c - cStart)*4 + r; 139 140 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 141 } 142 } 143 /* Split faces have 2 vertices and the same cells as the parent */ 144 for (f = fStart; f < fEnd; ++f) { 145 for (r = 0; r < 2; ++r) { 146 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 147 PetscInt size; 148 149 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 150 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 151 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 152 } 153 } 154 /* Interior faces have 2 vertices and 2 cells */ 155 for (c = cStart; c < cEnd; ++c) { 156 for (r = 0; r < 4; ++r) { 157 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + r; 158 159 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 160 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 161 } 162 } 163 /* Old vertices have identical supports */ 164 for (v = vStart; v < vEnd; ++v) { 165 const PetscInt newp = vStartNew + (v - vStart); 166 PetscInt size; 167 168 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 169 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 170 } 171 /* Face vertices have 2 + cells supports */ 172 for (f = fStart; f < fEnd; ++f) { 173 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 174 PetscInt size; 175 176 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 177 ierr = DMPlexSetSupportSize(rdm, newp, 2 + size);CHKERRQ(ierr); 178 } 179 /* Cell vertices have 4 supports */ 180 for (c = cStart; c < cEnd; ++c) { 181 const PetscInt newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 182 183 ierr = DMPlexSetSupportSize(rdm, newp, 4);CHKERRQ(ierr); 184 } 185 break; 186 case 3: 187 /* Hybrid 2D */ 188 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 189 cMax = PetscMin(cEnd, cMax); 190 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 191 fMax = PetscMin(fEnd, fMax); 192 ierr = DMPlexSetHybridBounds(rdm, cStartNew + (cMax - cStart)*4, fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 193 /* Interior cells have 3 faces */ 194 for (c = cStart; c < cMax; ++c) { 195 for (r = 0; r < 4; ++r) { 196 const PetscInt newp = cStartNew + (c - cStart)*4 + r; 197 198 ierr = DMPlexSetConeSize(rdm, newp, 3);CHKERRQ(ierr); 199 } 200 } 201 /* Hybrid cells have 4 faces */ 202 for (c = cMax; c < cEnd; ++c) { 203 for (r = 0; r < 2; ++r) { 204 const PetscInt newp = cStartNew + (cMax - cStart)*4 + (c - cMax)*2 + r; 205 206 ierr = DMPlexSetConeSize(rdm, newp, 4);CHKERRQ(ierr); 207 } 208 } 209 /* Interior split faces have 2 vertices and the same cells as the parent */ 210 for (f = fStart; f < fMax; ++f) { 211 for (r = 0; r < 2; ++r) { 212 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 213 PetscInt size; 214 215 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 216 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 217 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 218 } 219 } 220 /* Interior cell faces have 2 vertices and 2 cells */ 221 for (c = cStart; c < cMax; ++c) { 222 for (r = 0; r < 3; ++r) { 223 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + r; 224 225 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 226 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 227 } 228 } 229 /* Hybrid faces have 2 vertices and the same cells */ 230 for (f = fMax; f < fEnd; ++f) { 231 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (f - fMax); 232 PetscInt size; 233 234 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 235 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 236 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 237 } 238 /* Hybrid cell faces have 2 vertices and 2 cells */ 239 for (c = cMax; c < cEnd; ++c) { 240 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 241 242 ierr = DMPlexSetConeSize(rdm, newp, 2);CHKERRQ(ierr); 243 ierr = DMPlexSetSupportSize(rdm, newp, 2);CHKERRQ(ierr); 244 } 245 /* Old vertices have identical supports */ 246 for (v = vStart; v < vEnd; ++v) { 247 const PetscInt newp = vStartNew + (v - vStart); 248 PetscInt size; 249 250 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 251 ierr = DMPlexSetSupportSize(rdm, newp, size);CHKERRQ(ierr); 252 } 253 /* Face vertices have 2 + (2 interior, 1 hybrid) supports */ 254 for (f = fStart; f < fMax; ++f) { 255 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 256 const PetscInt *support; 257 PetscInt size, newSize = 2, s; 258 259 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 260 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 261 for (s = 0; s < size; ++s) { 262 if (support[s] >= cMax) newSize += 1; 263 else newSize += 2; 264 } 265 ierr = DMPlexSetSupportSize(rdm, newp, newSize);CHKERRQ(ierr); 266 } 267 break; 268 default: 269 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 270 } 271 PetscFunctionReturn(0); 272 } 273 274 #undef __FUNCT__ 275 #define __FUNCT__ "CellRefinerSetCones" 276 PetscErrorCode CellRefinerSetCones(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 277 { 278 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, r, p; 279 PetscInt maxSupportSize, *supportRef; 280 PetscErrorCode ierr; 281 282 PetscFunctionBegin; 283 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 284 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 285 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 286 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 287 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 288 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 289 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 290 ierr = GetDepthEnd_Private(depth, depthSize, &cEndNew, &fEndNew, &eEndNew, &vEndNew);CHKERRQ(ierr); 291 switch (refiner) { 292 case 1: 293 /* Simplicial 2D */ 294 /* 295 2 296 |\ 297 | \ 298 | \ 299 | \ 300 | C \ 301 | \ 302 | \ 303 2---1---1 304 |\ D / \ 305 | 2 0 \ 306 |A \ / B \ 307 0---0-------1 308 */ 309 /* All cells have 3 faces */ 310 for (c = cStart; c < cEnd; ++c) { 311 const PetscInt newp = cStartNew + (c - cStart)*4; 312 const PetscInt *cone, *ornt; 313 PetscInt coneNew[3], orntNew[3]; 314 315 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 316 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 317 /* A triangle */ 318 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 319 orntNew[0] = ornt[0]; 320 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 2; 321 orntNew[1] = -2; 322 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 323 orntNew[2] = ornt[2]; 324 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 325 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 326 #if 1 327 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); 328 for (p = 0; p < 3; ++p) { 329 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); 330 } 331 #endif 332 /* B triangle */ 333 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 334 orntNew[0] = ornt[0]; 335 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 336 orntNew[1] = ornt[1]; 337 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 0; 338 orntNew[2] = -2; 339 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 340 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 341 #if 1 342 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); 343 for (p = 0; p < 3; ++p) { 344 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); 345 } 346 #endif 347 /* C triangle */ 348 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 1; 349 orntNew[0] = -2; 350 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 351 orntNew[1] = ornt[1]; 352 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 353 orntNew[2] = ornt[2]; 354 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 355 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 356 #if 1 357 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); 358 for (p = 0; p < 3; ++p) { 359 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); 360 } 361 #endif 362 /* D triangle */ 363 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 0; 364 orntNew[0] = 0; 365 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 1; 366 orntNew[1] = 0; 367 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + 2; 368 orntNew[2] = 0; 369 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 370 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 371 #if 1 372 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); 373 for (p = 0; p < 3; ++p) { 374 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); 375 } 376 #endif 377 } 378 /* Split faces have 2 vertices and the same cells as the parent */ 379 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 380 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 381 for (f = fStart; f < fEnd; ++f) { 382 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 383 384 for (r = 0; r < 2; ++r) { 385 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 386 const PetscInt *cone, *support; 387 PetscInt coneNew[2], coneSize, c, supportSize, s; 388 389 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 390 coneNew[0] = vStartNew + (cone[0] - vStart); 391 coneNew[1] = vStartNew + (cone[1] - vStart); 392 coneNew[(r+1)%2] = newv; 393 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 394 #if 1 395 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 396 for (p = 0; p < 2; ++p) { 397 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); 398 } 399 #endif 400 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 401 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 402 for (s = 0; s < supportSize; ++s) { 403 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 404 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 405 for (c = 0; c < coneSize; ++c) { 406 if (cone[c] == f) break; 407 } 408 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (c+r)%3; 409 } 410 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 411 #if 1 412 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 413 for (p = 0; p < supportSize; ++p) { 414 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); 415 } 416 #endif 417 } 418 } 419 /* Interior faces have 2 vertices and 2 cells */ 420 for (c = cStart; c < cEnd; ++c) { 421 const PetscInt *cone; 422 423 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 424 for (r = 0; r < 3; ++r) { 425 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*3 + r; 426 PetscInt coneNew[2]; 427 PetscInt supportNew[2]; 428 429 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 430 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - fStart); 431 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 432 #if 1 433 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 434 for (p = 0; p < 2; ++p) { 435 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); 436 } 437 #endif 438 supportNew[0] = (c - cStart)*4 + (r+1)%3; 439 supportNew[1] = (c - cStart)*4 + 3; 440 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 441 #if 1 442 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 443 for (p = 0; p < 2; ++p) { 444 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); 445 } 446 #endif 447 } 448 } 449 /* Old vertices have identical supports */ 450 for (v = vStart; v < vEnd; ++v) { 451 const PetscInt newp = vStartNew + (v - vStart); 452 const PetscInt *support, *cone; 453 PetscInt size, s; 454 455 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 456 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 457 for (s = 0; s < size; ++s) { 458 PetscInt r = 0; 459 460 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 461 if (cone[1] == v) r = 1; 462 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 463 } 464 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 465 #if 1 466 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 467 for (p = 0; p < size; ++p) { 468 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); 469 } 470 #endif 471 } 472 /* Face vertices have 2 + cells*2 supports */ 473 for (f = fStart; f < fEnd; ++f) { 474 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 475 const PetscInt *cone, *support; 476 PetscInt size, s; 477 478 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 479 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 480 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 481 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 482 for (s = 0; s < size; ++s) { 483 PetscInt r = 0; 484 485 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 486 if (cone[1] == f) r = 1; 487 else if (cone[2] == f) r = 2; 488 supportRef[2+s*2+0] = fStartNew + (fEnd - fStart)*2 + (support[s] - cStart)*3 + (r+2)%3; 489 supportRef[2+s*2+1] = fStartNew + (fEnd - fStart)*2 + (support[s] - cStart)*3 + r; 490 } 491 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 492 #if 1 493 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 494 for (p = 0; p < 2+size*2; ++p) { 495 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); 496 } 497 #endif 498 } 499 ierr = PetscFree(supportRef);CHKERRQ(ierr); 500 break; 501 case 2: 502 /* Hex 2D */ 503 /* 504 3---------2---------2 505 | | | 506 | D 2 C | 507 | | | 508 3----3----0----1----1 509 | | | 510 | A 0 B | 511 | | | 512 0---------0---------1 513 */ 514 /* All cells have 4 faces */ 515 for (c = cStart; c < cEnd; ++c) { 516 const PetscInt newp = (c - cStart)*4; 517 const PetscInt *cone, *ornt; 518 PetscInt coneNew[4], orntNew[4]; 519 520 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 521 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 522 /* A quad */ 523 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 524 orntNew[0] = ornt[0]; 525 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 0; 526 orntNew[1] = 0; 527 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 3; 528 orntNew[2] = -2; 529 coneNew[3] = fStartNew + (cone[3] - fStart)*2 + (ornt[3] < 0 ? 0 : 1); 530 orntNew[3] = ornt[3]; 531 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 532 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 533 #if 1 534 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); 535 for (p = 0; p < 4; ++p) { 536 if ((coneNew[p] < fStartNew) || (coneNew[p] >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", coneNew[p], fStartNew, fEndNew); 537 } 538 #endif 539 /* B quad */ 540 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 541 orntNew[0] = ornt[0]; 542 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 543 orntNew[1] = ornt[1]; 544 coneNew[2] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 1; 545 orntNew[2] = 0; 546 coneNew[3] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 0; 547 orntNew[3] = -2; 548 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 549 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 550 #if 1 551 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); 552 for (p = 0; p < 4; ++p) { 553 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); 554 } 555 #endif 556 /* C quad */ 557 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 1; 558 orntNew[0] = -2; 559 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 560 orntNew[1] = ornt[1]; 561 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 562 orntNew[2] = ornt[2]; 563 coneNew[3] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 2; 564 orntNew[3] = 0; 565 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 566 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 567 #if 1 568 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); 569 for (p = 0; p < 4; ++p) { 570 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); 571 } 572 #endif 573 /* D quad */ 574 coneNew[0] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 3; 575 orntNew[0] = 0; 576 coneNew[1] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + 2; 577 orntNew[1] = -2; 578 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 579 orntNew[2] = ornt[2]; 580 coneNew[3] = fStartNew + (cone[3] - fStart)*2 + (ornt[3] < 0 ? 1 : 0); 581 orntNew[3] = ornt[3]; 582 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 583 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 584 #if 1 585 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); 586 for (p = 0; p < 4; ++p) { 587 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); 588 } 589 #endif 590 } 591 /* Split faces have 2 vertices and the same cells as the parent */ 592 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 593 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 594 for (f = fStart; f < fEnd; ++f) { 595 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 596 597 for (r = 0; r < 2; ++r) { 598 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 599 const PetscInt *cone, *support; 600 PetscInt coneNew[2], coneSize, c, supportSize, s; 601 602 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 603 coneNew[0] = vStartNew + (cone[0] - vStart); 604 coneNew[1] = vStartNew + (cone[1] - vStart); 605 coneNew[(r+1)%2] = newv; 606 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 607 #if 1 608 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 609 for (p = 0; p < 2; ++p) { 610 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); 611 } 612 #endif 613 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 614 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 615 for (s = 0; s < supportSize; ++s) { 616 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 617 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 618 for (c = 0; c < coneSize; ++c) { 619 if (cone[c] == f) break; 620 } 621 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (c+r)%4; 622 } 623 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 624 #if 1 625 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 626 for (p = 0; p < supportSize; ++p) { 627 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); 628 } 629 #endif 630 } 631 } 632 /* Interior faces have 2 vertices and 2 cells */ 633 for (c = cStart; c < cEnd; ++c) { 634 const PetscInt *cone; 635 PetscInt coneNew[2], supportNew[2]; 636 637 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 638 for (r = 0; r < 4; ++r) { 639 const PetscInt newp = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + r; 640 641 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 642 coneNew[1] = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 643 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 644 #if 1 645 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 646 for (p = 0; p < 2; ++p) { 647 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); 648 } 649 #endif 650 supportNew[0] = (c - cStart)*4 + r; 651 supportNew[1] = (c - cStart)*4 + (r+1)%4; 652 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 653 #if 1 654 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 655 for (p = 0; p < 2; ++p) { 656 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); 657 } 658 #endif 659 } 660 } 661 /* Old vertices have identical supports */ 662 for (v = vStart; v < vEnd; ++v) { 663 const PetscInt newp = vStartNew + (v - vStart); 664 const PetscInt *support, *cone; 665 PetscInt size, s; 666 667 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 668 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 669 for (s = 0; s < size; ++s) { 670 PetscInt r = 0; 671 672 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 673 if (cone[1] == v) r = 1; 674 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 675 } 676 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 677 #if 1 678 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 679 for (p = 0; p < size; ++p) { 680 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); 681 } 682 #endif 683 } 684 /* Face vertices have 2 + cells supports */ 685 for (f = fStart; f < fEnd; ++f) { 686 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 687 const PetscInt *cone, *support; 688 PetscInt size, s; 689 690 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 691 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 692 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 693 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 694 for (s = 0; s < size; ++s) { 695 PetscInt r = 0; 696 697 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 698 if (cone[1] == f) r = 1; 699 else if (cone[2] == f) r = 2; 700 else if (cone[3] == f) r = 3; 701 supportRef[2+s] = fStartNew + (fEnd - fStart)*2 + (support[s] - cStart)*4 + r; 702 } 703 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 704 #if 1 705 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 706 for (p = 0; p < 2+size; ++p) { 707 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); 708 } 709 #endif 710 } 711 /* Cell vertices have 4 supports */ 712 for (c = cStart; c < cEnd; ++c) { 713 const PetscInt newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 714 PetscInt supportNew[4]; 715 716 for (r = 0; r < 4; ++r) { 717 supportNew[r] = fStartNew + (fEnd - fStart)*2 + (c - cStart)*4 + r; 718 } 719 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 720 } 721 break; 722 case 3: 723 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 724 cMax = PetscMin(cEnd, cMax); 725 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 726 fMax = PetscMin(fEnd, fMax); 727 /* Interior cells have 3 faces */ 728 for (c = cStart; c < cMax; ++c) { 729 const PetscInt newp = cStartNew + (c - cStart)*4; 730 const PetscInt *cone, *ornt; 731 PetscInt coneNew[3], orntNew[3]; 732 733 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 734 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 735 /* A triangle */ 736 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 737 orntNew[0] = ornt[0]; 738 coneNew[1] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 2; 739 orntNew[1] = -2; 740 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 0 : 1); 741 orntNew[2] = ornt[2]; 742 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 743 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 744 #if 1 745 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); 746 for (p = 0; p < 3; ++p) { 747 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); 748 } 749 #endif 750 /* B triangle */ 751 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 752 orntNew[0] = ornt[0]; 753 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 754 orntNew[1] = ornt[1]; 755 coneNew[2] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 0; 756 orntNew[2] = -2; 757 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 758 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 759 #if 1 760 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); 761 for (p = 0; p < 3; ++p) { 762 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); 763 } 764 #endif 765 /* C triangle */ 766 coneNew[0] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 1; 767 orntNew[0] = -2; 768 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 769 orntNew[1] = ornt[1]; 770 coneNew[2] = fStartNew + (cone[2] - fStart)*2 + (ornt[2] < 0 ? 1 : 0); 771 orntNew[2] = ornt[2]; 772 ierr = DMPlexSetCone(rdm, newp+2, coneNew);CHKERRQ(ierr); 773 ierr = DMPlexSetConeOrientation(rdm, newp+2, orntNew);CHKERRQ(ierr); 774 #if 1 775 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); 776 for (p = 0; p < 3; ++p) { 777 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); 778 } 779 #endif 780 /* D triangle */ 781 coneNew[0] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 0; 782 orntNew[0] = 0; 783 coneNew[1] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 1; 784 orntNew[1] = 0; 785 coneNew[2] = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + 2; 786 orntNew[2] = 0; 787 ierr = DMPlexSetCone(rdm, newp+3, coneNew);CHKERRQ(ierr); 788 ierr = DMPlexSetConeOrientation(rdm, newp+3, orntNew);CHKERRQ(ierr); 789 #if 1 790 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); 791 for (p = 0; p < 3; ++p) { 792 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); 793 } 794 #endif 795 } 796 /* 797 2----3----3 798 | | 799 | B | 800 | | 801 0----4--- 1 802 | | 803 | A | 804 | | 805 0----2----1 806 */ 807 /* Hybrid cells have 4 faces */ 808 for (c = cMax; c < cEnd; ++c) { 809 const PetscInt newp = cStartNew + (cMax - cStart)*4 + (c - cMax)*2; 810 const PetscInt *cone, *ornt; 811 PetscInt coneNew[4], orntNew[4]; 812 813 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 814 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 815 /* A quad */ 816 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 1 : 0); 817 orntNew[0] = ornt[0]; 818 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 1 : 0); 819 orntNew[1] = ornt[1]; 820 coneNew[2] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (cone[2] - fMax); 821 orntNew[2] = 0; 822 coneNew[3] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 823 orntNew[3] = 0; 824 ierr = DMPlexSetCone(rdm, newp+0, coneNew);CHKERRQ(ierr); 825 ierr = DMPlexSetConeOrientation(rdm, newp+0, orntNew);CHKERRQ(ierr); 826 #if 1 827 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); 828 for (p = 0; p < 4; ++p) { 829 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); 830 } 831 #endif 832 /* B quad */ 833 coneNew[0] = fStartNew + (cone[0] - fStart)*2 + (ornt[0] < 0 ? 0 : 1); 834 orntNew[0] = ornt[0]; 835 coneNew[1] = fStartNew + (cone[1] - fStart)*2 + (ornt[1] < 0 ? 0 : 1); 836 orntNew[1] = ornt[1]; 837 coneNew[2] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 838 orntNew[2] = 0; 839 coneNew[3] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (cone[3] - fMax); 840 orntNew[3] = 0; 841 ierr = DMPlexSetCone(rdm, newp+1, coneNew);CHKERRQ(ierr); 842 ierr = DMPlexSetConeOrientation(rdm, newp+1, orntNew);CHKERRQ(ierr); 843 #if 1 844 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); 845 for (p = 0; p < 4; ++p) { 846 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); 847 } 848 #endif 849 } 850 /* Interior split faces have 2 vertices and the same cells as the parent */ 851 ierr = DMPlexGetMaxSizes(dm, NULL, &maxSupportSize);CHKERRQ(ierr); 852 ierr = PetscMalloc((2 + maxSupportSize*2) * sizeof(PetscInt), &supportRef);CHKERRQ(ierr); 853 for (f = fStart; f < fMax; ++f) { 854 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 855 856 for (r = 0; r < 2; ++r) { 857 const PetscInt newp = fStartNew + (f - fStart)*2 + r; 858 const PetscInt *cone, *support; 859 PetscInt coneNew[2], coneSize, c, supportSize, s; 860 861 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 862 coneNew[0] = vStartNew + (cone[0] - vStart); 863 coneNew[1] = vStartNew + (cone[1] - vStart); 864 coneNew[(r+1)%2] = newv; 865 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 866 #if 1 867 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 868 for (p = 0; p < 2; ++p) { 869 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); 870 } 871 #endif 872 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 873 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 874 for (s = 0; s < supportSize; ++s) { 875 if (support[s] >= cMax) { 876 supportRef[s] = cStartNew + (cMax - cStart)*4 + (support[s] - cMax)*2 + r; 877 } else { 878 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 879 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 880 for (c = 0; c < coneSize; ++c) { 881 if (cone[c] == f) break; 882 } 883 supportRef[s] = cStartNew + (support[s] - cStart)*4 + (c+r)%3; 884 } 885 } 886 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 887 #if 1 888 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 889 for (p = 0; p < supportSize; ++p) { 890 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); 891 } 892 #endif 893 } 894 } 895 /* Interior cell faces have 2 vertices and 2 cells */ 896 for (c = cStart; c < cMax; ++c) { 897 const PetscInt *cone; 898 899 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 900 for (r = 0; r < 3; ++r) { 901 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (c - cStart)*3 + r; 902 PetscInt coneNew[2]; 903 PetscInt supportNew[2]; 904 905 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[r] - fStart); 906 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[(r+1)%3] - fStart); 907 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 908 #if 1 909 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 910 for (p = 0; p < 2; ++p) { 911 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); 912 } 913 #endif 914 supportNew[0] = (c - cStart)*4 + (r+1)%3; 915 supportNew[1] = (c - cStart)*4 + 3; 916 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 917 #if 1 918 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 919 for (p = 0; p < 2; ++p) { 920 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); 921 } 922 #endif 923 } 924 } 925 /* Interior hybrid faces have 2 vertices and the same cells */ 926 for (f = fMax; f < fEnd; ++f) { 927 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (f - fMax); 928 const PetscInt *cone; 929 const PetscInt *support; 930 PetscInt coneNew[2]; 931 PetscInt supportNew[2]; 932 PetscInt size, s, r; 933 934 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 935 coneNew[0] = vStartNew + (cone[0] - vStart); 936 coneNew[1] = vStartNew + (cone[1] - vStart); 937 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 938 #if 1 939 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 940 for (p = 0; p < 2; ++p) { 941 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); 942 } 943 #endif 944 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 945 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 946 for (s = 0; s < size; ++s) { 947 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 948 for (r = 0; r < 2; ++r) { 949 if (cone[r+2] == f) break; 950 } 951 supportNew[s] = (cMax - cStart)*4 + (support[s] - cMax)*2 + r; 952 } 953 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 954 #if 1 955 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 956 for (p = 0; p < size; ++p) { 957 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); 958 } 959 #endif 960 } 961 /* Cell hybrid faces have 2 vertices and 2 cells */ 962 for (c = cMax; c < cEnd; ++c) { 963 const PetscInt newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (c - cMax); 964 const PetscInt *cone; 965 PetscInt coneNew[2]; 966 PetscInt supportNew[2]; 967 968 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 969 coneNew[0] = vStartNew + (vEnd - vStart) + (cone[0] - fStart); 970 coneNew[1] = vStartNew + (vEnd - vStart) + (cone[1] - fStart); 971 ierr = DMPlexSetCone(rdm, newp, coneNew);CHKERRQ(ierr); 972 #if 1 973 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 974 for (p = 0; p < 2; ++p) { 975 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); 976 } 977 #endif 978 supportNew[0] = (cMax - cStart)*4 + (c - cMax)*2 + 0; 979 supportNew[1] = (cMax - cStart)*4 + (c - cMax)*2 + 1; 980 ierr = DMPlexSetSupport(rdm, newp, supportNew);CHKERRQ(ierr); 981 #if 1 982 if ((newp < fStartNew) || (newp >= fEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a face [%d, %d)", newp, fStartNew, fEndNew); 983 for (p = 0; p < 2; ++p) { 984 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); 985 } 986 #endif 987 } 988 /* Old vertices have identical supports */ 989 for (v = vStart; v < vEnd; ++v) { 990 const PetscInt newp = vStartNew + (v - vStart); 991 const PetscInt *support, *cone; 992 PetscInt size, s; 993 994 ierr = DMPlexGetSupportSize(dm, v, &size);CHKERRQ(ierr); 995 ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr); 996 for (s = 0; s < size; ++s) { 997 if (support[s] >= fMax) { 998 supportRef[s] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (support[s] - fMax); 999 } else { 1000 PetscInt r = 0; 1001 1002 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1003 if (cone[1] == v) r = 1; 1004 supportRef[s] = fStartNew + (support[s] - fStart)*2 + r; 1005 } 1006 } 1007 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1008 #if 1 1009 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1010 for (p = 0; p < size; ++p) { 1011 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); 1012 } 1013 #endif 1014 } 1015 /* Face vertices have 2 + (2 interior, 1 hybrid) supports */ 1016 for (f = fStart; f < fMax; ++f) { 1017 const PetscInt newp = vStartNew + (vEnd - vStart) + (f - fStart); 1018 const PetscInt *cone, *support; 1019 PetscInt size, newSize = 2, s; 1020 1021 ierr = DMPlexGetSupportSize(dm, f, &size);CHKERRQ(ierr); 1022 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 1023 supportRef[0] = fStartNew + (f - fStart)*2 + 0; 1024 supportRef[1] = fStartNew + (f - fStart)*2 + 1; 1025 for (s = 0; s < size; ++s) { 1026 PetscInt r = 0; 1027 1028 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1029 if (support[s] >= cMax) { 1030 supportRef[newSize+0] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (fEnd - fMax) + (support[s] - cMax); 1031 1032 newSize += 1; 1033 } else { 1034 if (cone[1] == f) r = 1; 1035 else if (cone[2] == f) r = 2; 1036 supportRef[newSize+0] = fStartNew + (fMax - fStart)*2 + (support[s] - cStart)*3 + (r+2)%3; 1037 supportRef[newSize+1] = fStartNew + (fMax - fStart)*2 + (support[s] - cStart)*3 + r; 1038 1039 newSize += 2; 1040 } 1041 } 1042 ierr = DMPlexSetSupport(rdm, newp, supportRef);CHKERRQ(ierr); 1043 #if 1 1044 if ((newp < vStartNew) || (newp >= vEndNew)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d is not a vertex [%d, %d)", newp, vStartNew, vEndNew); 1045 for (p = 0; p < newSize; ++p) { 1046 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); 1047 } 1048 #endif 1049 } 1050 ierr = PetscFree(supportRef);CHKERRQ(ierr); 1051 break; 1052 default: 1053 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 1054 } 1055 PetscFunctionReturn(0); 1056 } 1057 1058 #undef __FUNCT__ 1059 #define __FUNCT__ "CellRefinerSetCoordinates" 1060 PetscErrorCode CellRefinerSetCoordinates(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 1061 { 1062 PetscSection coordSection, coordSectionNew; 1063 Vec coordinates, coordinatesNew; 1064 PetscScalar *coords, *coordsNew; 1065 PetscInt dim, depth, coordSizeNew, cStart, cEnd, c, vStart, vStartNew, vEnd, v, fStart, fEnd, fMax, f; 1066 PetscErrorCode ierr; 1067 1068 PetscFunctionBegin; 1069 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1070 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1071 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1072 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1073 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1074 ierr = DMPlexGetHybridBounds(dm, NULL, &fMax, NULL, NULL);CHKERRQ(ierr); 1075 ierr = GetDepthStart_Private(depth, depthSize, NULL, NULL, NULL, &vStartNew);CHKERRQ(ierr); 1076 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1077 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &coordSectionNew);CHKERRQ(ierr); 1078 ierr = PetscSectionSetNumFields(coordSectionNew, 1);CHKERRQ(ierr); 1079 ierr = PetscSectionSetFieldComponents(coordSectionNew, 0, dim);CHKERRQ(ierr); 1080 ierr = PetscSectionSetChart(coordSectionNew, vStartNew, vStartNew+depthSize[0]);CHKERRQ(ierr); 1081 if (fMax < 0) fMax = fEnd; 1082 switch (refiner) { 1083 case 1: 1084 case 2: 1085 case 3: 1086 /* Simplicial and Hex 2D */ 1087 /* All vertices have the dim coordinates */ 1088 for (v = vStartNew; v < vStartNew+depthSize[0]; ++v) { 1089 ierr = PetscSectionSetDof(coordSectionNew, v, dim);CHKERRQ(ierr); 1090 ierr = PetscSectionSetFieldDof(coordSectionNew, v, 0, dim);CHKERRQ(ierr); 1091 } 1092 ierr = PetscSectionSetUp(coordSectionNew);CHKERRQ(ierr); 1093 ierr = DMPlexSetCoordinateSection(rdm, coordSectionNew);CHKERRQ(ierr); 1094 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1095 ierr = PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);CHKERRQ(ierr); 1096 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinatesNew);CHKERRQ(ierr); 1097 ierr = PetscObjectSetName((PetscObject) coordinatesNew, "coordinates");CHKERRQ(ierr); 1098 ierr = VecSetSizes(coordinatesNew, coordSizeNew, PETSC_DETERMINE);CHKERRQ(ierr); 1099 ierr = VecSetFromOptions(coordinatesNew);CHKERRQ(ierr); 1100 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1101 ierr = VecGetArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 1102 /* Old vertices have the same coordinates */ 1103 for (v = vStart; v < vEnd; ++v) { 1104 const PetscInt newv = vStartNew + (v - vStart); 1105 PetscInt off, offnew, d; 1106 1107 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1108 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 1109 for (d = 0; d < dim; ++d) { 1110 coordsNew[offnew+d] = coords[off+d]; 1111 } 1112 } 1113 /* Face vertices have the average of endpoint coordinates */ 1114 for (f = fStart; f < fMax; ++f) { 1115 const PetscInt newv = vStartNew + (vEnd - vStart) + (f - fStart); 1116 const PetscInt *cone; 1117 PetscInt coneSize, offA, offB, offnew, d; 1118 1119 ierr = DMPlexGetConeSize(dm, f, &coneSize);CHKERRQ(ierr); 1120 if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Face %d cone should have two vertices, not %d", f, coneSize); 1121 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1122 ierr = PetscSectionGetOffset(coordSection, cone[0], &offA);CHKERRQ(ierr); 1123 ierr = PetscSectionGetOffset(coordSection, cone[1], &offB);CHKERRQ(ierr); 1124 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 1125 for (d = 0; d < dim; ++d) { 1126 coordsNew[offnew+d] = 0.5*(coords[offA+d] + coords[offB+d]); 1127 } 1128 } 1129 /* Just Hex 2D */ 1130 if (refiner == 2) { 1131 /* Cell vertices have the average of corner coordinates */ 1132 for (c = cStart; c < cEnd; ++c) { 1133 const PetscInt newv = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (c - cStart); 1134 PetscInt *cone = NULL; 1135 PetscInt closureSize, coneSize = 0, offA, offB, offC, offD, offnew, p, d; 1136 1137 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 1138 for (p = 0; p < closureSize*2; p += 2) { 1139 const PetscInt point = cone[p]; 1140 if ((point >= vStart) && (point < vEnd)) cone[coneSize++] = point; 1141 } 1142 if (coneSize != 4) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Quad %d cone should have four vertices, not %d", c, coneSize); 1143 ierr = PetscSectionGetOffset(coordSection, cone[0], &offA);CHKERRQ(ierr); 1144 ierr = PetscSectionGetOffset(coordSection, cone[1], &offB);CHKERRQ(ierr); 1145 ierr = PetscSectionGetOffset(coordSection, cone[2], &offC);CHKERRQ(ierr); 1146 ierr = PetscSectionGetOffset(coordSection, cone[3], &offD);CHKERRQ(ierr); 1147 ierr = PetscSectionGetOffset(coordSectionNew, newv, &offnew);CHKERRQ(ierr); 1148 for (d = 0; d < dim; ++d) { 1149 coordsNew[offnew+d] = 0.25*(coords[offA+d] + coords[offB+d] + coords[offC+d] + coords[offD+d]); 1150 } 1151 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &cone);CHKERRQ(ierr); 1152 } 1153 } 1154 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1155 ierr = VecRestoreArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 1156 ierr = DMSetCoordinatesLocal(rdm, coordinatesNew);CHKERRQ(ierr); 1157 ierr = VecDestroy(&coordinatesNew);CHKERRQ(ierr); 1158 ierr = PetscSectionDestroy(&coordSectionNew);CHKERRQ(ierr); 1159 break; 1160 default: 1161 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 1162 } 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "DMPlexCreateProcessSF" 1168 PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess) 1169 { 1170 PetscInt numRoots, numLeaves, l; 1171 const PetscInt *localPoints; 1172 const PetscSFNode *remotePoints; 1173 PetscInt *localPointsNew; 1174 PetscSFNode *remotePointsNew; 1175 PetscInt *ranks, *ranksNew; 1176 PetscErrorCode ierr; 1177 1178 PetscFunctionBegin; 1179 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1180 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranks);CHKERRQ(ierr); 1181 for (l = 0; l < numLeaves; ++l) { 1182 ranks[l] = remotePoints[l].rank; 1183 } 1184 ierr = PetscSortRemoveDupsInt(&numLeaves, ranks);CHKERRQ(ierr); 1185 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &ranksNew);CHKERRQ(ierr); 1186 ierr = PetscMalloc(numLeaves * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 1187 ierr = PetscMalloc(numLeaves * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 1188 for (l = 0; l < numLeaves; ++l) { 1189 ranksNew[l] = ranks[l]; 1190 localPointsNew[l] = l; 1191 remotePointsNew[l].index = 0; 1192 remotePointsNew[l].rank = ranksNew[l]; 1193 } 1194 ierr = PetscFree(ranks);CHKERRQ(ierr); 1195 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);CHKERRQ(ierr); 1196 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);CHKERRQ(ierr); 1197 ierr = PetscSFSetFromOptions(*sfProcess);CHKERRQ(ierr); 1198 ierr = PetscSFSetGraph(*sfProcess, 1, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1199 PetscFunctionReturn(0); 1200 } 1201 1202 #undef __FUNCT__ 1203 #define __FUNCT__ "CellRefinerCreateSF" 1204 PetscErrorCode CellRefinerCreateSF(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 1205 { 1206 PetscSF sf, sfNew, sfProcess; 1207 IS processRanks; 1208 MPI_Datatype depthType; 1209 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 1210 const PetscInt *localPoints, *neighbors; 1211 const PetscSFNode *remotePoints; 1212 PetscInt *localPointsNew; 1213 PetscSFNode *remotePointsNew; 1214 PetscInt *depthSizeOld, *rdepthSize, *rdepthSizeOld, *rdepthMaxOld, *rvStart, *rvStartNew, *reStart, *reStartNew, *rfStart, *rfStartNew, *rcStart, *rcStartNew; 1215 PetscInt depth, numNeighbors, pStartNew, pEndNew, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eStartNew, eEnd, eMax, r, n; 1216 PetscErrorCode ierr; 1217 1218 PetscFunctionBegin; 1219 ierr = DMPlexGetChart(rdm, &pStartNew, &pEndNew);CHKERRQ(ierr); 1220 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1221 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1222 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 1223 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1224 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1225 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 1226 ierr = GetDepthStart_Private(depth, depthSize, &cStartNew, &fStartNew, &eStartNew, &vStartNew);CHKERRQ(ierr); 1227 switch (refiner) { 1228 case 3: 1229 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 1230 cMax = PetscMin(cEnd, cMax); 1231 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 1232 fMax = PetscMin(fEnd, fMax); 1233 } 1234 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 1235 ierr = DMGetPointSF(rdm, &sfNew);CHKERRQ(ierr); 1236 /* Caculate size of new SF */ 1237 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1238 if (numRoots < 0) PetscFunctionReturn(0); 1239 for (l = 0; l < numLeaves; ++l) { 1240 const PetscInt p = localPoints[l]; 1241 1242 switch (refiner) { 1243 case 1: 1244 /* Simplicial 2D */ 1245 if ((p >= vStart) && (p < vEnd)) { 1246 /* Old vertices stay the same */ 1247 ++numLeavesNew; 1248 } else if ((p >= fStart) && (p < fEnd)) { 1249 /* Old faces add new faces and vertex */ 1250 numLeavesNew += 1 + 2; 1251 } else if ((p >= cStart) && (p < cEnd)) { 1252 /* Old cells add new cells and interior faces */ 1253 numLeavesNew += 4 + 3; 1254 } 1255 break; 1256 case 2: 1257 /* Hex 2D */ 1258 if ((p >= vStart) && (p < vEnd)) { 1259 /* Old vertices stay the same */ 1260 ++numLeavesNew; 1261 } else if ((p >= fStart) && (p < fEnd)) { 1262 /* Old faces add new faces and vertex */ 1263 numLeavesNew += 1 + 2; 1264 } else if ((p >= cStart) && (p < cEnd)) { 1265 /* Old cells add new cells and interior faces */ 1266 numLeavesNew += 4 + 4; 1267 } 1268 break; 1269 default: 1270 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 1271 } 1272 } 1273 /* Communicate depthSizes for each remote rank */ 1274 ierr = DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);CHKERRQ(ierr); 1275 ierr = ISGetLocalSize(processRanks, &numNeighbors);CHKERRQ(ierr); 1276 ierr = PetscMalloc5((depth+1)*numNeighbors,PetscInt,&rdepthSize,numNeighbors,PetscInt,&rvStartNew,numNeighbors,PetscInt,&reStartNew,numNeighbors,PetscInt,&rfStartNew,numNeighbors,PetscInt,&rcStartNew);CHKERRQ(ierr); 1277 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); 1278 ierr = MPI_Type_contiguous(depth+1, MPIU_INT, &depthType);CHKERRQ(ierr); 1279 ierr = MPI_Type_commit(&depthType);CHKERRQ(ierr); 1280 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 1281 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSize, rdepthSize);CHKERRQ(ierr); 1282 for (n = 0; n < numNeighbors; ++n) { 1283 ierr = GetDepthStart_Private(depth, &rdepthSize[n*(depth+1)], &rcStartNew[n], &rfStartNew[n], &reStartNew[n], &rvStartNew[n]);CHKERRQ(ierr); 1284 } 1285 depthSizeOld[depth] = cMax; 1286 depthSizeOld[0] = vMax; 1287 depthSizeOld[depth-1] = fMax; 1288 depthSizeOld[1] = eMax; 1289 1290 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 1291 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthMaxOld);CHKERRQ(ierr); 1292 1293 depthSizeOld[depth] = cEnd - cStart; 1294 depthSizeOld[0] = vEnd - vStart; 1295 depthSizeOld[depth-1] = fEnd - fStart; 1296 depthSizeOld[1] = eEnd - eStart; 1297 1298 ierr = PetscSFBcastBegin(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 1299 ierr = PetscSFBcastEnd(sfProcess, depthType, depthSizeOld, rdepthSizeOld);CHKERRQ(ierr); 1300 for (n = 0; n < numNeighbors; ++n) { 1301 ierr = GetDepthStart_Private(depth, &rdepthSizeOld[n*(depth+1)], &rcStart[n], &rfStart[n], &reStart[n], &rvStart[n]);CHKERRQ(ierr); 1302 } 1303 ierr = MPI_Type_free(&depthType);CHKERRQ(ierr); 1304 ierr = PetscSFDestroy(&sfProcess);CHKERRQ(ierr); 1305 /* Calculate new point SF */ 1306 ierr = PetscMalloc(numLeavesNew * sizeof(PetscInt), &localPointsNew);CHKERRQ(ierr); 1307 ierr = PetscMalloc(numLeavesNew * sizeof(PetscSFNode), &remotePointsNew);CHKERRQ(ierr); 1308 ierr = ISGetIndices(processRanks, &neighbors);CHKERRQ(ierr); 1309 for (l = 0, m = 0; l < numLeaves; ++l) { 1310 PetscInt p = localPoints[l]; 1311 PetscInt rp = remotePoints[l].index, n; 1312 PetscMPIInt rrank = remotePoints[l].rank; 1313 1314 ierr = PetscFindInt(rrank, numNeighbors, neighbors, &n);CHKERRQ(ierr); 1315 if (n < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %d", rrank); 1316 switch (refiner) { 1317 case 1: 1318 /* Simplicial 2D */ 1319 if ((p >= vStart) && (p < vEnd)) { 1320 /* Old vertices stay the same */ 1321 localPointsNew[m] = vStartNew + (p - vStart); 1322 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 1323 remotePointsNew[m].rank = rrank; 1324 ++m; 1325 } else if ((p >= fStart) && (p < fEnd)) { 1326 /* Old faces add new faces and vertex */ 1327 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 1328 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 1329 remotePointsNew[m].rank = rrank; 1330 ++m; 1331 for (r = 0; r < 2; ++r, ++m) { 1332 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 1333 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 1334 remotePointsNew[m].rank = rrank; 1335 } 1336 } else if ((p >= cStart) && (p < cEnd)) { 1337 /* Old cells add new cells and interior faces */ 1338 for (r = 0; r < 4; ++r, ++m) { 1339 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 1340 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 1341 remotePointsNew[m].rank = rrank; 1342 } 1343 for (r = 0; r < 3; ++r, ++m) { 1344 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 1345 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*3 + r; 1346 remotePointsNew[m].rank = rrank; 1347 } 1348 } 1349 break; 1350 case 2: 1351 /* Hex 2D */ 1352 if ((p >= vStart) && (p < vEnd)) { 1353 /* Old vertices stay the same */ 1354 localPointsNew[m] = vStartNew + (p - vStart); 1355 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 1356 remotePointsNew[m].rank = rrank; 1357 ++m; 1358 } else if ((p >= fStart) && (p < fEnd)) { 1359 /* Old faces add new faces and vertex */ 1360 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 1361 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 1362 remotePointsNew[m].rank = rrank; 1363 ++m; 1364 for (r = 0; r < 2; ++r, ++m) { 1365 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 1366 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 1367 remotePointsNew[m].rank = rrank; 1368 } 1369 } else if ((p >= cStart) && (p < cEnd)) { 1370 /* Old cells add new cells and interior faces */ 1371 for (r = 0; r < 4; ++r, ++m) { 1372 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 1373 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 1374 remotePointsNew[m].rank = rrank; 1375 } 1376 for (r = 0; r < 4; ++r, ++m) { 1377 localPointsNew[m] = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 1378 remotePointsNew[m].index = rfStartNew[n] + rdepthSizeOld[n*(depth+1)+depth-1]*2 + (rp - rcStart[n])*4 + r; 1379 remotePointsNew[m].rank = rrank; 1380 } 1381 } 1382 break; 1383 case 3: 1384 /* Hybrid simplicial 2D */ 1385 if ((p >= vStart) && (p < vEnd)) { 1386 /* Old vertices stay the same */ 1387 localPointsNew[m] = vStartNew + (p - vStart); 1388 remotePointsNew[m].index = rvStartNew[n] + (rp - rvStart[n]); 1389 remotePointsNew[m].rank = rrank; 1390 ++m; 1391 } else if ((p >= fStart) && (p < fMax)) { 1392 /* Old interior faces add new faces and vertex */ 1393 localPointsNew[m] = vStartNew + (vEnd - vStart) + (p - fStart); 1394 remotePointsNew[m].index = rvStartNew[n] + rdepthSizeOld[n*(depth+1)+0] + (rp - rfStart[n]); 1395 remotePointsNew[m].rank = rrank; 1396 ++m; 1397 for (r = 0; r < 2; ++r, ++m) { 1398 localPointsNew[m] = fStartNew + (p - fStart)*2 + r; 1399 remotePointsNew[m].index = rfStartNew[n] + (rp - rfStart[n])*2 + r; 1400 remotePointsNew[m].rank = rrank; 1401 } 1402 } else if ((p >= fMax) && (p < fEnd)) { 1403 /* Old hybrid faces stay the same */ 1404 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - fMax); 1405 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rdepthMaxOld[n*(depth+1)+depth-1]); 1406 remotePointsNew[m].rank = rrank; 1407 ++m; 1408 } else if ((p >= cStart) && (p < cMax)) { 1409 /* Old interior cells add new cells and interior faces */ 1410 for (r = 0; r < 4; ++r, ++m) { 1411 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 1412 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 1413 remotePointsNew[m].rank = rrank; 1414 } 1415 for (r = 0; r < 3; ++r, ++m) { 1416 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (p - cStart)*3 + r; 1417 remotePointsNew[m].index = rfStartNew[n] + (rdepthMaxOld[n*(depth+1)+depth-1] - rfStart[n])*2 + (rp - rcStart[n])*3 + r; 1418 remotePointsNew[m].rank = rrank; 1419 } 1420 } else if ((p >= cStart) && (p < cMax)) { 1421 /* Old hybrid cells add new cells and hybrid face */ 1422 for (r = 0; r < 2; ++r, ++m) { 1423 localPointsNew[m] = cStartNew + (p - cStart)*4 + r; 1424 remotePointsNew[m].index = rcStartNew[n] + (rp - rcStart[n])*4 + r; 1425 remotePointsNew[m].rank = rrank; 1426 } 1427 localPointsNew[m] = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 1428 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]); 1429 remotePointsNew[m].rank = rrank; 1430 ++m; 1431 } 1432 break; 1433 default: 1434 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 1435 } 1436 } 1437 ierr = ISRestoreIndices(processRanks, &neighbors);CHKERRQ(ierr); 1438 ierr = ISDestroy(&processRanks);CHKERRQ(ierr); 1439 ierr = PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1440 ierr = PetscFree5(rdepthSize,rvStartNew,reStartNew,rfStartNew,rcStartNew);CHKERRQ(ierr); 1441 ierr = PetscFree6(depthSizeOld,rdepthSizeOld,rvStart,reStart,rfStart,rcStart);CHKERRQ(ierr); 1442 PetscFunctionReturn(0); 1443 } 1444 1445 #undef __FUNCT__ 1446 #define __FUNCT__ "CellRefinerCreateLabels" 1447 PetscErrorCode CellRefinerCreateLabels(CellRefiner refiner, DM dm, PetscInt depthSize[], DM rdm) 1448 { 1449 PetscInt numLabels, l; 1450 PetscInt newp, cStart, cStartNew, cEnd, cMax, vStart, vStartNew, vEnd, vMax, fStart, fStartNew, fEnd, fMax, eStart, eEnd, eMax, r; 1451 PetscErrorCode ierr; 1452 1453 PetscFunctionBegin; 1454 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1455 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 1456 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1457 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1458 1459 cStartNew = 0; 1460 vStartNew = depthSize[2]; 1461 fStartNew = depthSize[2] + depthSize[0]; 1462 1463 ierr = DMPlexGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 1464 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);CHKERRQ(ierr); 1465 switch (refiner) { 1466 case 3: 1467 if (cMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No cell maximum specified in hybrid mesh"); 1468 cMax = PetscMin(cEnd, cMax); 1469 if (fMax < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No face maximum specified in hybrid mesh"); 1470 fMax = PetscMin(fEnd, fMax); 1471 } 1472 for (l = 0; l < numLabels; ++l) { 1473 DMLabel label, labelNew; 1474 const char *lname; 1475 PetscBool isDepth; 1476 IS valueIS; 1477 const PetscInt *values; 1478 PetscInt numValues, val; 1479 1480 ierr = DMPlexGetLabelName(dm, l, &lname);CHKERRQ(ierr); 1481 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 1482 if (isDepth) continue; 1483 ierr = DMPlexCreateLabel(rdm, lname);CHKERRQ(ierr); 1484 ierr = DMPlexGetLabel(dm, lname, &label);CHKERRQ(ierr); 1485 ierr = DMPlexGetLabel(rdm, lname, &labelNew);CHKERRQ(ierr); 1486 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 1487 ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr); 1488 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 1489 for (val = 0; val < numValues; ++val) { 1490 IS pointIS; 1491 const PetscInt *points; 1492 PetscInt numPoints, n; 1493 1494 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 1495 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 1496 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1497 for (n = 0; n < numPoints; ++n) { 1498 const PetscInt p = points[n]; 1499 switch (refiner) { 1500 case 1: 1501 /* Simplicial 2D */ 1502 if ((p >= vStart) && (p < vEnd)) { 1503 /* Old vertices stay the same */ 1504 newp = vStartNew + (p - vStart); 1505 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1506 } else if ((p >= fStart) && (p < fEnd)) { 1507 /* Old faces add new faces and vertex */ 1508 newp = vStartNew + (vEnd - vStart) + (p - fStart); 1509 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1510 for (r = 0; r < 2; ++r) { 1511 newp = fStartNew + (p - fStart)*2 + r; 1512 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1513 } 1514 } else if ((p >= cStart) && (p < cEnd)) { 1515 /* Old cells add new cells and interior faces */ 1516 for (r = 0; r < 4; ++r) { 1517 newp = cStartNew + (p - cStart)*4 + r; 1518 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1519 } 1520 for (r = 0; r < 3; ++r) { 1521 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 1522 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1523 } 1524 } 1525 break; 1526 case 2: 1527 /* Hex 2D */ 1528 if ((p >= vStart) && (p < vEnd)) { 1529 /* Old vertices stay the same */ 1530 newp = vStartNew + (p - vStart); 1531 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1532 } else if ((p >= fStart) && (p < fEnd)) { 1533 /* Old faces add new faces and vertex */ 1534 newp = vStartNew + (vEnd - vStart) + (p - fStart); 1535 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1536 for (r = 0; r < 2; ++r) { 1537 newp = fStartNew + (p - fStart)*2 + r; 1538 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1539 } 1540 } else if ((p >= cStart) && (p < cEnd)) { 1541 /* Old cells add new cells and interior faces and vertex */ 1542 for (r = 0; r < 4; ++r) { 1543 newp = cStartNew + (p - cStart)*4 + r; 1544 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1545 } 1546 for (r = 0; r < 4; ++r) { 1547 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*4 + r; 1548 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1549 } 1550 newp = vStartNew + (vEnd - vStart) + (fEnd - fStart) + (p - cStart); 1551 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1552 } 1553 break; 1554 case 3: 1555 /* Hybrid simplicial 2D */ 1556 if ((p >= vStart) && (p < vEnd)) { 1557 /* Old vertices stay the same */ 1558 newp = vStartNew + (p - vStart); 1559 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1560 } else if ((p >= fStart) && (p < fMax)) { 1561 /* Old interior faces add new faces and vertex */ 1562 newp = vStartNew + (vEnd - vStart) + (p - fStart); 1563 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1564 for (r = 0; r < 2; ++r) { 1565 newp = fStartNew + (p - fStart)*2 + r; 1566 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1567 } 1568 } else if ((p >= fMax) && (p < fEnd)) { 1569 /* Old hybrid faces stay the same */ 1570 newp = fStartNew + (fMax - fStart)*2 + (p - fMax); 1571 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1572 } else if ((p >= cStart) && (p < cMax)) { 1573 /* Old interior cells add new cells and interior faces */ 1574 for (r = 0; r < 4; ++r) { 1575 newp = cStartNew + (p - cStart)*4 + r; 1576 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1577 } 1578 for (r = 0; r < 3; ++r) { 1579 newp = fStartNew + (fEnd - fStart)*2 + (p - cStart)*3 + r; 1580 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1581 } 1582 } else if ((p >= cMax) && (p < cEnd)) { 1583 /* Old hybrid cells add new cells and hybrid face */ 1584 for (r = 0; r < 2; ++r) { 1585 newp = cStartNew + (cMax - cStart)*4 + (p - cMax)*2 + r; 1586 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1587 } 1588 newp = fStartNew + (fMax - fStart)*2 + (cMax - cStart)*3 + (p - cMax); 1589 ierr = DMLabelSetValue(labelNew, newp, values[val]);CHKERRQ(ierr); 1590 } 1591 break; 1592 default: 1593 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell refiner %d", refiner); 1594 } 1595 } 1596 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1597 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1598 } 1599 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 1600 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 1601 if (0) { 1602 ierr = PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);CHKERRQ(ierr); 1603 ierr = DMLabelView(labelNew, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1604 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1605 } 1606 } 1607 PetscFunctionReturn(0); 1608 } 1609 1610 #undef __FUNCT__ 1611 #define __FUNCT__ "DMPlexRefineUniform_Internal" 1612 /* This will only work for interpolated meshes */ 1613 PetscErrorCode DMPlexRefineUniform_Internal(DM dm, CellRefiner cellRefiner, DM *dmRefined) 1614 { 1615 DM rdm; 1616 PetscInt *depthSize; 1617 PetscInt dim, depth = 0, d, pStart = 0, pEnd = 0; 1618 PetscErrorCode ierr; 1619 1620 PetscFunctionBegin; 1621 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &rdm);CHKERRQ(ierr); 1622 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 1623 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1624 ierr = DMPlexSetDimension(rdm, dim);CHKERRQ(ierr); 1625 /* Calculate number of new points of each depth */ 1626 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1627 ierr = PetscMalloc((depth+1) * sizeof(PetscInt), &depthSize);CHKERRQ(ierr); 1628 ierr = PetscMemzero(depthSize, (depth+1) * sizeof(PetscInt));CHKERRQ(ierr); 1629 ierr = CellRefinerGetSizes(cellRefiner, dm, depthSize);CHKERRQ(ierr); 1630 /* Step 1: Set chart */ 1631 for (d = 0; d <= depth; ++d) pEnd += depthSize[d]; 1632 ierr = DMPlexSetChart(rdm, pStart, pEnd);CHKERRQ(ierr); 1633 /* Step 2: Set cone/support sizes */ 1634 ierr = CellRefinerSetConeSizes(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 1635 /* Step 3: Setup refined DM */ 1636 ierr = DMSetUp(rdm);CHKERRQ(ierr); 1637 /* Step 4: Set cones and supports */ 1638 ierr = CellRefinerSetCones(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 1639 /* Step 5: Stratify */ 1640 ierr = DMPlexStratify(rdm);CHKERRQ(ierr); 1641 /* Step 6: Set coordinates for vertices */ 1642 ierr = CellRefinerSetCoordinates(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 1643 /* Step 7: Create pointSF */ 1644 ierr = CellRefinerCreateSF(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 1645 /* Step 8: Create labels */ 1646 ierr = CellRefinerCreateLabels(cellRefiner, dm, depthSize, rdm);CHKERRQ(ierr); 1647 ierr = PetscFree(depthSize);CHKERRQ(ierr); 1648 1649 *dmRefined = rdm; 1650 PetscFunctionReturn(0); 1651 } 1652 1653 #undef __FUNCT__ 1654 #define __FUNCT__ "DMPlexSetRefinementUniform" 1655 PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform) 1656 { 1657 DM_Plex *mesh = (DM_Plex*) dm->data; 1658 1659 PetscFunctionBegin; 1660 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1661 mesh->refinementUniform = refinementUniform; 1662 PetscFunctionReturn(0); 1663 } 1664 1665 #undef __FUNCT__ 1666 #define __FUNCT__ "DMPlexGetRefinementUniform" 1667 PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform) 1668 { 1669 DM_Plex *mesh = (DM_Plex*) dm->data; 1670 1671 PetscFunctionBegin; 1672 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1673 PetscValidPointer(refinementUniform, 2); 1674 *refinementUniform = mesh->refinementUniform; 1675 PetscFunctionReturn(0); 1676 } 1677 1678 #undef __FUNCT__ 1679 #define __FUNCT__ "DMPlexSetRefinementLimit" 1680 PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit) 1681 { 1682 DM_Plex *mesh = (DM_Plex*) dm->data; 1683 1684 PetscFunctionBegin; 1685 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1686 mesh->refinementLimit = refinementLimit; 1687 PetscFunctionReturn(0); 1688 } 1689 1690 #undef __FUNCT__ 1691 #define __FUNCT__ "DMPlexGetRefinementLimit" 1692 PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit) 1693 { 1694 DM_Plex *mesh = (DM_Plex*) dm->data; 1695 1696 PetscFunctionBegin; 1697 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1698 PetscValidPointer(refinementLimit, 2); 1699 /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */ 1700 *refinementLimit = mesh->refinementLimit; 1701 PetscFunctionReturn(0); 1702 } 1703 1704 #undef __FUNCT__ 1705 #define __FUNCT__ "DMPlexGetCellRefiner_Internal" 1706 PetscErrorCode DMPlexGetCellRefiner_Internal(DM dm, CellRefiner *cellRefiner) 1707 { 1708 PetscInt dim, cStart, coneSize, cMax; 1709 PetscErrorCode ierr; 1710 1711 PetscFunctionBegin; 1712 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1713 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 1714 ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr); 1715 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1716 switch (dim) { 1717 case 2: 1718 switch (coneSize) { 1719 case 3: 1720 if (cMax >= 0) *cellRefiner = 3; /* Hybrid */ 1721 else *cellRefiner = 1; /* Triangular */ 1722 break; 1723 case 4: 1724 if (cMax >= 0) *cellRefiner = 4; /* Hybrid */ 1725 else *cellRefiner = 2; /* Quadrilateral */ 1726 break; 1727 default: 1728 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown coneSize %d in dimension %d for cell refiner", coneSize, dim); 1729 } 1730 break; 1731 default: 1732 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown dimension %d for cell refiner", dim); 1733 } 1734 PetscFunctionReturn(0); 1735 } 1736