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