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