Actual source code: plexrefsbr.c

  1: #include <petsc/private/dmplextransformimpl.h>
  2: #include <petscsf.h>

  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";

 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: };

 24: static PetscErrorCode PointQueueCreate(PetscInt size, PointQueue *queue)
 25: {
 26:   PointQueue     q;

 30:   if (size < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Queue size %D must be non-negative", size);
 31:   PetscCalloc1(1, &q);
 32:   q->size = size;
 33:   PetscMalloc1(q->size, &q->points);
 34:   q->num   = 0;
 35:   q->front = 0;
 36:   q->back  = q->size-1;
 37:   *queue = q;
 38:   return(0);
 39: }

 41: static PetscErrorCode PointQueueDestroy(PointQueue *queue)
 42: {
 43:   PointQueue     q = *queue;

 47:   PetscFree(q->points);
 48:   PetscFree(q);
 49:   *queue = NULL;
 50:   return(0);
 51: }

 53: static PetscErrorCode PointQueueEnsureSize(PointQueue queue)
 54: {

 58:   if (queue->num < queue->size) return(0);
 59:   queue->size *= 2;
 60:   PetscRealloc(queue->size * sizeof(PetscInt), &queue->points);
 61:   return(0);
 62: }

 64: static PetscErrorCode PointQueueEnqueue(PointQueue queue, PetscInt p)
 65: {

 69:   PointQueueEnsureSize(queue);
 70:   queue->back = (queue->back + 1) % queue->size;
 71:   queue->points[queue->back] = p;
 72:   ++queue->num;
 73:   return(0);
 74: }

 76: static PetscErrorCode PointQueueDequeue(PointQueue queue, PetscInt *p)
 77: {
 79:   if (!queue->num) SETERRQ(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:   return(0);
 84: }

 86: #if 0
 87: static PetscErrorCode PointQueueFront(PointQueue queue, PetscInt *p)
 88: {
 90:   if (!queue->num) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot get the front of an empty queue");
 91:   *p = queue->points[queue->front];
 92:   return(0);
 93: }

 95: static PetscErrorCode PointQueueBack(PointQueue queue, PetscInt *p)
 96: {
 98:   if (!queue->num) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot get the back of an empty queue");
 99:   *p = queue->points[queue->back];
100:   return(0);
101: }
102: #endif

104: PETSC_STATIC_INLINE PetscBool PointQueueEmpty(PointQueue queue)
105: {
106:   if (!queue->num) return PETSC_TRUE;
107:   return PETSC_FALSE;
108: }

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;

118:   DMPlexTransformGetDM(tr, &dm);
119:   PetscSectionGetOffset(sbr->secEdgeLen, edge, &off);
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;

128:     DMGetCoordinateDM(dm, &cdm);
129:     DMPlexGetCone(dm, edge, &cone);
130:     DMPlexGetConeSize(dm, edge, &coneSize);
131:     if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Edge %D cone size must be 2, not %D", edge, coneSize);
132:     DMGetCoordinateDim(dm, &cdim);
133:     DMGetCoordinatesLocal(dm, &coordsLocal);
134:     VecGetArrayRead(coordsLocal, &coords);
135:     DMPlexPointLocalRead(cdm, cone[0], coords, &cA);
136:     DMPlexPointLocalRead(cdm, cone[1], coords, &cB);
137:     sbr->edgeLen[off] = DMPlex_DistD_Internal(cdim, cA, cB);
138:     VecRestoreArrayRead(coordsLocal, &coords);
139:   }
140:   *len = sbr->edgeLen[off];
141:   return(0);
142: }

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;

153:   DMPlexTransformGetDM(tr, &dm);
154:   while (!PointQueueEmpty(queue)) {
155:     PetscInt        p = -1;
156:     const PetscInt *support;
157:     PetscInt        supportSize, s;

159:     PointQueueDequeue(queue, &p);
160:     DMPlexGetSupport(dm, p, &support);
161:     DMPlexGetSupportSize(dm, p, &supportSize);
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;

169:       DMLabelGetValue(sbr->splitPoints, cell, &cval);
170:       if (cval == 2) continue;
171:       DMPlexGetCone(dm, cell, &cone);
172:       DMPlexGetConeSize(dm, cell, &coneSize);
173:       SBRGetEdgeLen_Private(tr, cone[0], &maxlen);
174:       maxedge = cone[0];
175:       for (c = 1; c < coneSize; ++c) {
176:         SBRGetEdgeLen_Private(tr, cone[c], &len);
177:         if (len > maxlen) {maxlen = len; maxedge = cone[c];}
178:       }
179:       DMLabelGetValue(sbr->splitPoints, maxedge, &eval);
180:       if (eval != 1) {
181:         DMLabelSetValue(sbr->splitPoints, maxedge, 1);
182:         PointQueueEnqueue(queue, maxedge);
183:       }
184:       DMLabelSetValue(sbr->splitPoints, cell, 2);
185:     }
186:   }
187:   return(0);
188: }

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;

202:   DMPlexTransformGetDM(tr, &dm);
203:   /* Add in leaves */
204:   PetscSFGetGraph(pointSF, NULL, &Nl, &points, NULL);
205:   for (l = 0; l < Nl; ++l) {
206:     DMLabelGetValue(splitPoints, points[l], &val);
207:     if (val > 0) splitArray[points[l]] = val;
208:   }
209:   /* Add in shared roots */
210:   PetscSFComputeDegreeBegin(pointSF, &degree);
211:   PetscSFComputeDegreeEnd(pointSF, &degree);
212:   DMPlexGetChart(dm, &pStart, &pEnd);
213:   for (p = pStart; p < pEnd; ++p) {
214:     if (degree[p]) {
215:       DMLabelGetValue(splitPoints, p, &val);
216:       if (val > 0) splitArray[p] = val;
217:     }
218:   }
219:   return(0);
220: }

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;

234:   DMPlexTransformGetDM(tr, &dm);
235:   /* Read out leaves */
236:   PetscSFGetGraph(pointSF, NULL, &Nl, &points, NULL);
237:   for (l = 0; l < Nl; ++l) {
238:     const PetscInt p    = points[l];
239:     const PetscInt cval = splitArray[p];

241:     if (cval) {
242:       DMLabelGetValue(splitPoints, p, &val);
243:       if (val <= 0) {
244:         DMLabelSetValue(splitPoints, p, cval);
245:         PointQueueEnqueue(queue, p);
246:       }
247:     }
248:   }
249:   /* Read out shared roots */
250:   PetscSFComputeDegreeBegin(pointSF, &degree);
251:   PetscSFComputeDegreeEnd(pointSF, &degree);
252:   DMPlexGetChart(dm, &pStart, &pEnd);
253:   for (p = pStart; p < pEnd; ++p) {
254:     if (degree[p]) {
255:       const PetscInt cval = splitArray[p];

257:       if (cval) {
258:         DMLabelGetValue(splitPoints, p, &val);
259:         if (val <= 0) {
260:           DMLabelSetValue(splitPoints, p, cval);
261:           PointQueueEnqueue(queue, p);
262:         }
263:       }
264:     }
265:   }
266:   return(0);
267: }

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:

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;

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;

305:   DMPlexTransformGetDM(tr, &dm);
306:   DMLabelCreate(PETSC_COMM_SELF, "Split Points", &sbr->splitPoints);
307:   /* Create edge lengths */
308:   DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
309:   PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen);
310:   PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd);
311:   for (e = eStart; e < eEnd; ++e) {
312:     PetscSectionSetDof(sbr->secEdgeLen, e, 1);
313:   }
314:   PetscSectionSetUp(sbr->secEdgeLen);
315:   PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize);
316:   PetscCalloc1(edgeLenSize, &sbr->edgeLen);
317:   /* Add edges of cells that are marked for refinement to edge queue */
318:   DMPlexTransformGetActive(tr, &active);
319:   if (!active) SETERRQ(PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_WRONGSTATE, "DMPlexTransform must have an adaptation label in order to use SBR algorithm");
320:   PointQueueCreate(1024, &queue);
321:   DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS);
322:   DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc);
323:   if (refineIS) {ISGetIndices(refineIS, &refineCells);}
324:   for (c = 0; c < Nc; ++c) {
325:     const PetscInt cell = refineCells[c];
326:     PetscInt       depth;

328:     DMPlexGetPointDepth(dm, cell, &depth);
329:     if (depth == 1) {
330:       DMLabelSetValue(sbr->splitPoints, cell, 1);
331:       PointQueueEnqueue(queue, cell);
332:     } else {
333:       PetscInt *closure = NULL;
334:       PetscInt  Ncl, cl;

336:       DMLabelSetValue(sbr->splitPoints, cell, depth);
337:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
338:       for (cl = 0; cl < Ncl; cl += 2) {
339:         const PetscInt edge = closure[cl];

341:         if (edge >= eStart && edge < eEnd) {
342:           DMLabelSetValue(sbr->splitPoints, edge, 1);
343:           PointQueueEnqueue(queue, edge);
344:         }
345:       }
346:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
347:     }
348:   }
349:   if (refineIS) {ISRestoreIndices(refineIS, &refineCells);}
350:   ISDestroy(&refineIS);
351:   /* Setup communication */
352:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
353:   DMGetPointSF(dm, &pointSF);
354:   if (size > 1) {
355:     PetscInt pStart, pEnd;

357:     DMPlexGetChart(dm, &pStart, &pEnd);
358:     PetscCalloc1(pEnd-pStart, &sbr->splitArray);
359:   }
360:   /* While edge queue is not empty: */
361:   empty = PointQueueEmpty(queue);
362:   MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) dm));
363:   while (!empty) {
364:     SBRSplitLocalEdges_Private(tr, queue);
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.

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.

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:       SBRInitializeComm(tr, pointSF);
379:       PetscSFReduceBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX);
380:       PetscSFReduceEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX);
381:       PetscSFBcastBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE);
382:       PetscSFBcastEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE);
383:       SBRFinalizeComm(tr, pointSF, queue);
384:     }
385:     empty = PointQueueEmpty(queue);
386:     MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) dm));
387:   }
388:   PetscFree(sbr->splitArray);
389:   /* Calculate refineType for each cell */
390:   DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType);
391:   DMPlexGetChart(dm, &pStart, &pEnd);
392:   for (p = pStart; p < pEnd; ++p) {
393:     DMLabel        trType = tr->trType;
394:     DMPolytopeType ct;
395:     PetscInt       val;

397:     DMPlexGetCellType(dm, p, &ct);
398:     switch (ct) {
399:       case DM_POLYTOPE_POINT:
400:         DMLabelSetValue(trType, p, RT_VERTEX);break;
401:       case DM_POLYTOPE_SEGMENT:
402:         DMLabelGetValue(sbr->splitPoints, p, &val);
403:         if (val == 1) {DMLabelSetValue(trType, p, RT_EDGE_SPLIT);}
404:         else          {DMLabelSetValue(trType, p, RT_EDGE);}
405:         break;
406:       case DM_POLYTOPE_TRIANGLE:
407:         DMLabelGetValue(sbr->splitPoints, p, &val);
408:         if (val == 2) {
409:           const PetscInt *cone;
410:           PetscReal       lens[3];
411:           PetscInt        vals[3], i;

413:           DMPlexGetCone(dm, p, &cone);
414:           for (i = 0; i < 3; ++i) {
415:             DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i]);
416:             vals[i] = vals[i] < 0 ? 0 : vals[i];
417:             SBRGetEdgeLen_Private(tr, cone[i], &lens[i]);
418:           }
419:           if (vals[0] && vals[1] && vals[2]) {DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT);}
420:           else if (vals[0] && vals[1])       {DMLabelSetValue(trType, p, lens[0] > lens[1] ? RT_TRIANGLE_SPLIT_01 : RT_TRIANGLE_SPLIT_10);}
421:           else if (vals[1] && vals[2])       {DMLabelSetValue(trType, p, lens[1] > lens[2] ? RT_TRIANGLE_SPLIT_12 : RT_TRIANGLE_SPLIT_21);}
422:           else if (vals[2] && vals[0])       {DMLabelSetValue(trType, p, lens[2] > lens[0] ? RT_TRIANGLE_SPLIT_20 : RT_TRIANGLE_SPLIT_02);}
423:           else if (vals[0])                  {DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_0);}
424:           else if (vals[1])                  {DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_1);}
425:           else if (vals[2])                  {DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_2);}
426:           else SETERRQ4(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 {DMLabelSetValue(trType, p, RT_TRIANGLE);}
428:         break;
429:       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
430:     }
431:     DMLabelGetValue(sbr->splitPoints, p, &val);
432:   }
433:   /* Cleanup */
434:   PointQueueDestroy(&queue);
435:   return(0);
436: }

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;

444:   DMLabelGetValue(tr->trType, sp, &rt);
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:       DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew);
476:       break;
477:     default: DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
478:   }
479:   return(0);
480: }

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};

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:   return(0);
519: }

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};

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:   return(0);
650: }

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;

659:   if (p < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point argument is invalid");
660:   DMLabelGetValue(trType, p, &val);
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:       DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
674:       break;
675:     case DM_POLYTOPE_SEGMENT:
676:       if (val == RT_EDGE) {DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);}
677:       else                {DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt);}
678:       break;
679:     case DM_POLYTOPE_TRIANGLE:
680:       switch (val) {
681:         case RT_TRIANGLE_SPLIT_0: SBRGetTriangleSplitSingle(2, Nt, target, size, cone, ornt);break;
682:         case RT_TRIANGLE_SPLIT_1: SBRGetTriangleSplitSingle(0, Nt, target, size, cone, ornt);break;
683:         case RT_TRIANGLE_SPLIT_2: SBRGetTriangleSplitSingle(1, Nt, target, size, cone, ornt);break;
684:         case RT_TRIANGLE_SPLIT_21: SBRGetTriangleSplitDouble(-3, Nt, target, size, cone, ornt);break;
685:         case RT_TRIANGLE_SPLIT_10: SBRGetTriangleSplitDouble(-2, Nt, target, size, cone, ornt);break;
686:         case RT_TRIANGLE_SPLIT_02: SBRGetTriangleSplitDouble(-1, Nt, target, size, cone, ornt);break;
687:         case RT_TRIANGLE_SPLIT_12: SBRGetTriangleSplitDouble( 0, Nt, target, size, cone, ornt);break;
688:         case RT_TRIANGLE_SPLIT_20: SBRGetTriangleSplitDouble( 1, Nt, target, size, cone, ornt);break;
689:         case RT_TRIANGLE_SPLIT_01: SBRGetTriangleSplitDouble( 2, Nt, target, size, cone, ornt);break;
690:         case RT_TRIANGLE_SPLIT: DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt); break;
691:         default: DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
692:       }
693:       break;
694:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
695:   }
696:   return(0);
697: }

699: static PetscErrorCode DMPlexTransformSetFromOptions_SBR(PetscOptionItems *PetscOptionsObject, DMPlexTransform tr)
700: {
701:   PetscInt       cells[256], n = 256, i;
702:   PetscBool      flg;

707:   PetscOptionsHead(PetscOptionsObject,"DMPlex Options");
708:   PetscOptionsIntArray("-dm_plex_transform_sbr_ref_cell", "Mark cells for refinement", "", cells, &n, &flg);
709:   if (flg) {
710:     DMLabel active;

712:     DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active);
713:     for (i = 0; i < n; ++i) {DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE);}
714:     DMPlexTransformSetActive(tr, active);
715:     DMLabelDestroy(&active);
716:   }
717:   PetscOptionsTail();
718:   return(0);
719: }

721: static PetscErrorCode DMPlexTransformView_SBR(DMPlexTransform tr, PetscViewer viewer)
722: {
723:   PetscBool      isascii;

729:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
730:   if (isascii) {
731:     PetscViewerFormat format;
732:     const char       *name;

734:     PetscObjectGetName((PetscObject) tr, &name);
735:     PetscViewerASCIIPrintf(viewer, "SBR refinement %s\n", name ? name : "");
736:     PetscViewerGetFormat(viewer, &format);
737:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
738:       DMLabelView(tr->trType, viewer);
739:     }
740:   } else {
741:     SETERRQ1(PetscObjectComm((PetscObject) tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject) viewer)->type_name);
742:   }
743:   return(0);
744: }

746: static PetscErrorCode DMPlexTransformDestroy_SBR(DMPlexTransform tr)
747: {
748:   DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
749:   PetscErrorCode    ierr;

752:   PetscFree(sbr->edgeLen);
753:   PetscSectionDestroy(&sbr->secEdgeLen);
754:   DMLabelDestroy(&sbr->splitPoints);
755:   PetscFree(tr->data);
756:   return(0);
757: }

759: static PetscErrorCode DMPlexTransformInitialize_SBR(DMPlexTransform tr)
760: {
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:   return(0);
770: }

772: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform tr)
773: {
774:   DMPlexRefine_SBR *f;
775:   PetscErrorCode    ierr;

779:   PetscNewLog(tr, &f);
780:   tr->data = f;

782:   DMPlexTransformInitialize_SBR(tr);
783:   PetscCitationsRegister(SBRCitation, &SBRcite);
784:   return(0);
785: }