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