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