Actual source code: plextransform.c
1: #include <petsc/private/dmplextransformimpl.h>
3: #include <petsc/private/petscfeimpl.h>
5: PetscClassId DMPLEXTRANSFORM_CLASSID;
7: PetscFunctionList DMPlexTransformList = NULL;
8: PetscBool DMPlexTransformRegisterAllCalled = PETSC_FALSE;
10: /* Construct cell type order since we must loop over cell types in depth order */
11: static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
12: {
13: PetscInt *ctO, *ctOInv;
14: PetscInt c, d, off = 0;
16: PetscFunctionBegin;
17: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv));
18: for (d = 3; d >= dim; --d) {
19: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
20: if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d || c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue;
21: ctO[off++] = c;
22: }
23: }
24: if (dim != 0) {
25: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
26: if (DMPolytopeTypeGetDim((DMPolytopeType)c) != 0) continue;
27: ctO[off++] = c;
28: }
29: }
30: for (d = dim - 1; d > 0; --d) {
31: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
32: if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d || c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue;
33: ctO[off++] = c;
34: }
35: }
36: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
37: if (DMPolytopeTypeGetDim((DMPolytopeType)c) >= 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) continue;
38: ctO[off++] = c;
39: }
40: PetscCheck(off == DM_NUM_POLYTOPES + 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %" PetscInt_FMT " for cell type order", off);
41: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c;
42: *ctOrder = ctO;
43: *ctOrderInv = ctOInv;
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: /*@C
48: DMPlexTransformRegister - Adds a new transform component implementation
50: Not Collective
52: Input Parameters:
53: + name - The name of a new user-defined creation routine
54: - create_func - The creation routine
56: Example Usage:
57: .vb
58: DMPlexTransformRegister("my_transform", MyTransformCreate);
59: .ve
61: Then, your transform type can be chosen with the procedural interface via
62: .vb
63: DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
64: DMPlexTransformSetType(DMPlexTransform, "my_transform");
65: .ve
66: or at runtime via the option
67: .vb
68: -dm_plex_transform_type my_transform
69: .ve
71: Level: advanced
73: Note:
74: `DMPlexTransformRegister()` may be called multiple times to add several user-defined transforms
76: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()`
77: @*/
78: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
79: {
80: PetscFunctionBegin;
81: PetscCall(DMInitializePackage());
82: PetscCall(PetscFunctionListAdd(&DMPlexTransformList, name, create_func));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
87: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
88: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
89: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
90: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
91: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
92: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform);
93: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform);
95: /*@C
96: DMPlexTransformRegisterAll - Registers all of the transform components in the `DM` package.
98: Not Collective
100: Level: advanced
102: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
103: @*/
104: PetscErrorCode DMPlexTransformRegisterAll(void)
105: {
106: PetscFunctionBegin;
107: if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
108: DMPlexTransformRegisterAllCalled = PETSC_TRUE;
110: PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter));
111: PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular));
112: PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox));
113: PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld));
114: PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL));
115: PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR));
116: PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D));
117: PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDE, DMPlexTransformCreate_Extrude));
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: /*@C
122: DMPlexTransformRegisterDestroy - This function destroys the registered `DMPlexTransformType`. It is called from `PetscFinalize()`.
124: Level: developer
126: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()`
127: @*/
128: PetscErrorCode DMPlexTransformRegisterDestroy(void)
129: {
130: PetscFunctionBegin;
131: PetscCall(PetscFunctionListDestroy(&DMPlexTransformList));
132: DMPlexTransformRegisterAllCalled = PETSC_FALSE;
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: /*@
137: DMPlexTransformCreate - Creates an empty transform object. The type can then be set with `DMPlexTransformSetType()`.
139: Collective
141: Input Parameter:
142: . comm - The communicator for the transform object
144: Output Parameter:
145: . tr - The transform object
147: Level: beginner
149: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
150: @*/
151: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
152: {
153: DMPlexTransform t;
155: PetscFunctionBegin;
156: PetscAssertPointer(tr, 2);
157: *tr = NULL;
158: PetscCall(DMInitializePackage());
160: PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView));
161: t->setupcalled = PETSC_FALSE;
162: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom));
163: *tr = t;
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: /*@C
168: DMPlexTransformSetType - Sets the particular implementation for a transform.
170: Collective
172: Input Parameters:
173: + tr - The transform
174: - method - The name of the transform type
176: Options Database Key:
177: . -dm_plex_transform_type <type> - Sets the transform type; see `DMPlexTransformType`
179: Level: intermediate
181: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
182: @*/
183: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
184: {
185: PetscErrorCode (*r)(DMPlexTransform);
186: PetscBool match;
188: PetscFunctionBegin;
190: PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match));
191: if (match) PetscFunctionReturn(PETSC_SUCCESS);
193: PetscCall(DMPlexTransformRegisterAll());
194: PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r));
195: PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);
197: PetscTryTypeMethod(tr, destroy);
198: PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops)));
199: PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method));
200: PetscCall((*r)(tr));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*@C
205: DMPlexTransformGetType - Gets the type name (as a string) from the transform.
207: Not Collective
209: Input Parameter:
210: . tr - The `DMPlexTransform`
212: Output Parameter:
213: . type - The `DMPlexTransformType` name
215: Level: intermediate
217: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
218: @*/
219: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
220: {
221: PetscFunctionBegin;
223: PetscAssertPointer(type, 2);
224: PetscCall(DMPlexTransformRegisterAll());
225: *type = ((PetscObject)tr)->type_name;
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
230: {
231: PetscViewerFormat format;
233: PetscFunctionBegin;
234: PetscCall(PetscViewerGetFormat(v, &format));
235: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
236: const PetscInt *trTypes = NULL;
237: IS trIS;
238: PetscInt cols = 8;
239: PetscInt Nrt = 8, f, g;
241: if (tr->trType) PetscCall(DMLabelView(tr->trType, v));
242: PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n"));
243: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
244: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
245: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]));
246: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
247: PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n"));
248: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
249: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
250: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]));
251: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
253: if (tr->trType) {
254: PetscCall(DMLabelGetNumValues(tr->trType, &Nrt));
255: PetscCall(DMLabelGetValueIS(tr->trType, &trIS));
256: PetscCall(ISGetIndices(trIS, &trTypes));
257: }
258: PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n"));
259: PetscCall(PetscViewerASCIIPrintf(v, " "));
260: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
261: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
262: for (f = 0; f < Nrt; ++f) {
263: PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT " |", trTypes ? trTypes[f] : f));
264: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]));
265: PetscCall(PetscViewerASCIIPrintf(v, " |\n"));
266: }
267: if (tr->trType) {
268: PetscCall(ISRestoreIndices(trIS, &trTypes));
269: PetscCall(ISDestroy(&trIS));
270: }
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*@C
276: DMPlexTransformView - Views a `DMPlexTransform`
278: Collective
280: Input Parameters:
281: + tr - the `DMPlexTransform` object to view
282: - v - the viewer
284: Level: beginner
286: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
287: @*/
288: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
289: {
290: PetscBool isascii;
292: PetscFunctionBegin;
294: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v));
296: PetscCheckSameComm(tr, 1, v, 2);
297: PetscCall(PetscViewerCheckWritable(v));
298: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v));
299: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
300: if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v));
301: PetscTryTypeMethod(tr, view, v);
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*@
306: DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database
308: Collective
310: Input Parameter:
311: . tr - the `DMPlexTransform` object to set options for
313: Options Database Keys:
314: + -dm_plex_transform_type - Set the transform type, e.g. refine_regular
315: . -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point
316: - -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
318: Level: intermediate
320: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
321: @*/
322: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
323: {
324: char typeName[1024];
325: const char *defName = DMPLEXREFINEREGULAR;
326: PetscBool flg;
328: PetscFunctionBegin;
330: PetscObjectOptionsBegin((PetscObject)tr);
331: PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg));
332: if (flg) PetscCall(DMPlexTransformSetType(tr, typeName));
333: else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName));
334: PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &tr->labelMatchStrata, NULL));
335: 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));
336: PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
337: /* process any options handlers added with PetscObjectAddOptionsHandler() */
338: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject));
339: PetscOptionsEnd();
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*@C
344: DMPlexTransformDestroy - Destroys a `DMPlexTransform`
346: Collective
348: Input Parameter:
349: . tr - the transform object to destroy
351: Level: beginner
353: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
354: @*/
355: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
356: {
357: PetscInt c;
359: PetscFunctionBegin;
360: if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
362: if (--((PetscObject)*tr)->refct > 0) {
363: *tr = NULL;
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: PetscTryTypeMethod(*tr, destroy);
368: PetscCall(DMDestroy(&(*tr)->dm));
369: PetscCall(DMLabelDestroy(&(*tr)->active));
370: PetscCall(DMLabelDestroy(&(*tr)->trType));
371: PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld));
372: PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew));
373: PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew));
374: PetscCall(PetscFree((*tr)->offset));
375: PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd));
376: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
377: PetscCall(PetscFEDestroy(&(*tr)->coordFE[c]));
378: PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c]));
379: }
380: if ((*tr)->trVerts) {
381: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
382: DMPolytopeType *rct;
383: PetscInt *rsize, *rcone, *rornt, Nct, n, r;
385: if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) {
386: PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
387: for (n = 0; n < Nct; ++n) {
388: if (rct[n] == DM_POLYTOPE_POINT) continue;
389: for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r]));
390: PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]]));
391: }
392: }
393: PetscCall(PetscFree((*tr)->trSubVerts[c]));
394: PetscCall(PetscFree((*tr)->trVerts[c]));
395: }
396: }
397: PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts));
398: PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom));
399: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
400: PetscCall(PetscHeaderDestroy(tr));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
405: {
406: DMLabel trType = tr->trType;
407: PetscInt c, cN, *off;
409: PetscFunctionBegin;
410: if (trType) {
411: DM dm;
412: IS rtIS;
413: const PetscInt *reftypes;
414: PetscInt Nrt, r;
416: PetscCall(DMPlexTransformGetDM(tr, &dm));
417: PetscCall(DMLabelGetNumValues(trType, &Nrt));
418: PetscCall(DMLabelGetValueIS(trType, &rtIS));
419: PetscCall(ISGetIndices(rtIS, &reftypes));
420: PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off));
421: for (r = 0; r < Nrt; ++r) {
422: const PetscInt rt = reftypes[r];
423: IS rtIS;
424: const PetscInt *points;
425: DMPolytopeType ct;
426: PetscInt np, p;
428: PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
429: PetscCall(ISGetLocalSize(rtIS, &np));
430: PetscCall(ISGetIndices(rtIS, &points));
431: if (!np) continue;
432: p = points[0];
433: PetscCall(ISRestoreIndices(rtIS, &points));
434: PetscCall(ISDestroy(&rtIS));
435: PetscCall(DMPlexGetCellType(dm, p, &ct));
436: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
437: const DMPolytopeType ctNew = (DMPolytopeType)cN;
438: DMPolytopeType *rct;
439: PetscInt *rsize, *cone, *ornt;
440: PetscInt Nct, n, s;
442: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
443: off[r * DM_NUM_POLYTOPES + ctNew] = -1;
444: break;
445: }
446: off[r * DM_NUM_POLYTOPES + ctNew] = 0;
447: for (s = 0; s <= r; ++s) {
448: const PetscInt st = reftypes[s];
449: DMPolytopeType sct;
450: PetscInt q, qrt;
452: PetscCall(DMLabelGetStratumIS(trType, st, &rtIS));
453: PetscCall(ISGetLocalSize(rtIS, &np));
454: PetscCall(ISGetIndices(rtIS, &points));
455: if (!np) continue;
456: q = points[0];
457: PetscCall(ISRestoreIndices(rtIS, &points));
458: PetscCall(ISDestroy(&rtIS));
459: PetscCall(DMPlexGetCellType(dm, q, &sct));
460: PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt));
461: 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);
462: if (st == rt) {
463: for (n = 0; n < Nct; ++n)
464: if (rct[n] == ctNew) break;
465: if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
466: break;
467: }
468: for (n = 0; n < Nct; ++n) {
469: if (rct[n] == ctNew) {
470: PetscInt sn;
472: PetscCall(DMLabelGetStratumSize(trType, st, &sn));
473: off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
474: }
475: }
476: }
477: }
478: }
479: PetscCall(ISRestoreIndices(rtIS, &reftypes));
480: PetscCall(ISDestroy(&rtIS));
481: } else {
482: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off));
483: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
484: const DMPolytopeType ct = (DMPolytopeType)c;
485: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
486: const DMPolytopeType ctNew = (DMPolytopeType)cN;
487: DMPolytopeType *rct;
488: PetscInt *rsize, *cone, *ornt;
489: PetscInt Nct, n, i;
491: 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) {
492: off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
493: continue;
494: }
495: off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
496: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
497: const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i];
498: const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];
500: PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt));
501: if (ict == ct) {
502: for (n = 0; n < Nct; ++n)
503: if (rct[n] == ctNew) break;
504: if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
505: break;
506: }
507: for (n = 0; n < Nct; ++n)
508: if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
509: }
510: }
511: }
512: }
513: *offset = off;
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
518: {
519: DM dm;
520: DMPolytopeType ctCell;
521: PetscInt pStart, pEnd, p, c, celldim = 0;
523: PetscFunctionBegin;
525: if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
526: PetscTryTypeMethod(tr, setup);
527: PetscCall(DMPlexTransformGetDM(tr, &dm));
528: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
529: if (pEnd > pStart) {
530: PetscCall(DMPlexGetCellType(dm, 0, &ctCell));
531: } else {
532: PetscInt dim;
534: PetscCall(DMGetDimension(dm, &dim));
535: switch (dim) {
536: case 0:
537: ctCell = DM_POLYTOPE_POINT;
538: break;
539: case 1:
540: ctCell = DM_POLYTOPE_SEGMENT;
541: break;
542: case 2:
543: ctCell = DM_POLYTOPE_TRIANGLE;
544: break;
545: case 3:
546: ctCell = DM_POLYTOPE_TETRAHEDRON;
547: break;
548: default:
549: ctCell = DM_POLYTOPE_UNKNOWN;
550: }
551: }
552: PetscCall(DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
553: for (p = pStart; p < pEnd; ++p) {
554: DMPolytopeType ct;
555: DMPolytopeType *rct;
556: PetscInt *rsize, *cone, *ornt;
557: PetscInt Nct, n;
559: PetscCall(DMPlexGetCellType(dm, p, &ct));
560: 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);
561: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
562: for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
563: }
564: PetscCall(DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
565: /* Construct sizes and offsets for each cell type */
566: if (!tr->ctStart) {
567: PetscInt *ctS, *ctSN, *ctC, *ctCN;
569: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
570: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
571: for (p = pStart; p < pEnd; ++p) {
572: DMPolytopeType ct;
573: DMPolytopeType *rct;
574: PetscInt *rsize, *cone, *ornt;
575: PetscInt Nct, n;
577: PetscCall(DMPlexGetCellType(dm, p, &ct));
578: 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);
579: ++ctC[ct];
580: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
581: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
582: }
583: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
584: const PetscInt cto = tr->ctOrderOld[c];
585: const PetscInt cton = tr->ctOrderOld[c + 1];
586: const PetscInt ctn = tr->ctOrderNew[c];
587: const PetscInt ctnn = tr->ctOrderNew[c + 1];
589: ctS[cton] = ctS[cto] + ctC[cto];
590: ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
591: }
592: PetscCall(PetscFree2(ctC, ctCN));
593: tr->ctStart = ctS;
594: tr->ctStartNew = ctSN;
595: }
596: PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
597: // Compute depth information
598: tr->depth = -1;
599: for (c = 0; c < DM_NUM_POLYTOPES; ++c)
600: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
601: PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
602: for (PetscInt d = 0; d <= tr->depth; ++d) {
603: tr->depthStart[d] = PETSC_MAX_INT;
604: tr->depthEnd[d] = -1;
605: }
606: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
607: const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);
609: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
610: tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
611: tr->depthEnd[dep] = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
612: }
613: tr->setupcalled = PETSC_TRUE;
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
618: {
619: PetscFunctionBegin;
621: PetscAssertPointer(dm, 2);
622: *dm = tr->dm;
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
627: {
628: PetscFunctionBegin;
631: PetscCall(PetscObjectReference((PetscObject)dm));
632: PetscCall(DMDestroy(&tr->dm));
633: tr->dm = dm;
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
638: {
639: PetscFunctionBegin;
641: PetscAssertPointer(active, 2);
642: *active = tr->active;
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
647: {
648: PetscFunctionBegin;
651: PetscCall(PetscObjectReference((PetscObject)active));
652: PetscCall(DMLabelDestroy(&tr->active));
653: tr->active = active;
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
658: {
659: PetscFunctionBegin;
660: if (!tr->coordFE[ct]) {
661: PetscInt dim, cdim;
663: dim = DMPolytopeTypeGetDim(ct);
664: PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
665: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
666: {
667: PetscDualSpace dsp;
668: PetscQuadrature quad;
669: DM K;
670: PetscFEGeom *cg;
671: PetscScalar *Xq;
672: PetscReal *xq, *wq;
673: PetscInt Nq, q;
675: PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
676: PetscCall(PetscMalloc1(Nq * cdim, &xq));
677: for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
678: PetscCall(PetscMalloc1(Nq, &wq));
679: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
680: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
681: PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
682: PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));
684: PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
685: PetscCall(PetscDualSpaceGetDM(dsp, &K));
686: PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]));
687: cg = tr->refGeom[ct];
688: PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
689: PetscCall(PetscQuadratureDestroy(&quad));
690: }
691: }
692: *fe = tr->coordFE[ct];
693: PetscFunctionReturn(PETSC_SUCCESS);
694: }
696: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
697: {
698: PetscInt dim, cdim;
700: PetscFunctionBegin;
701: PetscCall(DMGetDimension(dm, &dim));
702: PetscCall(DMSetDimension(tdm, dim));
703: PetscCall(DMGetCoordinateDim(dm, &cdim));
704: PetscCall(DMSetCoordinateDim(tdm, cdim));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: /*@
709: DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`
711: Input Parameters:
712: + tr - The `DMPlexTransform` object
713: - dm - The original `DM`
715: Output Parameter:
716: . tdm - The transformed `DM`
718: Level: advanced
720: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
721: @*/
722: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
723: {
724: PetscFunctionBegin;
725: PetscUseTypeMethod(tr, setdimensions, dm, tdm);
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
730: {
731: PetscFunctionBegin;
732: if (pStart) *pStart = 0;
733: if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
734: PetscFunctionReturn(PETSC_SUCCESS);
735: }
737: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
738: {
739: PetscInt ctNew;
741: PetscFunctionBegin;
743: PetscAssertPointer(celltype, 3);
744: /* TODO Can do bisection since everything is sorted */
745: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
746: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
748: if (cell >= ctSN && cell < ctEN) break;
749: }
750: PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
751: *celltype = (DMPolytopeType)ctNew;
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
755: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
756: {
757: PetscFunctionBegin;
759: if (start) *start = tr->ctStartNew[celltype];
760: if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
761: PetscFunctionReturn(PETSC_SUCCESS);
762: }
764: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
765: {
766: PetscFunctionBegin;
768: *depth = tr->depth;
769: PetscFunctionReturn(PETSC_SUCCESS);
770: }
772: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
773: {
774: PetscFunctionBegin;
776: if (start) *start = tr->depthStart[depth];
777: if (end) *end = tr->depthEnd[depth];
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: /*@
782: DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.
784: Not Collective
786: Input Parameters:
787: + tr - The `DMPlexTransform`
788: . ct - The type of the original point which produces the new point
789: . ctNew - The type of the new point
790: . p - The original point which produces the new point
791: - r - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`
793: Output Parameter:
794: . pNew - The new point number
796: Level: developer
798: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
799: @*/
800: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
801: {
802: DMPolytopeType *rct;
803: PetscInt *rsize, *cone, *ornt;
804: PetscInt rt, Nct, n, off, rp;
805: DMLabel trType = tr->trType;
806: PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
807: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
808: PetscInt newp = ctSN, cind;
810: PetscFunctionBeginHot;
811: 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);
812: PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
813: if (trType) {
814: PetscCall(DMLabelGetValueIndex(trType, rt, &cind));
815: PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
816: 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);
817: } else {
818: cind = ct;
819: rp = p - ctS;
820: }
821: off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
822: 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);
823: newp += off;
824: for (n = 0; n < Nct; ++n) {
825: if (rct[n] == ctNew) {
826: if (rsize[n] && r >= rsize[n])
827: 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]);
828: newp += rp * rsize[n] + r;
829: break;
830: }
831: }
833: 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);
834: *pNew = newp;
835: PetscFunctionReturn(PETSC_SUCCESS);
836: }
838: /*@
839: DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.
841: Not Collective
843: Input Parameters:
844: + tr - The `DMPlexTransform`
845: - pNew - The new point number
847: Output Parameters:
848: + ct - The type of the original point which produces the new point
849: . ctNew - The type of the new point
850: . p - The original point which produces the new point
851: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
853: Level: developer
855: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
856: @*/
857: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
858: {
859: DMLabel trType = tr->trType;
860: DMPolytopeType *rct, ctN;
861: PetscInt *rsize, *cone, *ornt;
862: PetscInt rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
863: PetscInt offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;
865: PetscFunctionBegin;
866: PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
867: if (trType) {
868: DM dm;
869: IS rtIS;
870: const PetscInt *reftypes;
871: PetscInt Nrt, r, rtStart;
873: PetscCall(DMPlexTransformGetDM(tr, &dm));
874: PetscCall(DMLabelGetNumValues(trType, &Nrt));
875: PetscCall(DMLabelGetValueIS(trType, &rtIS));
876: PetscCall(ISGetIndices(rtIS, &reftypes));
877: for (r = 0; r < Nrt; ++r) {
878: const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];
880: if (tr->ctStartNew[ctN] + off > pNew) continue;
881: /* Check that any of this refinement type exist */
882: /* TODO Actually keep track of the number produced here instead */
883: if (off > offset) {
884: rt = reftypes[r];
885: offset = off;
886: }
887: }
888: PetscCall(ISRestoreIndices(rtIS, &reftypes));
889: PetscCall(ISDestroy(&rtIS));
890: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
891: /* TODO Map refinement types to cell types */
892: PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL));
893: PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
894: for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
895: PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
897: if ((rtStart >= ctS) && (rtStart < ctE)) break;
898: }
899: PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt);
900: } else {
901: for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
902: const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];
904: if (tr->ctStartNew[ctN] + off > pNew) continue;
905: if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
906: /* TODO Actually keep track of the number produced here instead */
907: if (off > offset) {
908: ctO = ctTmp;
909: offset = off;
910: }
911: }
912: rt = -1;
913: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
914: }
915: ctS = tr->ctStart[ctO];
916: ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
917: if (trType) {
918: for (rtS = ctS; rtS < ctE; ++rtS) {
919: PetscInt val;
920: PetscCall(DMLabelGetValue(trType, rtS, &val));
921: if (val == rt) break;
922: }
923: 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);
924: } else rtS = ctS;
925: PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
926: 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);
927: for (n = 0; n < Nct; ++n) {
928: if (rct[n] == ctN) {
929: PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c;
931: if (trType) {
932: for (c = ctS; c < ctE; ++c) {
933: PetscCall(DMLabelGetValue(trType, c, &val));
934: if (val == rt) {
935: if (tmp < rsize[n]) break;
936: tmp -= rsize[n];
937: }
938: }
939: PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew);
940: rp = c - ctS;
941: rO = tmp;
942: } else {
943: // This assumes that all points of type ctO transform the same way
944: rp = tmp / rsize[n];
945: rO = tmp % rsize[n];
946: }
947: break;
948: }
949: }
950: PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
951: pO = rp + ctS;
952: 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);
953: if (ct) *ct = (DMPolytopeType)ctO;
954: if (ctNew) *ctNew = (DMPolytopeType)ctN;
955: if (p) *p = pO;
956: if (r) *r = rO;
957: PetscFunctionReturn(PETSC_SUCCESS);
958: }
960: /*@
961: DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.
963: Input Parameters:
964: + tr - The `DMPlexTransform` object
965: . source - The source cell type
966: - p - The source point, which can also determine the refine type
968: Output Parameters:
969: + rt - The refine type for this point
970: . Nt - The number of types produced by this point
971: . target - An array of length `Nt` giving the types produced
972: . size - An array of length `Nt` giving the number of cells of each type produced
973: . cone - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
974: - ornt - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point
976: Level: advanced
978: Notes:
979: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
980: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
981: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
982: .vb
983: the number of cones to be taken, or 0 for the current cell
984: the cell cone point number at each level from which it is subdivided
985: the replica number r of the subdivision.
986: .ve
987: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
988: .vb
989: Nt = 2
990: target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
991: size = {1, 2}
992: cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
993: ornt = { 0, 0, 0, 0}
994: .ve
996: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
997: @*/
998: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
999: {
1000: PetscFunctionBegin;
1001: PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
1002: PetscFunctionReturn(PETSC_SUCCESS);
1003: }
1005: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1006: {
1007: PetscFunctionBegin;
1008: *rnew = r;
1009: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: /* Returns the same thing */
1014: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1015: {
1016: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
1017: static PetscInt vertexS[] = {1};
1018: static PetscInt vertexC[] = {0};
1019: static PetscInt vertexO[] = {0};
1020: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
1021: static PetscInt edgeS[] = {1};
1022: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1023: static PetscInt edgeO[] = {0, 0};
1024: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
1025: static PetscInt tedgeS[] = {1};
1026: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1027: static PetscInt tedgeO[] = {0, 0};
1028: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
1029: static PetscInt triS[] = {1};
1030: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1031: static PetscInt triO[] = {0, 0, 0};
1032: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
1033: static PetscInt quadS[] = {1};
1034: 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};
1035: static PetscInt quadO[] = {0, 0, 0, 0};
1036: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
1037: static PetscInt tquadS[] = {1};
1038: 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};
1039: static PetscInt tquadO[] = {0, 0, 0, 0};
1040: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
1041: static PetscInt tetS[] = {1};
1042: 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};
1043: static PetscInt tetO[] = {0, 0, 0, 0};
1044: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
1045: static PetscInt hexS[] = {1};
1046: 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};
1047: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
1048: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
1049: static PetscInt tripS[] = {1};
1050: 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};
1051: static PetscInt tripO[] = {0, 0, 0, 0, 0};
1052: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
1053: static PetscInt ttripS[] = {1};
1054: 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};
1055: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
1056: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
1057: static PetscInt tquadpS[] = {1};
1058: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
1059: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1060: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
1061: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
1062: static PetscInt pyrS[] = {1};
1063: 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};
1064: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
1066: PetscFunctionBegin;
1067: if (rt) *rt = 0;
1068: switch (source) {
1069: case DM_POLYTOPE_POINT:
1070: *Nt = 1;
1071: *target = vertexT;
1072: *size = vertexS;
1073: *cone = vertexC;
1074: *ornt = vertexO;
1075: break;
1076: case DM_POLYTOPE_SEGMENT:
1077: *Nt = 1;
1078: *target = edgeT;
1079: *size = edgeS;
1080: *cone = edgeC;
1081: *ornt = edgeO;
1082: break;
1083: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1084: *Nt = 1;
1085: *target = tedgeT;
1086: *size = tedgeS;
1087: *cone = tedgeC;
1088: *ornt = tedgeO;
1089: break;
1090: case DM_POLYTOPE_TRIANGLE:
1091: *Nt = 1;
1092: *target = triT;
1093: *size = triS;
1094: *cone = triC;
1095: *ornt = triO;
1096: break;
1097: case DM_POLYTOPE_QUADRILATERAL:
1098: *Nt = 1;
1099: *target = quadT;
1100: *size = quadS;
1101: *cone = quadC;
1102: *ornt = quadO;
1103: break;
1104: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1105: *Nt = 1;
1106: *target = tquadT;
1107: *size = tquadS;
1108: *cone = tquadC;
1109: *ornt = tquadO;
1110: break;
1111: case DM_POLYTOPE_TETRAHEDRON:
1112: *Nt = 1;
1113: *target = tetT;
1114: *size = tetS;
1115: *cone = tetC;
1116: *ornt = tetO;
1117: break;
1118: case DM_POLYTOPE_HEXAHEDRON:
1119: *Nt = 1;
1120: *target = hexT;
1121: *size = hexS;
1122: *cone = hexC;
1123: *ornt = hexO;
1124: break;
1125: case DM_POLYTOPE_TRI_PRISM:
1126: *Nt = 1;
1127: *target = tripT;
1128: *size = tripS;
1129: *cone = tripC;
1130: *ornt = tripO;
1131: break;
1132: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1133: *Nt = 1;
1134: *target = ttripT;
1135: *size = ttripS;
1136: *cone = ttripC;
1137: *ornt = ttripO;
1138: break;
1139: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1140: *Nt = 1;
1141: *target = tquadpT;
1142: *size = tquadpS;
1143: *cone = tquadpC;
1144: *ornt = tquadpO;
1145: break;
1146: case DM_POLYTOPE_PYRAMID:
1147: *Nt = 1;
1148: *target = pyrT;
1149: *size = pyrS;
1150: *cone = pyrC;
1151: *ornt = pyrO;
1152: break;
1153: default:
1154: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1155: }
1156: PetscFunctionReturn(PETSC_SUCCESS);
1157: }
1159: /*@
1160: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
1162: Not Collective
1164: Input Parameters:
1165: + tr - The `DMPlexTransform`
1166: . sct - The source point cell type, from whom the new cell is being produced
1167: . sp - The source point
1168: . so - The orientation of the source point in its enclosing parent
1169: . tct - The target point cell type
1170: . r - The replica number requested for the produced cell type
1171: - o - The orientation of the replica
1173: Output Parameters:
1174: + rnew - The replica number, given the orientation of the parent
1175: - onew - The replica orientation, given the orientation of the parent
1177: Level: advanced
1179: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1180: @*/
1181: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1182: {
1183: PetscFunctionBeginHot;
1184: PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1189: {
1190: DM dm;
1191: PetscInt pStart, pEnd, p, pNew;
1193: PetscFunctionBegin;
1194: PetscCall(DMPlexTransformGetDM(tr, &dm));
1195: /* Must create the celltype label here so that we do not automatically try to compute the types */
1196: PetscCall(DMCreateLabel(rdm, "celltype"));
1197: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1198: for (p = pStart; p < pEnd; ++p) {
1199: DMPolytopeType ct;
1200: DMPolytopeType *rct;
1201: PetscInt *rsize, *rcone, *rornt;
1202: PetscInt Nct, n, r;
1204: PetscCall(DMPlexGetCellType(dm, p, &ct));
1205: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1206: for (n = 0; n < Nct; ++n) {
1207: for (r = 0; r < rsize[n]; ++r) {
1208: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1209: PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
1210: PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
1211: }
1212: }
1213: }
1214: /* Let the DM know we have set all the cell types */
1215: {
1216: DMLabel ctLabel;
1217: DM_Plex *plex = (DM_Plex *)rdm->data;
1219: PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
1220: PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
1221: }
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1226: {
1227: DMPolytopeType ctNew;
1229: PetscFunctionBegin;
1231: PetscAssertPointer(coneSize, 3);
1232: PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
1233: *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType)ctNew);
1234: PetscFunctionReturn(PETSC_SUCCESS);
1235: }
1237: /* The orientation o is for the interior of the cell p */
1238: 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[])
1239: {
1240: DM dm;
1241: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1242: const PetscInt *cone;
1243: PetscInt c, coff = *coneoff, ooff = *orntoff;
1245: PetscFunctionBegin;
1246: PetscCall(DMPlexTransformGetDM(tr, &dm));
1247: PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
1248: for (c = 0; c < csizeNew; ++c) {
1249: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1250: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1251: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1252: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1253: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1254: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1255: const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1256: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1257: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1258: PetscInt lc;
1260: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1261: for (lc = 0; lc < fn; ++lc) {
1262: const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po);
1263: const PetscInt acp = rcone[coff++];
1264: const PetscInt pcp = parr[acp * 2];
1265: const PetscInt pco = parr[acp * 2 + 1];
1266: const PetscInt *ppornt;
1268: ppp = pp;
1269: pp = pcone[pcp];
1270: PetscCall(DMPlexGetCellType(dm, pp, &pct));
1271: // Restore the parent cone from the last iterate
1272: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
1273: PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
1274: PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
1275: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1276: PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
1277: }
1278: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
1279: pr = rcone[coff++];
1280: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1281: PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
1282: PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
1283: orntNew[c] = fo;
1284: }
1285: PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
1286: *coneoff = coff;
1287: *orntoff = ooff;
1288: PetscFunctionReturn(PETSC_SUCCESS);
1289: }
1291: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1292: {
1293: DM dm;
1294: DMPolytopeType ct;
1295: PetscInt *coneNew, *orntNew;
1296: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1298: PetscFunctionBegin;
1299: PetscCall(DMPlexTransformGetDM(tr, &dm));
1300: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1301: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1302: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1303: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1304: for (p = pStart; p < pEnd; ++p) {
1305: PetscInt coff, ooff;
1306: DMPolytopeType *rct;
1307: PetscInt *rsize, *rcone, *rornt;
1308: PetscInt Nct, n, r;
1310: PetscCall(DMPlexGetCellType(dm, p, &ct));
1311: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1312: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1313: const DMPolytopeType ctNew = rct[n];
1315: for (r = 0; r < rsize[n]; ++r) {
1316: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1317: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
1318: PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
1319: PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
1320: }
1321: }
1322: }
1323: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1324: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1325: PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
1326: PetscCall(DMPlexSymmetrize(rdm));
1327: PetscCall(DMPlexStratify(rdm));
1328: PetscFunctionReturn(PETSC_SUCCESS);
1329: }
1331: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1332: {
1333: DM dm;
1334: DMPolytopeType ct, qct;
1335: DMPolytopeType *rct;
1336: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1337: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1339: PetscFunctionBegin;
1341: PetscAssertPointer(cone, 4);
1342: PetscAssertPointer(ornt, 5);
1343: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1344: PetscCall(DMPlexTransformGetDM(tr, &dm));
1345: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1346: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1347: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1348: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1349: for (n = 0; n < Nct; ++n) {
1350: const DMPolytopeType ctNew = rct[n];
1351: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1352: PetscInt Nr = rsize[n], fn, c;
1354: if (ctNew == qct) Nr = r;
1355: for (nr = 0; nr < Nr; ++nr) {
1356: for (c = 0; c < csizeNew; ++c) {
1357: ++coff; /* Cell type of new cone point */
1358: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1359: coff += fn;
1360: ++coff; /* Replica number of new cone point */
1361: ++ooff; /* Orientation of new cone point */
1362: }
1363: }
1364: if (ctNew == qct) break;
1365: }
1366: PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1367: *cone = qcone;
1368: *ornt = qornt;
1369: PetscFunctionReturn(PETSC_SUCCESS);
1370: }
1372: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1373: {
1374: DM dm;
1375: DMPolytopeType ct, qct;
1376: DMPolytopeType *rct;
1377: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1378: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1380: PetscFunctionBegin;
1382: if (cone) PetscAssertPointer(cone, 3);
1383: if (ornt) PetscAssertPointer(ornt, 4);
1384: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1385: PetscCall(DMPlexTransformGetDM(tr, &dm));
1386: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1387: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1388: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1389: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1390: for (n = 0; n < Nct; ++n) {
1391: const DMPolytopeType ctNew = rct[n];
1392: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1393: PetscInt Nr = rsize[n], fn, c;
1395: if (ctNew == qct) Nr = r;
1396: for (nr = 0; nr < Nr; ++nr) {
1397: for (c = 0; c < csizeNew; ++c) {
1398: ++coff; /* Cell type of new cone point */
1399: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1400: coff += fn;
1401: ++coff; /* Replica number of new cone point */
1402: ++ooff; /* Orientation of new cone point */
1403: }
1404: }
1405: if (ctNew == qct) break;
1406: }
1407: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1408: if (cone) *cone = qcone;
1409: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1410: if (ornt) *ornt = qornt;
1411: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1412: PetscFunctionReturn(PETSC_SUCCESS);
1413: }
1415: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1416: {
1417: DM dm;
1419: PetscFunctionBegin;
1421: PetscCall(DMPlexTransformGetDM(tr, &dm));
1422: if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
1423: if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
1424: PetscFunctionReturn(PETSC_SUCCESS);
1425: }
1427: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1428: {
1429: PetscInt ict;
1431: PetscFunctionBegin;
1432: PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
1433: for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1434: const DMPolytopeType ct = (DMPolytopeType)ict;
1435: DMPlexTransform reftr;
1436: DM refdm, trdm;
1437: Vec coordinates;
1438: const PetscScalar *coords;
1439: DMPolytopeType *rct;
1440: PetscInt *rsize, *rcone, *rornt;
1441: PetscInt Nct, n, r, pNew;
1442: PetscInt trdim, vStart, vEnd, Nc;
1443: const PetscInt debug = 0;
1444: const char *typeName;
1446: /* Since points are 0-dimensional, coordinates make no sense */
1447: if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue;
1448: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
1449: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
1450: PetscCall(DMPlexTransformSetDM(reftr, refdm));
1451: PetscCall(DMPlexTransformGetType(tr, &typeName));
1452: PetscCall(DMPlexTransformSetType(reftr, typeName));
1453: PetscCall(DMPlexTransformSetUp(reftr));
1454: PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));
1456: PetscCall(DMGetDimension(trdm, &trdim));
1457: PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
1458: tr->trNv[ct] = vEnd - vStart;
1459: PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
1460: PetscCall(VecGetLocalSize(coordinates, &Nc));
1461: 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);
1462: PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
1463: PetscCall(VecGetArrayRead(coordinates, &coords));
1464: PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
1465: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1467: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
1468: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1469: for (n = 0; n < Nct; ++n) {
1470: /* Since points are 0-dimensional, coordinates make no sense */
1471: if (rct[n] == DM_POLYTOPE_POINT) continue;
1472: PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
1473: for (r = 0; r < rsize[n]; ++r) {
1474: PetscInt *closure = NULL;
1475: PetscInt clSize, cl, Nv = 0;
1477: PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
1478: PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
1479: PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1480: for (cl = 0; cl < clSize * 2; cl += 2) {
1481: const PetscInt sv = closure[cl];
1483: if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1484: }
1485: PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1486: 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]);
1487: }
1488: }
1489: if (debug) {
1490: DMPolytopeType *rct;
1491: PetscInt *rsize, *rcone, *rornt;
1492: PetscInt v, dE = trdim, d, off = 0;
1494: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
1495: for (v = 0; v < tr->trNv[ct]; ++v) {
1496: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1497: for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
1498: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1499: }
1501: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1502: for (n = 0; n < Nct; ++n) {
1503: if (rct[n] == DM_POLYTOPE_POINT) continue;
1504: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
1505: for (r = 0; r < rsize[n]; ++r) {
1506: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1507: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
1508: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1509: }
1510: }
1511: }
1512: PetscCall(DMDestroy(&refdm));
1513: PetscCall(DMDestroy(&trdm));
1514: PetscCall(DMPlexTransformDestroy(&reftr));
1515: }
1516: PetscFunctionReturn(PETSC_SUCCESS);
1517: }
1519: /*@C
1520: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
1522: Input Parameters:
1523: + tr - The `DMPlexTransform` object
1524: - ct - The cell type
1526: Output Parameters:
1527: + Nv - The number of transformed vertices in the closure of the reference cell of given type
1528: - trVerts - The coordinates of these vertices in the reference cell
1530: Level: developer
1532: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
1533: @*/
1534: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1535: {
1536: PetscFunctionBegin;
1537: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1538: if (Nv) *Nv = tr->trNv[ct];
1539: if (trVerts) *trVerts = tr->trVerts[ct];
1540: PetscFunctionReturn(PETSC_SUCCESS);
1541: }
1543: /*@C
1544: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
1546: Input Parameters:
1547: + tr - The `DMPlexTransform` object
1548: . ct - The cell type
1549: . rct - The subcell type
1550: - r - The subcell index
1552: Output Parameter:
1553: . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()`
1555: Level: developer
1557: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
1558: @*/
1559: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1560: {
1561: PetscFunctionBegin;
1562: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1563: PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
1564: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1565: PetscFunctionReturn(PETSC_SUCCESS);
1566: }
1568: /* Computes new vertex as the barycenter, or centroid */
1569: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1570: {
1571: PetscInt v, d;
1573: PetscFunctionBeginHot;
1574: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
1575: for (d = 0; d < dE; ++d) out[d] = 0.0;
1576: for (v = 0; v < Nv; ++v)
1577: for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1578: for (d = 0; d < dE; ++d) out[d] /= Nv;
1579: PetscFunctionReturn(PETSC_SUCCESS);
1580: }
1582: /*@
1583: DMPlexTransformMapCoordinates - Calculate new coordinates for produced points
1585: Not collective
1587: Input Parameters:
1588: + tr - The `DMPlexTransform`
1589: . pct - The cell type of the parent, from whom the new cell is being produced
1590: . ct - The type being produced
1591: . p - The original point
1592: . r - The replica number requested for the produced cell type
1593: . Nv - Number of vertices in the closure of the parent cell
1594: . dE - Spatial dimension
1595: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
1597: Output Parameter:
1598: . out - The coordinates of the new vertices
1600: Level: intermediate
1602: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
1603: @*/
1604: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1605: {
1606: PetscFunctionBeginHot;
1607: if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1608: PetscFunctionReturn(PETSC_SUCCESS);
1609: }
1611: /*
1612: DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label
1614: Not Collective
1616: Input Parameters:
1617: + tr - The `DMPlexTransform`
1618: . label - The label in the transformed mesh
1619: . pp - The parent point in the original mesh
1620: . pct - The cell type of the parent point
1621: . p - The point in the transformed mesh
1622: . ct - The cell type of the point
1623: . r - The replica number of the point
1624: - val - The label value of the parent point
1626: Level: developer
1628: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
1629: */
1630: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
1631: {
1632: PetscFunctionBeginHot;
1633: if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
1634: PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
1635: PetscFunctionReturn(PETSC_SUCCESS);
1636: }
1638: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1639: {
1640: DM dm;
1641: IS valueIS;
1642: const PetscInt *values;
1643: PetscInt defVal, Nv, val;
1645: PetscFunctionBegin;
1646: PetscCall(DMPlexTransformGetDM(tr, &dm));
1647: PetscCall(DMLabelGetDefaultValue(label, &defVal));
1648: PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
1649: PetscCall(DMLabelGetValueIS(label, &valueIS));
1650: PetscCall(ISGetLocalSize(valueIS, &Nv));
1651: PetscCall(ISGetIndices(valueIS, &values));
1652: for (val = 0; val < Nv; ++val) {
1653: IS pointIS;
1654: const PetscInt *points;
1655: PetscInt numPoints, p;
1657: /* Ensure refined label is created with same number of strata as
1658: * original (even if no entries here). */
1659: PetscCall(DMLabelAddStratum(labelNew, values[val]));
1660: PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
1661: PetscCall(ISGetLocalSize(pointIS, &numPoints));
1662: PetscCall(ISGetIndices(pointIS, &points));
1663: for (p = 0; p < numPoints; ++p) {
1664: const PetscInt point = points[p];
1665: DMPolytopeType ct;
1666: DMPolytopeType *rct;
1667: PetscInt *rsize, *rcone, *rornt;
1668: PetscInt Nct, n, r, pNew = 0;
1670: PetscCall(DMPlexGetCellType(dm, point, &ct));
1671: PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1672: for (n = 0; n < Nct; ++n) {
1673: for (r = 0; r < rsize[n]; ++r) {
1674: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
1675: PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
1676: }
1677: }
1678: }
1679: PetscCall(ISRestoreIndices(pointIS, &points));
1680: PetscCall(ISDestroy(&pointIS));
1681: }
1682: PetscCall(ISRestoreIndices(valueIS, &values));
1683: PetscCall(ISDestroy(&valueIS));
1684: PetscFunctionReturn(PETSC_SUCCESS);
1685: }
1687: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1688: {
1689: DM dm;
1690: PetscInt numLabels, l;
1692: PetscFunctionBegin;
1693: PetscCall(DMPlexTransformGetDM(tr, &dm));
1694: PetscCall(DMGetNumLabels(dm, &numLabels));
1695: for (l = 0; l < numLabels; ++l) {
1696: DMLabel label, labelNew;
1697: const char *lname;
1698: PetscBool isDepth, isCellType;
1700: PetscCall(DMGetLabelName(dm, l, &lname));
1701: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
1702: if (isDepth) continue;
1703: PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
1704: if (isCellType) continue;
1705: PetscCall(DMCreateLabel(rdm, lname));
1706: PetscCall(DMGetLabel(dm, lname, &label));
1707: PetscCall(DMGetLabel(rdm, lname, &labelNew));
1708: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1709: }
1710: PetscFunctionReturn(PETSC_SUCCESS);
1711: }
1713: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1714: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1715: {
1716: DM dm;
1717: PetscInt Nf, f, Nds, s;
1719: PetscFunctionBegin;
1720: PetscCall(DMPlexTransformGetDM(tr, &dm));
1721: PetscCall(DMGetNumFields(dm, &Nf));
1722: for (f = 0; f < Nf; ++f) {
1723: DMLabel label, labelNew;
1724: PetscObject obj;
1725: const char *lname;
1727: PetscCall(DMGetField(rdm, f, &label, &obj));
1728: if (!label) continue;
1729: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1730: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1731: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1732: PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
1733: PetscCall(DMLabelDestroy(&labelNew));
1734: }
1735: PetscCall(DMGetNumDS(dm, &Nds));
1736: for (s = 0; s < Nds; ++s) {
1737: DMLabel label, labelNew;
1738: const char *lname;
1740: PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
1741: if (!label) continue;
1742: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1743: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1744: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1745: PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
1746: PetscCall(DMLabelDestroy(&labelNew));
1747: }
1748: PetscFunctionReturn(PETSC_SUCCESS);
1749: }
1751: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1752: {
1753: DM dm;
1754: PetscSF sf, sfNew;
1755: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
1756: const PetscInt *localPoints;
1757: const PetscSFNode *remotePoints;
1758: PetscInt *localPointsNew;
1759: PetscSFNode *remotePointsNew;
1760: PetscInt pStartNew, pEndNew, pNew;
1761: /* Brute force algorithm */
1762: PetscSF rsf;
1763: PetscSection s;
1764: const PetscInt *rootdegree;
1765: PetscInt *rootPointsNew, *remoteOffsets;
1766: PetscInt numPointsNew, pStart, pEnd, p;
1768: PetscFunctionBegin;
1769: PetscCall(DMPlexTransformGetDM(tr, &dm));
1770: PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
1771: PetscCall(DMGetPointSF(dm, &sf));
1772: PetscCall(DMGetPointSF(rdm, &sfNew));
1773: /* Calculate size of new SF */
1774: PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
1775: if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
1776: for (l = 0; l < numLeaves; ++l) {
1777: const PetscInt p = localPoints[l];
1778: DMPolytopeType ct;
1779: DMPolytopeType *rct;
1780: PetscInt *rsize, *rcone, *rornt;
1781: PetscInt Nct, n;
1783: PetscCall(DMPlexGetCellType(dm, p, &ct));
1784: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1785: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
1786: }
1787: /* Send new root point numbers
1788: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1789: */
1790: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1791: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
1792: PetscCall(PetscSectionSetChart(s, pStart, pEnd));
1793: for (p = pStart; p < pEnd; ++p) {
1794: DMPolytopeType ct;
1795: DMPolytopeType *rct;
1796: PetscInt *rsize, *rcone, *rornt;
1797: PetscInt Nct, n;
1799: PetscCall(DMPlexGetCellType(dm, p, &ct));
1800: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1801: for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
1802: }
1803: PetscCall(PetscSectionSetUp(s));
1804: PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
1805: PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
1806: PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
1807: PetscCall(PetscFree(remoteOffsets));
1808: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1809: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1810: PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
1811: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1812: for (p = pStart; p < pEnd; ++p) {
1813: DMPolytopeType ct;
1814: DMPolytopeType *rct;
1815: PetscInt *rsize, *rcone, *rornt;
1816: PetscInt Nct, n, r, off;
1818: if (!rootdegree[p - pStart]) continue;
1819: PetscCall(PetscSectionGetOffset(s, p, &off));
1820: PetscCall(DMPlexGetCellType(dm, p, &ct));
1821: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1822: for (n = 0, m = 0; n < Nct; ++n) {
1823: for (r = 0; r < rsize[n]; ++r, ++m) {
1824: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1825: rootPointsNew[off + m] = pNew;
1826: }
1827: }
1828: }
1829: PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1830: PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1831: PetscCall(PetscSFDestroy(&rsf));
1832: PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
1833: PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
1834: for (l = 0, m = 0; l < numLeaves; ++l) {
1835: const PetscInt p = localPoints[l];
1836: DMPolytopeType ct;
1837: DMPolytopeType *rct;
1838: PetscInt *rsize, *rcone, *rornt;
1839: PetscInt Nct, n, r, q, off;
1841: PetscCall(PetscSectionGetOffset(s, p, &off));
1842: PetscCall(DMPlexGetCellType(dm, p, &ct));
1843: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1844: for (n = 0, q = 0; n < Nct; ++n) {
1845: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1846: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1847: localPointsNew[m] = pNew;
1848: remotePointsNew[m].index = rootPointsNew[off + q];
1849: remotePointsNew[m].rank = remotePoints[l].rank;
1850: }
1851: }
1852: }
1853: PetscCall(PetscSectionDestroy(&s));
1854: PetscCall(PetscFree(rootPointsNew));
1855: /* SF needs sorted leaves to correctly calculate Gather */
1856: {
1857: PetscSFNode *rp, *rtmp;
1858: PetscInt *lp, *idx, *ltmp, i;
1860: PetscCall(PetscMalloc1(numLeavesNew, &idx));
1861: PetscCall(PetscMalloc1(numLeavesNew, &lp));
1862: PetscCall(PetscMalloc1(numLeavesNew, &rp));
1863: for (i = 0; i < numLeavesNew; ++i) {
1864: 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);
1865: idx[i] = i;
1866: }
1867: PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
1868: for (i = 0; i < numLeavesNew; ++i) {
1869: lp[i] = localPointsNew[idx[i]];
1870: rp[i] = remotePointsNew[idx[i]];
1871: }
1872: ltmp = localPointsNew;
1873: localPointsNew = lp;
1874: rtmp = remotePointsNew;
1875: remotePointsNew = rp;
1876: PetscCall(PetscFree(idx));
1877: PetscCall(PetscFree(ltmp));
1878: PetscCall(PetscFree(rtmp));
1879: }
1880: PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
1881: PetscFunctionReturn(PETSC_SUCCESS);
1882: }
1884: /*
1885: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`.
1887: Not Collective
1889: Input Parameters:
1890: + tr - The `DMPlexTransform`
1891: . ct - The type of the parent cell
1892: . rct - The type of the produced cell
1893: . r - The index of the produced cell
1894: - x - The localized coordinates for the parent cell
1896: Output Parameter:
1897: . xr - The localized coordinates for the produced cell
1899: Level: developer
1901: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
1902: */
1903: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1904: {
1905: PetscFE fe = NULL;
1906: PetscInt cdim, v, *subcellV;
1908: PetscFunctionBegin;
1909: PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
1910: PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
1911: PetscCall(PetscFEGetNumComponents(fe, &cdim));
1912: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
1913: PetscFunctionReturn(PETSC_SUCCESS);
1914: }
1916: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1917: {
1918: DM dm, cdm, cdmCell, cdmNew, cdmCellNew;
1919: PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
1920: Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
1921: const PetscScalar *coords;
1922: PetscScalar *coordsNew;
1923: const PetscReal *maxCell, *Lstart, *L;
1924: PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1925: PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;
1927: PetscFunctionBegin;
1928: PetscCall(DMPlexTransformGetDM(tr, &dm));
1929: PetscCall(DMGetCoordinateDM(dm, &cdm));
1930: PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
1931: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
1932: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
1933: if (localized) {
1934: /* Localize coordinates of new vertices */
1935: localizeVertices = PETSC_TRUE;
1936: /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
1937: if (!maxCell) localizeCells = PETSC_TRUE;
1938: }
1939: PetscCall(DMGetCoordinateSection(dm, &coordSection));
1940: PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
1941: if (maxCell) {
1942: PetscReal maxCellNew[3];
1944: for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
1945: PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L));
1946: }
1947: PetscCall(DMGetCoordinateDim(rdm, &dE));
1948: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
1949: PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
1950: PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
1951: PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
1952: PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
1953: /* Localization should be inherited */
1954: /* Stefano calculates parent cells for each new cell for localization */
1955: /* Localized cells need coordinates of closure */
1956: for (v = vStartNew; v < vEndNew; ++v) {
1957: PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
1958: PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
1959: }
1960: PetscCall(PetscSectionSetUp(coordSectionNew));
1961: PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));
1963: if (localizeCells) {
1964: PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
1965: PetscCall(DMClone(cdmNew, &cdmCellNew));
1966: PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
1967: PetscCall(DMDestroy(&cdmCellNew));
1969: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
1970: PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
1971: PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
1972: PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
1973: PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));
1975: PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
1976: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1977: for (c = cStart; c < cEnd; ++c) {
1978: PetscInt dof;
1980: PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
1981: if (dof) {
1982: DMPolytopeType ct;
1983: DMPolytopeType *rct;
1984: PetscInt *rsize, *rcone, *rornt;
1985: PetscInt dim, cNew, Nct, n, r;
1987: PetscCall(DMPlexGetCellType(dm, c, &ct));
1988: dim = DMPolytopeTypeGetDim(ct);
1989: PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1990: /* This allows for different cell types */
1991: for (n = 0; n < Nct; ++n) {
1992: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
1993: for (r = 0; r < rsize[n]; ++r) {
1994: PetscInt *closure = NULL;
1995: PetscInt clSize, cl, Nv = 0;
1997: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
1998: PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
1999: for (cl = 0; cl < clSize * 2; cl += 2) {
2000: if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
2001: }
2002: PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2003: PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
2004: PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
2005: }
2006: }
2007: }
2008: }
2009: PetscCall(PetscSectionSetUp(coordSectionCellNew));
2010: PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
2011: }
2012: PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
2013: {
2014: VecType vtype;
2015: PetscInt coordSizeNew, bs;
2016: const char *name;
2018: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
2019: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
2020: PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
2021: PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
2022: PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
2023: PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
2024: PetscCall(VecGetBlockSize(coordsLocal, &bs));
2025: PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
2026: PetscCall(VecGetType(coordsLocal, &vtype));
2027: PetscCall(VecSetType(coordsLocalNew, vtype));
2028: }
2029: PetscCall(VecGetArrayRead(coordsLocal, &coords));
2030: PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
2031: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2032: /* First set coordinates for vertices */
2033: for (p = pStart; p < pEnd; ++p) {
2034: DMPolytopeType ct;
2035: DMPolytopeType *rct;
2036: PetscInt *rsize, *rcone, *rornt;
2037: PetscInt Nct, n, r;
2038: PetscBool hasVertex = PETSC_FALSE;
2040: PetscCall(DMPlexGetCellType(dm, p, &ct));
2041: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2042: for (n = 0; n < Nct; ++n) {
2043: if (rct[n] == DM_POLYTOPE_POINT) {
2044: hasVertex = PETSC_TRUE;
2045: break;
2046: }
2047: }
2048: if (hasVertex) {
2049: const PetscScalar *icoords = NULL;
2050: const PetscScalar *array = NULL;
2051: PetscScalar *pcoords = NULL;
2052: PetscBool isDG;
2053: PetscInt Nc, Nv, v, d;
2055: PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2057: icoords = pcoords;
2058: Nv = Nc / dEo;
2059: if (ct != DM_POLYTOPE_POINT) {
2060: if (localizeVertices && maxCell) {
2061: PetscScalar anchor[3];
2063: for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
2064: for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
2065: }
2066: }
2067: for (n = 0; n < Nct; ++n) {
2068: if (rct[n] != DM_POLYTOPE_POINT) continue;
2069: for (r = 0; r < rsize[n]; ++r) {
2070: PetscScalar vcoords[3];
2071: PetscInt vNew, off;
2073: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
2074: PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
2075: PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
2076: PetscCall(DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
2077: }
2078: }
2079: PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2080: }
2081: }
2082: PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
2083: PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
2084: PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
2085: PetscCall(VecDestroy(&coordsLocalNew));
2086: PetscCall(PetscSectionDestroy(&coordSectionNew));
2087: /* Then set coordinates for cells by localizing */
2088: if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
2089: else {
2090: VecType vtype;
2091: PetscInt coordSizeNew, bs;
2092: const char *name;
2094: PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
2095: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
2096: PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
2097: PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
2098: PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
2099: PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
2100: PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
2101: PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
2102: PetscCall(VecGetType(coordsLocalCell, &vtype));
2103: PetscCall(VecSetType(coordsLocalCellNew, vtype));
2104: PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
2105: PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));
2107: for (p = pStart; p < pEnd; ++p) {
2108: DMPolytopeType ct;
2109: DMPolytopeType *rct;
2110: PetscInt *rsize, *rcone, *rornt;
2111: PetscInt dof = 0, Nct, n, r;
2113: PetscCall(DMPlexGetCellType(dm, p, &ct));
2114: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2115: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
2116: if (dof) {
2117: const PetscScalar *pcoords;
2119: PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords));
2120: for (n = 0; n < Nct; ++n) {
2121: const PetscInt Nr = rsize[n];
2123: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2124: for (r = 0; r < Nr; ++r) {
2125: PetscInt pNew, offNew;
2127: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2128: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2129: cell to the ones it produces. */
2130: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2131: PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
2132: PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
2133: }
2134: }
2135: }
2136: }
2137: PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
2138: PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
2139: PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
2140: PetscCall(VecDestroy(&coordsLocalCellNew));
2141: PetscCall(PetscSectionDestroy(&coordSectionCellNew));
2142: }
2143: PetscFunctionReturn(PETSC_SUCCESS);
2144: }
2146: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
2147: {
2148: DM rdm;
2149: DMPlexInterpolatedFlag interp;
2150: PetscInt pStart, pEnd;
2152: PetscFunctionBegin;
2155: PetscAssertPointer(tdm, 3);
2156: PetscCall(PetscLogEventBegin(DMPLEX_Transform, tr, dm, 0, 0));
2157: PetscCall(DMPlexTransformSetDM(tr, dm));
2159: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
2160: PetscCall(DMSetType(rdm, DMPLEX));
2161: PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
2162: /* Calculate number of new points of each depth */
2163: PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
2164: PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
2165: /* Step 1: Set chart */
2166: PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
2167: PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
2168: /* Step 2: Set cone/support sizes (automatically stratifies) */
2169: PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
2170: /* Step 3: Setup refined DM */
2171: PetscCall(DMSetUp(rdm));
2172: /* Step 4: Set cones and supports (automatically symmetrizes) */
2173: PetscCall(DMPlexTransformSetCones(tr, rdm));
2174: /* Step 5: Create pointSF */
2175: PetscCall(DMPlexTransformCreateSF(tr, rdm));
2176: /* Step 6: Create labels */
2177: PetscCall(DMPlexTransformCreateLabels(tr, rdm));
2178: /* Step 7: Set coordinates */
2179: PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
2180: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm));
2181: // If the original DM was configured from options, the transformed DM should be as well
2182: rdm->setfromoptionscalled = dm->setfromoptionscalled;
2183: PetscCall(PetscLogEventEnd(DMPLEX_Transform, tr, dm, 0, 0));
2184: *tdm = rdm;
2185: PetscFunctionReturn(PETSC_SUCCESS);
2186: }
2188: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2189: {
2190: DMPlexTransform tr;
2191: DM cdm, rcdm;
2192: const char *prefix;
2194: PetscFunctionBegin;
2195: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
2196: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2197: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
2198: PetscCall(DMPlexTransformSetDM(tr, dm));
2199: PetscCall(DMPlexTransformSetFromOptions(tr));
2200: PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
2201: PetscCall(DMPlexTransformSetUp(tr));
2202: PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
2203: PetscCall(DMPlexTransformApply(tr, dm, rdm));
2204: PetscCall(DMCopyDisc(dm, *rdm));
2205: PetscCall(DMGetCoordinateDM(dm, &cdm));
2206: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
2207: PetscCall(DMCopyDisc(cdm, rcdm));
2208: PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
2209: PetscCall(DMCopyDisc(dm, *rdm));
2210: PetscCall(DMPlexTransformDestroy(&tr));
2211: ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2212: PetscFunctionReturn(PETSC_SUCCESS);
2213: }