Actual source code: plextransform.c
1: #include <petsc/private/dmplextransformimpl.h>
3: #include <petsc/private/petscfeimpl.h>
5: PetscClassId DMPLEXTRANSFORM_CLASSID;
7: PetscFunctionList DMPlexTransformList = NULL;
8: PetscBool DMPlexTransformRegisterAllCalled = PETSC_FALSE;
10: PetscLogEvent DMPLEXTRANSFORM_SetUp, DMPLEXTRANSFORM_Apply, DMPLEXTRANSFORM_SetConeSizes, DMPLEXTRANSFORM_SetCones, DMPLEXTRANSFORM_CreateSF, DMPLEXTRANSFORM_CreateLabels, DMPLEXTRANSFORM_SetCoordinates;
12: /* Construct cell type order since we must loop over cell types in the same dimensional order they are stored in the plex if dm != NULL
13: OR in standard plex ordering if dm == NULL */
14: static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DM dm, PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
15: {
16: PetscInt *ctO, *ctOInv;
17: PetscInt d, c, off = 0;
18: PetscInt dimOrder[5] = {3, 2, 1, 0, -1};
20: PetscFunctionBegin;
21: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv));
22: if (dm) { // Order the dimensions by their starting location
23: PetscInt hStart[4] = {-1, -1, -1, -1};
24: for (d = 0; d <= dim; ++d) PetscCall(DMPlexGetDepthStratum(dm, dim - d, &hStart[d], NULL));
25: PetscCall(PetscSortIntWithArray(dim + 1, hStart, &dimOrder[3 - dim]));
26: } else if (dim > 1) { // Standard plex ordering. dimOrder is in correct order if dim > 1
27: off = 4 - dim;
28: dimOrder[off++] = 0;
29: for (d = dim - 1; d > 0; --d) dimOrder[off++] = d;
30: }
32: off = 0;
33: for (d = 0; d < 5; ++d) {
34: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
35: if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue;
36: if (DMPolytopeTypeGetDim((DMPolytopeType)c) == dimOrder[d]) ctO[off++] = c;
37: }
38: }
39: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
40: if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) ctO[off++] = c;
41: }
42: ctO[off++] = DM_NUM_POLYTOPES;
43: PetscCheck(off == DM_NUM_POLYTOPES + 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %" PetscInt_FMT " for cell type order", off);
45: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c;
47: *ctOrder = ctO;
48: *ctOrderInv = ctOInv;
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: /*@C
53: DMPlexTransformRegister - Adds a new transform component implementation
55: Not Collective
57: Input Parameters:
58: + name - The name of a new user-defined creation routine
59: - create_func - The creation routine
61: Example Usage:
62: .vb
63: DMPlexTransformRegister("my_transform", MyTransformCreate);
64: .ve
66: Then, your transform type can be chosen with the procedural interface via
67: .vb
68: DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
69: DMPlexTransformSetType(DMPlexTransform, "my_transform");
70: .ve
71: or at runtime via the option
72: .vb
73: -dm_plex_transform_type my_transform
74: .ve
76: Level: advanced
78: Note:
79: `DMPlexTransformRegister()` may be called multiple times to add several user-defined transforms
81: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()`
82: @*/
83: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
84: {
85: PetscFunctionBegin;
86: PetscCall(DMInitializePackage());
87: PetscCall(PetscFunctionListAdd(&DMPlexTransformList, name, create_func));
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
92: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
93: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
94: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
95: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
96: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
97: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform);
98: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform);
99: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform);
101: /*@C
102: DMPlexTransformRegisterAll - Registers all of the transform components in the `DM` package.
104: Not Collective
106: Level: advanced
108: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
109: @*/
110: PetscErrorCode DMPlexTransformRegisterAll(void)
111: {
112: PetscFunctionBegin;
113: if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
114: DMPlexTransformRegisterAllCalled = PETSC_TRUE;
116: PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter));
117: PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular));
118: PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox));
119: PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld));
120: PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL));
121: PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR));
122: PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D));
123: PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDETYPE, DMPlexTransformCreate_Extrude));
124: PetscCall(DMPlexTransformRegister(DMPLEXCOHESIVEEXTRUDE, DMPlexTransformCreate_Cohesive));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /*@C
129: DMPlexTransformRegisterDestroy - This function destroys the registered `DMPlexTransformType`. It is called from `PetscFinalize()`.
131: Not collective
133: Level: developer
135: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()`
136: @*/
137: PetscErrorCode DMPlexTransformRegisterDestroy(void)
138: {
139: PetscFunctionBegin;
140: PetscCall(PetscFunctionListDestroy(&DMPlexTransformList));
141: DMPlexTransformRegisterAllCalled = PETSC_FALSE;
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: /*@
146: DMPlexTransformCreate - Creates an empty transform object. The type can then be set with `DMPlexTransformSetType()`.
148: Collective
150: Input Parameter:
151: . comm - The communicator for the transform object
153: Output Parameter:
154: . tr - The transform object
156: Level: beginner
158: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
159: @*/
160: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
161: {
162: DMPlexTransform t;
164: PetscFunctionBegin;
165: PetscAssertPointer(tr, 2);
166: *tr = NULL;
167: PetscCall(DMInitializePackage());
169: PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView));
170: t->setupcalled = PETSC_FALSE;
171: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom));
172: *tr = t;
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: /*@
177: DMPlexTransformSetType - Sets the particular implementation for a transform.
179: Collective
181: Input Parameters:
182: + tr - The transform
183: - method - The name of the transform type
185: Options Database Key:
186: . -dm_plex_transform_type <type> - Sets the transform type; see `DMPlexTransformType`
188: Level: intermediate
190: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
191: @*/
192: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
193: {
194: PetscErrorCode (*r)(DMPlexTransform);
195: PetscBool match;
197: PetscFunctionBegin;
199: PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match));
200: if (match) PetscFunctionReturn(PETSC_SUCCESS);
202: PetscCall(DMPlexTransformRegisterAll());
203: PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r));
204: PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);
206: PetscTryTypeMethod(tr, destroy);
207: PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops)));
208: PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method));
209: PetscCall((*r)(tr));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*@
214: DMPlexTransformGetType - Gets the type name (as a string) from the transform.
216: Not Collective
218: Input Parameter:
219: . tr - The `DMPlexTransform`
221: Output Parameter:
222: . type - The `DMPlexTransformType` name
224: Level: intermediate
226: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
227: @*/
228: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
229: {
230: PetscFunctionBegin;
232: PetscAssertPointer(type, 2);
233: PetscCall(DMPlexTransformRegisterAll());
234: *type = ((PetscObject)tr)->type_name;
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
239: {
240: PetscViewerFormat format;
242: PetscFunctionBegin;
243: PetscCall(PetscViewerGetFormat(v, &format));
244: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
245: const PetscInt *trTypes = NULL;
246: IS trIS;
247: PetscInt cols = 8;
248: PetscInt Nrt = 8, f, g;
250: if (tr->trType) PetscCall(DMLabelView(tr->trType, v));
251: PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n"));
252: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
253: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
254: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]));
255: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
256: PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n"));
257: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
258: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
259: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]));
260: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
262: if (tr->trType) {
263: PetscCall(DMLabelGetNumValues(tr->trType, &Nrt));
264: PetscCall(DMLabelGetValueIS(tr->trType, &trIS));
265: PetscCall(ISGetIndices(trIS, &trTypes));
266: }
267: PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n"));
268: PetscCall(PetscViewerASCIIPrintf(v, " "));
269: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
270: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
271: for (f = 0; f < Nrt; ++f) {
272: PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT " |", trTypes ? trTypes[f] : f));
273: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]));
274: PetscCall(PetscViewerASCIIPrintf(v, " |\n"));
275: }
276: if (tr->trType) {
277: PetscCall(ISRestoreIndices(trIS, &trTypes));
278: PetscCall(ISDestroy(&trIS));
279: }
280: }
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*@
285: DMPlexTransformView - Views a `DMPlexTransform`
287: Collective
289: Input Parameters:
290: + tr - the `DMPlexTransform` object to view
291: - v - the viewer
293: Level: beginner
295: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
296: @*/
297: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
298: {
299: PetscBool isascii;
301: PetscFunctionBegin;
303: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v));
305: PetscCheckSameComm(tr, 1, v, 2);
306: PetscCall(PetscViewerCheckWritable(v));
307: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v));
308: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
309: if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v));
310: PetscTryTypeMethod(tr, view, v);
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*@
315: DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database
317: Collective
319: Input Parameter:
320: . tr - the `DMPlexTransform` object to set options for
322: Options Database Keys:
323: + -dm_plex_transform_type - Set the transform type, e.g. refine_regular
324: . -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point
325: . -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
326: . -dm_plex_transform_active <name> - Name for active mesh label
327: - -dm_plex_transform_active_values <v0,v1,...> - Values in the active label
329: Level: intermediate
331: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
332: @*/
333: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
334: {
335: char typeName[1024], active[PETSC_MAX_PATH_LEN];
336: const char *defName = DMPLEXREFINEREGULAR;
337: PetscBool flg, match;
339: PetscFunctionBegin;
341: PetscObjectOptionsBegin((PetscObject)tr);
342: PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg));
343: if (flg) PetscCall(DMPlexTransformSetType(tr, typeName));
344: else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName));
345: PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &match, &flg));
346: if (flg) PetscCall(DMPlexTransformSetMatchStrata(tr, match));
347: 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));
348: PetscCall(PetscOptionsString("-dm_plex_transform_active", "Name for active mesh label", "DMPlexTransformSetActive", active, active, sizeof(active), &flg));
349: if (flg) {
350: DM dm;
351: DMLabel label;
352: PetscInt values[16];
353: PetscInt n = 16;
355: PetscCall(DMPlexTransformGetDM(tr, &dm));
356: PetscCall(DMGetLabel(dm, active, &label));
357: PetscCall(PetscOptionsIntArray("-dm_plex_transform_active_values", "The label values to be active", "DMPlexTransformSetActive", values, &n, &flg));
358: if (flg && n) {
359: DMLabel newlabel;
361: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Active", &newlabel));
362: for (PetscInt i = 0; i < n; ++i) {
363: IS is;
365: PetscCall(DMLabelGetStratumIS(label, values[i], &is));
366: PetscCall(DMLabelInsertIS(newlabel, is, values[i]));
367: PetscCall(ISDestroy(&is));
368: }
369: PetscCall(DMPlexTransformSetActive(tr, newlabel));
370: PetscCall(DMLabelDestroy(&newlabel));
371: } else {
372: PetscCall(DMPlexTransformSetActive(tr, label));
373: }
374: }
375: PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
376: /* process any options handlers added with PetscObjectAddOptionsHandler() */
377: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject));
378: PetscOptionsEnd();
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: /*@
383: DMPlexTransformDestroy - Destroys a `DMPlexTransform`
385: Collective
387: Input Parameter:
388: . tr - the transform object to destroy
390: Level: beginner
392: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
393: @*/
394: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
395: {
396: PetscInt c;
398: PetscFunctionBegin;
399: if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
401: if (--((PetscObject)*tr)->refct > 0) {
402: *tr = NULL;
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: PetscTryTypeMethod(*tr, destroy);
407: PetscCall(DMDestroy(&(*tr)->dm));
408: PetscCall(DMLabelDestroy(&(*tr)->active));
409: PetscCall(DMLabelDestroy(&(*tr)->trType));
410: PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld));
411: PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew));
412: PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew));
413: PetscCall(PetscFree((*tr)->offset));
414: PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd));
415: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
416: PetscCall(PetscFEDestroy(&(*tr)->coordFE[c]));
417: PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c]));
418: }
419: if ((*tr)->trVerts) {
420: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
421: DMPolytopeType *rct;
422: PetscInt *rsize, *rcone, *rornt, Nct, n, r;
424: if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) {
425: PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
426: for (n = 0; n < Nct; ++n) {
427: if (rct[n] == DM_POLYTOPE_POINT) continue;
428: for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r]));
429: PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]]));
430: }
431: }
432: PetscCall(PetscFree((*tr)->trSubVerts[c]));
433: PetscCall(PetscFree((*tr)->trVerts[c]));
434: }
435: }
436: PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts));
437: PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom));
438: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
439: PetscCall(PetscHeaderDestroy(tr));
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
444: {
445: DMLabel trType = tr->trType;
446: PetscInt c, cN, *off;
448: PetscFunctionBegin;
449: if (trType) {
450: DM dm;
451: IS rtIS;
452: const PetscInt *reftypes;
453: PetscInt Nrt, r;
455: PetscCall(DMPlexTransformGetDM(tr, &dm));
456: PetscCall(DMLabelGetNumValues(trType, &Nrt));
457: PetscCall(DMLabelGetValueIS(trType, &rtIS));
458: PetscCall(ISGetIndices(rtIS, &reftypes));
459: PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off));
460: for (r = 0; r < Nrt; ++r) {
461: const PetscInt rt = reftypes[r];
462: IS rtIS;
463: const PetscInt *points;
464: DMPolytopeType ct;
465: PetscInt np, p;
467: PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
468: PetscCall(ISGetLocalSize(rtIS, &np));
469: PetscCall(ISGetIndices(rtIS, &points));
470: if (!np) continue;
471: p = points[0];
472: PetscCall(ISRestoreIndices(rtIS, &points));
473: PetscCall(ISDestroy(&rtIS));
474: PetscCall(DMPlexGetCellType(dm, p, &ct));
475: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
476: const DMPolytopeType ctNew = (DMPolytopeType)cN;
477: DMPolytopeType *rct;
478: PetscInt *rsize, *cone, *ornt;
479: PetscInt Nct, n, s;
481: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
482: off[r * DM_NUM_POLYTOPES + ctNew] = -1;
483: break;
484: }
485: off[r * DM_NUM_POLYTOPES + ctNew] = 0;
486: for (s = 0; s <= r; ++s) {
487: const PetscInt st = reftypes[s];
488: DMPolytopeType sct;
489: PetscInt q, qrt;
491: PetscCall(DMLabelGetStratumIS(trType, st, &rtIS));
492: PetscCall(ISGetLocalSize(rtIS, &np));
493: PetscCall(ISGetIndices(rtIS, &points));
494: if (!np) continue;
495: q = points[0];
496: PetscCall(ISRestoreIndices(rtIS, &points));
497: PetscCall(ISDestroy(&rtIS));
498: PetscCall(DMPlexGetCellType(dm, q, &sct));
499: PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt));
500: 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);
501: if (st == rt) {
502: for (n = 0; n < Nct; ++n)
503: if (rct[n] == ctNew) break;
504: if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
505: break;
506: }
507: for (n = 0; n < Nct; ++n) {
508: if (rct[n] == ctNew) {
509: PetscInt sn;
511: PetscCall(DMLabelGetStratumSize(trType, st, &sn));
512: off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
513: }
514: }
515: }
516: }
517: }
518: PetscCall(ISRestoreIndices(rtIS, &reftypes));
519: PetscCall(ISDestroy(&rtIS));
520: } else {
521: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off));
522: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
523: const DMPolytopeType ct = (DMPolytopeType)c;
524: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
525: const DMPolytopeType ctNew = (DMPolytopeType)cN;
526: DMPolytopeType *rct;
527: PetscInt *rsize, *cone, *ornt;
528: PetscInt Nct, n, i;
530: 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) {
531: off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
532: continue;
533: }
534: off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
535: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
536: const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i];
537: const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];
539: PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt));
540: if (ict == ct) {
541: for (n = 0; n < Nct; ++n)
542: if (rct[n] == ctNew) break;
543: if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
544: break;
545: }
546: for (n = 0; n < Nct; ++n)
547: if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
548: }
549: }
550: }
551: }
552: *offset = off;
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: /*@
557: DMPlexTransformSetUp - Create the tables that drive the transform
559: Input Parameter:
560: . tr - The `DMPlexTransform` object
562: Level: intermediate
564: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
565: @*/
566: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
567: {
568: DMPolytopeType ctCell;
569: DM dm;
570: PetscInt pStart, pEnd, p, c, celldim = 0;
572: PetscFunctionBegin;
574: if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
575: PetscCall(DMPlexTransformGetDM(tr, &dm));
576: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0));
577: PetscTryTypeMethod(tr, setup);
578: PetscCall(DMSetSnapToGeomModel(dm, NULL));
579: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
581: if (pEnd > pStart) {
582: // Ignore cells hanging off of embedded surfaces
583: PetscInt c = pStart;
585: ctCell = DM_POLYTOPE_FV_GHOST;
586: while (DMPolytopeTypeGetDim(ctCell) < 0) PetscCall(DMPlexGetCellType(dm, c++, &ctCell));
587: } else {
588: PetscInt dim;
590: PetscCall(DMGetDimension(dm, &dim));
591: switch (dim) {
592: case 0:
593: ctCell = DM_POLYTOPE_POINT;
594: break;
595: case 1:
596: ctCell = DM_POLYTOPE_SEGMENT;
597: break;
598: case 2:
599: ctCell = DM_POLYTOPE_TRIANGLE;
600: break;
601: case 3:
602: ctCell = DM_POLYTOPE_TETRAHEDRON;
603: break;
604: default:
605: ctCell = DM_POLYTOPE_UNKNOWN;
606: }
607: }
608: PetscCall(DMPlexCreateCellTypeOrder_Internal(dm, DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
609: for (p = pStart; p < pEnd; ++p) {
610: DMPolytopeType ct;
611: DMPolytopeType *rct;
612: PetscInt *rsize, *cone, *ornt;
613: PetscInt Nct, n;
615: PetscCall(DMPlexGetCellType(dm, p, &ct));
616: 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);
617: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
618: for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
619: }
620: PetscCall(DMPlexCreateCellTypeOrder_Internal(NULL, celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
621: /* Construct sizes and offsets for each cell type */
622: if (!tr->ctStart) {
623: PetscInt *ctS, *ctSN, *ctC, *ctCN;
625: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
626: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
627: for (p = pStart; p < pEnd; ++p) {
628: DMPolytopeType ct;
629: DMPolytopeType *rct;
630: PetscInt *rsize, *cone, *ornt;
631: PetscInt Nct, n;
633: PetscCall(DMPlexGetCellType(dm, p, &ct));
634: 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);
635: ++ctC[ct];
636: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
637: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
638: }
639: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
640: const PetscInt cto = tr->ctOrderOld[c];
641: const PetscInt cton = tr->ctOrderOld[c + 1];
642: const PetscInt ctn = tr->ctOrderNew[c];
643: const PetscInt ctnn = tr->ctOrderNew[c + 1];
645: ctS[cton] = ctS[cto] + ctC[cto];
646: ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
647: }
648: PetscCall(PetscFree2(ctC, ctCN));
649: tr->ctStart = ctS;
650: tr->ctStartNew = ctSN;
651: }
652: PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
653: // Compute depth information
654: tr->depth = -1;
655: for (c = 0; c < DM_NUM_POLYTOPES; ++c)
656: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
657: PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
658: for (PetscInt d = 0; d <= tr->depth; ++d) {
659: tr->depthStart[d] = PETSC_INT_MAX;
660: tr->depthEnd[d] = -1;
661: }
662: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
663: const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);
665: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
666: tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
667: tr->depthEnd[dep] = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
668: }
669: tr->setupcalled = PETSC_TRUE;
670: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0));
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: /*@
675: DMPlexTransformGetDM - Get the base `DM` for the transform
677: Input Parameter:
678: . tr - The `DMPlexTransform` object
680: Output Parameter:
681: . dm - The original `DM` which will be transformed
683: Level: intermediate
685: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
686: @*/
687: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
688: {
689: PetscFunctionBegin;
691: PetscAssertPointer(dm, 2);
692: *dm = tr->dm;
693: PetscFunctionReturn(PETSC_SUCCESS);
694: }
696: /*@
697: DMPlexTransformSetDM - Set the base `DM` for the transform
699: Input Parameters:
700: + tr - The `DMPlexTransform` object
701: - dm - The original `DM` which will be transformed
703: Level: intermediate
705: Note:
706: The user does not typically call this, as it is called by `DMPlexTransformApply()`.
708: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
709: @*/
710: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
711: {
712: PetscFunctionBegin;
715: PetscCall(PetscObjectReference((PetscObject)dm));
716: PetscCall(DMDestroy(&tr->dm));
717: tr->dm = dm;
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: /*@
722: DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform
724: Input Parameter:
725: . tr - The `DMPlexTransform` object
727: Output Parameter:
728: . active - The `DMLabel` indicating which points will be transformed
730: Level: intermediate
732: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
733: @*/
734: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
735: {
736: PetscFunctionBegin;
738: PetscAssertPointer(active, 2);
739: *active = tr->active;
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@
744: DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform
746: Input Parameters:
747: + tr - The `DMPlexTransform` object
748: - active - The `DMLabel` indicating which points will be transformed
750: Level: intermediate
752: Note:
753: This only applies to transforms listed in [](plex_transform_table) that operate on a subset of the mesh.
755: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
756: @*/
757: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
758: {
759: PetscFunctionBegin;
762: PetscCall(PetscObjectReference((PetscObject)active));
763: PetscCall(DMLabelDestroy(&tr->active));
764: tr->active = active;
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: /*@
769: DMPlexTransformGetTransformTypes - Get the `DMLabel` marking the transform type of each point for the transform
771: Input Parameter:
772: . tr - The `DMPlexTransform` object
774: Output Parameter:
775: . trType - The `DMLabel` indicating the transform type for each point
777: Level: intermediate
779: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetTransformType()`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
780: @*/
781: PetscErrorCode DMPlexTransformGetTransformTypes(DMPlexTransform tr, DMLabel *trType)
782: {
783: PetscFunctionBegin;
785: PetscAssertPointer(trType, 2);
786: *trType = tr->trType;
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: /*@
791: DMPlexTransformSetTransformTypes - Set the `DMLabel` marking the transform type of each point for the transform
793: Input Parameters:
794: + tr - The `DMPlexTransform` object
795: - trType - The original `DM` which will be transformed
797: Level: intermediate
799: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetTransformTypes()`, `DMPlexTransformGetActive())`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
800: @*/
801: PetscErrorCode DMPlexTransformSetTransformTypes(DMPlexTransform tr, DMLabel trType)
802: {
803: PetscFunctionBegin;
806: PetscCall(PetscObjectReference((PetscObject)trType));
807: PetscCall(DMLabelDestroy(&tr->trType));
808: tr->trType = trType;
809: PetscFunctionReturn(PETSC_SUCCESS);
810: }
812: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
813: {
814: PetscFunctionBegin;
815: if (!tr->coordFE[ct]) {
816: PetscInt dim, cdim;
818: dim = DMPolytopeTypeGetDim(ct);
819: PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
820: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
821: {
822: PetscDualSpace dsp;
823: PetscQuadrature quad;
824: DM K;
825: PetscFEGeom *cg;
826: PetscScalar *Xq;
827: PetscReal *xq, *wq;
828: PetscInt Nq, q;
830: PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
831: PetscCall(PetscMalloc1(Nq * cdim, &xq));
832: for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
833: PetscCall(PetscMalloc1(Nq, &wq));
834: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
835: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
836: PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
837: PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));
839: PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
840: PetscCall(PetscDualSpaceGetDM(dsp, &K));
841: PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FEGEOM_BASIC, &tr->refGeom[ct]));
842: cg = tr->refGeom[ct];
843: PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
844: PetscCall(PetscQuadratureDestroy(&quad));
845: }
846: }
847: *fe = tr->coordFE[ct];
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
852: {
853: PetscInt dim, cdim;
855: PetscFunctionBegin;
856: PetscCall(DMGetDimension(dm, &dim));
857: PetscCall(DMSetDimension(tdm, dim));
858: PetscCall(DMGetCoordinateDim(dm, &cdim));
859: PetscCall(DMSetCoordinateDim(tdm, cdim));
860: PetscFunctionReturn(PETSC_SUCCESS);
861: }
863: /*@
864: DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`
866: Input Parameters:
867: + tr - The `DMPlexTransform` object
868: - dm - The original `DM`
870: Output Parameter:
871: . tdm - The transformed `DM`
873: Level: advanced
875: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
876: @*/
877: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
878: {
879: PetscFunctionBegin;
880: PetscUseTypeMethod(tr, setdimensions, dm, tdm);
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
885: {
886: PetscFunctionBegin;
887: if (pStart) *pStart = 0;
888: if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
893: {
894: PetscInt ctNew;
896: PetscFunctionBegin;
898: PetscAssertPointer(celltype, 3);
899: /* TODO Can do bisection since everything is sorted */
900: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
901: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
903: if (cell >= ctSN && cell < ctEN) break;
904: }
905: PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
906: *celltype = (DMPolytopeType)ctNew;
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
911: {
912: PetscFunctionBegin;
914: if (start) *start = tr->ctStartNew[celltype];
915: if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
920: {
921: PetscFunctionBegin;
923: *depth = tr->depth;
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
928: {
929: PetscFunctionBegin;
931: if (start) *start = tr->depthStart[depth];
932: if (end) *end = tr->depthEnd[depth];
933: PetscFunctionReturn(PETSC_SUCCESS);
934: }
936: /*@
937: DMPlexTransformGetMatchStrata - Get the flag which determines what points get added to the transformed labels
939: Not Collective
941: Input Parameter:
942: . tr - The `DMPlexTransform`
944: Output Parameter:
945: . match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels
947: Level: intermediate
949: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetMatchStrata()`, `DMPlexGetPointDepth()`
950: @*/
951: PetscErrorCode DMPlexTransformGetMatchStrata(DMPlexTransform tr, PetscBool *match)
952: {
953: PetscFunctionBegin;
955: PetscAssertPointer(match, 2);
956: *match = tr->labelMatchStrata;
957: PetscFunctionReturn(PETSC_SUCCESS);
958: }
960: /*@
961: DMPlexTransformSetMatchStrata - Set the flag which determines what points get added to the transformed labels
963: Not Collective
965: Input Parameters:
966: + tr - The `DMPlexTransform`
967: - match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels
969: Level: intermediate
971: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetMatchStrata()`, `DMPlexGetPointDepth()`
972: @*/
973: PetscErrorCode DMPlexTransformSetMatchStrata(DMPlexTransform tr, PetscBool match)
974: {
975: PetscFunctionBegin;
977: tr->labelMatchStrata = match;
978: PetscFunctionReturn(PETSC_SUCCESS);
979: }
981: /*@
982: DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.
984: Not Collective
986: Input Parameters:
987: + tr - The `DMPlexTransform`
988: . ct - The type of the original point which produces the new point
989: . ctNew - The type of the new point
990: . p - The original point which produces the new point
991: - r - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`
993: Output Parameter:
994: . pNew - The new point number
996: Level: developer
998: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
999: @*/
1000: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
1001: {
1002: DMPolytopeType *rct;
1003: PetscInt *rsize, *cone, *ornt;
1004: PetscInt rt, Nct, n, off, rp;
1005: DMLabel trType = tr->trType;
1006: PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
1007: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
1008: PetscInt newp = ctSN, cind;
1010: PetscFunctionBeginHot;
1011: 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);
1012: PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
1013: if (trType) {
1014: PetscCall(DMLabelGetValueIndex(trType, rt, &cind));
1015: PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
1016: 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);
1017: } else {
1018: cind = ct;
1019: rp = p - ctS;
1020: }
1021: off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
1022: 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);
1023: newp += off;
1024: for (n = 0; n < Nct; ++n) {
1025: if (rct[n] == ctNew) {
1026: if (rsize[n] && r >= rsize[n])
1027: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ") for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
1028: newp += rp * rsize[n] + r;
1029: break;
1030: }
1031: }
1033: 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);
1034: *pNew = newp;
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1038: /*@
1039: DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.
1041: Not Collective
1043: Input Parameters:
1044: + tr - The `DMPlexTransform`
1045: - pNew - The new point number
1047: Output Parameters:
1048: + ct - The type of the original point which produces the new point
1049: . ctNew - The type of the new point
1050: . p - The original point which produces the new point
1051: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
1053: Level: developer
1055: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
1056: @*/
1057: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
1058: {
1059: DMLabel trType = tr->trType;
1060: DMPolytopeType *rct, ctN;
1061: PetscInt *rsize, *cone, *ornt;
1062: PetscInt rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
1063: PetscInt offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;
1065: PetscFunctionBegin;
1066: PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
1067: if (trType) {
1068: DM dm;
1069: IS rtIS;
1070: const PetscInt *reftypes;
1071: PetscInt Nrt, r, rtStart;
1073: PetscCall(DMPlexTransformGetDM(tr, &dm));
1074: PetscCall(DMLabelGetNumValues(trType, &Nrt));
1075: PetscCall(DMLabelGetValueIS(trType, &rtIS));
1076: PetscCall(ISGetIndices(rtIS, &reftypes));
1077: for (r = 0; r < Nrt; ++r) {
1078: const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];
1080: if (tr->ctStartNew[ctN] + off > pNew) continue;
1081: /* Check that any of this refinement type exist */
1082: /* TODO Actually keep track of the number produced here instead */
1083: if (off > offset) {
1084: rt = reftypes[r];
1085: offset = off;
1086: }
1087: }
1088: PetscCall(ISRestoreIndices(rtIS, &reftypes));
1089: PetscCall(ISDestroy(&rtIS));
1090: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
1091: /* TODO Map refinement types to cell types */
1092: PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL));
1093: PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
1094: for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
1095: PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
1097: if ((rtStart >= ctS) && (rtStart < ctE)) break;
1098: }
1099: PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt);
1100: } else {
1101: for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
1102: const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];
1104: if (tr->ctStartNew[ctN] + off > pNew) continue;
1105: if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
1106: /* TODO Actually keep track of the number produced here instead */
1107: if (off > offset) {
1108: ctO = ctTmp;
1109: offset = off;
1110: }
1111: }
1112: rt = -1;
1113: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
1114: }
1115: ctS = tr->ctStart[ctO];
1116: ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
1117: if (trType) {
1118: for (rtS = ctS; rtS < ctE; ++rtS) {
1119: PetscInt val;
1120: PetscCall(DMLabelGetValue(trType, rtS, &val));
1121: if (val == rt) break;
1122: }
1123: PetscCheck(rtS < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find point of type %s with refine type %" PetscInt_FMT, DMPolytopeTypes[ctO], rt);
1124: } else rtS = ctS;
1125: PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
1126: 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);
1127: for (n = 0; n < Nct; ++n) {
1128: if (rct[n] == ctN) {
1129: PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c;
1131: if (trType) {
1132: for (c = ctS; c < ctE; ++c) {
1133: PetscCall(DMLabelGetValue(trType, c, &val));
1134: if (val == rt) {
1135: if (tmp < rsize[n]) break;
1136: tmp -= rsize[n];
1137: }
1138: }
1139: PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew);
1140: rp = c - ctS;
1141: rO = tmp;
1142: } else {
1143: // This assumes that all points of type ctO transform the same way
1144: rp = tmp / rsize[n];
1145: rO = tmp % rsize[n];
1146: }
1147: break;
1148: }
1149: }
1150: PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
1151: pO = rp + ctS;
1152: 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);
1153: if (ct) *ct = (DMPolytopeType)ctO;
1154: if (ctNew) *ctNew = ctN;
1155: if (p) *p = pO;
1156: if (r) *r = rO;
1157: PetscFunctionReturn(PETSC_SUCCESS);
1158: }
1160: /*@
1161: DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.
1163: Input Parameters:
1164: + tr - The `DMPlexTransform` object
1165: . source - The source cell type
1166: - p - The source point, which can also determine the refine type
1168: Output Parameters:
1169: + rt - The refine type for this point
1170: . Nt - The number of types produced by this point
1171: . target - An array of length `Nt` giving the types produced
1172: . size - An array of length `Nt` giving the number of cells of each type produced
1173: . cone - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
1174: - ornt - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point
1176: Level: advanced
1178: Notes:
1179: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
1180: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
1181: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
1182: .vb
1183: the number of cones to be taken, or 0 for the current cell
1184: the cell cone point number at each level from which it is subdivided
1185: the replica number r of the subdivision.
1186: .ve
1187: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
1188: .vb
1189: Nt = 2
1190: target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
1191: size = {1, 2}
1192: cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
1193: ornt = { 0, 0, 0, 0}
1194: .ve
1196: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
1197: @*/
1198: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1199: {
1200: PetscFunctionBegin;
1201: PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1206: {
1207: PetscFunctionBegin;
1208: *rnew = r;
1209: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: /* Returns the same thing */
1214: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1215: {
1216: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
1217: static PetscInt vertexS[] = {1};
1218: static PetscInt vertexC[] = {0};
1219: static PetscInt vertexO[] = {0};
1220: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
1221: static PetscInt edgeS[] = {1};
1222: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1223: static PetscInt edgeO[] = {0, 0};
1224: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
1225: static PetscInt tedgeS[] = {1};
1226: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1227: static PetscInt tedgeO[] = {0, 0};
1228: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
1229: static PetscInt triS[] = {1};
1230: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1231: static PetscInt triO[] = {0, 0, 0};
1232: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
1233: static PetscInt quadS[] = {1};
1234: 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};
1235: static PetscInt quadO[] = {0, 0, 0, 0};
1236: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
1237: static PetscInt tquadS[] = {1};
1238: 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};
1239: static PetscInt tquadO[] = {0, 0, 0, 0};
1240: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
1241: static PetscInt tetS[] = {1};
1242: 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};
1243: static PetscInt tetO[] = {0, 0, 0, 0};
1244: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
1245: static PetscInt hexS[] = {1};
1246: 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};
1247: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
1248: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
1249: static PetscInt tripS[] = {1};
1250: 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};
1251: static PetscInt tripO[] = {0, 0, 0, 0, 0};
1252: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
1253: static PetscInt ttripS[] = {1};
1254: 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};
1255: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
1256: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
1257: static PetscInt tquadpS[] = {1};
1258: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
1259: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1260: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
1261: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
1262: static PetscInt pyrS[] = {1};
1263: 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};
1264: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
1266: PetscFunctionBegin;
1267: if (rt) *rt = 0;
1268: switch (source) {
1269: case DM_POLYTOPE_POINT:
1270: *Nt = 1;
1271: *target = vertexT;
1272: *size = vertexS;
1273: *cone = vertexC;
1274: *ornt = vertexO;
1275: break;
1276: case DM_POLYTOPE_SEGMENT:
1277: *Nt = 1;
1278: *target = edgeT;
1279: *size = edgeS;
1280: *cone = edgeC;
1281: *ornt = edgeO;
1282: break;
1283: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1284: *Nt = 1;
1285: *target = tedgeT;
1286: *size = tedgeS;
1287: *cone = tedgeC;
1288: *ornt = tedgeO;
1289: break;
1290: case DM_POLYTOPE_TRIANGLE:
1291: *Nt = 1;
1292: *target = triT;
1293: *size = triS;
1294: *cone = triC;
1295: *ornt = triO;
1296: break;
1297: case DM_POLYTOPE_QUADRILATERAL:
1298: *Nt = 1;
1299: *target = quadT;
1300: *size = quadS;
1301: *cone = quadC;
1302: *ornt = quadO;
1303: break;
1304: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1305: *Nt = 1;
1306: *target = tquadT;
1307: *size = tquadS;
1308: *cone = tquadC;
1309: *ornt = tquadO;
1310: break;
1311: case DM_POLYTOPE_TETRAHEDRON:
1312: *Nt = 1;
1313: *target = tetT;
1314: *size = tetS;
1315: *cone = tetC;
1316: *ornt = tetO;
1317: break;
1318: case DM_POLYTOPE_HEXAHEDRON:
1319: *Nt = 1;
1320: *target = hexT;
1321: *size = hexS;
1322: *cone = hexC;
1323: *ornt = hexO;
1324: break;
1325: case DM_POLYTOPE_TRI_PRISM:
1326: *Nt = 1;
1327: *target = tripT;
1328: *size = tripS;
1329: *cone = tripC;
1330: *ornt = tripO;
1331: break;
1332: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1333: *Nt = 1;
1334: *target = ttripT;
1335: *size = ttripS;
1336: *cone = ttripC;
1337: *ornt = ttripO;
1338: break;
1339: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1340: *Nt = 1;
1341: *target = tquadpT;
1342: *size = tquadpS;
1343: *cone = tquadpC;
1344: *ornt = tquadpO;
1345: break;
1346: case DM_POLYTOPE_PYRAMID:
1347: *Nt = 1;
1348: *target = pyrT;
1349: *size = pyrS;
1350: *cone = pyrC;
1351: *ornt = pyrO;
1352: break;
1353: default:
1354: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1355: }
1356: PetscFunctionReturn(PETSC_SUCCESS);
1357: }
1359: /*@
1360: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
1362: Not Collective
1364: Input Parameters:
1365: + tr - The `DMPlexTransform`
1366: . sct - The source point cell type, from whom the new cell is being produced
1367: . sp - The source point
1368: . so - The orientation of the source point in its enclosing parent
1369: . tct - The target point cell type
1370: . r - The replica number requested for the produced cell type
1371: - o - The orientation of the replica
1373: Output Parameters:
1374: + rnew - The replica number, given the orientation of the parent
1375: - onew - The replica orientation, given the orientation of the parent
1377: Level: advanced
1379: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1380: @*/
1381: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1382: {
1383: PetscFunctionBeginHot;
1384: PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1385: PetscFunctionReturn(PETSC_SUCCESS);
1386: }
1388: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1389: {
1390: DM dm;
1391: PetscInt pStart, pEnd, pNew;
1393: PetscFunctionBegin;
1394: PetscCall(DMPlexTransformGetDM(tr, &dm));
1395: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0));
1396: /* Must create the celltype label here so that we do not automatically try to compute the types */
1397: PetscCall(DMCreateLabel(rdm, "celltype"));
1398: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1399: for (PetscInt p = pStart; p < pEnd; ++p) {
1400: DMPolytopeType ct;
1401: DMPolytopeType *rct;
1402: PetscInt *rsize, *rcone, *rornt;
1403: PetscInt Nct, n, r;
1405: PetscCall(DMPlexGetCellType(dm, p, &ct));
1406: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1407: for (n = 0; n < Nct; ++n) {
1408: for (r = 0; r < rsize[n]; ++r) {
1409: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1410: PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
1411: PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
1412: }
1413: }
1414: }
1415: /* Let the DM know we have set all the cell types */
1416: {
1417: DMLabel ctLabel;
1418: DM_Plex *plex = (DM_Plex *)rdm->data;
1420: PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
1421: PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
1422: }
1423: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0));
1424: PetscFunctionReturn(PETSC_SUCCESS);
1425: }
1427: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1428: {
1429: DMPolytopeType ctNew;
1431: PetscFunctionBegin;
1433: PetscAssertPointer(coneSize, 3);
1434: PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
1435: *coneSize = DMPolytopeTypeGetConeSize(ctNew);
1436: PetscFunctionReturn(PETSC_SUCCESS);
1437: }
1439: /* The orientation o is for the interior of the cell p */
1440: 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[])
1441: {
1442: DM dm;
1443: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1444: const PetscInt *cone;
1445: DMPolytopeType *newft = NULL;
1446: PetscInt c, coff = *coneoff, ooff = *orntoff;
1447: PetscInt dim, cr = 0, co = 0, nr, no;
1449: PetscFunctionBegin;
1450: PetscCall(DMPlexTransformGetDM(tr, &dm));
1451: PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
1452: // Check if we have to permute this cell
1453: PetscCall(DMGetDimension(dm, &dim));
1454: if (DMPolytopeTypeGetDim(ctNew) == dim && DMPolytopeTypeGetDim(ct) == dim - 1) {
1455: PetscCall(DMPlexTransformGetSubcellOrientation(tr, ct, p, o, ctNew, cr, co, &nr, &no));
1456: if (cr != nr || co != no) PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newft));
1457: }
1458: for (c = 0; c < csizeNew; ++c) {
1459: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1460: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1461: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1462: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1463: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1464: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1465: const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1466: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1467: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1468: PetscInt lc;
1470: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1471: for (lc = 0; lc < fn; ++lc) {
1472: const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po);
1473: const PetscInt acp = rcone[coff++];
1474: const PetscInt pcp = parr[acp * 2];
1475: const PetscInt pco = parr[acp * 2 + 1];
1476: const PetscInt *ppornt;
1478: ppp = pp;
1479: pp = pcone[pcp];
1480: PetscCall(DMPlexGetCellType(dm, pp, &pct));
1481: // Restore the parent cone from the last iterate
1482: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
1483: PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
1484: PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
1485: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1486: PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
1487: }
1488: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
1489: pr = rcone[coff++];
1490: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1491: PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
1492: PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
1493: orntNew[c] = fo;
1494: if (newft) newft[c] = ft;
1495: }
1496: PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
1497: if (newft) {
1498: const PetscInt *arr;
1499: PetscInt *newcone, *newornt;
1501: arr = DMPolytopeTypeGetArrangement(ctNew, no);
1502: PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newcone));
1503: PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newornt));
1504: for (PetscInt c = 0; c < csizeNew; ++c) {
1505: DMPolytopeType ft = newft[c];
1506: PetscInt nO;
1508: nO = DMPolytopeTypeGetNumArrangements(ft) / 2;
1509: newcone[c] = coneNew[arr[c * 2 + 0]];
1510: newornt[c] = DMPolytopeTypeComposeOrientation(ft, arr[c * 2 + 1], orntNew[arr[c * 2 + 0]]);
1511: 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]);
1512: }
1513: for (PetscInt c = 0; c < csizeNew; ++c) {
1514: coneNew[c] = newcone[c];
1515: orntNew[c] = newornt[c];
1516: }
1517: PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newcone));
1518: PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newornt));
1519: PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newft));
1520: }
1521: *coneoff = coff;
1522: *orntoff = ooff;
1523: PetscFunctionReturn(PETSC_SUCCESS);
1524: }
1526: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1527: {
1528: DM dm;
1529: DMPolytopeType ct;
1530: PetscInt *coneNew, *orntNew;
1531: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1533: PetscFunctionBegin;
1534: PetscCall(DMPlexTransformGetDM(tr, &dm));
1535: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0));
1536: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1537: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1538: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1539: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1540: for (p = pStart; p < pEnd; ++p) {
1541: PetscInt coff, ooff;
1542: DMPolytopeType *rct;
1543: PetscInt *rsize, *rcone, *rornt;
1544: PetscInt Nct, n, r;
1546: PetscCall(DMPlexGetCellType(dm, p, &ct));
1547: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1548: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1549: const DMPolytopeType ctNew = rct[n];
1551: for (r = 0; r < rsize[n]; ++r) {
1552: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1553: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
1554: PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
1555: PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
1556: }
1557: }
1558: }
1559: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1560: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1561: PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
1562: PetscCall(DMPlexSymmetrize(rdm));
1563: PetscCall(DMPlexStratify(rdm));
1564: PetscTryTypeMethod(tr, ordersupports, dm, rdm);
1565: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0));
1566: PetscFunctionReturn(PETSC_SUCCESS);
1567: }
1569: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1570: {
1571: DM dm;
1572: DMPolytopeType ct, qct;
1573: DMPolytopeType *rct;
1574: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1575: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1577: PetscFunctionBegin;
1579: PetscAssertPointer(cone, 4);
1580: PetscAssertPointer(ornt, 5);
1581: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1582: PetscCall(DMPlexTransformGetDM(tr, &dm));
1583: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1584: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1585: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1586: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1587: for (n = 0; n < Nct; ++n) {
1588: const DMPolytopeType ctNew = rct[n];
1589: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1590: PetscInt Nr = rsize[n], fn, c;
1592: if (ctNew == qct) Nr = r;
1593: for (nr = 0; nr < Nr; ++nr) {
1594: for (c = 0; c < csizeNew; ++c) {
1595: ++coff; /* Cell type of new cone point */
1596: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1597: coff += fn;
1598: ++coff; /* Replica number of new cone point */
1599: ++ooff; /* Orientation of new cone point */
1600: }
1601: }
1602: if (ctNew == qct) break;
1603: }
1604: PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1605: *cone = qcone;
1606: *ornt = qornt;
1607: PetscFunctionReturn(PETSC_SUCCESS);
1608: }
1610: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1611: {
1612: DM dm;
1613: DMPolytopeType ct, qct;
1614: DMPolytopeType *rct;
1615: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1616: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1618: PetscFunctionBegin;
1620: if (cone) PetscAssertPointer(cone, 3);
1621: if (ornt) PetscAssertPointer(ornt, 4);
1622: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1623: PetscCall(DMPlexTransformGetDM(tr, &dm));
1624: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1625: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1626: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1627: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1628: for (n = 0; n < Nct; ++n) {
1629: const DMPolytopeType ctNew = rct[n];
1630: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1631: PetscInt Nr = rsize[n], fn, c;
1633: if (ctNew == qct) Nr = r;
1634: for (nr = 0; nr < Nr; ++nr) {
1635: for (c = 0; c < csizeNew; ++c) {
1636: ++coff; /* Cell type of new cone point */
1637: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1638: coff += fn;
1639: ++coff; /* Replica number of new cone point */
1640: ++ooff; /* Orientation of new cone point */
1641: }
1642: }
1643: if (ctNew == qct) break;
1644: }
1645: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1646: if (cone) *cone = qcone;
1647: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1648: if (ornt) *ornt = qornt;
1649: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1650: PetscFunctionReturn(PETSC_SUCCESS);
1651: }
1653: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1654: {
1655: DM dm;
1657: PetscFunctionBegin;
1659: PetscCall(DMPlexTransformGetDM(tr, &dm));
1660: if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
1661: if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
1662: PetscFunctionReturn(PETSC_SUCCESS);
1663: }
1665: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1666: {
1667: PetscInt ict;
1669: PetscFunctionBegin;
1670: PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
1671: for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1672: const DMPolytopeType ct = (DMPolytopeType)ict;
1673: DMPlexTransform reftr;
1674: DM refdm, trdm;
1675: Vec coordinates;
1676: const PetscScalar *coords;
1677: DMPolytopeType *rct;
1678: PetscInt *rsize, *rcone, *rornt;
1679: PetscInt Nct, n, r, pNew = 0;
1680: PetscInt trdim, vStart, vEnd, Nc;
1681: const PetscInt debug = 0;
1682: const char *typeName;
1684: /* Since points are 0-dimensional, coordinates make no sense */
1685: if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue;
1686: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
1687: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
1688: PetscCall(DMPlexTransformSetDM(reftr, refdm));
1689: PetscCall(DMPlexTransformGetType(tr, &typeName));
1690: PetscCall(DMPlexTransformSetType(reftr, typeName));
1691: PetscCall(DMPlexTransformSetUp(reftr));
1692: PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));
1694: PetscCall(DMGetDimension(trdm, &trdim));
1695: PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
1696: tr->trNv[ct] = vEnd - vStart;
1697: PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
1698: PetscCall(VecGetLocalSize(coordinates, &Nc));
1699: 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);
1700: PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
1701: PetscCall(VecGetArrayRead(coordinates, &coords));
1702: PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
1703: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1705: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
1706: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1707: for (n = 0; n < Nct; ++n) {
1708: /* Since points are 0-dimensional, coordinates make no sense */
1709: if (rct[n] == DM_POLYTOPE_POINT) continue;
1710: PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
1711: for (r = 0; r < rsize[n]; ++r) {
1712: PetscInt *closure = NULL;
1713: PetscInt clSize, cl, Nv = 0;
1715: PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
1716: PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
1717: PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1718: for (cl = 0; cl < clSize * 2; cl += 2) {
1719: const PetscInt sv = closure[cl];
1721: if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1722: }
1723: PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1724: 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]);
1725: }
1726: }
1727: if (debug) {
1728: DMPolytopeType *rct;
1729: PetscInt *rsize, *rcone, *rornt;
1730: PetscInt v, dE = trdim, d, off = 0;
1732: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
1733: for (v = 0; v < tr->trNv[ct]; ++v) {
1734: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1735: for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
1736: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1737: }
1739: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1740: for (n = 0; n < Nct; ++n) {
1741: if (rct[n] == DM_POLYTOPE_POINT) continue;
1742: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
1743: for (r = 0; r < rsize[n]; ++r) {
1744: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1745: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
1746: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1747: }
1748: }
1749: }
1750: PetscCall(DMDestroy(&refdm));
1751: PetscCall(DMDestroy(&trdm));
1752: PetscCall(DMPlexTransformDestroy(&reftr));
1753: }
1754: PetscFunctionReturn(PETSC_SUCCESS);
1755: }
1757: /*@C
1758: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
1760: Input Parameters:
1761: + tr - The `DMPlexTransform` object
1762: - ct - The cell type
1764: Output Parameters:
1765: + Nv - The number of transformed vertices in the closure of the reference cell of given type
1766: - trVerts - The coordinates of these vertices in the reference cell
1768: Level: developer
1770: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
1771: @*/
1772: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1773: {
1774: PetscFunctionBegin;
1775: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1776: if (Nv) *Nv = tr->trNv[ct];
1777: if (trVerts) *trVerts = tr->trVerts[ct];
1778: PetscFunctionReturn(PETSC_SUCCESS);
1779: }
1781: /*@C
1782: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
1784: Input Parameters:
1785: + tr - The `DMPlexTransform` object
1786: . ct - The cell type
1787: . rct - The subcell type
1788: - r - The subcell index
1790: Output Parameter:
1791: . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()`
1793: Level: developer
1795: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
1796: @*/
1797: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1798: {
1799: PetscFunctionBegin;
1800: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1801: PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
1802: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1803: PetscFunctionReturn(PETSC_SUCCESS);
1804: }
1806: /* Computes new vertex as the barycenter, or centroid */
1807: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1808: {
1809: PetscInt v, d;
1811: PetscFunctionBeginHot;
1812: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
1813: for (d = 0; d < dE; ++d) out[d] = 0.0;
1814: for (v = 0; v < Nv; ++v)
1815: for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1816: for (d = 0; d < dE; ++d) out[d] /= Nv;
1817: PetscFunctionReturn(PETSC_SUCCESS);
1818: }
1820: /*@
1821: DMPlexTransformMapCoordinates - Calculate new coordinates for produced points
1823: Not collective
1825: Input Parameters:
1826: + tr - The `DMPlexTransform`
1827: . pct - The cell type of the parent, from whom the new cell is being produced
1828: . ct - The type being produced
1829: . p - The original point
1830: . r - The replica number requested for the produced cell type
1831: . Nv - Number of vertices in the closure of the parent cell
1832: . dE - Spatial dimension
1833: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
1835: Output Parameter:
1836: . out - The coordinates of the new vertices
1838: Level: intermediate
1840: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
1841: @*/
1842: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1843: {
1844: PetscFunctionBeginHot;
1845: if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1846: PetscFunctionReturn(PETSC_SUCCESS);
1847: }
1849: /*
1850: DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label
1852: Not Collective
1854: Input Parameters:
1855: + tr - The `DMPlexTransform`
1856: . label - The label in the transformed mesh
1857: . pp - The parent point in the original mesh
1858: . pct - The cell type of the parent point
1859: . p - The point in the transformed mesh
1860: . ct - The cell type of the point
1861: . r - The replica number of the point
1862: - val - The label value of the parent point
1864: Level: developer
1866: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
1867: */
1868: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
1869: {
1870: PetscFunctionBeginHot;
1871: if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
1872: PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
1873: PetscFunctionReturn(PETSC_SUCCESS);
1874: }
1876: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1877: {
1878: DM dm;
1879: IS valueIS;
1880: const PetscInt *values;
1881: PetscInt defVal, Nv, val;
1883: PetscFunctionBegin;
1884: PetscCall(DMPlexTransformGetDM(tr, &dm));
1885: PetscCall(DMLabelGetDefaultValue(label, &defVal));
1886: PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
1887: PetscCall(DMLabelGetValueIS(label, &valueIS));
1888: PetscCall(ISGetLocalSize(valueIS, &Nv));
1889: PetscCall(ISGetIndices(valueIS, &values));
1890: for (val = 0; val < Nv; ++val) {
1891: IS pointIS;
1892: const PetscInt *points;
1893: PetscInt numPoints, p;
1895: /* Ensure refined label is created with same number of strata as
1896: * original (even if no entries here). */
1897: PetscCall(DMLabelAddStratum(labelNew, values[val]));
1898: PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
1899: PetscCall(ISGetLocalSize(pointIS, &numPoints));
1900: PetscCall(ISGetIndices(pointIS, &points));
1901: for (p = 0; p < numPoints; ++p) {
1902: const PetscInt point = points[p];
1903: DMPolytopeType ct;
1904: DMPolytopeType *rct;
1905: PetscInt *rsize, *rcone, *rornt;
1906: PetscInt Nct, n, r, pNew = 0;
1908: PetscCall(DMPlexGetCellType(dm, point, &ct));
1909: PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1910: for (n = 0; n < Nct; ++n) {
1911: for (r = 0; r < rsize[n]; ++r) {
1912: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
1913: PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
1914: }
1915: }
1916: }
1917: PetscCall(ISRestoreIndices(pointIS, &points));
1918: PetscCall(ISDestroy(&pointIS));
1919: }
1920: PetscCall(ISRestoreIndices(valueIS, &values));
1921: PetscCall(ISDestroy(&valueIS));
1922: PetscFunctionReturn(PETSC_SUCCESS);
1923: }
1925: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1926: {
1927: DM dm;
1928: PetscInt numLabels, l;
1930: PetscFunctionBegin;
1931: PetscCall(DMPlexTransformGetDM(tr, &dm));
1932: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0));
1933: PetscCall(DMGetNumLabels(dm, &numLabels));
1934: for (l = 0; l < numLabels; ++l) {
1935: DMLabel label, labelNew;
1936: const char *lname;
1937: PetscBool isDepth, isCellType;
1939: PetscCall(DMGetLabelName(dm, l, &lname));
1940: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
1941: if (isDepth) continue;
1942: PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
1943: if (isCellType) continue;
1944: PetscCall(DMCreateLabel(rdm, lname));
1945: PetscCall(DMGetLabel(dm, lname, &label));
1946: PetscCall(DMGetLabel(rdm, lname, &labelNew));
1947: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1948: }
1949: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0));
1950: PetscFunctionReturn(PETSC_SUCCESS);
1951: }
1953: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1954: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1955: {
1956: DM dm;
1957: PetscInt Nf, f, Nds, s;
1959: PetscFunctionBegin;
1960: PetscCall(DMPlexTransformGetDM(tr, &dm));
1961: PetscCall(DMGetNumFields(dm, &Nf));
1962: for (f = 0; f < Nf; ++f) {
1963: DMLabel label, labelNew;
1964: PetscObject obj;
1965: const char *lname;
1967: PetscCall(DMGetField(rdm, f, &label, &obj));
1968: if (!label) continue;
1969: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1970: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1971: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1972: PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
1973: PetscCall(DMLabelDestroy(&labelNew));
1974: }
1975: PetscCall(DMGetNumDS(dm, &Nds));
1976: for (s = 0; s < Nds; ++s) {
1977: DMLabel label, labelNew;
1978: const char *lname;
1980: PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
1981: if (!label) continue;
1982: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1983: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1984: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1985: PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
1986: PetscCall(DMLabelDestroy(&labelNew));
1987: }
1988: PetscFunctionReturn(PETSC_SUCCESS);
1989: }
1991: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1992: {
1993: DM dm;
1994: PetscSF sf, sfNew;
1995: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
1996: const PetscInt *localPoints;
1997: const PetscSFNode *remotePoints;
1998: PetscInt *localPointsNew;
1999: PetscSFNode *remotePointsNew;
2000: PetscInt pStartNew, pEndNew, pNew;
2001: /* Brute force algorithm */
2002: PetscSF rsf;
2003: PetscSection s;
2004: const PetscInt *rootdegree;
2005: PetscInt *rootPointsNew, *remoteOffsets;
2006: PetscInt numPointsNew, pStart, pEnd, p;
2008: PetscFunctionBegin;
2009: PetscCall(DMPlexTransformGetDM(tr, &dm));
2010: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
2011: PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
2012: PetscCall(DMGetPointSF(dm, &sf));
2013: PetscCall(DMGetPointSF(rdm, &sfNew));
2014: /* Calculate size of new SF */
2015: PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
2016: if (numRoots < 0) {
2017: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
2018: PetscFunctionReturn(PETSC_SUCCESS);
2019: }
2020: for (l = 0; l < numLeaves; ++l) {
2021: const PetscInt p = localPoints[l];
2022: DMPolytopeType ct;
2023: DMPolytopeType *rct;
2024: PetscInt *rsize, *rcone, *rornt;
2025: PetscInt Nct, n;
2027: PetscCall(DMPlexGetCellType(dm, p, &ct));
2028: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2029: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
2030: }
2031: /* Send new root point numbers
2032: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
2033: */
2034: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2035: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
2036: PetscCall(PetscSectionSetChart(s, pStart, pEnd));
2037: for (p = pStart; p < pEnd; ++p) {
2038: DMPolytopeType ct;
2039: DMPolytopeType *rct;
2040: PetscInt *rsize, *rcone, *rornt;
2041: PetscInt Nct, n;
2043: PetscCall(DMPlexGetCellType(dm, p, &ct));
2044: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2045: for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
2046: }
2047: PetscCall(PetscSectionSetUp(s));
2048: PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
2049: PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
2050: PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
2051: PetscCall(PetscFree(remoteOffsets));
2052: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
2053: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
2054: PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
2055: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
2056: for (p = pStart; p < pEnd; ++p) {
2057: DMPolytopeType ct;
2058: DMPolytopeType *rct;
2059: PetscInt *rsize, *rcone, *rornt;
2060: PetscInt Nct, n, r, off;
2062: if (!rootdegree[p - pStart]) continue;
2063: PetscCall(PetscSectionGetOffset(s, p, &off));
2064: PetscCall(DMPlexGetCellType(dm, p, &ct));
2065: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2066: for (n = 0, m = 0; n < Nct; ++n) {
2067: for (r = 0; r < rsize[n]; ++r, ++m) {
2068: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2069: rootPointsNew[off + m] = pNew;
2070: }
2071: }
2072: }
2073: PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
2074: PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
2075: PetscCall(PetscSFDestroy(&rsf));
2076: PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
2077: PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
2078: for (l = 0, m = 0; l < numLeaves; ++l) {
2079: const PetscInt p = localPoints[l];
2080: DMPolytopeType ct;
2081: DMPolytopeType *rct;
2082: PetscInt *rsize, *rcone, *rornt;
2083: PetscInt Nct, n, r, q, off;
2085: PetscCall(PetscSectionGetOffset(s, p, &off));
2086: PetscCall(DMPlexGetCellType(dm, p, &ct));
2087: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2088: for (n = 0, q = 0; n < Nct; ++n) {
2089: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
2090: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2091: localPointsNew[m] = pNew;
2092: remotePointsNew[m].index = rootPointsNew[off + q];
2093: remotePointsNew[m].rank = remotePoints[l].rank;
2094: }
2095: }
2096: }
2097: PetscCall(PetscSectionDestroy(&s));
2098: PetscCall(PetscFree(rootPointsNew));
2099: /* SF needs sorted leaves to correctly calculate Gather */
2100: {
2101: PetscSFNode *rp, *rtmp;
2102: PetscInt *lp, *idx, *ltmp, i;
2104: PetscCall(PetscMalloc1(numLeavesNew, &idx));
2105: PetscCall(PetscMalloc1(numLeavesNew, &lp));
2106: PetscCall(PetscMalloc1(numLeavesNew, &rp));
2107: for (i = 0; i < numLeavesNew; ++i) {
2108: 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);
2109: idx[i] = i;
2110: }
2111: PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
2112: for (i = 0; i < numLeavesNew; ++i) {
2113: lp[i] = localPointsNew[idx[i]];
2114: rp[i] = remotePointsNew[idx[i]];
2115: }
2116: ltmp = localPointsNew;
2117: localPointsNew = lp;
2118: rtmp = remotePointsNew;
2119: remotePointsNew = rp;
2120: PetscCall(PetscFree(idx));
2121: PetscCall(PetscFree(ltmp));
2122: PetscCall(PetscFree(rtmp));
2123: }
2124: PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
2125: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
2126: PetscFunctionReturn(PETSC_SUCCESS);
2127: }
2129: /*
2130: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`.
2132: Not Collective
2134: Input Parameters:
2135: + tr - The `DMPlexTransform`
2136: . ct - The type of the parent cell
2137: . rct - The type of the produced cell
2138: . r - The index of the produced cell
2139: - x - The localized coordinates for the parent cell
2141: Output Parameter:
2142: . xr - The localized coordinates for the produced cell
2144: Level: developer
2146: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
2147: */
2148: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
2149: {
2150: PetscFE fe = NULL;
2151: PetscInt cdim, v, *subcellV;
2153: PetscFunctionBegin;
2154: PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
2155: PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
2156: PetscCall(PetscFEGetNumComponents(fe, &cdim));
2157: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
2158: PetscFunctionReturn(PETSC_SUCCESS);
2159: }
2161: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
2162: {
2163: DM dm, cdm, cdmCell, cdmNew, cdmCellNew;
2164: PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
2165: Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
2166: const PetscScalar *coords;
2167: PetscScalar *coordsNew;
2168: const PetscReal *maxCell, *Lstart, *L;
2169: PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
2170: PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;
2172: PetscFunctionBegin;
2173: // Need to clear the DMField for coordinates
2174: PetscCall(DMSetCoordinateField(rdm, NULL));
2175: PetscCall(DMPlexTransformGetDM(tr, &dm));
2176: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0));
2177: PetscCall(DMGetCoordinateDM(dm, &cdm));
2178: PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
2179: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
2180: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
2181: if (localized) {
2182: /* Localize coordinates of new vertices */
2183: localizeVertices = PETSC_TRUE;
2184: /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
2185: if (!maxCell) localizeCells = PETSC_TRUE;
2186: }
2187: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2188: PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
2189: if (maxCell) {
2190: PetscReal maxCellNew[3];
2192: for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
2193: PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L));
2194: }
2195: PetscCall(DMGetCoordinateDim(rdm, &dE));
2196: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
2197: PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
2198: PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
2199: PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
2200: PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
2201: /* Localization should be inherited */
2202: /* Stefano calculates parent cells for each new cell for localization */
2203: /* Localized cells need coordinates of closure */
2204: for (v = vStartNew; v < vEndNew; ++v) {
2205: PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
2206: PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
2207: }
2208: PetscCall(PetscSectionSetUp(coordSectionNew));
2209: PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));
2211: if (localizeCells) {
2212: PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
2213: PetscCall(DMClone(cdmNew, &cdmCellNew));
2214: PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
2215: PetscCall(DMDestroy(&cdmCellNew));
2217: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
2218: PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
2219: PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
2220: PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
2221: PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));
2223: PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
2224: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2225: for (c = cStart; c < cEnd; ++c) {
2226: PetscInt dof;
2228: PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
2229: if (dof) {
2230: DMPolytopeType ct;
2231: DMPolytopeType *rct;
2232: PetscInt *rsize, *rcone, *rornt;
2233: PetscInt dim, cNew, Nct, n, r;
2235: PetscCall(DMPlexGetCellType(dm, c, &ct));
2236: dim = DMPolytopeTypeGetDim(ct);
2237: PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2238: /* This allows for different cell types */
2239: for (n = 0; n < Nct; ++n) {
2240: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2241: for (r = 0; r < rsize[n]; ++r) {
2242: PetscInt *closure = NULL;
2243: PetscInt clSize, cl, Nv = 0;
2245: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
2246: PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2247: for (cl = 0; cl < clSize * 2; cl += 2) {
2248: if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
2249: }
2250: PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2251: PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
2252: PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
2253: }
2254: }
2255: }
2256: }
2257: PetscCall(PetscSectionSetUp(coordSectionCellNew));
2258: PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
2259: }
2260: PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
2261: {
2262: VecType vtype;
2263: PetscInt coordSizeNew, bs;
2264: const char *name;
2266: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
2267: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
2268: PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
2269: PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
2270: PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
2271: PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
2272: PetscCall(VecGetBlockSize(coordsLocal, &bs));
2273: PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
2274: PetscCall(VecGetType(coordsLocal, &vtype));
2275: PetscCall(VecSetType(coordsLocalNew, vtype));
2276: }
2277: PetscCall(VecGetArrayRead(coordsLocal, &coords));
2278: PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
2279: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2280: /* First set coordinates for vertices */
2281: for (p = pStart; p < pEnd; ++p) {
2282: DMPolytopeType ct;
2283: DMPolytopeType *rct;
2284: PetscInt *rsize, *rcone, *rornt;
2285: PetscInt Nct, n, r;
2286: PetscBool hasVertex = PETSC_FALSE;
2288: PetscCall(DMPlexGetCellType(dm, p, &ct));
2289: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2290: for (n = 0; n < Nct; ++n) {
2291: if (rct[n] == DM_POLYTOPE_POINT) {
2292: hasVertex = PETSC_TRUE;
2293: break;
2294: }
2295: }
2296: if (hasVertex) {
2297: const PetscScalar *icoords = NULL;
2298: const PetscScalar *array = NULL;
2299: PetscScalar *pcoords = NULL;
2300: PetscBool isDG;
2301: PetscInt Nc, Nv, v, d;
2303: PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2305: icoords = pcoords;
2306: Nv = Nc / dEo;
2307: if (ct != DM_POLYTOPE_POINT) {
2308: if (localizeVertices && maxCell) {
2309: PetscScalar anchor[3];
2311: for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
2312: for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
2313: }
2314: }
2315: for (n = 0; n < Nct; ++n) {
2316: if (rct[n] != DM_POLYTOPE_POINT) continue;
2317: for (r = 0; r < rsize[n]; ++r) {
2318: PetscScalar vcoords[3];
2319: PetscInt vNew, off;
2321: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
2322: PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
2323: PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
2324: PetscCall(DMSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
2325: }
2326: }
2327: PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2328: }
2329: }
2330: PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
2331: PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
2332: PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
2333: PetscCall(VecDestroy(&coordsLocalNew));
2334: PetscCall(PetscSectionDestroy(&coordSectionNew));
2335: /* Then set coordinates for cells by localizing */
2336: if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
2337: else {
2338: VecType vtype;
2339: PetscInt coordSizeNew, bs;
2340: const char *name;
2342: PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
2343: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
2344: PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
2345: PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
2346: PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
2347: PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
2348: PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
2349: PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
2350: PetscCall(VecGetType(coordsLocalCell, &vtype));
2351: PetscCall(VecSetType(coordsLocalCellNew, vtype));
2352: PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
2353: PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));
2355: for (p = pStart; p < pEnd; ++p) {
2356: DMPolytopeType ct;
2357: DMPolytopeType *rct;
2358: PetscInt *rsize, *rcone, *rornt;
2359: PetscInt dof = 0, Nct, n, r;
2361: PetscCall(DMPlexGetCellType(dm, p, &ct));
2362: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2363: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
2364: if (dof) {
2365: const PetscScalar *pcoords;
2367: PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords));
2368: for (n = 0; n < Nct; ++n) {
2369: const PetscInt Nr = rsize[n];
2371: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2372: for (r = 0; r < Nr; ++r) {
2373: PetscInt pNew, offNew;
2375: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2376: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2377: cell to the ones it produces. */
2378: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2379: PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
2380: PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
2381: }
2382: }
2383: }
2384: }
2385: PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
2386: PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
2387: PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
2388: PetscCall(VecDestroy(&coordsLocalCellNew));
2389: PetscCall(PetscSectionDestroy(&coordSectionCellNew));
2390: }
2391: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0));
2392: PetscFunctionReturn(PETSC_SUCCESS);
2393: }
2395: /*@
2396: DMPlexTransformApply - Execute the transformation, producing another `DM`
2398: Collective
2400: Input Parameters:
2401: + tr - The `DMPlexTransform` object
2402: - dm - The original `DM`
2404: Output Parameter:
2405: . tdm - The transformed `DM`
2407: Level: intermediate
2409: Options Database Keys:
2410: + -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point
2411: . -dm_plex_transform_label_replica_inc <num> - Increment for the label value to be multiplied by the replica number
2412: - -dm_plex_transform_active <name> - Name for active mesh label
2414: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()`
2415: @*/
2416: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
2417: {
2418: DM rdm;
2419: DMPlexInterpolatedFlag interp;
2420: PetscInt pStart, pEnd;
2422: PetscFunctionBegin;
2425: PetscAssertPointer(tdm, 3);
2426: PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0));
2427: PetscCall(DMPlexTransformSetDM(tr, dm));
2429: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
2430: PetscCall(DMSetType(rdm, DMPLEX));
2431: PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
2432: /* Calculate number of new points of each depth */
2433: PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
2434: PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
2435: /* Step 1: Set chart */
2436: PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
2437: PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
2438: /* Step 2: Set cone/support sizes (automatically stratifies) */
2439: PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
2440: /* Step 3: Setup refined DM */
2441: PetscCall(DMSetUp(rdm));
2442: /* Step 4: Set cones and supports (automatically symmetrizes) */
2443: PetscCall(DMPlexTransformSetCones(tr, rdm));
2444: /* Step 5: Create pointSF */
2445: PetscCall(DMPlexTransformCreateSF(tr, rdm));
2446: /* Step 6: Create labels */
2447: PetscCall(DMPlexTransformCreateLabels(tr, rdm));
2448: /* Step 7: Set coordinates */
2449: PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
2450: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm));
2451: // If the original DM was configured from options, the transformed DM should be as well
2452: rdm->setfromoptionscalled = dm->setfromoptionscalled;
2453: PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0));
2454: *tdm = rdm;
2455: PetscFunctionReturn(PETSC_SUCCESS);
2456: }
2458: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2459: {
2460: DMPlexTransform tr;
2461: DM cdm, rcdm;
2462: const char *prefix;
2463: PetscBool save;
2465: PetscFunctionBegin;
2466: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
2467: PetscCall(PetscObjectSetName((PetscObject)tr, "Adapt Label Transform"));
2468: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2469: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
2470: PetscCall(DMPlexTransformSetDM(tr, dm));
2471: PetscCall(DMPlexTransformSetFromOptions(tr));
2472: if (adaptLabel) PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
2473: PetscCall(DMPlexTransformSetUp(tr));
2474: PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
2475: PetscCall(DMPlexTransformApply(tr, dm, rdm));
2476: PetscCall(DMCopyDisc(dm, *rdm));
2477: PetscCall(DMGetCoordinateDM(dm, &cdm));
2478: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
2479: PetscCall(DMCopyDisc(cdm, rcdm));
2480: PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
2481: PetscCall(DMCopyDisc(dm, *rdm));
2482: PetscCall(DMPlexGetSaveTransform(dm, &save));
2483: if (save) PetscCall(DMPlexSetTransform(*rdm, tr));
2484: PetscCall(DMPlexTransformDestroy(&tr));
2485: ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2486: PetscFunctionReturn(PETSC_SUCCESS);
2487: }