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