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 || c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) 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 || c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue;
 33:       ctO[off++] = c;
 34:     }
 35:   }
 36:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 37:     if (DMPolytopeTypeGetDim((DMPolytopeType)c) >= 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) 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);
 94: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform);

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

 99:   Not Collective

101:   Level: advanced

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

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

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

126:   Not collective

128:   Level: developer

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

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

143:   Collective

145:   Input Parameter:
146: . comm - The communicator for the transform object

148:   Output Parameter:
149: . tr - The transform object

151:   Level: beginner

153: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
154: @*/
155: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
156: {
157:   DMPlexTransform t;

159:   PetscFunctionBegin;
160:   PetscAssertPointer(tr, 2);
161:   *tr = NULL;
162:   PetscCall(DMInitializePackage());

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

171: /*@
172:   DMPlexTransformSetType - Sets the particular implementation for a transform.

174:   Collective

176:   Input Parameters:
177: + tr     - The transform
178: - method - The name of the transform type

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

183:   Level: intermediate

185: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
186: @*/
187: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
188: {
189:   PetscErrorCode (*r)(DMPlexTransform);
190:   PetscBool match;

192:   PetscFunctionBegin;
194:   PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match));
195:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

197:   PetscCall(DMPlexTransformRegisterAll());
198:   PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r));
199:   PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);

201:   PetscTryTypeMethod(tr, destroy);
202:   PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops)));
203:   PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method));
204:   PetscCall((*r)(tr));
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

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

211:   Not Collective

213:   Input Parameter:
214: . tr - The `DMPlexTransform`

216:   Output Parameter:
217: . type - The `DMPlexTransformType` name

219:   Level: intermediate

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

233: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
234: {
235:   PetscViewerFormat format;

237:   PetscFunctionBegin;
238:   PetscCall(PetscViewerGetFormat(v, &format));
239:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
240:     const PetscInt *trTypes = NULL;
241:     IS              trIS;
242:     PetscInt        cols = 8;
243:     PetscInt        Nrt  = 8, f, g;

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

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

279: /*@C
280:   DMPlexTransformView - Views a `DMPlexTransform`

282:   Collective

284:   Input Parameters:
285: + tr - the `DMPlexTransform` object to view
286: - v  - the viewer

288:   Level: beginner

290: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
291: @*/
292: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
293: {
294:   PetscBool isascii;

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

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

312:   Collective

314:   Input Parameter:
315: . tr - the `DMPlexTransform` object to set options for

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

322:   Level: intermediate

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

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

347: /*@
348:   DMPlexTransformDestroy - Destroys a `DMPlexTransform`

350:   Collective

352:   Input Parameter:
353: . tr - the transform object to destroy

355:   Level: beginner

357: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
358: @*/
359: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
360: {
361:   PetscInt c;

363:   PetscFunctionBegin;
364:   if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
366:   if (--((PetscObject)*tr)->refct > 0) {
367:     *tr = NULL;
368:     PetscFunctionReturn(PETSC_SUCCESS);
369:   }

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

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

408: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
409: {
410:   DMLabel  trType = tr->trType;
411:   PetscInt c, cN, *off;

413:   PetscFunctionBegin;
414:   if (trType) {
415:     DM              dm;
416:     IS              rtIS;
417:     const PetscInt *reftypes;
418:     PetscInt        Nrt, r;

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

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

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

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

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

495:         if (DMPolytopeTypeGetDim(ct) < 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE || DMPolytopeTypeGetDim(ctNew) < 0 || ctNew == DM_POLYTOPE_UNKNOWN_CELL || ctNew == DM_POLYTOPE_UNKNOWN_FACE) {
496:           off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
497:           continue;
498:         }
499:         off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
500:         for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
501:           const DMPolytopeType ict  = (DMPolytopeType)ctOrderOld[i];
502:           const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];

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

521: /*@
522:   DMPlexTransformSetUp - Create the tables that drive the transform

524:   Input Parameter:
525: . tr - The `DMPlexTransform` object

527:   Level: intermediate

529: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
530: @*/
531: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
532: {
533:   DM             dm;
534:   DMPolytopeType ctCell;
535:   PetscInt       pStart, pEnd, p, c, celldim = 0;

537:   PetscFunctionBegin;
539:   if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
540:   PetscTryTypeMethod(tr, setup);
541:   PetscCall(DMPlexTransformGetDM(tr, &dm));
542:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
543:   if (pEnd > pStart) {
544:     PetscCall(DMPlexGetCellType(dm, 0, &ctCell));
545:   } else {
546:     PetscInt dim;

548:     PetscCall(DMGetDimension(dm, &dim));
549:     switch (dim) {
550:     case 0:
551:       ctCell = DM_POLYTOPE_POINT;
552:       break;
553:     case 1:
554:       ctCell = DM_POLYTOPE_SEGMENT;
555:       break;
556:     case 2:
557:       ctCell = DM_POLYTOPE_TRIANGLE;
558:       break;
559:     case 3:
560:       ctCell = DM_POLYTOPE_TETRAHEDRON;
561:       break;
562:     default:
563:       ctCell = DM_POLYTOPE_UNKNOWN;
564:     }
565:   }
566:   PetscCall(DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
567:   for (p = pStart; p < pEnd; ++p) {
568:     DMPolytopeType  ct;
569:     DMPolytopeType *rct;
570:     PetscInt       *rsize, *cone, *ornt;
571:     PetscInt        Nct, n;

573:     PetscCall(DMPlexGetCellType(dm, p, &ct));
574:     PetscCheck(ct != DM_POLYTOPE_UNKNOWN && ct != DM_POLYTOPE_UNKNOWN_CELL && ct != DM_POLYTOPE_UNKNOWN_FACE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p);
575:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
576:     for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
577:   }
578:   PetscCall(DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
579:   /* Construct sizes and offsets for each cell type */
580:   if (!tr->ctStart) {
581:     PetscInt *ctS, *ctSN, *ctC, *ctCN;

583:     PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
584:     PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
585:     for (p = pStart; p < pEnd; ++p) {
586:       DMPolytopeType  ct;
587:       DMPolytopeType *rct;
588:       PetscInt       *rsize, *cone, *ornt;
589:       PetscInt        Nct, n;

591:       PetscCall(DMPlexGetCellType(dm, p, &ct));
592:       PetscCheck(ct != DM_POLYTOPE_UNKNOWN && ct != DM_POLYTOPE_UNKNOWN_CELL && ct != DM_POLYTOPE_UNKNOWN_FACE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p);
593:       ++ctC[ct];
594:       PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
595:       for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
596:     }
597:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
598:       const PetscInt cto  = tr->ctOrderOld[c];
599:       const PetscInt cton = tr->ctOrderOld[c + 1];
600:       const PetscInt ctn  = tr->ctOrderNew[c];
601:       const PetscInt ctnn = tr->ctOrderNew[c + 1];

603:       ctS[cton]  = ctS[cto] + ctC[cto];
604:       ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
605:     }
606:     PetscCall(PetscFree2(ctC, ctCN));
607:     tr->ctStart    = ctS;
608:     tr->ctStartNew = ctSN;
609:   }
610:   PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
611:   // Compute depth information
612:   tr->depth = -1;
613:   for (c = 0; c < DM_NUM_POLYTOPES; ++c)
614:     if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
615:   PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
616:   for (PetscInt d = 0; d <= tr->depth; ++d) {
617:     tr->depthStart[d] = PETSC_MAX_INT;
618:     tr->depthEnd[d]   = -1;
619:   }
620:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
621:     const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);

623:     if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
624:     tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
625:     tr->depthEnd[dep]   = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
626:   }
627:   tr->setupcalled = PETSC_TRUE;
628:   PetscFunctionReturn(PETSC_SUCCESS);
629: }

631: /*@
632:   DMPlexTransformGetDM - Get the base `DM` for the transform

634:   Input Parameter:
635: . tr - The `DMPlexTransform` object

637:   Output Parameter:
638: . dm - The original `DM` which will be transformed

640:   Level: intermediate

642: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
643: @*/
644: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
645: {
646:   PetscFunctionBegin;
648:   PetscAssertPointer(dm, 2);
649:   *dm = tr->dm;
650:   PetscFunctionReturn(PETSC_SUCCESS);
651: }

653: /*@
654:   DMPlexTransformSetDM - Set the base `DM` for the transform

656:   Input Parameters:
657: + tr - The `DMPlexTransform` object
658: - dm - The original `DM` which will be transformed

660:   Level: intermediate

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

665: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
666: @*/
667: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
668: {
669:   PetscFunctionBegin;
672:   PetscCall(PetscObjectReference((PetscObject)dm));
673:   PetscCall(DMDestroy(&tr->dm));
674:   tr->dm = dm;
675:   PetscFunctionReturn(PETSC_SUCCESS);
676: }

678: /*@
679:   DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform

681:   Input Parameter:
682: . tr - The `DMPlexTransform` object

684:   Output Parameter:
685: . active - The `DMLabel` indicating which points will be transformed

687:   Level: intermediate

689: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
690: @*/
691: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
692: {
693:   PetscFunctionBegin;
695:   PetscAssertPointer(active, 2);
696:   *active = tr->active;
697:   PetscFunctionReturn(PETSC_SUCCESS);
698: }

700: /*@
701:   DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform

703:   Input Parameters:
704: + tr     - The `DMPlexTransform` object
705: - active - The original `DM` which will be transformed

707:   Level: intermediate

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

712: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
713: @*/
714: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
715: {
716:   PetscFunctionBegin;
719:   PetscCall(PetscObjectReference((PetscObject)active));
720:   PetscCall(DMLabelDestroy(&tr->active));
721:   tr->active = active;
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

725: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
726: {
727:   PetscFunctionBegin;
728:   if (!tr->coordFE[ct]) {
729:     PetscInt dim, cdim;

731:     dim = DMPolytopeTypeGetDim(ct);
732:     PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
733:     PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
734:     {
735:       PetscDualSpace  dsp;
736:       PetscQuadrature quad;
737:       DM              K;
738:       PetscFEGeom    *cg;
739:       PetscScalar    *Xq;
740:       PetscReal      *xq, *wq;
741:       PetscInt        Nq, q;

743:       PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
744:       PetscCall(PetscMalloc1(Nq * cdim, &xq));
745:       for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
746:       PetscCall(PetscMalloc1(Nq, &wq));
747:       for (q = 0; q < Nq; ++q) wq[q] = 1.0;
748:       PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
749:       PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
750:       PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));

752:       PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
753:       PetscCall(PetscDualSpaceGetDM(dsp, &K));
754:       PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]));
755:       cg = tr->refGeom[ct];
756:       PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
757:       PetscCall(PetscQuadratureDestroy(&quad));
758:     }
759:   }
760:   *fe = tr->coordFE[ct];
761:   PetscFunctionReturn(PETSC_SUCCESS);
762: }

764: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
765: {
766:   PetscInt dim, cdim;

768:   PetscFunctionBegin;
769:   PetscCall(DMGetDimension(dm, &dim));
770:   PetscCall(DMSetDimension(tdm, dim));
771:   PetscCall(DMGetCoordinateDim(dm, &cdim));
772:   PetscCall(DMSetCoordinateDim(tdm, cdim));
773:   PetscFunctionReturn(PETSC_SUCCESS);
774: }

776: /*@
777:   DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`

779:   Input Parameters:
780: + tr - The `DMPlexTransform` object
781: - dm - The original `DM`

783:   Output Parameter:
784: . tdm - The transformed `DM`

786:   Level: advanced

788: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
789: @*/
790: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
791: {
792:   PetscFunctionBegin;
793:   PetscUseTypeMethod(tr, setdimensions, dm, tdm);
794:   PetscFunctionReturn(PETSC_SUCCESS);
795: }

797: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
798: {
799:   PetscFunctionBegin;
800:   if (pStart) *pStart = 0;
801:   if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
802:   PetscFunctionReturn(PETSC_SUCCESS);
803: }

805: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
806: {
807:   PetscInt ctNew;

809:   PetscFunctionBegin;
811:   PetscAssertPointer(celltype, 3);
812:   /* TODO Can do bisection since everything is sorted */
813:   for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
814:     PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];

816:     if (cell >= ctSN && cell < ctEN) break;
817:   }
818:   PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
819:   *celltype = (DMPolytopeType)ctNew;
820:   PetscFunctionReturn(PETSC_SUCCESS);
821: }

823: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
824: {
825:   PetscFunctionBegin;
827:   if (start) *start = tr->ctStartNew[celltype];
828:   if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
829:   PetscFunctionReturn(PETSC_SUCCESS);
830: }

832: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
833: {
834:   PetscFunctionBegin;
836:   *depth = tr->depth;
837:   PetscFunctionReturn(PETSC_SUCCESS);
838: }

840: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
841: {
842:   PetscFunctionBegin;
844:   if (start) *start = tr->depthStart[depth];
845:   if (end) *end = tr->depthEnd[depth];
846:   PetscFunctionReturn(PETSC_SUCCESS);
847: }

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

852:   Not Collective

854:   Input Parameters:
855: + tr    - The `DMPlexTransform`
856: . ct    - The type of the original point which produces the new point
857: . ctNew - The type of the new point
858: . p     - The original point which produces the new point
859: - r     - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`

861:   Output Parameter:
862: . pNew - The new point number

864:   Level: developer

866: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
867: @*/
868: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
869: {
870:   DMPolytopeType *rct;
871:   PetscInt       *rsize, *cone, *ornt;
872:   PetscInt        rt, Nct, n, off, rp;
873:   DMLabel         trType = tr->trType;
874:   PetscInt        ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
875:   PetscInt        ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
876:   PetscInt        newp = ctSN, cind;

878:   PetscFunctionBeginHot;
879:   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);
880:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
881:   if (trType) {
882:     PetscCall(DMLabelGetValueIndex(trType, rt, &cind));
883:     PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
884:     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);
885:   } else {
886:     cind = ct;
887:     rp   = p - ctS;
888:   }
889:   off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
890:   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);
891:   newp += off;
892:   for (n = 0; n < Nct; ++n) {
893:     if (rct[n] == ctNew) {
894:       if (rsize[n] && r >= rsize[n])
895:         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]);
896:       newp += rp * rsize[n] + r;
897:       break;
898:     }
899:   }

901:   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);
902:   *pNew = newp;
903:   PetscFunctionReturn(PETSC_SUCCESS);
904: }

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

909:   Not Collective

911:   Input Parameters:
912: + tr   - The `DMPlexTransform`
913: - pNew - The new point number

915:   Output Parameters:
916: + ct    - The type of the original point which produces the new point
917: . ctNew - The type of the new point
918: . p     - The original point which produces the new point
919: - r     - The replica number of the new point, meaning it is the rth point of type ctNew produced from p

921:   Level: developer

923: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
924: @*/
925: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
926: {
927:   DMLabel         trType = tr->trType;
928:   DMPolytopeType *rct, ctN;
929:   PetscInt       *rsize, *cone, *ornt;
930:   PetscInt        rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
931:   PetscInt        offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;

933:   PetscFunctionBegin;
934:   PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
935:   if (trType) {
936:     DM              dm;
937:     IS              rtIS;
938:     const PetscInt *reftypes;
939:     PetscInt        Nrt, r, rtStart;

941:     PetscCall(DMPlexTransformGetDM(tr, &dm));
942:     PetscCall(DMLabelGetNumValues(trType, &Nrt));
943:     PetscCall(DMLabelGetValueIS(trType, &rtIS));
944:     PetscCall(ISGetIndices(rtIS, &reftypes));
945:     for (r = 0; r < Nrt; ++r) {
946:       const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];

948:       if (tr->ctStartNew[ctN] + off > pNew) continue;
949:       /* Check that any of this refinement type exist */
950:       /* TODO Actually keep track of the number produced here instead */
951:       if (off > offset) {
952:         rt     = reftypes[r];
953:         offset = off;
954:       }
955:     }
956:     PetscCall(ISRestoreIndices(rtIS, &reftypes));
957:     PetscCall(ISDestroy(&rtIS));
958:     PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
959:     /* TODO Map refinement types to cell types */
960:     PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL));
961:     PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
962:     for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
963:       PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];

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

972:       if (tr->ctStartNew[ctN] + off > pNew) continue;
973:       if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
974:       /* TODO Actually keep track of the number produced here instead */
975:       if (off > offset) {
976:         ctO    = ctTmp;
977:         offset = off;
978:       }
979:     }
980:     rt = -1;
981:     PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
982:   }
983:   ctS = tr->ctStart[ctO];
984:   ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
985:   if (trType) {
986:     for (rtS = ctS; rtS < ctE; ++rtS) {
987:       PetscInt val;
988:       PetscCall(DMLabelGetValue(trType, rtS, &val));
989:       if (val == rt) break;
990:     }
991:     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);
992:   } else rtS = ctS;
993:   PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
994:   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);
995:   for (n = 0; n < Nct; ++n) {
996:     if (rct[n] == ctN) {
997:       PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c;

999:       if (trType) {
1000:         for (c = ctS; c < ctE; ++c) {
1001:           PetscCall(DMLabelGetValue(trType, c, &val));
1002:           if (val == rt) {
1003:             if (tmp < rsize[n]) break;
1004:             tmp -= rsize[n];
1005:           }
1006:         }
1007:         PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew);
1008:         rp = c - ctS;
1009:         rO = tmp;
1010:       } else {
1011:         // This assumes that all points of type ctO transform the same way
1012:         rp = tmp / rsize[n];
1013:         rO = tmp % rsize[n];
1014:       }
1015:       break;
1016:     }
1017:   }
1018:   PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
1019:   pO = rp + ctS;
1020:   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);
1021:   if (ct) *ct = (DMPolytopeType)ctO;
1022:   if (ctNew) *ctNew = (DMPolytopeType)ctN;
1023:   if (p) *p = pO;
1024:   if (r) *r = rO;
1025:   PetscFunctionReturn(PETSC_SUCCESS);
1026: }

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

1031:   Input Parameters:
1032: + tr     - The `DMPlexTransform` object
1033: . source - The source cell type
1034: - p      - The source point, which can also determine the refine type

1036:   Output Parameters:
1037: + rt     - The refine type for this point
1038: . Nt     - The number of types produced by this point
1039: . target - An array of length `Nt` giving the types produced
1040: . size   - An array of length `Nt` giving the number of cells of each type produced
1041: . cone   - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
1042: - ornt   - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point

1044:   Level: advanced

1046:   Notes:
1047:   The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
1048:   need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
1049:   division (described in these outputs) of the cell in the original mesh. The point identifier is given by
1050: .vb
1051:    the number of cones to be taken, or 0 for the current cell
1052:    the cell cone point number at each level from which it is subdivided
1053:    the replica number r of the subdivision.
1054: .ve
1055:   The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
1056: .vb
1057:    Nt     = 2
1058:    target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
1059:    size   = {1, 2}
1060:    cone   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,  DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
1061:    ornt   = {                         0,                       0,                        0,                          0}
1062: .ve

1064: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
1065: @*/
1066: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1067: {
1068:   PetscFunctionBegin;
1069:   PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
1070:   PetscFunctionReturn(PETSC_SUCCESS);
1071: }

1073: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1074: {
1075:   PetscFunctionBegin;
1076:   *rnew = r;
1077:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
1078:   PetscFunctionReturn(PETSC_SUCCESS);
1079: }

1081: /* Returns the same thing */
1082: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1083: {
1084:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
1085:   static PetscInt       vertexS[] = {1};
1086:   static PetscInt       vertexC[] = {0};
1087:   static PetscInt       vertexO[] = {0};
1088:   static DMPolytopeType edgeT[]   = {DM_POLYTOPE_SEGMENT};
1089:   static PetscInt       edgeS[]   = {1};
1090:   static PetscInt       edgeC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1091:   static PetscInt       edgeO[]   = {0, 0};
1092:   static DMPolytopeType tedgeT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR};
1093:   static PetscInt       tedgeS[]  = {1};
1094:   static PetscInt       tedgeC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1095:   static PetscInt       tedgeO[]  = {0, 0};
1096:   static DMPolytopeType triT[]    = {DM_POLYTOPE_TRIANGLE};
1097:   static PetscInt       triS[]    = {1};
1098:   static PetscInt       triC[]    = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1099:   static PetscInt       triO[]    = {0, 0, 0};
1100:   static DMPolytopeType quadT[]   = {DM_POLYTOPE_QUADRILATERAL};
1101:   static PetscInt       quadS[]   = {1};
1102:   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};
1103:   static PetscInt       quadO[]   = {0, 0, 0, 0};
1104:   static DMPolytopeType tquadT[]  = {DM_POLYTOPE_SEG_PRISM_TENSOR};
1105:   static PetscInt       tquadS[]  = {1};
1106:   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};
1107:   static PetscInt       tquadO[]  = {0, 0, 0, 0};
1108:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_TETRAHEDRON};
1109:   static PetscInt       tetS[]    = {1};
1110:   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};
1111:   static PetscInt       tetO[]    = {0, 0, 0, 0};
1112:   static DMPolytopeType hexT[]    = {DM_POLYTOPE_HEXAHEDRON};
1113:   static PetscInt       hexS[]    = {1};
1114:   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};
1115:   static PetscInt       hexO[] = {0, 0, 0, 0, 0, 0};
1116:   static DMPolytopeType tripT[]   = {DM_POLYTOPE_TRI_PRISM};
1117:   static PetscInt       tripS[]   = {1};
1118:   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};
1119:   static PetscInt       tripO[]   = {0, 0, 0, 0, 0};
1120:   static DMPolytopeType ttripT[]  = {DM_POLYTOPE_TRI_PRISM_TENSOR};
1121:   static PetscInt       ttripS[]  = {1};
1122:   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};
1123:   static PetscInt       ttripO[]  = {0, 0, 0, 0, 0};
1124:   static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
1125:   static PetscInt       tquadpS[] = {1};
1126:   static PetscInt       tquadpC[] = {DM_POLYTOPE_QUADRILATERAL,    1, 0, 0, DM_POLYTOPE_QUADRILATERAL,    1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
1127:                                      DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1128:   static PetscInt       tquadpO[] = {0, 0, 0, 0, 0, 0};
1129:   static DMPolytopeType pyrT[]    = {DM_POLYTOPE_PYRAMID};
1130:   static PetscInt       pyrS[]    = {1};
1131:   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};
1132:   static PetscInt       pyrO[]    = {0, 0, 0, 0, 0};

1134:   PetscFunctionBegin;
1135:   if (rt) *rt = 0;
1136:   switch (source) {
1137:   case DM_POLYTOPE_POINT:
1138:     *Nt     = 1;
1139:     *target = vertexT;
1140:     *size   = vertexS;
1141:     *cone   = vertexC;
1142:     *ornt   = vertexO;
1143:     break;
1144:   case DM_POLYTOPE_SEGMENT:
1145:     *Nt     = 1;
1146:     *target = edgeT;
1147:     *size   = edgeS;
1148:     *cone   = edgeC;
1149:     *ornt   = edgeO;
1150:     break;
1151:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1152:     *Nt     = 1;
1153:     *target = tedgeT;
1154:     *size   = tedgeS;
1155:     *cone   = tedgeC;
1156:     *ornt   = tedgeO;
1157:     break;
1158:   case DM_POLYTOPE_TRIANGLE:
1159:     *Nt     = 1;
1160:     *target = triT;
1161:     *size   = triS;
1162:     *cone   = triC;
1163:     *ornt   = triO;
1164:     break;
1165:   case DM_POLYTOPE_QUADRILATERAL:
1166:     *Nt     = 1;
1167:     *target = quadT;
1168:     *size   = quadS;
1169:     *cone   = quadC;
1170:     *ornt   = quadO;
1171:     break;
1172:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1173:     *Nt     = 1;
1174:     *target = tquadT;
1175:     *size   = tquadS;
1176:     *cone   = tquadC;
1177:     *ornt   = tquadO;
1178:     break;
1179:   case DM_POLYTOPE_TETRAHEDRON:
1180:     *Nt     = 1;
1181:     *target = tetT;
1182:     *size   = tetS;
1183:     *cone   = tetC;
1184:     *ornt   = tetO;
1185:     break;
1186:   case DM_POLYTOPE_HEXAHEDRON:
1187:     *Nt     = 1;
1188:     *target = hexT;
1189:     *size   = hexS;
1190:     *cone   = hexC;
1191:     *ornt   = hexO;
1192:     break;
1193:   case DM_POLYTOPE_TRI_PRISM:
1194:     *Nt     = 1;
1195:     *target = tripT;
1196:     *size   = tripS;
1197:     *cone   = tripC;
1198:     *ornt   = tripO;
1199:     break;
1200:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
1201:     *Nt     = 1;
1202:     *target = ttripT;
1203:     *size   = ttripS;
1204:     *cone   = ttripC;
1205:     *ornt   = ttripO;
1206:     break;
1207:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1208:     *Nt     = 1;
1209:     *target = tquadpT;
1210:     *size   = tquadpS;
1211:     *cone   = tquadpC;
1212:     *ornt   = tquadpO;
1213:     break;
1214:   case DM_POLYTOPE_PYRAMID:
1215:     *Nt     = 1;
1216:     *target = pyrT;
1217:     *size   = pyrS;
1218:     *cone   = pyrC;
1219:     *ornt   = pyrO;
1220:     break;
1221:   default:
1222:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1223:   }
1224:   PetscFunctionReturn(PETSC_SUCCESS);
1225: }

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

1230:   Not Collective

1232:   Input Parameters:
1233: + tr  - The `DMPlexTransform`
1234: . sct - The source point cell type, from whom the new cell is being produced
1235: . sp  - The source point
1236: . so  - The orientation of the source point in its enclosing parent
1237: . tct - The target point cell type
1238: . r   - The replica number requested for the produced cell type
1239: - o   - The orientation of the replica

1241:   Output Parameters:
1242: + rnew - The replica number, given the orientation of the parent
1243: - onew - The replica orientation, given the orientation of the parent

1245:   Level: advanced

1247: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1248: @*/
1249: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1250: {
1251:   PetscFunctionBeginHot;
1252:   PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1253:   PetscFunctionReturn(PETSC_SUCCESS);
1254: }

1256: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1257: {
1258:   DM       dm;
1259:   PetscInt pStart, pEnd, p, pNew;

1261:   PetscFunctionBegin;
1262:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1263:   /* Must create the celltype label here so that we do not automatically try to compute the types */
1264:   PetscCall(DMCreateLabel(rdm, "celltype"));
1265:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1266:   for (p = pStart; p < pEnd; ++p) {
1267:     DMPolytopeType  ct;
1268:     DMPolytopeType *rct;
1269:     PetscInt       *rsize, *rcone, *rornt;
1270:     PetscInt        Nct, n, r;

1272:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1273:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1274:     for (n = 0; n < Nct; ++n) {
1275:       for (r = 0; r < rsize[n]; ++r) {
1276:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1277:         PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
1278:         PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
1279:       }
1280:     }
1281:   }
1282:   /* Let the DM know we have set all the cell types */
1283:   {
1284:     DMLabel  ctLabel;
1285:     DM_Plex *plex = (DM_Plex *)rdm->data;

1287:     PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
1288:     PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
1289:   }
1290:   PetscFunctionReturn(PETSC_SUCCESS);
1291: }

1293: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1294: {
1295:   DMPolytopeType ctNew;

1297:   PetscFunctionBegin;
1299:   PetscAssertPointer(coneSize, 3);
1300:   PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
1301:   *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType)ctNew);
1302:   PetscFunctionReturn(PETSC_SUCCESS);
1303: }

1305: /* The orientation o is for the interior of the cell p */
1306: 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[])
1307: {
1308:   DM              dm;
1309:   const PetscInt  csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1310:   const PetscInt *cone;
1311:   PetscInt        c, coff = *coneoff, ooff = *orntoff;

1313:   PetscFunctionBegin;
1314:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1315:   PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
1316:   for (c = 0; c < csizeNew; ++c) {
1317:     PetscInt             ppp   = -1;                            /* Parent Parent point: Parent of point pp */
1318:     PetscInt             pp    = p;                             /* Parent point: Point in the original mesh producing new cone point */
1319:     PetscInt             po    = 0;                             /* Orientation of parent point pp in parent parent point ppp */
1320:     DMPolytopeType       pct   = ct;                            /* Parent type: Cell type for parent of new cone point */
1321:     const PetscInt      *pcone = cone;                          /* Parent cone: Cone of parent point pp */
1322:     PetscInt             pr    = -1;                            /* Replica number of pp that produces new cone point  */
1323:     const DMPolytopeType ft    = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1324:     const PetscInt       fn    = rcone[coff++];                 /* Number of cones of p that need to be taken when producing new cone point */
1325:     PetscInt             fo    = rornt[ooff++];                 /* Orientation of new cone point in pNew */
1326:     PetscInt             lc;

1328:     /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1329:     for (lc = 0; lc < fn; ++lc) {
1330:       const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po);
1331:       const PetscInt  acp  = rcone[coff++];
1332:       const PetscInt  pcp  = parr[acp * 2];
1333:       const PetscInt  pco  = parr[acp * 2 + 1];
1334:       const PetscInt *ppornt;

1336:       ppp = pp;
1337:       pp  = pcone[pcp];
1338:       PetscCall(DMPlexGetCellType(dm, pp, &pct));
1339:       // Restore the parent cone from the last iterate
1340:       if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
1341:       PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
1342:       PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
1343:       po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1344:       PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
1345:     }
1346:     if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
1347:     pr = rcone[coff++];
1348:     /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1349:     PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
1350:     PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
1351:     orntNew[c] = fo;
1352:   }
1353:   PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
1354:   *coneoff = coff;
1355:   *orntoff = ooff;
1356:   PetscFunctionReturn(PETSC_SUCCESS);
1357: }

1359: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1360: {
1361:   DM             dm;
1362:   DMPolytopeType ct;
1363:   PetscInt      *coneNew, *orntNew;
1364:   PetscInt       maxConeSize = 0, pStart, pEnd, p, pNew;

1366:   PetscFunctionBegin;
1367:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1368:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1369:   PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1370:   PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1371:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1372:   for (p = pStart; p < pEnd; ++p) {
1373:     PetscInt        coff, ooff;
1374:     DMPolytopeType *rct;
1375:     PetscInt       *rsize, *rcone, *rornt;
1376:     PetscInt        Nct, n, r;

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

1383:       for (r = 0; r < rsize[n]; ++r) {
1384:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1385:         PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
1386:         PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
1387:         PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
1388:       }
1389:     }
1390:   }
1391:   PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1392:   PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1393:   PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
1394:   PetscCall(DMPlexSymmetrize(rdm));
1395:   PetscCall(DMPlexStratify(rdm));
1396:   PetscFunctionReturn(PETSC_SUCCESS);
1397: }

1399: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1400: {
1401:   DM              dm;
1402:   DMPolytopeType  ct, qct;
1403:   DMPolytopeType *rct;
1404:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
1405:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

1407:   PetscFunctionBegin;
1409:   PetscAssertPointer(cone, 4);
1410:   PetscAssertPointer(ornt, 5);
1411:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1412:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1413:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1414:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1415:   PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1416:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1417:   for (n = 0; n < Nct; ++n) {
1418:     const DMPolytopeType ctNew    = rct[n];
1419:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1420:     PetscInt             Nr       = rsize[n], fn, c;

1422:     if (ctNew == qct) Nr = r;
1423:     for (nr = 0; nr < Nr; ++nr) {
1424:       for (c = 0; c < csizeNew; ++c) {
1425:         ++coff;             /* Cell type of new cone point */
1426:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1427:         coff += fn;
1428:         ++coff; /* Replica number of new cone point */
1429:         ++ooff; /* Orientation of new cone point */
1430:       }
1431:     }
1432:     if (ctNew == qct) break;
1433:   }
1434:   PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1435:   *cone = qcone;
1436:   *ornt = qornt;
1437:   PetscFunctionReturn(PETSC_SUCCESS);
1438: }

1440: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1441: {
1442:   DM              dm;
1443:   DMPolytopeType  ct, qct;
1444:   DMPolytopeType *rct;
1445:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
1446:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

1448:   PetscFunctionBegin;
1450:   if (cone) PetscAssertPointer(cone, 3);
1451:   if (ornt) PetscAssertPointer(ornt, 4);
1452:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1453:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1454:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1455:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1456:   PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1457:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1458:   for (n = 0; n < Nct; ++n) {
1459:     const DMPolytopeType ctNew    = rct[n];
1460:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1461:     PetscInt             Nr       = rsize[n], fn, c;

1463:     if (ctNew == qct) Nr = r;
1464:     for (nr = 0; nr < Nr; ++nr) {
1465:       for (c = 0; c < csizeNew; ++c) {
1466:         ++coff;             /* Cell type of new cone point */
1467:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1468:         coff += fn;
1469:         ++coff; /* Replica number of new cone point */
1470:         ++ooff; /* Orientation of new cone point */
1471:       }
1472:     }
1473:     if (ctNew == qct) break;
1474:   }
1475:   PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1476:   if (cone) *cone = qcone;
1477:   else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1478:   if (ornt) *ornt = qornt;
1479:   else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1480:   PetscFunctionReturn(PETSC_SUCCESS);
1481: }

1483: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1484: {
1485:   DM dm;

1487:   PetscFunctionBegin;
1489:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1490:   if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
1491:   if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
1492:   PetscFunctionReturn(PETSC_SUCCESS);
1493: }

1495: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1496: {
1497:   PetscInt ict;

1499:   PetscFunctionBegin;
1500:   PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
1501:   for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1502:     const DMPolytopeType ct = (DMPolytopeType)ict;
1503:     DMPlexTransform      reftr;
1504:     DM                   refdm, trdm;
1505:     Vec                  coordinates;
1506:     const PetscScalar   *coords;
1507:     DMPolytopeType      *rct;
1508:     PetscInt            *rsize, *rcone, *rornt;
1509:     PetscInt             Nct, n, r, pNew;
1510:     PetscInt             trdim, vStart, vEnd, Nc;
1511:     const PetscInt       debug = 0;
1512:     const char          *typeName;

1514:     /* Since points are 0-dimensional, coordinates make no sense */
1515:     if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue;
1516:     PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
1517:     PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
1518:     PetscCall(DMPlexTransformSetDM(reftr, refdm));
1519:     PetscCall(DMPlexTransformGetType(tr, &typeName));
1520:     PetscCall(DMPlexTransformSetType(reftr, typeName));
1521:     PetscCall(DMPlexTransformSetUp(reftr));
1522:     PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));

1524:     PetscCall(DMGetDimension(trdm, &trdim));
1525:     PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
1526:     tr->trNv[ct] = vEnd - vStart;
1527:     PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
1528:     PetscCall(VecGetLocalSize(coordinates, &Nc));
1529:     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);
1530:     PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
1531:     PetscCall(VecGetArrayRead(coordinates, &coords));
1532:     PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
1533:     PetscCall(VecRestoreArrayRead(coordinates, &coords));

1535:     PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
1536:     PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1537:     for (n = 0; n < Nct; ++n) {
1538:       /* Since points are 0-dimensional, coordinates make no sense */
1539:       if (rct[n] == DM_POLYTOPE_POINT) continue;
1540:       PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
1541:       for (r = 0; r < rsize[n]; ++r) {
1542:         PetscInt *closure = NULL;
1543:         PetscInt  clSize, cl, Nv = 0;

1545:         PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
1546:         PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
1547:         PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1548:         for (cl = 0; cl < clSize * 2; cl += 2) {
1549:           const PetscInt sv = closure[cl];

1551:           if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1552:         }
1553:         PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1554:         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]);
1555:       }
1556:     }
1557:     if (debug) {
1558:       DMPolytopeType *rct;
1559:       PetscInt       *rsize, *rcone, *rornt;
1560:       PetscInt        v, dE = trdim, d, off = 0;

1562:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
1563:       for (v = 0; v < tr->trNv[ct]; ++v) {
1564:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  "));
1565:         for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
1566:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1567:       }

1569:       PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1570:       for (n = 0; n < Nct; ++n) {
1571:         if (rct[n] == DM_POLYTOPE_POINT) continue;
1572:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
1573:         for (r = 0; r < rsize[n]; ++r) {
1574:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  "));
1575:           for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
1576:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1577:         }
1578:       }
1579:     }
1580:     PetscCall(DMDestroy(&refdm));
1581:     PetscCall(DMDestroy(&trdm));
1582:     PetscCall(DMPlexTransformDestroy(&reftr));
1583:   }
1584:   PetscFunctionReturn(PETSC_SUCCESS);
1585: }

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

1590:   Input Parameters:
1591: + tr - The `DMPlexTransform` object
1592: - ct - The cell type

1594:   Output Parameters:
1595: + Nv      - The number of transformed vertices in the closure of the reference cell of given type
1596: - trVerts - The coordinates of these vertices in the reference cell

1598:   Level: developer

1600: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
1601: @*/
1602: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1603: {
1604:   PetscFunctionBegin;
1605:   if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1606:   if (Nv) *Nv = tr->trNv[ct];
1607:   if (trVerts) *trVerts = tr->trVerts[ct];
1608:   PetscFunctionReturn(PETSC_SUCCESS);
1609: }

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

1614:   Input Parameters:
1615: + tr  - The `DMPlexTransform` object
1616: . ct  - The cell type
1617: . rct - The subcell type
1618: - r   - The subcell index

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

1623:   Level: developer

1625: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
1626: @*/
1627: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1628: {
1629:   PetscFunctionBegin;
1630:   if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1631:   PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
1632:   if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1633:   PetscFunctionReturn(PETSC_SUCCESS);
1634: }

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

1641:   PetscFunctionBeginHot;
1642:   PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
1643:   for (d = 0; d < dE; ++d) out[d] = 0.0;
1644:   for (v = 0; v < Nv; ++v)
1645:     for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1646:   for (d = 0; d < dE; ++d) out[d] /= Nv;
1647:   PetscFunctionReturn(PETSC_SUCCESS);
1648: }

1650: /*@
1651:   DMPlexTransformMapCoordinates - Calculate new coordinates for produced points

1653:   Not collective

1655:   Input Parameters:
1656: + tr  - The `DMPlexTransform`
1657: . pct - The cell type of the parent, from whom the new cell is being produced
1658: . ct  - The type being produced
1659: . p   - The original point
1660: . r   - The replica number requested for the produced cell type
1661: . Nv  - Number of vertices in the closure of the parent cell
1662: . dE  - Spatial dimension
1663: - in  - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell

1665:   Output Parameter:
1666: . out - The coordinates of the new vertices

1668:   Level: intermediate

1670: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
1671: @*/
1672: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1673: {
1674:   PetscFunctionBeginHot;
1675:   if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1676:   PetscFunctionReturn(PETSC_SUCCESS);
1677: }

1679: /*
1680:   DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label

1682:   Not Collective

1684:   Input Parameters:
1685: + tr    - The `DMPlexTransform`
1686: . label - The label in the transformed mesh
1687: . pp    - The parent point in the original mesh
1688: . pct   - The cell type of the parent point
1689: . p     - The point in the transformed mesh
1690: . ct    - The cell type of the point
1691: . r     - The replica number of the point
1692: - val   - The label value of the parent point

1694:   Level: developer

1696: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
1697: */
1698: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
1699: {
1700:   PetscFunctionBeginHot;
1701:   if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
1702:   PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
1703:   PetscFunctionReturn(PETSC_SUCCESS);
1704: }

1706: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1707: {
1708:   DM              dm;
1709:   IS              valueIS;
1710:   const PetscInt *values;
1711:   PetscInt        defVal, Nv, val;

1713:   PetscFunctionBegin;
1714:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1715:   PetscCall(DMLabelGetDefaultValue(label, &defVal));
1716:   PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
1717:   PetscCall(DMLabelGetValueIS(label, &valueIS));
1718:   PetscCall(ISGetLocalSize(valueIS, &Nv));
1719:   PetscCall(ISGetIndices(valueIS, &values));
1720:   for (val = 0; val < Nv; ++val) {
1721:     IS              pointIS;
1722:     const PetscInt *points;
1723:     PetscInt        numPoints, p;

1725:     /* Ensure refined label is created with same number of strata as
1726:      * original (even if no entries here). */
1727:     PetscCall(DMLabelAddStratum(labelNew, values[val]));
1728:     PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
1729:     PetscCall(ISGetLocalSize(pointIS, &numPoints));
1730:     PetscCall(ISGetIndices(pointIS, &points));
1731:     for (p = 0; p < numPoints; ++p) {
1732:       const PetscInt  point = points[p];
1733:       DMPolytopeType  ct;
1734:       DMPolytopeType *rct;
1735:       PetscInt       *rsize, *rcone, *rornt;
1736:       PetscInt        Nct, n, r, pNew = 0;

1738:       PetscCall(DMPlexGetCellType(dm, point, &ct));
1739:       PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1740:       for (n = 0; n < Nct; ++n) {
1741:         for (r = 0; r < rsize[n]; ++r) {
1742:           PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
1743:           PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
1744:         }
1745:       }
1746:     }
1747:     PetscCall(ISRestoreIndices(pointIS, &points));
1748:     PetscCall(ISDestroy(&pointIS));
1749:   }
1750:   PetscCall(ISRestoreIndices(valueIS, &values));
1751:   PetscCall(ISDestroy(&valueIS));
1752:   PetscFunctionReturn(PETSC_SUCCESS);
1753: }

1755: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1756: {
1757:   DM       dm;
1758:   PetscInt numLabels, l;

1760:   PetscFunctionBegin;
1761:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1762:   PetscCall(DMGetNumLabels(dm, &numLabels));
1763:   for (l = 0; l < numLabels; ++l) {
1764:     DMLabel     label, labelNew;
1765:     const char *lname;
1766:     PetscBool   isDepth, isCellType;

1768:     PetscCall(DMGetLabelName(dm, l, &lname));
1769:     PetscCall(PetscStrcmp(lname, "depth", &isDepth));
1770:     if (isDepth) continue;
1771:     PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
1772:     if (isCellType) continue;
1773:     PetscCall(DMCreateLabel(rdm, lname));
1774:     PetscCall(DMGetLabel(dm, lname, &label));
1775:     PetscCall(DMGetLabel(rdm, lname, &labelNew));
1776:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
1777:   }
1778:   PetscFunctionReturn(PETSC_SUCCESS);
1779: }

1781: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1782: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1783: {
1784:   DM       dm;
1785:   PetscInt Nf, f, Nds, s;

1787:   PetscFunctionBegin;
1788:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1789:   PetscCall(DMGetNumFields(dm, &Nf));
1790:   for (f = 0; f < Nf; ++f) {
1791:     DMLabel     label, labelNew;
1792:     PetscObject obj;
1793:     const char *lname;

1795:     PetscCall(DMGetField(rdm, f, &label, &obj));
1796:     if (!label) continue;
1797:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1798:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1799:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
1800:     PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
1801:     PetscCall(DMLabelDestroy(&labelNew));
1802:   }
1803:   PetscCall(DMGetNumDS(dm, &Nds));
1804:   for (s = 0; s < Nds; ++s) {
1805:     DMLabel     label, labelNew;
1806:     const char *lname;

1808:     PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
1809:     if (!label) continue;
1810:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1811:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1812:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
1813:     PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
1814:     PetscCall(DMLabelDestroy(&labelNew));
1815:   }
1816:   PetscFunctionReturn(PETSC_SUCCESS);
1817: }

1819: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1820: {
1821:   DM                 dm;
1822:   PetscSF            sf, sfNew;
1823:   PetscInt           numRoots, numLeaves, numLeavesNew = 0, l, m;
1824:   const PetscInt    *localPoints;
1825:   const PetscSFNode *remotePoints;
1826:   PetscInt          *localPointsNew;
1827:   PetscSFNode       *remotePointsNew;
1828:   PetscInt           pStartNew, pEndNew, pNew;
1829:   /* Brute force algorithm */
1830:   PetscSF         rsf;
1831:   PetscSection    s;
1832:   const PetscInt *rootdegree;
1833:   PetscInt       *rootPointsNew, *remoteOffsets;
1834:   PetscInt        numPointsNew, pStart, pEnd, p;

1836:   PetscFunctionBegin;
1837:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1838:   PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
1839:   PetscCall(DMGetPointSF(dm, &sf));
1840:   PetscCall(DMGetPointSF(rdm, &sfNew));
1841:   /* Calculate size of new SF */
1842:   PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
1843:   if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
1844:   for (l = 0; l < numLeaves; ++l) {
1845:     const PetscInt  p = localPoints[l];
1846:     DMPolytopeType  ct;
1847:     DMPolytopeType *rct;
1848:     PetscInt       *rsize, *rcone, *rornt;
1849:     PetscInt        Nct, n;

1851:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1852:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1853:     for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
1854:   }
1855:   /* Send new root point numbers
1856:        It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1857:   */
1858:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1859:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
1860:   PetscCall(PetscSectionSetChart(s, pStart, pEnd));
1861:   for (p = pStart; p < pEnd; ++p) {
1862:     DMPolytopeType  ct;
1863:     DMPolytopeType *rct;
1864:     PetscInt       *rsize, *rcone, *rornt;
1865:     PetscInt        Nct, n;

1867:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1868:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1869:     for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
1870:   }
1871:   PetscCall(PetscSectionSetUp(s));
1872:   PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
1873:   PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
1874:   PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
1875:   PetscCall(PetscFree(remoteOffsets));
1876:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1877:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1878:   PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
1879:   for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1880:   for (p = pStart; p < pEnd; ++p) {
1881:     DMPolytopeType  ct;
1882:     DMPolytopeType *rct;
1883:     PetscInt       *rsize, *rcone, *rornt;
1884:     PetscInt        Nct, n, r, off;

1886:     if (!rootdegree[p - pStart]) continue;
1887:     PetscCall(PetscSectionGetOffset(s, p, &off));
1888:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1889:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1890:     for (n = 0, m = 0; n < Nct; ++n) {
1891:       for (r = 0; r < rsize[n]; ++r, ++m) {
1892:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1893:         rootPointsNew[off + m] = pNew;
1894:       }
1895:     }
1896:   }
1897:   PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1898:   PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1899:   PetscCall(PetscSFDestroy(&rsf));
1900:   PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
1901:   PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
1902:   for (l = 0, m = 0; l < numLeaves; ++l) {
1903:     const PetscInt  p = localPoints[l];
1904:     DMPolytopeType  ct;
1905:     DMPolytopeType *rct;
1906:     PetscInt       *rsize, *rcone, *rornt;
1907:     PetscInt        Nct, n, r, q, off;

1909:     PetscCall(PetscSectionGetOffset(s, p, &off));
1910:     PetscCall(DMPlexGetCellType(dm, p, &ct));
1911:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1912:     for (n = 0, q = 0; n < Nct; ++n) {
1913:       for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1914:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1915:         localPointsNew[m]        = pNew;
1916:         remotePointsNew[m].index = rootPointsNew[off + q];
1917:         remotePointsNew[m].rank  = remotePoints[l].rank;
1918:       }
1919:     }
1920:   }
1921:   PetscCall(PetscSectionDestroy(&s));
1922:   PetscCall(PetscFree(rootPointsNew));
1923:   /* SF needs sorted leaves to correctly calculate Gather */
1924:   {
1925:     PetscSFNode *rp, *rtmp;
1926:     PetscInt    *lp, *idx, *ltmp, i;

1928:     PetscCall(PetscMalloc1(numLeavesNew, &idx));
1929:     PetscCall(PetscMalloc1(numLeavesNew, &lp));
1930:     PetscCall(PetscMalloc1(numLeavesNew, &rp));
1931:     for (i = 0; i < numLeavesNew; ++i) {
1932:       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);
1933:       idx[i] = i;
1934:     }
1935:     PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
1936:     for (i = 0; i < numLeavesNew; ++i) {
1937:       lp[i] = localPointsNew[idx[i]];
1938:       rp[i] = remotePointsNew[idx[i]];
1939:     }
1940:     ltmp            = localPointsNew;
1941:     localPointsNew  = lp;
1942:     rtmp            = remotePointsNew;
1943:     remotePointsNew = rp;
1944:     PetscCall(PetscFree(idx));
1945:     PetscCall(PetscFree(ltmp));
1946:     PetscCall(PetscFree(rtmp));
1947:   }
1948:   PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
1949:   PetscFunctionReturn(PETSC_SUCCESS);
1950: }

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

1955:   Not Collective

1957:   Input Parameters:
1958: + tr  - The `DMPlexTransform`
1959: . ct  - The type of the parent cell
1960: . rct - The type of the produced cell
1961: . r   - The index of the produced cell
1962: - x   - The localized coordinates for the parent cell

1964:   Output Parameter:
1965: . xr  - The localized coordinates for the produced cell

1967:   Level: developer

1969: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
1970: */
1971: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1972: {
1973:   PetscFE  fe = NULL;
1974:   PetscInt cdim, v, *subcellV;

1976:   PetscFunctionBegin;
1977:   PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
1978:   PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
1979:   PetscCall(PetscFEGetNumComponents(fe, &cdim));
1980:   for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
1981:   PetscFunctionReturn(PETSC_SUCCESS);
1982: }

1984: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1985: {
1986:   DM                 dm, cdm, cdmCell, cdmNew, cdmCellNew;
1987:   PetscSection       coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
1988:   Vec                coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
1989:   const PetscScalar *coords;
1990:   PetscScalar       *coordsNew;
1991:   const PetscReal   *maxCell, *Lstart, *L;
1992:   PetscBool          localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1993:   PetscInt           dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;

1995:   PetscFunctionBegin;
1996:   PetscCall(DMPlexTransformGetDM(tr, &dm));
1997:   PetscCall(DMGetCoordinateDM(dm, &cdm));
1998:   PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
1999:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
2000:   PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
2001:   if (localized) {
2002:     /* Localize coordinates of new vertices */
2003:     localizeVertices = PETSC_TRUE;
2004:     /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
2005:     if (!maxCell) localizeCells = PETSC_TRUE;
2006:   }
2007:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
2008:   PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
2009:   if (maxCell) {
2010:     PetscReal maxCellNew[3];

2012:     for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
2013:     PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L));
2014:   }
2015:   PetscCall(DMGetCoordinateDim(rdm, &dE));
2016:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
2017:   PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
2018:   PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
2019:   PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
2020:   PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
2021:   /* Localization should be inherited */
2022:   /*   Stefano calculates parent cells for each new cell for localization */
2023:   /*   Localized cells need coordinates of closure */
2024:   for (v = vStartNew; v < vEndNew; ++v) {
2025:     PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
2026:     PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
2027:   }
2028:   PetscCall(PetscSectionSetUp(coordSectionNew));
2029:   PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));

2031:   if (localizeCells) {
2032:     PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
2033:     PetscCall(DMClone(cdmNew, &cdmCellNew));
2034:     PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
2035:     PetscCall(DMDestroy(&cdmCellNew));

2037:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
2038:     PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
2039:     PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
2040:     PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
2041:     PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));

2043:     PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
2044:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2045:     for (c = cStart; c < cEnd; ++c) {
2046:       PetscInt dof;

2048:       PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
2049:       if (dof) {
2050:         DMPolytopeType  ct;
2051:         DMPolytopeType *rct;
2052:         PetscInt       *rsize, *rcone, *rornt;
2053:         PetscInt        dim, cNew, Nct, n, r;

2055:         PetscCall(DMPlexGetCellType(dm, c, &ct));
2056:         dim = DMPolytopeTypeGetDim(ct);
2057:         PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2058:         /* This allows for different cell types */
2059:         for (n = 0; n < Nct; ++n) {
2060:           if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2061:           for (r = 0; r < rsize[n]; ++r) {
2062:             PetscInt *closure = NULL;
2063:             PetscInt  clSize, cl, Nv = 0;

2065:             PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
2066:             PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2067:             for (cl = 0; cl < clSize * 2; cl += 2) {
2068:               if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
2069:             }
2070:             PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2071:             PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
2072:             PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
2073:           }
2074:         }
2075:       }
2076:     }
2077:     PetscCall(PetscSectionSetUp(coordSectionCellNew));
2078:     PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
2079:   }
2080:   PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
2081:   {
2082:     VecType     vtype;
2083:     PetscInt    coordSizeNew, bs;
2084:     const char *name;

2086:     PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
2087:     PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
2088:     PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
2089:     PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
2090:     PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
2091:     PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
2092:     PetscCall(VecGetBlockSize(coordsLocal, &bs));
2093:     PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
2094:     PetscCall(VecGetType(coordsLocal, &vtype));
2095:     PetscCall(VecSetType(coordsLocalNew, vtype));
2096:   }
2097:   PetscCall(VecGetArrayRead(coordsLocal, &coords));
2098:   PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
2099:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2100:   /* First set coordinates for vertices */
2101:   for (p = pStart; p < pEnd; ++p) {
2102:     DMPolytopeType  ct;
2103:     DMPolytopeType *rct;
2104:     PetscInt       *rsize, *rcone, *rornt;
2105:     PetscInt        Nct, n, r;
2106:     PetscBool       hasVertex = PETSC_FALSE;

2108:     PetscCall(DMPlexGetCellType(dm, p, &ct));
2109:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2110:     for (n = 0; n < Nct; ++n) {
2111:       if (rct[n] == DM_POLYTOPE_POINT) {
2112:         hasVertex = PETSC_TRUE;
2113:         break;
2114:       }
2115:     }
2116:     if (hasVertex) {
2117:       const PetscScalar *icoords = NULL;
2118:       const PetscScalar *array   = NULL;
2119:       PetscScalar       *pcoords = NULL;
2120:       PetscBool          isDG;
2121:       PetscInt           Nc, Nv, v, d;

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

2125:       icoords = pcoords;
2126:       Nv      = Nc / dEo;
2127:       if (ct != DM_POLYTOPE_POINT) {
2128:         if (localizeVertices && maxCell) {
2129:           PetscScalar anchor[3];

2131:           for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
2132:           for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
2133:         }
2134:       }
2135:       for (n = 0; n < Nct; ++n) {
2136:         if (rct[n] != DM_POLYTOPE_POINT) continue;
2137:         for (r = 0; r < rsize[n]; ++r) {
2138:           PetscScalar vcoords[3];
2139:           PetscInt    vNew, off;

2141:           PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
2142:           PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
2143:           PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
2144:           PetscCall(DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
2145:         }
2146:       }
2147:       PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2148:     }
2149:   }
2150:   PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
2151:   PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
2152:   PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
2153:   PetscCall(VecDestroy(&coordsLocalNew));
2154:   PetscCall(PetscSectionDestroy(&coordSectionNew));
2155:   /* Then set coordinates for cells by localizing */
2156:   if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
2157:   else {
2158:     VecType     vtype;
2159:     PetscInt    coordSizeNew, bs;
2160:     const char *name;

2162:     PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
2163:     PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
2164:     PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
2165:     PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
2166:     PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
2167:     PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
2168:     PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
2169:     PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
2170:     PetscCall(VecGetType(coordsLocalCell, &vtype));
2171:     PetscCall(VecSetType(coordsLocalCellNew, vtype));
2172:     PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
2173:     PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));

2175:     for (p = pStart; p < pEnd; ++p) {
2176:       DMPolytopeType  ct;
2177:       DMPolytopeType *rct;
2178:       PetscInt       *rsize, *rcone, *rornt;
2179:       PetscInt        dof = 0, Nct, n, r;

2181:       PetscCall(DMPlexGetCellType(dm, p, &ct));
2182:       PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2183:       if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
2184:       if (dof) {
2185:         const PetscScalar *pcoords;

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

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

2195:             /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2196:                DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2197:                cell to the ones it produces. */
2198:             PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2199:             PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
2200:             PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
2201:           }
2202:         }
2203:       }
2204:     }
2205:     PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
2206:     PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
2207:     PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
2208:     PetscCall(VecDestroy(&coordsLocalCellNew));
2209:     PetscCall(PetscSectionDestroy(&coordSectionCellNew));
2210:   }
2211:   PetscFunctionReturn(PETSC_SUCCESS);
2212: }

2214: /*@
2215:   DMPlexTransformApply - Execute the transformation, producing another `DM`

2217:   Collective

2219:   Input Parameters:
2220: + tr - The `DMPlexTransform` object
2221: - dm - The original `DM`

2223:   Output Parameter:
2224: . tdm - The transformed `DM`

2226:   Level: intermediate

2228: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()`
2229: @*/
2230: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
2231: {
2232:   DM                     rdm;
2233:   DMPlexInterpolatedFlag interp;
2234:   PetscInt               pStart, pEnd;

2236:   PetscFunctionBegin;
2239:   PetscAssertPointer(tdm, 3);
2240:   PetscCall(PetscLogEventBegin(DMPLEX_Transform, tr, dm, 0, 0));
2241:   PetscCall(DMPlexTransformSetDM(tr, dm));

2243:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
2244:   PetscCall(DMSetType(rdm, DMPLEX));
2245:   PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
2246:   /* Calculate number of new points of each depth */
2247:   PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
2248:   PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
2249:   /* Step 1: Set chart */
2250:   PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
2251:   PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
2252:   /* Step 2: Set cone/support sizes (automatically stratifies) */
2253:   PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
2254:   /* Step 3: Setup refined DM */
2255:   PetscCall(DMSetUp(rdm));
2256:   /* Step 4: Set cones and supports (automatically symmetrizes) */
2257:   PetscCall(DMPlexTransformSetCones(tr, rdm));
2258:   /* Step 5: Create pointSF */
2259:   PetscCall(DMPlexTransformCreateSF(tr, rdm));
2260:   /* Step 6: Create labels */
2261:   PetscCall(DMPlexTransformCreateLabels(tr, rdm));
2262:   /* Step 7: Set coordinates */
2263:   PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
2264:   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm));
2265:   // If the original DM was configured from options, the transformed DM should be as well
2266:   rdm->setfromoptionscalled = dm->setfromoptionscalled;
2267:   PetscCall(PetscLogEventEnd(DMPLEX_Transform, tr, dm, 0, 0));
2268:   *tdm = rdm;
2269:   PetscFunctionReturn(PETSC_SUCCESS);
2270: }

2272: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2273: {
2274:   DMPlexTransform tr;
2275:   DM              cdm, rcdm;
2276:   const char     *prefix;

2278:   PetscFunctionBegin;
2279:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
2280:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2281:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
2282:   PetscCall(DMPlexTransformSetDM(tr, dm));
2283:   PetscCall(DMPlexTransformSetFromOptions(tr));
2284:   PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
2285:   PetscCall(DMPlexTransformSetUp(tr));
2286:   PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
2287:   PetscCall(DMPlexTransformApply(tr, dm, rdm));
2288:   PetscCall(DMCopyDisc(dm, *rdm));
2289:   PetscCall(DMGetCoordinateDM(dm, &cdm));
2290:   PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
2291:   PetscCall(DMCopyDisc(cdm, rcdm));
2292:   PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
2293:   PetscCall(DMCopyDisc(dm, *rdm));
2294:   PetscCall(DMPlexTransformDestroy(&tr));
2295:   ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2296:   PetscFunctionReturn(PETSC_SUCCESS);
2297: }