Actual source code: plextransform.c

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

  3: #include <petsc/private/petscfeimpl.h>

  5: PetscClassId DMPLEXTRANSFORM_CLASSID;

  7: PetscFunctionList DMPlexTransformList              = NULL;
  8: PetscBool         DMPlexTransformRegisterAllCalled = PETSC_FALSE;

 10: /* Construct cell type order since we must loop over cell types in depth order */
 11: static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
 12: {
 13:   PetscInt *ctO, *ctOInv;
 14:   PetscInt  c, d, off = 0;

 16:   PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv);
 17:   for (d = 3; d >= dim; --d) {
 18:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 19:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d) continue;
 20:       ctO[off++] = c;
 21:     }
 22:   }
 23:   if (dim != 0) {
 24:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 25:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) != 0) continue;
 26:       ctO[off++] = c;
 27:     }
 28:   }
 29:   for (d = dim - 1; d > 0; --d) {
 30:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 31:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d) continue;
 32:       ctO[off++] = c;
 33:     }
 34:   }
 35:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 36:     if (DMPolytopeTypeGetDim((DMPolytopeType)c) >= 0) continue;
 37:     ctO[off++] = c;
 38:   }
 40:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c;
 41:   *ctOrder    = ctO;
 42:   *ctOrderInv = ctOInv;
 43:   return 0;
 44: }

 46: /*@C
 47:   DMPlexTransformRegister - Adds a new transform component implementation

 49:   Not Collective

 51:   Input Parameters:
 52: + name        - The name of a new user-defined creation routine
 53: - create_func - The creation routine itself

 55:   Notes:
 56:   DMPlexTransformRegister() may be called multiple times to add several user-defined transforms

 58:   Sample usage:
 59: .vb
 60:   DMPlexTransformRegister("my_transform", MyTransformCreate);
 61: .ve

 63:   Then, your transform type can be chosen with the procedural interface via
 64: .vb
 65:   DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
 66:   DMPlexTransformSetType(DMPlexTransform, "my_transform");
 67: .ve
 68:   or at runtime via the option
 69: .vb
 70:   -dm_plex_transform_type my_transform
 71: .ve

 73:   Level: advanced

 75: .seealso: `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()`
 76: @*/
 77: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
 78: {
 79:   DMInitializePackage();
 80:   PetscFunctionListAdd(&DMPlexTransformList, name, create_func);
 81:   return 0;
 82: }

 84: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
 85: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
 86: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
 87: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
 88: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
 89: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
 90: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform);
 91: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform);

 93: /*@C
 94:   DMPlexTransformRegisterAll - Registers all of the transform components in the DM package.

 96:   Not Collective

 98:   Level: advanced

100: .seealso: `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
101: @*/
102: PetscErrorCode DMPlexTransformRegisterAll(void)
103: {
104:   if (DMPlexTransformRegisterAllCalled) return 0;
105:   DMPlexTransformRegisterAllCalled = PETSC_TRUE;

107:   DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter);
108:   DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular);
109:   DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox);
110:   DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld);
111:   DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL);
112:   DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR);
113:   DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D);
114:   DMPlexTransformRegister(DMPLEXEXTRUDE, DMPlexTransformCreate_Extrude);
115:   return 0;
116: }

118: /*@C
119:   DMPlexTransformRegisterDestroy - This function destroys the . It is called from PetscFinalize().

121:   Level: developer

123: .seealso: `PetscInitialize()`
124: @*/
125: PetscErrorCode DMPlexTransformRegisterDestroy(void)
126: {
127:   PetscFunctionListDestroy(&DMPlexTransformList);
128:   DMPlexTransformRegisterAllCalled = PETSC_FALSE;
129:   return 0;
130: }

132: /*@
133:   DMPlexTransformCreate - Creates an empty transform object. The type can then be set with DMPlexTransformSetType().

135:   Collective

137:   Input Parameter:
138: . comm - The communicator for the transform object

140:   Output Parameter:
141: . dm - The transform object

143:   Level: beginner

145: .seealso: `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
146: @*/
147: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
148: {
149:   DMPlexTransform t;

152:   *tr = NULL;
153:   DMInitializePackage();

155:   PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView);
156:   t->setupcalled = PETSC_FALSE;
157:   PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom);
158:   *tr = t;
159:   return 0;
160: }

162: /*@C
163:   DMSetType - Sets the particular implementation for a transform.

165:   Collective on tr

167:   Input Parameters:
168: + tr     - The transform
169: - method - The name of the transform type

171:   Options Database Key:
172: . -dm_plex_transform_type <type> - Sets the transform type; use -help for a list of available types

174:   Notes:
175:   See "petsc/include/petscdmplextransform.h" for available transform types

177:   Level: intermediate

179: .seealso: `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
180: @*/
181: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
182: {
183:   PetscErrorCode (*r)(DMPlexTransform);
184:   PetscBool match;

187:   PetscObjectTypeCompare((PetscObject)tr, method, &match);
188:   if (match) return 0;

190:   DMPlexTransformRegisterAll();
191:   PetscFunctionListFind(DMPlexTransformList, method, &r);

194:   PetscTryTypeMethod(tr, destroy);
195:   PetscMemzero(tr->ops, sizeof(*tr->ops));
196:   PetscObjectChangeTypeName((PetscObject)tr, method);
197:   (*r)(tr);
198:   return 0;
199: }

201: /*@C
202:   DMPlexTransformGetType - Gets the type name (as a string) from the transform.

204:   Not Collective

206:   Input Parameter:
207: . tr  - The DMPlexTransform

209:   Output Parameter:
210: . type - The DMPlexTransform type name

212:   Level: intermediate

214: .seealso: `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
215: @*/
216: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
217: {
220:   DMPlexTransformRegisterAll();
221:   *type = ((PetscObject)tr)->type_name;
222:   return 0;
223: }

225: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
226: {
227:   PetscViewerFormat format;

229:   PetscViewerGetFormat(v, &format);
230:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
231:     const PetscInt *trTypes = NULL;
232:     IS              trIS;
233:     PetscInt        cols = 8;
234:     PetscInt        Nrt  = 8, f, g;

236:     PetscViewerASCIIPrintf(v, "Source Starts\n");
237:     for (g = 0; g <= cols; ++g) PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]);
238:     PetscViewerASCIIPrintf(v, "\n");
239:     for (f = 0; f <= cols; ++f) PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]);
240:     PetscViewerASCIIPrintf(v, "\n");
241:     PetscViewerASCIIPrintf(v, "Target Starts\n");
242:     for (g = 0; g <= cols; ++g) PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]);
243:     PetscViewerASCIIPrintf(v, "\n");
244:     for (f = 0; f <= cols; ++f) PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]);
245:     PetscViewerASCIIPrintf(v, "\n");

247:     if (tr->trType) {
248:       DMLabelGetNumValues(tr->trType, &Nrt);
249:       DMLabelGetValueIS(tr->trType, &trIS);
250:       ISGetIndices(trIS, &trTypes);
251:     }
252:     PetscViewerASCIIPrintf(v, "Offsets\n");
253:     PetscViewerASCIIPrintf(v, "     ");
254:     for (g = 0; g < cols; ++g) PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]);
255:     PetscViewerASCIIPrintf(v, "\n");
256:     for (f = 0; f < Nrt; ++f) {
257:       PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT "  |", trTypes ? trTypes[f] : f);
258:       for (g = 0; g < cols; ++g) PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]);
259:       PetscViewerASCIIPrintf(v, " |\n");
260:     }
261:     if (trTypes) {
262:       ISGetIndices(trIS, &trTypes);
263:       ISDestroy(&trIS);
264:     }
265:   }
266:   return 0;
267: }

269: /*@C
270:   DMPlexTransformView - Views a DMPlexTransform

272:   Collective on tr

274:   Input Parameters:
275: + tr - the DMPlexTransform object to view
276: - v  - the viewer

278:   Level: beginner

280: .seealso `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
281: @*/
282: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
283: {
284:   PetscBool isascii;

287:   if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v);
290:   PetscViewerCheckWritable(v);
291:   PetscObjectPrintClassNamePrefixType((PetscObject)tr, v);
292:   PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii);
293:   if (isascii) DMPlexTransformView_Ascii(tr, v);
294:   PetscTryTypeMethod(tr, view, v);
295:   return 0;
296: }

298: /*@
299:   DMPlexTransformSetFromOptions - Sets parameters in a transform from the options database

301:   Collective on tr

303:   Input Parameter:
304: . tr - the DMPlexTransform object to set options for

306:   Options Database:
307: . -dm_plex_transform_type - Set the transform type, e.g. refine_regular

309:   Level: intermediate

311: .seealso `DMPlexTransformView()`, `DMPlexTransformCreate()`
312: @*/
313: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
314: {
315:   char        typeName[1024];
316:   const char *defName = DMPLEXREFINEREGULAR;
317:   PetscBool   flg;

320:   PetscObjectOptionsBegin((PetscObject)tr);
321:   PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg);
322:   if (flg) DMPlexTransformSetType(tr, typeName);
323:   else if (!((PetscObject)tr)->type_name) DMPlexTransformSetType(tr, defName);
324:   PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
325:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
326:   PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject);
327:   PetscOptionsEnd();
328:   return 0;
329: }

331: /*@C
332:   DMPlexTransformDestroy - Destroys a transform.

334:   Collective on tr

336:   Input Parameter:
337: . tr - the transform object to destroy

339:   Level: beginner

341: .seealso `DMPlexTransformView()`, `DMPlexTransformCreate()`
342: @*/
343: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
344: {
345:   PetscInt c;

347:   if (!*tr) return 0;
349:   if (--((PetscObject)(*tr))->refct > 0) {
350:     *tr = NULL;
351:     return 0;
352:   }

354:   if ((*tr)->ops->destroy) (*(*tr)->ops->destroy)(*tr);
355:   DMDestroy(&(*tr)->dm);
356:   DMLabelDestroy(&(*tr)->active);
357:   DMLabelDestroy(&(*tr)->trType);
358:   PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld);
359:   PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew);
360:   PetscFree2((*tr)->ctStart, (*tr)->ctStartNew);
361:   PetscFree((*tr)->offset);
362:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
363:     PetscFEDestroy(&(*tr)->coordFE[c]);
364:     PetscFEGeomDestroy(&(*tr)->refGeom[c]);
365:   }
366:   if ((*tr)->trVerts) {
367:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
368:       DMPolytopeType *rct;
369:       PetscInt       *rsize, *rcone, *rornt, Nct, n, r;

371:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0) {
372:         DMPlexTransformCellTransform((*tr), (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
373:         for (n = 0; n < Nct; ++n) {
374:           if (rct[n] == DM_POLYTOPE_POINT) continue;
375:           for (r = 0; r < rsize[n]; ++r) PetscFree((*tr)->trSubVerts[c][rct[n]][r]);
376:           PetscFree((*tr)->trSubVerts[c][rct[n]]);
377:         }
378:       }
379:       PetscFree((*tr)->trSubVerts[c]);
380:       PetscFree((*tr)->trVerts[c]);
381:     }
382:   }
383:   PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts);
384:   PetscFree2((*tr)->coordFE, (*tr)->refGeom);
385:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
386:   PetscHeaderDestroy(tr);
387:   return 0;
388: }

390: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
391: {
392:   DMLabel  trType = tr->trType;
393:   PetscInt c, cN, *off;

395:   if (trType) {
396:     DM              dm;
397:     IS              rtIS;
398:     const PetscInt *reftypes;
399:     PetscInt        Nrt, r;

401:     DMPlexTransformGetDM(tr, &dm);
402:     DMLabelGetNumValues(trType, &Nrt);
403:     DMLabelGetValueIS(trType, &rtIS);
404:     ISGetIndices(rtIS, &reftypes);
405:     PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off);
406:     for (r = 0; r < Nrt; ++r) {
407:       const PetscInt  rt = reftypes[r];
408:       IS              rtIS;
409:       const PetscInt *points;
410:       DMPolytopeType  ct;
411:       PetscInt        p;

413:       DMLabelGetStratumIS(trType, rt, &rtIS);
414:       ISGetIndices(rtIS, &points);
415:       p = points[0];
416:       ISRestoreIndices(rtIS, &points);
417:       ISDestroy(&rtIS);
418:       DMPlexGetCellType(dm, p, &ct);
419:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
420:         const DMPolytopeType ctNew = (DMPolytopeType)cN;
421:         DMPolytopeType      *rct;
422:         PetscInt            *rsize, *cone, *ornt;
423:         PetscInt             Nct, n, s;

425:         if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
426:           off[r * DM_NUM_POLYTOPES + ctNew] = -1;
427:           break;
428:         }
429:         off[r * DM_NUM_POLYTOPES + ctNew] = 0;
430:         for (s = 0; s <= r; ++s) {
431:           const PetscInt st = reftypes[s];
432:           DMPolytopeType sct;
433:           PetscInt       q, qrt;

435:           DMLabelGetStratumIS(trType, st, &rtIS);
436:           ISGetIndices(rtIS, &points);
437:           q = points[0];
438:           ISRestoreIndices(rtIS, &points);
439:           ISDestroy(&rtIS);
440:           DMPlexGetCellType(dm, q, &sct);
441:           DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt);
443:           if (st == rt) {
444:             for (n = 0; n < Nct; ++n)
445:               if (rct[n] == ctNew) break;
446:             if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
447:             break;
448:           }
449:           for (n = 0; n < Nct; ++n) {
450:             if (rct[n] == ctNew) {
451:               PetscInt sn;

453:               DMLabelGetStratumSize(trType, st, &sn);
454:               off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
455:             }
456:           }
457:         }
458:       }
459:     }
460:     ISRestoreIndices(rtIS, &reftypes);
461:     ISDestroy(&rtIS);
462:   } else {
463:     PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off);
464:     for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
465:       const DMPolytopeType ct = (DMPolytopeType)c;
466:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
467:         const DMPolytopeType ctNew = (DMPolytopeType)cN;
468:         DMPolytopeType      *rct;
469:         PetscInt            *rsize, *cone, *ornt;
470:         PetscInt             Nct, n, i;

472:         if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
473:           off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
474:           continue;
475:         }
476:         off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
477:         for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
478:           const DMPolytopeType ict  = (DMPolytopeType)ctOrderOld[i];
479:           const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];

481:           DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt);
482:           if (ict == ct) {
483:             for (n = 0; n < Nct; ++n)
484:               if (rct[n] == ctNew) break;
485:             if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
486:             break;
487:           }
488:           for (n = 0; n < Nct; ++n)
489:             if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
490:         }
491:       }
492:     }
493:   }
494:   *offset = off;
495:   return 0;
496: }

498: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
499: {
500:   DM             dm;
501:   DMPolytopeType ctCell;
502:   PetscInt       pStart, pEnd, p, c, celldim = 0;

505:   if (tr->setupcalled) return 0;
506:   PetscTryTypeMethod(tr, setup);
507:   DMPlexTransformGetDM(tr, &dm);
508:   DMPlexGetChart(dm, &pStart, &pEnd);
509:   if (pEnd > pStart) {
510:     DMPlexGetCellType(dm, 0, &ctCell);
511:   } else {
512:     PetscInt dim;

514:     DMGetDimension(dm, &dim);
515:     switch (dim) {
516:     case 0:
517:       ctCell = DM_POLYTOPE_POINT;
518:       break;
519:     case 1:
520:       ctCell = DM_POLYTOPE_SEGMENT;
521:       break;
522:     case 2:
523:       ctCell = DM_POLYTOPE_TRIANGLE;
524:       break;
525:     case 3:
526:       ctCell = DM_POLYTOPE_TETRAHEDRON;
527:       break;
528:     default:
529:       ctCell = DM_POLYTOPE_UNKNOWN;
530:     }
531:   }
532:   DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld);
533:   for (p = pStart; p < pEnd; ++p) {
534:     DMPolytopeType  ct;
535:     DMPolytopeType *rct;
536:     PetscInt       *rsize, *cone, *ornt;
537:     PetscInt        Nct, n;

539:     DMPlexGetCellType(dm, p, &ct);
541:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt);
542:     for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
543:   }
544:   DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew);
545:   /* Construct sizes and offsets for each cell type */
546:   if (!tr->ctStart) {
547:     PetscInt *ctS, *ctSN, *ctC, *ctCN;

549:     PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN);
550:     PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN);
551:     for (p = pStart; p < pEnd; ++p) {
552:       DMPolytopeType  ct;
553:       DMPolytopeType *rct;
554:       PetscInt       *rsize, *cone, *ornt;
555:       PetscInt        Nct, n;

557:       DMPlexGetCellType(dm, p, &ct);
559:       ++ctC[ct];
560:       DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt);
561:       for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
562:     }
563:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
564:       const PetscInt cto  = tr->ctOrderOld[c];
565:       const PetscInt cton = tr->ctOrderOld[c + 1];
566:       const PetscInt ctn  = tr->ctOrderNew[c];
567:       const PetscInt ctnn = tr->ctOrderNew[c + 1];

569:       ctS[cton]  = ctS[cto] + ctC[cto];
570:       ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
571:     }
572:     PetscFree2(ctC, ctCN);
573:     tr->ctStart    = ctS;
574:     tr->ctStartNew = ctSN;
575:   }
576:   DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset);
577:   tr->setupcalled = PETSC_TRUE;
578:   return 0;
579: }

581: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
582: {
585:   *dm = tr->dm;
586:   return 0;
587: }

589: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
590: {
593:   PetscObjectReference((PetscObject)dm);
594:   DMDestroy(&tr->dm);
595:   tr->dm = dm;
596:   return 0;
597: }

599: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
600: {
603:   *active = tr->active;
604:   return 0;
605: }

607: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
608: {
611:   PetscObjectReference((PetscObject)active);
612:   DMLabelDestroy(&tr->active);
613:   tr->active = active;
614:   return 0;
615: }

617: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
618: {
619:   if (!tr->coordFE[ct]) {
620:     PetscInt dim, cdim;

622:     dim = DMPolytopeTypeGetDim(ct);
623:     DMGetCoordinateDim(tr->dm, &cdim);
624:     PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]);
625:     {
626:       PetscDualSpace  dsp;
627:       PetscQuadrature quad;
628:       DM              K;
629:       PetscFEGeom    *cg;
630:       PetscScalar    *Xq;
631:       PetscReal      *xq, *wq;
632:       PetscInt        Nq, q;

634:       DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq);
635:       PetscMalloc1(Nq * cdim, &xq);
636:       for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
637:       PetscMalloc1(Nq, &wq);
638:       for (q = 0; q < Nq; ++q) wq[q] = 1.0;
639:       PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
640:       PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);
641:       PetscFESetQuadrature(tr->coordFE[ct], quad);

643:       PetscFEGetDualSpace(tr->coordFE[ct], &dsp);
644:       PetscDualSpaceGetDM(dsp, &K);
645:       PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]);
646:       cg = tr->refGeom[ct];
647:       DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
648:       PetscQuadratureDestroy(&quad);
649:     }
650:   }
651:   *fe = tr->coordFE[ct];
652:   return 0;
653: }

655: /*@
656:   DMPlexTransformSetDimensions - Set the dimensions for the transformed DM

658:   Input Parameters:
659: + tr - The DMPlexTransform object
660: - dm - The original DM

662:   Output Parameter:
663: . tdm - The transformed DM

665:   Level: advanced

667: .seealso: `DMPlexTransformApply()`, `DMPlexTransformCreate()`
668: @*/
669: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
670: {
671:   if (tr->ops->setdimensions) PetscUseTypeMethod(tr, setdimensions, dm, tdm);
672:   else {
673:     PetscInt dim, cdim;

675:     DMGetDimension(dm, &dim);
676:     DMSetDimension(tdm, dim);
677:     DMGetCoordinateDim(dm, &cdim);
678:     DMSetCoordinateDim(tdm, cdim);
679:   }
680:   return 0;
681: }

683: /*@
684:   DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.

686:   Not collective

688:   Input Parameters:
689: + tr    - The DMPlexTransform
690: . ct    - The type of the original point which produces the new point
691: . ctNew - The type of the new point
692: . p     - The original point which produces the new point
693: - r     - The replica number of the new point, meaning it is the rth point of type ctNew produced from p

695:   Output Parameters:
696: . pNew  - The new point number

698:   Level: developer

700: .seealso: `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
701: @*/
702: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
703: {
704:   DMPolytopeType *rct;
705:   PetscInt       *rsize, *cone, *ornt;
706:   PetscInt        rt, Nct, n, off, rp;
707:   DMLabel         trType = tr->trType;
708:   PetscInt        ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
709:   PetscInt        ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
710:   PetscInt        newp = ctSN, cind;

714:   DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt);
715:   if (trType) {
716:     DMLabelGetValueIndex(trType, rt, &cind);
717:     DMLabelGetStratumPointIndex(trType, rt, p, &rp);
719:   } else {
720:     cind = ct;
721:     rp   = p - ctS;
722:   }
723:   off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
725:   newp += off;
726:   for (n = 0; n < Nct; ++n) {
727:     if (rct[n] == ctNew) {
728:       if (rsize[n] && r >= rsize[n])
729:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ") for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
730:       newp += rp * rsize[n] + r;
731:       break;
732:     }
733:   }

736:   *pNew = newp;
737:   return 0;
738: }

740: /*@
741:   DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.

743:   Not collective

745:   Input Parameters:
746: + tr    - The DMPlexTransform
747: - pNew  - The new point number

749:   Output Parameters:
750: + ct    - The type of the original point which produces the new point
751: . ctNew - The type of the new point
752: . p     - The original point which produces the new point
753: - r     - The replica number of the new point, meaning it is the rth point of type ctNew produced from p

755:   Level: developer

757: .seealso: `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
758: @*/
759: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
760: {
761:   DMLabel         trType = tr->trType;
762:   DMPolytopeType *rct;
763:   PetscInt       *rsize, *cone, *ornt;
764:   PetscInt        rt, Nct, n, rp = 0, rO = 0, pO;
765:   PetscInt        offset = -1, ctS, ctE, ctO = 0, ctN, ctTmp;

767:   for (ctN = 0; ctN < DM_NUM_POLYTOPES; ++ctN) {
768:     PetscInt ctSN = tr->ctStartNew[ctN], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctN] + 1]];

770:     if ((pNew >= ctSN) && (pNew < ctEN)) break;
771:   }
773:   if (trType) {
774:     DM              dm;
775:     IS              rtIS;
776:     const PetscInt *reftypes;
777:     PetscInt        Nrt, r, rtStart;

779:     DMPlexTransformGetDM(tr, &dm);
780:     DMLabelGetNumValues(trType, &Nrt);
781:     DMLabelGetValueIS(trType, &rtIS);
782:     ISGetIndices(rtIS, &reftypes);
783:     for (r = 0; r < Nrt; ++r) {
784:       const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];

786:       if (tr->ctStartNew[ctN] + off > pNew) continue;
787:       /* Check that any of this refinement type exist */
788:       /* TODO Actually keep track of the number produced here instead */
789:       if (off > offset) {
790:         rt     = reftypes[r];
791:         offset = off;
792:       }
793:     }
794:     ISRestoreIndices(rtIS, &reftypes);
795:     ISDestroy(&rtIS);
797:     /* TODO Map refinement types to cell types */
798:     DMLabelGetStratumBounds(trType, rt, &rtStart, NULL);
800:     for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
801:       PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];

803:       if ((rtStart >= ctS) && (rtStart < ctE)) break;
804:     }
806:   } else {
807:     for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
808:       const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];

810:       if (tr->ctStartNew[ctN] + off > pNew) continue;
811:       if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
812:       /* TODO Actually keep track of the number produced here instead */
813:       if (off > offset) {
814:         ctO    = ctTmp;
815:         offset = off;
816:       }
817:     }
819:   }
820:   ctS = tr->ctStart[ctO];
821:   ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
822:   DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, ctS, &rt, &Nct, &rct, &rsize, &cone, &ornt);
823:   for (n = 0; n < Nct; ++n) {
824:     if ((PetscInt)rct[n] == ctN) {
825:       PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset;

827:       rp = tmp / rsize[n];
828:       rO = tmp % rsize[n];
829:       break;
830:     }
831:   }
833:   pO = rp + ctS;
835:   if (ct) *ct = (DMPolytopeType)ctO;
836:   if (ctNew) *ctNew = (DMPolytopeType)ctN;
837:   if (p) *p = pO;
838:   if (r) *r = rO;
839:   return 0;
840: }

842: /*@
843:   DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.

845:   Input Parameters:
846: + tr     - The DMPlexTransform object
847: . source - The source cell type
848: - p      - The source point, which can also determine the refine type

850:   Output Parameters:
851: + rt     - The refine type for this point
852: . Nt     - The number of types produced by this point
853: . target - An array of length Nt giving the types produced
854: . size   - An array of length Nt giving the number of cells of each type produced
855: . cone   - An array of length Nt*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
856: - ornt   - An array of length Nt*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point

858:   Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
859:   need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
860:   division (described in these outputs) of the cell in the original mesh. The point identifier is given by
861: $   the number of cones to be taken, or 0 for the current cell
862: $   the cell cone point number at each level from which it is subdivided
863: $   the replica number r of the subdivision.
864: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
865: $   Nt     = 2
866: $   target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
867: $   size   = {1, 2}
868: $   cone   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,  DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
869: $   ornt   = {                         0,                       0,                        0,                          0}

871:   Level: advanced

873: .seealso: `DMPlexTransformApply()`, `DMPlexTransformCreate()`
874: @*/
875: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
876: {
877:   PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
878:   return 0;
879: }

881: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
882: {
883:   *rnew = r;
884:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
885:   return 0;
886: }

888: /* Returns the same thing */
889: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
890: {
891:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
892:   static PetscInt       vertexS[] = {1};
893:   static PetscInt       vertexC[] = {0};
894:   static PetscInt       vertexO[] = {0};
895:   static DMPolytopeType edgeT[]   = {DM_POLYTOPE_SEGMENT};
896:   static PetscInt       edgeS[]   = {1};
897:   static PetscInt       edgeC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
898:   static PetscInt       edgeO[]   = {0, 0};
899:   static DMPolytopeType tedgeT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR};
900:   static PetscInt       tedgeS[]  = {1};
901:   static PetscInt       tedgeC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
902:   static PetscInt       tedgeO[]  = {0, 0};
903:   static DMPolytopeType triT[]    = {DM_POLYTOPE_TRIANGLE};
904:   static PetscInt       triS[]    = {1};
905:   static PetscInt       triC[]    = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
906:   static PetscInt       triO[]    = {0, 0, 0};
907:   static DMPolytopeType quadT[]   = {DM_POLYTOPE_QUADRILATERAL};
908:   static PetscInt       quadS[]   = {1};
909:   static PetscInt       quadC[]   = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0};
910:   static PetscInt       quadO[]   = {0, 0, 0, 0};
911:   static DMPolytopeType tquadT[]  = {DM_POLYTOPE_SEG_PRISM_TENSOR};
912:   static PetscInt       tquadS[]  = {1};
913:   static PetscInt       tquadC[]  = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0};
914:   static PetscInt       tquadO[]  = {0, 0, 0, 0};
915:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_TETRAHEDRON};
916:   static PetscInt       tetS[]    = {1};
917:   static PetscInt       tetC[]    = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0};
918:   static PetscInt       tetO[]    = {0, 0, 0, 0};
919:   static DMPolytopeType hexT[]    = {DM_POLYTOPE_HEXAHEDRON};
920:   static PetscInt       hexS[]    = {1};
921:   static PetscInt       hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0};
922:   static PetscInt       hexO[] = {0, 0, 0, 0, 0, 0};
923:   static DMPolytopeType tripT[]   = {DM_POLYTOPE_TRI_PRISM};
924:   static PetscInt       tripS[]   = {1};
925:   static PetscInt       tripC[]   = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
926:   static PetscInt       tripO[]   = {0, 0, 0, 0, 0};
927:   static DMPolytopeType ttripT[]  = {DM_POLYTOPE_TRI_PRISM_TENSOR};
928:   static PetscInt       ttripS[]  = {1};
929:   static PetscInt       ttripC[]  = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
930:   static PetscInt       ttripO[]  = {0, 0, 0, 0, 0};
931:   static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
932:   static PetscInt       tquadpS[] = {1};
933:   static PetscInt       tquadpC[] = {DM_POLYTOPE_QUADRILATERAL,    1, 0, 0, DM_POLYTOPE_QUADRILATERAL,    1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
934:                                      DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
935:   static PetscInt       tquadpO[] = {0, 0, 0, 0, 0, 0};
936:   static DMPolytopeType pyrT[]    = {DM_POLYTOPE_PYRAMID};
937:   static PetscInt       pyrS[]    = {1};
938:   static PetscInt       pyrC[]    = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0};
939:   static PetscInt       pyrO[]    = {0, 0, 0, 0, 0};

941:   if (rt) *rt = 0;
942:   switch (source) {
943:   case DM_POLYTOPE_POINT:
944:     *Nt     = 1;
945:     *target = vertexT;
946:     *size   = vertexS;
947:     *cone   = vertexC;
948:     *ornt   = vertexO;
949:     break;
950:   case DM_POLYTOPE_SEGMENT:
951:     *Nt     = 1;
952:     *target = edgeT;
953:     *size   = edgeS;
954:     *cone   = edgeC;
955:     *ornt   = edgeO;
956:     break;
957:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
958:     *Nt     = 1;
959:     *target = tedgeT;
960:     *size   = tedgeS;
961:     *cone   = tedgeC;
962:     *ornt   = tedgeO;
963:     break;
964:   case DM_POLYTOPE_TRIANGLE:
965:     *Nt     = 1;
966:     *target = triT;
967:     *size   = triS;
968:     *cone   = triC;
969:     *ornt   = triO;
970:     break;
971:   case DM_POLYTOPE_QUADRILATERAL:
972:     *Nt     = 1;
973:     *target = quadT;
974:     *size   = quadS;
975:     *cone   = quadC;
976:     *ornt   = quadO;
977:     break;
978:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
979:     *Nt     = 1;
980:     *target = tquadT;
981:     *size   = tquadS;
982:     *cone   = tquadC;
983:     *ornt   = tquadO;
984:     break;
985:   case DM_POLYTOPE_TETRAHEDRON:
986:     *Nt     = 1;
987:     *target = tetT;
988:     *size   = tetS;
989:     *cone   = tetC;
990:     *ornt   = tetO;
991:     break;
992:   case DM_POLYTOPE_HEXAHEDRON:
993:     *Nt     = 1;
994:     *target = hexT;
995:     *size   = hexS;
996:     *cone   = hexC;
997:     *ornt   = hexO;
998:     break;
999:   case DM_POLYTOPE_TRI_PRISM:
1000:     *Nt     = 1;
1001:     *target = tripT;
1002:     *size   = tripS;
1003:     *cone   = tripC;
1004:     *ornt   = tripO;
1005:     break;
1006:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
1007:     *Nt     = 1;
1008:     *target = ttripT;
1009:     *size   = ttripS;
1010:     *cone   = ttripC;
1011:     *ornt   = ttripO;
1012:     break;
1013:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1014:     *Nt     = 1;
1015:     *target = tquadpT;
1016:     *size   = tquadpS;
1017:     *cone   = tquadpC;
1018:     *ornt   = tquadpO;
1019:     break;
1020:   case DM_POLYTOPE_PYRAMID:
1021:     *Nt     = 1;
1022:     *target = pyrT;
1023:     *size   = pyrS;
1024:     *cone   = pyrC;
1025:     *ornt   = pyrO;
1026:     break;
1027:   default:
1028:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1029:   }
1030:   return 0;
1031: }

1033: /*@
1034:   DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point

1036:   Not collective

1038:   Input Parameters:
1039: + tr  - The DMPlexTransform
1040: . sct - The source point cell type, from whom the new cell is being produced
1041: . sp  - The source point
1042: . so  - The orientation of the source point in its enclosing parent
1043: . tct - The target point cell type
1044: . r   - The replica number requested for the produced cell type
1045: - o   - The orientation of the replica

1047:   Output Parameters:
1048: + rnew - The replica number, given the orientation of the parent
1049: - onew - The replica orientation, given the orientation of the parent

1051:   Level: advanced

1053: .seealso: `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1054: @*/
1055: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1056: {
1058:   PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1059:   return 0;
1060: }

1062: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1063: {
1064:   DM       dm;
1065:   PetscInt pStart, pEnd, p, pNew;

1067:   DMPlexTransformGetDM(tr, &dm);
1068:   /* Must create the celltype label here so that we do not automatically try to compute the types */
1069:   DMCreateLabel(rdm, "celltype");
1070:   DMPlexGetChart(dm, &pStart, &pEnd);
1071:   for (p = pStart; p < pEnd; ++p) {
1072:     DMPolytopeType  ct;
1073:     DMPolytopeType *rct;
1074:     PetscInt       *rsize, *rcone, *rornt;
1075:     PetscInt        Nct, n, r;

1077:     DMPlexGetCellType(dm, p, &ct);
1078:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1079:     for (n = 0; n < Nct; ++n) {
1080:       for (r = 0; r < rsize[n]; ++r) {
1081:         DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1082:         DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
1083:         DMPlexSetCellType(rdm, pNew, rct[n]);
1084:       }
1085:     }
1086:   }
1087:   /* Let the DM know we have set all the cell types */
1088:   {
1089:     DMLabel  ctLabel;
1090:     DM_Plex *plex = (DM_Plex *)rdm->data;

1092:     DMPlexGetCellTypeLabel(rdm, &ctLabel);
1093:     PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState);
1094:   }
1095:   return 0;
1096: }

1098: #if 0
1099: static PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1100: {
1101:   PetscInt ctNew;

1105:   /* TODO Can do bisection since everything is sorted */
1106:   for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
1107:     PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew]+1]];

1109:     if (q >= ctSN && q < ctEN) break;
1110:   }
1112:   *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType) ctNew);
1113:   return 0;
1114: }
1115: #endif

1117: /* The orientation o is for the interior of the cell p */
1118: static PetscErrorCode DMPlexTransformGetCone_Internal(DMPlexTransform tr, PetscInt p, PetscInt o, DMPolytopeType ct, DMPolytopeType ctNew, const PetscInt rcone[], PetscInt *coneoff, const PetscInt rornt[], PetscInt *orntoff, PetscInt coneNew[], PetscInt orntNew[])
1119: {
1120:   DM              dm;
1121:   const PetscInt  csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1122:   const PetscInt *cone;
1123:   PetscInt        c, coff = *coneoff, ooff = *orntoff;

1125:   DMPlexTransformGetDM(tr, &dm);
1126:   DMPlexGetCone(dm, p, &cone);
1127:   for (c = 0; c < csizeNew; ++c) {
1128:     PetscInt             ppp   = -1;                            /* Parent Parent point: Parent of point pp */
1129:     PetscInt             pp    = p;                             /* Parent point: Point in the original mesh producing new cone point */
1130:     PetscInt             po    = 0;                             /* Orientation of parent point pp in parent parent point ppp */
1131:     DMPolytopeType       pct   = ct;                            /* Parent type: Cell type for parent of new cone point */
1132:     const PetscInt      *pcone = cone;                          /* Parent cone: Cone of parent point pp */
1133:     PetscInt             pr    = -1;                            /* Replica number of pp that produces new cone point  */
1134:     const DMPolytopeType ft    = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1135:     const PetscInt       fn    = rcone[coff++];                 /* Number of cones of p that need to be taken when producing new cone point */
1136:     PetscInt             fo    = rornt[ooff++];                 /* Orientation of new cone point in pNew */
1137:     PetscInt             lc;

1139:     /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1140:     for (lc = 0; lc < fn; ++lc) {
1141:       const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po);
1142:       const PetscInt  acp  = rcone[coff++];
1143:       const PetscInt  pcp  = parr[acp * 2];
1144:       const PetscInt  pco  = parr[acp * 2 + 1];
1145:       const PetscInt *ppornt;

1147:       ppp = pp;
1148:       pp  = pcone[pcp];
1149:       DMPlexGetCellType(dm, pp, &pct);
1150:       DMPlexGetCone(dm, pp, &pcone);
1151:       DMPlexGetConeOrientation(dm, ppp, &ppornt);
1152:       po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1153:     }
1154:     pr = rcone[coff++];
1155:     /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1156:     DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo);
1157:     DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]);
1158:     orntNew[c] = fo;
1159:   }
1160:   *coneoff = coff;
1161:   *orntoff = ooff;
1162:   return 0;
1163: }

1165: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1166: {
1167:   DM             dm;
1168:   DMPolytopeType ct;
1169:   PetscInt      *coneNew, *orntNew;
1170:   PetscInt       maxConeSize = 0, pStart, pEnd, p, pNew;

1172:   DMPlexTransformGetDM(tr, &dm);
1173:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1174:   DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1175:   DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1176:   DMPlexGetChart(dm, &pStart, &pEnd);
1177:   for (p = pStart; p < pEnd; ++p) {
1178:     const PetscInt *cone, *ornt;
1179:     PetscInt        coff, ooff;
1180:     DMPolytopeType *rct;
1181:     PetscInt       *rsize, *rcone, *rornt;
1182:     PetscInt        Nct, n, r;

1184:     DMPlexGetCellType(dm, p, &ct);
1185:     DMPlexGetCone(dm, p, &cone);
1186:     DMPlexGetConeOrientation(dm, p, &ornt);
1187:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1188:     for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1189:       const DMPolytopeType ctNew = rct[n];

1191:       for (r = 0; r < rsize[n]; ++r) {
1192:         DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1193:         DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew);
1194:         DMPlexSetCone(rdm, pNew, coneNew);
1195:         DMPlexSetConeOrientation(rdm, pNew, orntNew);
1196:       }
1197:     }
1198:   }
1199:   DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1200:   DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1201:   DMViewFromOptions(rdm, NULL, "-rdm_view");
1202:   DMPlexSymmetrize(rdm);
1203:   DMPlexStratify(rdm);
1204:   return 0;
1205: }

1207: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1208: {
1209:   DM              dm;
1210:   DMPolytopeType  ct, qct;
1211:   DMPolytopeType *rct;
1212:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
1213:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

1218:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1219:   DMPlexTransformGetDM(tr, &dm);
1220:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1221:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1222:   DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1223:   DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1224:   for (n = 0; n < Nct; ++n) {
1225:     const DMPolytopeType ctNew    = rct[n];
1226:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1227:     PetscInt             Nr       = rsize[n], fn, c;

1229:     if (ctNew == qct) Nr = r;
1230:     for (nr = 0; nr < Nr; ++nr) {
1231:       for (c = 0; c < csizeNew; ++c) {
1232:         ++coff;             /* Cell type of new cone point */
1233:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1234:         coff += fn;
1235:         ++coff; /* Replica number of new cone point */
1236:         ++ooff; /* Orientation of new cone point */
1237:       }
1238:     }
1239:     if (ctNew == qct) break;
1240:   }
1241:   DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1242:   *cone = qcone;
1243:   *ornt = qornt;
1244:   return 0;
1245: }

1247: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1248: {
1249:   DM              dm;
1250:   DMPolytopeType  ct, qct;
1251:   DMPolytopeType *rct;
1252:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
1253:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

1257:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1258:   DMPlexTransformGetDM(tr, &dm);
1259:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1260:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1261:   DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1262:   DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1263:   for (n = 0; n < Nct; ++n) {
1264:     const DMPolytopeType ctNew    = rct[n];
1265:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1266:     PetscInt             Nr       = rsize[n], fn, c;

1268:     if (ctNew == qct) Nr = r;
1269:     for (nr = 0; nr < Nr; ++nr) {
1270:       for (c = 0; c < csizeNew; ++c) {
1271:         ++coff;             /* Cell type of new cone point */
1272:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1273:         coff += fn;
1274:         ++coff; /* Replica number of new cone point */
1275:         ++ooff; /* Orientation of new cone point */
1276:       }
1277:     }
1278:     if (ctNew == qct) break;
1279:   }
1280:   DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1281:   *cone = qcone;
1282:   *ornt = qornt;
1283:   return 0;
1284: }

1286: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1287: {
1288:   DM dm;

1292:   DMPlexTransformGetDM(tr, &dm);
1293:   DMRestoreWorkArray(dm, 0, MPIU_INT, cone);
1294:   DMRestoreWorkArray(dm, 0, MPIU_INT, ornt);
1295:   return 0;
1296: }

1298: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1299: {
1300:   PetscInt ict;

1302:   PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts);
1303:   for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1304:     const DMPolytopeType ct = (DMPolytopeType)ict;
1305:     DMPlexTransform      reftr;
1306:     DM                   refdm, trdm;
1307:     Vec                  coordinates;
1308:     const PetscScalar   *coords;
1309:     DMPolytopeType      *rct;
1310:     PetscInt            *rsize, *rcone, *rornt;
1311:     PetscInt             Nct, n, r, pNew;
1312:     PetscInt             trdim, vStart, vEnd, Nc;
1313:     const PetscInt       debug = 0;
1314:     const char          *typeName;

1316:     /* Since points are 0-dimensional, coordinates make no sense */
1317:     if (DMPolytopeTypeGetDim(ct) <= 0) continue;
1318:     DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm);
1319:     DMPlexTransformCreate(PETSC_COMM_SELF, &reftr);
1320:     DMPlexTransformSetDM(reftr, refdm);
1321:     DMPlexTransformGetType(tr, &typeName);
1322:     DMPlexTransformSetType(reftr, typeName);
1323:     DMPlexTransformSetUp(reftr);
1324:     DMPlexTransformApply(reftr, refdm, &trdm);

1326:     DMGetDimension(trdm, &trdim);
1327:     DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd);
1328:     tr->trNv[ct] = vEnd - vStart;
1329:     DMGetCoordinatesLocal(trdm, &coordinates);
1330:     VecGetLocalSize(coordinates, &Nc);
1332:     PetscCalloc1(Nc, &tr->trVerts[ct]);
1333:     VecGetArrayRead(coordinates, &coords);
1334:     PetscArraycpy(tr->trVerts[ct], coords, Nc);
1335:     VecRestoreArrayRead(coordinates, &coords);

1337:     PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]);
1338:     DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1339:     for (n = 0; n < Nct; ++n) {
1340:       /* Since points are 0-dimensional, coordinates make no sense */
1341:       if (rct[n] == DM_POLYTOPE_POINT) continue;
1342:       PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]);
1343:       for (r = 0; r < rsize[n]; ++r) {
1344:         PetscInt *closure = NULL;
1345:         PetscInt  clSize, cl, Nv = 0;

1347:         PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]);
1348:         DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew);
1349:         DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1350:         for (cl = 0; cl < clSize * 2; cl += 2) {
1351:           const PetscInt sv = closure[cl];

1353:           if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1354:         }
1355:         DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1357:       }
1358:     }
1359:     if (debug) {
1360:       DMPolytopeType *rct;
1361:       PetscInt       *rsize, *rcone, *rornt;
1362:       PetscInt        v, dE = trdim, d, off = 0;

1364:       PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]);
1365:       for (v = 0; v < tr->trNv[ct]; ++v) {
1366:         PetscPrintf(PETSC_COMM_SELF, "  ");
1367:         for (d = 0; d < dE; ++d) PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]));
1368:         PetscPrintf(PETSC_COMM_SELF, "\n");
1369:       }

1371:       DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1372:       for (n = 0; n < Nct; ++n) {
1373:         if (rct[n] == DM_POLYTOPE_POINT) continue;
1374:         PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]);
1375:         for (r = 0; r < rsize[n]; ++r) {
1376:           PetscPrintf(PETSC_COMM_SELF, "  ");
1377:           for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]);
1378:           PetscPrintf(PETSC_COMM_SELF, "\n");
1379:         }
1380:       }
1381:     }
1382:     DMDestroy(&refdm);
1383:     DMDestroy(&trdm);
1384:     DMPlexTransformDestroy(&reftr);
1385:   }
1386:   return 0;
1387: }

1389: /*
1390:   DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type

1392:   Input Parameters:
1393: + tr - The DMPlexTransform object
1394: - ct - The cell type

1396:   Output Parameters:
1397: + Nv      - The number of transformed vertices in the closure of the reference cell of given type
1398: - trVerts - The coordinates of these vertices in the reference cell

1400:   Level: developer

1402: .seealso: `DMPlexTransformGetSubcellVertices()`
1403: */
1404: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1405: {
1406:   if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1407:   if (Nv) *Nv = tr->trNv[ct];
1408:   if (trVerts) *trVerts = tr->trVerts[ct];
1409:   return 0;
1410: }

1412: /*
1413:   DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type

1415:   Input Parameters:
1416: + tr  - The DMPlexTransform object
1417: . ct  - The cell type
1418: . rct - The subcell type
1419: - r   - The subcell index

1421:   Output Parameter:
1422: . subVerts - The indices of these vertices in the set of vertices returned by DMPlexTransformGetCellVertices()

1424:   Level: developer

1426: .seealso: `DMPlexTransformGetCellVertices()`
1427: */
1428: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1429: {
1430:   if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1432:   if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1433:   return 0;
1434: }

1436: /* Computes new vertex as the barycenter, or centroid */
1437: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1438: {
1439:   PetscInt v, d;

1443:   for (d = 0; d < dE; ++d) out[d] = 0.0;
1444:   for (v = 0; v < Nv; ++v)
1445:     for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1446:   for (d = 0; d < dE; ++d) out[d] /= Nv;
1447:   return 0;
1448: }

1450: /*@
1451:   DMPlexTransformMapCoordinates - Calculate new coordinates for produced points

1453:   Not collective

1455:   Input Parameters:
1456: + cr   - The DMPlexCellRefiner
1457: . pct  - The cell type of the parent, from whom the new cell is being produced
1458: . ct   - The type being produced
1459: . p    - The original point
1460: . r    - The replica number requested for the produced cell type
1461: . Nv   - Number of vertices in the closure of the parent cell
1462: . dE   - Spatial dimension
1463: - in   - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell

1465:   Output Parameter:
1466: . out - The coordinates of the new vertices

1468:   Level: intermediate

1470: .seealso: `DMPlexTransformApply()`
1471: @*/
1472: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1473: {
1475:   PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1476:   return 0;
1477: }

1479: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1480: {
1481:   DM              dm;
1482:   IS              valueIS;
1483:   const PetscInt *values;
1484:   PetscInt        defVal, Nv, val;

1486:   DMPlexTransformGetDM(tr, &dm);
1487:   DMLabelGetDefaultValue(label, &defVal);
1488:   DMLabelSetDefaultValue(labelNew, defVal);
1489:   DMLabelGetValueIS(label, &valueIS);
1490:   ISGetLocalSize(valueIS, &Nv);
1491:   ISGetIndices(valueIS, &values);
1492:   for (val = 0; val < Nv; ++val) {
1493:     IS              pointIS;
1494:     const PetscInt *points;
1495:     PetscInt        numPoints, p;

1497:     /* Ensure refined label is created with same number of strata as
1498:      * original (even if no entries here). */
1499:     DMLabelAddStratum(labelNew, values[val]);
1500:     DMLabelGetStratumIS(label, values[val], &pointIS);
1501:     ISGetLocalSize(pointIS, &numPoints);
1502:     ISGetIndices(pointIS, &points);
1503:     for (p = 0; p < numPoints; ++p) {
1504:       const PetscInt  point = points[p];
1505:       DMPolytopeType  ct;
1506:       DMPolytopeType *rct;
1507:       PetscInt       *rsize, *rcone, *rornt;
1508:       PetscInt        Nct, n, r, pNew = 0;

1510:       DMPlexGetCellType(dm, point, &ct);
1511:       DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1512:       for (n = 0; n < Nct; ++n) {
1513:         for (r = 0; r < rsize[n]; ++r) {
1514:           DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew);
1515:           DMLabelSetValue(labelNew, pNew, values[val]);
1516:         }
1517:       }
1518:     }
1519:     ISRestoreIndices(pointIS, &points);
1520:     ISDestroy(&pointIS);
1521:   }
1522:   ISRestoreIndices(valueIS, &values);
1523:   ISDestroy(&valueIS);
1524:   return 0;
1525: }

1527: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1528: {
1529:   DM       dm;
1530:   PetscInt numLabels, l;

1532:   DMPlexTransformGetDM(tr, &dm);
1533:   DMGetNumLabels(dm, &numLabels);
1534:   for (l = 0; l < numLabels; ++l) {
1535:     DMLabel     label, labelNew;
1536:     const char *lname;
1537:     PetscBool   isDepth, isCellType;

1539:     DMGetLabelName(dm, l, &lname);
1540:     PetscStrcmp(lname, "depth", &isDepth);
1541:     if (isDepth) continue;
1542:     PetscStrcmp(lname, "celltype", &isCellType);
1543:     if (isCellType) continue;
1544:     DMCreateLabel(rdm, lname);
1545:     DMGetLabel(dm, lname, &label);
1546:     DMGetLabel(rdm, lname, &labelNew);
1547:     RefineLabel_Internal(tr, label, labelNew);
1548:   }
1549:   return 0;
1550: }

1552: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1553: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1554: {
1555:   DM       dm;
1556:   PetscInt Nf, f, Nds, s;

1558:   DMPlexTransformGetDM(tr, &dm);
1559:   DMGetNumFields(dm, &Nf);
1560:   for (f = 0; f < Nf; ++f) {
1561:     DMLabel     label, labelNew;
1562:     PetscObject obj;
1563:     const char *lname;

1565:     DMGetField(rdm, f, &label, &obj);
1566:     if (!label) continue;
1567:     PetscObjectGetName((PetscObject)label, &lname);
1568:     DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1569:     RefineLabel_Internal(tr, label, labelNew);
1570:     DMSetField_Internal(rdm, f, labelNew, obj);
1571:     DMLabelDestroy(&labelNew);
1572:   }
1573:   DMGetNumDS(dm, &Nds);
1574:   for (s = 0; s < Nds; ++s) {
1575:     DMLabel     label, labelNew;
1576:     const char *lname;

1578:     DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
1579:     if (!label) continue;
1580:     PetscObjectGetName((PetscObject)label, &lname);
1581:     DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1582:     RefineLabel_Internal(tr, label, labelNew);
1583:     DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
1584:     DMLabelDestroy(&labelNew);
1585:   }
1586:   return 0;
1587: }

1589: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1590: {
1591:   DM                 dm;
1592:   PetscSF            sf, sfNew;
1593:   PetscInt           numRoots, numLeaves, numLeavesNew = 0, l, m;
1594:   const PetscInt    *localPoints;
1595:   const PetscSFNode *remotePoints;
1596:   PetscInt          *localPointsNew;
1597:   PetscSFNode       *remotePointsNew;
1598:   PetscInt           pStartNew, pEndNew, pNew;
1599:   /* Brute force algorithm */
1600:   PetscSF         rsf;
1601:   PetscSection    s;
1602:   const PetscInt *rootdegree;
1603:   PetscInt       *rootPointsNew, *remoteOffsets;
1604:   PetscInt        numPointsNew, pStart, pEnd, p;

1606:   DMPlexTransformGetDM(tr, &dm);
1607:   DMPlexGetChart(rdm, &pStartNew, &pEndNew);
1608:   DMGetPointSF(dm, &sf);
1609:   DMGetPointSF(rdm, &sfNew);
1610:   /* Calculate size of new SF */
1611:   PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
1612:   if (numRoots < 0) return 0;
1613:   for (l = 0; l < numLeaves; ++l) {
1614:     const PetscInt  p = localPoints[l];
1615:     DMPolytopeType  ct;
1616:     DMPolytopeType *rct;
1617:     PetscInt       *rsize, *rcone, *rornt;
1618:     PetscInt        Nct, n;

1620:     DMPlexGetCellType(dm, p, &ct);
1621:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1622:     for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
1623:   }
1624:   /* Send new root point numbers
1625:        It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1626:   */
1627:   DMPlexGetChart(dm, &pStart, &pEnd);
1628:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s);
1629:   PetscSectionSetChart(s, pStart, pEnd);
1630:   for (p = pStart; p < pEnd; ++p) {
1631:     DMPolytopeType  ct;
1632:     DMPolytopeType *rct;
1633:     PetscInt       *rsize, *rcone, *rornt;
1634:     PetscInt        Nct, n;

1636:     DMPlexGetCellType(dm, p, &ct);
1637:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1638:     for (n = 0; n < Nct; ++n) PetscSectionAddDof(s, p, rsize[n]);
1639:   }
1640:   PetscSectionSetUp(s);
1641:   PetscSectionGetStorageSize(s, &numPointsNew);
1642:   PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets);
1643:   PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf);
1644:   PetscFree(remoteOffsets);
1645:   PetscSFComputeDegreeBegin(sf, &rootdegree);
1646:   PetscSFComputeDegreeEnd(sf, &rootdegree);
1647:   PetscMalloc1(numPointsNew, &rootPointsNew);
1648:   for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1649:   for (p = pStart; p < pEnd; ++p) {
1650:     DMPolytopeType  ct;
1651:     DMPolytopeType *rct;
1652:     PetscInt       *rsize, *rcone, *rornt;
1653:     PetscInt        Nct, n, r, off;

1655:     if (!rootdegree[p - pStart]) continue;
1656:     PetscSectionGetOffset(s, p, &off);
1657:     DMPlexGetCellType(dm, p, &ct);
1658:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1659:     for (n = 0, m = 0; n < Nct; ++n) {
1660:       for (r = 0; r < rsize[n]; ++r, ++m) {
1661:         DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1662:         rootPointsNew[off + m] = pNew;
1663:       }
1664:     }
1665:   }
1666:   PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE);
1667:   PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE);
1668:   PetscSFDestroy(&rsf);
1669:   PetscMalloc1(numLeavesNew, &localPointsNew);
1670:   PetscMalloc1(numLeavesNew, &remotePointsNew);
1671:   for (l = 0, m = 0; l < numLeaves; ++l) {
1672:     const PetscInt  p = localPoints[l];
1673:     DMPolytopeType  ct;
1674:     DMPolytopeType *rct;
1675:     PetscInt       *rsize, *rcone, *rornt;
1676:     PetscInt        Nct, n, r, q, off;

1678:     PetscSectionGetOffset(s, p, &off);
1679:     DMPlexGetCellType(dm, p, &ct);
1680:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1681:     for (n = 0, q = 0; n < Nct; ++n) {
1682:       for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1683:         DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1684:         localPointsNew[m]        = pNew;
1685:         remotePointsNew[m].index = rootPointsNew[off + q];
1686:         remotePointsNew[m].rank  = remotePoints[l].rank;
1687:       }
1688:     }
1689:   }
1690:   PetscSectionDestroy(&s);
1691:   PetscFree(rootPointsNew);
1692:   /* SF needs sorted leaves to correctly calculate Gather */
1693:   {
1694:     PetscSFNode *rp, *rtmp;
1695:     PetscInt    *lp, *idx, *ltmp, i;

1697:     PetscMalloc1(numLeavesNew, &idx);
1698:     PetscMalloc1(numLeavesNew, &lp);
1699:     PetscMalloc1(numLeavesNew, &rp);
1700:     for (i = 0; i < numLeavesNew; ++i) {
1702:       idx[i] = i;
1703:     }
1704:     PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
1705:     for (i = 0; i < numLeavesNew; ++i) {
1706:       lp[i] = localPointsNew[idx[i]];
1707:       rp[i] = remotePointsNew[idx[i]];
1708:     }
1709:     ltmp            = localPointsNew;
1710:     localPointsNew  = lp;
1711:     rtmp            = remotePointsNew;
1712:     remotePointsNew = rp;
1713:     PetscFree(idx);
1714:     PetscFree(ltmp);
1715:     PetscFree(rtmp);
1716:   }
1717:   PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1718:   return 0;
1719: }

1721: /*
1722:   DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.

1724:   Not collective

1726:   Input Parameters:
1727: + cr  - The DMPlexCellRefiner
1728: . ct  - The type of the parent cell
1729: . rct - The type of the produced cell
1730: . r   - The index of the produced cell
1731: - x   - The localized coordinates for the parent cell

1733:   Output Parameter:
1734: . xr  - The localized coordinates for the produced cell

1736:   Level: developer

1738: .seealso: `DMPlexCellRefinerSetCoordinates()`
1739: */
1740: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1741: {
1742:   PetscFE  fe = NULL;
1743:   PetscInt cdim, v, *subcellV;

1745:   DMPlexTransformGetCoordinateFE(tr, ct, &fe);
1746:   DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV);
1747:   PetscFEGetNumComponents(fe, &cdim);
1748:   for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]);
1749:   return 0;
1750: }

1752: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1753: {
1754:   DM                 dm, cdm, cdmCell, cdmNew, cdmCellNew;
1755:   PetscSection       coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
1756:   Vec                coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
1757:   const PetscScalar *coords;
1758:   PetscScalar       *coordsNew;
1759:   const PetscReal   *maxCell, *Lstart, *L;
1760:   PetscBool          localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1761:   PetscInt           dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;

1763:   DMPlexTransformGetDM(tr, &dm);
1764:   DMGetCoordinateDM(dm, &cdm);
1765:   DMGetCellCoordinateDM(dm, &cdmCell);
1766:   DMGetCoordinatesLocalized(dm, &localized);
1767:   DMGetPeriodicity(dm, &maxCell, &Lstart, &L);
1768:   if (localized) {
1769:     /* Localize coordinates of new vertices */
1770:     localizeVertices = PETSC_TRUE;
1771:     /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
1772:     if (!maxCell) localizeCells = PETSC_TRUE;
1773:   }
1774:   DMGetCoordinateSection(dm, &coordSection);
1775:   PetscSectionGetFieldComponents(coordSection, 0, &dEo);
1776:   if (maxCell) {
1777:     PetscReal maxCellNew[3];

1779:     for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
1780:     DMSetPeriodicity(rdm, maxCellNew, Lstart, L);
1781:   }
1782:   DMGetCoordinateDim(rdm, &dE);
1783:   PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew);
1784:   PetscSectionSetNumFields(coordSectionNew, 1);
1785:   PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
1786:   DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
1787:   PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);
1788:   /* Localization should be inherited */
1789:   /*   Stefano calculates parent cells for each new cell for localization */
1790:   /*   Localized cells need coordinates of closure */
1791:   for (v = vStartNew; v < vEndNew; ++v) {
1792:     PetscSectionSetDof(coordSectionNew, v, dE);
1793:     PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
1794:   }
1795:   PetscSectionSetUp(coordSectionNew);
1796:   DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);

1798:   if (localizeCells) {
1799:     DMGetCoordinateDM(rdm, &cdmNew);
1800:     DMClone(cdmNew, &cdmCellNew);
1801:     DMSetCellCoordinateDM(rdm, cdmCellNew);
1802:     DMDestroy(&cdmCellNew);

1804:     PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew);
1805:     PetscSectionSetNumFields(coordSectionCellNew, 1);
1806:     PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE);
1807:     DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew);
1808:     PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew);

1810:     DMGetCellCoordinateSection(dm, &coordSectionCell);
1811:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1812:     for (c = cStart; c < cEnd; ++c) {
1813:       PetscInt dof;

1815:       PetscSectionGetDof(coordSectionCell, c, &dof);
1816:       if (dof) {
1817:         DMPolytopeType  ct;
1818:         DMPolytopeType *rct;
1819:         PetscInt       *rsize, *rcone, *rornt;
1820:         PetscInt        dim, cNew, Nct, n, r;

1822:         DMPlexGetCellType(dm, c, &ct);
1823:         dim = DMPolytopeTypeGetDim(ct);
1824:         DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1825:         /* This allows for different cell types */
1826:         for (n = 0; n < Nct; ++n) {
1827:           if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
1828:           for (r = 0; r < rsize[n]; ++r) {
1829:             PetscInt *closure = NULL;
1830:             PetscInt  clSize, cl, Nv = 0;

1832:             DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew);
1833:             DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1834:             for (cl = 0; cl < clSize * 2; cl += 2) {
1835:               if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
1836:             }
1837:             DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1838:             PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE);
1839:             PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE);
1840:           }
1841:         }
1842:       }
1843:     }
1844:     PetscSectionSetUp(coordSectionCellNew);
1845:     DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew);
1846:   }
1847:   DMViewFromOptions(dm, NULL, "-coarse_dm_view");
1848:   {
1849:     VecType     vtype;
1850:     PetscInt    coordSizeNew, bs;
1851:     const char *name;

1853:     DMGetCoordinatesLocal(dm, &coordsLocal);
1854:     VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
1855:     PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
1856:     VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
1857:     PetscObjectGetName((PetscObject)coordsLocal, &name);
1858:     PetscObjectSetName((PetscObject)coordsLocalNew, name);
1859:     VecGetBlockSize(coordsLocal, &bs);
1860:     VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE);
1861:     VecGetType(coordsLocal, &vtype);
1862:     VecSetType(coordsLocalNew, vtype);
1863:   }
1864:   VecGetArrayRead(coordsLocal, &coords);
1865:   VecGetArray(coordsLocalNew, &coordsNew);
1866:   DMPlexGetChart(dm, &pStart, &pEnd);
1867:   /* First set coordinates for vertices */
1868:   for (p = pStart; p < pEnd; ++p) {
1869:     DMPolytopeType  ct;
1870:     DMPolytopeType *rct;
1871:     PetscInt       *rsize, *rcone, *rornt;
1872:     PetscInt        Nct, n, r;
1873:     PetscBool       hasVertex = PETSC_FALSE;

1875:     DMPlexGetCellType(dm, p, &ct);
1876:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1877:     for (n = 0; n < Nct; ++n) {
1878:       if (rct[n] == DM_POLYTOPE_POINT) {
1879:         hasVertex = PETSC_TRUE;
1880:         break;
1881:       }
1882:     }
1883:     if (hasVertex) {
1884:       const PetscScalar *icoords = NULL;
1885:       const PetscScalar *array   = NULL;
1886:       PetscScalar       *pcoords = NULL;
1887:       PetscBool          isDG;
1888:       PetscInt           Nc, Nv, v, d;

1890:       DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords);

1892:       icoords = pcoords;
1893:       Nv      = Nc / dEo;
1894:       if (ct != DM_POLYTOPE_POINT) {
1895:         if (localizeVertices && maxCell) {
1896:           PetscScalar anchor[3];

1898:           for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
1899:           for (v = 0; v < Nv; ++v) DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]);
1900:         }
1901:       }
1902:       for (n = 0; n < Nct; ++n) {
1903:         if (rct[n] != DM_POLYTOPE_POINT) continue;
1904:         for (r = 0; r < rsize[n]; ++r) {
1905:           PetscScalar vcoords[3];
1906:           PetscInt    vNew, off;

1908:           DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew);
1909:           PetscSectionGetOffset(coordSectionNew, vNew, &off);
1910:           DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords);
1911:           DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]);
1912:         }
1913:       }
1914:       DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords);
1915:     }
1916:   }
1917:   VecRestoreArrayRead(coordsLocal, &coords);
1918:   VecRestoreArray(coordsLocalNew, &coordsNew);
1919:   DMSetCoordinatesLocal(rdm, coordsLocalNew);
1920:   VecDestroy(&coordsLocalNew);
1921:   PetscSectionDestroy(&coordSectionNew);
1922:   /* Then set coordinates for cells by localizing */
1923:   if (!localizeCells) DMLocalizeCoordinates(rdm);
1924:   else {
1925:     VecType     vtype;
1926:     PetscInt    coordSizeNew, bs;
1927:     const char *name;

1929:     DMGetCellCoordinatesLocal(dm, &coordsLocalCell);
1930:     VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew);
1931:     PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew);
1932:     VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE);
1933:     PetscObjectGetName((PetscObject)coordsLocalCell, &name);
1934:     PetscObjectSetName((PetscObject)coordsLocalCellNew, name);
1935:     VecGetBlockSize(coordsLocalCell, &bs);
1936:     VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE);
1937:     VecGetType(coordsLocalCell, &vtype);
1938:     VecSetType(coordsLocalCellNew, vtype);
1939:     VecGetArrayRead(coordsLocalCell, &coords);
1940:     VecGetArray(coordsLocalCellNew, &coordsNew);

1942:     for (p = pStart; p < pEnd; ++p) {
1943:       DMPolytopeType  ct;
1944:       DMPolytopeType *rct;
1945:       PetscInt       *rsize, *rcone, *rornt;
1946:       PetscInt        dof = 0, Nct, n, r;

1948:       DMPlexGetCellType(dm, p, &ct);
1949:       DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1950:       if (p >= cStart && p < cEnd) PetscSectionGetDof(coordSectionCell, p, &dof);
1951:       if (dof) {
1952:         const PetscScalar *pcoords;

1954:         DMPlexPointLocalRead(cdmCell, p, coords, &pcoords);
1955:         for (n = 0; n < Nct; ++n) {
1956:           const PetscInt Nr = rsize[n];

1958:           if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
1959:           for (r = 0; r < Nr; ++r) {
1960:             PetscInt pNew, offNew;

1962:             /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
1963:                DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
1964:                cell to the ones it produces. */
1965:             DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1966:             PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew);
1967:             DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
1968:           }
1969:         }
1970:       }
1971:     }
1972:     VecRestoreArrayRead(coordsLocalCell, &coords);
1973:     VecRestoreArray(coordsLocalCellNew, &coordsNew);
1974:     DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew);
1975:     VecDestroy(&coordsLocalCellNew);
1976:     PetscSectionDestroy(&coordSectionCellNew);
1977:   }
1978:   return 0;
1979: }

1981: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
1982: {
1983:   DM                     rdm;
1984:   DMPlexInterpolatedFlag interp;

1989:   DMPlexTransformSetDM(tr, dm);

1991:   DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
1992:   DMSetType(rdm, DMPLEX);
1993:   DMPlexTransformSetDimensions(tr, dm, rdm);
1994:   /* Calculate number of new points of each depth */
1995:   DMPlexIsInterpolatedCollective(dm, &interp);
1997:   /* Step 1: Set chart */
1998:   DMPlexSetChart(rdm, 0, tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]]);
1999:   /* Step 2: Set cone/support sizes (automatically stratifies) */
2000:   DMPlexTransformSetConeSizes(tr, rdm);
2001:   /* Step 3: Setup refined DM */
2002:   DMSetUp(rdm);
2003:   /* Step 4: Set cones and supports (automatically symmetrizes) */
2004:   DMPlexTransformSetCones(tr, rdm);
2005:   /* Step 5: Create pointSF */
2006:   DMPlexTransformCreateSF(tr, rdm);
2007:   /* Step 6: Create labels */
2008:   DMPlexTransformCreateLabels(tr, rdm);
2009:   /* Step 7: Set coordinates */
2010:   DMPlexTransformSetCoordinates(tr, rdm);
2011:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm);
2012:   *tdm = rdm;
2013:   return 0;
2014: }

2016: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2017: {
2018:   DMPlexTransform tr;
2019:   DM              cdm, rcdm;
2020:   const char     *prefix;

2022:   DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr);
2023:   PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
2024:   PetscObjectSetOptionsPrefix((PetscObject)tr, prefix);
2025:   DMPlexTransformSetDM(tr, dm);
2026:   DMPlexTransformSetFromOptions(tr);
2027:   DMPlexTransformSetActive(tr, adaptLabel);
2028:   DMPlexTransformSetUp(tr);
2029:   PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view");
2030:   DMPlexTransformApply(tr, dm, rdm);
2031:   DMCopyDisc(dm, *rdm);
2032:   DMGetCoordinateDM(dm, &cdm);
2033:   DMGetCoordinateDM(*rdm, &rcdm);
2034:   DMCopyDisc(cdm, rcdm);
2035:   DMPlexTransformCreateDiscLabels(tr, *rdm);
2036:   DMCopyDisc(dm, *rdm);
2037:   DMPlexTransformDestroy(&tr);
2038:   ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2039:   return 0;
2040: }