1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/ 2 #include <petscsf.h> 3 4 PetscBool SBRcite = PETSC_FALSE; 5 const char SBRCitation[] = "@article{PlazaCarey2000,\n" 6 " title = {Local refinement of simplicial grids based on the skeleton},\n" 7 " journal = {Applied Numerical Mathematics},\n" 8 " author = {A. Plaza and Graham F. Carey},\n" 9 " volume = {32},\n" 10 " number = {3},\n" 11 " pages = {195--218},\n" 12 " doi = {10.1016/S0168-9274(99)00022-7},\n" 13 " year = {2000}\n}\n"; 14 15 static PetscErrorCode SBRGetEdgeLen_Private(DMPlexTransform tr, PetscInt edge, PetscReal *len) { 16 DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data; 17 DM dm; 18 PetscInt off; 19 20 PetscFunctionBeginHot; 21 PetscCall(DMPlexTransformGetDM(tr, &dm)); 22 PetscCall(PetscSectionGetOffset(sbr->secEdgeLen, edge, &off)); 23 if (sbr->edgeLen[off] <= 0.0) { 24 DM cdm; 25 Vec coordsLocal; 26 const PetscScalar *coords; 27 const PetscInt *cone; 28 PetscScalar *cA, *cB; 29 PetscInt coneSize, cdim; 30 31 PetscCall(DMGetCoordinateDM(dm, &cdm)); 32 PetscCall(DMPlexGetCone(dm, edge, &cone)); 33 PetscCall(DMPlexGetConeSize(dm, edge, &coneSize)); 34 PetscCheck(coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Edge %" PetscInt_FMT " cone size must be 2, not %" PetscInt_FMT, edge, coneSize); 35 PetscCall(DMGetCoordinateDim(dm, &cdim)); 36 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 37 PetscCall(VecGetArrayRead(coordsLocal, &coords)); 38 PetscCall(DMPlexPointLocalRead(cdm, cone[0], coords, &cA)); 39 PetscCall(DMPlexPointLocalRead(cdm, cone[1], coords, &cB)); 40 sbr->edgeLen[off] = DMPlex_DistD_Internal(cdim, cA, cB); 41 PetscCall(VecRestoreArrayRead(coordsLocal, &coords)); 42 } 43 *len = sbr->edgeLen[off]; 44 PetscFunctionReturn(0); 45 } 46 47 /* Mark local edges that should be split */ 48 /* TODO This will not work in 3D */ 49 static PetscErrorCode SBRSplitLocalEdges_Private(DMPlexTransform tr, DMPlexPointQueue queue) { 50 DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data; 51 DM dm; 52 53 PetscFunctionBegin; 54 PetscCall(DMPlexTransformGetDM(tr, &dm)); 55 while (!DMPlexPointQueueEmpty(queue)) { 56 PetscInt p = -1; 57 const PetscInt *support; 58 PetscInt supportSize, s; 59 60 PetscCall(DMPlexPointQueueDequeue(queue, &p)); 61 PetscCall(DMPlexGetSupport(dm, p, &support)); 62 PetscCall(DMPlexGetSupportSize(dm, p, &supportSize)); 63 for (s = 0; s < supportSize; ++s) { 64 const PetscInt cell = support[s]; 65 const PetscInt *cone; 66 PetscInt coneSize, c; 67 PetscInt cval, eval, maxedge; 68 PetscReal len, maxlen; 69 70 PetscCall(DMLabelGetValue(sbr->splitPoints, cell, &cval)); 71 if (cval == 2) continue; 72 PetscCall(DMPlexGetCone(dm, cell, &cone)); 73 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 74 PetscCall(SBRGetEdgeLen_Private(tr, cone[0], &maxlen)); 75 maxedge = cone[0]; 76 for (c = 1; c < coneSize; ++c) { 77 PetscCall(SBRGetEdgeLen_Private(tr, cone[c], &len)); 78 if (len > maxlen) { 79 maxlen = len; 80 maxedge = cone[c]; 81 } 82 } 83 PetscCall(DMLabelGetValue(sbr->splitPoints, maxedge, &eval)); 84 if (eval != 1) { 85 PetscCall(DMLabelSetValue(sbr->splitPoints, maxedge, 1)); 86 PetscCall(DMPlexPointQueueEnqueue(queue, maxedge)); 87 } 88 PetscCall(DMLabelSetValue(sbr->splitPoints, cell, 2)); 89 } 90 } 91 PetscFunctionReturn(0); 92 } 93 94 static PetscErrorCode splitPoint(PETSC_UNUSED DMLabel label, PetscInt p, PETSC_UNUSED PetscInt val, void *ctx) { 95 DMPlexPointQueue queue = (DMPlexPointQueue)ctx; 96 97 PetscFunctionBegin; 98 PetscCall(DMPlexPointQueueEnqueue(queue, p)); 99 PetscFunctionReturn(0); 100 } 101 102 /* 103 The 'splitPoints' label marks mesh points to be divided. It marks edges with 1, triangles with 2, and tetrahedra with 3. 104 Then the refinement type is calculated as follows: 105 106 vertex: 0 107 edge unsplit: 1 108 edge split: 2 109 triangle unsplit: 3 110 triangle split all edges: 4 111 triangle split edges 0 1: 5 112 triangle split edges 1 0: 6 113 triangle split edges 1 2: 7 114 triangle split edges 2 1: 8 115 triangle split edges 2 0: 9 116 triangle split edges 0 2: 10 117 triangle split edge 0: 11 118 triangle split edge 1: 12 119 triangle split edge 2: 13 120 */ 121 typedef enum { 122 RT_VERTEX, 123 RT_EDGE, 124 RT_EDGE_SPLIT, 125 RT_TRIANGLE, 126 RT_TRIANGLE_SPLIT, 127 RT_TRIANGLE_SPLIT_01, 128 RT_TRIANGLE_SPLIT_10, 129 RT_TRIANGLE_SPLIT_12, 130 RT_TRIANGLE_SPLIT_21, 131 RT_TRIANGLE_SPLIT_20, 132 RT_TRIANGLE_SPLIT_02, 133 RT_TRIANGLE_SPLIT_0, 134 RT_TRIANGLE_SPLIT_1, 135 RT_TRIANGLE_SPLIT_2 136 } RefinementType; 137 138 static PetscErrorCode DMPlexTransformSetUp_SBR(DMPlexTransform tr) { 139 DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data; 140 DM dm; 141 DMLabel active; 142 PetscSF pointSF; 143 DMPlexPointQueue queue = NULL; 144 IS refineIS; 145 const PetscInt *refineCells; 146 PetscInt pStart, pEnd, p, eStart, eEnd, e, edgeLenSize, Nc, c; 147 PetscBool empty; 148 149 PetscFunctionBegin; 150 PetscCall(DMPlexTransformGetDM(tr, &dm)); 151 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Split Points", &sbr->splitPoints)); 152 /* Create edge lengths */ 153 PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd)); 154 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen)); 155 PetscCall(PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd)); 156 for (e = eStart; e < eEnd; ++e) PetscCall(PetscSectionSetDof(sbr->secEdgeLen, e, 1)); 157 PetscCall(PetscSectionSetUp(sbr->secEdgeLen)); 158 PetscCall(PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize)); 159 PetscCall(PetscCalloc1(edgeLenSize, &sbr->edgeLen)); 160 /* Add edges of cells that are marked for refinement to edge queue */ 161 PetscCall(DMPlexTransformGetActive(tr, &active)); 162 PetscCheck(active, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONGSTATE, "DMPlexTransform must have an adaptation label in order to use SBR algorithm"); 163 PetscCall(DMPlexPointQueueCreate(1024, &queue)); 164 PetscCall(DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS)); 165 PetscCall(DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc)); 166 if (refineIS) PetscCall(ISGetIndices(refineIS, &refineCells)); 167 for (c = 0; c < Nc; ++c) { 168 const PetscInt cell = refineCells[c]; 169 PetscInt depth; 170 171 PetscCall(DMPlexGetPointDepth(dm, cell, &depth)); 172 if (depth == 1) { 173 PetscCall(DMLabelSetValue(sbr->splitPoints, cell, 1)); 174 PetscCall(DMPlexPointQueueEnqueue(queue, cell)); 175 } else { 176 PetscInt *closure = NULL; 177 PetscInt Ncl, cl; 178 179 PetscCall(DMLabelSetValue(sbr->splitPoints, cell, depth)); 180 PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure)); 181 for (cl = 0; cl < Ncl; cl += 2) { 182 const PetscInt edge = closure[cl]; 183 184 if (edge >= eStart && edge < eEnd) { 185 PetscCall(DMLabelSetValue(sbr->splitPoints, edge, 1)); 186 PetscCall(DMPlexPointQueueEnqueue(queue, edge)); 187 } 188 } 189 PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure)); 190 } 191 } 192 if (refineIS) PetscCall(ISRestoreIndices(refineIS, &refineCells)); 193 PetscCall(ISDestroy(&refineIS)); 194 /* Setup communication */ 195 PetscCall(DMGetPointSF(dm, &pointSF)); 196 PetscCall(DMLabelPropagateBegin(sbr->splitPoints, pointSF)); 197 /* While edge queue is not empty: */ 198 PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty)); 199 while (!empty) { 200 PetscCall(SBRSplitLocalEdges_Private(tr, queue)); 201 /* Communicate marked edges 202 An easy implementation is to allocate an array the size of the number of points. We put the splitPoints marks into the 203 array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent. 204 205 TODO: We could use in-place communication with a different SF 206 We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has 207 already been marked. If not, it might have been handled on the process in this round, but we add it anyway. 208 209 In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new 210 values will have 1+0=1 and old values will have 1+1=2. Loop over these, resetting the values to 1, and adding any new 211 edge to the queue. 212 */ 213 PetscCall(DMLabelPropagatePush(sbr->splitPoints, pointSF, splitPoint, queue)); 214 PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty)); 215 } 216 PetscCall(DMLabelPropagateEnd(sbr->splitPoints, pointSF)); 217 /* Calculate refineType for each cell */ 218 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType)); 219 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 220 for (p = pStart; p < pEnd; ++p) { 221 DMLabel trType = tr->trType; 222 DMPolytopeType ct; 223 PetscInt val; 224 225 PetscCall(DMPlexGetCellType(dm, p, &ct)); 226 switch (ct) { 227 case DM_POLYTOPE_POINT: PetscCall(DMLabelSetValue(trType, p, RT_VERTEX)); break; 228 case DM_POLYTOPE_SEGMENT: 229 PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val)); 230 if (val == 1) PetscCall(DMLabelSetValue(trType, p, RT_EDGE_SPLIT)); 231 else PetscCall(DMLabelSetValue(trType, p, RT_EDGE)); 232 break; 233 case DM_POLYTOPE_TRIANGLE: 234 PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val)); 235 if (val == 2) { 236 const PetscInt *cone; 237 PetscReal lens[3]; 238 PetscInt vals[3], i; 239 240 PetscCall(DMPlexGetCone(dm, p, &cone)); 241 for (i = 0; i < 3; ++i) { 242 PetscCall(DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i])); 243 vals[i] = vals[i] < 0 ? 0 : vals[i]; 244 PetscCall(SBRGetEdgeLen_Private(tr, cone[i], &lens[i])); 245 } 246 if (vals[0] && vals[1] && vals[2]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT)); 247 else if (vals[0] && vals[1]) PetscCall(DMLabelSetValue(trType, p, lens[0] > lens[1] ? RT_TRIANGLE_SPLIT_01 : RT_TRIANGLE_SPLIT_10)); 248 else if (vals[1] && vals[2]) PetscCall(DMLabelSetValue(trType, p, lens[1] > lens[2] ? RT_TRIANGLE_SPLIT_12 : RT_TRIANGLE_SPLIT_21)); 249 else if (vals[2] && vals[0]) PetscCall(DMLabelSetValue(trType, p, lens[2] > lens[0] ? RT_TRIANGLE_SPLIT_20 : RT_TRIANGLE_SPLIT_02)); 250 else if (vals[0]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_0)); 251 else if (vals[1]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_1)); 252 else if (vals[2]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_2)); 253 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT " does not fit any refinement type (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")", p, vals[0], vals[1], vals[2]); 254 } else PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE)); 255 break; 256 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]); 257 } 258 PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val)); 259 } 260 /* Cleanup */ 261 PetscCall(DMPlexPointQueueDestroy(&queue)); 262 PetscFunctionReturn(0); 263 } 264 265 static PetscErrorCode DMPlexTransformGetSubcellOrientation_SBR(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) { 266 PetscInt rt; 267 268 PetscFunctionBeginHot; 269 PetscCall(DMLabelGetValue(tr->trType, sp, &rt)); 270 *rnew = r; 271 *onew = o; 272 switch (rt) { 273 case RT_TRIANGLE_SPLIT_01: 274 case RT_TRIANGLE_SPLIT_10: 275 case RT_TRIANGLE_SPLIT_12: 276 case RT_TRIANGLE_SPLIT_21: 277 case RT_TRIANGLE_SPLIT_20: 278 case RT_TRIANGLE_SPLIT_02: 279 switch (tct) { 280 case DM_POLYTOPE_SEGMENT: break; 281 case DM_POLYTOPE_TRIANGLE: break; 282 default: break; 283 } 284 break; 285 case RT_TRIANGLE_SPLIT_0: 286 case RT_TRIANGLE_SPLIT_1: 287 case RT_TRIANGLE_SPLIT_2: 288 switch (tct) { 289 case DM_POLYTOPE_SEGMENT: break; 290 case DM_POLYTOPE_TRIANGLE: 291 *onew = so < 0 ? -(o + 1) : o; 292 *rnew = so < 0 ? (r + 1) % 2 : r; 293 break; 294 default: break; 295 } 296 break; 297 case RT_EDGE_SPLIT: 298 case RT_TRIANGLE_SPLIT: PetscCall(DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew)); break; 299 default: PetscCall(DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew)); 300 } 301 PetscFunctionReturn(0); 302 } 303 304 /* Add 1 edge inside this triangle, making 2 new triangles. 305 2 306 |\ 307 | \ 308 | \ 309 | \ 310 | 1 311 | \ 312 | B \ 313 2 1 314 | / \ 315 | ____/ 0 316 |/ A \ 317 0-----0-----1 318 */ 319 static PetscErrorCode SBRGetTriangleSplitSingle(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) { 320 const PetscInt *arr = DMPolytopeTypeGetArrangment(DM_POLYTOPE_TRIANGLE, o); 321 static DMPolytopeType triT1[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE}; 322 static PetscInt triS1[] = {1, 2}; 323 static PetscInt triC1[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, 324 DM_POLYTOPE_SEGMENT, 0, 0}; 325 static PetscInt triO1[] = {0, 0, 0, 0, -1, 0, 0, 0}; 326 327 PetscFunctionBeginHot; 328 /* To get the other divisions, we reorient the triangle */ 329 triC1[2] = arr[0 * 2]; 330 triC1[7] = arr[1 * 2]; 331 triC1[11] = arr[0 * 2]; 332 triC1[15] = arr[1 * 2]; 333 triC1[22] = arr[1 * 2]; 334 triC1[26] = arr[2 * 2]; 335 *Nt = 2; 336 *target = triT1; 337 *size = triS1; 338 *cone = triC1; 339 *ornt = triO1; 340 PetscFunctionReturn(0); 341 } 342 343 /* Add 2 edges inside this triangle, making 3 new triangles. 344 RT_TRIANGLE_SPLIT_12 345 2 346 |\ 347 | \ 348 | \ 349 0 \ 350 | 1 351 | \ 352 | B \ 353 2-------1 354 | C / \ 355 1 ____/ 0 356 |/ A \ 357 0-----0-----1 358 RT_TRIANGLE_SPLIT_10 359 2 360 |\ 361 | \ 362 | \ 363 0 \ 364 | 1 365 | \ 366 | A \ 367 2 1 368 | /|\ 369 1 ____/ / 0 370 |/ C / B \ 371 0-----0-----1 372 RT_TRIANGLE_SPLIT_20 373 2 374 |\ 375 | \ 376 | \ 377 0 \ 378 | \ 379 | \ 380 | \ 381 2 A 1 382 |\ \ 383 1 ---\ \ 384 |B \_C----\\ 385 0-----0-----1 386 RT_TRIANGLE_SPLIT_21 387 2 388 |\ 389 | \ 390 | \ 391 0 \ 392 | \ 393 | B \ 394 | \ 395 2-------1 396 |\ C \ 397 1 ---\ \ 398 | A ----\\ 399 0-----0-----1 400 RT_TRIANGLE_SPLIT_01 401 2 402 |\ 403 |\\ 404 || \ 405 | \ \ 406 | | \ 407 | | \ 408 | | \ 409 2 \ C 1 410 | A | / \ 411 | | |B \ 412 | \/ \ 413 0-----0-----1 414 RT_TRIANGLE_SPLIT_02 415 2 416 |\ 417 |\\ 418 || \ 419 | \ \ 420 | | \ 421 | | \ 422 | | \ 423 2 C \ 1 424 |\ | \ 425 | \__| A \ 426 | B \\ \ 427 0-----0-----1 428 */ 429 static PetscErrorCode SBRGetTriangleSplitDouble(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) { 430 PetscInt e0, e1; 431 const PetscInt *arr = DMPolytopeTypeGetArrangment(DM_POLYTOPE_TRIANGLE, o); 432 static DMPolytopeType triT2[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE}; 433 static PetscInt triS2[] = {2, 3}; 434 static PetscInt triC2[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1}; 435 static PetscInt triO2[] = {0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0}; 436 437 PetscFunctionBeginHot; 438 /* To get the other divisions, we reorient the triangle */ 439 triC2[2] = arr[0 * 2]; 440 triC2[3] = arr[0 * 2 + 1] ? 1 : 0; 441 triC2[7] = arr[1 * 2]; 442 triC2[11] = arr[1 * 2]; 443 triC2[15] = arr[2 * 2]; 444 /* Swap the first two edges if the triangle is reversed */ 445 e0 = o < 0 ? 23 : 19; 446 e1 = o < 0 ? 19 : 23; 447 triC2[e0] = arr[0 * 2]; 448 triC2[e0 + 1] = 0; 449 triC2[e1] = arr[1 * 2]; 450 triC2[e1 + 1] = o < 0 ? 1 : 0; 451 triO2[6] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]); 452 /* Swap the first two edges if the triangle is reversed */ 453 e0 = o < 0 ? 34 : 30; 454 e1 = o < 0 ? 30 : 34; 455 triC2[e0] = arr[1 * 2]; 456 triC2[e0 + 1] = o < 0 ? 0 : 1; 457 triC2[e1] = arr[2 * 2]; 458 triC2[e1 + 1] = o < 0 ? 1 : 0; 459 triO2[9] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]); 460 /* Swap the last two edges if the triangle is reversed */ 461 triC2[41] = arr[2 * 2]; 462 triC2[42] = o < 0 ? 0 : 1; 463 triC2[45] = o < 0 ? 1 : 0; 464 triC2[48] = o < 0 ? 0 : 1; 465 triO2[11] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[1 * 2 + 1]); 466 triO2[12] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[2 * 2 + 1]); 467 *Nt = 2; 468 *target = triT2; 469 *size = triS2; 470 *cone = triC2; 471 *ornt = triO2; 472 PetscFunctionReturn(0); 473 } 474 475 static PetscErrorCode DMPlexTransformCellTransform_SBR(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) { 476 DMLabel trType = tr->trType; 477 PetscInt val; 478 479 PetscFunctionBeginHot; 480 PetscCheck(p >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point argument is invalid"); 481 PetscCall(DMLabelGetValue(trType, p, &val)); 482 if (rt) *rt = val; 483 switch (source) { 484 case DM_POLYTOPE_POINT: 485 case DM_POLYTOPE_POINT_PRISM_TENSOR: 486 case DM_POLYTOPE_QUADRILATERAL: 487 case DM_POLYTOPE_SEG_PRISM_TENSOR: 488 case DM_POLYTOPE_TETRAHEDRON: 489 case DM_POLYTOPE_HEXAHEDRON: 490 case DM_POLYTOPE_TRI_PRISM: 491 case DM_POLYTOPE_TRI_PRISM_TENSOR: 492 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 493 case DM_POLYTOPE_PYRAMID: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt)); break; 494 case DM_POLYTOPE_SEGMENT: 495 if (val == RT_EDGE) PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt)); 496 else PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt)); 497 break; 498 case DM_POLYTOPE_TRIANGLE: 499 switch (val) { 500 case RT_TRIANGLE_SPLIT_0: PetscCall(SBRGetTriangleSplitSingle(2, Nt, target, size, cone, ornt)); break; 501 case RT_TRIANGLE_SPLIT_1: PetscCall(SBRGetTriangleSplitSingle(0, Nt, target, size, cone, ornt)); break; 502 case RT_TRIANGLE_SPLIT_2: PetscCall(SBRGetTriangleSplitSingle(1, Nt, target, size, cone, ornt)); break; 503 case RT_TRIANGLE_SPLIT_21: PetscCall(SBRGetTriangleSplitDouble(-3, Nt, target, size, cone, ornt)); break; 504 case RT_TRIANGLE_SPLIT_10: PetscCall(SBRGetTriangleSplitDouble(-2, Nt, target, size, cone, ornt)); break; 505 case RT_TRIANGLE_SPLIT_02: PetscCall(SBRGetTriangleSplitDouble(-1, Nt, target, size, cone, ornt)); break; 506 case RT_TRIANGLE_SPLIT_12: PetscCall(SBRGetTriangleSplitDouble(0, Nt, target, size, cone, ornt)); break; 507 case RT_TRIANGLE_SPLIT_20: PetscCall(SBRGetTriangleSplitDouble(1, Nt, target, size, cone, ornt)); break; 508 case RT_TRIANGLE_SPLIT_01: PetscCall(SBRGetTriangleSplitDouble(2, Nt, target, size, cone, ornt)); break; 509 case RT_TRIANGLE_SPLIT: PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt)); break; 510 default: PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt)); 511 } 512 break; 513 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 514 } 515 PetscFunctionReturn(0); 516 } 517 518 static PetscErrorCode DMPlexTransformSetFromOptions_SBR(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject) { 519 PetscInt cells[256], n = 256, i; 520 PetscBool flg; 521 522 PetscFunctionBegin; 523 PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Options"); 524 PetscCall(PetscOptionsIntArray("-dm_plex_transform_sbr_ref_cell", "Mark cells for refinement", "", cells, &n, &flg)); 525 if (flg) { 526 DMLabel active; 527 528 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active)); 529 for (i = 0; i < n; ++i) PetscCall(DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE)); 530 PetscCall(DMPlexTransformSetActive(tr, active)); 531 PetscCall(DMLabelDestroy(&active)); 532 } 533 PetscOptionsHeadEnd(); 534 PetscFunctionReturn(0); 535 } 536 537 static PetscErrorCode DMPlexTransformView_SBR(DMPlexTransform tr, PetscViewer viewer) { 538 PetscBool isascii; 539 540 PetscFunctionBegin; 541 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 542 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 543 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 544 if (isascii) { 545 PetscViewerFormat format; 546 const char *name; 547 548 PetscCall(PetscObjectGetName((PetscObject)tr, &name)); 549 PetscCall(PetscViewerASCIIPrintf(viewer, "SBR refinement %s\n", name ? name : "")); 550 PetscCall(PetscViewerGetFormat(viewer, &format)); 551 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(DMLabelView(tr->trType, viewer)); 552 } else { 553 SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name); 554 } 555 PetscFunctionReturn(0); 556 } 557 558 static PetscErrorCode DMPlexTransformDestroy_SBR(DMPlexTransform tr) { 559 DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data; 560 561 PetscFunctionBegin; 562 PetscCall(PetscFree(sbr->edgeLen)); 563 PetscCall(PetscSectionDestroy(&sbr->secEdgeLen)); 564 PetscCall(DMLabelDestroy(&sbr->splitPoints)); 565 PetscCall(PetscFree(tr->data)); 566 PetscFunctionReturn(0); 567 } 568 569 static PetscErrorCode DMPlexTransformInitialize_SBR(DMPlexTransform tr) { 570 PetscFunctionBegin; 571 tr->ops->view = DMPlexTransformView_SBR; 572 tr->ops->setfromoptions = DMPlexTransformSetFromOptions_SBR; 573 tr->ops->setup = DMPlexTransformSetUp_SBR; 574 tr->ops->destroy = DMPlexTransformDestroy_SBR; 575 tr->ops->celltransform = DMPlexTransformCellTransform_SBR; 576 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_SBR; 577 tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal; 578 PetscFunctionReturn(0); 579 } 580 581 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform tr) { 582 DMPlexRefine_SBR *f; 583 584 PetscFunctionBegin; 585 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 586 PetscCall(PetscNewLog(tr, &f)); 587 tr->data = f; 588 589 PetscCall(DMPlexTransformInitialize_SBR(tr)); 590 PetscCall(PetscCitationsRegister(SBRCitation, &SBRcite)); 591 PetscFunctionReturn(0); 592 } 593