Actual source code: plextransform.c
1: #include <petsc/private/dmplextransformimpl.h>
3: #include <petsc/private/petscfeimpl.h>
4: #include <petsc/private/hashmapi.h>
5: #include <petsc/private/dmlabelimpl.h>
7: PetscClassId DMPLEXTRANSFORM_CLASSID;
9: PetscFunctionList DMPlexTransformList = NULL;
10: PetscBool DMPlexTransformRegisterAllCalled = PETSC_FALSE;
12: PetscLogEvent DMPLEXTRANSFORM_SetUp, DMPLEXTRANSFORM_Apply, DMPLEXTRANSFORM_SetConeSizes, DMPLEXTRANSFORM_SetCones, DMPLEXTRANSFORM_CreateSF, DMPLEXTRANSFORM_CreateLabels, DMPLEXTRANSFORM_SetCoordinates, DMPLEXTRANSFORM_Check;
14: /* 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
15: OR in standard plex ordering if dm == NULL */
16: static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DM dm, PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
17: {
18: PetscInt *ctO, *ctOInv;
19: PetscInt d, c, off = 0;
20: PetscInt dimOrder[5] = {3, 2, 1, 0, -1};
22: PetscFunctionBegin;
23: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv));
24: if (dm) { // Order the dimensions by their starting location
25: PetscInt hStart[4] = {-1, -1, -1, -1};
26: for (d = 0; d <= dim; ++d) PetscCall(DMPlexGetDepthStratum(dm, dim - d, &hStart[d], NULL));
27: PetscCall(PetscSortIntWithArray(dim + 1, hStart, &dimOrder[3 - dim]));
28: } else if (dim > 1) { // Standard plex ordering. dimOrder is in correct order if dim > 1
29: off = 4 - dim;
30: dimOrder[off++] = 0;
31: for (d = dim - 1; d > 0; --d) dimOrder[off++] = d;
32: }
34: off = 0;
35: for (d = 0; d < 5; ++d) {
36: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
37: if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue;
38: if (DMPolytopeTypeGetDim((DMPolytopeType)c) == dimOrder[d]) ctO[off++] = c;
39: }
40: }
41: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
42: if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) ctO[off++] = c;
43: }
44: ctO[off++] = DM_NUM_POLYTOPES;
45: PetscCheck(off == DM_NUM_POLYTOPES + 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %" PetscInt_FMT " for cell type order", off);
47: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c;
49: *ctOrder = ctO;
50: *ctOrderInv = ctOInv;
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /*@
55: DMPlexTransformRegister - Adds a new transform component implementation
57: Not Collective
59: Input Parameters:
60: + name - The name of a new user-defined creation routine
61: - create_func - The creation routine
63: Example Usage:
64: .vb
65: DMPlexTransformRegister("my_transform", MyTransformCreate);
66: .ve
68: Then, your transform type can be chosen with the procedural interface via
69: .vb
70: DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
71: DMPlexTransformSetType(DMPlexTransform, "my_transform");
72: .ve
73: or at runtime via the option
74: .vb
75: -dm_plex_transform_type my_transform
76: .ve
78: Level: advanced
80: Note:
81: `DMPlexTransformRegister()` may be called multiple times to add several user-defined transforms
83: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()`
84: @*/
85: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
86: {
87: PetscFunctionBegin;
88: PetscCall(DMInitializePackage());
89: PetscCall(PetscFunctionListAdd(&DMPlexTransformList, name, create_func));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
94: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
95: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
96: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToSimplex(DMPlexTransform);
97: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
98: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
99: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
100: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform);
101: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform);
102: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform);
104: /*@
105: DMPlexTransformRegisterAll - Registers all of the transform components in the `DM` package.
107: Not Collective
109: Level: advanced
111: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
112: @*/
113: PetscErrorCode DMPlexTransformRegisterAll(void)
114: {
115: PetscFunctionBegin;
116: if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
117: DMPlexTransformRegisterAllCalled = PETSC_TRUE;
119: PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter));
120: PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular));
121: PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox));
122: PetscCall(DMPlexTransformRegister(DMPLEXREFINETOSIMPLEX, DMPlexTransformCreate_ToSimplex));
123: PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld));
124: PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL));
125: PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR));
126: PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D));
127: PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDETYPE, DMPlexTransformCreate_Extrude));
128: PetscCall(DMPlexTransformRegister(DMPLEXCOHESIVEEXTRUDE, DMPlexTransformCreate_Cohesive));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: /*@
133: DMPlexTransformRegisterDestroy - This function destroys the registered `DMPlexTransformType`. It is called from `PetscFinalize()`.
135: Not collective
137: Level: developer
139: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()`
140: @*/
141: PetscErrorCode DMPlexTransformRegisterDestroy(void)
142: {
143: PetscFunctionBegin;
144: PetscCall(PetscFunctionListDestroy(&DMPlexTransformList));
145: DMPlexTransformRegisterAllCalled = PETSC_FALSE;
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: DMPlexTransformCreate - Creates an empty transform object. The type can then be set with `DMPlexTransformSetType()`.
152: Collective
154: Input Parameter:
155: . comm - The communicator for the transform object
157: Output Parameter:
158: . tr - The transform object
160: Level: beginner
162: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
163: @*/
164: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
165: {
166: DMPlexTransform t;
168: PetscFunctionBegin;
169: PetscAssertPointer(tr, 2);
170: *tr = NULL;
171: PetscCall(DMInitializePackage());
173: PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView));
174: t->setupcalled = PETSC_FALSE;
175: t->redFactor = 2.0;
176: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom));
177: *tr = t;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: /*@
182: DMPlexTransformSetType - Sets the particular implementation for a transform.
184: Collective
186: Input Parameters:
187: + tr - The transform
188: - method - The name of the transform type
190: Options Database Key:
191: . -dm_plex_transform_type type - Sets the transform type; see `DMPlexTransformType`
193: Level: intermediate
195: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
196: @*/
197: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
198: {
199: PetscErrorCode (*r)(DMPlexTransform);
200: PetscBool match;
202: PetscFunctionBegin;
204: PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match));
205: if (match) PetscFunctionReturn(PETSC_SUCCESS);
207: PetscCall(DMPlexTransformRegisterAll());
208: PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r));
209: PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);
211: PetscTryTypeMethod(tr, destroy);
212: PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops)));
213: PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method));
214: PetscCall((*r)(tr));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*@
219: DMPlexTransformGetType - Gets the type name (as a string) from the transform.
221: Not Collective
223: Input Parameter:
224: . tr - The `DMPlexTransform`
226: Output Parameter:
227: . type - The `DMPlexTransformType` name
229: Level: intermediate
231: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
232: @*/
233: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
234: {
235: PetscFunctionBegin;
237: PetscAssertPointer(type, 2);
238: PetscCall(DMPlexTransformRegisterAll());
239: *type = ((PetscObject)tr)->type_name;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
244: {
245: PetscViewerFormat format;
247: PetscFunctionBegin;
248: PetscCall(PetscViewerGetFormat(v, &format));
249: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
250: const PetscInt *trTypes = NULL;
251: IS trIS;
252: PetscInt cols = 8;
253: PetscInt Nrt = 8, f, g;
254: PetscMPIInt size, rank;
256: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)tr), &rank));
257: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)tr), &size));
258: PetscCall(PetscViewerASCIIPushSynchronized(v));
259: if (tr->trType) PetscCall(DMLabelView(tr->trType, v));
260: if (size > 1) PetscCall(PetscViewerASCIISynchronizedPrintf(v, "Process: %d\n", rank));
261: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "Source Starts\n"));
262: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIISynchronizedPrintf(v, " %14s", DMPolytopeTypes[g]));
263: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "\n"));
264: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIISynchronizedPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]));
265: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "\n"));
266: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "Target Starts\n"));
267: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIISynchronizedPrintf(v, " %14s", DMPolytopeTypes[g]));
268: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "\n"));
269: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIISynchronizedPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]));
270: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "\n"));
272: if (tr->trType) {
273: PetscCall(DMLabelGetNumValues(tr->trType, &Nrt));
274: PetscCall(DMLabelGetValueIS(tr->trType, &trIS));
275: PetscCall(ISGetIndices(trIS, &trTypes));
276: }
277: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "Offsets\n"));
278: PetscCall(PetscViewerASCIISynchronizedPrintf(v, " "));
279: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIISynchronizedPrintf(v, " %14s", DMPolytopeTypes[g]));
280: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "\n"));
281: for (f = 0; f < Nrt; ++f) {
282: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "%2" PetscInt_FMT " |", trTypes ? trTypes[f] : f));
283: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIISynchronizedPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]));
284: PetscCall(PetscViewerASCIISynchronizedPrintf(v, " |\n"));
285: }
286: if (tr->trType) {
287: PetscCall(ISRestoreIndices(trIS, &trTypes));
288: PetscCall(ISDestroy(&trIS));
289: }
290: PetscCall(PetscViewerFlush(v));
291: }
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: /*@
296: DMPlexTransformView - Views a `DMPlexTransform`
298: Collective
300: Input Parameters:
301: + tr - the `DMPlexTransform` object to view
302: - v - the viewer
304: Level: beginner
306: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
307: @*/
308: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
309: {
310: PetscBool isascii;
312: PetscFunctionBegin;
314: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v));
316: PetscCheckSameComm(tr, 1, v, 2);
317: PetscCall(PetscViewerCheckWritable(v));
318: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v));
319: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
320: if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v));
321: PetscTryTypeMethod(tr, view, v);
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*@
326: DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database
328: Collective
330: Input Parameter:
331: . tr - the `DMPlexTransform` object to set options for
333: Options Database Keys:
334: + -dm_plex_transform_type type - Set the transform type, e.g. refine_regular
335: . -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point
336: . -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
337: . -dm_plex_transform_active name - Name for active mesh label
338: - -dm_plex_transform_active_values v0,v1,... - Values in the active label
340: Level: intermediate
342: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
343: @*/
344: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
345: {
346: char typeName[1024], active[PETSC_MAX_PATH_LEN];
347: const char *defName = DMPLEXREFINEREGULAR;
348: PetscBool flg, match;
350: PetscFunctionBegin;
352: PetscObjectOptionsBegin((PetscObject)tr);
353: PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, sizeof(typeName), &flg));
354: if (flg) PetscCall(DMPlexTransformSetType(tr, typeName));
355: else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName));
356: PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &match, &flg));
357: if (flg) PetscCall(DMPlexTransformSetMatchStrata(tr, match));
358: 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));
359: PetscCall(PetscOptionsString("-dm_plex_transform_active", "Name for active mesh label", "DMPlexTransformSetActive", active, active, sizeof(active), &flg));
360: if (flg) {
361: DM dm;
362: DMLabel label;
363: PetscInt values[16];
364: PetscInt n = 16;
366: PetscCall(DMPlexTransformGetDM(tr, &dm));
367: PetscCall(DMGetLabel(dm, active, &label));
368: PetscCall(PetscOptionsIntArray("-dm_plex_transform_active_values", "The label values to be active", "DMPlexTransformSetActive", values, &n, &flg));
369: if (flg && n) {
370: DMLabel newlabel;
372: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Active", &newlabel));
373: for (PetscInt i = 0; i < n; ++i) {
374: IS is;
376: PetscCall(DMLabelGetStratumIS(label, values[i], &is));
377: PetscCall(DMLabelInsertIS(newlabel, is, values[i]));
378: PetscCall(ISDestroy(&is));
379: }
380: PetscCall(DMPlexTransformSetActive(tr, newlabel));
381: PetscCall(DMLabelDestroy(&newlabel));
382: } else {
383: PetscCall(DMPlexTransformSetActive(tr, label));
384: }
385: }
386: PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
387: /* process any options handlers added with PetscObjectAddOptionsHandler() */
388: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject));
389: PetscOptionsEnd();
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: /*@
394: DMPlexTransformDestroy - Destroys a `DMPlexTransform`
396: Collective
398: Input Parameter:
399: . tr - the transform object to destroy
401: Level: beginner
403: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
404: @*/
405: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
406: {
407: PetscInt c;
409: PetscFunctionBegin;
410: if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
412: if (--((PetscObject)*tr)->refct > 0) {
413: *tr = NULL;
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: PetscTryTypeMethod(*tr, destroy);
418: PetscCall(DMDestroy(&(*tr)->dm));
419: PetscCall(DMLabelDestroy(&(*tr)->active));
420: PetscCall(DMLabelDestroy(&(*tr)->trType));
421: PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld));
422: PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew));
423: PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew));
424: PetscCall(PetscFree((*tr)->offset));
425: PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd));
426: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
427: PetscCall(PetscFEDestroy(&(*tr)->coordFE[c]));
428: PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c]));
429: }
430: if ((*tr)->trVerts) {
431: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
432: DMPolytopeType *rct;
433: PetscInt *rsize, *rcone, *rornt, Nct, n, r;
435: if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) {
436: PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
437: for (n = 0; n < Nct; ++n) {
438: if (rct[n] == DM_POLYTOPE_POINT) continue;
439: for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r]));
440: PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]]));
441: }
442: }
443: PetscCall(PetscFree((*tr)->trSubVerts[c]));
444: PetscCall(PetscFree((*tr)->trVerts[c]));
445: }
446: }
447: PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts));
448: PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom));
449: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
450: PetscCall(PetscHeaderDestroy(tr));
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
454: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
455: {
456: DMLabel trType = tr->trType;
457: PetscInt c, cN, *off;
459: PetscFunctionBegin;
460: if (trType) {
461: DM dm;
462: IS rtIS;
463: const PetscInt *reftypes;
464: PetscInt Nrt;
466: PetscCall(DMPlexTransformGetDM(tr, &dm));
467: PetscCall(DMLabelGetNumValues(trType, &Nrt));
468: PetscCall(DMLabelGetValueIS(trType, &rtIS));
469: PetscCall(ISGetIndices(rtIS, &reftypes));
470: PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off));
471: for (PetscInt r = 0; r < Nrt; ++r) {
472: const PetscInt rt = reftypes[r];
473: IS rtIS;
474: const PetscInt *points;
475: DMPolytopeType ct;
476: PetscInt np, p;
478: PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
479: PetscCall(ISGetLocalSize(rtIS, &np));
480: PetscCall(ISGetIndices(rtIS, &points));
481: if (!np) continue;
482: p = points[0];
483: PetscCall(ISRestoreIndices(rtIS, &points));
484: PetscCall(ISDestroy(&rtIS));
485: PetscCall(DMPlexGetCellType(dm, p, &ct));
486: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
487: const DMPolytopeType ctNew = (DMPolytopeType)cN;
488: DMPolytopeType *rct;
489: PetscInt *rsize, *cone, *ornt;
490: PetscInt Nct, n, s;
492: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
493: off[r * DM_NUM_POLYTOPES + ctNew] = -1;
494: break;
495: }
496: off[r * DM_NUM_POLYTOPES + ctNew] = 0;
497: for (s = 0; s <= r; ++s) {
498: const PetscInt st = reftypes[s];
499: DMPolytopeType sct;
500: PetscInt q, qrt;
502: PetscCall(DMLabelGetStratumIS(trType, st, &rtIS));
503: PetscCall(ISGetLocalSize(rtIS, &np));
504: PetscCall(ISGetIndices(rtIS, &points));
505: if (!np) continue;
506: q = points[0];
507: PetscCall(ISRestoreIndices(rtIS, &points));
508: PetscCall(ISDestroy(&rtIS));
509: PetscCall(DMPlexGetCellType(dm, q, &sct));
510: PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt));
511: 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);
512: if (st == rt) {
513: for (n = 0; n < Nct; ++n)
514: if (rct[n] == ctNew) break;
515: if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
516: break;
517: }
518: for (n = 0; n < Nct; ++n) {
519: if (rct[n] == ctNew) {
520: PetscInt sn;
522: PetscCall(DMLabelGetStratumSize(trType, st, &sn));
523: off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
524: }
525: }
526: }
527: }
528: }
529: PetscCall(ISRestoreIndices(rtIS, &reftypes));
530: PetscCall(ISDestroy(&rtIS));
531: } else {
532: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off));
533: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
534: const DMPolytopeType ct = (DMPolytopeType)c;
535: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
536: const DMPolytopeType ctNew = (DMPolytopeType)cN;
537: DMPolytopeType *rct;
538: PetscInt *rsize, *cone, *ornt;
539: PetscInt Nct, n, i;
541: 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) {
542: off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
543: continue;
544: }
545: off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
546: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
547: const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i];
548: const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];
550: PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt));
551: if (ict == ct) {
552: for (n = 0; n < Nct; ++n)
553: if (rct[n] == ctNew) break;
554: if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
555: break;
556: }
557: for (n = 0; n < Nct; ++n)
558: if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
559: }
560: }
561: }
562: }
563: *offset = off;
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: /*@
568: DMPlexTransformSetUp - Create the tables that drive the transform
570: Input Parameter:
571: . tr - The `DMPlexTransform` object
573: Level: intermediate
575: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
576: @*/
577: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
578: {
579: DMPolytopeType ctCell;
580: DM dm;
581: PetscInt pStart, pEnd, p, c, celldim = 0;
583: PetscFunctionBegin;
585: if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
586: PetscCall(DMPlexTransformGetDM(tr, &dm));
587: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0));
588: PetscTryTypeMethod(tr, setup);
589: PetscCall(DMSetSnapToGeomModel(dm, NULL));
590: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
592: if (pEnd > pStart) {
593: // Ignore cells hanging off of embedded surfaces
594: PetscInt c = pStart;
596: ctCell = DM_POLYTOPE_FV_GHOST;
597: while (DMPolytopeTypeGetDim(ctCell) < 0) PetscCall(DMPlexGetCellType(dm, c++, &ctCell));
598: } else {
599: PetscInt dim;
601: PetscCall(DMGetDimension(dm, &dim));
602: switch (dim) {
603: case 0:
604: ctCell = DM_POLYTOPE_POINT;
605: break;
606: case 1:
607: ctCell = DM_POLYTOPE_SEGMENT;
608: break;
609: case 2:
610: ctCell = DM_POLYTOPE_TRIANGLE;
611: break;
612: case 3:
613: ctCell = DM_POLYTOPE_TETRAHEDRON;
614: break;
615: default:
616: ctCell = DM_POLYTOPE_UNKNOWN;
617: }
618: }
619: PetscCall(DMPlexCreateCellTypeOrder_Internal(dm, DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
620: for (p = pStart; p < pEnd; ++p) {
621: DMPolytopeType ct;
622: DMPolytopeType *rct;
623: PetscInt *rsize, *cone, *ornt;
624: PetscInt Nct;
626: PetscCall(DMPlexGetCellType(dm, p, &ct));
627: 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);
628: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
629: for (PetscInt n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
630: }
631: PetscCall(DMPlexCreateCellTypeOrder_Internal(NULL, celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
632: /* Construct sizes and offsets for each cell type */
633: if (!tr->ctStart) {
634: PetscInt *ctS, *ctSN, *ctC, *ctCN;
636: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
637: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
638: for (p = pStart; p < pEnd; ++p) {
639: DMPolytopeType ct;
640: DMPolytopeType *rct;
641: PetscInt *rsize, *cone, *ornt;
642: PetscInt Nct;
644: PetscCall(DMPlexGetCellType(dm, p, &ct));
645: 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);
646: ++ctC[ct];
647: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
648: for (PetscInt n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
649: }
650: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
651: const PetscInt cto = tr->ctOrderOld[c];
652: const PetscInt cton = tr->ctOrderOld[c + 1];
653: const PetscInt ctn = tr->ctOrderNew[c];
654: const PetscInt ctnn = tr->ctOrderNew[c + 1];
656: ctS[cton] = ctS[cto] + ctC[cto];
657: ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
658: }
659: PetscCall(PetscFree2(ctC, ctCN));
660: tr->ctStart = ctS;
661: tr->ctStartNew = ctSN;
662: }
663: PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
664: // Compute depth information
665: tr->depth = -1;
666: for (c = 0; c < DM_NUM_POLYTOPES; ++c)
667: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
668: PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
669: for (PetscInt d = 0; d <= tr->depth; ++d) {
670: tr->depthStart[d] = PETSC_INT_MAX;
671: tr->depthEnd[d] = -1;
672: }
673: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
674: const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);
676: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
677: tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
678: tr->depthEnd[dep] = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
679: }
680: tr->setupcalled = PETSC_TRUE;
681: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0));
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: /*@
686: DMPlexTransformGetDM - Get the base `DM` for the transform
688: Input Parameter:
689: . tr - The `DMPlexTransform` object
691: Output Parameter:
692: . dm - The original `DM` which will be transformed
694: Level: intermediate
696: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
697: @*/
698: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
699: {
700: PetscFunctionBegin;
702: PetscAssertPointer(dm, 2);
703: *dm = tr->dm;
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: /*@
708: DMPlexTransformSetDM - Set the base `DM` for the transform
710: Input Parameters:
711: + tr - The `DMPlexTransform` object
712: - dm - The original `DM` which will be transformed
714: Level: intermediate
716: Note:
717: The user does not typically call this, as it is called by `DMPlexTransformApply()`.
719: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
720: @*/
721: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
722: {
723: PetscFunctionBegin;
725: if (dm) {
727: PetscCall(PetscObjectReference((PetscObject)dm));
728: }
729: PetscCall(DMDestroy(&tr->dm));
730: tr->dm = dm;
731: PetscFunctionReturn(PETSC_SUCCESS);
732: }
734: /*@
735: DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform
737: Input Parameter:
738: . tr - The `DMPlexTransform` object
740: Output Parameter:
741: . active - The `DMLabel` indicating which points will be transformed
743: Level: intermediate
745: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
746: @*/
747: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
748: {
749: PetscFunctionBegin;
751: PetscAssertPointer(active, 2);
752: *active = tr->active;
753: PetscFunctionReturn(PETSC_SUCCESS);
754: }
756: /*@
757: DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform
759: Input Parameters:
760: + tr - The `DMPlexTransform` object
761: - active - The `DMLabel` indicating which points will be transformed
763: Level: intermediate
765: Note:
766: This only applies to transforms listed in [](plex_transform_table) that operate on a subset of the mesh.
768: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
769: @*/
770: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
771: {
772: PetscFunctionBegin;
775: PetscCall(PetscObjectReference((PetscObject)active));
776: PetscCall(DMLabelDestroy(&tr->active));
777: tr->active = active;
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: /*@
782: DMPlexTransformGetTransformTypes - Get the `DMLabel` marking the transform type of each point for the transform
784: Input Parameter:
785: . tr - The `DMPlexTransform` object
787: Output Parameter:
788: . trType - The `DMLabel` indicating the transform type for each point
790: Level: intermediate
792: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetTransformType()`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
793: @*/
794: PetscErrorCode DMPlexTransformGetTransformTypes(DMPlexTransform tr, DMLabel *trType)
795: {
796: PetscFunctionBegin;
798: PetscAssertPointer(trType, 2);
799: *trType = tr->trType;
800: PetscFunctionReturn(PETSC_SUCCESS);
801: }
803: /*@
804: DMPlexTransformSetTransformTypes - Set the `DMLabel` marking the transform type of each point for the transform
806: Input Parameters:
807: + tr - The `DMPlexTransform` object
808: - trType - The original `DM` which will be transformed
810: Level: intermediate
812: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetTransformTypes()`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
813: @*/
814: PetscErrorCode DMPlexTransformSetTransformTypes(DMPlexTransform tr, DMLabel trType)
815: {
816: PetscFunctionBegin;
819: PetscCall(PetscObjectReference((PetscObject)trType));
820: PetscCall(DMLabelDestroy(&tr->trType));
821: tr->trType = trType;
822: PetscFunctionReturn(PETSC_SUCCESS);
823: }
825: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
826: {
827: PetscFunctionBegin;
828: if (!tr->coordFE[ct]) {
829: PetscInt dim, cdim;
831: dim = DMPolytopeTypeGetDim(ct);
832: PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
833: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
834: {
835: PetscDualSpace dsp;
836: PetscQuadrature quad;
837: DM K;
838: PetscFEGeom *cg;
839: PetscScalar *Xq;
840: PetscReal *xq, *wq;
841: PetscInt Nq;
843: PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
844: PetscCall(PetscMalloc1(Nq * cdim, &xq));
845: for (PetscInt q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
846: PetscCall(PetscMalloc1(Nq, &wq));
847: for (PetscInt q = 0; q < Nq; ++q) wq[q] = 1.0;
848: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
849: PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
850: PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));
852: PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
853: PetscCall(PetscDualSpaceGetDM(dsp, &K));
854: PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FEGEOM_BASIC, &tr->refGeom[ct]));
855: cg = tr->refGeom[ct];
856: PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
857: PetscCall(PetscQuadratureDestroy(&quad));
858: }
859: }
860: *fe = tr->coordFE[ct];
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
865: {
866: PetscInt dim, cdim;
868: PetscFunctionBegin;
869: PetscCall(DMGetDimension(dm, &dim));
870: PetscCall(DMSetDimension(tdm, dim));
871: PetscCall(DMGetCoordinateDim(dm, &cdim));
872: PetscCall(DMSetCoordinateDim(tdm, cdim));
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: /*@
877: DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`
879: Input Parameters:
880: + tr - The `DMPlexTransform` object
881: - dm - The original `DM`
883: Output Parameter:
884: . trdm - The transformed `DM`
886: Level: advanced
888: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
889: @*/
890: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM trdm)
891: {
892: PetscFunctionBegin;
893: PetscUseTypeMethod(tr, setdimensions, dm, trdm);
894: PetscFunctionReturn(PETSC_SUCCESS);
895: }
897: /*@
898: DMPlexTransformGetChart - Get the chart `[pStart, pEnd)` for the points produced by the transform
900: Not Collective
902: Input Parameter:
903: . tr - The `DMPlexTransform`
905: Output Parameters:
906: + pStart - The first point in the transformed mesh, or `NULL` if not needed
907: - pEnd - One past the last point in the transformed mesh, or `NULL` if not needed
909: Level: developer
911: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformGetCellType()`, `DMPlexTransformGetCellTypeStratum()`
912: @*/
913: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
914: {
915: PetscFunctionBegin;
916: if (pStart) *pStart = 0;
917: if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: /*@
922: DMPlexTransformGetCellType - Return the cell type for a point in the transformed mesh
924: Not Collective
926: Input Parameters:
927: + tr - The `DMPlexTransform`
928: - cell - The point number in the transformed mesh
930: Output Parameter:
931: . celltype - The `DMPolytopeType` of the point
933: Level: developer
935: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetChart()`, `DMPlexTransformGetCellTypeStratum()`
936: @*/
937: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
938: {
939: PetscInt ctNew;
941: PetscFunctionBegin;
943: PetscAssertPointer(celltype, 3);
944: /* TODO Can do bisection since everything is sorted */
945: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
946: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
948: if (cell >= ctSN && cell < ctEN) break;
949: }
950: PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
951: *celltype = (DMPolytopeType)ctNew;
952: PetscFunctionReturn(PETSC_SUCCESS);
953: }
955: /*@
956: DMPlexTransformGetCellTypeStratum - Return the point range for a given cell type in the transformed mesh
958: Not Collective
960: Input Parameters:
961: + tr - The `DMPlexTransform`
962: - celltype - The `DMPolytopeType` of the requested stratum
964: Output Parameters:
965: + start - The first point of the stratum, or `NULL` if not needed
966: - end - One past the last point of the stratum, or `NULL` if not needed
968: Level: developer
970: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellType()`, `DMPlexTransformGetChart()`, `DMPlexGetDepthStratum()`
971: @*/
972: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
973: {
974: PetscFunctionBegin;
976: if (start) *start = tr->ctStartNew[celltype];
977: if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
978: PetscFunctionReturn(PETSC_SUCCESS);
979: }
981: /*@
982: DMPlexTransformGetDepth - Return the topological depth of the transformed mesh
984: Not Collective
986: Input Parameter:
987: . tr - The `DMPlexTransform`
989: Output Parameter:
990: . depth - The depth of the transformed mesh
992: Level: developer
994: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDepthStratum()`, `DMPlexGetDepth()`
995: @*/
996: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
997: {
998: PetscFunctionBegin;
1000: *depth = tr->depth;
1001: PetscFunctionReturn(PETSC_SUCCESS);
1002: }
1004: /*@
1005: DMPlexTransformGetDepthStratum - Return the point range for a given depth in the transformed mesh
1007: Not Collective
1009: Input Parameters:
1010: + tr - The `DMPlexTransform`
1011: - depth - The requested depth in the transformed mesh
1013: Output Parameters:
1014: + start - The first point at the given depth, or `NULL` if not needed
1015: - end - One past the last point at the given depth, or `NULL` if not needed
1017: Level: developer
1019: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDepth()`, `DMPlexGetDepthStratum()`
1020: @*/
1021: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
1022: {
1023: PetscFunctionBegin;
1025: if (start) *start = tr->depthStart[depth];
1026: if (end) *end = tr->depthEnd[depth];
1027: PetscFunctionReturn(PETSC_SUCCESS);
1028: }
1030: /*@
1031: DMPlexTransformGetMatchStrata - Get the flag which determines what points get added to the transformed labels
1033: Not Collective
1035: Input Parameter:
1036: . tr - The `DMPlexTransform`
1038: Output Parameter:
1039: . match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels
1041: Level: intermediate
1043: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetMatchStrata()`, `DMPlexGetPointDepth()`
1044: @*/
1045: PetscErrorCode DMPlexTransformGetMatchStrata(DMPlexTransform tr, PetscBool *match)
1046: {
1047: PetscFunctionBegin;
1049: PetscAssertPointer(match, 2);
1050: *match = tr->labelMatchStrata;
1051: PetscFunctionReturn(PETSC_SUCCESS);
1052: }
1054: /*@
1055: DMPlexTransformSetMatchStrata - Set the flag which determines what points get added to the transformed labels
1057: Not Collective
1059: Input Parameters:
1060: + tr - The `DMPlexTransform`
1061: - match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels
1063: Level: intermediate
1065: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetMatchStrata()`, `DMPlexGetPointDepth()`
1066: @*/
1067: PetscErrorCode DMPlexTransformSetMatchStrata(DMPlexTransform tr, PetscBool match)
1068: {
1069: PetscFunctionBegin;
1071: tr->labelMatchStrata = match;
1072: PetscFunctionReturn(PETSC_SUCCESS);
1073: }
1075: /*@
1076: DMPlexTransformCheck - Verify that the given `DM`, produced by this `DMPlexTransform`, is valid
1078: Input Parameters:
1079: + tr - The `DMPlexTransform` object
1080: - dm - The `DM` to check
1082: Level: advanced
1084: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
1085: @*/
1086: PetscErrorCode DMPlexTransformCheck(DMPlexTransform tr, DM dm)
1087: {
1088: PetscFunctionBegin;
1089: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_Check, tr, dm, 0, 0));
1090: PetscTryTypeMethod(tr, check, dm);
1091: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_Check, tr, dm, 0, 0));
1092: PetscFunctionReturn(PETSC_SUCCESS);
1093: }
1095: /*@
1096: DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.
1098: Not Collective
1100: Input Parameters:
1101: + tr - The `DMPlexTransform`
1102: . ct - The type of the original point which produces the new point
1103: . ctNew - The type of the new point
1104: . p - The original point which produces the new point
1105: - r - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`
1107: Output Parameter:
1108: . pNew - The new point number
1110: Level: developer
1112: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
1113: @*/
1114: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
1115: {
1116: DMPolytopeType *rct;
1117: PetscInt *rsize, *cone, *ornt;
1118: PetscInt rt, Nct, n, off, rp;
1119: DMLabel trType = tr->trType;
1120: PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
1121: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
1122: PetscInt newp = ctSN, cind;
1124: PetscFunctionBeginHot;
1125: 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);
1126: PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
1127: if (trType) {
1128: PetscCall(DMLabelLookupStratum(trType, rt, &cind));
1129: PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
1130: 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);
1131: } else {
1132: cind = ct;
1133: rp = p - ctS;
1134: }
1135: off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
1136: 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);
1137: newp += off;
1138: for (n = 0; n < Nct; ++n) {
1139: if (rct[n] == ctNew) {
1140: PetscCheck(!rsize[n] || r < rsize[n], 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]);
1141: newp += rp * rsize[n] + r;
1142: if (!(newp >= ctSN && newp <= ctEN)) {
1143: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Problem with point %" PetscInt_FMT " %s replica %" PetscInt_FMT "\n", p, DMPolytopeTypes[ct], r));
1144: PetscCall(PetscPrintf(PETSC_COMM_SELF, " n %" PetscInt_FMT " rsize %" PetscInt_FMT " rt %" PetscInt_FMT " cind %" PetscInt_FMT " rp %" PetscInt_FMT "\n", n, rsize[n], rt, cind, rp));
1145: }
1146: break;
1147: }
1148: }
1150: 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);
1151: *pNew = newp;
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: /*@
1156: DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.
1158: Not Collective
1160: Input Parameters:
1161: + tr - The `DMPlexTransform`
1162: - pNew - The new point number
1164: Output Parameters:
1165: + ct - The type of the original point which produces the new point
1166: . ctNew - The type of the new point
1167: . p - The original point which produces the new point
1168: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
1170: Level: developer
1172: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
1173: @*/
1174: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
1175: {
1176: DMLabel trType = tr->trType;
1177: DMPolytopeType *rct, ctN;
1178: PetscInt *rsize, *cone, *ornt;
1179: PetscInt rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
1180: PetscInt offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;
1182: PetscFunctionBegin;
1183: PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
1184: if (trType) {
1185: DM dm;
1186: IS rtIS;
1187: const PetscInt *reftypes;
1188: PetscInt Nrt, r;
1190: PetscCall(DMPlexTransformGetDM(tr, &dm));
1191: PetscCall(DMLabelGetNumValues(trType, &Nrt));
1192: PetscCall(DMLabelGetValueIS(trType, &rtIS));
1193: PetscCall(ISGetIndices(rtIS, &reftypes));
1194: for (r = 0; r < Nrt; ++r) {
1195: const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];
1197: if (tr->ctStartNew[ctN] + off > pNew) continue;
1198: /* Check that any of this refinement type exist */
1199: /* TODO Actually keep track of the number produced here instead */
1200: if (off > offset) {
1201: rt = reftypes[r];
1202: offset = off;
1203: }
1204: }
1205: PetscCall(ISRestoreIndices(rtIS, &reftypes));
1206: PetscCall(ISDestroy(&rtIS));
1207: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
1208: /* TODO Map refinement types to cell types */
1209: PetscCall(DMLabelGetStratumBounds(trType, rt, &rtS, NULL));
1210: PetscCheck(rtS >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
1211: for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
1212: PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
1214: if (rtS >= ctS && rtS < ctE) break;
1215: }
1216: PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt);
1217: } else {
1218: for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
1219: const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];
1221: if (tr->ctStartNew[ctN] + off > pNew) continue;
1222: if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
1223: /* TODO Actually keep track of the number produced here instead */
1224: if (off > offset) {
1225: ctO = ctTmp;
1226: offset = off;
1227: }
1228: }
1229: rt = -1;
1230: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
1231: }
1232: ctS = tr->ctStart[ctO];
1233: ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
1234: if (!trType) rtS = ctS;
1235: PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
1236: 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);
1237: for (n = 0; n < Nct; ++n) {
1238: if (rct[n] == ctN) {
1239: PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, c;
1241: if (trType) {
1242: IS rtIS;
1243: const PetscInt *points;
1244: const PetscInt idx = tmp / rsize[n];
1245: PetscInt pStart, pEnd;
1247: PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
1248: PetscCall(ISGetPointRange(rtIS, &pStart, &pEnd, &points));
1249: PetscCheck(idx < pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could not be found due to invalid index %" PetscInt_FMT " not in range [0, %" PetscInt_FMT ")", pNew, idx, pEnd - pStart);
1250: c = points ? points[idx] : pStart + idx;
1251: PetscCall(ISRestorePointRange(rtIS, &pStart, &pEnd, &points));
1252: PetscCall(ISDestroy(&rtIS));
1253: rp = c - ctS;
1254: rO = tmp % rsize[n];
1255: } else {
1256: // This assumes that all points of type ctO transform the same way
1257: rp = tmp / rsize[n];
1258: rO = tmp % rsize[n];
1259: }
1260: break;
1261: }
1262: }
1263: PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
1264: pO = rp + ctS;
1265: 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);
1266: if (ct) *ct = (DMPolytopeType)ctO;
1267: if (ctNew) *ctNew = ctN;
1268: if (p) *p = pO;
1269: if (r) *r = rO;
1270: PetscFunctionReturn(PETSC_SUCCESS);
1271: }
1273: /*@
1274: DMPlexTransformCreateSplitCellLabel - Mark the cells of a transformed mesh whose source cell was genuinely split
1276: Not Collective
1278: Input Parameters:
1279: + tr - The `DMPlexTransform` that produced `dm`
1280: - dm - The transformed `DM`
1282: Output Parameter:
1283: . label - A `DMLabel` marking with 1 the cells of `dm` whose source cell produced more than one cell
1285: Level: advanced
1287: Notes:
1288: For `DMPLEXREFINESBR`, the label also marks cells whose source cell was not flagged in the active label,
1289: since the transform splits further cells to restore conformity.
1291: The transform is available only if saved with `DMPlexSetSaveTransform()`. Use `DMPlexGetTransform()` to get it.
1293: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformGetActive()`, `DMPlexGetTransform()`, `DMPlexLabelComplete()`
1294: @*/
1295: PetscErrorCode DMPlexTransformCreateSplitCellLabel(DMPlexTransform tr, DM dm, DMLabel *label)
1296: {
1297: PetscHMapI children;
1298: PetscInt *source;
1299: PetscInt cStart, cEnd;
1301: PetscFunctionBegin;
1304: PetscAssertPointer(label, 3);
1305: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1306: /* Count the cells each source point produced, rather than asking the transform how it refines a cell type, so that
1307: this does not have to know how any particular transform reports its refinement types */
1308: PetscCall(PetscHMapICreate(&children));
1309: PetscCall(PetscMalloc1(cEnd - cStart, &source));
1310: for (PetscInt c = cStart; c < cEnd; ++c) {
1311: PetscInt p, n = 0;
1313: PetscCall(DMPlexTransformGetSourcePoint(tr, c, NULL, NULL, &p, NULL));
1314: PetscCall(PetscHMapIGetWithDefault(children, p, 0, &n));
1315: PetscCall(PetscHMapISet(children, p, n + 1));
1316: source[c - cStart] = p;
1317: }
1318: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "split cells", label));
1319: for (PetscInt c = cStart; c < cEnd; ++c) {
1320: PetscInt n = 0;
1322: PetscCall(PetscHMapIGetWithDefault(children, source[c - cStart], 0, &n));
1323: if (n > 1) PetscCall(DMLabelSetValue(*label, c, 1));
1324: }
1325: PetscCall(PetscFree(source));
1326: PetscCall(PetscHMapIDestroy(&children));
1327: PetscFunctionReturn(PETSC_SUCCESS);
1328: }
1330: /*@
1331: DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.
1333: Input Parameters:
1334: + tr - The `DMPlexTransform` object
1335: . source - The source cell type
1336: - p - The source point, which can also determine the refine type
1338: Output Parameters:
1339: + rt - The refine type for this point
1340: . Nt - The number of types produced by this point
1341: . target - An array of length `Nt` giving the types produced
1342: . size - An array of length `Nt` giving the number of cells of each type produced
1343: . cone - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
1344: - ornt - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point
1346: Level: advanced
1348: Notes:
1349: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
1350: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
1351: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
1352: .vb
1353: the number of cones to be taken, or 0 for the current cell
1354: the cell cone point number at each level from which it is subdivided
1355: the replica number r of the subdivision.
1356: .ve
1357: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
1358: .vb
1359: Nt = 2
1360: target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
1361: size = {1, 2}
1362: cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
1363: ornt = { 0, 0, 0, 0}
1364: .ve
1366: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
1367: @*/
1368: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1369: {
1370: PetscFunctionBegin;
1371: PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
1372: PetscFunctionReturn(PETSC_SUCCESS);
1373: }
1375: /*@
1376: DMPlexTransformGetSubcellOrientationIdentity - Default `getsubcellorientation` implementation for transforms that reproduce the input mesh
1378: Not Collective
1380: Input Parameters:
1381: + tr - The `DMPlexTransform`
1382: . sct - The source point cell type
1383: . sp - The source point
1384: . so - The orientation of the source point in its enclosing parent
1385: . tct - The target point cell type
1386: . r - The replica number requested for the produced cell type
1387: - o - The orientation of the replica
1389: Output Parameters:
1390: + rnew - The replica number, given the orientation of the parent (returns `r`)
1391: - onew - The replica orientation composed with the source orientation
1393: Level: developer
1395: Note:
1396: This is the identity variant used by transforms such as the "identity" refiner where each source
1397: point produces itself, so the replica number is unchanged and the returned orientation is simply
1398: `o` composed with `so` via `DMPolytopeTypeComposeOrientation()`.
1400: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetSubcellOrientation()`, `DMPlexTransformCellTransformIdentity()`, `DMPolytopeTypeComposeOrientation()`
1401: @*/
1402: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1403: {
1404: PetscFunctionBegin;
1405: *rnew = r;
1406: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
1407: PetscFunctionReturn(PETSC_SUCCESS);
1408: }
1410: /*@
1411: DMPlexTransformCellTransformIdentity - Default `celltransform` implementation for transforms that reproduce the input mesh
1413: Not Collective
1415: Input Parameters:
1416: + tr - The `DMPlexTransform`
1417: . source - The cell type of the source point
1418: - p - The source point
1420: Output Parameters:
1421: + rt - Refinement type of the source point (set to 0), or `NULL`
1422: . Nt - Number of target cell types produced (always 1)
1423: . target - Array of produced cell types (a single-element array containing `source`)
1424: . size - Array of replica counts for each produced type (a single-element array containing 1)
1425: . cone - Cone description used by `DMPlexTransformGetCone()`; encodes that the replica takes the entire parent cone
1426: - ornt - Orientation array associated with `cone`; all zero for identity
1428: Level: developer
1430: Note:
1431: This routine returns statically allocated arrays describing an identity refinement for each supported
1432: `DMPolytopeType`; every source point produces a single replica of the same type with unchanged cone
1433: and orientation.
1435: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformGetSubcellOrientationIdentity()`
1436: @*/
1437: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1438: {
1439: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
1440: static PetscInt vertexS[] = {1};
1441: static PetscInt vertexC[] = {0};
1442: static PetscInt vertexO[] = {0};
1443: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
1444: static PetscInt edgeS[] = {1};
1445: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1446: static PetscInt edgeO[] = {0, 0};
1447: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
1448: static PetscInt tedgeS[] = {1};
1449: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1450: static PetscInt tedgeO[] = {0, 0};
1451: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
1452: static PetscInt triS[] = {1};
1453: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1454: static PetscInt triO[] = {0, 0, 0};
1455: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
1456: static PetscInt quadS[] = {1};
1457: 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};
1458: static PetscInt quadO[] = {0, 0, 0, 0};
1459: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
1460: static PetscInt tquadS[] = {1};
1461: 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};
1462: static PetscInt tquadO[] = {0, 0, 0, 0};
1463: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
1464: static PetscInt tetS[] = {1};
1465: 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};
1466: static PetscInt tetO[] = {0, 0, 0, 0};
1467: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
1468: static PetscInt hexS[] = {1};
1469: 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};
1470: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
1471: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
1472: static PetscInt tripS[] = {1};
1473: 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};
1474: static PetscInt tripO[] = {0, 0, 0, 0, 0};
1475: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
1476: static PetscInt ttripS[] = {1};
1477: 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};
1478: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
1479: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
1480: static PetscInt tquadpS[] = {1};
1481: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
1482: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1483: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
1484: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
1485: static PetscInt pyrS[] = {1};
1486: 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};
1487: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
1489: PetscFunctionBegin;
1490: if (rt) *rt = 0;
1491: switch (source) {
1492: case DM_POLYTOPE_POINT:
1493: *Nt = 1;
1494: *target = vertexT;
1495: *size = vertexS;
1496: *cone = vertexC;
1497: *ornt = vertexO;
1498: break;
1499: case DM_POLYTOPE_SEGMENT:
1500: *Nt = 1;
1501: *target = edgeT;
1502: *size = edgeS;
1503: *cone = edgeC;
1504: *ornt = edgeO;
1505: break;
1506: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1507: *Nt = 1;
1508: *target = tedgeT;
1509: *size = tedgeS;
1510: *cone = tedgeC;
1511: *ornt = tedgeO;
1512: break;
1513: case DM_POLYTOPE_TRIANGLE:
1514: *Nt = 1;
1515: *target = triT;
1516: *size = triS;
1517: *cone = triC;
1518: *ornt = triO;
1519: break;
1520: case DM_POLYTOPE_QUADRILATERAL:
1521: *Nt = 1;
1522: *target = quadT;
1523: *size = quadS;
1524: *cone = quadC;
1525: *ornt = quadO;
1526: break;
1527: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1528: *Nt = 1;
1529: *target = tquadT;
1530: *size = tquadS;
1531: *cone = tquadC;
1532: *ornt = tquadO;
1533: break;
1534: case DM_POLYTOPE_TETRAHEDRON:
1535: *Nt = 1;
1536: *target = tetT;
1537: *size = tetS;
1538: *cone = tetC;
1539: *ornt = tetO;
1540: break;
1541: case DM_POLYTOPE_HEXAHEDRON:
1542: *Nt = 1;
1543: *target = hexT;
1544: *size = hexS;
1545: *cone = hexC;
1546: *ornt = hexO;
1547: break;
1548: case DM_POLYTOPE_TRI_PRISM:
1549: *Nt = 1;
1550: *target = tripT;
1551: *size = tripS;
1552: *cone = tripC;
1553: *ornt = tripO;
1554: break;
1555: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1556: *Nt = 1;
1557: *target = ttripT;
1558: *size = ttripS;
1559: *cone = ttripC;
1560: *ornt = ttripO;
1561: break;
1562: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1563: *Nt = 1;
1564: *target = tquadpT;
1565: *size = tquadpS;
1566: *cone = tquadpC;
1567: *ornt = tquadpO;
1568: break;
1569: case DM_POLYTOPE_PYRAMID:
1570: *Nt = 1;
1571: *target = pyrT;
1572: *size = pyrS;
1573: *cone = pyrC;
1574: *ornt = pyrO;
1575: break;
1576: default:
1577: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1578: }
1579: PetscFunctionReturn(PETSC_SUCCESS);
1580: }
1582: /*@
1583: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
1585: Not Collective
1587: Input Parameters:
1588: + tr - The `DMPlexTransform`
1589: . sct - The source point cell type, from whom the new cell is being produced
1590: . sp - The source point
1591: . so - The orientation of the source point in its enclosing parent
1592: . tct - The target point cell type
1593: . r - The replica number requested for the produced cell type
1594: - o - The orientation of the replica
1596: Output Parameters:
1597: + rnew - The replica number, given the orientation of the parent
1598: - onew - The replica orientation, given the orientation of the parent
1600: Level: advanced
1602: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1603: @*/
1604: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1605: {
1606: PetscFunctionBeginHot;
1607: PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1608: PetscFunctionReturn(PETSC_SUCCESS);
1609: }
1611: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1612: {
1613: DM dm;
1614: PetscInt pStart, pEnd, pNew;
1616: PetscFunctionBegin;
1617: PetscCall(DMPlexTransformGetDM(tr, &dm));
1618: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0));
1619: /* Must create the celltype label here so that we do not automatically try to compute the types */
1620: PetscCall(DMCreateLabel(rdm, "celltype"));
1621: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1622: for (PetscInt p = pStart; p < pEnd; ++p) {
1623: DMPolytopeType ct;
1624: DMPolytopeType *rct;
1625: PetscInt *rsize, *rcone, *rornt;
1626: PetscInt Nct, n, r;
1628: PetscCall(DMPlexGetCellType(dm, p, &ct));
1629: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1630: for (n = 0; n < Nct; ++n) {
1631: for (r = 0; r < rsize[n]; ++r) {
1632: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1633: PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
1634: PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
1635: }
1636: }
1637: }
1638: /* Let the DM know we have set all the cell types */
1639: {
1640: DMLabel ctLabel;
1641: DM_Plex *plex = (DM_Plex *)rdm->data;
1643: PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
1644: PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
1645: }
1646: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0));
1647: PetscFunctionReturn(PETSC_SUCCESS);
1648: }
1650: /*@
1651: DMPlexTransformGetConeSize - Return the cone size of a point in the transformed mesh
1653: Not Collective
1655: Input Parameters:
1656: + tr - The `DMPlexTransform`
1657: - q - The point number in the transformed mesh
1659: Output Parameter:
1660: . coneSize - The number of points in the cone of `q`
1662: Level: developer
1664: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetCone()`, `DMPlexTransformGetCellType()`, `DMPlexGetConeSize()`
1665: @*/
1666: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1667: {
1668: DMPolytopeType ctNew;
1670: PetscFunctionBegin;
1672: PetscAssertPointer(coneSize, 3);
1673: PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
1674: *coneSize = DMPolytopeTypeGetConeSize(ctNew);
1675: PetscFunctionReturn(PETSC_SUCCESS);
1676: }
1678: /* The orientation o is for the interior of the cell p */
1679: 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[])
1680: {
1681: DM dm;
1682: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1683: const PetscInt *cone;
1684: DMPolytopeType *newft = NULL;
1685: PetscInt c, coff = *coneoff, ooff = *orntoff;
1686: PetscInt dim, cr = 0, co = 0, nr, no;
1688: PetscFunctionBegin;
1689: PetscCall(DMPlexTransformGetDM(tr, &dm));
1690: PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
1691: // Check if we have to permute this cell
1692: PetscCall(DMGetDimension(dm, &dim));
1693: if (DMPolytopeTypeGetDim(ctNew) == dim && DMPolytopeTypeGetDim(ct) == dim - 1) {
1694: PetscCall(DMPlexTransformGetSubcellOrientation(tr, ct, p, o, ctNew, cr, co, &nr, &no));
1695: if (cr != nr || co != no) PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newft));
1696: }
1697: for (c = 0; c < csizeNew; ++c) {
1698: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1699: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1700: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1701: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1702: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1703: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1704: const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1705: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1706: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1707: PetscInt lc;
1709: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1710: for (lc = 0; lc < fn; ++lc) {
1711: const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po);
1712: const PetscInt acp = rcone[coff++];
1713: const PetscInt pcp = parr[acp * 2];
1714: const PetscInt pco = parr[acp * 2 + 1];
1715: const PetscInt *ppornt;
1717: ppp = pp;
1718: pp = pcone[pcp];
1719: PetscCall(DMPlexGetCellType(dm, pp, &pct));
1720: // Restore the parent cone from the last iterate
1721: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
1722: PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
1723: PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
1724: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1725: PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
1726: }
1727: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
1728: pr = rcone[coff++];
1729: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1730: PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
1731: PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
1732: orntNew[c] = fo;
1733: if (newft) newft[c] = ft;
1734: }
1735: PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
1736: if (newft) {
1737: const PetscInt *arr;
1738: PetscInt *newcone, *newornt;
1740: arr = DMPolytopeTypeGetArrangement(ctNew, no);
1741: PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newcone));
1742: PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newornt));
1743: for (PetscInt c = 0; c < csizeNew; ++c) {
1744: DMPolytopeType ft = newft[c];
1745: PetscInt nO;
1747: nO = DMPolytopeTypeGetNumArrangements(ft) / 2;
1748: newcone[c] = coneNew[arr[c * 2 + 0]];
1749: newornt[c] = DMPolytopeTypeComposeOrientation(ft, arr[c * 2 + 1], orntNew[arr[c * 2 + 0]]);
1750: 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]);
1751: }
1752: for (PetscInt c = 0; c < csizeNew; ++c) {
1753: coneNew[c] = newcone[c];
1754: orntNew[c] = newornt[c];
1755: }
1756: PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newcone));
1757: PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newornt));
1758: PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newft));
1759: }
1760: *coneoff = coff;
1761: *orntoff = ooff;
1762: PetscFunctionReturn(PETSC_SUCCESS);
1763: }
1765: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1766: {
1767: DM dm;
1768: DMPolytopeType ct;
1769: PetscInt *coneNew, *orntNew;
1770: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1772: PetscFunctionBegin;
1773: PetscCall(DMPlexTransformGetDM(tr, &dm));
1774: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0));
1775: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1776: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1777: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1778: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1779: for (p = pStart; p < pEnd; ++p) {
1780: PetscInt coff, ooff;
1781: DMPolytopeType *rct;
1782: PetscInt *rsize, *rcone, *rornt;
1783: PetscInt Nct, n, r;
1785: PetscCall(DMPlexGetCellType(dm, p, &ct));
1786: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1787: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1788: const DMPolytopeType ctNew = rct[n];
1790: for (r = 0; r < rsize[n]; ++r) {
1791: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1792: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
1793: PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
1794: PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
1795: }
1796: }
1797: }
1798: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1799: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1800: PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
1801: PetscCall(DMPlexSymmetrize(rdm));
1802: PetscCall(DMPlexStratify(rdm));
1803: PetscCall(DMPlexTransformOrderSupports(tr, dm, rdm));
1804: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0));
1805: PetscFunctionReturn(PETSC_SUCCESS);
1806: }
1808: /*@
1809: DMPlexTransformGetConeOriented - Return the cone of a point in the transformed mesh, computed using a specified parent orientation
1811: Not Collective
1813: Input Parameters:
1814: + tr - The `DMPlexTransform`
1815: . q - The point number in the transformed mesh
1816: - po - The orientation of the parent cell in the original mesh to use when producing the cone
1818: Output Parameters:
1819: + cone - The cone points, obtained from an internal work array
1820: - ornt - The orientations of the cone points, obtained from an internal work array
1822: Level: developer
1824: Note:
1825: Both `cone` and `ornt` are returned in work arrays that must be released with `DMPlexTransformRestoreCone()`.
1827: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetCone()`, `DMPlexTransformRestoreCone()`, `DMPlexTransformGetConeSize()`
1828: @*/
1829: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1830: {
1831: DM dm;
1832: DMPolytopeType ct, qct;
1833: DMPolytopeType *rct;
1834: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1835: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1837: PetscFunctionBegin;
1839: PetscAssertPointer(cone, 4);
1840: PetscAssertPointer(ornt, 5);
1841: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1842: PetscCall(DMPlexTransformGetDM(tr, &dm));
1843: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1844: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1845: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1846: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1847: for (n = 0; n < Nct; ++n) {
1848: const DMPolytopeType ctNew = rct[n];
1849: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1850: PetscInt Nr = rsize[n], fn, c;
1852: if (ctNew == qct) Nr = r;
1853: for (nr = 0; nr < Nr; ++nr) {
1854: for (c = 0; c < csizeNew; ++c) {
1855: ++coff; /* Cell type of new cone point */
1856: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1857: coff += fn;
1858: ++coff; /* Replica number of new cone point */
1859: ++ooff; /* Orientation of new cone point */
1860: }
1861: }
1862: if (ctNew == qct) break;
1863: }
1864: PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1865: *cone = qcone;
1866: *ornt = qornt;
1867: PetscFunctionReturn(PETSC_SUCCESS);
1868: }
1870: /*@
1871: DMPlexTransformGetCone - Return the cone of a point in the transformed mesh
1873: Not Collective
1875: Input Parameters:
1876: + tr - The `DMPlexTransform`
1877: - q - The point number in the transformed mesh
1879: Output Parameters:
1880: + cone - The cone points, obtained from an internal work array, or `NULL` if not requested
1881: - ornt - The orientations of the cone points, obtained from an internal work array, or `NULL` if not requested
1883: Level: developer
1885: Note:
1886: Any non-`NULL` output must be released with `DMPlexTransformRestoreCone()`.
1888: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRestoreCone()`, `DMPlexTransformGetConeOriented()`, `DMPlexTransformGetConeSize()`, `DMPlexGetCone()`
1889: @*/
1890: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1891: {
1892: DM dm;
1893: DMPolytopeType ct, qct;
1894: DMPolytopeType *rct;
1895: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1896: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1898: PetscFunctionBegin;
1900: if (cone) PetscAssertPointer(cone, 3);
1901: if (ornt) PetscAssertPointer(ornt, 4);
1902: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1903: PetscCall(DMPlexTransformGetDM(tr, &dm));
1904: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1905: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1906: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1907: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1908: for (n = 0; n < Nct; ++n) {
1909: const DMPolytopeType ctNew = rct[n];
1910: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1911: PetscInt Nr = rsize[n], fn, c;
1913: if (ctNew == qct) Nr = r;
1914: for (nr = 0; nr < Nr; ++nr) {
1915: for (c = 0; c < csizeNew; ++c) {
1916: ++coff; /* Cell type of new cone point */
1917: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1918: coff += fn;
1919: ++coff; /* Replica number of new cone point */
1920: ++ooff; /* Orientation of new cone point */
1921: }
1922: }
1923: if (ctNew == qct) break;
1924: }
1925: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1926: if (cone) *cone = qcone;
1927: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1928: if (ornt) *ornt = qornt;
1929: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1930: PetscFunctionReturn(PETSC_SUCCESS);
1931: }
1933: /*@
1934: DMPlexTransformRestoreCone - Return the work arrays produced by `DMPlexTransformGetCone()` or `DMPlexTransformGetConeOriented()`
1936: Not Collective
1938: Input Parameters:
1939: + tr - The `DMPlexTransform`
1940: . q - The point number in the transformed mesh
1941: . cone - The cone points to release, or `NULL`
1942: - ornt - The orientations to release, or `NULL`
1944: Level: developer
1946: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetCone()`, `DMPlexTransformGetConeOriented()`
1947: @*/
1948: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1949: {
1950: DM dm;
1952: PetscFunctionBegin;
1954: PetscCall(DMPlexTransformGetDM(tr, &dm));
1955: if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
1956: if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
1957: PetscFunctionReturn(PETSC_SUCCESS);
1958: }
1960: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1961: {
1962: PetscFunctionBegin;
1963: PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
1964: for (PetscInt ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1965: const DMPolytopeType ct = (DMPolytopeType)ict;
1966: DMPlexTransform reftr;
1967: DM refdm, trdm;
1968: Vec coordinates;
1969: const PetscScalar *coords;
1970: DMPolytopeType *rct;
1971: PetscInt *rsize, *rcone, *rornt;
1972: PetscInt Nct, n, r, pNew = 0;
1973: PetscInt trdim, vStart, vEnd, Nc;
1974: const PetscInt debug = 0;
1975: const char *typeName;
1977: /* Since points are 0-dimensional, coordinates make no sense */
1978: if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue;
1979: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
1980: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
1981: PetscCall(DMPlexTransformSetDM(reftr, refdm));
1982: PetscCall(DMPlexTransformGetType(tr, &typeName));
1983: PetscCall(DMPlexTransformSetType(reftr, typeName));
1984: PetscCall(DMPlexTransformSetUp(reftr));
1985: PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));
1987: PetscCall(DMGetDimension(trdm, &trdim));
1988: PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
1989: tr->trNv[ct] = vEnd - vStart;
1990: PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
1991: PetscCall(VecGetLocalSize(coordinates, &Nc));
1992: 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);
1993: PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
1994: PetscCall(VecGetArrayRead(coordinates, &coords));
1995: PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
1996: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1998: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
1999: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2000: for (n = 0; n < Nct; ++n) {
2001: /* Since points are 0-dimensional, coordinates make no sense */
2002: if (rct[n] == DM_POLYTOPE_POINT) continue;
2003: PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
2004: for (r = 0; r < rsize[n]; ++r) {
2005: PetscInt *closure = NULL;
2006: PetscInt clSize, cl, Nv = 0;
2008: PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
2009: PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
2010: PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
2011: for (cl = 0; cl < clSize * 2; cl += 2) {
2012: const PetscInt sv = closure[cl];
2014: if (sv >= vStart && sv < vEnd) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
2015: }
2016: PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
2017: 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]);
2018: }
2019: }
2020: if (debug) {
2021: DMPolytopeType *rct;
2022: PetscInt *rsize, *rcone, *rornt;
2023: PetscInt v, dE = trdim, d, off = 0;
2025: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
2026: for (v = 0; v < tr->trNv[ct]; ++v) {
2027: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
2028: for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
2029: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
2030: }
2032: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2033: for (n = 0; n < Nct; ++n) {
2034: if (rct[n] == DM_POLYTOPE_POINT) continue;
2035: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
2036: for (r = 0; r < rsize[n]; ++r) {
2037: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
2038: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
2039: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
2040: }
2041: }
2042: }
2043: PetscCall(DMDestroy(&refdm));
2044: PetscCall(DMDestroy(&trdm));
2045: PetscCall(DMPlexTransformDestroy(&reftr));
2046: }
2047: PetscFunctionReturn(PETSC_SUCCESS);
2048: }
2050: /*@
2051: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
2053: Input Parameters:
2054: + tr - The `DMPlexTransform` object
2055: - ct - The cell type
2057: Output Parameters:
2058: + Nv - The number of transformed vertices in the closure of the reference cell of given type
2059: - trVerts - The coordinates of these vertices in the reference cell
2061: Level: developer
2063: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
2064: @*/
2065: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
2066: {
2067: PetscFunctionBegin;
2068: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
2069: if (Nv) *Nv = tr->trNv[ct];
2070: if (trVerts) *trVerts = tr->trVerts[ct];
2071: PetscFunctionReturn(PETSC_SUCCESS);
2072: }
2074: /*@
2075: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
2077: Input Parameters:
2078: + tr - The `DMPlexTransform` object
2079: . ct - The cell type
2080: . rct - The subcell type
2081: - r - The subcell index
2083: Output Parameter:
2084: . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()`
2086: Level: developer
2088: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
2089: @*/
2090: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
2091: {
2092: PetscFunctionBegin;
2093: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
2094: PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
2095: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
2096: PetscFunctionReturn(PETSC_SUCCESS);
2097: }
2099: /*@
2100: DMPlexTransformOrderSupports - Reorder newly introduced point supports
2102: Collective
2104: Input Parameters:
2105: + tr - The `DMPlexTransform`
2106: . dm - The original `DM`
2107: - trdm - The transformed `DM` which is reordered
2109: Level: intermediate
2111: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
2112: @*/
2113: PetscErrorCode DMPlexTransformOrderSupports(DMPlexTransform tr, DM dm, DM trdm)
2114: {
2115: PetscFunctionBegin;
2119: PetscTryTypeMethod(tr, ordersupports, dm, trdm);
2120: PetscFunctionReturn(PETSC_SUCCESS);
2121: }
2123: /* Computes new vertex as the barycenter, or centroid */
2124: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
2125: {
2126: PetscInt v, d;
2128: PetscFunctionBeginHot;
2129: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
2130: for (d = 0; d < dE; ++d) out[d] = 0.0;
2131: for (v = 0; v < Nv; ++v)
2132: for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
2133: for (d = 0; d < dE; ++d) out[d] /= Nv;
2134: PetscFunctionReturn(PETSC_SUCCESS);
2135: }
2137: /*@
2138: DMPlexTransformMapCoordinates - Calculate new coordinates for produced points
2140: Not collective
2142: Input Parameters:
2143: + tr - The `DMPlexTransform`
2144: . pct - The cell type of the parent, from whom the new cell is being produced
2145: . ct - The type being produced
2146: . p - The original point
2147: . r - The replica number requested for the produced cell type
2148: . Nv - Number of vertices in the closure of the parent cell
2149: . dE - Spatial dimension
2150: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
2152: Output Parameter:
2153: . out - The coordinates of the new vertices
2155: Level: intermediate
2157: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
2158: @*/
2159: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
2160: {
2161: PetscFunctionBeginHot;
2162: if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
2163: PetscFunctionReturn(PETSC_SUCCESS);
2164: }
2166: /*
2167: DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label
2169: Not Collective
2171: Input Parameters:
2172: + tr - The `DMPlexTransform`
2173: . label - The label in the transformed mesh
2174: . pp - The parent point in the original mesh
2175: . pct - The cell type of the parent point
2176: . p - The point in the transformed mesh
2177: . ct - The cell type of the point
2178: . r - The replica number of the point
2179: - val - The label value of the parent point
2181: Level: developer
2183: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
2184: */
2185: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
2186: {
2187: PetscFunctionBeginHot;
2188: if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
2189: PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
2190: PetscFunctionReturn(PETSC_SUCCESS);
2191: }
2193: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
2194: {
2195: DM dm;
2196: IS valueIS;
2197: const PetscInt *values;
2198: PetscInt defVal, Nv, val;
2200: PetscFunctionBegin;
2201: PetscCall(DMPlexTransformGetDM(tr, &dm));
2202: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2203: PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
2204: PetscCall(DMLabelGetValueIS(label, &valueIS));
2205: PetscCall(ISGetLocalSize(valueIS, &Nv));
2206: PetscCall(ISGetIndices(valueIS, &values));
2207: for (val = 0; val < Nv; ++val) {
2208: IS pointIS;
2209: const PetscInt *points;
2210: PetscInt numPoints, p;
2212: /* Ensure refined label is created with same number of strata as
2213: * original (even if no entries here). */
2214: PetscCall(DMLabelAddStratum(labelNew, values[val]));
2215: PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
2216: PetscCall(ISGetLocalSize(pointIS, &numPoints));
2217: PetscCall(ISGetIndices(pointIS, &points));
2218: for (p = 0; p < numPoints; ++p) {
2219: const PetscInt point = points[p];
2220: DMPolytopeType ct;
2221: DMPolytopeType *rct;
2222: PetscInt *rsize, *rcone, *rornt;
2223: PetscInt Nct, n, r, pNew = 0;
2225: PetscCall(DMPlexGetCellType(dm, point, &ct));
2226: PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2227: for (n = 0; n < Nct; ++n) {
2228: for (r = 0; r < rsize[n]; ++r) {
2229: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
2230: PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
2231: }
2232: }
2233: }
2234: PetscCall(ISRestoreIndices(pointIS, &points));
2235: PetscCall(ISDestroy(&pointIS));
2236: }
2237: PetscCall(ISRestoreIndices(valueIS, &values));
2238: PetscCall(ISDestroy(&valueIS));
2239: PetscFunctionReturn(PETSC_SUCCESS);
2240: }
2242: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
2243: {
2244: DM dm;
2245: PetscInt numLabels, l;
2247: PetscFunctionBegin;
2248: PetscCall(DMPlexTransformGetDM(tr, &dm));
2249: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0));
2250: PetscCall(DMGetNumLabels(dm, &numLabels));
2251: for (l = 0; l < numLabels; ++l) {
2252: DMLabel label, labelNew;
2253: const char *lname;
2254: PetscBool isDepth, isCellType;
2256: PetscCall(DMGetLabelName(dm, l, &lname));
2257: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
2258: if (isDepth) continue;
2259: PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
2260: if (isCellType) continue;
2261: PetscCall(DMCreateLabel(rdm, lname));
2262: PetscCall(DMGetLabel(dm, lname, &label));
2263: PetscCall(DMGetLabel(rdm, lname, &labelNew));
2264: PetscCall(RefineLabel_Internal(tr, label, labelNew));
2265: }
2266: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0));
2267: PetscFunctionReturn(PETSC_SUCCESS);
2268: }
2270: /*@
2271: DMPlexTransformCreateDiscLabels - Refine the labels which define field and discrete system regions on the transformed `DM`
2273: Not Collective
2275: Input Parameters:
2276: + tr - The `DMPlexTransform`
2277: - rdm - The refined `DM` produced by the transform
2279: Level: developer
2281: Note:
2282: Region labels attached to fields (see `DMSetField()`) and to discrete systems (see `DMSetRegionNumDS()`) are
2283: not automatically included in the list of `DM` labels, so this routine walks each field and each `PetscDS`
2284: and updates the labels to refer to the refined points.
2286: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMSetField()`, `DMSetRegionNumDS()`
2287: @*/
2288: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
2289: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
2290: {
2291: DM dm;
2292: PetscInt Nf, f, Nds, s;
2294: PetscFunctionBegin;
2295: PetscCall(DMPlexTransformGetDM(tr, &dm));
2296: PetscCall(DMGetNumFields(dm, &Nf));
2297: for (f = 0; f < Nf; ++f) {
2298: DMLabel label, labelNew;
2299: PetscObject obj;
2300: const char *lname;
2302: PetscCall(DMGetField(rdm, f, &label, &obj));
2303: if (!label) continue;
2304: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2305: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
2306: PetscCall(RefineLabel_Internal(tr, label, labelNew));
2307: PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
2308: PetscCall(DMLabelDestroy(&labelNew));
2309: }
2310: PetscCall(DMGetNumDS(dm, &Nds));
2311: for (s = 0; s < Nds; ++s) {
2312: DMLabel label, labelNew;
2313: const char *lname;
2315: PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
2316: if (!label) continue;
2317: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2318: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
2319: PetscCall(RefineLabel_Internal(tr, label, labelNew));
2320: PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
2321: PetscCall(DMLabelDestroy(&labelNew));
2322: }
2323: PetscFunctionReturn(PETSC_SUCCESS);
2324: }
2326: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
2327: {
2328: DM dm;
2329: PetscSF sf, sfNew;
2330: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
2331: const PetscInt *localPoints;
2332: const PetscSFNode *remotePoints;
2333: PetscInt *localPointsNew;
2334: PetscSFNode *remotePointsNew;
2335: PetscInt pStartNew, pEndNew, pNew;
2336: /* Brute force algorithm */
2337: PetscSF rsf;
2338: PetscSection s;
2339: const PetscInt *rootdegree;
2340: PetscInt *rootPointsNew, *remoteOffsets;
2341: PetscInt numPointsNew, pStart, pEnd, p;
2343: PetscFunctionBegin;
2344: PetscCall(DMPlexTransformGetDM(tr, &dm));
2345: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
2346: PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
2347: PetscCall(DMGetPointSF(dm, &sf));
2348: PetscCall(DMGetPointSF(rdm, &sfNew));
2349: /* Calculate size of new SF */
2350: PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
2351: if (numRoots < 0) {
2352: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
2353: PetscFunctionReturn(PETSC_SUCCESS);
2354: }
2355: for (l = 0; l < numLeaves; ++l) {
2356: const PetscInt p = localPoints ? localPoints[l] : l;
2357: DMPolytopeType ct;
2358: DMPolytopeType *rct;
2359: PetscInt *rsize, *rcone, *rornt;
2360: PetscInt Nct, n;
2362: PetscCall(DMPlexGetCellType(dm, p, &ct));
2363: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2364: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
2365: }
2366: /* Send new root point numbers
2367: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
2368: */
2369: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2370: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
2371: PetscCall(PetscSectionSetChart(s, pStart, pEnd));
2372: for (p = pStart; p < pEnd; ++p) {
2373: DMPolytopeType ct;
2374: DMPolytopeType *rct;
2375: PetscInt *rsize, *rcone, *rornt;
2376: PetscInt Nct, n;
2378: PetscCall(DMPlexGetCellType(dm, p, &ct));
2379: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2380: for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
2381: }
2382: PetscCall(PetscSectionSetUp(s));
2383: PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
2384: PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
2385: PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
2386: PetscCall(PetscFree(remoteOffsets));
2387: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
2388: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
2389: PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
2390: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
2391: for (p = pStart; p < pEnd; ++p) {
2392: DMPolytopeType ct;
2393: DMPolytopeType *rct;
2394: PetscInt *rsize, *rcone, *rornt;
2395: PetscInt Nct, n, r, off;
2397: if (!rootdegree[p - pStart]) continue;
2398: PetscCall(PetscSectionGetOffset(s, p, &off));
2399: PetscCall(DMPlexGetCellType(dm, p, &ct));
2400: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2401: for (n = 0, m = 0; n < Nct; ++n) {
2402: for (r = 0; r < rsize[n]; ++r, ++m) {
2403: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2404: rootPointsNew[off + m] = pNew;
2405: }
2406: }
2407: }
2408: PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
2409: PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
2410: PetscCall(PetscSFDestroy(&rsf));
2411: PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
2412: PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
2413: for (l = 0, m = 0; l < numLeaves; ++l) {
2414: const PetscInt p = localPoints ? localPoints[l] : l;
2415: DMPolytopeType ct;
2416: DMPolytopeType *rct;
2417: PetscInt *rsize, *rcone, *rornt;
2418: PetscInt Nct, n, r, q, off;
2420: PetscCall(PetscSectionGetOffset(s, p, &off));
2421: PetscCall(DMPlexGetCellType(dm, p, &ct));
2422: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2423: for (n = 0, q = 0; n < Nct; ++n) {
2424: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
2425: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2426: localPointsNew[m] = pNew;
2427: remotePointsNew[m].index = rootPointsNew[off + q];
2428: remotePointsNew[m].rank = remotePoints[l].rank;
2429: }
2430: }
2431: }
2432: PetscCall(PetscSectionDestroy(&s));
2433: PetscCall(PetscFree(rootPointsNew));
2434: /* SF needs sorted leaves to correctly calculate Gather */
2435: {
2436: PetscSFNode *rp, *rtmp;
2437: PetscInt *lp, *idx, *ltmp, i;
2439: PetscCall(PetscMalloc1(numLeavesNew, &idx));
2440: PetscCall(PetscMalloc1(numLeavesNew, &lp));
2441: PetscCall(PetscMalloc1(numLeavesNew, &rp));
2442: for (i = 0; i < numLeavesNew; ++i) {
2443: 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);
2444: idx[i] = i;
2445: }
2446: PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
2447: for (i = 0; i < numLeavesNew; ++i) {
2448: lp[i] = localPointsNew[idx[i]];
2449: rp[i] = remotePointsNew[idx[i]];
2450: }
2451: ltmp = localPointsNew;
2452: localPointsNew = lp;
2453: rtmp = remotePointsNew;
2454: remotePointsNew = rp;
2455: PetscCall(PetscFree(idx));
2456: PetscCall(PetscFree(ltmp));
2457: PetscCall(PetscFree(rtmp));
2458: }
2459: PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
2460: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
2461: if (PetscDefined(USE_DEBUG)) {
2462: PetscInt overlap;
2464: // Need to set overlap because some transforms put cells in the overlap
2465: PetscCall(DMPlexGetOverlap(rdm, &overlap));
2466: PetscCall(DMPlexSetOverlap(rdm, NULL, 1));
2467: PetscCall(DMPlexCheckPointSF(rdm, sfNew, PETSC_FALSE));
2468: PetscCall(DMPlexSetOverlap(rdm, NULL, overlap));
2469: }
2470: PetscFunctionReturn(PETSC_SUCCESS);
2471: }
2473: /*
2474: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`.
2476: Not Collective
2478: Input Parameters:
2479: + tr - The `DMPlexTransform`
2480: . ct - The type of the parent cell
2481: . rct - The type of the produced cell
2482: . r - The index of the produced cell
2483: - x - The localized coordinates for the parent cell
2485: Output Parameter:
2486: . xr - The localized coordinates for the produced cell
2488: Level: developer
2490: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
2491: */
2492: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
2493: {
2494: PetscFE fe = NULL;
2495: PetscInt cdim, v, *subcellV;
2497: PetscFunctionBegin;
2498: PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
2499: PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
2500: PetscCall(PetscFEGetNumComponents(fe, &cdim));
2501: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
2502: PetscFunctionReturn(PETSC_SUCCESS);
2503: }
2505: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
2506: {
2507: DM dm, cdm, cdmCell, cdmNew, cdmCellNew;
2508: PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
2509: Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
2510: const PetscScalar *coords;
2511: PetscScalar *coordsNew;
2512: const PetscReal *maxCell, *Lstart, *L;
2513: PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE, sparseLocalize;
2514: PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;
2516: PetscFunctionBegin;
2517: // Need to clear the DMField for coordinates
2518: PetscCall(DMSetCoordinateField(rdm, NULL));
2519: PetscCall(DMPlexTransformGetDM(tr, &dm));
2520: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0));
2521: PetscCall(DMGetCoordinateDM(dm, &cdm));
2522: PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
2523: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
2524: PetscCall(DMGetSparseLocalize(dm, &sparseLocalize));
2525: PetscCall(DMSetSparseLocalize(rdm, sparseLocalize));
2526: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
2527: if (localized) {
2528: /* Localize coordinates of new vertices */
2529: localizeVertices = PETSC_TRUE;
2530: /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
2531: if (!maxCell) localizeCells = PETSC_TRUE;
2532: }
2533: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2534: PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
2535: PetscCall(DMGetCoordinateDim(rdm, &dE));
2536: if (maxCell) {
2537: PetscReal *LstartNew, *LNew, *maxCellNew;
2539: PetscCall(PetscMalloc3(dE, &LstartNew, dE, &LNew, dE, &maxCellNew));
2540: for (d = 0; d < dEo; ++d) {
2541: LstartNew[d] = Lstart[d];
2542: LNew[d] = L[d];
2543: maxCellNew[d] = maxCell[d] / tr->redFactor;
2544: }
2545: for (d = dEo; d < dE; ++d) {
2546: LstartNew[d] = 0.;
2547: LNew[d] = -1.;
2548: maxCellNew[d] = -1.;
2549: }
2550: PetscCall(DMSetPeriodicity(rdm, maxCellNew, LstartNew, LNew));
2551: PetscCall(PetscFree3(LstartNew, LNew, maxCellNew));
2552: }
2553: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
2554: PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
2555: PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
2556: PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
2557: PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
2558: /* Localization should be inherited */
2559: /* Stefano calculates parent cells for each new cell for localization */
2560: /* Localized cells need coordinates of closure */
2561: for (v = vStartNew; v < vEndNew; ++v) {
2562: PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
2563: PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
2564: }
2565: PetscCall(PetscSectionSetUp(coordSectionNew));
2566: PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));
2568: if (localizeCells) {
2569: PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
2570: PetscCall(DMClone(cdmNew, &cdmCellNew));
2571: PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
2572: PetscCall(DMDestroy(&cdmCellNew));
2574: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
2575: PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
2576: PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
2577: PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
2578: PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));
2580: PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
2581: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2582: for (c = cStart; c < cEnd; ++c) {
2583: PetscInt dof;
2585: PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
2586: if (dof) {
2587: DMPolytopeType ct;
2588: DMPolytopeType *rct;
2589: PetscInt *rsize, *rcone, *rornt;
2590: PetscInt dim, cNew, Nct, n, r;
2592: PetscCall(DMPlexGetCellType(dm, c, &ct));
2593: dim = DMPolytopeTypeGetDim(ct);
2594: PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2595: /* This allows for different cell types */
2596: for (n = 0; n < Nct; ++n) {
2597: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2598: for (r = 0; r < rsize[n]; ++r) {
2599: PetscInt *closure = NULL;
2600: PetscInt clSize, cl, Nv = 0;
2602: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
2603: PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2604: for (cl = 0; cl < clSize * 2; cl += 2) {
2605: if (closure[cl] >= vStartNew && closure[cl] < vEndNew) ++Nv;
2606: }
2607: PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2608: PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
2609: PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
2610: }
2611: }
2612: }
2613: }
2614: PetscCall(PetscSectionSetUp(coordSectionCellNew));
2615: PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
2616: }
2617: PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
2618: {
2619: VecType vtype;
2620: PetscInt coordSizeNew, bs;
2621: const char *name;
2623: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
2624: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
2625: PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
2626: PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
2627: PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
2628: PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
2629: PetscCall(VecGetBlockSize(coordsLocal, &bs));
2630: PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
2631: PetscCall(VecGetType(coordsLocal, &vtype));
2632: PetscCall(VecSetType(coordsLocalNew, vtype));
2633: }
2634: PetscCall(VecGetArrayRead(coordsLocal, &coords));
2635: PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
2636: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2637: /* First set coordinates for vertices */
2638: for (p = pStart; p < pEnd; ++p) {
2639: DMPolytopeType ct;
2640: DMPolytopeType *rct;
2641: PetscInt *rsize, *rcone, *rornt;
2642: PetscInt Nct, n, r;
2643: PetscBool hasVertex = PETSC_FALSE;
2645: PetscCall(DMPlexGetCellType(dm, p, &ct));
2646: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2647: for (n = 0; n < Nct; ++n) {
2648: if (rct[n] == DM_POLYTOPE_POINT) {
2649: hasVertex = PETSC_TRUE;
2650: break;
2651: }
2652: }
2653: if (hasVertex) {
2654: const PetscScalar *icoords = NULL;
2655: const PetscScalar *array = NULL;
2656: PetscScalar *pcoords = NULL;
2657: PetscBool isDG;
2658: PetscInt Nc, Nv, v, d;
2660: PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2662: icoords = pcoords;
2663: Nv = Nc / dEo;
2664: if (ct != DM_POLYTOPE_POINT) {
2665: if (localizeVertices && maxCell) {
2666: PetscScalar anchor[3];
2668: for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
2669: for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
2670: }
2671: }
2672: for (n = 0; n < Nct; ++n) {
2673: if (rct[n] != DM_POLYTOPE_POINT) continue;
2674: for (r = 0; r < rsize[n]; ++r) {
2675: PetscScalar vcoords[3];
2676: PetscInt vNew, off;
2678: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
2679: PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
2680: PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
2681: PetscCall(DMSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
2682: }
2683: }
2684: PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2685: }
2686: }
2687: PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
2688: PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
2689: PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
2690: PetscCall(VecDestroy(&coordsLocalNew));
2691: PetscCall(PetscSectionDestroy(&coordSectionNew));
2692: /* Then set coordinates for cells by localizing */
2693: if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
2694: else {
2695: VecType vtype;
2696: PetscInt coordSizeNew, bs;
2697: const char *name;
2699: PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
2700: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
2701: PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
2702: PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
2703: PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
2704: PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
2705: PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
2706: PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
2707: PetscCall(VecGetType(coordsLocalCell, &vtype));
2708: PetscCall(VecSetType(coordsLocalCellNew, vtype));
2709: PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
2710: PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));
2712: for (p = pStart; p < pEnd; ++p) {
2713: DMPolytopeType ct;
2714: DMPolytopeType *rct;
2715: PetscInt *rsize, *rcone, *rornt;
2716: PetscInt dof = 0, Nct, n, r;
2718: PetscCall(DMPlexGetCellType(dm, p, &ct));
2719: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2720: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
2721: if (dof) {
2722: const PetscScalar *pcoords;
2724: PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords));
2725: for (n = 0; n < Nct; ++n) {
2726: const PetscInt Nr = rsize[n];
2728: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2729: for (r = 0; r < Nr; ++r) {
2730: PetscInt pNew, offNew;
2732: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2733: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2734: cell to the ones it produces. */
2735: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2736: PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
2737: PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
2738: }
2739: }
2740: }
2741: }
2742: PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
2743: PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
2744: PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
2745: PetscCall(VecDestroy(&coordsLocalCellNew));
2746: PetscCall(PetscSectionDestroy(&coordSectionCellNew));
2747: }
2748: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0));
2749: PetscFunctionReturn(PETSC_SUCCESS);
2750: }
2752: /*@
2753: DMPlexTransformApply - Execute the transformation, producing another `DM`
2755: Collective
2757: Input Parameters:
2758: + tr - The `DMPlexTransform` object
2759: - dm - The original `DM`
2761: Output Parameter:
2762: . trdm - The transformed `DM`
2764: Level: intermediate
2766: Options Database Keys:
2767: + -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point
2768: . -dm_plex_transform_label_replica_inc num - Increment for the label value to be multiplied by the replica number
2769: - -dm_plex_transform_active name - Name for active mesh label
2771: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()`
2772: @*/
2773: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *trdm)
2774: {
2775: DM rdm;
2776: DMPlexInterpolatedFlag interp;
2777: PetscInt pStart, pEnd;
2779: PetscFunctionBegin;
2782: PetscAssertPointer(trdm, 3);
2783: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0));
2784: PetscCall(DMPlexTransformSetDM(tr, dm));
2786: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
2787: PetscCall(DMSetType(rdm, DMPLEX));
2788: PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
2789: /* Calculate number of new points of each depth */
2790: PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
2791: PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
2792: /* Step 1: Set chart */
2793: PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
2794: PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
2795: /* Step 2: Set cone/support sizes (automatically stratifies) */
2796: PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
2797: /* Step 3: Setup refined DM */
2798: PetscCall(DMSetUp(rdm));
2799: /* Step 4: Set cones and supports (automatically symmetrizes) */
2800: PetscCall(DMPlexTransformSetCones(tr, rdm));
2801: /* Step 5: Create pointSF */
2802: PetscCall(DMPlexTransformCreateSF(tr, rdm));
2803: /* Step 6: Create labels */
2804: PetscCall(DMPlexTransformCreateLabels(tr, rdm));
2805: /* Step 7: Set coordinates */
2806: PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
2807: // Do not copy periodicity, which was handled in DMPlexTransformSetCoordinates()
2808: PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_TRUE, rdm));
2809: // If the original DM was configured from options, the transformed DM should be as well
2810: rdm->setfromoptionscalled = dm->setfromoptionscalled;
2811: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0));
2812: *trdm = rdm;
2813: PetscFunctionReturn(PETSC_SUCCESS);
2814: }
2816: /*@
2817: DMPlexTransformAdaptLabel - Adapt a `DMPLEX` using a `DMPlexTransform` driven by a `DMLabel` marking cells to be refined or coarsened.
2819: Collective
2821: Input Parameters:
2822: + dm - the input `DMPLEX`
2823: . metric - unused; present to conform to the `DMAdaptor` label-based interface
2824: . adaptLabel - a `DMLabel` marking cells with `DM_ADAPT_REFINE`, `DM_ADAPT_COARSEN`, etc.
2825: - rgLabel - unused region-tag label; present to conform to the `DMAdaptor` interface
2827: Output Parameter:
2828: . rdm - the adapted `DMPLEX`
2830: Level: developer
2832: Note:
2833: This routine is registered as the "cellrefiner" adaptor by `DMGenerateRegisterAll()` and is invoked through `DMAdaptLabel()`.
2835: .seealso: `DMPLEX`, `DMPlexTransform`, `DMAdaptLabel()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`, `DMLabel`
2836: @*/
2837: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2838: {
2839: DMPlexTransform tr;
2840: DM cdm, rcdm;
2841: const char *prefix;
2842: PetscBool save;
2844: PetscFunctionBegin;
2845: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
2846: PetscCall(PetscObjectSetName((PetscObject)tr, "Adapt Label Transform"));
2847: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2848: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
2849: PetscCall(DMPlexTransformSetDM(tr, dm));
2850: PetscCall(DMPlexTransformSetFromOptions(tr));
2851: if (adaptLabel) PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
2852: PetscCall(DMPlexTransformSetUp(tr));
2853: PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
2854: PetscCall(DMPlexTransformApply(tr, dm, rdm));
2855: PetscCall(DMCopyDisc(dm, *rdm));
2856: PetscCall(DMGetCoordinateDM(dm, &cdm));
2857: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
2858: PetscCall(DMCopyDisc(cdm, rcdm));
2859: PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
2860: PetscCall(DMCopyDisc(dm, *rdm));
2861: PetscCall(DMPlexGetSaveTransform(dm, &save));
2862: if (save) PetscCall(DMPlexSetTransform(*rdm, tr));
2863: PetscCall(DMPlexTransformDestroy(&tr));
2864: ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2865: PetscFunctionReturn(PETSC_SUCCESS);
2866: }