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;

 29:   PetscCalloc1(1, &q);
 30:   q->size = size;
 31:   PetscMalloc1(q->size, &q->points);
 32:   q->num   = 0;
 33:   q->front = 0;
 34:   q->back  = q->size-1;
 35:   *queue = q;
 36:   return 0;
 37: }

 39: static PetscErrorCode PointQueueDestroy(PointQueue *queue)
 40: {
 41:   PointQueue     q = *queue;

 43:   PetscFree(q->points);
 44:   PetscFree(q);
 45:   *queue = NULL;
 46:   return 0;
 47: }

 49: static PetscErrorCode PointQueueEnsureSize(PointQueue queue)
 50: {
 51:   if (queue->num < queue->size) return 0;
 52:   queue->size *= 2;
 53:   PetscRealloc(queue->size * sizeof(PetscInt), &queue->points);
 54:   return 0;
 55: }

 57: static PetscErrorCode PointQueueEnqueue(PointQueue queue, PetscInt p)
 58: {
 59:   PointQueueEnsureSize(queue);
 60:   queue->back = (queue->back + 1) % queue->size;
 61:   queue->points[queue->back] = p;
 62:   ++queue->num;
 63:   return 0;
 64: }

 66: static PetscErrorCode PointQueueDequeue(PointQueue queue, PetscInt *p)
 67: {
 69:   *p = queue->points[queue->front];
 70:   queue->front = (queue->front + 1) % queue->size;
 71:   --queue->num;
 72:   return 0;
 73: }

 75: #if 0
 76: static PetscErrorCode PointQueueFront(PointQueue queue, PetscInt *p)
 77: {
 79:   *p = queue->points[queue->front];
 80:   return 0;
 81: }

 83: static PetscErrorCode PointQueueBack(PointQueue queue, PetscInt *p)
 84: {
 86:   *p = queue->points[queue->back];
 87:   return 0;
 88: }
 89: #endif

 91: static inline PetscBool PointQueueEmpty(PointQueue queue)
 92: {
 93:   if (!queue->num) return PETSC_TRUE;
 94:   return PETSC_FALSE;
 95: }

 97: static PetscErrorCode SBRGetEdgeLen_Private(DMPlexTransform tr, PetscInt edge, PetscReal *len)
 98: {
 99:   DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
100:   DM                dm;
101:   PetscInt          off;

104:   DMPlexTransformGetDM(tr, &dm);
105:   PetscSectionGetOffset(sbr->secEdgeLen, edge, &off);
106:   if (sbr->edgeLen[off] <= 0.0) {
107:     DM                 cdm;
108:     Vec                coordsLocal;
109:     const PetscScalar *coords;
110:     const PetscInt    *cone;
111:     PetscScalar       *cA, *cB;
112:     PetscInt           coneSize, cdim;

114:     DMGetCoordinateDM(dm, &cdm);
115:     DMPlexGetCone(dm, edge, &cone);
116:     DMPlexGetConeSize(dm, edge, &coneSize);
118:     DMGetCoordinateDim(dm, &cdim);
119:     DMGetCoordinatesLocal(dm, &coordsLocal);
120:     VecGetArrayRead(coordsLocal, &coords);
121:     DMPlexPointLocalRead(cdm, cone[0], coords, &cA);
122:     DMPlexPointLocalRead(cdm, cone[1], coords, &cB);
123:     sbr->edgeLen[off] = DMPlex_DistD_Internal(cdim, cA, cB);
124:     VecRestoreArrayRead(coordsLocal, &coords);
125:   }
126:   *len = sbr->edgeLen[off];
127:   return 0;
128: }

130: /* Mark local edges that should be split */
131: /* TODO This will not work in 3D */
132: static PetscErrorCode SBRSplitLocalEdges_Private(DMPlexTransform tr, PointQueue queue)
133: {
134:   DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
135:   DM                dm;

137:   DMPlexTransformGetDM(tr, &dm);
138:   while (!PointQueueEmpty(queue)) {
139:     PetscInt        p = -1;
140:     const PetscInt *support;
141:     PetscInt        supportSize, s;

143:     PointQueueDequeue(queue, &p);
144:     DMPlexGetSupport(dm, p, &support);
145:     DMPlexGetSupportSize(dm, p, &supportSize);
146:     for (s = 0; s < supportSize; ++s) {
147:       const PetscInt  cell = support[s];
148:       const PetscInt *cone;
149:       PetscInt        coneSize, c;
150:       PetscInt        cval, eval, maxedge;
151:       PetscReal       len, maxlen;

153:       DMLabelGetValue(sbr->splitPoints, cell, &cval);
154:       if (cval == 2) continue;
155:       DMPlexGetCone(dm, cell, &cone);
156:       DMPlexGetConeSize(dm, cell, &coneSize);
157:       SBRGetEdgeLen_Private(tr, cone[0], &maxlen);
158:       maxedge = cone[0];
159:       for (c = 1; c < coneSize; ++c) {
160:         SBRGetEdgeLen_Private(tr, cone[c], &len);
161:         if (len > maxlen) {maxlen = len; maxedge = cone[c];}
162:       }
163:       DMLabelGetValue(sbr->splitPoints, maxedge, &eval);
164:       if (eval != 1) {
165:         DMLabelSetValue(sbr->splitPoints, maxedge, 1);
166:         PointQueueEnqueue(queue, maxedge);
167:       }
168:       DMLabelSetValue(sbr->splitPoints, cell, 2);
169:     }
170:   }
171:   return 0;
172: }

174: static PetscErrorCode SBRInitializeComm(DMPlexTransform tr, PetscSF pointSF)
175: {
176:   DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
177:   DM                dm;
178:   DMLabel           splitPoints = sbr->splitPoints;
179:   PetscInt         *splitArray  = sbr->splitArray;
180:   const PetscInt   *degree;
181:   const PetscInt   *points;
182:   PetscInt          Nl, l, pStart, pEnd, p, val;

184:   DMPlexTransformGetDM(tr, &dm);
185:   /* Add in leaves */
186:   PetscSFGetGraph(pointSF, NULL, &Nl, &points, NULL);
187:   for (l = 0; l < Nl; ++l) {
188:     DMLabelGetValue(splitPoints, points[l], &val);
189:     if (val > 0) splitArray[points[l]] = val;
190:   }
191:   /* Add in shared roots */
192:   PetscSFComputeDegreeBegin(pointSF, &degree);
193:   PetscSFComputeDegreeEnd(pointSF, &degree);
194:   DMPlexGetChart(dm, &pStart, &pEnd);
195:   for (p = pStart; p < pEnd; ++p) {
196:     if (degree[p]) {
197:       DMLabelGetValue(splitPoints, p, &val);
198:       if (val > 0) splitArray[p] = val;
199:     }
200:   }
201:   return 0;
202: }

204: static PetscErrorCode SBRFinalizeComm(DMPlexTransform tr, PetscSF pointSF, PointQueue queue)
205: {
206:   DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
207:   DM                dm;
208:   DMLabel           splitPoints = sbr->splitPoints;
209:   PetscInt         *splitArray  = sbr->splitArray;
210:   const PetscInt   *degree;
211:   const PetscInt   *points;
212:   PetscInt          Nl, l, pStart, pEnd, p, val;

214:   DMPlexTransformGetDM(tr, &dm);
215:   /* Read out leaves */
216:   PetscSFGetGraph(pointSF, NULL, &Nl, &points, NULL);
217:   for (l = 0; l < Nl; ++l) {
218:     const PetscInt p    = points[l];
219:     const PetscInt cval = splitArray[p];

221:     if (cval) {
222:       DMLabelGetValue(splitPoints, p, &val);
223:       if (val <= 0) {
224:         DMLabelSetValue(splitPoints, p, cval);
225:         PointQueueEnqueue(queue, p);
226:       }
227:     }
228:   }
229:   /* Read out shared roots */
230:   PetscSFComputeDegreeBegin(pointSF, &degree);
231:   PetscSFComputeDegreeEnd(pointSF, &degree);
232:   DMPlexGetChart(dm, &pStart, &pEnd);
233:   for (p = pStart; p < pEnd; ++p) {
234:     if (degree[p]) {
235:       const PetscInt cval = splitArray[p];

237:       if (cval) {
238:         DMLabelGetValue(splitPoints, p, &val);
239:         if (val <= 0) {
240:           DMLabelSetValue(splitPoints, p, cval);
241:           PointQueueEnqueue(queue, p);
242:         }
243:       }
244:     }
245:   }
246:   return 0;
247: }

249: /*
250:   The 'splitPoints' label marks mesh points to be divided. It marks edges with 1, triangles with 2, and tetrahedra with 3.
251:   Then the refinement type is calculated as follows:

253:     vertex:                   0
254:     edge unsplit:             1
255:     edge split:               2
256:     triangle unsplit:         3
257:     triangle split all edges: 4
258:     triangle split edges 0 1: 5
259:     triangle split edges 1 0: 6
260:     triangle split edges 1 2: 7
261:     triangle split edges 2 1: 8
262:     triangle split edges 2 0: 9
263:     triangle split edges 0 2: 10
264:     triangle split edge 0:    11
265:     triangle split edge 1:    12
266:     triangle split edge 2:    13
267: */
268: 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;

270: static PetscErrorCode DMPlexTransformSetUp_SBR(DMPlexTransform tr)
271: {
272:   DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
273:   DM                dm;
274:   DMLabel           active;
275:   PetscSF           pointSF;
276:   PointQueue        queue = NULL;
277:   IS                refineIS;
278:   const PetscInt   *refineCells;
279:   PetscMPIInt       size;
280:   PetscInt          pStart, pEnd, p, eStart, eEnd, e, edgeLenSize, Nc, c;
281:   PetscBool         empty;

283:   DMPlexTransformGetDM(tr, &dm);
284:   DMLabelCreate(PETSC_COMM_SELF, "Split Points", &sbr->splitPoints);
285:   /* Create edge lengths */
286:   DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
287:   PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen);
288:   PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd);
289:   for (e = eStart; e < eEnd; ++e) {
290:     PetscSectionSetDof(sbr->secEdgeLen, e, 1);
291:   }
292:   PetscSectionSetUp(sbr->secEdgeLen);
293:   PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize);
294:   PetscCalloc1(edgeLenSize, &sbr->edgeLen);
295:   /* Add edges of cells that are marked for refinement to edge queue */
296:   DMPlexTransformGetActive(tr, &active);
298:   PointQueueCreate(1024, &queue);
299:   DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS);
300:   DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc);
301:   if (refineIS) ISGetIndices(refineIS, &refineCells);
302:   for (c = 0; c < Nc; ++c) {
303:     const PetscInt cell = refineCells[c];
304:     PetscInt       depth;

306:     DMPlexGetPointDepth(dm, cell, &depth);
307:     if (depth == 1) {
308:       DMLabelSetValue(sbr->splitPoints, cell, 1);
309:       PointQueueEnqueue(queue, cell);
310:     } else {
311:       PetscInt *closure = NULL;
312:       PetscInt  Ncl, cl;

314:       DMLabelSetValue(sbr->splitPoints, cell, depth);
315:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
316:       for (cl = 0; cl < Ncl; cl += 2) {
317:         const PetscInt edge = closure[cl];

319:         if (edge >= eStart && edge < eEnd) {
320:           DMLabelSetValue(sbr->splitPoints, edge, 1);
321:           PointQueueEnqueue(queue, edge);
322:         }
323:       }
324:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
325:     }
326:   }
327:   if (refineIS) ISRestoreIndices(refineIS, &refineCells);
328:   ISDestroy(&refineIS);
329:   /* Setup communication */
330:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
331:   DMGetPointSF(dm, &pointSF);
332:   if (size > 1) {
333:     PetscInt pStart, pEnd;

335:     DMPlexGetChart(dm, &pStart, &pEnd);
336:     PetscCalloc1(pEnd-pStart, &sbr->splitArray);
337:   }
338:   /* While edge queue is not empty: */
339:   empty = PointQueueEmpty(queue);
340:   MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) dm));
341:   while (!empty) {
342:     SBRSplitLocalEdges_Private(tr, queue);
343:     /* Communicate marked edges
344:          An easy implementation is to allocate an array the size of the number of points. We put the splitPoints marks into the
345:          array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.

347:          TODO: We could use in-place communication with a different SF
348:            We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
349:            already been marked. If not, it might have been handled on the process in this round, but we add it anyway.

351:            In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
352:            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
353:            edge to the queue.
354:     */
355:     if (size > 1) {
356:       SBRInitializeComm(tr, pointSF);
357:       PetscSFReduceBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX);
358:       PetscSFReduceEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX);
359:       PetscSFBcastBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE);
360:       PetscSFBcastEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE);
361:       SBRFinalizeComm(tr, pointSF, queue);
362:     }
363:     empty = PointQueueEmpty(queue);
364:     MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) dm));
365:   }
366:   PetscFree(sbr->splitArray);
367:   /* Calculate refineType for each cell */
368:   DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType);
369:   DMPlexGetChart(dm, &pStart, &pEnd);
370:   for (p = pStart; p < pEnd; ++p) {
371:     DMLabel        trType = tr->trType;
372:     DMPolytopeType ct;
373:     PetscInt       val;

375:     DMPlexGetCellType(dm, p, &ct);
376:     switch (ct) {
377:       case DM_POLYTOPE_POINT:
378:         DMLabelSetValue(trType, p, RT_VERTEX);break;
379:       case DM_POLYTOPE_SEGMENT:
380:         DMLabelGetValue(sbr->splitPoints, p, &val);
381:         if (val == 1) DMLabelSetValue(trType, p, RT_EDGE_SPLIT);
382:         else          DMLabelSetValue(trType, p, RT_EDGE);
383:         break;
384:       case DM_POLYTOPE_TRIANGLE:
385:         DMLabelGetValue(sbr->splitPoints, p, &val);
386:         if (val == 2) {
387:           const PetscInt *cone;
388:           PetscReal       lens[3];
389:           PetscInt        vals[3], i;

391:           DMPlexGetCone(dm, p, &cone);
392:           for (i = 0; i < 3; ++i) {
393:             DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i]);
394:             vals[i] = vals[i] < 0 ? 0 : vals[i];
395:             SBRGetEdgeLen_Private(tr, cone[i], &lens[i]);
396:           }
397:           if (vals[0] && vals[1] && vals[2]) DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT);
398:           else if (vals[0] && vals[1])       DMLabelSetValue(trType, p, lens[0] > lens[1] ? RT_TRIANGLE_SPLIT_01 : RT_TRIANGLE_SPLIT_10);
399:           else if (vals[1] && vals[2])       DMLabelSetValue(trType, p, lens[1] > lens[2] ? RT_TRIANGLE_SPLIT_12 : RT_TRIANGLE_SPLIT_21);
400:           else if (vals[2] && vals[0])       DMLabelSetValue(trType, p, lens[2] > lens[0] ? RT_TRIANGLE_SPLIT_20 : RT_TRIANGLE_SPLIT_02);
401:           else if (vals[0])                  DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_0);
402:           else if (vals[1])                  DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_1);
403:           else if (vals[2])                  DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_2);
404:           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]);
405:         } else DMLabelSetValue(trType, p, RT_TRIANGLE);
406:         break;
407:       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
408:     }
409:     DMLabelGetValue(sbr->splitPoints, p, &val);
410:   }
411:   /* Cleanup */
412:   PointQueueDestroy(&queue);
413:   return 0;
414: }

416: static PetscErrorCode DMPlexTransformGetSubcellOrientation_SBR(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
417: {
418:   PetscInt         rt;

421:   DMLabelGetValue(tr->trType, sp, &rt);
422:   *rnew = r;
423:   *onew = o;
424:   switch (rt) {
425:     case RT_TRIANGLE_SPLIT_01:
426:     case RT_TRIANGLE_SPLIT_10:
427:     case RT_TRIANGLE_SPLIT_12:
428:     case RT_TRIANGLE_SPLIT_21:
429:     case RT_TRIANGLE_SPLIT_20:
430:     case RT_TRIANGLE_SPLIT_02:
431:       switch (tct) {
432:         case DM_POLYTOPE_SEGMENT:  break;
433:         case DM_POLYTOPE_TRIANGLE: break;
434:         default: break;
435:       }
436:       break;
437:     case RT_TRIANGLE_SPLIT_0:
438:     case RT_TRIANGLE_SPLIT_1:
439:     case RT_TRIANGLE_SPLIT_2:
440:       switch (tct) {
441:         case DM_POLYTOPE_SEGMENT:
442:           break;
443:         case DM_POLYTOPE_TRIANGLE:
444:           *onew = so < 0 ? -(o+1)  : o;
445:           *rnew = so < 0 ? (r+1)%2 : r;
446:           break;
447:         default: break;
448:       }
449:       break;
450:     case RT_EDGE_SPLIT:
451:     case RT_TRIANGLE_SPLIT:
452:       DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew);
453:       break;
454:     default: DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
455:   }
456:   return 0;
457: }

459: /* Add 1 edge inside this triangle, making 2 new triangles.
460:  2
461:  |\
462:  | \
463:  |  \
464:  |   \
465:  |    1
466:  |     \
467:  |  B   \
468:  2       1
469:  |      / \
470:  | ____/   0
471:  |/    A    \
472:  0-----0-----1
473: */
474: static PetscErrorCode SBRGetTriangleSplitSingle(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
475: {
476:   const PetscInt       *arr     = DMPolytopeTypeGetArrangment(DM_POLYTOPE_TRIANGLE, o);
477:   static DMPolytopeType triT1[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
478:   static PetscInt       triS1[] = {1, 2};
479:   static PetscInt       triC1[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
480:                                    DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0,    0,
481:                                    DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0,    0};
482:   static PetscInt       triO1[] = {0, 0,
483:                                    0,  0, -1,
484:                                    0,  0,  0};

487:   /* To get the other divisions, we reorient the triangle */
488:   triC1[2]  = arr[0*2];
489:   triC1[7]  = arr[1*2];
490:   triC1[11] = arr[0*2];
491:   triC1[15] = arr[1*2];
492:   triC1[22] = arr[1*2];
493:   triC1[26] = arr[2*2];
494:   *Nt = 2; *target = triT1; *size = triS1; *cone = triC1; *ornt = triO1;
495:   return 0;
496: }

498: /* Add 2 edges inside this triangle, making 3 new triangles.
499:  RT_TRIANGLE_SPLIT_12
500:  2
501:  |\
502:  | \
503:  |  \
504:  0   \
505:  |    1
506:  |     \
507:  |  B   \
508:  2-------1
509:  |   C  / \
510:  1 ____/   0
511:  |/    A    \
512:  0-----0-----1
513:  RT_TRIANGLE_SPLIT_10
514:  2
515:  |\
516:  | \
517:  |  \
518:  0   \
519:  |    1
520:  |     \
521:  |  A   \
522:  2       1
523:  |      /|\
524:  1 ____/ / 0
525:  |/ C   / B \
526:  0-----0-----1
527:  RT_TRIANGLE_SPLIT_20
528:  2
529:  |\
530:  | \
531:  |  \
532:  0   \
533:  |    \
534:  |     \
535:  |      \
536:  2   A   1
537:  |\       \
538:  1 ---\    \
539:  |B \_C----\\
540:  0-----0-----1
541:  RT_TRIANGLE_SPLIT_21
542:  2
543:  |\
544:  | \
545:  |  \
546:  0   \
547:  |    \
548:  |  B  \
549:  |      \
550:  2-------1
551:  |\     C \
552:  1 ---\    \
553:  |  A  ----\\
554:  0-----0-----1
555:  RT_TRIANGLE_SPLIT_01
556:  2
557:  |\
558:  |\\
559:  || \
560:  | \ \
561:  |  | \
562:  |  |  \
563:  |  |   \
564:  2   \ C 1
565:  |  A | / \
566:  |    | |B \
567:  |     \/   \
568:  0-----0-----1
569:  RT_TRIANGLE_SPLIT_02
570:  2
571:  |\
572:  |\\
573:  || \
574:  | \ \
575:  |  | \
576:  |  |  \
577:  |  |   \
578:  2 C \   1
579:  |\   |   \
580:  | \__|  A \
581:  | B  \\    \
582:  0-----0-----1
583: */
584: static PetscErrorCode SBRGetTriangleSplitDouble(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
585: {
586:   PetscInt              e0, e1;
587:   const PetscInt       *arr     = DMPolytopeTypeGetArrangment(DM_POLYTOPE_TRIANGLE, o);
588:   static DMPolytopeType triT2[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
589:   static PetscInt       triS2[] = {2, 3};
590:   static PetscInt       triC2[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
591:                                    DM_POLYTOPE_POINT, 1, 1,    0, DM_POLYTOPE_POINT, 1, 2, 0,
592:                                    DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0,    0,
593:                                    DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0,    1,
594:                                    DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 0,    0, DM_POLYTOPE_SEGMENT, 0,    1};
595:   static PetscInt       triO2[] = {0, 0,
596:                                    0, 0,
597:                                    0,  0, -1,
598:                                    0,  0, -1,
599:                                    0,  0,  0};

602:   /* To get the other divisions, we reorient the triangle */
603:   triC2[2]  = arr[0*2]; triC2[3] = arr[0*2+1] ? 1 : 0;
604:   triC2[7]  = arr[1*2];
605:   triC2[11] = arr[1*2];
606:   triC2[15] = arr[2*2];
607:   /* Swap the first two edges if the triangle is reversed */
608:   e0 = o < 0 ? 23: 19;
609:   e1 = o < 0 ? 19: 23;
610:   triC2[e0] = arr[0*2]; triC2[e0+1] = 0;
611:   triC2[e1] = arr[1*2]; triC2[e1+1] = o < 0 ? 1 : 0;
612:   triO2[6]  = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2*2+1]);
613:   /* Swap the first two edges if the triangle is reversed */
614:   e0 = o < 0 ? 34: 30;
615:   e1 = o < 0 ? 30: 34;
616:   triC2[e0] = arr[1*2]; triC2[e0+1] = o < 0 ? 0 : 1;
617:   triC2[e1] = arr[2*2]; triC2[e1+1] = o < 0 ? 1 : 0;
618:   triO2[9]  = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2*2+1]);
619:   /* Swap the last two edges if the triangle is reversed */
620:   triC2[41] = arr[2*2]; triC2[42] = o < 0 ? 0 : 1;
621:   triC2[45] = o < 0 ? 1 : 0;
622:   triC2[48] = o < 0 ? 0 : 1;
623:   triO2[11] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[1*2+1]);
624:   triO2[12] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[2*2+1]);
625:   *Nt = 2; *target = triT2; *size = triS2; *cone = triC2; *ornt = triO2;
626:   return 0;
627: }

629: static PetscErrorCode DMPlexTransformCellTransform_SBR(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
630: {
631:   DMLabel        trType = tr->trType;
632:   PetscInt       val;

636:   DMLabelGetValue(trType, p, &val);
637:   if (rt) *rt = val;
638:   switch (source) {
639:     case DM_POLYTOPE_POINT:
640:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
641:     case DM_POLYTOPE_QUADRILATERAL:
642:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
643:     case DM_POLYTOPE_TETRAHEDRON:
644:     case DM_POLYTOPE_HEXAHEDRON:
645:     case DM_POLYTOPE_TRI_PRISM:
646:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
647:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
648:     case DM_POLYTOPE_PYRAMID:
649:       DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
650:       break;
651:     case DM_POLYTOPE_SEGMENT:
652:       if (val == RT_EDGE) DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
653:       else                DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt);
654:       break;
655:     case DM_POLYTOPE_TRIANGLE:
656:       switch (val) {
657:         case RT_TRIANGLE_SPLIT_0: SBRGetTriangleSplitSingle(2, Nt, target, size, cone, ornt);break;
658:         case RT_TRIANGLE_SPLIT_1: SBRGetTriangleSplitSingle(0, Nt, target, size, cone, ornt);break;
659:         case RT_TRIANGLE_SPLIT_2: SBRGetTriangleSplitSingle(1, Nt, target, size, cone, ornt);break;
660:         case RT_TRIANGLE_SPLIT_21: SBRGetTriangleSplitDouble(-3, Nt, target, size, cone, ornt);break;
661:         case RT_TRIANGLE_SPLIT_10: SBRGetTriangleSplitDouble(-2, Nt, target, size, cone, ornt);break;
662:         case RT_TRIANGLE_SPLIT_02: SBRGetTriangleSplitDouble(-1, Nt, target, size, cone, ornt);break;
663:         case RT_TRIANGLE_SPLIT_12: SBRGetTriangleSplitDouble( 0, Nt, target, size, cone, ornt);break;
664:         case RT_TRIANGLE_SPLIT_20: SBRGetTriangleSplitDouble( 1, Nt, target, size, cone, ornt);break;
665:         case RT_TRIANGLE_SPLIT_01: SBRGetTriangleSplitDouble( 2, Nt, target, size, cone, ornt);break;
666:         case RT_TRIANGLE_SPLIT: DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt); break;
667:         default: DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
668:       }
669:       break;
670:     default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
671:   }
672:   return 0;
673: }

675: static PetscErrorCode DMPlexTransformSetFromOptions_SBR(PetscOptionItems *PetscOptionsObject, DMPlexTransform tr)
676: {
677:   PetscInt       cells[256], n = 256, i;
678:   PetscBool      flg;

681:   PetscOptionsHead(PetscOptionsObject,"DMPlex Options");
682:   PetscOptionsIntArray("-dm_plex_transform_sbr_ref_cell", "Mark cells for refinement", "", cells, &n, &flg);
683:   if (flg) {
684:     DMLabel active;

686:     DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active);
687:     for (i = 0; i < n; ++i) DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE);
688:     DMPlexTransformSetActive(tr, active);
689:     DMLabelDestroy(&active);
690:   }
691:   PetscOptionsTail();
692:   return 0;
693: }

695: static PetscErrorCode DMPlexTransformView_SBR(DMPlexTransform tr, PetscViewer viewer)
696: {
697:   PetscBool      isascii;

701:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
702:   if (isascii) {
703:     PetscViewerFormat format;
704:     const char       *name;

706:     PetscObjectGetName((PetscObject) tr, &name);
707:     PetscViewerASCIIPrintf(viewer, "SBR refinement %s\n", name ? name : "");
708:     PetscViewerGetFormat(viewer, &format);
709:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
710:       DMLabelView(tr->trType, viewer);
711:     }
712:   } else {
713:     SETERRQ(PetscObjectComm((PetscObject) tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject) viewer)->type_name);
714:   }
715:   return 0;
716: }

718: static PetscErrorCode DMPlexTransformDestroy_SBR(DMPlexTransform tr)
719: {
720:   DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;

722:   PetscFree(sbr->edgeLen);
723:   PetscSectionDestroy(&sbr->secEdgeLen);
724:   DMLabelDestroy(&sbr->splitPoints);
725:   PetscFree(tr->data);
726:   return 0;
727: }

729: static PetscErrorCode DMPlexTransformInitialize_SBR(DMPlexTransform tr)
730: {
731:   tr->ops->view           = DMPlexTransformView_SBR;
732:   tr->ops->setfromoptions = DMPlexTransformSetFromOptions_SBR;
733:   tr->ops->setup          = DMPlexTransformSetUp_SBR;
734:   tr->ops->destroy        = DMPlexTransformDestroy_SBR;
735:   tr->ops->celltransform  = DMPlexTransformCellTransform_SBR;
736:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_SBR;
737:   tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
738:   return 0;
739: }

741: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform tr)
742: {
743:   DMPlexRefine_SBR *f;

746:   PetscNewLog(tr, &f);
747:   tr->data = f;

749:   DMPlexTransformInitialize_SBR(tr);
750:   PetscCitationsRegister(SBRCitation, &SBRcite);
751:   return 0;
752: }