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