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