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:   Sample usage:
 56: .vb
 57:   DMPlexTransformRegister("my_transform", MyTransformCreate);
 58: .ve

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

 70:   Level: advanced

 72:   Note:
 73:   `DMPlexTransformRegister()` may be called multiple times to add several user-defined transforms

 75: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `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 registered `DMPlexTransformType`. It is called from `PetscFinalize()`.

121:   Level: developer

123: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `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:   Level: intermediate

176:   Notes:
177:   See "petsc/include/petscdmplextransform.h" for available transform types

179: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `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 `DMPlexTransformType` name

212:   Level: intermediate

214: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `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 Key:
307: . -dm_plex_transform_type - Set the transform type, e.g. refine_regular

309:   Level: intermediate

311: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `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 `DMPlexTransform`

334:   Collective on tr

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

339:   Level: beginner

341: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `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:   Level: advanced

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

878: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
879: @*/
880: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
881: {
882:   PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
883:   return 0;
884: }

886: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
887: {
888:   *rnew = r;
889:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
890:   return 0;
891: }

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

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

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

1041:   Not Collective

1043:   Input Parameters:
1044: + tr  - The `DMPlexTransform`
1045: . sct - The source point cell type, from whom the new cell is being produced
1046: . sp  - The source point
1047: . so  - The orientation of the source point in its enclosing parent
1048: . tct - The target point cell type
1049: . r   - The replica number requested for the produced cell type
1050: - o   - The orientation of the replica

1052:   Output Parameters:
1053: + rnew - The replica number, given the orientation of the parent
1054: - onew - The replica orientation, given the orientation of the parent

1056:   Level: advanced

1058: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1059: @*/
1060: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1061: {
1063:   PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1064:   return 0;
1065: }

1067: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1068: {
1069:   DM       dm;
1070:   PetscInt pStart, pEnd, p, pNew;

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

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

1097:     DMPlexGetCellTypeLabel(rdm, &ctLabel);
1098:     PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState);
1099:   }
1100:   return 0;
1101: }

1103: #if 0
1104: static PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1105: {
1106:   PetscInt ctNew;

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

1114:     if (q >= ctSN && q < ctEN) break;
1115:   }
1117:   *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType) ctNew);
1118:   return 0;
1119: }
1120: #endif

1122: /* The orientation o is for the interior of the cell p */
1123: 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[])
1124: {
1125:   DM              dm;
1126:   const PetscInt  csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1127:   const PetscInt *cone;
1128:   PetscInt        c, coff = *coneoff, ooff = *orntoff;

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

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

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

1170: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1171: {
1172:   DM             dm;
1173:   DMPolytopeType ct;
1174:   PetscInt      *coneNew, *orntNew;
1175:   PetscInt       maxConeSize = 0, pStart, pEnd, p, pNew;

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

1189:     DMPlexGetCellType(dm, p, &ct);
1190:     DMPlexGetCone(dm, p, &cone);
1191:     DMPlexGetConeOrientation(dm, p, &ornt);
1192:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1193:     for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1194:       const DMPolytopeType ctNew = rct[n];

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

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

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

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

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

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

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

1291: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1292: {
1293:   DM dm;

1297:   DMPlexTransformGetDM(tr, &dm);
1298:   DMRestoreWorkArray(dm, 0, MPIU_INT, cone);
1299:   DMRestoreWorkArray(dm, 0, MPIU_INT, ornt);
1300:   return 0;
1301: }

1303: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1304: {
1305:   PetscInt ict;

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

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

1331:     DMGetDimension(trdm, &trdim);
1332:     DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd);
1333:     tr->trNv[ct] = vEnd - vStart;
1334:     DMGetCoordinatesLocal(trdm, &coordinates);
1335:     VecGetLocalSize(coordinates, &Nc);
1337:     PetscCalloc1(Nc, &tr->trVerts[ct]);
1338:     VecGetArrayRead(coordinates, &coords);
1339:     PetscArraycpy(tr->trVerts[ct], coords, Nc);
1340:     VecRestoreArrayRead(coordinates, &coords);

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

1352:         PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]);
1353:         DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew);
1354:         DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1355:         for (cl = 0; cl < clSize * 2; cl += 2) {
1356:           const PetscInt sv = closure[cl];

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

1369:       PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]);
1370:       for (v = 0; v < tr->trNv[ct]; ++v) {
1371:         PetscPrintf(PETSC_COMM_SELF, "  ");
1372:         for (d = 0; d < dE; ++d) PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]));
1373:         PetscPrintf(PETSC_COMM_SELF, "\n");
1374:       }

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

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

1397:   Input Parameters:
1398: + tr - The DMPlexTransform object
1399: - ct - The cell type

1401:   Output Parameters:
1402: + Nv      - The number of transformed vertices in the closure of the reference cell of given type
1403: - trVerts - The coordinates of these vertices in the reference cell

1405:   Level: developer

1407: .seealso: `DMPlexTransformGetSubcellVertices()`
1408: */
1409: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1410: {
1411:   if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1412:   if (Nv) *Nv = tr->trNv[ct];
1413:   if (trVerts) *trVerts = tr->trVerts[ct];
1414:   return 0;
1415: }

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

1420:   Input Parameters:
1421: + tr  - The DMPlexTransform object
1422: . ct  - The cell type
1423: . rct - The subcell type
1424: - r   - The subcell index

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

1429:   Level: developer

1431: .seealso: `DMPlexTransformGetCellVertices()`
1432: */
1433: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1434: {
1435:   if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1437:   if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1438:   return 0;
1439: }

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

1448:   for (d = 0; d < dE; ++d) out[d] = 0.0;
1449:   for (v = 0; v < Nv; ++v)
1450:     for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1451:   for (d = 0; d < dE; ++d) out[d] /= Nv;
1452:   return 0;
1453: }

1455: /*@
1456:   DMPlexTransformMapCoordinates - Calculate new coordinates for produced points

1458:   Not collective

1460:   Input Parameters:
1461: + tr   - The `DMPlexTransform`
1462: . pct  - The cell type of the parent, from whom the new cell is being produced
1463: . ct   - The type being produced
1464: . p    - The original point
1465: . r    - The replica number requested for the produced cell type
1466: . Nv   - Number of vertices in the closure of the parent cell
1467: . dE   - Spatial dimension
1468: - in   - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell

1470:   Output Parameter:
1471: . out - The coordinates of the new vertices

1473:   Level: intermediate

1475: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransform`, `DMPlexTransformApply()`
1476: @*/
1477: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1478: {
1480:   PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1481:   return 0;
1482: }

1484: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1485: {
1486:   DM              dm;
1487:   IS              valueIS;
1488:   const PetscInt *values;
1489:   PetscInt        defVal, Nv, val;

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

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

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

1532: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1533: {
1534:   DM       dm;
1535:   PetscInt numLabels, l;

1537:   DMPlexTransformGetDM(tr, &dm);
1538:   DMGetNumLabels(dm, &numLabels);
1539:   for (l = 0; l < numLabels; ++l) {
1540:     DMLabel     label, labelNew;
1541:     const char *lname;
1542:     PetscBool   isDepth, isCellType;

1544:     DMGetLabelName(dm, l, &lname);
1545:     PetscStrcmp(lname, "depth", &isDepth);
1546:     if (isDepth) continue;
1547:     PetscStrcmp(lname, "celltype", &isCellType);
1548:     if (isCellType) continue;
1549:     DMCreateLabel(rdm, lname);
1550:     DMGetLabel(dm, lname, &label);
1551:     DMGetLabel(rdm, lname, &labelNew);
1552:     RefineLabel_Internal(tr, label, labelNew);
1553:   }
1554:   return 0;
1555: }

1557: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1558: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1559: {
1560:   DM       dm;
1561:   PetscInt Nf, f, Nds, s;

1563:   DMPlexTransformGetDM(tr, &dm);
1564:   DMGetNumFields(dm, &Nf);
1565:   for (f = 0; f < Nf; ++f) {
1566:     DMLabel     label, labelNew;
1567:     PetscObject obj;
1568:     const char *lname;

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

1583:     DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
1584:     if (!label) continue;
1585:     PetscObjectGetName((PetscObject)label, &lname);
1586:     DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1587:     RefineLabel_Internal(tr, label, labelNew);
1588:     DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
1589:     DMLabelDestroy(&labelNew);
1590:   }
1591:   return 0;
1592: }

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

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

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

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

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

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

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

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

1729:   Not collective

1731:   Input Parameters:
1732: + cr  - The DMPlexCellRefiner
1733: . ct  - The type of the parent cell
1734: . rct - The type of the produced cell
1735: . r   - The index of the produced cell
1736: - x   - The localized coordinates for the parent cell

1738:   Output Parameter:
1739: . xr  - The localized coordinates for the produced cell

1741:   Level: developer

1743: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
1744: @*/
1745: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1746: {
1747:   PetscFE  fe = NULL;
1748:   PetscInt cdim, v, *subcellV;

1750:   DMPlexTransformGetCoordinateFE(tr, ct, &fe);
1751:   DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV);
1752:   PetscFEGetNumComponents(fe, &cdim);
1753:   for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]);
1754:   return 0;
1755: }

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

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

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

1803:   if (localizeCells) {
1804:     DMGetCoordinateDM(rdm, &cdmNew);
1805:     DMClone(cdmNew, &cdmCellNew);
1806:     DMSetCellCoordinateDM(rdm, cdmCellNew);
1807:     DMDestroy(&cdmCellNew);

1809:     PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew);
1810:     PetscSectionSetNumFields(coordSectionCellNew, 1);
1811:     PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE);
1812:     DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew);
1813:     PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew);

1815:     DMGetCellCoordinateSection(dm, &coordSectionCell);
1816:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1817:     for (c = cStart; c < cEnd; ++c) {
1818:       PetscInt dof;

1820:       PetscSectionGetDof(coordSectionCell, c, &dof);
1821:       if (dof) {
1822:         DMPolytopeType  ct;
1823:         DMPolytopeType *rct;
1824:         PetscInt       *rsize, *rcone, *rornt;
1825:         PetscInt        dim, cNew, Nct, n, r;

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

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

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

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

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

1897:       icoords = pcoords;
1898:       Nv      = Nc / dEo;
1899:       if (ct != DM_POLYTOPE_POINT) {
1900:         if (localizeVertices && maxCell) {
1901:           PetscScalar anchor[3];

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

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

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

1947:     for (p = pStart; p < pEnd; ++p) {
1948:       DMPolytopeType  ct;
1949:       DMPolytopeType *rct;
1950:       PetscInt       *rsize, *rcone, *rornt;
1951:       PetscInt        dof = 0, Nct, n, r;

1953:       DMPlexGetCellType(dm, p, &ct);
1954:       DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1955:       if (p >= cStart && p < cEnd) PetscSectionGetDof(coordSectionCell, p, &dof);
1956:       if (dof) {
1957:         const PetscScalar *pcoords;

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

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

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

1986: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
1987: {
1988:   DM                     rdm;
1989:   DMPlexInterpolatedFlag interp;

1994:   DMPlexTransformSetDM(tr, dm);

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

2021: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2022: {
2023:   DMPlexTransform tr;
2024:   DM              cdm, rcdm;
2025:   const char     *prefix;

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