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 CHKERRQ(PetscCalloc1(1, &q)); 31 q->size = size; 32 CHKERRQ(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 CHKERRQ(PetscFree(q->points)); 46 CHKERRQ(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 CHKERRQ(PetscRealloc(queue->size * sizeof(PetscInt), &queue->points)); 57 PetscFunctionReturn(0); 58 } 59 60 static PetscErrorCode PointQueueEnqueue(PointQueue queue, PetscInt p) 61 { 62 PetscFunctionBegin; 63 CHKERRQ(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 CHKERRQ(DMPlexTransformGetDM(tr, &dm)); 112 CHKERRQ(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 CHKERRQ(DMGetCoordinateDM(dm, &cdm)); 122 CHKERRQ(DMPlexGetCone(dm, edge, &cone)); 123 CHKERRQ(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 CHKERRQ(DMGetCoordinateDim(dm, &cdim)); 126 CHKERRQ(DMGetCoordinatesLocal(dm, &coordsLocal)); 127 CHKERRQ(VecGetArrayRead(coordsLocal, &coords)); 128 CHKERRQ(DMPlexPointLocalRead(cdm, cone[0], coords, &cA)); 129 CHKERRQ(DMPlexPointLocalRead(cdm, cone[1], coords, &cB)); 130 sbr->edgeLen[off] = DMPlex_DistD_Internal(cdim, cA, cB); 131 CHKERRQ(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 CHKERRQ(DMPlexTransformGetDM(tr, &dm)); 146 while (!PointQueueEmpty(queue)) { 147 PetscInt p = -1; 148 const PetscInt *support; 149 PetscInt supportSize, s; 150 151 CHKERRQ(PointQueueDequeue(queue, &p)); 152 CHKERRQ(DMPlexGetSupport(dm, p, &support)); 153 CHKERRQ(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 CHKERRQ(DMLabelGetValue(sbr->splitPoints, cell, &cval)); 162 if (cval == 2) continue; 163 CHKERRQ(DMPlexGetCone(dm, cell, &cone)); 164 CHKERRQ(DMPlexGetConeSize(dm, cell, &coneSize)); 165 CHKERRQ(SBRGetEdgeLen_Private(tr, cone[0], &maxlen)); 166 maxedge = cone[0]; 167 for (c = 1; c < coneSize; ++c) { 168 CHKERRQ(SBRGetEdgeLen_Private(tr, cone[c], &len)); 169 if (len > maxlen) {maxlen = len; maxedge = cone[c];} 170 } 171 CHKERRQ(DMLabelGetValue(sbr->splitPoints, maxedge, &eval)); 172 if (eval != 1) { 173 CHKERRQ(DMLabelSetValue(sbr->splitPoints, maxedge, 1)); 174 CHKERRQ(PointQueueEnqueue(queue, maxedge)); 175 } 176 CHKERRQ(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 CHKERRQ(DMPlexTransformGetDM(tr, &dm)); 194 /* Add in leaves */ 195 CHKERRQ(PetscSFGetGraph(pointSF, NULL, &Nl, &points, NULL)); 196 for (l = 0; l < Nl; ++l) { 197 CHKERRQ(DMLabelGetValue(splitPoints, points[l], &val)); 198 if (val > 0) splitArray[points[l]] = val; 199 } 200 /* Add in shared roots */ 201 CHKERRQ(PetscSFComputeDegreeBegin(pointSF, °ree)); 202 CHKERRQ(PetscSFComputeDegreeEnd(pointSF, °ree)); 203 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 204 for (p = pStart; p < pEnd; ++p) { 205 if (degree[p]) { 206 CHKERRQ(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 CHKERRQ(DMPlexTransformGetDM(tr, &dm)); 225 /* Read out leaves */ 226 CHKERRQ(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 CHKERRQ(DMLabelGetValue(splitPoints, p, &val)); 233 if (val <= 0) { 234 CHKERRQ(DMLabelSetValue(splitPoints, p, cval)); 235 CHKERRQ(PointQueueEnqueue(queue, p)); 236 } 237 } 238 } 239 /* Read out shared roots */ 240 CHKERRQ(PetscSFComputeDegreeBegin(pointSF, °ree)); 241 CHKERRQ(PetscSFComputeDegreeEnd(pointSF, °ree)); 242 CHKERRQ(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 CHKERRQ(DMLabelGetValue(splitPoints, p, &val)); 249 if (val <= 0) { 250 CHKERRQ(DMLabelSetValue(splitPoints, p, cval)); 251 CHKERRQ(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 CHKERRQ(DMPlexTransformGetDM(tr, &dm)); 295 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "Split Points", &sbr->splitPoints)); 296 /* Create edge lengths */ 297 CHKERRQ(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd)); 298 CHKERRQ(PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen)); 299 CHKERRQ(PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd)); 300 for (e = eStart; e < eEnd; ++e) { 301 CHKERRQ(PetscSectionSetDof(sbr->secEdgeLen, e, 1)); 302 } 303 CHKERRQ(PetscSectionSetUp(sbr->secEdgeLen)); 304 CHKERRQ(PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize)); 305 CHKERRQ(PetscCalloc1(edgeLenSize, &sbr->edgeLen)); 306 /* Add edges of cells that are marked for refinement to edge queue */ 307 CHKERRQ(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 CHKERRQ(PointQueueCreate(1024, &queue)); 310 CHKERRQ(DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS)); 311 CHKERRQ(DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc)); 312 if (refineIS) CHKERRQ(ISGetIndices(refineIS, &refineCells)); 313 for (c = 0; c < Nc; ++c) { 314 const PetscInt cell = refineCells[c]; 315 PetscInt depth; 316 317 CHKERRQ(DMPlexGetPointDepth(dm, cell, &depth)); 318 if (depth == 1) { 319 CHKERRQ(DMLabelSetValue(sbr->splitPoints, cell, 1)); 320 CHKERRQ(PointQueueEnqueue(queue, cell)); 321 } else { 322 PetscInt *closure = NULL; 323 PetscInt Ncl, cl; 324 325 CHKERRQ(DMLabelSetValue(sbr->splitPoints, cell, depth)); 326 CHKERRQ(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 CHKERRQ(DMLabelSetValue(sbr->splitPoints, edge, 1)); 332 CHKERRQ(PointQueueEnqueue(queue, edge)); 333 } 334 } 335 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure)); 336 } 337 } 338 if (refineIS) CHKERRQ(ISRestoreIndices(refineIS, &refineCells)); 339 CHKERRQ(ISDestroy(&refineIS)); 340 /* Setup communication */ 341 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 342 CHKERRQ(DMGetPointSF(dm, &pointSF)); 343 if (size > 1) { 344 PetscInt pStart, pEnd; 345 346 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 347 CHKERRQ(PetscCalloc1(pEnd-pStart, &sbr->splitArray)); 348 } 349 /* While edge queue is not empty: */ 350 empty = PointQueueEmpty(queue); 351 CHKERRMPI(MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) dm))); 352 while (!empty) { 353 CHKERRQ(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 CHKERRQ(SBRInitializeComm(tr, pointSF)); 368 CHKERRQ(PetscSFReduceBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX)); 369 CHKERRQ(PetscSFReduceEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX)); 370 CHKERRQ(PetscSFBcastBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE)); 371 CHKERRQ(PetscSFBcastEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE)); 372 CHKERRQ(SBRFinalizeComm(tr, pointSF, queue)); 373 } 374 empty = PointQueueEmpty(queue); 375 CHKERRMPI(MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) dm))); 376 } 377 CHKERRQ(PetscFree(sbr->splitArray)); 378 /* Calculate refineType for each cell */ 379 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType)); 380 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 381 for (p = pStart; p < pEnd; ++p) { 382 DMLabel trType = tr->trType; 383 DMPolytopeType ct; 384 PetscInt val; 385 386 CHKERRQ(DMPlexGetCellType(dm, p, &ct)); 387 switch (ct) { 388 case DM_POLYTOPE_POINT: 389 CHKERRQ(DMLabelSetValue(trType, p, RT_VERTEX));break; 390 case DM_POLYTOPE_SEGMENT: 391 CHKERRQ(DMLabelGetValue(sbr->splitPoints, p, &val)); 392 if (val == 1) CHKERRQ(DMLabelSetValue(trType, p, RT_EDGE_SPLIT)); 393 else CHKERRQ(DMLabelSetValue(trType, p, RT_EDGE)); 394 break; 395 case DM_POLYTOPE_TRIANGLE: 396 CHKERRQ(DMLabelGetValue(sbr->splitPoints, p, &val)); 397 if (val == 2) { 398 const PetscInt *cone; 399 PetscReal lens[3]; 400 PetscInt vals[3], i; 401 402 CHKERRQ(DMPlexGetCone(dm, p, &cone)); 403 for (i = 0; i < 3; ++i) { 404 CHKERRQ(DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i])); 405 vals[i] = vals[i] < 0 ? 0 : vals[i]; 406 CHKERRQ(SBRGetEdgeLen_Private(tr, cone[i], &lens[i])); 407 } 408 if (vals[0] && vals[1] && vals[2]) CHKERRQ(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT)); 409 else if (vals[0] && vals[1]) CHKERRQ(DMLabelSetValue(trType, p, lens[0] > lens[1] ? RT_TRIANGLE_SPLIT_01 : RT_TRIANGLE_SPLIT_10)); 410 else if (vals[1] && vals[2]) CHKERRQ(DMLabelSetValue(trType, p, lens[1] > lens[2] ? RT_TRIANGLE_SPLIT_12 : RT_TRIANGLE_SPLIT_21)); 411 else if (vals[2] && vals[0]) CHKERRQ(DMLabelSetValue(trType, p, lens[2] > lens[0] ? RT_TRIANGLE_SPLIT_20 : RT_TRIANGLE_SPLIT_02)); 412 else if (vals[0]) CHKERRQ(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_0)); 413 else if (vals[1]) CHKERRQ(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_1)); 414 else if (vals[2]) CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelGetValue(sbr->splitPoints, p, &val)); 421 } 422 /* Cleanup */ 423 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew)); 464 break; 465 default: CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt)); 661 break; 662 case DM_POLYTOPE_SEGMENT: 663 if (val == RT_EDGE) CHKERRQ(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt)); 664 else CHKERRQ(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: CHKERRQ(SBRGetTriangleSplitSingle(2, Nt, target, size, cone, ornt));break; 669 case RT_TRIANGLE_SPLIT_1: CHKERRQ(SBRGetTriangleSplitSingle(0, Nt, target, size, cone, ornt));break; 670 case RT_TRIANGLE_SPLIT_2: CHKERRQ(SBRGetTriangleSplitSingle(1, Nt, target, size, cone, ornt));break; 671 case RT_TRIANGLE_SPLIT_21: CHKERRQ(SBRGetTriangleSplitDouble(-3, Nt, target, size, cone, ornt));break; 672 case RT_TRIANGLE_SPLIT_10: CHKERRQ(SBRGetTriangleSplitDouble(-2, Nt, target, size, cone, ornt));break; 673 case RT_TRIANGLE_SPLIT_02: CHKERRQ(SBRGetTriangleSplitDouble(-1, Nt, target, size, cone, ornt));break; 674 case RT_TRIANGLE_SPLIT_12: CHKERRQ(SBRGetTriangleSplitDouble( 0, Nt, target, size, cone, ornt));break; 675 case RT_TRIANGLE_SPLIT_20: CHKERRQ(SBRGetTriangleSplitDouble( 1, Nt, target, size, cone, ornt));break; 676 case RT_TRIANGLE_SPLIT_01: CHKERRQ(SBRGetTriangleSplitDouble( 2, Nt, target, size, cone, ornt));break; 677 case RT_TRIANGLE_SPLIT: CHKERRQ(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt)); break; 678 default: CHKERRQ(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 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"DMPlex Options")); 694 CHKERRQ(PetscOptionsIntArray("-dm_plex_transform_sbr_ref_cell", "Mark cells for refinement", "", cells, &n, &flg)); 695 if (flg) { 696 DMLabel active; 697 698 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active)); 699 for (i = 0; i < n; ++i) CHKERRQ(DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE)); 700 CHKERRQ(DMPlexTransformSetActive(tr, active)); 701 CHKERRQ(DMLabelDestroy(&active)); 702 } 703 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii)); 715 if (isascii) { 716 PetscViewerFormat format; 717 const char *name; 718 719 CHKERRQ(PetscObjectGetName((PetscObject) tr, &name)); 720 CHKERRQ(PetscViewerASCIIPrintf(viewer, "SBR refinement %s\n", name ? name : "")); 721 CHKERRQ(PetscViewerGetFormat(viewer, &format)); 722 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 723 CHKERRQ(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 CHKERRQ(PetscFree(sbr->edgeLen)); 737 CHKERRQ(PetscSectionDestroy(&sbr->secEdgeLen)); 738 CHKERRQ(DMLabelDestroy(&sbr->splitPoints)); 739 CHKERRQ(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 CHKERRQ(PetscNewLog(tr, &f)); 763 tr->data = f; 764 765 CHKERRQ(DMPlexTransformInitialize_SBR(tr)); 766 CHKERRQ(PetscCitationsRegister(SBRCitation, &SBRcite)); 767 PetscFunctionReturn(0); 768 } 769