Actual source code: plextransform.c

  1: #include "petsc/private/petscimpl.h"
  2: #include "petscdmplex.h"
  3: #include "petscdmplextransform.h"
  4: #include "petscdmplextransformtypes.h"
  5: #include "petscerror.h"
  6: #include "petscsystypes.h"
  7: #include <petsc/private/dmplextransformimpl.h>

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

 11: PetscClassId DMPLEXTRANSFORM_CLASSID;

 13: PetscFunctionList DMPlexTransformList              = NULL;
 14: PetscBool         DMPlexTransformRegisterAllCalled = PETSC_FALSE;

 16: /* Construct cell type order since we must loop over cell types in depth order */
 17: static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
 18: {
 19:   PetscInt *ctO, *ctOInv;
 20:   PetscInt  c, d, off = 0;

 22:   PetscFunctionBegin;
 23:   PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv));
 24:   for (d = 3; d >= dim; --d) {
 25:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 26:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d || c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue;
 27:       ctO[off++] = c;
 28:     }
 29:   }
 30:   if (dim != 0) {
 31:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 32:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) != 0) continue;
 33:       ctO[off++] = c;
 34:     }
 35:   }
 36:   for (d = dim - 1; d > 0; --d) {
 37:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 38:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d || c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue;
 39:       ctO[off++] = c;
 40:     }
 41:   }
 42:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 43:     if (DMPolytopeTypeGetDim((DMPolytopeType)c) >= 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) continue;
 44:     ctO[off++] = c;
 45:   }
 46:   PetscCheck(off == DM_NUM_POLYTOPES + 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %" PetscInt_FMT " for cell type order", off);
 47:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c;
 48:   *ctOrder    = ctO;
 49:   *ctOrderInv = ctOInv;
 50:   PetscFunctionReturn(PETSC_SUCCESS);
 51: }

 53: /*@C
 54:   DMPlexTransformRegister - Adds a new transform component implementation

 56:   Not Collective

 58:   Input Parameters:
 59: + name        - The name of a new user-defined creation routine
 60: - create_func - The creation routine

 62:   Example Usage:
 63: .vb
 64:   DMPlexTransformRegister("my_transform", MyTransformCreate);
 65: .ve

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

 77:   Level: advanced

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

 82: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()`
 83: @*/
 84: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
 85: {
 86:   PetscFunctionBegin;
 87:   PetscCall(DMInitializePackage());
 88:   PetscCall(PetscFunctionListAdd(&DMPlexTransformList, name, create_func));
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
 93: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
 94: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
 95: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
 96: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
 97: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
 98: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform);
 99: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform);
100: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform);

102: /*@C
103:   DMPlexTransformRegisterAll - Registers all of the transform components in the `DM` package.

105:   Not Collective

107:   Level: advanced

109: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
110: @*/
111: PetscErrorCode DMPlexTransformRegisterAll(void)
112: {
113:   PetscFunctionBegin;
114:   if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
115:   DMPlexTransformRegisterAllCalled = PETSC_TRUE;

117:   PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter));
118:   PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular));
119:   PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox));
120:   PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld));
121:   PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL));
122:   PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR));
123:   PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D));
124:   PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDE, DMPlexTransformCreate_Extrude));
125:   PetscCall(DMPlexTransformRegister(DMPLEXCOHESIVEEXTRUDE, DMPlexTransformCreate_Cohesive));
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: /*@C
130:   DMPlexTransformRegisterDestroy - This function destroys the registered `DMPlexTransformType`. It is called from `PetscFinalize()`.

132:   Not collective

134:   Level: developer

136: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()`
137: @*/
138: PetscErrorCode DMPlexTransformRegisterDestroy(void)
139: {
140:   PetscFunctionBegin;
141:   PetscCall(PetscFunctionListDestroy(&DMPlexTransformList));
142:   DMPlexTransformRegisterAllCalled = PETSC_FALSE;
143:   PetscFunctionReturn(PETSC_SUCCESS);
144: }

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

149:   Collective

151:   Input Parameter:
152: . comm - The communicator for the transform object

154:   Output Parameter:
155: . tr - The transform object

157:   Level: beginner

159: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
160: @*/
161: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
162: {
163:   DMPlexTransform t;

165:   PetscFunctionBegin;
166:   PetscAssertPointer(tr, 2);
167:   *tr = NULL;
168:   PetscCall(DMInitializePackage());

170:   PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView));
171:   t->setupcalled = PETSC_FALSE;
172:   PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom));
173:   *tr = t;
174:   PetscFunctionReturn(PETSC_SUCCESS);
175: }

177: /*@
178:   DMPlexTransformSetType - Sets the particular implementation for a transform.

180:   Collective

182:   Input Parameters:
183: + tr     - The transform
184: - method - The name of the transform type

186:   Options Database Key:
187: . -dm_plex_transform_type <type> - Sets the transform type; see `DMPlexTransformType`

189:   Level: intermediate

191: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
192: @*/
193: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
194: {
195:   PetscErrorCode (*r)(DMPlexTransform);
196:   PetscBool match;

198:   PetscFunctionBegin;
200:   PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match));
201:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

203:   PetscCall(DMPlexTransformRegisterAll());
204:   PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r));
205:   PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);

207:   PetscTryTypeMethod(tr, destroy);
208:   PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops)));
209:   PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method));
210:   PetscCall((*r)(tr));
211:   PetscFunctionReturn(PETSC_SUCCESS);
212: }

214: /*@
215:   DMPlexTransformGetType - Gets the type name (as a string) from the transform.

217:   Not Collective

219:   Input Parameter:
220: . tr - The `DMPlexTransform`

222:   Output Parameter:
223: . type - The `DMPlexTransformType` name

225:   Level: intermediate

227: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
228: @*/
229: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
230: {
231:   PetscFunctionBegin;
233:   PetscAssertPointer(type, 2);
234:   PetscCall(DMPlexTransformRegisterAll());
235:   *type = ((PetscObject)tr)->type_name;
236:   PetscFunctionReturn(PETSC_SUCCESS);
237: }

239: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
240: {
241:   PetscViewerFormat format;

243:   PetscFunctionBegin;
244:   PetscCall(PetscViewerGetFormat(v, &format));
245:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
246:     const PetscInt *trTypes = NULL;
247:     IS              trIS;
248:     PetscInt        cols = 8;
249:     PetscInt        Nrt  = 8, f, g;

251:     if (tr->trType) PetscCall(DMLabelView(tr->trType, v));
252:     PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n"));
253:     for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
254:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
255:     for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]));
256:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
257:     PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n"));
258:     for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
259:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
260:     for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]));
261:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));

263:     if (tr->trType) {
264:       PetscCall(DMLabelGetNumValues(tr->trType, &Nrt));
265:       PetscCall(DMLabelGetValueIS(tr->trType, &trIS));
266:       PetscCall(ISGetIndices(trIS, &trTypes));
267:     }
268:     PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n"));
269:     PetscCall(PetscViewerASCIIPrintf(v, "     "));
270:     for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
271:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
272:     for (f = 0; f < Nrt; ++f) {
273:       PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT "  |", trTypes ? trTypes[f] : f));
274:       for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]));
275:       PetscCall(PetscViewerASCIIPrintf(v, " |\n"));
276:     }
277:     if (tr->trType) {
278:       PetscCall(ISRestoreIndices(trIS, &trTypes));
279:       PetscCall(ISDestroy(&trIS));
280:     }
281:   }
282:   PetscFunctionReturn(PETSC_SUCCESS);
283: }

285: /*@
286:   DMPlexTransformView - Views a `DMPlexTransform`

288:   Collective

290:   Input Parameters:
291: + tr - the `DMPlexTransform` object to view
292: - v  - the viewer

294:   Level: beginner

296: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
297: @*/
298: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
299: {
300:   PetscBool isascii;

302:   PetscFunctionBegin;
304:   if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v));
306:   PetscCheckSameComm(tr, 1, v, 2);
307:   PetscCall(PetscViewerCheckWritable(v));
308:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v));
309:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
310:   if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v));
311:   PetscTryTypeMethod(tr, view, v);
312:   PetscFunctionReturn(PETSC_SUCCESS);
313: }

315: /*@
316:   DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database

318:   Collective

320:   Input Parameter:
321: . tr - the `DMPlexTransform` object to set options for

323:   Options Database Keys:
324: + -dm_plex_transform_type                      - Set the transform type, e.g. refine_regular
325: . -dm_plex_transform_label_match_strata        - Only label points of the same stratum as the producing point
326: . -dm_plex_transform_label_replica_inc <inc>   - Increment for the label value to be multiplied by the replica number, so that the new label value is oldValue + r * inc
327: . -dm_plex_transform_active <name>             - Name for active mesh label
328: - -dm_plex_transform_active_values <v0,v1,...> - Values in the active label

330:   Level: intermediate

332: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
333: @*/
334: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
335: {
336:   char        typeName[1024], active[PETSC_MAX_PATH_LEN];
337:   const char *defName = DMPLEXREFINEREGULAR;
338:   PetscBool   flg, match;

340:   PetscFunctionBegin;
342:   PetscObjectOptionsBegin((PetscObject)tr);
343:   PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg));
344:   if (flg) PetscCall(DMPlexTransformSetType(tr, typeName));
345:   else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName));
346:   PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &match, &flg));
347:   if (flg) PetscCall(DMPlexTransformSetMatchStrata(tr, match));
348:   PetscCall(PetscOptionsInt("-dm_plex_transform_label_replica_inc", "Increment for the label value to be multiplied by the replica number", "", tr->labelReplicaInc, &tr->labelReplicaInc, NULL));
349:   PetscCall(PetscOptionsString("-dm_plex_transform_active", "Name for active mesh label", "DMPlexTransformSetActive", active, active, sizeof(active), &flg));
350:   if (flg) {
351:     DM       dm;
352:     DMLabel  label;
353:     PetscInt values[16];
354:     PetscInt n = 16;

356:     PetscCall(DMPlexTransformGetDM(tr, &dm));
357:     PetscCall(DMGetLabel(dm, active, &label));
358:     PetscCall(PetscOptionsIntArray("-dm_plex_transform_active_values", "The label values to be active", "DMPlexTransformSetActive", values, &n, &flg));
359:     if (flg && n) {
360:       DMLabel newlabel;

362:       PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Active", &newlabel));
363:       for (PetscInt i = 0; i < n; ++i) {
364:         IS is;

366:         PetscCall(DMLabelGetStratumIS(label, values[i], &is));
367:         PetscCall(DMLabelInsertIS(newlabel, is, values[i]));
368:         PetscCall(ISDestroy(&is));
369:       }
370:       PetscCall(DMPlexTransformSetActive(tr, newlabel));
371:       PetscCall(DMLabelDestroy(&newlabel));
372:     } else {
373:       PetscCall(DMPlexTransformSetActive(tr, label));
374:     }
375:   }
376:   PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
377:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
378:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject));
379:   PetscOptionsEnd();
380:   PetscFunctionReturn(PETSC_SUCCESS);
381: }

383: /*@
384:   DMPlexTransformDestroy - Destroys a `DMPlexTransform`

386:   Collective

388:   Input Parameter:
389: . tr - the transform object to destroy

391:   Level: beginner

393: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
394: @*/
395: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
396: {
397:   PetscInt c;

399:   PetscFunctionBegin;
400:   if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
402:   if (--((PetscObject)*tr)->refct > 0) {
403:     *tr = NULL;
404:     PetscFunctionReturn(PETSC_SUCCESS);
405:   }

407:   PetscTryTypeMethod(*tr, destroy);
408:   PetscCall(DMDestroy(&(*tr)->dm));
409:   PetscCall(DMLabelDestroy(&(*tr)->active));
410:   PetscCall(DMLabelDestroy(&(*tr)->trType));
411:   PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld));
412:   PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew));
413:   PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew));
414:   PetscCall(PetscFree((*tr)->offset));
415:   PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd));
416:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
417:     PetscCall(PetscFEDestroy(&(*tr)->coordFE[c]));
418:     PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c]));
419:   }
420:   if ((*tr)->trVerts) {
421:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
422:       DMPolytopeType *rct;
423:       PetscInt       *rsize, *rcone, *rornt, Nct, n, r;

425:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) {
426:         PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
427:         for (n = 0; n < Nct; ++n) {
428:           if (rct[n] == DM_POLYTOPE_POINT) continue;
429:           for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r]));
430:           PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]]));
431:         }
432:       }
433:       PetscCall(PetscFree((*tr)->trSubVerts[c]));
434:       PetscCall(PetscFree((*tr)->trVerts[c]));
435:     }
436:   }
437:   PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts));
438:   PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom));
439:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
440:   PetscCall(PetscHeaderDestroy(tr));
441:   PetscFunctionReturn(PETSC_SUCCESS);
442: }

444: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
445: {
446:   DMLabel  trType = tr->trType;
447:   PetscInt c, cN, *off;

449:   PetscFunctionBegin;
450:   if (trType) {
451:     DM              dm;
452:     IS              rtIS;
453:     const PetscInt *reftypes;
454:     PetscInt        Nrt, r;

456:     PetscCall(DMPlexTransformGetDM(tr, &dm));
457:     PetscCall(DMLabelGetNumValues(trType, &Nrt));
458:     PetscCall(DMLabelGetValueIS(trType, &rtIS));
459:     PetscCall(ISGetIndices(rtIS, &reftypes));
460:     PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off));
461:     for (r = 0; r < Nrt; ++r) {
462:       const PetscInt  rt = reftypes[r];
463:       IS              rtIS;
464:       const PetscInt *points;
465:       DMPolytopeType  ct;
466:       PetscInt        np, p;

468:       PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
469:       PetscCall(ISGetLocalSize(rtIS, &np));
470:       PetscCall(ISGetIndices(rtIS, &points));
471:       if (!np) continue;
472:       p = points[0];
473:       PetscCall(ISRestoreIndices(rtIS, &points));
474:       PetscCall(ISDestroy(&rtIS));
475:       PetscCall(DMPlexGetCellType(dm, p, &ct));
476:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
477:         const DMPolytopeType ctNew = (DMPolytopeType)cN;
478:         DMPolytopeType      *rct;
479:         PetscInt            *rsize, *cone, *ornt;
480:         PetscInt             Nct, n, s;

482:         if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
483:           off[r * DM_NUM_POLYTOPES + ctNew] = -1;
484:           break;
485:         }
486:         off[r * DM_NUM_POLYTOPES + ctNew] = 0;
487:         for (s = 0; s <= r; ++s) {
488:           const PetscInt st = reftypes[s];
489:           DMPolytopeType sct;
490:           PetscInt       q, qrt;

492:           PetscCall(DMLabelGetStratumIS(trType, st, &rtIS));
493:           PetscCall(ISGetLocalSize(rtIS, &np));
494:           PetscCall(ISGetIndices(rtIS, &points));
495:           if (!np) continue;
496:           q = points[0];
497:           PetscCall(ISRestoreIndices(rtIS, &points));
498:           PetscCall(ISDestroy(&rtIS));
499:           PetscCall(DMPlexGetCellType(dm, q, &sct));
500:           PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt));
501:           PetscCheck(st == qrt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Refine type %" PetscInt_FMT " of point %" PetscInt_FMT " does not match predicted type %" PetscInt_FMT, qrt, q, st);
502:           if (st == rt) {
503:             for (n = 0; n < Nct; ++n)
504:               if (rct[n] == ctNew) break;
505:             if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
506:             break;
507:           }
508:           for (n = 0; n < Nct; ++n) {
509:             if (rct[n] == ctNew) {
510:               PetscInt sn;

512:               PetscCall(DMLabelGetStratumSize(trType, st, &sn));
513:               off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
514:             }
515:           }
516:         }
517:       }
518:     }
519:     PetscCall(ISRestoreIndices(rtIS, &reftypes));
520:     PetscCall(ISDestroy(&rtIS));
521:   } else {
522:     PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off));
523:     for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
524:       const DMPolytopeType ct = (DMPolytopeType)c;
525:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
526:         const DMPolytopeType ctNew = (DMPolytopeType)cN;
527:         DMPolytopeType      *rct;
528:         PetscInt            *rsize, *cone, *ornt;
529:         PetscInt             Nct, n, i;

531:         if (DMPolytopeTypeGetDim(ct) < 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE || DMPolytopeTypeGetDim(ctNew) < 0 || ctNew == DM_POLYTOPE_UNKNOWN_CELL || ctNew == DM_POLYTOPE_UNKNOWN_FACE) {
532:           off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
533:           continue;
534:         }
535:         off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
536:         for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
537:           const DMPolytopeType ict  = (DMPolytopeType)ctOrderOld[i];
538:           const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];

540:           PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt));
541:           if (ict == ct) {
542:             for (n = 0; n < Nct; ++n)
543:               if (rct[n] == ctNew) break;
544:             if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
545:             break;
546:           }
547:           for (n = 0; n < Nct; ++n)
548:             if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
549:         }
550:       }
551:     }
552:   }
553:   *offset = off;
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }

557: /*@
558:   DMPlexTransformSetUp - Create the tables that drive the transform

560:   Input Parameter:
561: . tr - The `DMPlexTransform` object

563:   Level: intermediate

565: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
566: @*/
567: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
568: {
569:   DM             dm;
570:   DMPolytopeType ctCell;
571:   PetscInt       pStart, pEnd, p, c, celldim = 0;

573:   PetscFunctionBegin;
575:   if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
576:   PetscTryTypeMethod(tr, setup);
577:   PetscCall(DMPlexTransformGetDM(tr, &dm));
578:   PetscCall(DMSetSnapToGeomModel(dm, NULL));
579:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
580:   if (pEnd > pStart) {
581:     // Ignore cells hanging off of embedded surfaces
582:     PetscInt c = pStart;

584:     ctCell = DM_POLYTOPE_FV_GHOST;
585:     while (DMPolytopeTypeGetDim(ctCell) < 0) PetscCall(DMPlexGetCellType(dm, c++, &ctCell));
586:   } else {
587:     PetscInt dim;

589:     PetscCall(DMGetDimension(dm, &dim));
590:     switch (dim) {
591:     case 0:
592:       ctCell = DM_POLYTOPE_POINT;
593:       break;
594:     case 1:
595:       ctCell = DM_POLYTOPE_SEGMENT;
596:       break;
597:     case 2:
598:       ctCell = DM_POLYTOPE_TRIANGLE;
599:       break;
600:     case 3:
601:       ctCell = DM_POLYTOPE_TETRAHEDRON;
602:       break;
603:     default:
604:       ctCell = DM_POLYTOPE_UNKNOWN;
605:     }
606:   }
607:   PetscCall(DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
608:   for (p = pStart; p < pEnd; ++p) {
609:     DMPolytopeType  ct;
610:     DMPolytopeType *rct;
611:     PetscInt       *rsize, *cone, *ornt;
612:     PetscInt        Nct, n;

614:     PetscCall(DMPlexGetCellType(dm, p, &ct));
615:     PetscCheck(ct != DM_POLYTOPE_UNKNOWN && ct != DM_POLYTOPE_UNKNOWN_CELL && ct != DM_POLYTOPE_UNKNOWN_FACE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p);
616:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
617:     for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
618:   }
619:   PetscCall(DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
620:   /* Construct sizes and offsets for each cell type */
621:   if (!tr->ctStart) {
622:     PetscInt *ctS, *ctSN, *ctC, *ctCN;

624:     PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
625:     PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
626:     for (p = pStart; p < pEnd; ++p) {
627:       DMPolytopeType  ct;
628:       DMPolytopeType *rct;
629:       PetscInt       *rsize, *cone, *ornt;
630:       PetscInt        Nct, n;

632:       PetscCall(DMPlexGetCellType(dm, p, &ct));
633:       PetscCheck(ct != DM_POLYTOPE_UNKNOWN && ct != DM_POLYTOPE_UNKNOWN_CELL && ct != DM_POLYTOPE_UNKNOWN_FACE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p);
634:       ++ctC[ct];
635:       PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
636:       for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
637:     }
638:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
639:       const PetscInt cto  = tr->ctOrderOld[c];
640:       const PetscInt cton = tr->ctOrderOld[c + 1];
641:       const PetscInt ctn  = tr->ctOrderNew[c];
642:       const PetscInt ctnn = tr->ctOrderNew[c + 1];

644:       ctS[cton]  = ctS[cto] + ctC[cto];
645:       ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
646:     }
647:     PetscCall(PetscFree2(ctC, ctCN));
648:     tr->ctStart    = ctS;
649:     tr->ctStartNew = ctSN;
650:   }
651:   PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
652:   // Compute depth information
653:   tr->depth = -1;
654:   for (c = 0; c < DM_NUM_POLYTOPES; ++c)
655:     if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
656:   PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
657:   for (PetscInt d = 0; d <= tr->depth; ++d) {
658:     tr->depthStart[d] = PETSC_INT_MAX;
659:     tr->depthEnd[d]   = -1;
660:   }
661:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
662:     const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);

664:     if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
665:     tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
666:     tr->depthEnd[dep]   = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
667:   }
668:   tr->setupcalled = PETSC_TRUE;
669:   PetscFunctionReturn(PETSC_SUCCESS);
670: }

672: /*@
673:   DMPlexTransformGetDM - Get the base `DM` for the transform

675:   Input Parameter:
676: . tr - The `DMPlexTransform` object

678:   Output Parameter:
679: . dm - The original `DM` which will be transformed

681:   Level: intermediate

683: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
684: @*/
685: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
686: {
687:   PetscFunctionBegin;
689:   PetscAssertPointer(dm, 2);
690:   *dm = tr->dm;
691:   PetscFunctionReturn(PETSC_SUCCESS);
692: }

694: /*@
695:   DMPlexTransformSetDM - Set the base `DM` for the transform

697:   Input Parameters:
698: + tr - The `DMPlexTransform` object
699: - dm - The original `DM` which will be transformed

701:   Level: intermediate

703:   Note:
704:   The user does not typically call this, as it is called by `DMPlexTransformApply()`.

706: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
707: @*/
708: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
709: {
710:   PetscFunctionBegin;
713:   PetscCall(PetscObjectReference((PetscObject)dm));
714:   PetscCall(DMDestroy(&tr->dm));
715:   tr->dm = dm;
716:   PetscFunctionReturn(PETSC_SUCCESS);
717: }

719: /*@
720:   DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform

722:   Input Parameter:
723: . tr - The `DMPlexTransform` object

725:   Output Parameter:
726: . active - The `DMLabel` indicating which points will be transformed

728:   Level: intermediate

730: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
731: @*/
732: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
733: {
734:   PetscFunctionBegin;
736:   PetscAssertPointer(active, 2);
737:   *active = tr->active;
738:   PetscFunctionReturn(PETSC_SUCCESS);
739: }

741: /*@
742:   DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform

744:   Input Parameters:
745: + tr     - The `DMPlexTransform` object
746: - active - The original `DM` which will be transformed

748:   Level: intermediate

750:   Note:
751:   This only applies to transforms that can operator on a subset of the mesh, listed in [](plex_transform_table).

753: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
754: @*/
755: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
756: {
757:   PetscFunctionBegin;
760:   PetscCall(PetscObjectReference((PetscObject)active));
761:   PetscCall(DMLabelDestroy(&tr->active));
762:   tr->active = active;
763:   PetscFunctionReturn(PETSC_SUCCESS);
764: }

766: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
767: {
768:   PetscFunctionBegin;
769:   if (!tr->coordFE[ct]) {
770:     PetscInt dim, cdim;

772:     dim = DMPolytopeTypeGetDim(ct);
773:     PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
774:     PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
775:     {
776:       PetscDualSpace  dsp;
777:       PetscQuadrature quad;
778:       DM              K;
779:       PetscFEGeom    *cg;
780:       PetscScalar    *Xq;
781:       PetscReal      *xq, *wq;
782:       PetscInt        Nq, q;

784:       PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
785:       PetscCall(PetscMalloc1(Nq * cdim, &xq));
786:       for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
787:       PetscCall(PetscMalloc1(Nq, &wq));
788:       for (q = 0; q < Nq; ++q) wq[q] = 1.0;
789:       PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
790:       PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
791:       PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));

793:       PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
794:       PetscCall(PetscDualSpaceGetDM(dsp, &K));
795:       PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]));
796:       cg = tr->refGeom[ct];
797:       PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
798:       PetscCall(PetscQuadratureDestroy(&quad));
799:     }
800:   }
801:   *fe = tr->coordFE[ct];
802:   PetscFunctionReturn(PETSC_SUCCESS);
803: }

805: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
806: {
807:   PetscInt dim, cdim;

809:   PetscFunctionBegin;
810:   PetscCall(DMGetDimension(dm, &dim));
811:   PetscCall(DMSetDimension(tdm, dim));
812:   PetscCall(DMGetCoordinateDim(dm, &cdim));
813:   PetscCall(DMSetCoordinateDim(tdm, cdim));
814:   PetscFunctionReturn(PETSC_SUCCESS);
815: }

817: /*@
818:   DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`

820:   Input Parameters:
821: + tr - The `DMPlexTransform` object
822: - dm - The original `DM`

824:   Output Parameter:
825: . tdm - The transformed `DM`

827:   Level: advanced

829: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
830: @*/
831: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
832: {
833:   PetscFunctionBegin;
834:   PetscUseTypeMethod(tr, setdimensions, dm, tdm);
835:   PetscFunctionReturn(PETSC_SUCCESS);
836: }

838: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
839: {
840:   PetscFunctionBegin;
841:   if (pStart) *pStart = 0;
842:   if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
843:   PetscFunctionReturn(PETSC_SUCCESS);
844: }

846: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
847: {
848:   PetscInt ctNew;

850:   PetscFunctionBegin;
852:   PetscAssertPointer(celltype, 3);
853:   /* TODO Can do bisection since everything is sorted */
854:   for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
855:     PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];

857:     if (cell >= ctSN && cell < ctEN) break;
858:   }
859:   PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
860:   *celltype = (DMPolytopeType)ctNew;
861:   PetscFunctionReturn(PETSC_SUCCESS);
862: }

864: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
865: {
866:   PetscFunctionBegin;
868:   if (start) *start = tr->ctStartNew[celltype];
869:   if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
870:   PetscFunctionReturn(PETSC_SUCCESS);
871: }

873: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
874: {
875:   PetscFunctionBegin;
877:   *depth = tr->depth;
878:   PetscFunctionReturn(PETSC_SUCCESS);
879: }

881: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
882: {
883:   PetscFunctionBegin;
885:   if (start) *start = tr->depthStart[depth];
886:   if (end) *end = tr->depthEnd[depth];
887:   PetscFunctionReturn(PETSC_SUCCESS);
888: }

890: /*@
891:   DMPlexTransformGetMatchStrata - Get the flag which determines what points get added to the transformed labels

893:   Not Collective

895:   Input Parameter:
896: . tr - The `DMPlexTransform`

898:   Output Parameter:
899: . match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels

901:   Level: intermediate

903: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetMatchStrata()`, `DMPlexGetPointDepth()`
904: @*/
905: PetscErrorCode DMPlexTransformGetMatchStrata(DMPlexTransform tr, PetscBool *match)
906: {
907:   PetscFunctionBegin;
909:   PetscAssertPointer(match, 2);
910:   *match = tr->labelMatchStrata;
911:   PetscFunctionReturn(PETSC_SUCCESS);
912: }

914: /*@
915:   DMPlexTransformSetMatchStrata - Set the flag which determines what points get added to the transformed labels

917:   Not Collective

919:   Input Parameters:
920: + tr    - The `DMPlexTransform`
921: - match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels

923:   Level: intermediate

925: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetMatchStrata()`, `DMPlexGetPointDepth()`
926: @*/
927: PetscErrorCode DMPlexTransformSetMatchStrata(DMPlexTransform tr, PetscBool match)
928: {
929:   PetscFunctionBegin;
931:   tr->labelMatchStrata = match;
932:   PetscFunctionReturn(PETSC_SUCCESS);
933: }

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

938:   Not Collective

940:   Input Parameters:
941: + tr    - The `DMPlexTransform`
942: . ct    - The type of the original point which produces the new point
943: . ctNew - The type of the new point
944: . p     - The original point which produces the new point
945: - r     - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`

947:   Output Parameter:
948: . pNew - The new point number

950:   Level: developer

952: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
953: @*/
954: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
955: {
956:   DMPolytopeType *rct;
957:   PetscInt       *rsize, *cone, *ornt;
958:   PetscInt        rt, Nct, n, off, rp;
959:   DMLabel         trType = tr->trType;
960:   PetscInt        ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
961:   PetscInt        ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
962:   PetscInt        newp = ctSN, cind;

964:   PetscFunctionBeginHot;
965:   PetscCheck(p >= ctS && p < ctE, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " is not a %s [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, DMPolytopeTypes[ct], ctS, ctE);
966:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
967:   if (trType) {
968:     PetscCall(DMLabelGetValueIndex(trType, rt, &cind));
969:     PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
970:     PetscCheck(rp >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s point %" PetscInt_FMT " does not have refine type %" PetscInt_FMT, DMPolytopeTypes[ct], p, rt);
971:   } else {
972:     cind = ct;
973:     rp   = p - ctS;
974:   }
975:   off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
976:   PetscCheck(off >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s (%" PetscInt_FMT ") of point %" PetscInt_FMT " does not produce type %s for transform %s", DMPolytopeTypes[ct], rt, p, DMPolytopeTypes[ctNew], tr->hdr.type_name);
977:   newp += off;
978:   for (n = 0; n < Nct; ++n) {
979:     if (rct[n] == ctNew) {
980:       if (rsize[n] && r >= rsize[n])
981:         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]);
982:       newp += rp * rsize[n] + r;
983:       break;
984:     }
985:   }

987:   PetscCheck(!(newp < ctSN) && !(newp >= ctEN), PETSC_COMM_SELF, PETSC_ERR_PLIB, "New point %" PetscInt_FMT " is not a %s [%" PetscInt_FMT ", %" PetscInt_FMT ")", newp, DMPolytopeTypes[ctNew], ctSN, ctEN);
988:   *pNew = newp;
989:   PetscFunctionReturn(PETSC_SUCCESS);
990: }

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

995:   Not Collective

997:   Input Parameters:
998: + tr   - The `DMPlexTransform`
999: - pNew - The new point number

1001:   Output Parameters:
1002: + ct    - The type of the original point which produces the new point
1003: . ctNew - The type of the new point
1004: . p     - The original point which produces the new point
1005: - r     - The replica number of the new point, meaning it is the rth point of type ctNew produced from p

1007:   Level: developer

1009: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
1010: @*/
1011: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
1012: {
1013:   DMLabel         trType = tr->trType;
1014:   DMPolytopeType *rct, ctN;
1015:   PetscInt       *rsize, *cone, *ornt;
1016:   PetscInt        rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
1017:   PetscInt        offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;

1019:   PetscFunctionBegin;
1020:   PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
1021:   if (trType) {
1022:     DM              dm;
1023:     IS              rtIS;
1024:     const PetscInt *reftypes;
1025:     PetscInt        Nrt, r, rtStart;

1027:     PetscCall(DMPlexTransformGetDM(tr, &dm));
1028:     PetscCall(DMLabelGetNumValues(trType, &Nrt));
1029:     PetscCall(DMLabelGetValueIS(trType, &rtIS));
1030:     PetscCall(ISGetIndices(rtIS, &reftypes));
1031:     for (r = 0; r < Nrt; ++r) {
1032:       const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];

1034:       if (tr->ctStartNew[ctN] + off > pNew) continue;
1035:       /* Check that any of this refinement type exist */
1036:       /* TODO Actually keep track of the number produced here instead */
1037:       if (off > offset) {
1038:         rt     = reftypes[r];
1039:         offset = off;
1040:       }
1041:     }
1042:     PetscCall(ISRestoreIndices(rtIS, &reftypes));
1043:     PetscCall(ISDestroy(&rtIS));
1044:     PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
1045:     /* TODO Map refinement types to cell types */
1046:     PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL));
1047:     PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
1048:     for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
1049:       PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];

1051:       if ((rtStart >= ctS) && (rtStart < ctE)) break;
1052:     }
1053:     PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt);
1054:   } else {
1055:     for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
1056:       const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];

1058:       if (tr->ctStartNew[ctN] + off > pNew) continue;
1059:       if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
1060:       /* TODO Actually keep track of the number produced here instead */
1061:       if (off > offset) {
1062:         ctO    = ctTmp;
1063:         offset = off;
1064:       }
1065:     }
1066:     rt = -1;
1067:     PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
1068:   }
1069:   ctS = tr->ctStart[ctO];
1070:   ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
1071:   if (trType) {
1072:     for (rtS = ctS; rtS < ctE; ++rtS) {
1073:       PetscInt val;
1074:       PetscCall(DMLabelGetValue(trType, rtS, &val));
1075:       if (val == rt) break;
1076:     }
1077:     PetscCheck(rtS < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find point of type %s with refine type %" PetscInt_FMT, DMPolytopeTypes[ctO], rt);
1078:   } else rtS = ctS;
1079:   PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
1080:   PetscCheck(!trType || rt == rtTmp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " has refine type %" PetscInt_FMT " != %" PetscInt_FMT " refine type which produced point %" PetscInt_FMT, rtS, rtTmp, rt, pNew);
1081:   for (n = 0; n < Nct; ++n) {
1082:     if (rct[n] == ctN) {
1083:       PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c;

1085:       if (trType) {
1086:         for (c = ctS; c < ctE; ++c) {
1087:           PetscCall(DMLabelGetValue(trType, c, &val));
1088:           if (val == rt) {
1089:             if (tmp < rsize[n]) break;
1090:             tmp -= rsize[n];
1091:           }
1092:         }
1093:         PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew);
1094:         rp = c - ctS;
1095:         rO = tmp;
1096:       } else {
1097:         // This assumes that all points of type ctO transform the same way
1098:         rp = tmp / rsize[n];
1099:         rO = tmp % rsize[n];
1100:       }
1101:       break;
1102:     }
1103:   }
1104:   PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
1105:   pO = rp + ctS;
1106:   PetscCheck(!(pO < ctS) && !(pO >= ctE), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Source point %" PetscInt_FMT " is not a %s [%" PetscInt_FMT ", %" PetscInt_FMT ")", pO, DMPolytopeTypes[ctO], ctS, ctE);
1107:   if (ct) *ct = (DMPolytopeType)ctO;
1108:   if (ctNew) *ctNew = ctN;
1109:   if (p) *p = pO;
1110:   if (r) *r = rO;
1111:   PetscFunctionReturn(PETSC_SUCCESS);
1112: }

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

1117:   Input Parameters:
1118: + tr     - The `DMPlexTransform` object
1119: . source - The source cell type
1120: - p      - The source point, which can also determine the refine type

1122:   Output Parameters:
1123: + rt     - The refine type for this point
1124: . Nt     - The number of types produced by this point
1125: . target - An array of length `Nt` giving the types produced
1126: . size   - An array of length `Nt` giving the number of cells of each type produced
1127: . cone   - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
1128: - ornt   - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point

1130:   Level: advanced

1132:   Notes:
1133:   The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
1134:   need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
1135:   division (described in these outputs) of the cell in the original mesh. The point identifier is given by
1136: .vb
1137:    the number of cones to be taken, or 0 for the current cell
1138:    the cell cone point number at each level from which it is subdivided
1139:    the replica number r of the subdivision.
1140: .ve
1141:   The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
1142: .vb
1143:    Nt     = 2
1144:    target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
1145:    size   = {1, 2}
1146:    cone   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,  DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
1147:    ornt   = {                         0,                       0,                        0,                          0}
1148: .ve

1150: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
1151: @*/
1152: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1153: {
1154:   PetscFunctionBegin;
1155:   PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
1156:   PetscFunctionReturn(PETSC_SUCCESS);
1157: }

1159: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1160: {
1161:   PetscFunctionBegin;
1162:   *rnew = r;
1163:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
1164:   PetscFunctionReturn(PETSC_SUCCESS);
1165: }

1167: /* Returns the same thing */
1168: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1169: {
1170:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
1171:   static PetscInt       vertexS[] = {1};
1172:   static PetscInt       vertexC[] = {0};
1173:   static PetscInt       vertexO[] = {0};
1174:   static DMPolytopeType edgeT[]   = {DM_POLYTOPE_SEGMENT};
1175:   static PetscInt       edgeS[]   = {1};
1176:   static PetscInt       edgeC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1177:   static PetscInt       edgeO[]   = {0, 0};
1178:   static DMPolytopeType tedgeT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR};
1179:   static PetscInt       tedgeS[]  = {1};
1180:   static PetscInt       tedgeC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1181:   static PetscInt       tedgeO[]  = {0, 0};
1182:   static DMPolytopeType triT[]    = {DM_POLYTOPE_TRIANGLE};
1183:   static PetscInt       triS[]    = {1};
1184:   static PetscInt       triC[]    = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1185:   static PetscInt       triO[]    = {0, 0, 0};
1186:   static DMPolytopeType quadT[]   = {DM_POLYTOPE_QUADRILATERAL};
1187:   static PetscInt       quadS[]   = {1};
1188:   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};
1189:   static PetscInt       quadO[]   = {0, 0, 0, 0};
1190:   static DMPolytopeType tquadT[]  = {DM_POLYTOPE_SEG_PRISM_TENSOR};
1191:   static PetscInt       tquadS[]  = {1};
1192:   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};
1193:   static PetscInt       tquadO[]  = {0, 0, 0, 0};
1194:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_TETRAHEDRON};
1195:   static PetscInt       tetS[]    = {1};
1196:   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};
1197:   static PetscInt       tetO[]    = {0, 0, 0, 0};
1198:   static DMPolytopeType hexT[]    = {DM_POLYTOPE_HEXAHEDRON};
1199:   static PetscInt       hexS[]    = {1};
1200:   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};
1201:   static PetscInt       hexO[] = {0, 0, 0, 0, 0, 0};
1202:   static DMPolytopeType tripT[]   = {DM_POLYTOPE_TRI_PRISM};
1203:   static PetscInt       tripS[]   = {1};
1204:   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};
1205:   static PetscInt       tripO[]   = {0, 0, 0, 0, 0};
1206:   static DMPolytopeType ttripT[]  = {DM_POLYTOPE_TRI_PRISM_TENSOR};
1207:   static PetscInt       ttripS[]  = {1};
1208:   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};
1209:   static PetscInt       ttripO[]  = {0, 0, 0, 0, 0};
1210:   static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
1211:   static PetscInt       tquadpS[] = {1};
1212:   static PetscInt       tquadpC[] = {DM_POLYTOPE_QUADRILATERAL,    1, 0, 0, DM_POLYTOPE_QUADRILATERAL,    1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
1213:                                      DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1214:   static PetscInt       tquadpO[] = {0, 0, 0, 0, 0, 0};
1215:   static DMPolytopeType pyrT[]    = {DM_POLYTOPE_PYRAMID};
1216:   static PetscInt       pyrS[]    = {1};
1217:   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};
1218:   static PetscInt       pyrO[]    = {0, 0, 0, 0, 0};

1220:   PetscFunctionBegin;
1221:   if (rt) *rt = 0;
1222:   switch (source) {
1223:   case DM_POLYTOPE_POINT:
1224:     *Nt     = 1;
1225:     *target = vertexT;
1226:     *size   = vertexS;
1227:     *cone   = vertexC;
1228:     *ornt   = vertexO;
1229:     break;
1230:   case DM_POLYTOPE_SEGMENT:
1231:     *Nt     = 1;
1232:     *target = edgeT;
1233:     *size   = edgeS;
1234:     *cone   = edgeC;
1235:     *ornt   = edgeO;
1236:     break;
1237:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1238:     *Nt     = 1;
1239:     *target = tedgeT;
1240:     *size   = tedgeS;
1241:     *cone   = tedgeC;
1242:     *ornt   = tedgeO;
1243:     break;
1244:   case DM_POLYTOPE_TRIANGLE:
1245:     *Nt     = 1;
1246:     *target = triT;
1247:     *size   = triS;
1248:     *cone   = triC;
1249:     *ornt   = triO;
1250:     break;
1251:   case DM_POLYTOPE_QUADRILATERAL:
1252:     *Nt     = 1;
1253:     *target = quadT;
1254:     *size   = quadS;
1255:     *cone   = quadC;
1256:     *ornt   = quadO;
1257:     break;
1258:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1259:     *Nt     = 1;
1260:     *target = tquadT;
1261:     *size   = tquadS;
1262:     *cone   = tquadC;
1263:     *ornt   = tquadO;
1264:     break;
1265:   case DM_POLYTOPE_TETRAHEDRON:
1266:     *Nt     = 1;
1267:     *target = tetT;
1268:     *size   = tetS;
1269:     *cone   = tetC;
1270:     *ornt   = tetO;
1271:     break;
1272:   case DM_POLYTOPE_HEXAHEDRON:
1273:     *Nt     = 1;
1274:     *target = hexT;
1275:     *size   = hexS;
1276:     *cone   = hexC;
1277:     *ornt   = hexO;
1278:     break;
1279:   case DM_POLYTOPE_TRI_PRISM:
1280:     *Nt     = 1;
1281:     *target = tripT;
1282:     *size   = tripS;
1283:     *cone   = tripC;
1284:     *ornt   = tripO;
1285:     break;
1286:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
1287:     *Nt     = 1;
1288:     *target = ttripT;
1289:     *size   = ttripS;
1290:     *cone   = ttripC;
1291:     *ornt   = ttripO;
1292:     break;
1293:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1294:     *Nt     = 1;
1295:     *target = tquadpT;
1296:     *size   = tquadpS;
1297:     *cone   = tquadpC;
1298:     *ornt   = tquadpO;
1299:     break;
1300:   case DM_POLYTOPE_PYRAMID:
1301:     *Nt     = 1;
1302:     *target = pyrT;
1303:     *size   = pyrS;
1304:     *cone   = pyrC;
1305:     *ornt   = pyrO;
1306:     break;
1307:   default:
1308:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1309:   }
1310:   PetscFunctionReturn(PETSC_SUCCESS);
1311: }

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

1316:   Not Collective

1318:   Input Parameters:
1319: + tr  - The `DMPlexTransform`
1320: . sct - The source point cell type, from whom the new cell is being produced
1321: . sp  - The source point
1322: . so  - The orientation of the source point in its enclosing parent
1323: . tct - The target point cell type
1324: . r   - The replica number requested for the produced cell type
1325: - o   - The orientation of the replica

1327:   Output Parameters:
1328: + rnew - The replica number, given the orientation of the parent
1329: - onew - The replica orientation, given the orientation of the parent

1331:   Level: advanced

1333: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1334: @*/
1335: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1336: {
1337:   PetscFunctionBeginHot;
1338:   PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1339:   PetscFunctionReturn(PETSC_SUCCESS);
1340: }

1342: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1343: {
1344:   DM       dm;
1345:   PetscInt pStart, pEnd, p, pNew;

1347:   PetscFunctionBegin;
1348:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1349:   /* Must create the celltype label here so that we do not automatically try to compute the types */
1350:   PetscCall(DMCreateLabel(rdm, "celltype"));
1351:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1352:   for (p = pStart; p < pEnd; ++p) {
1353:     DMPolytopeType  ct;
1354:     DMPolytopeType *rct;
1355:     PetscInt       *rsize, *rcone, *rornt;
1356:     PetscInt        Nct, n, r;

1358:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1359:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1360:     for (n = 0; n < Nct; ++n) {
1361:       for (r = 0; r < rsize[n]; ++r) {
1362:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1363:         PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
1364:         PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
1365:       }
1366:     }
1367:   }
1368:   /* Let the DM know we have set all the cell types */
1369:   {
1370:     DMLabel  ctLabel;
1371:     DM_Plex *plex = (DM_Plex *)rdm->data;

1373:     PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
1374:     PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
1375:   }
1376:   PetscFunctionReturn(PETSC_SUCCESS);
1377: }

1379: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1380: {
1381:   DMPolytopeType ctNew;

1383:   PetscFunctionBegin;
1385:   PetscAssertPointer(coneSize, 3);
1386:   PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
1387:   *coneSize = DMPolytopeTypeGetConeSize(ctNew);
1388:   PetscFunctionReturn(PETSC_SUCCESS);
1389: }

1391: /* The orientation o is for the interior of the cell p */
1392: 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[])
1393: {
1394:   DM              dm;
1395:   const PetscInt  csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1396:   const PetscInt *cone;
1397:   PetscInt        c, coff = *coneoff, ooff = *orntoff;

1399:   PetscFunctionBegin;
1400:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1401:   PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
1402:   for (c = 0; c < csizeNew; ++c) {
1403:     PetscInt             ppp   = -1;                            /* Parent Parent point: Parent of point pp */
1404:     PetscInt             pp    = p;                             /* Parent point: Point in the original mesh producing new cone point */
1405:     PetscInt             po    = 0;                             /* Orientation of parent point pp in parent parent point ppp */
1406:     DMPolytopeType       pct   = ct;                            /* Parent type: Cell type for parent of new cone point */
1407:     const PetscInt      *pcone = cone;                          /* Parent cone: Cone of parent point pp */
1408:     PetscInt             pr    = -1;                            /* Replica number of pp that produces new cone point  */
1409:     const DMPolytopeType ft    = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1410:     const PetscInt       fn    = rcone[coff++];                 /* Number of cones of p that need to be taken when producing new cone point */
1411:     PetscInt             fo    = rornt[ooff++];                 /* Orientation of new cone point in pNew */
1412:     PetscInt             lc;

1414:     /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1415:     for (lc = 0; lc < fn; ++lc) {
1416:       const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po);
1417:       const PetscInt  acp  = rcone[coff++];
1418:       const PetscInt  pcp  = parr[acp * 2];
1419:       const PetscInt  pco  = parr[acp * 2 + 1];
1420:       const PetscInt *ppornt;

1422:       ppp = pp;
1423:       pp  = pcone[pcp];
1424:       PetscCall(DMPlexGetCellType(dm, pp, &pct));
1425:       // Restore the parent cone from the last iterate
1426:       if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
1427:       PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
1428:       PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
1429:       po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1430:       PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
1431:     }
1432:     if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
1433:     pr = rcone[coff++];
1434:     /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1435:     PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
1436:     PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
1437:     orntNew[c] = fo;
1438:   }
1439:   PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
1440:   *coneoff = coff;
1441:   *orntoff = ooff;
1442:   PetscFunctionReturn(PETSC_SUCCESS);
1443: }

1445: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1446: {
1447:   DM             dm;
1448:   DMPolytopeType ct;
1449:   PetscInt      *coneNew, *orntNew;
1450:   PetscInt       maxConeSize = 0, pStart, pEnd, p, pNew;

1452:   PetscFunctionBegin;
1453:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1454:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1455:   PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1456:   PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1457:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1458:   for (p = pStart; p < pEnd; ++p) {
1459:     PetscInt        coff, ooff;
1460:     DMPolytopeType *rct;
1461:     PetscInt       *rsize, *rcone, *rornt;
1462:     PetscInt        Nct, n, r;

1464:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1465:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1466:     for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1467:       const DMPolytopeType ctNew = rct[n];

1469:       for (r = 0; r < rsize[n]; ++r) {
1470:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1471:         PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
1472:         PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
1473:         PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
1474:       }
1475:     }
1476:   }
1477:   PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1478:   PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1479:   PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
1480:   PetscCall(DMPlexSymmetrize(rdm));
1481:   PetscCall(DMPlexStratify(rdm));
1482:   PetscFunctionReturn(PETSC_SUCCESS);
1483: }

1485: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1486: {
1487:   DM              dm;
1488:   DMPolytopeType  ct, qct;
1489:   DMPolytopeType *rct;
1490:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
1491:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

1493:   PetscFunctionBegin;
1495:   PetscAssertPointer(cone, 4);
1496:   PetscAssertPointer(ornt, 5);
1497:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1498:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1499:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1500:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1501:   PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1502:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1503:   for (n = 0; n < Nct; ++n) {
1504:     const DMPolytopeType ctNew    = rct[n];
1505:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1506:     PetscInt             Nr       = rsize[n], fn, c;

1508:     if (ctNew == qct) Nr = r;
1509:     for (nr = 0; nr < Nr; ++nr) {
1510:       for (c = 0; c < csizeNew; ++c) {
1511:         ++coff;             /* Cell type of new cone point */
1512:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1513:         coff += fn;
1514:         ++coff; /* Replica number of new cone point */
1515:         ++ooff; /* Orientation of new cone point */
1516:       }
1517:     }
1518:     if (ctNew == qct) break;
1519:   }
1520:   PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1521:   *cone = qcone;
1522:   *ornt = qornt;
1523:   PetscFunctionReturn(PETSC_SUCCESS);
1524: }

1526: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1527: {
1528:   DM              dm;
1529:   DMPolytopeType  ct, qct;
1530:   DMPolytopeType *rct;
1531:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
1532:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

1534:   PetscFunctionBegin;
1536:   if (cone) PetscAssertPointer(cone, 3);
1537:   if (ornt) PetscAssertPointer(ornt, 4);
1538:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1539:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1540:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1541:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1542:   PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1543:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1544:   for (n = 0; n < Nct; ++n) {
1545:     const DMPolytopeType ctNew    = rct[n];
1546:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1547:     PetscInt             Nr       = rsize[n], fn, c;

1549:     if (ctNew == qct) Nr = r;
1550:     for (nr = 0; nr < Nr; ++nr) {
1551:       for (c = 0; c < csizeNew; ++c) {
1552:         ++coff;             /* Cell type of new cone point */
1553:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1554:         coff += fn;
1555:         ++coff; /* Replica number of new cone point */
1556:         ++ooff; /* Orientation of new cone point */
1557:       }
1558:     }
1559:     if (ctNew == qct) break;
1560:   }
1561:   PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1562:   if (cone) *cone = qcone;
1563:   else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1564:   if (ornt) *ornt = qornt;
1565:   else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1566:   PetscFunctionReturn(PETSC_SUCCESS);
1567: }

1569: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1570: {
1571:   DM dm;

1573:   PetscFunctionBegin;
1575:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1576:   if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
1577:   if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
1578:   PetscFunctionReturn(PETSC_SUCCESS);
1579: }

1581: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1582: {
1583:   PetscInt ict;

1585:   PetscFunctionBegin;
1586:   PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
1587:   for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1588:     const DMPolytopeType ct = (DMPolytopeType)ict;
1589:     DMPlexTransform      reftr;
1590:     DM                   refdm, trdm;
1591:     Vec                  coordinates;
1592:     const PetscScalar   *coords;
1593:     DMPolytopeType      *rct;
1594:     PetscInt            *rsize, *rcone, *rornt;
1595:     PetscInt             Nct, n, r, pNew = 0;
1596:     PetscInt             trdim, vStart, vEnd, Nc;
1597:     const PetscInt       debug = 0;
1598:     const char          *typeName;

1600:     /* Since points are 0-dimensional, coordinates make no sense */
1601:     if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue;
1602:     PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
1603:     PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
1604:     PetscCall(DMPlexTransformSetDM(reftr, refdm));
1605:     PetscCall(DMPlexTransformGetType(tr, &typeName));
1606:     PetscCall(DMPlexTransformSetType(reftr, typeName));
1607:     PetscCall(DMPlexTransformSetUp(reftr));
1608:     PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));

1610:     PetscCall(DMGetDimension(trdm, &trdim));
1611:     PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
1612:     tr->trNv[ct] = vEnd - vStart;
1613:     PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
1614:     PetscCall(VecGetLocalSize(coordinates, &Nc));
1615:     PetscCheck(tr->trNv[ct] * trdim == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell type %s, transformed coordinate size %" PetscInt_FMT " != %" PetscInt_FMT " size of coordinate storage", DMPolytopeTypes[ct], tr->trNv[ct] * trdim, Nc);
1616:     PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
1617:     PetscCall(VecGetArrayRead(coordinates, &coords));
1618:     PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
1619:     PetscCall(VecRestoreArrayRead(coordinates, &coords));

1621:     PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
1622:     PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1623:     for (n = 0; n < Nct; ++n) {
1624:       /* Since points are 0-dimensional, coordinates make no sense */
1625:       if (rct[n] == DM_POLYTOPE_POINT) continue;
1626:       PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
1627:       for (r = 0; r < rsize[n]; ++r) {
1628:         PetscInt *closure = NULL;
1629:         PetscInt  clSize, cl, Nv = 0;

1631:         PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
1632:         PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
1633:         PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1634:         for (cl = 0; cl < clSize * 2; cl += 2) {
1635:           const PetscInt sv = closure[cl];

1637:           if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1638:         }
1639:         PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1640:         PetscCheck(Nv == DMPolytopeTypeGetNumVertices(rct[n]), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of vertices %" PetscInt_FMT " != %" PetscInt_FMT " for %s subcell %" PetscInt_FMT " from cell %s", Nv, DMPolytopeTypeGetNumVertices(rct[n]), DMPolytopeTypes[rct[n]], r, DMPolytopeTypes[ct]);
1641:       }
1642:     }
1643:     if (debug) {
1644:       DMPolytopeType *rct;
1645:       PetscInt       *rsize, *rcone, *rornt;
1646:       PetscInt        v, dE = trdim, d, off = 0;

1648:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
1649:       for (v = 0; v < tr->trNv[ct]; ++v) {
1650:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  "));
1651:         for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
1652:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1653:       }

1655:       PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1656:       for (n = 0; n < Nct; ++n) {
1657:         if (rct[n] == DM_POLYTOPE_POINT) continue;
1658:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
1659:         for (r = 0; r < rsize[n]; ++r) {
1660:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  "));
1661:           for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
1662:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1663:         }
1664:       }
1665:     }
1666:     PetscCall(DMDestroy(&refdm));
1667:     PetscCall(DMDestroy(&trdm));
1668:     PetscCall(DMPlexTransformDestroy(&reftr));
1669:   }
1670:   PetscFunctionReturn(PETSC_SUCCESS);
1671: }

1673: /*@C
1674:   DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type

1676:   Input Parameters:
1677: + tr - The `DMPlexTransform` object
1678: - ct - The cell type

1680:   Output Parameters:
1681: + Nv      - The number of transformed vertices in the closure of the reference cell of given type
1682: - trVerts - The coordinates of these vertices in the reference cell

1684:   Level: developer

1686: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
1687: @*/
1688: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1689: {
1690:   PetscFunctionBegin;
1691:   if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1692:   if (Nv) *Nv = tr->trNv[ct];
1693:   if (trVerts) *trVerts = tr->trVerts[ct];
1694:   PetscFunctionReturn(PETSC_SUCCESS);
1695: }

1697: /*@C
1698:   DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type

1700:   Input Parameters:
1701: + tr  - The `DMPlexTransform` object
1702: . ct  - The cell type
1703: . rct - The subcell type
1704: - r   - The subcell index

1706:   Output Parameter:
1707: . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()`

1709:   Level: developer

1711: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
1712: @*/
1713: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1714: {
1715:   PetscFunctionBegin;
1716:   if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1717:   PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
1718:   if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1719:   PetscFunctionReturn(PETSC_SUCCESS);
1720: }

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

1727:   PetscFunctionBeginHot;
1728:   PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
1729:   for (d = 0; d < dE; ++d) out[d] = 0.0;
1730:   for (v = 0; v < Nv; ++v)
1731:     for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1732:   for (d = 0; d < dE; ++d) out[d] /= Nv;
1733:   PetscFunctionReturn(PETSC_SUCCESS);
1734: }

1736: /*@
1737:   DMPlexTransformMapCoordinates - Calculate new coordinates for produced points

1739:   Not collective

1741:   Input Parameters:
1742: + tr  - The `DMPlexTransform`
1743: . pct - The cell type of the parent, from whom the new cell is being produced
1744: . ct  - The type being produced
1745: . p   - The original point
1746: . r   - The replica number requested for the produced cell type
1747: . Nv  - Number of vertices in the closure of the parent cell
1748: . dE  - Spatial dimension
1749: - in  - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell

1751:   Output Parameter:
1752: . out - The coordinates of the new vertices

1754:   Level: intermediate

1756: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
1757: @*/
1758: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1759: {
1760:   PetscFunctionBeginHot;
1761:   if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1762:   PetscFunctionReturn(PETSC_SUCCESS);
1763: }

1765: /*
1766:   DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label

1768:   Not Collective

1770:   Input Parameters:
1771: + tr    - The `DMPlexTransform`
1772: . label - The label in the transformed mesh
1773: . pp    - The parent point in the original mesh
1774: . pct   - The cell type of the parent point
1775: . p     - The point in the transformed mesh
1776: . ct    - The cell type of the point
1777: . r     - The replica number of the point
1778: - val   - The label value of the parent point

1780:   Level: developer

1782: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
1783: */
1784: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
1785: {
1786:   PetscFunctionBeginHot;
1787:   if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
1788:   PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
1789:   PetscFunctionReturn(PETSC_SUCCESS);
1790: }

1792: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1793: {
1794:   DM              dm;
1795:   IS              valueIS;
1796:   const PetscInt *values;
1797:   PetscInt        defVal, Nv, val;

1799:   PetscFunctionBegin;
1800:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1801:   PetscCall(DMLabelGetDefaultValue(label, &defVal));
1802:   PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
1803:   PetscCall(DMLabelGetValueIS(label, &valueIS));
1804:   PetscCall(ISGetLocalSize(valueIS, &Nv));
1805:   PetscCall(ISGetIndices(valueIS, &values));
1806:   for (val = 0; val < Nv; ++val) {
1807:     IS              pointIS;
1808:     const PetscInt *points;
1809:     PetscInt        numPoints, p;

1811:     /* Ensure refined label is created with same number of strata as
1812:      * original (even if no entries here). */
1813:     PetscCall(DMLabelAddStratum(labelNew, values[val]));
1814:     PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
1815:     PetscCall(ISGetLocalSize(pointIS, &numPoints));
1816:     PetscCall(ISGetIndices(pointIS, &points));
1817:     for (p = 0; p < numPoints; ++p) {
1818:       const PetscInt  point = points[p];
1819:       DMPolytopeType  ct;
1820:       DMPolytopeType *rct;
1821:       PetscInt       *rsize, *rcone, *rornt;
1822:       PetscInt        Nct, n, r, pNew = 0;

1824:       PetscCall(DMPlexGetCellType(dm, point, &ct));
1825:       PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1826:       for (n = 0; n < Nct; ++n) {
1827:         for (r = 0; r < rsize[n]; ++r) {
1828:           PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
1829:           PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
1830:         }
1831:       }
1832:     }
1833:     PetscCall(ISRestoreIndices(pointIS, &points));
1834:     PetscCall(ISDestroy(&pointIS));
1835:   }
1836:   PetscCall(ISRestoreIndices(valueIS, &values));
1837:   PetscCall(ISDestroy(&valueIS));
1838:   PetscFunctionReturn(PETSC_SUCCESS);
1839: }

1841: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1842: {
1843:   DM       dm;
1844:   PetscInt numLabels, l;

1846:   PetscFunctionBegin;
1847:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1848:   PetscCall(DMGetNumLabels(dm, &numLabels));
1849:   for (l = 0; l < numLabels; ++l) {
1850:     DMLabel     label, labelNew;
1851:     const char *lname;
1852:     PetscBool   isDepth, isCellType;

1854:     PetscCall(DMGetLabelName(dm, l, &lname));
1855:     PetscCall(PetscStrcmp(lname, "depth", &isDepth));
1856:     if (isDepth) continue;
1857:     PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
1858:     if (isCellType) continue;
1859:     PetscCall(DMCreateLabel(rdm, lname));
1860:     PetscCall(DMGetLabel(dm, lname, &label));
1861:     PetscCall(DMGetLabel(rdm, lname, &labelNew));
1862:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
1863:   }
1864:   PetscFunctionReturn(PETSC_SUCCESS);
1865: }

1867: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1868: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1869: {
1870:   DM       dm;
1871:   PetscInt Nf, f, Nds, s;

1873:   PetscFunctionBegin;
1874:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1875:   PetscCall(DMGetNumFields(dm, &Nf));
1876:   for (f = 0; f < Nf; ++f) {
1877:     DMLabel     label, labelNew;
1878:     PetscObject obj;
1879:     const char *lname;

1881:     PetscCall(DMGetField(rdm, f, &label, &obj));
1882:     if (!label) continue;
1883:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1884:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1885:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
1886:     PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
1887:     PetscCall(DMLabelDestroy(&labelNew));
1888:   }
1889:   PetscCall(DMGetNumDS(dm, &Nds));
1890:   for (s = 0; s < Nds; ++s) {
1891:     DMLabel     label, labelNew;
1892:     const char *lname;

1894:     PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
1895:     if (!label) continue;
1896:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1897:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1898:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
1899:     PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
1900:     PetscCall(DMLabelDestroy(&labelNew));
1901:   }
1902:   PetscFunctionReturn(PETSC_SUCCESS);
1903: }

1905: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1906: {
1907:   DM                 dm;
1908:   PetscSF            sf, sfNew;
1909:   PetscInt           numRoots, numLeaves, numLeavesNew = 0, l, m;
1910:   const PetscInt    *localPoints;
1911:   const PetscSFNode *remotePoints;
1912:   PetscInt          *localPointsNew;
1913:   PetscSFNode       *remotePointsNew;
1914:   PetscInt           pStartNew, pEndNew, pNew;
1915:   /* Brute force algorithm */
1916:   PetscSF         rsf;
1917:   PetscSection    s;
1918:   const PetscInt *rootdegree;
1919:   PetscInt       *rootPointsNew, *remoteOffsets;
1920:   PetscInt        numPointsNew, pStart, pEnd, p;

1922:   PetscFunctionBegin;
1923:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1924:   PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
1925:   PetscCall(DMGetPointSF(dm, &sf));
1926:   PetscCall(DMGetPointSF(rdm, &sfNew));
1927:   /* Calculate size of new SF */
1928:   PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
1929:   if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
1930:   for (l = 0; l < numLeaves; ++l) {
1931:     const PetscInt  p = localPoints[l];
1932:     DMPolytopeType  ct;
1933:     DMPolytopeType *rct;
1934:     PetscInt       *rsize, *rcone, *rornt;
1935:     PetscInt        Nct, n;

1937:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1938:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1939:     for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
1940:   }
1941:   /* Send new root point numbers
1942:        It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1943:   */
1944:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1945:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
1946:   PetscCall(PetscSectionSetChart(s, pStart, pEnd));
1947:   for (p = pStart; p < pEnd; ++p) {
1948:     DMPolytopeType  ct;
1949:     DMPolytopeType *rct;
1950:     PetscInt       *rsize, *rcone, *rornt;
1951:     PetscInt        Nct, n;

1953:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1954:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1955:     for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
1956:   }
1957:   PetscCall(PetscSectionSetUp(s));
1958:   PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
1959:   PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
1960:   PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
1961:   PetscCall(PetscFree(remoteOffsets));
1962:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1963:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1964:   PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
1965:   for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1966:   for (p = pStart; p < pEnd; ++p) {
1967:     DMPolytopeType  ct;
1968:     DMPolytopeType *rct;
1969:     PetscInt       *rsize, *rcone, *rornt;
1970:     PetscInt        Nct, n, r, off;

1972:     if (!rootdegree[p - pStart]) continue;
1973:     PetscCall(PetscSectionGetOffset(s, p, &off));
1974:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1975:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1976:     for (n = 0, m = 0; n < Nct; ++n) {
1977:       for (r = 0; r < rsize[n]; ++r, ++m) {
1978:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1979:         rootPointsNew[off + m] = pNew;
1980:       }
1981:     }
1982:   }
1983:   PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1984:   PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1985:   PetscCall(PetscSFDestroy(&rsf));
1986:   PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
1987:   PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
1988:   for (l = 0, m = 0; l < numLeaves; ++l) {
1989:     const PetscInt  p = localPoints[l];
1990:     DMPolytopeType  ct;
1991:     DMPolytopeType *rct;
1992:     PetscInt       *rsize, *rcone, *rornt;
1993:     PetscInt        Nct, n, r, q, off;

1995:     PetscCall(PetscSectionGetOffset(s, p, &off));
1996:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1997:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1998:     for (n = 0, q = 0; n < Nct; ++n) {
1999:       for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
2000:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2001:         localPointsNew[m]        = pNew;
2002:         remotePointsNew[m].index = rootPointsNew[off + q];
2003:         remotePointsNew[m].rank  = remotePoints[l].rank;
2004:       }
2005:     }
2006:   }
2007:   PetscCall(PetscSectionDestroy(&s));
2008:   PetscCall(PetscFree(rootPointsNew));
2009:   /* SF needs sorted leaves to correctly calculate Gather */
2010:   {
2011:     PetscSFNode *rp, *rtmp;
2012:     PetscInt    *lp, *idx, *ltmp, i;

2014:     PetscCall(PetscMalloc1(numLeavesNew, &idx));
2015:     PetscCall(PetscMalloc1(numLeavesNew, &lp));
2016:     PetscCall(PetscMalloc1(numLeavesNew, &rp));
2017:     for (i = 0; i < numLeavesNew; ++i) {
2018:       PetscCheck(!(localPointsNew[i] < pStartNew) && !(localPointsNew[i] >= pEndNew), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local SF point %" PetscInt_FMT " (%" PetscInt_FMT ") not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", localPointsNew[i], i, pStartNew, pEndNew);
2019:       idx[i] = i;
2020:     }
2021:     PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
2022:     for (i = 0; i < numLeavesNew; ++i) {
2023:       lp[i] = localPointsNew[idx[i]];
2024:       rp[i] = remotePointsNew[idx[i]];
2025:     }
2026:     ltmp            = localPointsNew;
2027:     localPointsNew  = lp;
2028:     rtmp            = remotePointsNew;
2029:     remotePointsNew = rp;
2030:     PetscCall(PetscFree(idx));
2031:     PetscCall(PetscFree(ltmp));
2032:     PetscCall(PetscFree(rtmp));
2033:   }
2034:   PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
2035:   PetscFunctionReturn(PETSC_SUCCESS);
2036: }

2038: /*
2039:   DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`.

2041:   Not Collective

2043:   Input Parameters:
2044: + tr  - The `DMPlexTransform`
2045: . ct  - The type of the parent cell
2046: . rct - The type of the produced cell
2047: . r   - The index of the produced cell
2048: - x   - The localized coordinates for the parent cell

2050:   Output Parameter:
2051: . xr  - The localized coordinates for the produced cell

2053:   Level: developer

2055: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
2056: */
2057: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
2058: {
2059:   PetscFE  fe = NULL;
2060:   PetscInt cdim, v, *subcellV;

2062:   PetscFunctionBegin;
2063:   PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
2064:   PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
2065:   PetscCall(PetscFEGetNumComponents(fe, &cdim));
2066:   for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
2067:   PetscFunctionReturn(PETSC_SUCCESS);
2068: }

2070: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
2071: {
2072:   DM                 dm, cdm, cdmCell, cdmNew, cdmCellNew;
2073:   PetscSection       coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
2074:   Vec                coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
2075:   const PetscScalar *coords;
2076:   PetscScalar       *coordsNew;
2077:   const PetscReal   *maxCell, *Lstart, *L;
2078:   PetscBool          localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
2079:   PetscInt           dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;

2081:   PetscFunctionBegin;
2082:   // Need to clear the DMField for coordinates
2083:   PetscCall(DMSetCoordinateField(rdm, NULL));
2084:   PetscCall(DMPlexTransformGetDM(tr, &dm));
2085:   PetscCall(DMGetCoordinateDM(dm, &cdm));
2086:   PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
2087:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
2088:   PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
2089:   if (localized) {
2090:     /* Localize coordinates of new vertices */
2091:     localizeVertices = PETSC_TRUE;
2092:     /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
2093:     if (!maxCell) localizeCells = PETSC_TRUE;
2094:   }
2095:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
2096:   PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
2097:   if (maxCell) {
2098:     PetscReal maxCellNew[3];

2100:     for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
2101:     PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L));
2102:   }
2103:   PetscCall(DMGetCoordinateDim(rdm, &dE));
2104:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
2105:   PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
2106:   PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
2107:   PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
2108:   PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
2109:   /* Localization should be inherited */
2110:   /*   Stefano calculates parent cells for each new cell for localization */
2111:   /*   Localized cells need coordinates of closure */
2112:   for (v = vStartNew; v < vEndNew; ++v) {
2113:     PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
2114:     PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
2115:   }
2116:   PetscCall(PetscSectionSetUp(coordSectionNew));
2117:   PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));

2119:   if (localizeCells) {
2120:     PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
2121:     PetscCall(DMClone(cdmNew, &cdmCellNew));
2122:     PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
2123:     PetscCall(DMDestroy(&cdmCellNew));

2125:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
2126:     PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
2127:     PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
2128:     PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
2129:     PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));

2131:     PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
2132:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2133:     for (c = cStart; c < cEnd; ++c) {
2134:       PetscInt dof;

2136:       PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
2137:       if (dof) {
2138:         DMPolytopeType  ct;
2139:         DMPolytopeType *rct;
2140:         PetscInt       *rsize, *rcone, *rornt;
2141:         PetscInt        dim, cNew, Nct, n, r;

2143:         PetscCall(DMPlexGetCellType(dm, c, &ct));
2144:         dim = DMPolytopeTypeGetDim(ct);
2145:         PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2146:         /* This allows for different cell types */
2147:         for (n = 0; n < Nct; ++n) {
2148:           if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2149:           for (r = 0; r < rsize[n]; ++r) {
2150:             PetscInt *closure = NULL;
2151:             PetscInt  clSize, cl, Nv = 0;

2153:             PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
2154:             PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2155:             for (cl = 0; cl < clSize * 2; cl += 2) {
2156:               if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
2157:             }
2158:             PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2159:             PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
2160:             PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
2161:           }
2162:         }
2163:       }
2164:     }
2165:     PetscCall(PetscSectionSetUp(coordSectionCellNew));
2166:     PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
2167:   }
2168:   PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
2169:   {
2170:     VecType     vtype;
2171:     PetscInt    coordSizeNew, bs;
2172:     const char *name;

2174:     PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
2175:     PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
2176:     PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
2177:     PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
2178:     PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
2179:     PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
2180:     PetscCall(VecGetBlockSize(coordsLocal, &bs));
2181:     PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
2182:     PetscCall(VecGetType(coordsLocal, &vtype));
2183:     PetscCall(VecSetType(coordsLocalNew, vtype));
2184:   }
2185:   PetscCall(VecGetArrayRead(coordsLocal, &coords));
2186:   PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
2187:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2188:   /* First set coordinates for vertices */
2189:   for (p = pStart; p < pEnd; ++p) {
2190:     DMPolytopeType  ct;
2191:     DMPolytopeType *rct;
2192:     PetscInt       *rsize, *rcone, *rornt;
2193:     PetscInt        Nct, n, r;
2194:     PetscBool       hasVertex = PETSC_FALSE;

2196:     PetscCall(DMPlexGetCellType(dm, p, &ct));
2197:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2198:     for (n = 0; n < Nct; ++n) {
2199:       if (rct[n] == DM_POLYTOPE_POINT) {
2200:         hasVertex = PETSC_TRUE;
2201:         break;
2202:       }
2203:     }
2204:     if (hasVertex) {
2205:       const PetscScalar *icoords = NULL;
2206:       const PetscScalar *array   = NULL;
2207:       PetscScalar       *pcoords = NULL;
2208:       PetscBool          isDG;
2209:       PetscInt           Nc, Nv, v, d;

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

2213:       icoords = pcoords;
2214:       Nv      = Nc / dEo;
2215:       if (ct != DM_POLYTOPE_POINT) {
2216:         if (localizeVertices && maxCell) {
2217:           PetscScalar anchor[3];

2219:           for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
2220:           for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
2221:         }
2222:       }
2223:       for (n = 0; n < Nct; ++n) {
2224:         if (rct[n] != DM_POLYTOPE_POINT) continue;
2225:         for (r = 0; r < rsize[n]; ++r) {
2226:           PetscScalar vcoords[3];
2227:           PetscInt    vNew, off;

2229:           PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
2230:           PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
2231:           PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
2232:           PetscCall(DMSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
2233:         }
2234:       }
2235:       PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2236:     }
2237:   }
2238:   PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
2239:   PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
2240:   PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
2241:   PetscCall(VecDestroy(&coordsLocalNew));
2242:   PetscCall(PetscSectionDestroy(&coordSectionNew));
2243:   /* Then set coordinates for cells by localizing */
2244:   if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
2245:   else {
2246:     VecType     vtype;
2247:     PetscInt    coordSizeNew, bs;
2248:     const char *name;

2250:     PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
2251:     PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
2252:     PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
2253:     PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
2254:     PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
2255:     PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
2256:     PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
2257:     PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
2258:     PetscCall(VecGetType(coordsLocalCell, &vtype));
2259:     PetscCall(VecSetType(coordsLocalCellNew, vtype));
2260:     PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
2261:     PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));

2263:     for (p = pStart; p < pEnd; ++p) {
2264:       DMPolytopeType  ct;
2265:       DMPolytopeType *rct;
2266:       PetscInt       *rsize, *rcone, *rornt;
2267:       PetscInt        dof = 0, Nct, n, r;

2269:       PetscCall(DMPlexGetCellType(dm, p, &ct));
2270:       PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2271:       if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
2272:       if (dof) {
2273:         const PetscScalar *pcoords;

2275:         PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords));
2276:         for (n = 0; n < Nct; ++n) {
2277:           const PetscInt Nr = rsize[n];

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

2283:             /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2284:                DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2285:                cell to the ones it produces. */
2286:             PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2287:             PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
2288:             PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
2289:           }
2290:         }
2291:       }
2292:     }
2293:     PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
2294:     PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
2295:     PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
2296:     PetscCall(VecDestroy(&coordsLocalCellNew));
2297:     PetscCall(PetscSectionDestroy(&coordSectionCellNew));
2298:   }
2299:   PetscFunctionReturn(PETSC_SUCCESS);
2300: }

2302: /*@
2303:   DMPlexTransformApply - Execute the transformation, producing another `DM`

2305:   Collective

2307:   Input Parameters:
2308: + tr - The `DMPlexTransform` object
2309: - dm - The original `DM`

2311:   Output Parameter:
2312: . tdm - The transformed `DM`

2314:   Level: intermediate

2316:   Options Database Keys:
2317: + -dm_plex_transform_label_match_strata      - Only label points of the same stratum as the producing point
2318: . -dm_plex_transform_label_replica_inc <num> - Increment for the label value to be multiplied by the replica number
2319: - -dm_plex_transform_active <name>           - Name for active mesh label

2321: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()`
2322: @*/
2323: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
2324: {
2325:   DM                     rdm;
2326:   DMPlexInterpolatedFlag interp;
2327:   PetscInt               pStart, pEnd;

2329:   PetscFunctionBegin;
2332:   PetscAssertPointer(tdm, 3);
2333:   PetscCall(PetscLogEventBegin(DMPLEX_Transform, tr, dm, 0, 0));
2334:   PetscCall(DMPlexTransformSetDM(tr, dm));

2336:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
2337:   PetscCall(DMSetType(rdm, DMPLEX));
2338:   PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
2339:   /* Calculate number of new points of each depth */
2340:   PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
2341:   PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
2342:   /* Step 1: Set chart */
2343:   PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
2344:   PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
2345:   /* Step 2: Set cone/support sizes (automatically stratifies) */
2346:   PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
2347:   /* Step 3: Setup refined DM */
2348:   PetscCall(DMSetUp(rdm));
2349:   /* Step 4: Set cones and supports (automatically symmetrizes) */
2350:   PetscCall(DMPlexTransformSetCones(tr, rdm));
2351:   /* Step 5: Create pointSF */
2352:   PetscCall(DMPlexTransformCreateSF(tr, rdm));
2353:   /* Step 6: Create labels */
2354:   PetscCall(DMPlexTransformCreateLabels(tr, rdm));
2355:   /* Step 7: Set coordinates */
2356:   PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
2357:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm));
2358:   // If the original DM was configured from options, the transformed DM should be as well
2359:   rdm->setfromoptionscalled = dm->setfromoptionscalled;
2360:   PetscCall(PetscLogEventEnd(DMPLEX_Transform, tr, dm, 0, 0));
2361:   *tdm = rdm;
2362:   PetscFunctionReturn(PETSC_SUCCESS);
2363: }

2365: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2366: {
2367:   DMPlexTransform tr;
2368:   DM              cdm, rcdm;
2369:   const char     *prefix;

2371:   PetscFunctionBegin;
2372:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
2373:   PetscCall(PetscObjectSetName((PetscObject)tr, "Adapt Label Transform"));
2374:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2375:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
2376:   PetscCall(DMPlexTransformSetDM(tr, dm));
2377:   PetscCall(DMPlexTransformSetFromOptions(tr));
2378:   PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
2379:   PetscCall(DMPlexTransformSetUp(tr));
2380:   PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
2381:   PetscCall(DMPlexTransformApply(tr, dm, rdm));
2382:   PetscCall(DMCopyDisc(dm, *rdm));
2383:   PetscCall(DMGetCoordinateDM(dm, &cdm));
2384:   PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
2385:   PetscCall(DMCopyDisc(cdm, rcdm));
2386:   PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
2387:   PetscCall(DMCopyDisc(dm, *rdm));
2388:   PetscCall(DMPlexTransformDestroy(&tr));
2389:   ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2390:   PetscFunctionReturn(PETSC_SUCCESS);
2391: }