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: PetscLogEvent DMPLEXTRANSFORM_SetUp, DMPLEXTRANSFORM_Apply, DMPLEXTRANSFORM_SetConeSizes, DMPLEXTRANSFORM_SetCones, DMPLEXTRANSFORM_CreateSF, DMPLEXTRANSFORM_CreateLabels, DMPLEXTRANSFORM_SetCoordinates;
12: /* Construct cell type order since we must loop over cell types in the same dimensional order they are stored in the plex if dm != NULL
13: OR in standard plex ordering if dm == NULL */
14: static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DM dm, PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
15: {
16: PetscInt *ctO, *ctOInv;
17: PetscInt d, c, off = 0;
18: PetscInt dimOrder[5] = {3, 2, 1, 0, -1};
20: PetscFunctionBegin;
21: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv));
22: if (dm) { // Order the dimensions by their starting location
23: PetscInt hStart[4] = {-1, -1, -1, -1};
24: for (d = 0; d <= dim; ++d) PetscCall(DMPlexGetDepthStratum(dm, dim - d, &hStart[d], NULL));
25: PetscCall(PetscSortIntWithArray(dim + 1, hStart, &dimOrder[3 - dim]));
26: } else if (dim > 1) { // Standard plex ordering. dimOrder is in correct order if dim > 1
27: off = 4 - dim;
28: dimOrder[off++] = 0;
29: for (d = dim - 1; d > 0; --d) dimOrder[off++] = d;
30: }
32: off = 0;
33: for (d = 0; d < 5; ++d) {
34: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
35: if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue;
36: if (DMPolytopeTypeGetDim((DMPolytopeType)c) == dimOrder[d]) ctO[off++] = c;
37: }
38: }
39: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
40: if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) ctO[off++] = c;
41: }
42: ctO[off++] = DM_NUM_POLYTOPES;
43: PetscCheck(off == DM_NUM_POLYTOPES + 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %" PetscInt_FMT " for cell type order", off);
45: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c;
47: *ctOrder = ctO;
48: *ctOrderInv = ctOInv;
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: /*@C
53: DMPlexTransformRegister - Adds a new transform component implementation
55: Not Collective
57: Input Parameters:
58: + name - The name of a new user-defined creation routine
59: - create_func - The creation routine
61: Example Usage:
62: .vb
63: DMPlexTransformRegister("my_transform", MyTransformCreate);
64: .ve
66: Then, your transform type can be chosen with the procedural interface via
67: .vb
68: DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
69: DMPlexTransformSetType(DMPlexTransform, "my_transform");
70: .ve
71: or at runtime via the option
72: .vb
73: -dm_plex_transform_type my_transform
74: .ve
76: Level: advanced
78: Note:
79: `DMPlexTransformRegister()` may be called multiple times to add several user-defined transforms
81: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()`
82: @*/
83: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
84: {
85: PetscFunctionBegin;
86: PetscCall(DMInitializePackage());
87: PetscCall(PetscFunctionListAdd(&DMPlexTransformList, name, create_func));
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
92: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
93: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
94: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToSimplex(DMPlexTransform);
95: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
96: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
97: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
98: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform);
99: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform);
100: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform);
102: /*@C
103: DMPlexTransformRegisterAll - Registers all of the transform components in the `DM` package.
105: Not Collective
107: Level: advanced
109: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
110: @*/
111: PetscErrorCode DMPlexTransformRegisterAll(void)
112: {
113: PetscFunctionBegin;
114: if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
115: DMPlexTransformRegisterAllCalled = PETSC_TRUE;
117: PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter));
118: PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular));
119: PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox));
120: PetscCall(DMPlexTransformRegister(DMPLEXREFINETOSIMPLEX, DMPlexTransformCreate_ToSimplex));
121: PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld));
122: PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL));
123: PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR));
124: PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D));
125: PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDETYPE, DMPlexTransformCreate_Extrude));
126: PetscCall(DMPlexTransformRegister(DMPLEXCOHESIVEEXTRUDE, DMPlexTransformCreate_Cohesive));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: /*@C
131: DMPlexTransformRegisterDestroy - This function destroys the registered `DMPlexTransformType`. It is called from `PetscFinalize()`.
133: Not collective
135: Level: developer
137: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()`
138: @*/
139: PetscErrorCode DMPlexTransformRegisterDestroy(void)
140: {
141: PetscFunctionBegin;
142: PetscCall(PetscFunctionListDestroy(&DMPlexTransformList));
143: DMPlexTransformRegisterAllCalled = PETSC_FALSE;
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /*@
148: DMPlexTransformCreate - Creates an empty transform object. The type can then be set with `DMPlexTransformSetType()`.
150: Collective
152: Input Parameter:
153: . comm - The communicator for the transform object
155: Output Parameter:
156: . tr - The transform object
158: Level: beginner
160: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
161: @*/
162: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
163: {
164: DMPlexTransform t;
166: PetscFunctionBegin;
167: PetscAssertPointer(tr, 2);
168: *tr = NULL;
169: PetscCall(DMInitializePackage());
171: PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView));
172: t->setupcalled = PETSC_FALSE;
173: t->redFactor = 2.0;
174: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom));
175: *tr = t;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*@
180: DMPlexTransformSetType - Sets the particular implementation for a transform.
182: Collective
184: Input Parameters:
185: + tr - The transform
186: - method - The name of the transform type
188: Options Database Key:
189: . -dm_plex_transform_type type - Sets the transform type; see `DMPlexTransformType`
191: Level: intermediate
193: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
194: @*/
195: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
196: {
197: PetscErrorCode (*r)(DMPlexTransform);
198: PetscBool match;
200: PetscFunctionBegin;
202: PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match));
203: if (match) PetscFunctionReturn(PETSC_SUCCESS);
205: PetscCall(DMPlexTransformRegisterAll());
206: PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r));
207: PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);
209: PetscTryTypeMethod(tr, destroy);
210: PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops)));
211: PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method));
212: PetscCall((*r)(tr));
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: /*@
217: DMPlexTransformGetType - Gets the type name (as a string) from the transform.
219: Not Collective
221: Input Parameter:
222: . tr - The `DMPlexTransform`
224: Output Parameter:
225: . type - The `DMPlexTransformType` name
227: Level: intermediate
229: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
230: @*/
231: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
232: {
233: PetscFunctionBegin;
235: PetscAssertPointer(type, 2);
236: PetscCall(DMPlexTransformRegisterAll());
237: *type = ((PetscObject)tr)->type_name;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
242: {
243: PetscViewerFormat format;
245: PetscFunctionBegin;
246: PetscCall(PetscViewerGetFormat(v, &format));
247: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
248: const PetscInt *trTypes = NULL;
249: IS trIS;
250: PetscInt cols = 8;
251: PetscInt Nrt = 8, f, g;
253: if (tr->trType) PetscCall(DMLabelView(tr->trType, v));
254: PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n"));
255: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
256: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
257: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]));
258: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
259: PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n"));
260: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
261: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
262: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]));
263: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
265: if (tr->trType) {
266: PetscCall(DMLabelGetNumValues(tr->trType, &Nrt));
267: PetscCall(DMLabelGetValueIS(tr->trType, &trIS));
268: PetscCall(ISGetIndices(trIS, &trTypes));
269: }
270: PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n"));
271: PetscCall(PetscViewerASCIIPrintf(v, " "));
272: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
273: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
274: for (f = 0; f < Nrt; ++f) {
275: PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT " |", trTypes ? trTypes[f] : f));
276: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]));
277: PetscCall(PetscViewerASCIIPrintf(v, " |\n"));
278: }
279: if (tr->trType) {
280: PetscCall(ISRestoreIndices(trIS, &trTypes));
281: PetscCall(ISDestroy(&trIS));
282: }
283: }
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*@
288: DMPlexTransformView - Views a `DMPlexTransform`
290: Collective
292: Input Parameters:
293: + tr - the `DMPlexTransform` object to view
294: - v - the viewer
296: Level: beginner
298: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
299: @*/
300: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
301: {
302: PetscBool isascii;
304: PetscFunctionBegin;
306: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v));
308: PetscCheckSameComm(tr, 1, v, 2);
309: PetscCall(PetscViewerCheckWritable(v));
310: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v));
311: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
312: if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v));
313: PetscTryTypeMethod(tr, view, v);
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: /*@
318: DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database
320: Collective
322: Input Parameter:
323: . tr - the `DMPlexTransform` object to set options for
325: Options Database Keys:
326: + -dm_plex_transform_type type - Set the transform type, e.g. refine_regular
327: . -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point
328: . -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
329: . -dm_plex_transform_active name - Name for active mesh label
330: - -dm_plex_transform_active_values v0,v1,... - Values in the active label
332: Level: intermediate
334: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
335: @*/
336: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
337: {
338: char typeName[1024], active[PETSC_MAX_PATH_LEN];
339: const char *defName = DMPLEXREFINEREGULAR;
340: PetscBool flg, match;
342: PetscFunctionBegin;
344: PetscObjectOptionsBegin((PetscObject)tr);
345: PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg));
346: if (flg) PetscCall(DMPlexTransformSetType(tr, typeName));
347: else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName));
348: PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &match, &flg));
349: if (flg) PetscCall(DMPlexTransformSetMatchStrata(tr, match));
350: 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));
351: PetscCall(PetscOptionsString("-dm_plex_transform_active", "Name for active mesh label", "DMPlexTransformSetActive", active, active, sizeof(active), &flg));
352: if (flg) {
353: DM dm;
354: DMLabel label;
355: PetscInt values[16];
356: PetscInt n = 16;
358: PetscCall(DMPlexTransformGetDM(tr, &dm));
359: PetscCall(DMGetLabel(dm, active, &label));
360: PetscCall(PetscOptionsIntArray("-dm_plex_transform_active_values", "The label values to be active", "DMPlexTransformSetActive", values, &n, &flg));
361: if (flg && n) {
362: DMLabel newlabel;
364: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Active", &newlabel));
365: for (PetscInt i = 0; i < n; ++i) {
366: IS is;
368: PetscCall(DMLabelGetStratumIS(label, values[i], &is));
369: PetscCall(DMLabelInsertIS(newlabel, is, values[i]));
370: PetscCall(ISDestroy(&is));
371: }
372: PetscCall(DMPlexTransformSetActive(tr, newlabel));
373: PetscCall(DMLabelDestroy(&newlabel));
374: } else {
375: PetscCall(DMPlexTransformSetActive(tr, label));
376: }
377: }
378: PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
379: /* process any options handlers added with PetscObjectAddOptionsHandler() */
380: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject));
381: PetscOptionsEnd();
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*@
386: DMPlexTransformDestroy - Destroys a `DMPlexTransform`
388: Collective
390: Input Parameter:
391: . tr - the transform object to destroy
393: Level: beginner
395: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
396: @*/
397: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
398: {
399: PetscInt c;
401: PetscFunctionBegin;
402: if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
404: if (--((PetscObject)*tr)->refct > 0) {
405: *tr = NULL;
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: PetscTryTypeMethod(*tr, destroy);
410: PetscCall(DMDestroy(&(*tr)->dm));
411: PetscCall(DMLabelDestroy(&(*tr)->active));
412: PetscCall(DMLabelDestroy(&(*tr)->trType));
413: PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld));
414: PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew));
415: PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew));
416: PetscCall(PetscFree((*tr)->offset));
417: PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd));
418: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
419: PetscCall(PetscFEDestroy(&(*tr)->coordFE[c]));
420: PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c]));
421: }
422: if ((*tr)->trVerts) {
423: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
424: DMPolytopeType *rct;
425: PetscInt *rsize, *rcone, *rornt, Nct, n, r;
427: if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) {
428: PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
429: for (n = 0; n < Nct; ++n) {
430: if (rct[n] == DM_POLYTOPE_POINT) continue;
431: for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r]));
432: PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]]));
433: }
434: }
435: PetscCall(PetscFree((*tr)->trSubVerts[c]));
436: PetscCall(PetscFree((*tr)->trVerts[c]));
437: }
438: }
439: PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts));
440: PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom));
441: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
442: PetscCall(PetscHeaderDestroy(tr));
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
447: {
448: DMLabel trType = tr->trType;
449: PetscInt c, cN, *off;
451: PetscFunctionBegin;
452: if (trType) {
453: DM dm;
454: IS rtIS;
455: const PetscInt *reftypes;
456: PetscInt Nrt, r;
458: PetscCall(DMPlexTransformGetDM(tr, &dm));
459: PetscCall(DMLabelGetNumValues(trType, &Nrt));
460: PetscCall(DMLabelGetValueIS(trType, &rtIS));
461: PetscCall(ISGetIndices(rtIS, &reftypes));
462: PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off));
463: for (r = 0; r < Nrt; ++r) {
464: const PetscInt rt = reftypes[r];
465: IS rtIS;
466: const PetscInt *points;
467: DMPolytopeType ct;
468: PetscInt np, p;
470: PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
471: PetscCall(ISGetLocalSize(rtIS, &np));
472: PetscCall(ISGetIndices(rtIS, &points));
473: if (!np) continue;
474: p = points[0];
475: PetscCall(ISRestoreIndices(rtIS, &points));
476: PetscCall(ISDestroy(&rtIS));
477: PetscCall(DMPlexGetCellType(dm, p, &ct));
478: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
479: const DMPolytopeType ctNew = (DMPolytopeType)cN;
480: DMPolytopeType *rct;
481: PetscInt *rsize, *cone, *ornt;
482: PetscInt Nct, n, s;
484: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
485: off[r * DM_NUM_POLYTOPES + ctNew] = -1;
486: break;
487: }
488: off[r * DM_NUM_POLYTOPES + ctNew] = 0;
489: for (s = 0; s <= r; ++s) {
490: const PetscInt st = reftypes[s];
491: DMPolytopeType sct;
492: PetscInt q, qrt;
494: PetscCall(DMLabelGetStratumIS(trType, st, &rtIS));
495: PetscCall(ISGetLocalSize(rtIS, &np));
496: PetscCall(ISGetIndices(rtIS, &points));
497: if (!np) continue;
498: q = points[0];
499: PetscCall(ISRestoreIndices(rtIS, &points));
500: PetscCall(ISDestroy(&rtIS));
501: PetscCall(DMPlexGetCellType(dm, q, &sct));
502: PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt));
503: 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);
504: if (st == rt) {
505: for (n = 0; n < Nct; ++n)
506: if (rct[n] == ctNew) break;
507: if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
508: break;
509: }
510: for (n = 0; n < Nct; ++n) {
511: if (rct[n] == ctNew) {
512: PetscInt sn;
514: PetscCall(DMLabelGetStratumSize(trType, st, &sn));
515: off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
516: }
517: }
518: }
519: }
520: }
521: PetscCall(ISRestoreIndices(rtIS, &reftypes));
522: PetscCall(ISDestroy(&rtIS));
523: } else {
524: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off));
525: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
526: const DMPolytopeType ct = (DMPolytopeType)c;
527: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
528: const DMPolytopeType ctNew = (DMPolytopeType)cN;
529: DMPolytopeType *rct;
530: PetscInt *rsize, *cone, *ornt;
531: PetscInt Nct, n, i;
533: 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) {
534: off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
535: continue;
536: }
537: off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
538: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
539: const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i];
540: const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];
542: PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt));
543: if (ict == ct) {
544: for (n = 0; n < Nct; ++n)
545: if (rct[n] == ctNew) break;
546: if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
547: break;
548: }
549: for (n = 0; n < Nct; ++n)
550: if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
551: }
552: }
553: }
554: }
555: *offset = off;
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: /*@
560: DMPlexTransformSetUp - Create the tables that drive the transform
562: Input Parameter:
563: . tr - The `DMPlexTransform` object
565: Level: intermediate
567: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
568: @*/
569: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
570: {
571: DMPolytopeType ctCell;
572: DM dm;
573: PetscInt pStart, pEnd, p, c, celldim = 0;
575: PetscFunctionBegin;
577: if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
578: PetscCall(DMPlexTransformGetDM(tr, &dm));
579: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0));
580: PetscTryTypeMethod(tr, setup);
581: PetscCall(DMSetSnapToGeomModel(dm, NULL));
582: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
584: if (pEnd > pStart) {
585: // Ignore cells hanging off of embedded surfaces
586: PetscInt c = pStart;
588: ctCell = DM_POLYTOPE_FV_GHOST;
589: while (DMPolytopeTypeGetDim(ctCell) < 0) PetscCall(DMPlexGetCellType(dm, c++, &ctCell));
590: } else {
591: PetscInt dim;
593: PetscCall(DMGetDimension(dm, &dim));
594: switch (dim) {
595: case 0:
596: ctCell = DM_POLYTOPE_POINT;
597: break;
598: case 1:
599: ctCell = DM_POLYTOPE_SEGMENT;
600: break;
601: case 2:
602: ctCell = DM_POLYTOPE_TRIANGLE;
603: break;
604: case 3:
605: ctCell = DM_POLYTOPE_TETRAHEDRON;
606: break;
607: default:
608: ctCell = DM_POLYTOPE_UNKNOWN;
609: }
610: }
611: PetscCall(DMPlexCreateCellTypeOrder_Internal(dm, DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
612: for (p = pStart; p < pEnd; ++p) {
613: DMPolytopeType ct;
614: DMPolytopeType *rct;
615: PetscInt *rsize, *cone, *ornt;
616: PetscInt Nct, n;
618: PetscCall(DMPlexGetCellType(dm, p, &ct));
619: 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);
620: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
621: for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
622: }
623: PetscCall(DMPlexCreateCellTypeOrder_Internal(NULL, celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
624: /* Construct sizes and offsets for each cell type */
625: if (!tr->ctStart) {
626: PetscInt *ctS, *ctSN, *ctC, *ctCN;
628: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
629: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
630: for (p = pStart; p < pEnd; ++p) {
631: DMPolytopeType ct;
632: DMPolytopeType *rct;
633: PetscInt *rsize, *cone, *ornt;
634: PetscInt Nct, n;
636: PetscCall(DMPlexGetCellType(dm, p, &ct));
637: 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);
638: ++ctC[ct];
639: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
640: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
641: }
642: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
643: const PetscInt cto = tr->ctOrderOld[c];
644: const PetscInt cton = tr->ctOrderOld[c + 1];
645: const PetscInt ctn = tr->ctOrderNew[c];
646: const PetscInt ctnn = tr->ctOrderNew[c + 1];
648: ctS[cton] = ctS[cto] + ctC[cto];
649: ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
650: }
651: PetscCall(PetscFree2(ctC, ctCN));
652: tr->ctStart = ctS;
653: tr->ctStartNew = ctSN;
654: }
655: PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
656: // Compute depth information
657: tr->depth = -1;
658: for (c = 0; c < DM_NUM_POLYTOPES; ++c)
659: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
660: PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
661: for (PetscInt d = 0; d <= tr->depth; ++d) {
662: tr->depthStart[d] = PETSC_INT_MAX;
663: tr->depthEnd[d] = -1;
664: }
665: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
666: const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);
668: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
669: tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
670: tr->depthEnd[dep] = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
671: }
672: tr->setupcalled = PETSC_TRUE;
673: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0));
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: /*@
678: DMPlexTransformGetDM - Get the base `DM` for the transform
680: Input Parameter:
681: . tr - The `DMPlexTransform` object
683: Output Parameter:
684: . dm - The original `DM` which will be transformed
686: Level: intermediate
688: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
689: @*/
690: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
691: {
692: PetscFunctionBegin;
694: PetscAssertPointer(dm, 2);
695: *dm = tr->dm;
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: /*@
700: DMPlexTransformSetDM - Set the base `DM` for the transform
702: Input Parameters:
703: + tr - The `DMPlexTransform` object
704: - dm - The original `DM` which will be transformed
706: Level: intermediate
708: Note:
709: The user does not typically call this, as it is called by `DMPlexTransformApply()`.
711: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
712: @*/
713: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
714: {
715: PetscFunctionBegin;
718: PetscCall(PetscObjectReference((PetscObject)dm));
719: PetscCall(DMDestroy(&tr->dm));
720: tr->dm = dm;
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: /*@
725: DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform
727: Input Parameter:
728: . tr - The `DMPlexTransform` object
730: Output Parameter:
731: . active - The `DMLabel` indicating which points will be transformed
733: Level: intermediate
735: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
736: @*/
737: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
738: {
739: PetscFunctionBegin;
741: PetscAssertPointer(active, 2);
742: *active = tr->active;
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: /*@
747: DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform
749: Input Parameters:
750: + tr - The `DMPlexTransform` object
751: - active - The `DMLabel` indicating which points will be transformed
753: Level: intermediate
755: Note:
756: This only applies to transforms listed in [](plex_transform_table) that operate on a subset of the mesh.
758: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
759: @*/
760: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
761: {
762: PetscFunctionBegin;
765: PetscCall(PetscObjectReference((PetscObject)active));
766: PetscCall(DMLabelDestroy(&tr->active));
767: tr->active = active;
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: /*@
772: DMPlexTransformGetTransformTypes - Get the `DMLabel` marking the transform type of each point for the transform
774: Input Parameter:
775: . tr - The `DMPlexTransform` object
777: Output Parameter:
778: . trType - The `DMLabel` indicating the transform type for each point
780: Level: intermediate
782: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetTransformType()`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
783: @*/
784: PetscErrorCode DMPlexTransformGetTransformTypes(DMPlexTransform tr, DMLabel *trType)
785: {
786: PetscFunctionBegin;
788: PetscAssertPointer(trType, 2);
789: *trType = tr->trType;
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: /*@
794: DMPlexTransformSetTransformTypes - Set the `DMLabel` marking the transform type of each point for the transform
796: Input Parameters:
797: + tr - The `DMPlexTransform` object
798: - trType - The original `DM` which will be transformed
800: Level: intermediate
802: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetTransformTypes()`, `DMPlexTransformGetActive())`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
803: @*/
804: PetscErrorCode DMPlexTransformSetTransformTypes(DMPlexTransform tr, DMLabel trType)
805: {
806: PetscFunctionBegin;
809: PetscCall(PetscObjectReference((PetscObject)trType));
810: PetscCall(DMLabelDestroy(&tr->trType));
811: tr->trType = trType;
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
816: {
817: PetscFunctionBegin;
818: if (!tr->coordFE[ct]) {
819: PetscInt dim, cdim;
821: dim = DMPolytopeTypeGetDim(ct);
822: PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
823: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
824: {
825: PetscDualSpace dsp;
826: PetscQuadrature quad;
827: DM K;
828: PetscFEGeom *cg;
829: PetscScalar *Xq;
830: PetscReal *xq, *wq;
831: PetscInt Nq, q;
833: PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
834: PetscCall(PetscMalloc1(Nq * cdim, &xq));
835: for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
836: PetscCall(PetscMalloc1(Nq, &wq));
837: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
838: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
839: PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
840: PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));
842: PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
843: PetscCall(PetscDualSpaceGetDM(dsp, &K));
844: PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FEGEOM_BASIC, &tr->refGeom[ct]));
845: cg = tr->refGeom[ct];
846: PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
847: PetscCall(PetscQuadratureDestroy(&quad));
848: }
849: }
850: *fe = tr->coordFE[ct];
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
855: {
856: PetscInt dim, cdim;
858: PetscFunctionBegin;
859: PetscCall(DMGetDimension(dm, &dim));
860: PetscCall(DMSetDimension(tdm, dim));
861: PetscCall(DMGetCoordinateDim(dm, &cdim));
862: PetscCall(DMSetCoordinateDim(tdm, cdim));
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: /*@
867: DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`
869: Input Parameters:
870: + tr - The `DMPlexTransform` object
871: - dm - The original `DM`
873: Output Parameter:
874: . trdm - The transformed `DM`
876: Level: advanced
878: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
879: @*/
880: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM trdm)
881: {
882: PetscFunctionBegin;
883: PetscUseTypeMethod(tr, setdimensions, dm, trdm);
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: /*@
888: DMPlexTransformGetChart - Get the chart `[pStart, pEnd)` for the points produced by the transform
890: Not Collective
892: Input Parameter:
893: . tr - The `DMPlexTransform`
895: Output Parameters:
896: + pStart - The first point in the transformed mesh, or `NULL` if not needed
897: - pEnd - One past the last point in the transformed mesh, or `NULL` if not needed
899: Level: developer
901: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformGetCellType()`, `DMPlexTransformGetCellTypeStratum()`
902: @*/
903: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
904: {
905: PetscFunctionBegin;
906: if (pStart) *pStart = 0;
907: if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
908: PetscFunctionReturn(PETSC_SUCCESS);
909: }
911: /*@
912: DMPlexTransformGetCellType - Return the cell type for a point in the transformed mesh
914: Not Collective
916: Input Parameters:
917: + tr - The `DMPlexTransform`
918: - cell - The point number in the transformed mesh
920: Output Parameter:
921: . celltype - The `DMPolytopeType` of the point
923: Level: developer
925: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetChart()`, `DMPlexTransformGetCellTypeStratum()`
926: @*/
927: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
928: {
929: PetscInt ctNew;
931: PetscFunctionBegin;
933: PetscAssertPointer(celltype, 3);
934: /* TODO Can do bisection since everything is sorted */
935: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
936: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
938: if (cell >= ctSN && cell < ctEN) break;
939: }
940: PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
941: *celltype = (DMPolytopeType)ctNew;
942: PetscFunctionReturn(PETSC_SUCCESS);
943: }
945: /*@
946: DMPlexTransformGetCellTypeStratum - Return the point range for a given cell type in the transformed mesh
948: Not Collective
950: Input Parameters:
951: + tr - The `DMPlexTransform`
952: - celltype - The `DMPolytopeType` of the requested stratum
954: Output Parameters:
955: + start - The first point of the stratum, or `NULL` if not needed
956: - end - One past the last point of the stratum, or `NULL` if not needed
958: Level: developer
960: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellType()`, `DMPlexTransformGetChart()`, `DMPlexGetDepthStratum()`
961: @*/
962: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
963: {
964: PetscFunctionBegin;
966: if (start) *start = tr->ctStartNew[celltype];
967: if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@
972: DMPlexTransformGetDepth - Return the topological depth of the transformed mesh
974: Not Collective
976: Input Parameter:
977: . tr - The `DMPlexTransform`
979: Output Parameter:
980: . depth - The depth of the transformed mesh
982: Level: developer
984: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDepthStratum()`, `DMPlexGetDepth()`
985: @*/
986: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
987: {
988: PetscFunctionBegin;
990: *depth = tr->depth;
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: /*@
995: DMPlexTransformGetDepthStratum - Return the point range for a given depth in the transformed mesh
997: Not Collective
999: Input Parameters:
1000: + tr - The `DMPlexTransform`
1001: - depth - The requested depth in the transformed mesh
1003: Output Parameters:
1004: + start - The first point at the given depth, or `NULL` if not needed
1005: - end - One past the last point at the given depth, or `NULL` if not needed
1007: Level: developer
1009: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDepth()`, `DMPlexGetDepthStratum()`
1010: @*/
1011: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
1012: {
1013: PetscFunctionBegin;
1015: if (start) *start = tr->depthStart[depth];
1016: if (end) *end = tr->depthEnd[depth];
1017: PetscFunctionReturn(PETSC_SUCCESS);
1018: }
1020: /*@
1021: DMPlexTransformGetMatchStrata - Get the flag which determines what points get added to the transformed labels
1023: Not Collective
1025: Input Parameter:
1026: . tr - The `DMPlexTransform`
1028: Output Parameter:
1029: . match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels
1031: Level: intermediate
1033: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetMatchStrata()`, `DMPlexGetPointDepth()`
1034: @*/
1035: PetscErrorCode DMPlexTransformGetMatchStrata(DMPlexTransform tr, PetscBool *match)
1036: {
1037: PetscFunctionBegin;
1039: PetscAssertPointer(match, 2);
1040: *match = tr->labelMatchStrata;
1041: PetscFunctionReturn(PETSC_SUCCESS);
1042: }
1044: /*@
1045: DMPlexTransformSetMatchStrata - Set the flag which determines what points get added to the transformed labels
1047: Not Collective
1049: Input Parameters:
1050: + tr - The `DMPlexTransform`
1051: - match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels
1053: Level: intermediate
1055: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetMatchStrata()`, `DMPlexGetPointDepth()`
1056: @*/
1057: PetscErrorCode DMPlexTransformSetMatchStrata(DMPlexTransform tr, PetscBool match)
1058: {
1059: PetscFunctionBegin;
1061: tr->labelMatchStrata = match;
1062: PetscFunctionReturn(PETSC_SUCCESS);
1063: }
1065: /*@
1066: DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.
1068: Not Collective
1070: Input Parameters:
1071: + tr - The `DMPlexTransform`
1072: . ct - The type of the original point which produces the new point
1073: . ctNew - The type of the new point
1074: . p - The original point which produces the new point
1075: - r - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`
1077: Output Parameter:
1078: . pNew - The new point number
1080: Level: developer
1082: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
1083: @*/
1084: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
1085: {
1086: DMPolytopeType *rct;
1087: PetscInt *rsize, *cone, *ornt;
1088: PetscInt rt, Nct, n, off, rp;
1089: DMLabel trType = tr->trType;
1090: PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
1091: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
1092: PetscInt newp = ctSN, cind;
1094: PetscFunctionBeginHot;
1095: 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);
1096: PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
1097: if (trType) {
1098: PetscCall(DMLabelGetValueIndex(trType, rt, &cind));
1099: PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
1100: 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);
1101: } else {
1102: cind = ct;
1103: rp = p - ctS;
1104: }
1105: off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
1106: 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);
1107: newp += off;
1108: for (n = 0; n < Nct; ++n) {
1109: if (rct[n] == ctNew) {
1110: if (rsize[n] && r >= rsize[n])
1111: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %" PetscInt_FMT " for point %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ") for subcell type %s in cell type %s", r, p, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
1112: newp += rp * rsize[n] + r;
1113: break;
1114: }
1115: }
1117: 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);
1118: *pNew = newp;
1119: PetscFunctionReturn(PETSC_SUCCESS);
1120: }
1122: /*@
1123: DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.
1125: Not Collective
1127: Input Parameters:
1128: + tr - The `DMPlexTransform`
1129: - pNew - The new point number
1131: Output Parameters:
1132: + ct - The type of the original point which produces the new point
1133: . ctNew - The type of the new point
1134: . p - The original point which produces the new point
1135: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
1137: Level: developer
1139: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
1140: @*/
1141: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
1142: {
1143: DMLabel trType = tr->trType;
1144: DMPolytopeType *rct, ctN;
1145: PetscInt *rsize, *cone, *ornt;
1146: PetscInt rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
1147: PetscInt offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;
1149: PetscFunctionBegin;
1150: PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
1151: if (trType) {
1152: DM dm;
1153: IS rtIS;
1154: const PetscInt *reftypes;
1155: PetscInt Nrt, r, rtStart;
1157: PetscCall(DMPlexTransformGetDM(tr, &dm));
1158: PetscCall(DMLabelGetNumValues(trType, &Nrt));
1159: PetscCall(DMLabelGetValueIS(trType, &rtIS));
1160: PetscCall(ISGetIndices(rtIS, &reftypes));
1161: for (r = 0; r < Nrt; ++r) {
1162: const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];
1164: if (tr->ctStartNew[ctN] + off > pNew) continue;
1165: /* Check that any of this refinement type exist */
1166: /* TODO Actually keep track of the number produced here instead */
1167: if (off > offset) {
1168: rt = reftypes[r];
1169: offset = off;
1170: }
1171: }
1172: PetscCall(ISRestoreIndices(rtIS, &reftypes));
1173: PetscCall(ISDestroy(&rtIS));
1174: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
1175: /* TODO Map refinement types to cell types */
1176: PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL));
1177: PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
1178: for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
1179: PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
1181: if ((rtStart >= ctS) && (rtStart < ctE)) break;
1182: }
1183: PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt);
1184: } else {
1185: for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
1186: const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];
1188: if (tr->ctStartNew[ctN] + off > pNew) continue;
1189: if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
1190: /* TODO Actually keep track of the number produced here instead */
1191: if (off > offset) {
1192: ctO = ctTmp;
1193: offset = off;
1194: }
1195: }
1196: rt = -1;
1197: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
1198: }
1199: ctS = tr->ctStart[ctO];
1200: ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
1201: if (trType) {
1202: for (rtS = ctS; rtS < ctE; ++rtS) {
1203: PetscInt val;
1204: PetscCall(DMLabelGetValue(trType, rtS, &val));
1205: if (val == rt) break;
1206: }
1207: 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);
1208: } else rtS = ctS;
1209: PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
1210: 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);
1211: for (n = 0; n < Nct; ++n) {
1212: if (rct[n] == ctN) {
1213: PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c;
1215: if (trType) {
1216: for (c = ctS; c < ctE; ++c) {
1217: PetscCall(DMLabelGetValue(trType, c, &val));
1218: if (val == rt) {
1219: if (tmp < rsize[n]) break;
1220: tmp -= rsize[n];
1221: }
1222: }
1223: PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew);
1224: rp = c - ctS;
1225: rO = tmp;
1226: } else {
1227: // This assumes that all points of type ctO transform the same way
1228: rp = tmp / rsize[n];
1229: rO = tmp % rsize[n];
1230: }
1231: break;
1232: }
1233: }
1234: PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
1235: pO = rp + ctS;
1236: 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);
1237: if (ct) *ct = (DMPolytopeType)ctO;
1238: if (ctNew) *ctNew = ctN;
1239: if (p) *p = pO;
1240: if (r) *r = rO;
1241: PetscFunctionReturn(PETSC_SUCCESS);
1242: }
1244: /*@
1245: DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.
1247: Input Parameters:
1248: + tr - The `DMPlexTransform` object
1249: . source - The source cell type
1250: - p - The source point, which can also determine the refine type
1252: Output Parameters:
1253: + rt - The refine type for this point
1254: . Nt - The number of types produced by this point
1255: . target - An array of length `Nt` giving the types produced
1256: . size - An array of length `Nt` giving the number of cells of each type produced
1257: . cone - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
1258: - ornt - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point
1260: Level: advanced
1262: Notes:
1263: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
1264: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
1265: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
1266: .vb
1267: the number of cones to be taken, or 0 for the current cell
1268: the cell cone point number at each level from which it is subdivided
1269: the replica number r of the subdivision.
1270: .ve
1271: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
1272: .vb
1273: Nt = 2
1274: target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
1275: size = {1, 2}
1276: cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
1277: ornt = { 0, 0, 0, 0}
1278: .ve
1280: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
1281: @*/
1282: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1283: {
1284: PetscFunctionBegin;
1285: PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
1286: PetscFunctionReturn(PETSC_SUCCESS);
1287: }
1289: /*@
1290: DMPlexTransformGetSubcellOrientationIdentity - Default `getsubcellorientation` implementation for transforms that reproduce the input mesh
1292: Not Collective
1294: Input Parameters:
1295: + tr - The `DMPlexTransform`
1296: . sct - The source point cell type
1297: . sp - The source point
1298: . so - The orientation of the source point in its enclosing parent
1299: . tct - The target point cell type
1300: . r - The replica number requested for the produced cell type
1301: - o - The orientation of the replica
1303: Output Parameters:
1304: + rnew - The replica number, given the orientation of the parent (returns `r`)
1305: - onew - The replica orientation composed with the source orientation
1307: Level: developer
1309: Note:
1310: This is the identity variant used by transforms such as the "identity" refiner where each source
1311: point produces itself, so the replica number is unchanged and the returned orientation is simply
1312: `o` composed with `so` via `DMPolytopeTypeComposeOrientation()`.
1314: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetSubcellOrientation()`, `DMPlexTransformCellTransformIdentity()`, `DMPolytopeTypeComposeOrientation()`
1315: @*/
1316: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1317: {
1318: PetscFunctionBegin;
1319: *rnew = r;
1320: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
1321: PetscFunctionReturn(PETSC_SUCCESS);
1322: }
1324: /*@
1325: DMPlexTransformCellTransformIdentity - Default `celltransform` implementation for transforms that reproduce the input mesh
1327: Not Collective
1329: Input Parameters:
1330: + tr - The `DMPlexTransform`
1331: . source - The cell type of the source point
1332: - p - The source point
1334: Output Parameters:
1335: + rt - Refinement type of the source point (set to 0), or `NULL`
1336: . Nt - Number of target cell types produced (always 1)
1337: . target - Array of produced cell types (a single-element array containing `source`)
1338: . size - Array of replica counts for each produced type (a single-element array containing 1)
1339: . cone - Cone description used by `DMPlexTransformGetCone()`; encodes that the replica takes the entire parent cone
1340: - ornt - Orientation array associated with `cone`; all zero for identity
1342: Level: developer
1344: Note:
1345: This routine returns statically allocated arrays describing an identity refinement for each supported
1346: `DMPolytopeType`; every source point produces a single replica of the same type with unchanged cone
1347: and orientation.
1349: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformGetSubcellOrientationIdentity()`
1350: @*/
1351: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1352: {
1353: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
1354: static PetscInt vertexS[] = {1};
1355: static PetscInt vertexC[] = {0};
1356: static PetscInt vertexO[] = {0};
1357: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
1358: static PetscInt edgeS[] = {1};
1359: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1360: static PetscInt edgeO[] = {0, 0};
1361: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
1362: static PetscInt tedgeS[] = {1};
1363: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1364: static PetscInt tedgeO[] = {0, 0};
1365: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
1366: static PetscInt triS[] = {1};
1367: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1368: static PetscInt triO[] = {0, 0, 0};
1369: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
1370: static PetscInt quadS[] = {1};
1371: 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};
1372: static PetscInt quadO[] = {0, 0, 0, 0};
1373: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
1374: static PetscInt tquadS[] = {1};
1375: 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};
1376: static PetscInt tquadO[] = {0, 0, 0, 0};
1377: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
1378: static PetscInt tetS[] = {1};
1379: 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};
1380: static PetscInt tetO[] = {0, 0, 0, 0};
1381: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
1382: static PetscInt hexS[] = {1};
1383: 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};
1384: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
1385: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
1386: static PetscInt tripS[] = {1};
1387: 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};
1388: static PetscInt tripO[] = {0, 0, 0, 0, 0};
1389: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
1390: static PetscInt ttripS[] = {1};
1391: 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};
1392: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
1393: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
1394: static PetscInt tquadpS[] = {1};
1395: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
1396: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1397: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
1398: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
1399: static PetscInt pyrS[] = {1};
1400: 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};
1401: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
1403: PetscFunctionBegin;
1404: if (rt) *rt = 0;
1405: switch (source) {
1406: case DM_POLYTOPE_POINT:
1407: *Nt = 1;
1408: *target = vertexT;
1409: *size = vertexS;
1410: *cone = vertexC;
1411: *ornt = vertexO;
1412: break;
1413: case DM_POLYTOPE_SEGMENT:
1414: *Nt = 1;
1415: *target = edgeT;
1416: *size = edgeS;
1417: *cone = edgeC;
1418: *ornt = edgeO;
1419: break;
1420: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1421: *Nt = 1;
1422: *target = tedgeT;
1423: *size = tedgeS;
1424: *cone = tedgeC;
1425: *ornt = tedgeO;
1426: break;
1427: case DM_POLYTOPE_TRIANGLE:
1428: *Nt = 1;
1429: *target = triT;
1430: *size = triS;
1431: *cone = triC;
1432: *ornt = triO;
1433: break;
1434: case DM_POLYTOPE_QUADRILATERAL:
1435: *Nt = 1;
1436: *target = quadT;
1437: *size = quadS;
1438: *cone = quadC;
1439: *ornt = quadO;
1440: break;
1441: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1442: *Nt = 1;
1443: *target = tquadT;
1444: *size = tquadS;
1445: *cone = tquadC;
1446: *ornt = tquadO;
1447: break;
1448: case DM_POLYTOPE_TETRAHEDRON:
1449: *Nt = 1;
1450: *target = tetT;
1451: *size = tetS;
1452: *cone = tetC;
1453: *ornt = tetO;
1454: break;
1455: case DM_POLYTOPE_HEXAHEDRON:
1456: *Nt = 1;
1457: *target = hexT;
1458: *size = hexS;
1459: *cone = hexC;
1460: *ornt = hexO;
1461: break;
1462: case DM_POLYTOPE_TRI_PRISM:
1463: *Nt = 1;
1464: *target = tripT;
1465: *size = tripS;
1466: *cone = tripC;
1467: *ornt = tripO;
1468: break;
1469: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1470: *Nt = 1;
1471: *target = ttripT;
1472: *size = ttripS;
1473: *cone = ttripC;
1474: *ornt = ttripO;
1475: break;
1476: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1477: *Nt = 1;
1478: *target = tquadpT;
1479: *size = tquadpS;
1480: *cone = tquadpC;
1481: *ornt = tquadpO;
1482: break;
1483: case DM_POLYTOPE_PYRAMID:
1484: *Nt = 1;
1485: *target = pyrT;
1486: *size = pyrS;
1487: *cone = pyrC;
1488: *ornt = pyrO;
1489: break;
1490: default:
1491: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1492: }
1493: PetscFunctionReturn(PETSC_SUCCESS);
1494: }
1496: /*@
1497: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
1499: Not Collective
1501: Input Parameters:
1502: + tr - The `DMPlexTransform`
1503: . sct - The source point cell type, from whom the new cell is being produced
1504: . sp - The source point
1505: . so - The orientation of the source point in its enclosing parent
1506: . tct - The target point cell type
1507: . r - The replica number requested for the produced cell type
1508: - o - The orientation of the replica
1510: Output Parameters:
1511: + rnew - The replica number, given the orientation of the parent
1512: - onew - The replica orientation, given the orientation of the parent
1514: Level: advanced
1516: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1517: @*/
1518: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1519: {
1520: PetscFunctionBeginHot;
1521: PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1522: PetscFunctionReturn(PETSC_SUCCESS);
1523: }
1525: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1526: {
1527: DM dm;
1528: PetscInt pStart, pEnd, pNew;
1530: PetscFunctionBegin;
1531: PetscCall(DMPlexTransformGetDM(tr, &dm));
1532: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0));
1533: /* Must create the celltype label here so that we do not automatically try to compute the types */
1534: PetscCall(DMCreateLabel(rdm, "celltype"));
1535: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1536: for (PetscInt p = pStart; p < pEnd; ++p) {
1537: DMPolytopeType ct;
1538: DMPolytopeType *rct;
1539: PetscInt *rsize, *rcone, *rornt;
1540: PetscInt Nct, n, r;
1542: PetscCall(DMPlexGetCellType(dm, p, &ct));
1543: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1544: for (n = 0; n < Nct; ++n) {
1545: for (r = 0; r < rsize[n]; ++r) {
1546: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1547: PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
1548: PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
1549: }
1550: }
1551: }
1552: /* Let the DM know we have set all the cell types */
1553: {
1554: DMLabel ctLabel;
1555: DM_Plex *plex = (DM_Plex *)rdm->data;
1557: PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
1558: PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
1559: }
1560: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0));
1561: PetscFunctionReturn(PETSC_SUCCESS);
1562: }
1564: /*@
1565: DMPlexTransformGetConeSize - Return the cone size of a point in the transformed mesh
1567: Not Collective
1569: Input Parameters:
1570: + tr - The `DMPlexTransform`
1571: - q - The point number in the transformed mesh
1573: Output Parameter:
1574: . coneSize - The number of points in the cone of `q`
1576: Level: developer
1578: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetCone()`, `DMPlexTransformGetCellType()`, `DMPlexGetConeSize()`
1579: @*/
1580: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1581: {
1582: DMPolytopeType ctNew;
1584: PetscFunctionBegin;
1586: PetscAssertPointer(coneSize, 3);
1587: PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
1588: *coneSize = DMPolytopeTypeGetConeSize(ctNew);
1589: PetscFunctionReturn(PETSC_SUCCESS);
1590: }
1592: /* The orientation o is for the interior of the cell p */
1593: 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[])
1594: {
1595: DM dm;
1596: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1597: const PetscInt *cone;
1598: DMPolytopeType *newft = NULL;
1599: PetscInt c, coff = *coneoff, ooff = *orntoff;
1600: PetscInt dim, cr = 0, co = 0, nr, no;
1602: PetscFunctionBegin;
1603: PetscCall(DMPlexTransformGetDM(tr, &dm));
1604: PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
1605: // Check if we have to permute this cell
1606: PetscCall(DMGetDimension(dm, &dim));
1607: if (DMPolytopeTypeGetDim(ctNew) == dim && DMPolytopeTypeGetDim(ct) == dim - 1) {
1608: PetscCall(DMPlexTransformGetSubcellOrientation(tr, ct, p, o, ctNew, cr, co, &nr, &no));
1609: if (cr != nr || co != no) PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newft));
1610: }
1611: for (c = 0; c < csizeNew; ++c) {
1612: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1613: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1614: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1615: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1616: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1617: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1618: const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1619: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1620: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1621: PetscInt lc;
1623: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1624: for (lc = 0; lc < fn; ++lc) {
1625: const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po);
1626: const PetscInt acp = rcone[coff++];
1627: const PetscInt pcp = parr[acp * 2];
1628: const PetscInt pco = parr[acp * 2 + 1];
1629: const PetscInt *ppornt;
1631: ppp = pp;
1632: pp = pcone[pcp];
1633: PetscCall(DMPlexGetCellType(dm, pp, &pct));
1634: // Restore the parent cone from the last iterate
1635: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
1636: PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
1637: PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
1638: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1639: PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
1640: }
1641: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
1642: pr = rcone[coff++];
1643: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1644: PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
1645: PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
1646: orntNew[c] = fo;
1647: if (newft) newft[c] = ft;
1648: }
1649: PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
1650: if (newft) {
1651: const PetscInt *arr;
1652: PetscInt *newcone, *newornt;
1654: arr = DMPolytopeTypeGetArrangement(ctNew, no);
1655: PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newcone));
1656: PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newornt));
1657: for (PetscInt c = 0; c < csizeNew; ++c) {
1658: DMPolytopeType ft = newft[c];
1659: PetscInt nO;
1661: nO = DMPolytopeTypeGetNumArrangements(ft) / 2;
1662: newcone[c] = coneNew[arr[c * 2 + 0]];
1663: newornt[c] = DMPolytopeTypeComposeOrientation(ft, arr[c * 2 + 1], orntNew[arr[c * 2 + 0]]);
1664: PetscCheck(!newornt[c] || !(newornt[c] >= nO || newornt[c] < -nO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid orientation %" PetscInt_FMT " not in [%" PetscInt_FMT ",%" PetscInt_FMT ") for %s %" PetscInt_FMT, newornt[c], -nO, nO, DMPolytopeTypes[ft], coneNew[c]);
1665: }
1666: for (PetscInt c = 0; c < csizeNew; ++c) {
1667: coneNew[c] = newcone[c];
1668: orntNew[c] = newornt[c];
1669: }
1670: PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newcone));
1671: PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newornt));
1672: PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newft));
1673: }
1674: *coneoff = coff;
1675: *orntoff = ooff;
1676: PetscFunctionReturn(PETSC_SUCCESS);
1677: }
1679: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1680: {
1681: DM dm;
1682: DMPolytopeType ct;
1683: PetscInt *coneNew, *orntNew;
1684: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1686: PetscFunctionBegin;
1687: PetscCall(DMPlexTransformGetDM(tr, &dm));
1688: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0));
1689: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1690: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1691: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1692: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1693: for (p = pStart; p < pEnd; ++p) {
1694: PetscInt coff, ooff;
1695: DMPolytopeType *rct;
1696: PetscInt *rsize, *rcone, *rornt;
1697: PetscInt Nct, n, r;
1699: PetscCall(DMPlexGetCellType(dm, p, &ct));
1700: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1701: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1702: const DMPolytopeType ctNew = rct[n];
1704: for (r = 0; r < rsize[n]; ++r) {
1705: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1706: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
1707: PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
1708: PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
1709: }
1710: }
1711: }
1712: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1713: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1714: PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
1715: PetscCall(DMPlexSymmetrize(rdm));
1716: PetscCall(DMPlexStratify(rdm));
1717: PetscTryTypeMethod(tr, ordersupports, dm, rdm);
1718: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0));
1719: PetscFunctionReturn(PETSC_SUCCESS);
1720: }
1722: /*@
1723: DMPlexTransformGetConeOriented - Return the cone of a point in the transformed mesh, computed using a specified parent orientation
1725: Not Collective
1727: Input Parameters:
1728: + tr - The `DMPlexTransform`
1729: . q - The point number in the transformed mesh
1730: - po - The orientation of the parent cell in the original mesh to use when producing the cone
1732: Output Parameters:
1733: + cone - The cone points, obtained from an internal work array
1734: - ornt - The orientations of the cone points, obtained from an internal work array
1736: Level: developer
1738: Note:
1739: Both `cone` and `ornt` are returned in work arrays that must be released with `DMPlexTransformRestoreCone()`.
1741: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetCone()`, `DMPlexTransformRestoreCone()`, `DMPlexTransformGetConeSize()`
1742: @*/
1743: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1744: {
1745: DM dm;
1746: DMPolytopeType ct, qct;
1747: DMPolytopeType *rct;
1748: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1749: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1751: PetscFunctionBegin;
1753: PetscAssertPointer(cone, 4);
1754: PetscAssertPointer(ornt, 5);
1755: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1756: PetscCall(DMPlexTransformGetDM(tr, &dm));
1757: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1758: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1759: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1760: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1761: for (n = 0; n < Nct; ++n) {
1762: const DMPolytopeType ctNew = rct[n];
1763: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1764: PetscInt Nr = rsize[n], fn, c;
1766: if (ctNew == qct) Nr = r;
1767: for (nr = 0; nr < Nr; ++nr) {
1768: for (c = 0; c < csizeNew; ++c) {
1769: ++coff; /* Cell type of new cone point */
1770: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1771: coff += fn;
1772: ++coff; /* Replica number of new cone point */
1773: ++ooff; /* Orientation of new cone point */
1774: }
1775: }
1776: if (ctNew == qct) break;
1777: }
1778: PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1779: *cone = qcone;
1780: *ornt = qornt;
1781: PetscFunctionReturn(PETSC_SUCCESS);
1782: }
1784: /*@
1785: DMPlexTransformGetCone - Return the cone of a point in the transformed mesh
1787: Not Collective
1789: Input Parameters:
1790: + tr - The `DMPlexTransform`
1791: - q - The point number in the transformed mesh
1793: Output Parameters:
1794: + cone - The cone points, obtained from an internal work array, or `NULL` if not requested
1795: - ornt - The orientations of the cone points, obtained from an internal work array, or `NULL` if not requested
1797: Level: developer
1799: Note:
1800: Any non-`NULL` output must be released with `DMPlexTransformRestoreCone()`.
1802: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRestoreCone()`, `DMPlexTransformGetConeOriented()`, `DMPlexTransformGetConeSize()`, `DMPlexGetCone()`
1803: @*/
1804: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1805: {
1806: DM dm;
1807: DMPolytopeType ct, qct;
1808: DMPolytopeType *rct;
1809: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1810: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1812: PetscFunctionBegin;
1814: if (cone) PetscAssertPointer(cone, 3);
1815: if (ornt) PetscAssertPointer(ornt, 4);
1816: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1817: PetscCall(DMPlexTransformGetDM(tr, &dm));
1818: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1819: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1820: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1821: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1822: for (n = 0; n < Nct; ++n) {
1823: const DMPolytopeType ctNew = rct[n];
1824: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1825: PetscInt Nr = rsize[n], fn, c;
1827: if (ctNew == qct) Nr = r;
1828: for (nr = 0; nr < Nr; ++nr) {
1829: for (c = 0; c < csizeNew; ++c) {
1830: ++coff; /* Cell type of new cone point */
1831: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1832: coff += fn;
1833: ++coff; /* Replica number of new cone point */
1834: ++ooff; /* Orientation of new cone point */
1835: }
1836: }
1837: if (ctNew == qct) break;
1838: }
1839: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1840: if (cone) *cone = qcone;
1841: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1842: if (ornt) *ornt = qornt;
1843: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1844: PetscFunctionReturn(PETSC_SUCCESS);
1845: }
1847: /*@
1848: DMPlexTransformRestoreCone - Return the work arrays produced by `DMPlexTransformGetCone()` or `DMPlexTransformGetConeOriented()`
1850: Not Collective
1852: Input Parameters:
1853: + tr - The `DMPlexTransform`
1854: . q - The point number in the transformed mesh
1855: . cone - The cone points to release, or `NULL`
1856: - ornt - The orientations to release, or `NULL`
1858: Level: developer
1860: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetCone()`, `DMPlexTransformGetConeOriented()`
1861: @*/
1862: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1863: {
1864: DM dm;
1866: PetscFunctionBegin;
1868: PetscCall(DMPlexTransformGetDM(tr, &dm));
1869: if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
1870: if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
1871: PetscFunctionReturn(PETSC_SUCCESS);
1872: }
1874: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1875: {
1876: PetscInt ict;
1878: PetscFunctionBegin;
1879: PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
1880: for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1881: const DMPolytopeType ct = (DMPolytopeType)ict;
1882: DMPlexTransform reftr;
1883: DM refdm, trdm;
1884: Vec coordinates;
1885: const PetscScalar *coords;
1886: DMPolytopeType *rct;
1887: PetscInt *rsize, *rcone, *rornt;
1888: PetscInt Nct, n, r, pNew = 0;
1889: PetscInt trdim, vStart, vEnd, Nc;
1890: const PetscInt debug = 0;
1891: const char *typeName;
1893: /* Since points are 0-dimensional, coordinates make no sense */
1894: if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue;
1895: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
1896: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
1897: PetscCall(DMPlexTransformSetDM(reftr, refdm));
1898: PetscCall(DMPlexTransformGetType(tr, &typeName));
1899: PetscCall(DMPlexTransformSetType(reftr, typeName));
1900: PetscCall(DMPlexTransformSetUp(reftr));
1901: PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));
1903: PetscCall(DMGetDimension(trdm, &trdim));
1904: PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
1905: tr->trNv[ct] = vEnd - vStart;
1906: PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
1907: PetscCall(VecGetLocalSize(coordinates, &Nc));
1908: 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);
1909: PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
1910: PetscCall(VecGetArrayRead(coordinates, &coords));
1911: PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
1912: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1914: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
1915: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1916: for (n = 0; n < Nct; ++n) {
1917: /* Since points are 0-dimensional, coordinates make no sense */
1918: if (rct[n] == DM_POLYTOPE_POINT) continue;
1919: PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
1920: for (r = 0; r < rsize[n]; ++r) {
1921: PetscInt *closure = NULL;
1922: PetscInt clSize, cl, Nv = 0;
1924: PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
1925: PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
1926: PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1927: for (cl = 0; cl < clSize * 2; cl += 2) {
1928: const PetscInt sv = closure[cl];
1930: if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1931: }
1932: PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1933: 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]);
1934: }
1935: }
1936: if (debug) {
1937: DMPolytopeType *rct;
1938: PetscInt *rsize, *rcone, *rornt;
1939: PetscInt v, dE = trdim, d, off = 0;
1941: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
1942: for (v = 0; v < tr->trNv[ct]; ++v) {
1943: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1944: for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
1945: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1946: }
1948: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1949: for (n = 0; n < Nct; ++n) {
1950: if (rct[n] == DM_POLYTOPE_POINT) continue;
1951: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
1952: for (r = 0; r < rsize[n]; ++r) {
1953: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1954: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
1955: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1956: }
1957: }
1958: }
1959: PetscCall(DMDestroy(&refdm));
1960: PetscCall(DMDestroy(&trdm));
1961: PetscCall(DMPlexTransformDestroy(&reftr));
1962: }
1963: PetscFunctionReturn(PETSC_SUCCESS);
1964: }
1966: /*@C
1967: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
1969: Input Parameters:
1970: + tr - The `DMPlexTransform` object
1971: - ct - The cell type
1973: Output Parameters:
1974: + Nv - The number of transformed vertices in the closure of the reference cell of given type
1975: - trVerts - The coordinates of these vertices in the reference cell
1977: Level: developer
1979: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
1980: @*/
1981: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1982: {
1983: PetscFunctionBegin;
1984: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1985: if (Nv) *Nv = tr->trNv[ct];
1986: if (trVerts) *trVerts = tr->trVerts[ct];
1987: PetscFunctionReturn(PETSC_SUCCESS);
1988: }
1990: /*@C
1991: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
1993: Input Parameters:
1994: + tr - The `DMPlexTransform` object
1995: . ct - The cell type
1996: . rct - The subcell type
1997: - r - The subcell index
1999: Output Parameter:
2000: . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()`
2002: Level: developer
2004: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
2005: @*/
2006: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
2007: {
2008: PetscFunctionBegin;
2009: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
2010: PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
2011: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
2012: PetscFunctionReturn(PETSC_SUCCESS);
2013: }
2015: /* Computes new vertex as the barycenter, or centroid */
2016: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
2017: {
2018: PetscInt v, d;
2020: PetscFunctionBeginHot;
2021: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
2022: for (d = 0; d < dE; ++d) out[d] = 0.0;
2023: for (v = 0; v < Nv; ++v)
2024: for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
2025: for (d = 0; d < dE; ++d) out[d] /= Nv;
2026: PetscFunctionReturn(PETSC_SUCCESS);
2027: }
2029: /*@
2030: DMPlexTransformMapCoordinates - Calculate new coordinates for produced points
2032: Not collective
2034: Input Parameters:
2035: + tr - The `DMPlexTransform`
2036: . pct - The cell type of the parent, from whom the new cell is being produced
2037: . ct - The type being produced
2038: . p - The original point
2039: . r - The replica number requested for the produced cell type
2040: . Nv - Number of vertices in the closure of the parent cell
2041: . dE - Spatial dimension
2042: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
2044: Output Parameter:
2045: . out - The coordinates of the new vertices
2047: Level: intermediate
2049: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
2050: @*/
2051: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
2052: {
2053: PetscFunctionBeginHot;
2054: if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
2055: PetscFunctionReturn(PETSC_SUCCESS);
2056: }
2058: /*
2059: DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label
2061: Not Collective
2063: Input Parameters:
2064: + tr - The `DMPlexTransform`
2065: . label - The label in the transformed mesh
2066: . pp - The parent point in the original mesh
2067: . pct - The cell type of the parent point
2068: . p - The point in the transformed mesh
2069: . ct - The cell type of the point
2070: . r - The replica number of the point
2071: - val - The label value of the parent point
2073: Level: developer
2075: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
2076: */
2077: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
2078: {
2079: PetscFunctionBeginHot;
2080: if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
2081: PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
2082: PetscFunctionReturn(PETSC_SUCCESS);
2083: }
2085: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
2086: {
2087: DM dm;
2088: IS valueIS;
2089: const PetscInt *values;
2090: PetscInt defVal, Nv, val;
2092: PetscFunctionBegin;
2093: PetscCall(DMPlexTransformGetDM(tr, &dm));
2094: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2095: PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
2096: PetscCall(DMLabelGetValueIS(label, &valueIS));
2097: PetscCall(ISGetLocalSize(valueIS, &Nv));
2098: PetscCall(ISGetIndices(valueIS, &values));
2099: for (val = 0; val < Nv; ++val) {
2100: IS pointIS;
2101: const PetscInt *points;
2102: PetscInt numPoints, p;
2104: /* Ensure refined label is created with same number of strata as
2105: * original (even if no entries here). */
2106: PetscCall(DMLabelAddStratum(labelNew, values[val]));
2107: PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
2108: PetscCall(ISGetLocalSize(pointIS, &numPoints));
2109: PetscCall(ISGetIndices(pointIS, &points));
2110: for (p = 0; p < numPoints; ++p) {
2111: const PetscInt point = points[p];
2112: DMPolytopeType ct;
2113: DMPolytopeType *rct;
2114: PetscInt *rsize, *rcone, *rornt;
2115: PetscInt Nct, n, r, pNew = 0;
2117: PetscCall(DMPlexGetCellType(dm, point, &ct));
2118: PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2119: for (n = 0; n < Nct; ++n) {
2120: for (r = 0; r < rsize[n]; ++r) {
2121: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
2122: PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
2123: }
2124: }
2125: }
2126: PetscCall(ISRestoreIndices(pointIS, &points));
2127: PetscCall(ISDestroy(&pointIS));
2128: }
2129: PetscCall(ISRestoreIndices(valueIS, &values));
2130: PetscCall(ISDestroy(&valueIS));
2131: PetscFunctionReturn(PETSC_SUCCESS);
2132: }
2134: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
2135: {
2136: DM dm;
2137: PetscInt numLabels, l;
2139: PetscFunctionBegin;
2140: PetscCall(DMPlexTransformGetDM(tr, &dm));
2141: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0));
2142: PetscCall(DMGetNumLabels(dm, &numLabels));
2143: for (l = 0; l < numLabels; ++l) {
2144: DMLabel label, labelNew;
2145: const char *lname;
2146: PetscBool isDepth, isCellType;
2148: PetscCall(DMGetLabelName(dm, l, &lname));
2149: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
2150: if (isDepth) continue;
2151: PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
2152: if (isCellType) continue;
2153: PetscCall(DMCreateLabel(rdm, lname));
2154: PetscCall(DMGetLabel(dm, lname, &label));
2155: PetscCall(DMGetLabel(rdm, lname, &labelNew));
2156: PetscCall(RefineLabel_Internal(tr, label, labelNew));
2157: }
2158: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0));
2159: PetscFunctionReturn(PETSC_SUCCESS);
2160: }
2162: /*@
2163: DMPlexTransformCreateDiscLabels - Refine the labels which define field and discrete system regions on the transformed `DM`
2165: Not Collective
2167: Input Parameters:
2168: + tr - The `DMPlexTransform`
2169: - rdm - The refined `DM` produced by the transform
2171: Level: developer
2173: Note:
2174: Region labels attached to fields (see `DMSetField()`) and to discrete systems (see `DMSetRegionNumDS()`) are
2175: not automatically included in the list of `DM` labels, so this routine walks each field and each `PetscDS`
2176: and updates the labels to refer to the refined points.
2178: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMSetField()`, `DMSetRegionNumDS()`
2179: @*/
2180: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
2181: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
2182: {
2183: DM dm;
2184: PetscInt Nf, f, Nds, s;
2186: PetscFunctionBegin;
2187: PetscCall(DMPlexTransformGetDM(tr, &dm));
2188: PetscCall(DMGetNumFields(dm, &Nf));
2189: for (f = 0; f < Nf; ++f) {
2190: DMLabel label, labelNew;
2191: PetscObject obj;
2192: const char *lname;
2194: PetscCall(DMGetField(rdm, f, &label, &obj));
2195: if (!label) continue;
2196: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2197: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
2198: PetscCall(RefineLabel_Internal(tr, label, labelNew));
2199: PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
2200: PetscCall(DMLabelDestroy(&labelNew));
2201: }
2202: PetscCall(DMGetNumDS(dm, &Nds));
2203: for (s = 0; s < Nds; ++s) {
2204: DMLabel label, labelNew;
2205: const char *lname;
2207: PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
2208: if (!label) continue;
2209: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2210: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
2211: PetscCall(RefineLabel_Internal(tr, label, labelNew));
2212: PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
2213: PetscCall(DMLabelDestroy(&labelNew));
2214: }
2215: PetscFunctionReturn(PETSC_SUCCESS);
2216: }
2218: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
2219: {
2220: DM dm;
2221: PetscSF sf, sfNew;
2222: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
2223: const PetscInt *localPoints;
2224: const PetscSFNode *remotePoints;
2225: PetscInt *localPointsNew;
2226: PetscSFNode *remotePointsNew;
2227: PetscInt pStartNew, pEndNew, pNew;
2228: /* Brute force algorithm */
2229: PetscSF rsf;
2230: PetscSection s;
2231: const PetscInt *rootdegree;
2232: PetscInt *rootPointsNew, *remoteOffsets;
2233: PetscInt numPointsNew, pStart, pEnd, p;
2235: PetscFunctionBegin;
2236: PetscCall(DMPlexTransformGetDM(tr, &dm));
2237: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
2238: PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
2239: PetscCall(DMGetPointSF(dm, &sf));
2240: PetscCall(DMGetPointSF(rdm, &sfNew));
2241: /* Calculate size of new SF */
2242: PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
2243: if (numRoots < 0) {
2244: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
2245: PetscFunctionReturn(PETSC_SUCCESS);
2246: }
2247: for (l = 0; l < numLeaves; ++l) {
2248: const PetscInt p = localPoints[l];
2249: DMPolytopeType ct;
2250: DMPolytopeType *rct;
2251: PetscInt *rsize, *rcone, *rornt;
2252: PetscInt Nct, n;
2254: PetscCall(DMPlexGetCellType(dm, p, &ct));
2255: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2256: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
2257: }
2258: /* Send new root point numbers
2259: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
2260: */
2261: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2262: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
2263: PetscCall(PetscSectionSetChart(s, pStart, pEnd));
2264: for (p = pStart; p < pEnd; ++p) {
2265: DMPolytopeType ct;
2266: DMPolytopeType *rct;
2267: PetscInt *rsize, *rcone, *rornt;
2268: PetscInt Nct, n;
2270: PetscCall(DMPlexGetCellType(dm, p, &ct));
2271: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2272: for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
2273: }
2274: PetscCall(PetscSectionSetUp(s));
2275: PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
2276: PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
2277: PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
2278: PetscCall(PetscFree(remoteOffsets));
2279: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
2280: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
2281: PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
2282: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
2283: for (p = pStart; p < pEnd; ++p) {
2284: DMPolytopeType ct;
2285: DMPolytopeType *rct;
2286: PetscInt *rsize, *rcone, *rornt;
2287: PetscInt Nct, n, r, off;
2289: if (!rootdegree[p - pStart]) continue;
2290: PetscCall(PetscSectionGetOffset(s, p, &off));
2291: PetscCall(DMPlexGetCellType(dm, p, &ct));
2292: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2293: for (n = 0, m = 0; n < Nct; ++n) {
2294: for (r = 0; r < rsize[n]; ++r, ++m) {
2295: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2296: rootPointsNew[off + m] = pNew;
2297: }
2298: }
2299: }
2300: PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
2301: PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
2302: PetscCall(PetscSFDestroy(&rsf));
2303: PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
2304: PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
2305: for (l = 0, m = 0; l < numLeaves; ++l) {
2306: const PetscInt p = localPoints[l];
2307: DMPolytopeType ct;
2308: DMPolytopeType *rct;
2309: PetscInt *rsize, *rcone, *rornt;
2310: PetscInt Nct, n, r, q, off;
2312: PetscCall(PetscSectionGetOffset(s, p, &off));
2313: PetscCall(DMPlexGetCellType(dm, p, &ct));
2314: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2315: for (n = 0, q = 0; n < Nct; ++n) {
2316: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
2317: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2318: localPointsNew[m] = pNew;
2319: remotePointsNew[m].index = rootPointsNew[off + q];
2320: remotePointsNew[m].rank = remotePoints[l].rank;
2321: }
2322: }
2323: }
2324: PetscCall(PetscSectionDestroy(&s));
2325: PetscCall(PetscFree(rootPointsNew));
2326: /* SF needs sorted leaves to correctly calculate Gather */
2327: {
2328: PetscSFNode *rp, *rtmp;
2329: PetscInt *lp, *idx, *ltmp, i;
2331: PetscCall(PetscMalloc1(numLeavesNew, &idx));
2332: PetscCall(PetscMalloc1(numLeavesNew, &lp));
2333: PetscCall(PetscMalloc1(numLeavesNew, &rp));
2334: for (i = 0; i < numLeavesNew; ++i) {
2335: 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);
2336: idx[i] = i;
2337: }
2338: PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
2339: for (i = 0; i < numLeavesNew; ++i) {
2340: lp[i] = localPointsNew[idx[i]];
2341: rp[i] = remotePointsNew[idx[i]];
2342: }
2343: ltmp = localPointsNew;
2344: localPointsNew = lp;
2345: rtmp = remotePointsNew;
2346: remotePointsNew = rp;
2347: PetscCall(PetscFree(idx));
2348: PetscCall(PetscFree(ltmp));
2349: PetscCall(PetscFree(rtmp));
2350: }
2351: PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
2352: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
2353: PetscFunctionReturn(PETSC_SUCCESS);
2354: }
2356: /*
2357: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`.
2359: Not Collective
2361: Input Parameters:
2362: + tr - The `DMPlexTransform`
2363: . ct - The type of the parent cell
2364: . rct - The type of the produced cell
2365: . r - The index of the produced cell
2366: - x - The localized coordinates for the parent cell
2368: Output Parameter:
2369: . xr - The localized coordinates for the produced cell
2371: Level: developer
2373: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
2374: */
2375: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
2376: {
2377: PetscFE fe = NULL;
2378: PetscInt cdim, v, *subcellV;
2380: PetscFunctionBegin;
2381: PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
2382: PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
2383: PetscCall(PetscFEGetNumComponents(fe, &cdim));
2384: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
2385: PetscFunctionReturn(PETSC_SUCCESS);
2386: }
2388: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
2389: {
2390: DM dm, cdm, cdmCell, cdmNew, cdmCellNew;
2391: PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
2392: Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
2393: const PetscScalar *coords;
2394: PetscScalar *coordsNew;
2395: const PetscReal *maxCell, *Lstart, *L;
2396: PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE, sparseLocalize;
2397: PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;
2399: PetscFunctionBegin;
2400: // Need to clear the DMField for coordinates
2401: PetscCall(DMSetCoordinateField(rdm, NULL));
2402: PetscCall(DMPlexTransformGetDM(tr, &dm));
2403: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0));
2404: PetscCall(DMGetCoordinateDM(dm, &cdm));
2405: PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
2406: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
2407: PetscCall(DMGetSparseLocalize(dm, &sparseLocalize));
2408: PetscCall(DMSetSparseLocalize(rdm, sparseLocalize));
2409: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
2410: if (localized) {
2411: /* Localize coordinates of new vertices */
2412: localizeVertices = PETSC_TRUE;
2413: /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
2414: if (!maxCell) localizeCells = PETSC_TRUE;
2415: }
2416: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2417: PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
2418: PetscCall(DMGetCoordinateDim(rdm, &dE));
2419: if (maxCell) {
2420: PetscReal *LstartNew, *LNew, *maxCellNew;
2422: PetscCall(PetscMalloc3(dE, &LstartNew, dE, &LNew, dE, &maxCellNew));
2423: for (d = 0; d < dEo; ++d) {
2424: LstartNew[d] = Lstart[d];
2425: LNew[d] = L[d];
2426: maxCellNew[d] = maxCell[d] / tr->redFactor;
2427: }
2428: for (d = dEo; d < dE; ++d) {
2429: LstartNew[d] = 0.;
2430: LNew[d] = -1.;
2431: maxCellNew[d] = -1.;
2432: }
2433: PetscCall(DMSetPeriodicity(rdm, maxCellNew, LstartNew, LNew));
2434: PetscCall(PetscFree3(LstartNew, LNew, maxCellNew));
2435: }
2436: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
2437: PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
2438: PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
2439: PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
2440: PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
2441: /* Localization should be inherited */
2442: /* Stefano calculates parent cells for each new cell for localization */
2443: /* Localized cells need coordinates of closure */
2444: for (v = vStartNew; v < vEndNew; ++v) {
2445: PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
2446: PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
2447: }
2448: PetscCall(PetscSectionSetUp(coordSectionNew));
2449: PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));
2451: if (localizeCells) {
2452: PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
2453: PetscCall(DMClone(cdmNew, &cdmCellNew));
2454: PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
2455: PetscCall(DMDestroy(&cdmCellNew));
2457: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
2458: PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
2459: PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
2460: PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
2461: PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));
2463: PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
2464: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2465: for (c = cStart; c < cEnd; ++c) {
2466: PetscInt dof;
2468: PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
2469: if (dof) {
2470: DMPolytopeType ct;
2471: DMPolytopeType *rct;
2472: PetscInt *rsize, *rcone, *rornt;
2473: PetscInt dim, cNew, Nct, n, r;
2475: PetscCall(DMPlexGetCellType(dm, c, &ct));
2476: dim = DMPolytopeTypeGetDim(ct);
2477: PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2478: /* This allows for different cell types */
2479: for (n = 0; n < Nct; ++n) {
2480: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2481: for (r = 0; r < rsize[n]; ++r) {
2482: PetscInt *closure = NULL;
2483: PetscInt clSize, cl, Nv = 0;
2485: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
2486: PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2487: for (cl = 0; cl < clSize * 2; cl += 2) {
2488: if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
2489: }
2490: PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2491: PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
2492: PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
2493: }
2494: }
2495: }
2496: }
2497: PetscCall(PetscSectionSetUp(coordSectionCellNew));
2498: PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
2499: }
2500: PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
2501: {
2502: VecType vtype;
2503: PetscInt coordSizeNew, bs;
2504: const char *name;
2506: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
2507: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
2508: PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
2509: PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
2510: PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
2511: PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
2512: PetscCall(VecGetBlockSize(coordsLocal, &bs));
2513: PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
2514: PetscCall(VecGetType(coordsLocal, &vtype));
2515: PetscCall(VecSetType(coordsLocalNew, vtype));
2516: }
2517: PetscCall(VecGetArrayRead(coordsLocal, &coords));
2518: PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
2519: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2520: /* First set coordinates for vertices */
2521: for (p = pStart; p < pEnd; ++p) {
2522: DMPolytopeType ct;
2523: DMPolytopeType *rct;
2524: PetscInt *rsize, *rcone, *rornt;
2525: PetscInt Nct, n, r;
2526: PetscBool hasVertex = PETSC_FALSE;
2528: PetscCall(DMPlexGetCellType(dm, p, &ct));
2529: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2530: for (n = 0; n < Nct; ++n) {
2531: if (rct[n] == DM_POLYTOPE_POINT) {
2532: hasVertex = PETSC_TRUE;
2533: break;
2534: }
2535: }
2536: if (hasVertex) {
2537: const PetscScalar *icoords = NULL;
2538: const PetscScalar *array = NULL;
2539: PetscScalar *pcoords = NULL;
2540: PetscBool isDG;
2541: PetscInt Nc, Nv, v, d;
2543: PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2545: icoords = pcoords;
2546: Nv = Nc / dEo;
2547: if (ct != DM_POLYTOPE_POINT) {
2548: if (localizeVertices && maxCell) {
2549: PetscScalar anchor[3];
2551: for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
2552: for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
2553: }
2554: }
2555: for (n = 0; n < Nct; ++n) {
2556: if (rct[n] != DM_POLYTOPE_POINT) continue;
2557: for (r = 0; r < rsize[n]; ++r) {
2558: PetscScalar vcoords[3];
2559: PetscInt vNew, off;
2561: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
2562: PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
2563: PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
2564: PetscCall(DMSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
2565: }
2566: }
2567: PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2568: }
2569: }
2570: PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
2571: PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
2572: PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
2573: PetscCall(VecDestroy(&coordsLocalNew));
2574: PetscCall(PetscSectionDestroy(&coordSectionNew));
2575: /* Then set coordinates for cells by localizing */
2576: if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
2577: else {
2578: VecType vtype;
2579: PetscInt coordSizeNew, bs;
2580: const char *name;
2582: PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
2583: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
2584: PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
2585: PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
2586: PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
2587: PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
2588: PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
2589: PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
2590: PetscCall(VecGetType(coordsLocalCell, &vtype));
2591: PetscCall(VecSetType(coordsLocalCellNew, vtype));
2592: PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
2593: PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));
2595: for (p = pStart; p < pEnd; ++p) {
2596: DMPolytopeType ct;
2597: DMPolytopeType *rct;
2598: PetscInt *rsize, *rcone, *rornt;
2599: PetscInt dof = 0, Nct, n, r;
2601: PetscCall(DMPlexGetCellType(dm, p, &ct));
2602: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2603: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
2604: if (dof) {
2605: const PetscScalar *pcoords;
2607: PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords));
2608: for (n = 0; n < Nct; ++n) {
2609: const PetscInt Nr = rsize[n];
2611: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2612: for (r = 0; r < Nr; ++r) {
2613: PetscInt pNew, offNew;
2615: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2616: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2617: cell to the ones it produces. */
2618: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2619: PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
2620: PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
2621: }
2622: }
2623: }
2624: }
2625: PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
2626: PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
2627: PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
2628: PetscCall(VecDestroy(&coordsLocalCellNew));
2629: PetscCall(PetscSectionDestroy(&coordSectionCellNew));
2630: }
2631: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0));
2632: PetscFunctionReturn(PETSC_SUCCESS);
2633: }
2635: /*@
2636: DMPlexTransformApply - Execute the transformation, producing another `DM`
2638: Collective
2640: Input Parameters:
2641: + tr - The `DMPlexTransform` object
2642: - dm - The original `DM`
2644: Output Parameter:
2645: . trdm - The transformed `DM`
2647: Level: intermediate
2649: Options Database Keys:
2650: + -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point
2651: . -dm_plex_transform_label_replica_inc num - Increment for the label value to be multiplied by the replica number
2652: - -dm_plex_transform_active name - Name for active mesh label
2654: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()`
2655: @*/
2656: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *trdm)
2657: {
2658: DM rdm;
2659: DMPlexInterpolatedFlag interp;
2660: PetscInt pStart, pEnd;
2662: PetscFunctionBegin;
2665: PetscAssertPointer(trdm, 3);
2666: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0));
2667: PetscCall(DMPlexTransformSetDM(tr, dm));
2669: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
2670: PetscCall(DMSetType(rdm, DMPLEX));
2671: PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
2672: /* Calculate number of new points of each depth */
2673: PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
2674: PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
2675: /* Step 1: Set chart */
2676: PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
2677: PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
2678: /* Step 2: Set cone/support sizes (automatically stratifies) */
2679: PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
2680: /* Step 3: Setup refined DM */
2681: PetscCall(DMSetUp(rdm));
2682: /* Step 4: Set cones and supports (automatically symmetrizes) */
2683: PetscCall(DMPlexTransformSetCones(tr, rdm));
2684: /* Step 5: Create pointSF */
2685: PetscCall(DMPlexTransformCreateSF(tr, rdm));
2686: /* Step 6: Create labels */
2687: PetscCall(DMPlexTransformCreateLabels(tr, rdm));
2688: /* Step 7: Set coordinates */
2689: PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
2690: // Do not copy periodicity, which was handled in DMPlexTransformSetCoordinates()
2691: PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_TRUE, rdm));
2692: // If the original DM was configured from options, the transformed DM should be as well
2693: rdm->setfromoptionscalled = dm->setfromoptionscalled;
2694: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0));
2695: *trdm = rdm;
2696: PetscFunctionReturn(PETSC_SUCCESS);
2697: }
2699: /*@
2700: DMPlexTransformAdaptLabel - Adapt a `DMPLEX` using a `DMPlexTransform` driven by a `DMLabel` marking cells to be refined or coarsened.
2702: Collective
2704: Input Parameters:
2705: + dm - the input `DMPLEX`
2706: . metric - unused; present to conform to the `DMAdaptor` label-based interface
2707: . adaptLabel - a `DMLabel` marking cells with `DM_ADAPT_REFINE`, `DM_ADAPT_COARSEN`, etc.
2708: - rgLabel - unused region-tag label; present to conform to the `DMAdaptor` interface
2710: Output Parameter:
2711: . rdm - the adapted `DMPLEX`
2713: Level: developer
2715: Note:
2716: This routine is registered as the "cellrefiner" adaptor by `DMGenerateRegisterAll()` and is invoked through `DMAdaptLabel()`.
2718: .seealso: `DMPLEX`, `DMPlexTransform`, `DMAdaptLabel()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`, `DMLabel`
2719: @*/
2720: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2721: {
2722: DMPlexTransform tr;
2723: DM cdm, rcdm;
2724: const char *prefix;
2725: PetscBool save;
2727: PetscFunctionBegin;
2728: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
2729: PetscCall(PetscObjectSetName((PetscObject)tr, "Adapt Label Transform"));
2730: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2731: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
2732: PetscCall(DMPlexTransformSetDM(tr, dm));
2733: PetscCall(DMPlexTransformSetFromOptions(tr));
2734: if (adaptLabel) PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
2735: PetscCall(DMPlexTransformSetUp(tr));
2736: PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
2737: PetscCall(DMPlexTransformApply(tr, dm, rdm));
2738: PetscCall(DMCopyDisc(dm, *rdm));
2739: PetscCall(DMGetCoordinateDM(dm, &cdm));
2740: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
2741: PetscCall(DMCopyDisc(cdm, rcdm));
2742: PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
2743: PetscCall(DMCopyDisc(dm, *rdm));
2744: PetscCall(DMPlexGetSaveTransform(dm, &save));
2745: if (save) PetscCall(DMPlexSetTransform(*rdm, tr));
2746: PetscCall(DMPlexTransformDestroy(&tr));
2747: ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2748: PetscFunctionReturn(PETSC_SUCCESS);
2749: }