Actual source code: plextransform.c

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

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

  5: PetscClassId DMPLEXTRANSFORM_CLASSID;

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

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

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

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

 50:   Not Collective

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

 56:   Example Usage:
 57: .vb
 58:   DMPlexTransformRegister("my_transform", MyTransformCreate);
 59: .ve

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

 71:   Level: advanced

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

 76: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()`
 77: @*/
 78: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
 79: {
 80:   PetscFunctionBegin;
 81:   PetscCall(DMInitializePackage());
 82:   PetscCall(PetscFunctionListAdd(&DMPlexTransformList, name, create_func));
 83:   PetscFunctionReturn(PETSC_SUCCESS);
 84: }

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

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

 98:   Not Collective

100:   Level: advanced

102: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
103: @*/
104: PetscErrorCode DMPlexTransformRegisterAll(void)
105: {
106:   PetscFunctionBegin;
107:   if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
108:   DMPlexTransformRegisterAllCalled = PETSC_TRUE;

110:   PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter));
111:   PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular));
112:   PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox));
113:   PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld));
114:   PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL));
115:   PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR));
116:   PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D));
117:   PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDE, DMPlexTransformCreate_Extrude));
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

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

124:   Level: developer

126: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()`
127: @*/
128: PetscErrorCode DMPlexTransformRegisterDestroy(void)
129: {
130:   PetscFunctionBegin;
131:   PetscCall(PetscFunctionListDestroy(&DMPlexTransformList));
132:   DMPlexTransformRegisterAllCalled = PETSC_FALSE;
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

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

139:   Collective

141:   Input Parameter:
142: . comm - The communicator for the transform object

144:   Output Parameter:
145: . tr - The transform object

147:   Level: beginner

149: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
150: @*/
151: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
152: {
153:   DMPlexTransform t;

155:   PetscFunctionBegin;
156:   PetscAssertPointer(tr, 2);
157:   *tr = NULL;
158:   PetscCall(DMInitializePackage());

160:   PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView));
161:   t->setupcalled = PETSC_FALSE;
162:   PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom));
163:   *tr = t;
164:   PetscFunctionReturn(PETSC_SUCCESS);
165: }

167: /*@C
168:   DMPlexTransformSetType - Sets the particular implementation for a transform.

170:   Collective

172:   Input Parameters:
173: + tr     - The transform
174: - method - The name of the transform type

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

179:   Level: intermediate

181: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
182: @*/
183: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
184: {
185:   PetscErrorCode (*r)(DMPlexTransform);
186:   PetscBool match;

188:   PetscFunctionBegin;
190:   PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match));
191:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

193:   PetscCall(DMPlexTransformRegisterAll());
194:   PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r));
195:   PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);

197:   PetscTryTypeMethod(tr, destroy);
198:   PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops)));
199:   PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method));
200:   PetscCall((*r)(tr));
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

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

207:   Not Collective

209:   Input Parameter:
210: . tr - The `DMPlexTransform`

212:   Output Parameter:
213: . type - The `DMPlexTransformType` name

215:   Level: intermediate

217: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
218: @*/
219: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
220: {
221:   PetscFunctionBegin;
223:   PetscAssertPointer(type, 2);
224:   PetscCall(DMPlexTransformRegisterAll());
225:   *type = ((PetscObject)tr)->type_name;
226:   PetscFunctionReturn(PETSC_SUCCESS);
227: }

229: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
230: {
231:   PetscViewerFormat format;

233:   PetscFunctionBegin;
234:   PetscCall(PetscViewerGetFormat(v, &format));
235:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
236:     const PetscInt *trTypes = NULL;
237:     IS              trIS;
238:     PetscInt        cols = 8;
239:     PetscInt        Nrt  = 8, f, g;

241:     if (tr->trType) PetscCall(DMLabelView(tr->trType, v));
242:     PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n"));
243:     for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
244:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
245:     for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]));
246:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
247:     PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n"));
248:     for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
249:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
250:     for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]));
251:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));

253:     if (tr->trType) {
254:       PetscCall(DMLabelGetNumValues(tr->trType, &Nrt));
255:       PetscCall(DMLabelGetValueIS(tr->trType, &trIS));
256:       PetscCall(ISGetIndices(trIS, &trTypes));
257:     }
258:     PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n"));
259:     PetscCall(PetscViewerASCIIPrintf(v, "     "));
260:     for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
261:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
262:     for (f = 0; f < Nrt; ++f) {
263:       PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT "  |", trTypes ? trTypes[f] : f));
264:       for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]));
265:       PetscCall(PetscViewerASCIIPrintf(v, " |\n"));
266:     }
267:     if (tr->trType) {
268:       PetscCall(ISRestoreIndices(trIS, &trTypes));
269:       PetscCall(ISDestroy(&trIS));
270:     }
271:   }
272:   PetscFunctionReturn(PETSC_SUCCESS);
273: }

275: /*@C
276:   DMPlexTransformView - Views a `DMPlexTransform`

278:   Collective

280:   Input Parameters:
281: + tr - the `DMPlexTransform` object to view
282: - v  - the viewer

284:   Level: beginner

286: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
287: @*/
288: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
289: {
290:   PetscBool isascii;

292:   PetscFunctionBegin;
294:   if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v));
296:   PetscCheckSameComm(tr, 1, v, 2);
297:   PetscCall(PetscViewerCheckWritable(v));
298:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v));
299:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
300:   if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v));
301:   PetscTryTypeMethod(tr, view, v);
302:   PetscFunctionReturn(PETSC_SUCCESS);
303: }

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

308:   Collective

310:   Input Parameter:
311: . tr - the `DMPlexTransform` object to set options for

313:   Options Database Keys:
314: + -dm_plex_transform_type                    - Set the transform type, e.g. refine_regular
315: . -dm_plex_transform_label_match_strata      - Only label points of the same stratum as the producing point
316: - -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

318:   Level: intermediate

320: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
321: @*/
322: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
323: {
324:   char        typeName[1024];
325:   const char *defName = DMPLEXREFINEREGULAR;
326:   PetscBool   flg;

328:   PetscFunctionBegin;
330:   PetscObjectOptionsBegin((PetscObject)tr);
331:   PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg));
332:   if (flg) PetscCall(DMPlexTransformSetType(tr, typeName));
333:   else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName));
334:   PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &tr->labelMatchStrata, NULL));
335:   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));
336:   PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
337:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
338:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject));
339:   PetscOptionsEnd();
340:   PetscFunctionReturn(PETSC_SUCCESS);
341: }

343: /*@C
344:   DMPlexTransformDestroy - Destroys a `DMPlexTransform`

346:   Collective

348:   Input Parameter:
349: . tr - the transform object to destroy

351:   Level: beginner

353: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
354: @*/
355: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
356: {
357:   PetscInt c;

359:   PetscFunctionBegin;
360:   if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
362:   if (--((PetscObject)(*tr))->refct > 0) {
363:     *tr = NULL;
364:     PetscFunctionReturn(PETSC_SUCCESS);
365:   }

367:   if ((*tr)->ops->destroy) PetscCall((*(*tr)->ops->destroy)(*tr));
368:   PetscCall(DMDestroy(&(*tr)->dm));
369:   PetscCall(DMLabelDestroy(&(*tr)->active));
370:   PetscCall(DMLabelDestroy(&(*tr)->trType));
371:   PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld));
372:   PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew));
373:   PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew));
374:   PetscCall(PetscFree((*tr)->offset));
375:   PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd));
376:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
377:     PetscCall(PetscFEDestroy(&(*tr)->coordFE[c]));
378:     PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c]));
379:   }
380:   if ((*tr)->trVerts) {
381:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
382:       DMPolytopeType *rct;
383:       PetscInt       *rsize, *rcone, *rornt, Nct, n, r;

385:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0) {
386:         PetscCall(DMPlexTransformCellTransform((*tr), (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
387:         for (n = 0; n < Nct; ++n) {
388:           if (rct[n] == DM_POLYTOPE_POINT) continue;
389:           for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r]));
390:           PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]]));
391:         }
392:       }
393:       PetscCall(PetscFree((*tr)->trSubVerts[c]));
394:       PetscCall(PetscFree((*tr)->trVerts[c]));
395:     }
396:   }
397:   PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts));
398:   PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom));
399:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
400:   PetscCall(PetscHeaderDestroy(tr));
401:   PetscFunctionReturn(PETSC_SUCCESS);
402: }

404: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
405: {
406:   DMLabel  trType = tr->trType;
407:   PetscInt c, cN, *off;

409:   PetscFunctionBegin;
410:   if (trType) {
411:     DM              dm;
412:     IS              rtIS;
413:     const PetscInt *reftypes;
414:     PetscInt        Nrt, r;

416:     PetscCall(DMPlexTransformGetDM(tr, &dm));
417:     PetscCall(DMLabelGetNumValues(trType, &Nrt));
418:     PetscCall(DMLabelGetValueIS(trType, &rtIS));
419:     PetscCall(ISGetIndices(rtIS, &reftypes));
420:     PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off));
421:     for (r = 0; r < Nrt; ++r) {
422:       const PetscInt  rt = reftypes[r];
423:       IS              rtIS;
424:       const PetscInt *points;
425:       DMPolytopeType  ct;
426:       PetscInt        np, p;

428:       PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
429:       PetscCall(ISGetLocalSize(rtIS, &np));
430:       PetscCall(ISGetIndices(rtIS, &points));
431:       if (!np) continue;
432:       p = points[0];
433:       PetscCall(ISRestoreIndices(rtIS, &points));
434:       PetscCall(ISDestroy(&rtIS));
435:       PetscCall(DMPlexGetCellType(dm, p, &ct));
436:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
437:         const DMPolytopeType ctNew = (DMPolytopeType)cN;
438:         DMPolytopeType      *rct;
439:         PetscInt            *rsize, *cone, *ornt;
440:         PetscInt             Nct, n, s;

442:         if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
443:           off[r * DM_NUM_POLYTOPES + ctNew] = -1;
444:           break;
445:         }
446:         off[r * DM_NUM_POLYTOPES + ctNew] = 0;
447:         for (s = 0; s <= r; ++s) {
448:           const PetscInt st = reftypes[s];
449:           DMPolytopeType sct;
450:           PetscInt       q, qrt;

452:           PetscCall(DMLabelGetStratumIS(trType, st, &rtIS));
453:           PetscCall(ISGetLocalSize(rtIS, &np));
454:           PetscCall(ISGetIndices(rtIS, &points));
455:           if (!np) continue;
456:           q = points[0];
457:           PetscCall(ISRestoreIndices(rtIS, &points));
458:           PetscCall(ISDestroy(&rtIS));
459:           PetscCall(DMPlexGetCellType(dm, q, &sct));
460:           PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt));
461:           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);
462:           if (st == rt) {
463:             for (n = 0; n < Nct; ++n)
464:               if (rct[n] == ctNew) break;
465:             if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
466:             break;
467:           }
468:           for (n = 0; n < Nct; ++n) {
469:             if (rct[n] == ctNew) {
470:               PetscInt sn;

472:               PetscCall(DMLabelGetStratumSize(trType, st, &sn));
473:               off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
474:             }
475:           }
476:         }
477:       }
478:     }
479:     PetscCall(ISRestoreIndices(rtIS, &reftypes));
480:     PetscCall(ISDestroy(&rtIS));
481:   } else {
482:     PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off));
483:     for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
484:       const DMPolytopeType ct = (DMPolytopeType)c;
485:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
486:         const DMPolytopeType ctNew = (DMPolytopeType)cN;
487:         DMPolytopeType      *rct;
488:         PetscInt            *rsize, *cone, *ornt;
489:         PetscInt             Nct, n, i;

491:         if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
492:           off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
493:           continue;
494:         }
495:         off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
496:         for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
497:           const DMPolytopeType ict  = (DMPolytopeType)ctOrderOld[i];
498:           const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];

500:           PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt));
501:           if (ict == ct) {
502:             for (n = 0; n < Nct; ++n)
503:               if (rct[n] == ctNew) break;
504:             if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
505:             break;
506:           }
507:           for (n = 0; n < Nct; ++n)
508:             if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
509:         }
510:       }
511:     }
512:   }
513:   *offset = off;
514:   PetscFunctionReturn(PETSC_SUCCESS);
515: }

517: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
518: {
519:   DM             dm;
520:   DMPolytopeType ctCell;
521:   PetscInt       pStart, pEnd, p, c, celldim = 0;

523:   PetscFunctionBegin;
525:   if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
526:   PetscTryTypeMethod(tr, setup);
527:   PetscCall(DMPlexTransformGetDM(tr, &dm));
528:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
529:   if (pEnd > pStart) {
530:     PetscCall(DMPlexGetCellType(dm, 0, &ctCell));
531:   } else {
532:     PetscInt dim;

534:     PetscCall(DMGetDimension(dm, &dim));
535:     switch (dim) {
536:     case 0:
537:       ctCell = DM_POLYTOPE_POINT;
538:       break;
539:     case 1:
540:       ctCell = DM_POLYTOPE_SEGMENT;
541:       break;
542:     case 2:
543:       ctCell = DM_POLYTOPE_TRIANGLE;
544:       break;
545:     case 3:
546:       ctCell = DM_POLYTOPE_TETRAHEDRON;
547:       break;
548:     default:
549:       ctCell = DM_POLYTOPE_UNKNOWN;
550:     }
551:   }
552:   PetscCall(DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
553:   for (p = pStart; p < pEnd; ++p) {
554:     DMPolytopeType  ct;
555:     DMPolytopeType *rct;
556:     PetscInt       *rsize, *cone, *ornt;
557:     PetscInt        Nct, n;

559:     PetscCall(DMPlexGetCellType(dm, p, &ct));
560:     PetscCheck(ct != DM_POLYTOPE_UNKNOWN, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p);
561:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
562:     for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
563:   }
564:   PetscCall(DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
565:   /* Construct sizes and offsets for each cell type */
566:   if (!tr->ctStart) {
567:     PetscInt *ctS, *ctSN, *ctC, *ctCN;

569:     PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
570:     PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
571:     for (p = pStart; p < pEnd; ++p) {
572:       DMPolytopeType  ct;
573:       DMPolytopeType *rct;
574:       PetscInt       *rsize, *cone, *ornt;
575:       PetscInt        Nct, n;

577:       PetscCall(DMPlexGetCellType(dm, p, &ct));
578:       PetscCheck(ct != DM_POLYTOPE_UNKNOWN, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p);
579:       ++ctC[ct];
580:       PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
581:       for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
582:     }
583:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
584:       const PetscInt cto  = tr->ctOrderOld[c];
585:       const PetscInt cton = tr->ctOrderOld[c + 1];
586:       const PetscInt ctn  = tr->ctOrderNew[c];
587:       const PetscInt ctnn = tr->ctOrderNew[c + 1];

589:       ctS[cton]  = ctS[cto] + ctC[cto];
590:       ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
591:     }
592:     PetscCall(PetscFree2(ctC, ctCN));
593:     tr->ctStart    = ctS;
594:     tr->ctStartNew = ctSN;
595:   }
596:   PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
597:   // Compute depth information
598:   tr->depth = -1;
599:   for (c = 0; c < DM_NUM_POLYTOPES; ++c)
600:     if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
601:   PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
602:   for (PetscInt d = 0; d <= tr->depth; ++d) {
603:     tr->depthStart[d] = PETSC_MAX_INT;
604:     tr->depthEnd[d]   = -1;
605:   }
606:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
607:     const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);

609:     if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
610:     tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
611:     tr->depthEnd[dep]   = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
612:   }
613:   tr->setupcalled = PETSC_TRUE;
614:   PetscFunctionReturn(PETSC_SUCCESS);
615: }

617: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
618: {
619:   PetscFunctionBegin;
621:   PetscAssertPointer(dm, 2);
622:   *dm = tr->dm;
623:   PetscFunctionReturn(PETSC_SUCCESS);
624: }

626: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
627: {
628:   PetscFunctionBegin;
631:   PetscCall(PetscObjectReference((PetscObject)dm));
632:   PetscCall(DMDestroy(&tr->dm));
633:   tr->dm = dm;
634:   PetscFunctionReturn(PETSC_SUCCESS);
635: }

637: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
638: {
639:   PetscFunctionBegin;
641:   PetscAssertPointer(active, 2);
642:   *active = tr->active;
643:   PetscFunctionReturn(PETSC_SUCCESS);
644: }

646: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
647: {
648:   PetscFunctionBegin;
651:   PetscCall(PetscObjectReference((PetscObject)active));
652:   PetscCall(DMLabelDestroy(&tr->active));
653:   tr->active = active;
654:   PetscFunctionReturn(PETSC_SUCCESS);
655: }

657: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
658: {
659:   PetscFunctionBegin;
660:   if (!tr->coordFE[ct]) {
661:     PetscInt dim, cdim;

663:     dim = DMPolytopeTypeGetDim(ct);
664:     PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
665:     PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
666:     {
667:       PetscDualSpace  dsp;
668:       PetscQuadrature quad;
669:       DM              K;
670:       PetscFEGeom    *cg;
671:       PetscScalar    *Xq;
672:       PetscReal      *xq, *wq;
673:       PetscInt        Nq, q;

675:       PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
676:       PetscCall(PetscMalloc1(Nq * cdim, &xq));
677:       for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
678:       PetscCall(PetscMalloc1(Nq, &wq));
679:       for (q = 0; q < Nq; ++q) wq[q] = 1.0;
680:       PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
681:       PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
682:       PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));

684:       PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
685:       PetscCall(PetscDualSpaceGetDM(dsp, &K));
686:       PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]));
687:       cg = tr->refGeom[ct];
688:       PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
689:       PetscCall(PetscQuadratureDestroy(&quad));
690:     }
691:   }
692:   *fe = tr->coordFE[ct];
693:   PetscFunctionReturn(PETSC_SUCCESS);
694: }

696: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
697: {
698:   PetscInt dim, cdim;

700:   PetscFunctionBegin;
701:   PetscCall(DMGetDimension(dm, &dim));
702:   PetscCall(DMSetDimension(tdm, dim));
703:   PetscCall(DMGetCoordinateDim(dm, &cdim));
704:   PetscCall(DMSetCoordinateDim(tdm, cdim));
705:   PetscFunctionReturn(PETSC_SUCCESS);
706: }

708: /*@
709:   DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`

711:   Input Parameters:
712: + tr - The `DMPlexTransform` object
713: - dm - The original `DM`

715:   Output Parameter:
716: . tdm - The transformed `DM`

718:   Level: advanced

720: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
721: @*/
722: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
723: {
724:   PetscFunctionBegin;
725:   PetscUseTypeMethod(tr, setdimensions, dm, tdm);
726:   PetscFunctionReturn(PETSC_SUCCESS);
727: }

729: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
730: {
731:   PetscFunctionBegin;
732:   if (pStart) *pStart = 0;
733:   if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
734:   PetscFunctionReturn(PETSC_SUCCESS);
735: }

737: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
738: {
739:   PetscInt ctNew;

741:   PetscFunctionBegin;
743:   PetscAssertPointer(celltype, 3);
744:   /* TODO Can do bisection since everything is sorted */
745:   for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
746:     PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];

748:     if (cell >= ctSN && cell < ctEN) break;
749:   }
750:   PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
751:   *celltype = (DMPolytopeType)ctNew;
752:   PetscFunctionReturn(PETSC_SUCCESS);
753: }

755: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
756: {
757:   PetscFunctionBegin;
759:   if (start) *start = tr->ctStartNew[celltype];
760:   if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
761:   PetscFunctionReturn(PETSC_SUCCESS);
762: }

764: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
765: {
766:   PetscFunctionBegin;
768:   *depth = tr->depth;
769:   PetscFunctionReturn(PETSC_SUCCESS);
770: }

772: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
773: {
774:   PetscFunctionBegin;
776:   if (start) *start = tr->depthStart[depth];
777:   if (end) *end = tr->depthEnd[depth];
778:   PetscFunctionReturn(PETSC_SUCCESS);
779: }

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

784:   Not Collective

786:   Input Parameters:
787: + tr    - The `DMPlexTransform`
788: . ct    - The type of the original point which produces the new point
789: . ctNew - The type of the new point
790: . p     - The original point which produces the new point
791: - r     - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`

793:   Output Parameter:
794: . pNew - The new point number

796:   Level: developer

798: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
799: @*/
800: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
801: {
802:   DMPolytopeType *rct;
803:   PetscInt       *rsize, *cone, *ornt;
804:   PetscInt        rt, Nct, n, off, rp;
805:   DMLabel         trType = tr->trType;
806:   PetscInt        ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
807:   PetscInt        ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
808:   PetscInt        newp = ctSN, cind;

810:   PetscFunctionBeginHot;
811:   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);
812:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
813:   if (trType) {
814:     PetscCall(DMLabelGetValueIndex(trType, rt, &cind));
815:     PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
816:     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);
817:   } else {
818:     cind = ct;
819:     rp   = p - ctS;
820:   }
821:   off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
822:   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);
823:   newp += off;
824:   for (n = 0; n < Nct; ++n) {
825:     if (rct[n] == ctNew) {
826:       if (rsize[n] && r >= rsize[n])
827:         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]);
828:       newp += rp * rsize[n] + r;
829:       break;
830:     }
831:   }

833:   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);
834:   *pNew = newp;
835:   PetscFunctionReturn(PETSC_SUCCESS);
836: }

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

841:   Not Collective

843:   Input Parameters:
844: + tr   - The `DMPlexTransform`
845: - pNew - The new point number

847:   Output Parameters:
848: + ct    - The type of the original point which produces the new point
849: . ctNew - The type of the new point
850: . p     - The original point which produces the new point
851: - r     - The replica number of the new point, meaning it is the rth point of type ctNew produced from p

853:   Level: developer

855: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
856: @*/
857: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
858: {
859:   DMLabel         trType = tr->trType;
860:   DMPolytopeType *rct, ctN;
861:   PetscInt       *rsize, *cone, *ornt;
862:   PetscInt        rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
863:   PetscInt        offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;

865:   PetscFunctionBegin;
866:   PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
867:   if (trType) {
868:     DM              dm;
869:     IS              rtIS;
870:     const PetscInt *reftypes;
871:     PetscInt        Nrt, r, rtStart;

873:     PetscCall(DMPlexTransformGetDM(tr, &dm));
874:     PetscCall(DMLabelGetNumValues(trType, &Nrt));
875:     PetscCall(DMLabelGetValueIS(trType, &rtIS));
876:     PetscCall(ISGetIndices(rtIS, &reftypes));
877:     for (r = 0; r < Nrt; ++r) {
878:       const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];

880:       if (tr->ctStartNew[ctN] + off > pNew) continue;
881:       /* Check that any of this refinement type exist */
882:       /* TODO Actually keep track of the number produced here instead */
883:       if (off > offset) {
884:         rt     = reftypes[r];
885:         offset = off;
886:       }
887:     }
888:     PetscCall(ISRestoreIndices(rtIS, &reftypes));
889:     PetscCall(ISDestroy(&rtIS));
890:     PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
891:     /* TODO Map refinement types to cell types */
892:     PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL));
893:     PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
894:     for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
895:       PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];

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

904:       if (tr->ctStartNew[ctN] + off > pNew) continue;
905:       if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
906:       /* TODO Actually keep track of the number produced here instead */
907:       if (off > offset) {
908:         ctO    = ctTmp;
909:         offset = off;
910:       }
911:     }
912:     rt = -1;
913:     PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
914:   }
915:   ctS = tr->ctStart[ctO];
916:   ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
917:   if (trType) {
918:     for (rtS = ctS; rtS < ctE; ++rtS) {
919:       PetscInt val;
920:       PetscCall(DMLabelGetValue(trType, rtS, &val));
921:       if (val == rt) break;
922:     }
923:     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);
924:   } else rtS = ctS;
925:   PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
926:   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);
927:   for (n = 0; n < Nct; ++n) {
928:     if (rct[n] == ctN) {
929:       PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c;

931:       if (trType) {
932:         for (c = ctS; c < ctE; ++c) {
933:           PetscCall(DMLabelGetValue(trType, c, &val));
934:           if (val == rt) {
935:             if (tmp < rsize[n]) break;
936:             tmp -= rsize[n];
937:           }
938:         }
939:         PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew);
940:         rp = c - ctS;
941:         rO = tmp;
942:       } else {
943:         // This assumes that all points of type ctO transform the same way
944:         rp = tmp / rsize[n];
945:         rO = tmp % rsize[n];
946:       }
947:       break;
948:     }
949:   }
950:   PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
951:   pO = rp + ctS;
952:   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);
953:   if (ct) *ct = (DMPolytopeType)ctO;
954:   if (ctNew) *ctNew = (DMPolytopeType)ctN;
955:   if (p) *p = pO;
956:   if (r) *r = rO;
957:   PetscFunctionReturn(PETSC_SUCCESS);
958: }

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

963:   Input Parameters:
964: + tr     - The `DMPlexTransform` object
965: . source - The source cell type
966: - p      - The source point, which can also determine the refine type

968:   Output Parameters:
969: + rt     - The refine type for this point
970: . Nt     - The number of types produced by this point
971: . target - An array of length `Nt` giving the types produced
972: . size   - An array of length `Nt` giving the number of cells of each type produced
973: . cone   - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
974: - ornt   - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point

976:   Level: advanced

978:   Notes:
979:   The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
980:   need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
981:   division (described in these outputs) of the cell in the original mesh. The point identifier is given by
982: .vb
983:    the number of cones to be taken, or 0 for the current cell
984:    the cell cone point number at each level from which it is subdivided
985:    the replica number r of the subdivision.
986: .ve
987:   The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
988: .vb
989:    Nt     = 2
990:    target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
991:    size   = {1, 2}
992:    cone   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,  DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
993:    ornt   = {                         0,                       0,                        0,                          0}
994: .ve

996: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
997: @*/
998: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
999: {
1000:   PetscFunctionBegin;
1001:   PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
1002:   PetscFunctionReturn(PETSC_SUCCESS);
1003: }

1005: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1006: {
1007:   PetscFunctionBegin;
1008:   *rnew = r;
1009:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
1010:   PetscFunctionReturn(PETSC_SUCCESS);
1011: }

1013: /* Returns the same thing */
1014: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1015: {
1016:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
1017:   static PetscInt       vertexS[] = {1};
1018:   static PetscInt       vertexC[] = {0};
1019:   static PetscInt       vertexO[] = {0};
1020:   static DMPolytopeType edgeT[]   = {DM_POLYTOPE_SEGMENT};
1021:   static PetscInt       edgeS[]   = {1};
1022:   static PetscInt       edgeC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1023:   static PetscInt       edgeO[]   = {0, 0};
1024:   static DMPolytopeType tedgeT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR};
1025:   static PetscInt       tedgeS[]  = {1};
1026:   static PetscInt       tedgeC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1027:   static PetscInt       tedgeO[]  = {0, 0};
1028:   static DMPolytopeType triT[]    = {DM_POLYTOPE_TRIANGLE};
1029:   static PetscInt       triS[]    = {1};
1030:   static PetscInt       triC[]    = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1031:   static PetscInt       triO[]    = {0, 0, 0};
1032:   static DMPolytopeType quadT[]   = {DM_POLYTOPE_QUADRILATERAL};
1033:   static PetscInt       quadS[]   = {1};
1034:   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};
1035:   static PetscInt       quadO[]   = {0, 0, 0, 0};
1036:   static DMPolytopeType tquadT[]  = {DM_POLYTOPE_SEG_PRISM_TENSOR};
1037:   static PetscInt       tquadS[]  = {1};
1038:   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};
1039:   static PetscInt       tquadO[]  = {0, 0, 0, 0};
1040:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_TETRAHEDRON};
1041:   static PetscInt       tetS[]    = {1};
1042:   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};
1043:   static PetscInt       tetO[]    = {0, 0, 0, 0};
1044:   static DMPolytopeType hexT[]    = {DM_POLYTOPE_HEXAHEDRON};
1045:   static PetscInt       hexS[]    = {1};
1046:   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};
1047:   static PetscInt       hexO[] = {0, 0, 0, 0, 0, 0};
1048:   static DMPolytopeType tripT[]   = {DM_POLYTOPE_TRI_PRISM};
1049:   static PetscInt       tripS[]   = {1};
1050:   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};
1051:   static PetscInt       tripO[]   = {0, 0, 0, 0, 0};
1052:   static DMPolytopeType ttripT[]  = {DM_POLYTOPE_TRI_PRISM_TENSOR};
1053:   static PetscInt       ttripS[]  = {1};
1054:   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};
1055:   static PetscInt       ttripO[]  = {0, 0, 0, 0, 0};
1056:   static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
1057:   static PetscInt       tquadpS[] = {1};
1058:   static PetscInt       tquadpC[] = {DM_POLYTOPE_QUADRILATERAL,    1, 0, 0, DM_POLYTOPE_QUADRILATERAL,    1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
1059:                                      DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1060:   static PetscInt       tquadpO[] = {0, 0, 0, 0, 0, 0};
1061:   static DMPolytopeType pyrT[]    = {DM_POLYTOPE_PYRAMID};
1062:   static PetscInt       pyrS[]    = {1};
1063:   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};
1064:   static PetscInt       pyrO[]    = {0, 0, 0, 0, 0};

1066:   PetscFunctionBegin;
1067:   if (rt) *rt = 0;
1068:   switch (source) {
1069:   case DM_POLYTOPE_POINT:
1070:     *Nt     = 1;
1071:     *target = vertexT;
1072:     *size   = vertexS;
1073:     *cone   = vertexC;
1074:     *ornt   = vertexO;
1075:     break;
1076:   case DM_POLYTOPE_SEGMENT:
1077:     *Nt     = 1;
1078:     *target = edgeT;
1079:     *size   = edgeS;
1080:     *cone   = edgeC;
1081:     *ornt   = edgeO;
1082:     break;
1083:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1084:     *Nt     = 1;
1085:     *target = tedgeT;
1086:     *size   = tedgeS;
1087:     *cone   = tedgeC;
1088:     *ornt   = tedgeO;
1089:     break;
1090:   case DM_POLYTOPE_TRIANGLE:
1091:     *Nt     = 1;
1092:     *target = triT;
1093:     *size   = triS;
1094:     *cone   = triC;
1095:     *ornt   = triO;
1096:     break;
1097:   case DM_POLYTOPE_QUADRILATERAL:
1098:     *Nt     = 1;
1099:     *target = quadT;
1100:     *size   = quadS;
1101:     *cone   = quadC;
1102:     *ornt   = quadO;
1103:     break;
1104:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1105:     *Nt     = 1;
1106:     *target = tquadT;
1107:     *size   = tquadS;
1108:     *cone   = tquadC;
1109:     *ornt   = tquadO;
1110:     break;
1111:   case DM_POLYTOPE_TETRAHEDRON:
1112:     *Nt     = 1;
1113:     *target = tetT;
1114:     *size   = tetS;
1115:     *cone   = tetC;
1116:     *ornt   = tetO;
1117:     break;
1118:   case DM_POLYTOPE_HEXAHEDRON:
1119:     *Nt     = 1;
1120:     *target = hexT;
1121:     *size   = hexS;
1122:     *cone   = hexC;
1123:     *ornt   = hexO;
1124:     break;
1125:   case DM_POLYTOPE_TRI_PRISM:
1126:     *Nt     = 1;
1127:     *target = tripT;
1128:     *size   = tripS;
1129:     *cone   = tripC;
1130:     *ornt   = tripO;
1131:     break;
1132:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
1133:     *Nt     = 1;
1134:     *target = ttripT;
1135:     *size   = ttripS;
1136:     *cone   = ttripC;
1137:     *ornt   = ttripO;
1138:     break;
1139:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1140:     *Nt     = 1;
1141:     *target = tquadpT;
1142:     *size   = tquadpS;
1143:     *cone   = tquadpC;
1144:     *ornt   = tquadpO;
1145:     break;
1146:   case DM_POLYTOPE_PYRAMID:
1147:     *Nt     = 1;
1148:     *target = pyrT;
1149:     *size   = pyrS;
1150:     *cone   = pyrC;
1151:     *ornt   = pyrO;
1152:     break;
1153:   default:
1154:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1155:   }
1156:   PetscFunctionReturn(PETSC_SUCCESS);
1157: }

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

1162:   Not Collective

1164:   Input Parameters:
1165: + tr  - The `DMPlexTransform`
1166: . sct - The source point cell type, from whom the new cell is being produced
1167: . sp  - The source point
1168: . so  - The orientation of the source point in its enclosing parent
1169: . tct - The target point cell type
1170: . r   - The replica number requested for the produced cell type
1171: - o   - The orientation of the replica

1173:   Output Parameters:
1174: + rnew - The replica number, given the orientation of the parent
1175: - onew - The replica orientation, given the orientation of the parent

1177:   Level: advanced

1179: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1180: @*/
1181: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1182: {
1183:   PetscFunctionBeginHot;
1184:   PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1185:   PetscFunctionReturn(PETSC_SUCCESS);
1186: }

1188: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1189: {
1190:   DM       dm;
1191:   PetscInt pStart, pEnd, p, pNew;

1193:   PetscFunctionBegin;
1194:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1195:   /* Must create the celltype label here so that we do not automatically try to compute the types */
1196:   PetscCall(DMCreateLabel(rdm, "celltype"));
1197:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1198:   for (p = pStart; p < pEnd; ++p) {
1199:     DMPolytopeType  ct;
1200:     DMPolytopeType *rct;
1201:     PetscInt       *rsize, *rcone, *rornt;
1202:     PetscInt        Nct, n, r;

1204:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1205:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1206:     for (n = 0; n < Nct; ++n) {
1207:       for (r = 0; r < rsize[n]; ++r) {
1208:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1209:         PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
1210:         PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
1211:       }
1212:     }
1213:   }
1214:   /* Let the DM know we have set all the cell types */
1215:   {
1216:     DMLabel  ctLabel;
1217:     DM_Plex *plex = (DM_Plex *)rdm->data;

1219:     PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
1220:     PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
1221:   }
1222:   PetscFunctionReturn(PETSC_SUCCESS);
1223: }

1225: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1226: {
1227:   DMPolytopeType ctNew;

1229:   PetscFunctionBegin;
1231:   PetscAssertPointer(coneSize, 3);
1232:   PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
1233:   *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType)ctNew);
1234:   PetscFunctionReturn(PETSC_SUCCESS);
1235: }

1237: /* The orientation o is for the interior of the cell p */
1238: 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[])
1239: {
1240:   DM              dm;
1241:   const PetscInt  csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1242:   const PetscInt *cone;
1243:   PetscInt        c, coff = *coneoff, ooff = *orntoff;

1245:   PetscFunctionBegin;
1246:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1247:   PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
1248:   for (c = 0; c < csizeNew; ++c) {
1249:     PetscInt             ppp   = -1;                            /* Parent Parent point: Parent of point pp */
1250:     PetscInt             pp    = p;                             /* Parent point: Point in the original mesh producing new cone point */
1251:     PetscInt             po    = 0;                             /* Orientation of parent point pp in parent parent point ppp */
1252:     DMPolytopeType       pct   = ct;                            /* Parent type: Cell type for parent of new cone point */
1253:     const PetscInt      *pcone = cone;                          /* Parent cone: Cone of parent point pp */
1254:     PetscInt             pr    = -1;                            /* Replica number of pp that produces new cone point  */
1255:     const DMPolytopeType ft    = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1256:     const PetscInt       fn    = rcone[coff++];                 /* Number of cones of p that need to be taken when producing new cone point */
1257:     PetscInt             fo    = rornt[ooff++];                 /* Orientation of new cone point in pNew */
1258:     PetscInt             lc;

1260:     /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1261:     for (lc = 0; lc < fn; ++lc) {
1262:       const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po);
1263:       const PetscInt  acp  = rcone[coff++];
1264:       const PetscInt  pcp  = parr[acp * 2];
1265:       const PetscInt  pco  = parr[acp * 2 + 1];
1266:       const PetscInt *ppornt;

1268:       ppp = pp;
1269:       pp  = pcone[pcp];
1270:       PetscCall(DMPlexGetCellType(dm, pp, &pct));
1271:       // Restore the parent cone from the last iterate
1272:       if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
1273:       PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
1274:       PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
1275:       po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1276:       PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
1277:     }
1278:     if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
1279:     pr = rcone[coff++];
1280:     /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1281:     PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
1282:     PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
1283:     orntNew[c] = fo;
1284:   }
1285:   PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
1286:   *coneoff = coff;
1287:   *orntoff = ooff;
1288:   PetscFunctionReturn(PETSC_SUCCESS);
1289: }

1291: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1292: {
1293:   DM             dm;
1294:   DMPolytopeType ct;
1295:   PetscInt      *coneNew, *orntNew;
1296:   PetscInt       maxConeSize = 0, pStart, pEnd, p, pNew;

1298:   PetscFunctionBegin;
1299:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1300:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1301:   PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1302:   PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1303:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1304:   for (p = pStart; p < pEnd; ++p) {
1305:     PetscInt        coff, ooff;
1306:     DMPolytopeType *rct;
1307:     PetscInt       *rsize, *rcone, *rornt;
1308:     PetscInt        Nct, n, r;

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

1315:       for (r = 0; r < rsize[n]; ++r) {
1316:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1317:         PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
1318:         PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
1319:         PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
1320:       }
1321:     }
1322:   }
1323:   PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1324:   PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1325:   PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
1326:   PetscCall(DMPlexSymmetrize(rdm));
1327:   PetscCall(DMPlexStratify(rdm));
1328:   PetscFunctionReturn(PETSC_SUCCESS);
1329: }

1331: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1332: {
1333:   DM              dm;
1334:   DMPolytopeType  ct, qct;
1335:   DMPolytopeType *rct;
1336:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
1337:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

1339:   PetscFunctionBegin;
1341:   PetscAssertPointer(cone, 4);
1342:   PetscAssertPointer(ornt, 5);
1343:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1344:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1345:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1346:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1347:   PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1348:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1349:   for (n = 0; n < Nct; ++n) {
1350:     const DMPolytopeType ctNew    = rct[n];
1351:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1352:     PetscInt             Nr       = rsize[n], fn, c;

1354:     if (ctNew == qct) Nr = r;
1355:     for (nr = 0; nr < Nr; ++nr) {
1356:       for (c = 0; c < csizeNew; ++c) {
1357:         ++coff;             /* Cell type of new cone point */
1358:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1359:         coff += fn;
1360:         ++coff; /* Replica number of new cone point */
1361:         ++ooff; /* Orientation of new cone point */
1362:       }
1363:     }
1364:     if (ctNew == qct) break;
1365:   }
1366:   PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1367:   *cone = qcone;
1368:   *ornt = qornt;
1369:   PetscFunctionReturn(PETSC_SUCCESS);
1370: }

1372: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1373: {
1374:   DM              dm;
1375:   DMPolytopeType  ct, qct;
1376:   DMPolytopeType *rct;
1377:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
1378:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

1380:   PetscFunctionBegin;
1382:   if (cone) PetscAssertPointer(cone, 3);
1383:   if (ornt) PetscAssertPointer(ornt, 4);
1384:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1385:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1386:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1387:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1388:   PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1389:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1390:   for (n = 0; n < Nct; ++n) {
1391:     const DMPolytopeType ctNew    = rct[n];
1392:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1393:     PetscInt             Nr       = rsize[n], fn, c;

1395:     if (ctNew == qct) Nr = r;
1396:     for (nr = 0; nr < Nr; ++nr) {
1397:       for (c = 0; c < csizeNew; ++c) {
1398:         ++coff;             /* Cell type of new cone point */
1399:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1400:         coff += fn;
1401:         ++coff; /* Replica number of new cone point */
1402:         ++ooff; /* Orientation of new cone point */
1403:       }
1404:     }
1405:     if (ctNew == qct) break;
1406:   }
1407:   PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1408:   if (cone) *cone = qcone;
1409:   else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1410:   if (ornt) *ornt = qornt;
1411:   else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1412:   PetscFunctionReturn(PETSC_SUCCESS);
1413: }

1415: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1416: {
1417:   DM dm;

1419:   PetscFunctionBegin;
1421:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1422:   if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
1423:   if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
1424:   PetscFunctionReturn(PETSC_SUCCESS);
1425: }

1427: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1428: {
1429:   PetscInt ict;

1431:   PetscFunctionBegin;
1432:   PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
1433:   for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1434:     const DMPolytopeType ct = (DMPolytopeType)ict;
1435:     DMPlexTransform      reftr;
1436:     DM                   refdm, trdm;
1437:     Vec                  coordinates;
1438:     const PetscScalar   *coords;
1439:     DMPolytopeType      *rct;
1440:     PetscInt            *rsize, *rcone, *rornt;
1441:     PetscInt             Nct, n, r, pNew;
1442:     PetscInt             trdim, vStart, vEnd, Nc;
1443:     const PetscInt       debug = 0;
1444:     const char          *typeName;

1446:     /* Since points are 0-dimensional, coordinates make no sense */
1447:     if (DMPolytopeTypeGetDim(ct) <= 0) continue;
1448:     PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
1449:     PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
1450:     PetscCall(DMPlexTransformSetDM(reftr, refdm));
1451:     PetscCall(DMPlexTransformGetType(tr, &typeName));
1452:     PetscCall(DMPlexTransformSetType(reftr, typeName));
1453:     PetscCall(DMPlexTransformSetUp(reftr));
1454:     PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));

1456:     PetscCall(DMGetDimension(trdm, &trdim));
1457:     PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
1458:     tr->trNv[ct] = vEnd - vStart;
1459:     PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
1460:     PetscCall(VecGetLocalSize(coordinates, &Nc));
1461:     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);
1462:     PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
1463:     PetscCall(VecGetArrayRead(coordinates, &coords));
1464:     PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
1465:     PetscCall(VecRestoreArrayRead(coordinates, &coords));

1467:     PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
1468:     PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1469:     for (n = 0; n < Nct; ++n) {
1470:       /* Since points are 0-dimensional, coordinates make no sense */
1471:       if (rct[n] == DM_POLYTOPE_POINT) continue;
1472:       PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
1473:       for (r = 0; r < rsize[n]; ++r) {
1474:         PetscInt *closure = NULL;
1475:         PetscInt  clSize, cl, Nv = 0;

1477:         PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
1478:         PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
1479:         PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1480:         for (cl = 0; cl < clSize * 2; cl += 2) {
1481:           const PetscInt sv = closure[cl];

1483:           if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1484:         }
1485:         PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1486:         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]);
1487:       }
1488:     }
1489:     if (debug) {
1490:       DMPolytopeType *rct;
1491:       PetscInt       *rsize, *rcone, *rornt;
1492:       PetscInt        v, dE = trdim, d, off = 0;

1494:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
1495:       for (v = 0; v < tr->trNv[ct]; ++v) {
1496:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  "));
1497:         for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
1498:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1499:       }

1501:       PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1502:       for (n = 0; n < Nct; ++n) {
1503:         if (rct[n] == DM_POLYTOPE_POINT) continue;
1504:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
1505:         for (r = 0; r < rsize[n]; ++r) {
1506:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  "));
1507:           for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
1508:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1509:         }
1510:       }
1511:     }
1512:     PetscCall(DMDestroy(&refdm));
1513:     PetscCall(DMDestroy(&trdm));
1514:     PetscCall(DMPlexTransformDestroy(&reftr));
1515:   }
1516:   PetscFunctionReturn(PETSC_SUCCESS);
1517: }

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

1522:   Input Parameters:
1523: + tr - The `DMPlexTransform` object
1524: - ct - The cell type

1526:   Output Parameters:
1527: + Nv      - The number of transformed vertices in the closure of the reference cell of given type
1528: - trVerts - The coordinates of these vertices in the reference cell

1530:   Level: developer

1532: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
1533: @*/
1534: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1535: {
1536:   PetscFunctionBegin;
1537:   if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1538:   if (Nv) *Nv = tr->trNv[ct];
1539:   if (trVerts) *trVerts = tr->trVerts[ct];
1540:   PetscFunctionReturn(PETSC_SUCCESS);
1541: }

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

1546:   Input Parameters:
1547: + tr  - The `DMPlexTransform` object
1548: . ct  - The cell type
1549: . rct - The subcell type
1550: - r   - The subcell index

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

1555:   Level: developer

1557: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
1558: @*/
1559: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1560: {
1561:   PetscFunctionBegin;
1562:   if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1563:   PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
1564:   if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1565:   PetscFunctionReturn(PETSC_SUCCESS);
1566: }

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

1573:   PetscFunctionBeginHot;
1574:   PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
1575:   for (d = 0; d < dE; ++d) out[d] = 0.0;
1576:   for (v = 0; v < Nv; ++v)
1577:     for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1578:   for (d = 0; d < dE; ++d) out[d] /= Nv;
1579:   PetscFunctionReturn(PETSC_SUCCESS);
1580: }

1582: /*@
1583:   DMPlexTransformMapCoordinates - Calculate new coordinates for produced points

1585:   Not collective

1587:   Input Parameters:
1588: + tr  - The `DMPlexTransform`
1589: . pct - The cell type of the parent, from whom the new cell is being produced
1590: . ct  - The type being produced
1591: . p   - The original point
1592: . r   - The replica number requested for the produced cell type
1593: . Nv  - Number of vertices in the closure of the parent cell
1594: . dE  - Spatial dimension
1595: - in  - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell

1597:   Output Parameter:
1598: . out - The coordinates of the new vertices

1600:   Level: intermediate

1602: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
1603: @*/
1604: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1605: {
1606:   PetscFunctionBeginHot;
1607:   if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1608:   PetscFunctionReturn(PETSC_SUCCESS);
1609: }

1611: /*
1612:   DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label

1614:   Not Collective

1616:   Input Parameters:
1617: + tr    - The `DMPlexTransform`
1618: . label - The label in the transformed mesh
1619: . pp    - The parent point in the original mesh
1620: . pct   - The cell type of the parent point
1621: . p     - The point in the transformed mesh
1622: . ct    - The cell type of the point
1623: . r     - The replica number of the point
1624: - val   - The label value of the parent point

1626:   Level: developer

1628: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
1629: */
1630: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
1631: {
1632:   PetscFunctionBeginHot;
1633:   if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
1634:   PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
1635:   PetscFunctionReturn(PETSC_SUCCESS);
1636: }

1638: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1639: {
1640:   DM              dm;
1641:   IS              valueIS;
1642:   const PetscInt *values;
1643:   PetscInt        defVal, Nv, val;

1645:   PetscFunctionBegin;
1646:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1647:   PetscCall(DMLabelGetDefaultValue(label, &defVal));
1648:   PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
1649:   PetscCall(DMLabelGetValueIS(label, &valueIS));
1650:   PetscCall(ISGetLocalSize(valueIS, &Nv));
1651:   PetscCall(ISGetIndices(valueIS, &values));
1652:   for (val = 0; val < Nv; ++val) {
1653:     IS              pointIS;
1654:     const PetscInt *points;
1655:     PetscInt        numPoints, p;

1657:     /* Ensure refined label is created with same number of strata as
1658:      * original (even if no entries here). */
1659:     PetscCall(DMLabelAddStratum(labelNew, values[val]));
1660:     PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
1661:     PetscCall(ISGetLocalSize(pointIS, &numPoints));
1662:     PetscCall(ISGetIndices(pointIS, &points));
1663:     for (p = 0; p < numPoints; ++p) {
1664:       const PetscInt  point = points[p];
1665:       DMPolytopeType  ct;
1666:       DMPolytopeType *rct;
1667:       PetscInt       *rsize, *rcone, *rornt;
1668:       PetscInt        Nct, n, r, pNew = 0;

1670:       PetscCall(DMPlexGetCellType(dm, point, &ct));
1671:       PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1672:       for (n = 0; n < Nct; ++n) {
1673:         for (r = 0; r < rsize[n]; ++r) {
1674:           PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
1675:           PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
1676:         }
1677:       }
1678:     }
1679:     PetscCall(ISRestoreIndices(pointIS, &points));
1680:     PetscCall(ISDestroy(&pointIS));
1681:   }
1682:   PetscCall(ISRestoreIndices(valueIS, &values));
1683:   PetscCall(ISDestroy(&valueIS));
1684:   PetscFunctionReturn(PETSC_SUCCESS);
1685: }

1687: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1688: {
1689:   DM       dm;
1690:   PetscInt numLabels, l;

1692:   PetscFunctionBegin;
1693:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1694:   PetscCall(DMGetNumLabels(dm, &numLabels));
1695:   for (l = 0; l < numLabels; ++l) {
1696:     DMLabel     label, labelNew;
1697:     const char *lname;
1698:     PetscBool   isDepth, isCellType;

1700:     PetscCall(DMGetLabelName(dm, l, &lname));
1701:     PetscCall(PetscStrcmp(lname, "depth", &isDepth));
1702:     if (isDepth) continue;
1703:     PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
1704:     if (isCellType) continue;
1705:     PetscCall(DMCreateLabel(rdm, lname));
1706:     PetscCall(DMGetLabel(dm, lname, &label));
1707:     PetscCall(DMGetLabel(rdm, lname, &labelNew));
1708:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
1709:   }
1710:   PetscFunctionReturn(PETSC_SUCCESS);
1711: }

1713: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1714: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1715: {
1716:   DM       dm;
1717:   PetscInt Nf, f, Nds, s;

1719:   PetscFunctionBegin;
1720:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1721:   PetscCall(DMGetNumFields(dm, &Nf));
1722:   for (f = 0; f < Nf; ++f) {
1723:     DMLabel     label, labelNew;
1724:     PetscObject obj;
1725:     const char *lname;

1727:     PetscCall(DMGetField(rdm, f, &label, &obj));
1728:     if (!label) continue;
1729:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1730:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1731:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
1732:     PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
1733:     PetscCall(DMLabelDestroy(&labelNew));
1734:   }
1735:   PetscCall(DMGetNumDS(dm, &Nds));
1736:   for (s = 0; s < Nds; ++s) {
1737:     DMLabel     label, labelNew;
1738:     const char *lname;

1740:     PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
1741:     if (!label) continue;
1742:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1743:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1744:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
1745:     PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
1746:     PetscCall(DMLabelDestroy(&labelNew));
1747:   }
1748:   PetscFunctionReturn(PETSC_SUCCESS);
1749: }

1751: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1752: {
1753:   DM                 dm;
1754:   PetscSF            sf, sfNew;
1755:   PetscInt           numRoots, numLeaves, numLeavesNew = 0, l, m;
1756:   const PetscInt    *localPoints;
1757:   const PetscSFNode *remotePoints;
1758:   PetscInt          *localPointsNew;
1759:   PetscSFNode       *remotePointsNew;
1760:   PetscInt           pStartNew, pEndNew, pNew;
1761:   /* Brute force algorithm */
1762:   PetscSF         rsf;
1763:   PetscSection    s;
1764:   const PetscInt *rootdegree;
1765:   PetscInt       *rootPointsNew, *remoteOffsets;
1766:   PetscInt        numPointsNew, pStart, pEnd, p;

1768:   PetscFunctionBegin;
1769:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1770:   PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
1771:   PetscCall(DMGetPointSF(dm, &sf));
1772:   PetscCall(DMGetPointSF(rdm, &sfNew));
1773:   /* Calculate size of new SF */
1774:   PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
1775:   if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
1776:   for (l = 0; l < numLeaves; ++l) {
1777:     const PetscInt  p = localPoints[l];
1778:     DMPolytopeType  ct;
1779:     DMPolytopeType *rct;
1780:     PetscInt       *rsize, *rcone, *rornt;
1781:     PetscInt        Nct, n;

1783:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1784:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1785:     for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
1786:   }
1787:   /* Send new root point numbers
1788:        It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1789:   */
1790:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1791:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
1792:   PetscCall(PetscSectionSetChart(s, pStart, pEnd));
1793:   for (p = pStart; p < pEnd; ++p) {
1794:     DMPolytopeType  ct;
1795:     DMPolytopeType *rct;
1796:     PetscInt       *rsize, *rcone, *rornt;
1797:     PetscInt        Nct, n;

1799:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1800:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1801:     for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
1802:   }
1803:   PetscCall(PetscSectionSetUp(s));
1804:   PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
1805:   PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
1806:   PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
1807:   PetscCall(PetscFree(remoteOffsets));
1808:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1809:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1810:   PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
1811:   for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1812:   for (p = pStart; p < pEnd; ++p) {
1813:     DMPolytopeType  ct;
1814:     DMPolytopeType *rct;
1815:     PetscInt       *rsize, *rcone, *rornt;
1816:     PetscInt        Nct, n, r, off;

1818:     if (!rootdegree[p - pStart]) continue;
1819:     PetscCall(PetscSectionGetOffset(s, p, &off));
1820:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1821:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1822:     for (n = 0, m = 0; n < Nct; ++n) {
1823:       for (r = 0; r < rsize[n]; ++r, ++m) {
1824:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1825:         rootPointsNew[off + m] = pNew;
1826:       }
1827:     }
1828:   }
1829:   PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1830:   PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1831:   PetscCall(PetscSFDestroy(&rsf));
1832:   PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
1833:   PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
1834:   for (l = 0, m = 0; l < numLeaves; ++l) {
1835:     const PetscInt  p = localPoints[l];
1836:     DMPolytopeType  ct;
1837:     DMPolytopeType *rct;
1838:     PetscInt       *rsize, *rcone, *rornt;
1839:     PetscInt        Nct, n, r, q, off;

1841:     PetscCall(PetscSectionGetOffset(s, p, &off));
1842:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1843:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1844:     for (n = 0, q = 0; n < Nct; ++n) {
1845:       for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1846:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1847:         localPointsNew[m]        = pNew;
1848:         remotePointsNew[m].index = rootPointsNew[off + q];
1849:         remotePointsNew[m].rank  = remotePoints[l].rank;
1850:       }
1851:     }
1852:   }
1853:   PetscCall(PetscSectionDestroy(&s));
1854:   PetscCall(PetscFree(rootPointsNew));
1855:   /* SF needs sorted leaves to correctly calculate Gather */
1856:   {
1857:     PetscSFNode *rp, *rtmp;
1858:     PetscInt    *lp, *idx, *ltmp, i;

1860:     PetscCall(PetscMalloc1(numLeavesNew, &idx));
1861:     PetscCall(PetscMalloc1(numLeavesNew, &lp));
1862:     PetscCall(PetscMalloc1(numLeavesNew, &rp));
1863:     for (i = 0; i < numLeavesNew; ++i) {
1864:       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);
1865:       idx[i] = i;
1866:     }
1867:     PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
1868:     for (i = 0; i < numLeavesNew; ++i) {
1869:       lp[i] = localPointsNew[idx[i]];
1870:       rp[i] = remotePointsNew[idx[i]];
1871:     }
1872:     ltmp            = localPointsNew;
1873:     localPointsNew  = lp;
1874:     rtmp            = remotePointsNew;
1875:     remotePointsNew = rp;
1876:     PetscCall(PetscFree(idx));
1877:     PetscCall(PetscFree(ltmp));
1878:     PetscCall(PetscFree(rtmp));
1879:   }
1880:   PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
1881:   PetscFunctionReturn(PETSC_SUCCESS);
1882: }

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

1887:   Not Collective

1889:   Input Parameters:
1890: + tr  - The `DMPlexTransform`
1891: . ct  - The type of the parent cell
1892: . rct - The type of the produced cell
1893: . r   - The index of the produced cell
1894: - x   - The localized coordinates for the parent cell

1896:   Output Parameter:
1897: . xr  - The localized coordinates for the produced cell

1899:   Level: developer

1901: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
1902: */
1903: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1904: {
1905:   PetscFE  fe = NULL;
1906:   PetscInt cdim, v, *subcellV;

1908:   PetscFunctionBegin;
1909:   PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
1910:   PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
1911:   PetscCall(PetscFEGetNumComponents(fe, &cdim));
1912:   for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
1913:   PetscFunctionReturn(PETSC_SUCCESS);
1914: }

1916: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1917: {
1918:   DM                 dm, cdm, cdmCell, cdmNew, cdmCellNew;
1919:   PetscSection       coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
1920:   Vec                coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
1921:   const PetscScalar *coords;
1922:   PetscScalar       *coordsNew;
1923:   const PetscReal   *maxCell, *Lstart, *L;
1924:   PetscBool          localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1925:   PetscInt           dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;

1927:   PetscFunctionBegin;
1928:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1929:   PetscCall(DMGetCoordinateDM(dm, &cdm));
1930:   PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
1931:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
1932:   PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
1933:   if (localized) {
1934:     /* Localize coordinates of new vertices */
1935:     localizeVertices = PETSC_TRUE;
1936:     /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
1937:     if (!maxCell) localizeCells = PETSC_TRUE;
1938:   }
1939:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
1940:   PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
1941:   if (maxCell) {
1942:     PetscReal maxCellNew[3];

1944:     for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
1945:     PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L));
1946:   }
1947:   PetscCall(DMGetCoordinateDim(rdm, &dE));
1948:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
1949:   PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
1950:   PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
1951:   PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
1952:   PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
1953:   /* Localization should be inherited */
1954:   /*   Stefano calculates parent cells for each new cell for localization */
1955:   /*   Localized cells need coordinates of closure */
1956:   for (v = vStartNew; v < vEndNew; ++v) {
1957:     PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
1958:     PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
1959:   }
1960:   PetscCall(PetscSectionSetUp(coordSectionNew));
1961:   PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));

1963:   if (localizeCells) {
1964:     PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
1965:     PetscCall(DMClone(cdmNew, &cdmCellNew));
1966:     PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
1967:     PetscCall(DMDestroy(&cdmCellNew));

1969:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
1970:     PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
1971:     PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
1972:     PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
1973:     PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));

1975:     PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
1976:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1977:     for (c = cStart; c < cEnd; ++c) {
1978:       PetscInt dof;

1980:       PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
1981:       if (dof) {
1982:         DMPolytopeType  ct;
1983:         DMPolytopeType *rct;
1984:         PetscInt       *rsize, *rcone, *rornt;
1985:         PetscInt        dim, cNew, Nct, n, r;

1987:         PetscCall(DMPlexGetCellType(dm, c, &ct));
1988:         dim = DMPolytopeTypeGetDim(ct);
1989:         PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1990:         /* This allows for different cell types */
1991:         for (n = 0; n < Nct; ++n) {
1992:           if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
1993:           for (r = 0; r < rsize[n]; ++r) {
1994:             PetscInt *closure = NULL;
1995:             PetscInt  clSize, cl, Nv = 0;

1997:             PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
1998:             PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
1999:             for (cl = 0; cl < clSize * 2; cl += 2) {
2000:               if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
2001:             }
2002:             PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2003:             PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
2004:             PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
2005:           }
2006:         }
2007:       }
2008:     }
2009:     PetscCall(PetscSectionSetUp(coordSectionCellNew));
2010:     PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
2011:   }
2012:   PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
2013:   {
2014:     VecType     vtype;
2015:     PetscInt    coordSizeNew, bs;
2016:     const char *name;

2018:     PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
2019:     PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
2020:     PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
2021:     PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
2022:     PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
2023:     PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
2024:     PetscCall(VecGetBlockSize(coordsLocal, &bs));
2025:     PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
2026:     PetscCall(VecGetType(coordsLocal, &vtype));
2027:     PetscCall(VecSetType(coordsLocalNew, vtype));
2028:   }
2029:   PetscCall(VecGetArrayRead(coordsLocal, &coords));
2030:   PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
2031:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2032:   /* First set coordinates for vertices */
2033:   for (p = pStart; p < pEnd; ++p) {
2034:     DMPolytopeType  ct;
2035:     DMPolytopeType *rct;
2036:     PetscInt       *rsize, *rcone, *rornt;
2037:     PetscInt        Nct, n, r;
2038:     PetscBool       hasVertex = PETSC_FALSE;

2040:     PetscCall(DMPlexGetCellType(dm, p, &ct));
2041:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2042:     for (n = 0; n < Nct; ++n) {
2043:       if (rct[n] == DM_POLYTOPE_POINT) {
2044:         hasVertex = PETSC_TRUE;
2045:         break;
2046:       }
2047:     }
2048:     if (hasVertex) {
2049:       const PetscScalar *icoords = NULL;
2050:       const PetscScalar *array   = NULL;
2051:       PetscScalar       *pcoords = NULL;
2052:       PetscBool          isDG;
2053:       PetscInt           Nc, Nv, v, d;

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

2057:       icoords = pcoords;
2058:       Nv      = Nc / dEo;
2059:       if (ct != DM_POLYTOPE_POINT) {
2060:         if (localizeVertices && maxCell) {
2061:           PetscScalar anchor[3];

2063:           for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
2064:           for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
2065:         }
2066:       }
2067:       for (n = 0; n < Nct; ++n) {
2068:         if (rct[n] != DM_POLYTOPE_POINT) continue;
2069:         for (r = 0; r < rsize[n]; ++r) {
2070:           PetscScalar vcoords[3];
2071:           PetscInt    vNew, off;

2073:           PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
2074:           PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
2075:           PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
2076:           PetscCall(DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
2077:         }
2078:       }
2079:       PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2080:     }
2081:   }
2082:   PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
2083:   PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
2084:   PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
2085:   PetscCall(VecDestroy(&coordsLocalNew));
2086:   PetscCall(PetscSectionDestroy(&coordSectionNew));
2087:   /* Then set coordinates for cells by localizing */
2088:   if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
2089:   else {
2090:     VecType     vtype;
2091:     PetscInt    coordSizeNew, bs;
2092:     const char *name;

2094:     PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
2095:     PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
2096:     PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
2097:     PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
2098:     PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
2099:     PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
2100:     PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
2101:     PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
2102:     PetscCall(VecGetType(coordsLocalCell, &vtype));
2103:     PetscCall(VecSetType(coordsLocalCellNew, vtype));
2104:     PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
2105:     PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));

2107:     for (p = pStart; p < pEnd; ++p) {
2108:       DMPolytopeType  ct;
2109:       DMPolytopeType *rct;
2110:       PetscInt       *rsize, *rcone, *rornt;
2111:       PetscInt        dof = 0, Nct, n, r;

2113:       PetscCall(DMPlexGetCellType(dm, p, &ct));
2114:       PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2115:       if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
2116:       if (dof) {
2117:         const PetscScalar *pcoords;

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

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

2127:             /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2128:                DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2129:                cell to the ones it produces. */
2130:             PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2131:             PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
2132:             PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
2133:           }
2134:         }
2135:       }
2136:     }
2137:     PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
2138:     PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
2139:     PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
2140:     PetscCall(VecDestroy(&coordsLocalCellNew));
2141:     PetscCall(PetscSectionDestroy(&coordSectionCellNew));
2142:   }
2143:   PetscFunctionReturn(PETSC_SUCCESS);
2144: }

2146: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
2147: {
2148:   DM                     rdm;
2149:   DMPlexInterpolatedFlag interp;
2150:   PetscInt               pStart, pEnd;

2152:   PetscFunctionBegin;
2155:   PetscAssertPointer(tdm, 3);
2156:   PetscCall(PetscLogEventBegin(DMPLEX_Transform, tr, dm, 0, 0));
2157:   PetscCall(DMPlexTransformSetDM(tr, dm));

2159:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
2160:   PetscCall(DMSetType(rdm, DMPLEX));
2161:   PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
2162:   /* Calculate number of new points of each depth */
2163:   PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
2164:   PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
2165:   /* Step 1: Set chart */
2166:   PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
2167:   PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
2168:   /* Step 2: Set cone/support sizes (automatically stratifies) */
2169:   PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
2170:   /* Step 3: Setup refined DM */
2171:   PetscCall(DMSetUp(rdm));
2172:   /* Step 4: Set cones and supports (automatically symmetrizes) */
2173:   PetscCall(DMPlexTransformSetCones(tr, rdm));
2174:   /* Step 5: Create pointSF */
2175:   PetscCall(DMPlexTransformCreateSF(tr, rdm));
2176:   /* Step 6: Create labels */
2177:   PetscCall(DMPlexTransformCreateLabels(tr, rdm));
2178:   /* Step 7: Set coordinates */
2179:   PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
2180:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm));
2181:   // If the original DM was configured from options, the transformed DM should be as well
2182:   rdm->setfromoptionscalled = dm->setfromoptionscalled;
2183:   PetscCall(PetscLogEventEnd(DMPLEX_Transform, tr, dm, 0, 0));
2184:   *tdm = rdm;
2185:   PetscFunctionReturn(PETSC_SUCCESS);
2186: }

2188: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2189: {
2190:   DMPlexTransform tr;
2191:   DM              cdm, rcdm;
2192:   const char     *prefix;

2194:   PetscFunctionBegin;
2195:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
2196:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2197:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
2198:   PetscCall(DMPlexTransformSetDM(tr, dm));
2199:   PetscCall(DMPlexTransformSetFromOptions(tr));
2200:   PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
2201:   PetscCall(DMPlexTransformSetUp(tr));
2202:   PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
2203:   PetscCall(DMPlexTransformApply(tr, dm, rdm));
2204:   PetscCall(DMCopyDisc(dm, *rdm));
2205:   PetscCall(DMGetCoordinateDM(dm, &cdm));
2206:   PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
2207:   PetscCall(DMCopyDisc(cdm, rcdm));
2208:   PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
2209:   PetscCall(DMCopyDisc(dm, *rdm));
2210:   PetscCall(DMPlexTransformDestroy(&tr));
2211:   ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2212:   PetscFunctionReturn(PETSC_SUCCESS);
2213: }