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);
94: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform);
96: /*@C
97: DMPlexTransformRegisterAll - Registers all of the transform components in the `DM` package.
99: Not Collective
101: Level: advanced
103: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
104: @*/
105: PetscErrorCode DMPlexTransformRegisterAll(void)
106: {
107: PetscFunctionBegin;
108: if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
109: DMPlexTransformRegisterAllCalled = PETSC_TRUE;
111: PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter));
112: PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular));
113: PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox));
114: PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld));
115: PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL));
116: PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR));
117: PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D));
118: PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDE, DMPlexTransformCreate_Extrude));
119: PetscCall(DMPlexTransformRegister(DMPLEXCOHESIVEEXTRUDE, DMPlexTransformCreate_Cohesive));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: /*@C
124: DMPlexTransformRegisterDestroy - This function destroys the registered `DMPlexTransformType`. It is called from `PetscFinalize()`.
126: Not collective
128: Level: developer
130: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()`
131: @*/
132: PetscErrorCode DMPlexTransformRegisterDestroy(void)
133: {
134: PetscFunctionBegin;
135: PetscCall(PetscFunctionListDestroy(&DMPlexTransformList));
136: DMPlexTransformRegisterAllCalled = PETSC_FALSE;
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /*@
141: DMPlexTransformCreate - Creates an empty transform object. The type can then be set with `DMPlexTransformSetType()`.
143: Collective
145: Input Parameter:
146: . comm - The communicator for the transform object
148: Output Parameter:
149: . tr - The transform object
151: Level: beginner
153: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
154: @*/
155: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
156: {
157: DMPlexTransform t;
159: PetscFunctionBegin;
160: PetscAssertPointer(tr, 2);
161: *tr = NULL;
162: PetscCall(DMInitializePackage());
164: PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView));
165: t->setupcalled = PETSC_FALSE;
166: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom));
167: *tr = t;
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: /*@
172: DMPlexTransformSetType - Sets the particular implementation for a transform.
174: Collective
176: Input Parameters:
177: + tr - The transform
178: - method - The name of the transform type
180: Options Database Key:
181: . -dm_plex_transform_type <type> - Sets the transform type; see `DMPlexTransformType`
183: Level: intermediate
185: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
186: @*/
187: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
188: {
189: PetscErrorCode (*r)(DMPlexTransform);
190: PetscBool match;
192: PetscFunctionBegin;
194: PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match));
195: if (match) PetscFunctionReturn(PETSC_SUCCESS);
197: PetscCall(DMPlexTransformRegisterAll());
198: PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r));
199: PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);
201: PetscTryTypeMethod(tr, destroy);
202: PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops)));
203: PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method));
204: PetscCall((*r)(tr));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*@
209: DMPlexTransformGetType - Gets the type name (as a string) from the transform.
211: Not Collective
213: Input Parameter:
214: . tr - The `DMPlexTransform`
216: Output Parameter:
217: . type - The `DMPlexTransformType` name
219: Level: intermediate
221: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
222: @*/
223: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
224: {
225: PetscFunctionBegin;
227: PetscAssertPointer(type, 2);
228: PetscCall(DMPlexTransformRegisterAll());
229: *type = ((PetscObject)tr)->type_name;
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
234: {
235: PetscViewerFormat format;
237: PetscFunctionBegin;
238: PetscCall(PetscViewerGetFormat(v, &format));
239: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
240: const PetscInt *trTypes = NULL;
241: IS trIS;
242: PetscInt cols = 8;
243: PetscInt Nrt = 8, f, g;
245: if (tr->trType) PetscCall(DMLabelView(tr->trType, v));
246: PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n"));
247: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
248: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
249: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]));
250: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
251: PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n"));
252: for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
253: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
254: for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]));
255: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
257: if (tr->trType) {
258: PetscCall(DMLabelGetNumValues(tr->trType, &Nrt));
259: PetscCall(DMLabelGetValueIS(tr->trType, &trIS));
260: PetscCall(ISGetIndices(trIS, &trTypes));
261: }
262: PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n"));
263: PetscCall(PetscViewerASCIIPrintf(v, " "));
264: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
265: PetscCall(PetscViewerASCIIPrintf(v, "\n"));
266: for (f = 0; f < Nrt; ++f) {
267: PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT " |", trTypes ? trTypes[f] : f));
268: for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]));
269: PetscCall(PetscViewerASCIIPrintf(v, " |\n"));
270: }
271: if (tr->trType) {
272: PetscCall(ISRestoreIndices(trIS, &trTypes));
273: PetscCall(ISDestroy(&trIS));
274: }
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@
280: DMPlexTransformView - Views a `DMPlexTransform`
282: Collective
284: Input Parameters:
285: + tr - the `DMPlexTransform` object to view
286: - v - the viewer
288: Level: beginner
290: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
291: @*/
292: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
293: {
294: PetscBool isascii;
296: PetscFunctionBegin;
298: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v));
300: PetscCheckSameComm(tr, 1, v, 2);
301: PetscCall(PetscViewerCheckWritable(v));
302: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v));
303: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
304: if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v));
305: PetscTryTypeMethod(tr, view, v);
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*@
310: DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database
312: Collective
314: Input Parameter:
315: . tr - the `DMPlexTransform` object to set options for
317: Options Database Keys:
318: + -dm_plex_transform_type - Set the transform type, e.g. refine_regular
319: . -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point
320: - -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
322: Level: intermediate
324: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
325: @*/
326: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
327: {
328: char typeName[1024], active[PETSC_MAX_PATH_LEN];
329: const char *defName = DMPLEXREFINEREGULAR;
330: PetscBool flg;
332: PetscFunctionBegin;
334: PetscObjectOptionsBegin((PetscObject)tr);
335: PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg));
336: if (flg) PetscCall(DMPlexTransformSetType(tr, typeName));
337: else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName));
338: PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &tr->labelMatchStrata, NULL));
339: 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));
340: PetscCall(PetscOptionsString("-dm_plex_transform_active", "Name for active mesh label", "DMPlexTransformSetActive", active, active, sizeof(active), &flg));
341: if (flg) {
342: DM dm;
343: DMLabel label;
345: PetscCall(DMPlexTransformGetDM(tr, &dm));
346: PetscCall(DMGetLabel(dm, active, &label));
347: PetscCall(DMPlexTransformSetActive(tr, label));
348: }
349: PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
350: /* process any options handlers added with PetscObjectAddOptionsHandler() */
351: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject));
352: PetscOptionsEnd();
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@
357: DMPlexTransformDestroy - Destroys a `DMPlexTransform`
359: Collective
361: Input Parameter:
362: . tr - the transform object to destroy
364: Level: beginner
366: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
367: @*/
368: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
369: {
370: PetscInt c;
372: PetscFunctionBegin;
373: if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
375: if (--((PetscObject)*tr)->refct > 0) {
376: *tr = NULL;
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: PetscTryTypeMethod(*tr, destroy);
381: PetscCall(DMDestroy(&(*tr)->dm));
382: PetscCall(DMLabelDestroy(&(*tr)->active));
383: PetscCall(DMLabelDestroy(&(*tr)->trType));
384: PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld));
385: PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew));
386: PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew));
387: PetscCall(PetscFree((*tr)->offset));
388: PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd));
389: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
390: PetscCall(PetscFEDestroy(&(*tr)->coordFE[c]));
391: PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c]));
392: }
393: if ((*tr)->trVerts) {
394: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
395: DMPolytopeType *rct;
396: PetscInt *rsize, *rcone, *rornt, Nct, n, r;
398: if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) {
399: PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
400: for (n = 0; n < Nct; ++n) {
401: if (rct[n] == DM_POLYTOPE_POINT) continue;
402: for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r]));
403: PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]]));
404: }
405: }
406: PetscCall(PetscFree((*tr)->trSubVerts[c]));
407: PetscCall(PetscFree((*tr)->trVerts[c]));
408: }
409: }
410: PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts));
411: PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom));
412: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
413: PetscCall(PetscHeaderDestroy(tr));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
418: {
419: DMLabel trType = tr->trType;
420: PetscInt c, cN, *off;
422: PetscFunctionBegin;
423: if (trType) {
424: DM dm;
425: IS rtIS;
426: const PetscInt *reftypes;
427: PetscInt Nrt, r;
429: PetscCall(DMPlexTransformGetDM(tr, &dm));
430: PetscCall(DMLabelGetNumValues(trType, &Nrt));
431: PetscCall(DMLabelGetValueIS(trType, &rtIS));
432: PetscCall(ISGetIndices(rtIS, &reftypes));
433: PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off));
434: for (r = 0; r < Nrt; ++r) {
435: const PetscInt rt = reftypes[r];
436: IS rtIS;
437: const PetscInt *points;
438: DMPolytopeType ct;
439: PetscInt np, p;
441: PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
442: PetscCall(ISGetLocalSize(rtIS, &np));
443: PetscCall(ISGetIndices(rtIS, &points));
444: if (!np) continue;
445: p = points[0];
446: PetscCall(ISRestoreIndices(rtIS, &points));
447: PetscCall(ISDestroy(&rtIS));
448: PetscCall(DMPlexGetCellType(dm, p, &ct));
449: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
450: const DMPolytopeType ctNew = (DMPolytopeType)cN;
451: DMPolytopeType *rct;
452: PetscInt *rsize, *cone, *ornt;
453: PetscInt Nct, n, s;
455: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
456: off[r * DM_NUM_POLYTOPES + ctNew] = -1;
457: break;
458: }
459: off[r * DM_NUM_POLYTOPES + ctNew] = 0;
460: for (s = 0; s <= r; ++s) {
461: const PetscInt st = reftypes[s];
462: DMPolytopeType sct;
463: PetscInt q, qrt;
465: PetscCall(DMLabelGetStratumIS(trType, st, &rtIS));
466: PetscCall(ISGetLocalSize(rtIS, &np));
467: PetscCall(ISGetIndices(rtIS, &points));
468: if (!np) continue;
469: q = points[0];
470: PetscCall(ISRestoreIndices(rtIS, &points));
471: PetscCall(ISDestroy(&rtIS));
472: PetscCall(DMPlexGetCellType(dm, q, &sct));
473: PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt));
474: 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);
475: if (st == rt) {
476: for (n = 0; n < Nct; ++n)
477: if (rct[n] == ctNew) break;
478: if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
479: break;
480: }
481: for (n = 0; n < Nct; ++n) {
482: if (rct[n] == ctNew) {
483: PetscInt sn;
485: PetscCall(DMLabelGetStratumSize(trType, st, &sn));
486: off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
487: }
488: }
489: }
490: }
491: }
492: PetscCall(ISRestoreIndices(rtIS, &reftypes));
493: PetscCall(ISDestroy(&rtIS));
494: } else {
495: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off));
496: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
497: const DMPolytopeType ct = (DMPolytopeType)c;
498: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
499: const DMPolytopeType ctNew = (DMPolytopeType)cN;
500: DMPolytopeType *rct;
501: PetscInt *rsize, *cone, *ornt;
502: PetscInt Nct, n, i;
504: 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) {
505: off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
506: continue;
507: }
508: off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
509: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
510: const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i];
511: const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];
513: PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt));
514: if (ict == ct) {
515: for (n = 0; n < Nct; ++n)
516: if (rct[n] == ctNew) break;
517: if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
518: break;
519: }
520: for (n = 0; n < Nct; ++n)
521: if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
522: }
523: }
524: }
525: }
526: *offset = off;
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /*@
531: DMPlexTransformSetUp - Create the tables that drive the transform
533: Input Parameter:
534: . tr - The `DMPlexTransform` object
536: Level: intermediate
538: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
539: @*/
540: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
541: {
542: DM dm;
543: DMPolytopeType ctCell;
544: PetscInt pStart, pEnd, p, c, celldim = 0;
546: PetscFunctionBegin;
548: if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
549: PetscTryTypeMethod(tr, setup);
550: PetscCall(DMPlexTransformGetDM(tr, &dm));
551: PetscCall(DMSetSnapToGeomModel(dm, NULL));
552: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
553: if (pEnd > pStart) {
554: // Ignore cells hanging off of embedded surfaces
555: PetscInt c = pStart;
557: ctCell = DM_POLYTOPE_FV_GHOST;
558: while (DMPolytopeTypeGetDim(ctCell) < 0) PetscCall(DMPlexGetCellType(dm, c++, &ctCell));
559: } else {
560: PetscInt dim;
562: PetscCall(DMGetDimension(dm, &dim));
563: switch (dim) {
564: case 0:
565: ctCell = DM_POLYTOPE_POINT;
566: break;
567: case 1:
568: ctCell = DM_POLYTOPE_SEGMENT;
569: break;
570: case 2:
571: ctCell = DM_POLYTOPE_TRIANGLE;
572: break;
573: case 3:
574: ctCell = DM_POLYTOPE_TETRAHEDRON;
575: break;
576: default:
577: ctCell = DM_POLYTOPE_UNKNOWN;
578: }
579: }
580: PetscCall(DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
581: for (p = pStart; p < pEnd; ++p) {
582: DMPolytopeType ct;
583: DMPolytopeType *rct;
584: PetscInt *rsize, *cone, *ornt;
585: PetscInt Nct, n;
587: PetscCall(DMPlexGetCellType(dm, p, &ct));
588: 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);
589: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
590: for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
591: }
592: PetscCall(DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
593: /* Construct sizes and offsets for each cell type */
594: if (!tr->ctStart) {
595: PetscInt *ctS, *ctSN, *ctC, *ctCN;
597: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
598: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
599: for (p = pStart; p < pEnd; ++p) {
600: DMPolytopeType ct;
601: DMPolytopeType *rct;
602: PetscInt *rsize, *cone, *ornt;
603: PetscInt Nct, n;
605: PetscCall(DMPlexGetCellType(dm, p, &ct));
606: 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);
607: ++ctC[ct];
608: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
609: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
610: }
611: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
612: const PetscInt cto = tr->ctOrderOld[c];
613: const PetscInt cton = tr->ctOrderOld[c + 1];
614: const PetscInt ctn = tr->ctOrderNew[c];
615: const PetscInt ctnn = tr->ctOrderNew[c + 1];
617: ctS[cton] = ctS[cto] + ctC[cto];
618: ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
619: }
620: PetscCall(PetscFree2(ctC, ctCN));
621: tr->ctStart = ctS;
622: tr->ctStartNew = ctSN;
623: }
624: PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
625: // Compute depth information
626: tr->depth = -1;
627: for (c = 0; c < DM_NUM_POLYTOPES; ++c)
628: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
629: PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
630: for (PetscInt d = 0; d <= tr->depth; ++d) {
631: tr->depthStart[d] = PETSC_INT_MAX;
632: tr->depthEnd[d] = -1;
633: }
634: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
635: const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);
637: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
638: tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
639: tr->depthEnd[dep] = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
640: }
641: tr->setupcalled = PETSC_TRUE;
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: /*@
646: DMPlexTransformGetDM - Get the base `DM` for the transform
648: Input Parameter:
649: . tr - The `DMPlexTransform` object
651: Output Parameter:
652: . dm - The original `DM` which will be transformed
654: Level: intermediate
656: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
657: @*/
658: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
659: {
660: PetscFunctionBegin;
662: PetscAssertPointer(dm, 2);
663: *dm = tr->dm;
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: /*@
668: DMPlexTransformSetDM - Set the base `DM` for the transform
670: Input Parameters:
671: + tr - The `DMPlexTransform` object
672: - dm - The original `DM` which will be transformed
674: Level: intermediate
676: Note:
677: The user does not typically call this, as it is called by `DMPlexTransformApply()`.
679: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
680: @*/
681: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
682: {
683: PetscFunctionBegin;
686: PetscCall(PetscObjectReference((PetscObject)dm));
687: PetscCall(DMDestroy(&tr->dm));
688: tr->dm = dm;
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: /*@
693: DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform
695: Input Parameter:
696: . tr - The `DMPlexTransform` object
698: Output Parameter:
699: . active - The `DMLabel` indicating which points will be transformed
701: Level: intermediate
703: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
704: @*/
705: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
706: {
707: PetscFunctionBegin;
709: PetscAssertPointer(active, 2);
710: *active = tr->active;
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: /*@
715: DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform
717: Input Parameters:
718: + tr - The `DMPlexTransform` object
719: - active - The original `DM` which will be transformed
721: Level: intermediate
723: Note:
724: This only applies to transforms that can operator on a subset of the mesh, listed in [](plex_transform_table).
726: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
727: @*/
728: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
729: {
730: PetscFunctionBegin;
733: PetscCall(PetscObjectReference((PetscObject)active));
734: PetscCall(DMLabelDestroy(&tr->active));
735: tr->active = active;
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
740: {
741: PetscFunctionBegin;
742: if (!tr->coordFE[ct]) {
743: PetscInt dim, cdim;
745: dim = DMPolytopeTypeGetDim(ct);
746: PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
747: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
748: {
749: PetscDualSpace dsp;
750: PetscQuadrature quad;
751: DM K;
752: PetscFEGeom *cg;
753: PetscScalar *Xq;
754: PetscReal *xq, *wq;
755: PetscInt Nq, q;
757: PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
758: PetscCall(PetscMalloc1(Nq * cdim, &xq));
759: for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
760: PetscCall(PetscMalloc1(Nq, &wq));
761: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
762: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
763: PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
764: PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));
766: PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
767: PetscCall(PetscDualSpaceGetDM(dsp, &K));
768: PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]));
769: cg = tr->refGeom[ct];
770: PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
771: PetscCall(PetscQuadratureDestroy(&quad));
772: }
773: }
774: *fe = tr->coordFE[ct];
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
779: {
780: PetscInt dim, cdim;
782: PetscFunctionBegin;
783: PetscCall(DMGetDimension(dm, &dim));
784: PetscCall(DMSetDimension(tdm, dim));
785: PetscCall(DMGetCoordinateDim(dm, &cdim));
786: PetscCall(DMSetCoordinateDim(tdm, cdim));
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: /*@
791: DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`
793: Input Parameters:
794: + tr - The `DMPlexTransform` object
795: - dm - The original `DM`
797: Output Parameter:
798: . tdm - The transformed `DM`
800: Level: advanced
802: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
803: @*/
804: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
805: {
806: PetscFunctionBegin;
807: PetscUseTypeMethod(tr, setdimensions, dm, tdm);
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
812: {
813: PetscFunctionBegin;
814: if (pStart) *pStart = 0;
815: if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
820: {
821: PetscInt ctNew;
823: PetscFunctionBegin;
825: PetscAssertPointer(celltype, 3);
826: /* TODO Can do bisection since everything is sorted */
827: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
828: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
830: if (cell >= ctSN && cell < ctEN) break;
831: }
832: PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
833: *celltype = (DMPolytopeType)ctNew;
834: PetscFunctionReturn(PETSC_SUCCESS);
835: }
837: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
838: {
839: PetscFunctionBegin;
841: if (start) *start = tr->ctStartNew[celltype];
842: if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
843: PetscFunctionReturn(PETSC_SUCCESS);
844: }
846: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
847: {
848: PetscFunctionBegin;
850: *depth = tr->depth;
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
855: {
856: PetscFunctionBegin;
858: if (start) *start = tr->depthStart[depth];
859: if (end) *end = tr->depthEnd[depth];
860: PetscFunctionReturn(PETSC_SUCCESS);
861: }
863: /*@
864: DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.
866: Not Collective
868: Input Parameters:
869: + tr - The `DMPlexTransform`
870: . ct - The type of the original point which produces the new point
871: . ctNew - The type of the new point
872: . p - The original point which produces the new point
873: - r - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`
875: Output Parameter:
876: . pNew - The new point number
878: Level: developer
880: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
881: @*/
882: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
883: {
884: DMPolytopeType *rct;
885: PetscInt *rsize, *cone, *ornt;
886: PetscInt rt, Nct, n, off, rp;
887: DMLabel trType = tr->trType;
888: PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
889: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
890: PetscInt newp = ctSN, cind;
892: PetscFunctionBeginHot;
893: 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);
894: PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
895: if (trType) {
896: PetscCall(DMLabelGetValueIndex(trType, rt, &cind));
897: PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
898: 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);
899: } else {
900: cind = ct;
901: rp = p - ctS;
902: }
903: off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
904: 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);
905: newp += off;
906: for (n = 0; n < Nct; ++n) {
907: if (rct[n] == ctNew) {
908: if (rsize[n] && r >= rsize[n])
909: 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]);
910: newp += rp * rsize[n] + r;
911: break;
912: }
913: }
915: 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);
916: *pNew = newp;
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: /*@
921: DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.
923: Not Collective
925: Input Parameters:
926: + tr - The `DMPlexTransform`
927: - pNew - The new point number
929: Output Parameters:
930: + ct - The type of the original point which produces the new point
931: . ctNew - The type of the new point
932: . p - The original point which produces the new point
933: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
935: Level: developer
937: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
938: @*/
939: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
940: {
941: DMLabel trType = tr->trType;
942: DMPolytopeType *rct, ctN;
943: PetscInt *rsize, *cone, *ornt;
944: PetscInt rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
945: PetscInt offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;
947: PetscFunctionBegin;
948: PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
949: if (trType) {
950: DM dm;
951: IS rtIS;
952: const PetscInt *reftypes;
953: PetscInt Nrt, r, rtStart;
955: PetscCall(DMPlexTransformGetDM(tr, &dm));
956: PetscCall(DMLabelGetNumValues(trType, &Nrt));
957: PetscCall(DMLabelGetValueIS(trType, &rtIS));
958: PetscCall(ISGetIndices(rtIS, &reftypes));
959: for (r = 0; r < Nrt; ++r) {
960: const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];
962: if (tr->ctStartNew[ctN] + off > pNew) continue;
963: /* Check that any of this refinement type exist */
964: /* TODO Actually keep track of the number produced here instead */
965: if (off > offset) {
966: rt = reftypes[r];
967: offset = off;
968: }
969: }
970: PetscCall(ISRestoreIndices(rtIS, &reftypes));
971: PetscCall(ISDestroy(&rtIS));
972: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
973: /* TODO Map refinement types to cell types */
974: PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL));
975: PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
976: for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
977: PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
979: if ((rtStart >= ctS) && (rtStart < ctE)) break;
980: }
981: PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt);
982: } else {
983: for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
984: const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];
986: if (tr->ctStartNew[ctN] + off > pNew) continue;
987: if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
988: /* TODO Actually keep track of the number produced here instead */
989: if (off > offset) {
990: ctO = ctTmp;
991: offset = off;
992: }
993: }
994: rt = -1;
995: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
996: }
997: ctS = tr->ctStart[ctO];
998: ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
999: if (trType) {
1000: for (rtS = ctS; rtS < ctE; ++rtS) {
1001: PetscInt val;
1002: PetscCall(DMLabelGetValue(trType, rtS, &val));
1003: if (val == rt) break;
1004: }
1005: 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);
1006: } else rtS = ctS;
1007: PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
1008: 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);
1009: for (n = 0; n < Nct; ++n) {
1010: if (rct[n] == ctN) {
1011: PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c;
1013: if (trType) {
1014: for (c = ctS; c < ctE; ++c) {
1015: PetscCall(DMLabelGetValue(trType, c, &val));
1016: if (val == rt) {
1017: if (tmp < rsize[n]) break;
1018: tmp -= rsize[n];
1019: }
1020: }
1021: PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew);
1022: rp = c - ctS;
1023: rO = tmp;
1024: } else {
1025: // This assumes that all points of type ctO transform the same way
1026: rp = tmp / rsize[n];
1027: rO = tmp % rsize[n];
1028: }
1029: break;
1030: }
1031: }
1032: PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
1033: pO = rp + ctS;
1034: 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);
1035: if (ct) *ct = (DMPolytopeType)ctO;
1036: if (ctNew) *ctNew = (DMPolytopeType)ctN;
1037: if (p) *p = pO;
1038: if (r) *r = rO;
1039: PetscFunctionReturn(PETSC_SUCCESS);
1040: }
1042: /*@
1043: DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.
1045: Input Parameters:
1046: + tr - The `DMPlexTransform` object
1047: . source - The source cell type
1048: - p - The source point, which can also determine the refine type
1050: Output Parameters:
1051: + rt - The refine type for this point
1052: . Nt - The number of types produced by this point
1053: . target - An array of length `Nt` giving the types produced
1054: . size - An array of length `Nt` giving the number of cells of each type produced
1055: . cone - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
1056: - ornt - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point
1058: Level: advanced
1060: Notes:
1061: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
1062: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
1063: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
1064: .vb
1065: the number of cones to be taken, or 0 for the current cell
1066: the cell cone point number at each level from which it is subdivided
1067: the replica number r of the subdivision.
1068: .ve
1069: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
1070: .vb
1071: Nt = 2
1072: target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
1073: size = {1, 2}
1074: cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
1075: ornt = { 0, 0, 0, 0}
1076: .ve
1078: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
1079: @*/
1080: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1081: {
1082: PetscFunctionBegin;
1083: PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
1084: PetscFunctionReturn(PETSC_SUCCESS);
1085: }
1087: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1088: {
1089: PetscFunctionBegin;
1090: *rnew = r;
1091: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
1092: PetscFunctionReturn(PETSC_SUCCESS);
1093: }
1095: /* Returns the same thing */
1096: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1097: {
1098: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
1099: static PetscInt vertexS[] = {1};
1100: static PetscInt vertexC[] = {0};
1101: static PetscInt vertexO[] = {0};
1102: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
1103: static PetscInt edgeS[] = {1};
1104: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1105: static PetscInt edgeO[] = {0, 0};
1106: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
1107: static PetscInt tedgeS[] = {1};
1108: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1109: static PetscInt tedgeO[] = {0, 0};
1110: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
1111: static PetscInt triS[] = {1};
1112: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1113: static PetscInt triO[] = {0, 0, 0};
1114: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
1115: static PetscInt quadS[] = {1};
1116: 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};
1117: static PetscInt quadO[] = {0, 0, 0, 0};
1118: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
1119: static PetscInt tquadS[] = {1};
1120: 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};
1121: static PetscInt tquadO[] = {0, 0, 0, 0};
1122: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
1123: static PetscInt tetS[] = {1};
1124: 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};
1125: static PetscInt tetO[] = {0, 0, 0, 0};
1126: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
1127: static PetscInt hexS[] = {1};
1128: 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};
1129: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
1130: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
1131: static PetscInt tripS[] = {1};
1132: 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};
1133: static PetscInt tripO[] = {0, 0, 0, 0, 0};
1134: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
1135: static PetscInt ttripS[] = {1};
1136: 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};
1137: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
1138: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
1139: static PetscInt tquadpS[] = {1};
1140: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
1141: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1142: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
1143: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
1144: static PetscInt pyrS[] = {1};
1145: 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};
1146: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
1148: PetscFunctionBegin;
1149: if (rt) *rt = 0;
1150: switch (source) {
1151: case DM_POLYTOPE_POINT:
1152: *Nt = 1;
1153: *target = vertexT;
1154: *size = vertexS;
1155: *cone = vertexC;
1156: *ornt = vertexO;
1157: break;
1158: case DM_POLYTOPE_SEGMENT:
1159: *Nt = 1;
1160: *target = edgeT;
1161: *size = edgeS;
1162: *cone = edgeC;
1163: *ornt = edgeO;
1164: break;
1165: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1166: *Nt = 1;
1167: *target = tedgeT;
1168: *size = tedgeS;
1169: *cone = tedgeC;
1170: *ornt = tedgeO;
1171: break;
1172: case DM_POLYTOPE_TRIANGLE:
1173: *Nt = 1;
1174: *target = triT;
1175: *size = triS;
1176: *cone = triC;
1177: *ornt = triO;
1178: break;
1179: case DM_POLYTOPE_QUADRILATERAL:
1180: *Nt = 1;
1181: *target = quadT;
1182: *size = quadS;
1183: *cone = quadC;
1184: *ornt = quadO;
1185: break;
1186: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1187: *Nt = 1;
1188: *target = tquadT;
1189: *size = tquadS;
1190: *cone = tquadC;
1191: *ornt = tquadO;
1192: break;
1193: case DM_POLYTOPE_TETRAHEDRON:
1194: *Nt = 1;
1195: *target = tetT;
1196: *size = tetS;
1197: *cone = tetC;
1198: *ornt = tetO;
1199: break;
1200: case DM_POLYTOPE_HEXAHEDRON:
1201: *Nt = 1;
1202: *target = hexT;
1203: *size = hexS;
1204: *cone = hexC;
1205: *ornt = hexO;
1206: break;
1207: case DM_POLYTOPE_TRI_PRISM:
1208: *Nt = 1;
1209: *target = tripT;
1210: *size = tripS;
1211: *cone = tripC;
1212: *ornt = tripO;
1213: break;
1214: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1215: *Nt = 1;
1216: *target = ttripT;
1217: *size = ttripS;
1218: *cone = ttripC;
1219: *ornt = ttripO;
1220: break;
1221: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1222: *Nt = 1;
1223: *target = tquadpT;
1224: *size = tquadpS;
1225: *cone = tquadpC;
1226: *ornt = tquadpO;
1227: break;
1228: case DM_POLYTOPE_PYRAMID:
1229: *Nt = 1;
1230: *target = pyrT;
1231: *size = pyrS;
1232: *cone = pyrC;
1233: *ornt = pyrO;
1234: break;
1235: default:
1236: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1237: }
1238: PetscFunctionReturn(PETSC_SUCCESS);
1239: }
1241: /*@
1242: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
1244: Not Collective
1246: Input Parameters:
1247: + tr - The `DMPlexTransform`
1248: . sct - The source point cell type, from whom the new cell is being produced
1249: . sp - The source point
1250: . so - The orientation of the source point in its enclosing parent
1251: . tct - The target point cell type
1252: . r - The replica number requested for the produced cell type
1253: - o - The orientation of the replica
1255: Output Parameters:
1256: + rnew - The replica number, given the orientation of the parent
1257: - onew - The replica orientation, given the orientation of the parent
1259: Level: advanced
1261: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1262: @*/
1263: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1264: {
1265: PetscFunctionBeginHot;
1266: PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1271: {
1272: DM dm;
1273: PetscInt pStart, pEnd, p, pNew;
1275: PetscFunctionBegin;
1276: PetscCall(DMPlexTransformGetDM(tr, &dm));
1277: /* Must create the celltype label here so that we do not automatically try to compute the types */
1278: PetscCall(DMCreateLabel(rdm, "celltype"));
1279: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1280: for (p = pStart; p < pEnd; ++p) {
1281: DMPolytopeType ct;
1282: DMPolytopeType *rct;
1283: PetscInt *rsize, *rcone, *rornt;
1284: PetscInt Nct, n, r;
1286: PetscCall(DMPlexGetCellType(dm, p, &ct));
1287: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1288: for (n = 0; n < Nct; ++n) {
1289: for (r = 0; r < rsize[n]; ++r) {
1290: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1291: PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
1292: PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
1293: }
1294: }
1295: }
1296: /* Let the DM know we have set all the cell types */
1297: {
1298: DMLabel ctLabel;
1299: DM_Plex *plex = (DM_Plex *)rdm->data;
1301: PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
1302: PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
1303: }
1304: PetscFunctionReturn(PETSC_SUCCESS);
1305: }
1307: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1308: {
1309: DMPolytopeType ctNew;
1311: PetscFunctionBegin;
1313: PetscAssertPointer(coneSize, 3);
1314: PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
1315: *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType)ctNew);
1316: PetscFunctionReturn(PETSC_SUCCESS);
1317: }
1319: /* The orientation o is for the interior of the cell p */
1320: 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[])
1321: {
1322: DM dm;
1323: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1324: const PetscInt *cone;
1325: PetscInt c, coff = *coneoff, ooff = *orntoff;
1327: PetscFunctionBegin;
1328: PetscCall(DMPlexTransformGetDM(tr, &dm));
1329: PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
1330: for (c = 0; c < csizeNew; ++c) {
1331: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1332: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1333: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1334: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1335: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1336: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1337: const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1338: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1339: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1340: PetscInt lc;
1342: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1343: for (lc = 0; lc < fn; ++lc) {
1344: const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po);
1345: const PetscInt acp = rcone[coff++];
1346: const PetscInt pcp = parr[acp * 2];
1347: const PetscInt pco = parr[acp * 2 + 1];
1348: const PetscInt *ppornt;
1350: ppp = pp;
1351: pp = pcone[pcp];
1352: PetscCall(DMPlexGetCellType(dm, pp, &pct));
1353: // Restore the parent cone from the last iterate
1354: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
1355: PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
1356: PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
1357: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1358: PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
1359: }
1360: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
1361: pr = rcone[coff++];
1362: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1363: PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
1364: PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
1365: orntNew[c] = fo;
1366: }
1367: PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
1368: *coneoff = coff;
1369: *orntoff = ooff;
1370: PetscFunctionReturn(PETSC_SUCCESS);
1371: }
1373: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1374: {
1375: DM dm;
1376: DMPolytopeType ct;
1377: PetscInt *coneNew, *orntNew;
1378: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1380: PetscFunctionBegin;
1381: PetscCall(DMPlexTransformGetDM(tr, &dm));
1382: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1383: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1384: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1385: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1386: for (p = pStart; p < pEnd; ++p) {
1387: PetscInt coff, ooff;
1388: DMPolytopeType *rct;
1389: PetscInt *rsize, *rcone, *rornt;
1390: PetscInt Nct, n, r;
1392: PetscCall(DMPlexGetCellType(dm, p, &ct));
1393: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1394: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1395: const DMPolytopeType ctNew = rct[n];
1397: for (r = 0; r < rsize[n]; ++r) {
1398: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1399: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
1400: PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
1401: PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
1402: }
1403: }
1404: }
1405: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1406: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1407: PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
1408: PetscCall(DMPlexSymmetrize(rdm));
1409: PetscCall(DMPlexStratify(rdm));
1410: PetscFunctionReturn(PETSC_SUCCESS);
1411: }
1413: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1414: {
1415: DM dm;
1416: DMPolytopeType ct, qct;
1417: DMPolytopeType *rct;
1418: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1419: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1421: PetscFunctionBegin;
1423: PetscAssertPointer(cone, 4);
1424: PetscAssertPointer(ornt, 5);
1425: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1426: PetscCall(DMPlexTransformGetDM(tr, &dm));
1427: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1428: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1429: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1430: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1431: for (n = 0; n < Nct; ++n) {
1432: const DMPolytopeType ctNew = rct[n];
1433: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1434: PetscInt Nr = rsize[n], fn, c;
1436: if (ctNew == qct) Nr = r;
1437: for (nr = 0; nr < Nr; ++nr) {
1438: for (c = 0; c < csizeNew; ++c) {
1439: ++coff; /* Cell type of new cone point */
1440: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1441: coff += fn;
1442: ++coff; /* Replica number of new cone point */
1443: ++ooff; /* Orientation of new cone point */
1444: }
1445: }
1446: if (ctNew == qct) break;
1447: }
1448: PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1449: *cone = qcone;
1450: *ornt = qornt;
1451: PetscFunctionReturn(PETSC_SUCCESS);
1452: }
1454: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1455: {
1456: DM dm;
1457: DMPolytopeType ct, qct;
1458: DMPolytopeType *rct;
1459: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1460: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1462: PetscFunctionBegin;
1464: if (cone) PetscAssertPointer(cone, 3);
1465: if (ornt) PetscAssertPointer(ornt, 4);
1466: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1467: PetscCall(DMPlexTransformGetDM(tr, &dm));
1468: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1469: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1470: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1471: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1472: for (n = 0; n < Nct; ++n) {
1473: const DMPolytopeType ctNew = rct[n];
1474: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1475: PetscInt Nr = rsize[n], fn, c;
1477: if (ctNew == qct) Nr = r;
1478: for (nr = 0; nr < Nr; ++nr) {
1479: for (c = 0; c < csizeNew; ++c) {
1480: ++coff; /* Cell type of new cone point */
1481: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1482: coff += fn;
1483: ++coff; /* Replica number of new cone point */
1484: ++ooff; /* Orientation of new cone point */
1485: }
1486: }
1487: if (ctNew == qct) break;
1488: }
1489: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1490: if (cone) *cone = qcone;
1491: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1492: if (ornt) *ornt = qornt;
1493: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1494: PetscFunctionReturn(PETSC_SUCCESS);
1495: }
1497: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1498: {
1499: DM dm;
1501: PetscFunctionBegin;
1503: PetscCall(DMPlexTransformGetDM(tr, &dm));
1504: if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
1505: if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
1506: PetscFunctionReturn(PETSC_SUCCESS);
1507: }
1509: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1510: {
1511: PetscInt ict;
1513: PetscFunctionBegin;
1514: PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
1515: for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1516: const DMPolytopeType ct = (DMPolytopeType)ict;
1517: DMPlexTransform reftr;
1518: DM refdm, trdm;
1519: Vec coordinates;
1520: const PetscScalar *coords;
1521: DMPolytopeType *rct;
1522: PetscInt *rsize, *rcone, *rornt;
1523: PetscInt Nct, n, r, pNew = 0;
1524: PetscInt trdim, vStart, vEnd, Nc;
1525: const PetscInt debug = 0;
1526: const char *typeName;
1528: /* Since points are 0-dimensional, coordinates make no sense */
1529: if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue;
1530: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
1531: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
1532: PetscCall(DMPlexTransformSetDM(reftr, refdm));
1533: PetscCall(DMPlexTransformGetType(tr, &typeName));
1534: PetscCall(DMPlexTransformSetType(reftr, typeName));
1535: PetscCall(DMPlexTransformSetUp(reftr));
1536: PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));
1538: PetscCall(DMGetDimension(trdm, &trdim));
1539: PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
1540: tr->trNv[ct] = vEnd - vStart;
1541: PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
1542: PetscCall(VecGetLocalSize(coordinates, &Nc));
1543: 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);
1544: PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
1545: PetscCall(VecGetArrayRead(coordinates, &coords));
1546: PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
1547: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1549: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
1550: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1551: for (n = 0; n < Nct; ++n) {
1552: /* Since points are 0-dimensional, coordinates make no sense */
1553: if (rct[n] == DM_POLYTOPE_POINT) continue;
1554: PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
1555: for (r = 0; r < rsize[n]; ++r) {
1556: PetscInt *closure = NULL;
1557: PetscInt clSize, cl, Nv = 0;
1559: PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
1560: PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
1561: PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1562: for (cl = 0; cl < clSize * 2; cl += 2) {
1563: const PetscInt sv = closure[cl];
1565: if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1566: }
1567: PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1568: 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]);
1569: }
1570: }
1571: if (debug) {
1572: DMPolytopeType *rct;
1573: PetscInt *rsize, *rcone, *rornt;
1574: PetscInt v, dE = trdim, d, off = 0;
1576: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
1577: for (v = 0; v < tr->trNv[ct]; ++v) {
1578: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1579: for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
1580: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1581: }
1583: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1584: for (n = 0; n < Nct; ++n) {
1585: if (rct[n] == DM_POLYTOPE_POINT) continue;
1586: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
1587: for (r = 0; r < rsize[n]; ++r) {
1588: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1589: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
1590: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1591: }
1592: }
1593: }
1594: PetscCall(DMDestroy(&refdm));
1595: PetscCall(DMDestroy(&trdm));
1596: PetscCall(DMPlexTransformDestroy(&reftr));
1597: }
1598: PetscFunctionReturn(PETSC_SUCCESS);
1599: }
1601: /*@C
1602: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
1604: Input Parameters:
1605: + tr - The `DMPlexTransform` object
1606: - ct - The cell type
1608: Output Parameters:
1609: + Nv - The number of transformed vertices in the closure of the reference cell of given type
1610: - trVerts - The coordinates of these vertices in the reference cell
1612: Level: developer
1614: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
1615: @*/
1616: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1617: {
1618: PetscFunctionBegin;
1619: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1620: if (Nv) *Nv = tr->trNv[ct];
1621: if (trVerts) *trVerts = tr->trVerts[ct];
1622: PetscFunctionReturn(PETSC_SUCCESS);
1623: }
1625: /*@C
1626: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
1628: Input Parameters:
1629: + tr - The `DMPlexTransform` object
1630: . ct - The cell type
1631: . rct - The subcell type
1632: - r - The subcell index
1634: Output Parameter:
1635: . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()`
1637: Level: developer
1639: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
1640: @*/
1641: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1642: {
1643: PetscFunctionBegin;
1644: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1645: PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
1646: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1647: PetscFunctionReturn(PETSC_SUCCESS);
1648: }
1650: /* Computes new vertex as the barycenter, or centroid */
1651: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1652: {
1653: PetscInt v, d;
1655: PetscFunctionBeginHot;
1656: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
1657: for (d = 0; d < dE; ++d) out[d] = 0.0;
1658: for (v = 0; v < Nv; ++v)
1659: for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1660: for (d = 0; d < dE; ++d) out[d] /= Nv;
1661: PetscFunctionReturn(PETSC_SUCCESS);
1662: }
1664: /*@
1665: DMPlexTransformMapCoordinates - Calculate new coordinates for produced points
1667: Not collective
1669: Input Parameters:
1670: + tr - The `DMPlexTransform`
1671: . pct - The cell type of the parent, from whom the new cell is being produced
1672: . ct - The type being produced
1673: . p - The original point
1674: . r - The replica number requested for the produced cell type
1675: . Nv - Number of vertices in the closure of the parent cell
1676: . dE - Spatial dimension
1677: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
1679: Output Parameter:
1680: . out - The coordinates of the new vertices
1682: Level: intermediate
1684: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
1685: @*/
1686: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1687: {
1688: PetscFunctionBeginHot;
1689: if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1690: PetscFunctionReturn(PETSC_SUCCESS);
1691: }
1693: /*
1694: DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label
1696: Not Collective
1698: Input Parameters:
1699: + tr - The `DMPlexTransform`
1700: . label - The label in the transformed mesh
1701: . pp - The parent point in the original mesh
1702: . pct - The cell type of the parent point
1703: . p - The point in the transformed mesh
1704: . ct - The cell type of the point
1705: . r - The replica number of the point
1706: - val - The label value of the parent point
1708: Level: developer
1710: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
1711: */
1712: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
1713: {
1714: PetscFunctionBeginHot;
1715: if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
1716: PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
1717: PetscFunctionReturn(PETSC_SUCCESS);
1718: }
1720: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1721: {
1722: DM dm;
1723: IS valueIS;
1724: const PetscInt *values;
1725: PetscInt defVal, Nv, val;
1727: PetscFunctionBegin;
1728: PetscCall(DMPlexTransformGetDM(tr, &dm));
1729: PetscCall(DMLabelGetDefaultValue(label, &defVal));
1730: PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
1731: PetscCall(DMLabelGetValueIS(label, &valueIS));
1732: PetscCall(ISGetLocalSize(valueIS, &Nv));
1733: PetscCall(ISGetIndices(valueIS, &values));
1734: for (val = 0; val < Nv; ++val) {
1735: IS pointIS;
1736: const PetscInt *points;
1737: PetscInt numPoints, p;
1739: /* Ensure refined label is created with same number of strata as
1740: * original (even if no entries here). */
1741: PetscCall(DMLabelAddStratum(labelNew, values[val]));
1742: PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
1743: PetscCall(ISGetLocalSize(pointIS, &numPoints));
1744: PetscCall(ISGetIndices(pointIS, &points));
1745: for (p = 0; p < numPoints; ++p) {
1746: const PetscInt point = points[p];
1747: DMPolytopeType ct;
1748: DMPolytopeType *rct;
1749: PetscInt *rsize, *rcone, *rornt;
1750: PetscInt Nct, n, r, pNew = 0;
1752: PetscCall(DMPlexGetCellType(dm, point, &ct));
1753: PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1754: for (n = 0; n < Nct; ++n) {
1755: for (r = 0; r < rsize[n]; ++r) {
1756: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
1757: PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
1758: }
1759: }
1760: }
1761: PetscCall(ISRestoreIndices(pointIS, &points));
1762: PetscCall(ISDestroy(&pointIS));
1763: }
1764: PetscCall(ISRestoreIndices(valueIS, &values));
1765: PetscCall(ISDestroy(&valueIS));
1766: PetscFunctionReturn(PETSC_SUCCESS);
1767: }
1769: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1770: {
1771: DM dm;
1772: PetscInt numLabels, l;
1774: PetscFunctionBegin;
1775: PetscCall(DMPlexTransformGetDM(tr, &dm));
1776: PetscCall(DMGetNumLabels(dm, &numLabels));
1777: for (l = 0; l < numLabels; ++l) {
1778: DMLabel label, labelNew;
1779: const char *lname;
1780: PetscBool isDepth, isCellType;
1782: PetscCall(DMGetLabelName(dm, l, &lname));
1783: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
1784: if (isDepth) continue;
1785: PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
1786: if (isCellType) continue;
1787: PetscCall(DMCreateLabel(rdm, lname));
1788: PetscCall(DMGetLabel(dm, lname, &label));
1789: PetscCall(DMGetLabel(rdm, lname, &labelNew));
1790: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1791: }
1792: PetscFunctionReturn(PETSC_SUCCESS);
1793: }
1795: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1796: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1797: {
1798: DM dm;
1799: PetscInt Nf, f, Nds, s;
1801: PetscFunctionBegin;
1802: PetscCall(DMPlexTransformGetDM(tr, &dm));
1803: PetscCall(DMGetNumFields(dm, &Nf));
1804: for (f = 0; f < Nf; ++f) {
1805: DMLabel label, labelNew;
1806: PetscObject obj;
1807: const char *lname;
1809: PetscCall(DMGetField(rdm, f, &label, &obj));
1810: if (!label) continue;
1811: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1812: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1813: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1814: PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
1815: PetscCall(DMLabelDestroy(&labelNew));
1816: }
1817: PetscCall(DMGetNumDS(dm, &Nds));
1818: for (s = 0; s < Nds; ++s) {
1819: DMLabel label, labelNew;
1820: const char *lname;
1822: PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
1823: if (!label) continue;
1824: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1825: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1826: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1827: PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
1828: PetscCall(DMLabelDestroy(&labelNew));
1829: }
1830: PetscFunctionReturn(PETSC_SUCCESS);
1831: }
1833: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1834: {
1835: DM dm;
1836: PetscSF sf, sfNew;
1837: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
1838: const PetscInt *localPoints;
1839: const PetscSFNode *remotePoints;
1840: PetscInt *localPointsNew;
1841: PetscSFNode *remotePointsNew;
1842: PetscInt pStartNew, pEndNew, pNew;
1843: /* Brute force algorithm */
1844: PetscSF rsf;
1845: PetscSection s;
1846: const PetscInt *rootdegree;
1847: PetscInt *rootPointsNew, *remoteOffsets;
1848: PetscInt numPointsNew, pStart, pEnd, p;
1850: PetscFunctionBegin;
1851: PetscCall(DMPlexTransformGetDM(tr, &dm));
1852: PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
1853: PetscCall(DMGetPointSF(dm, &sf));
1854: PetscCall(DMGetPointSF(rdm, &sfNew));
1855: /* Calculate size of new SF */
1856: PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
1857: if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
1858: for (l = 0; l < numLeaves; ++l) {
1859: const PetscInt p = localPoints[l];
1860: DMPolytopeType ct;
1861: DMPolytopeType *rct;
1862: PetscInt *rsize, *rcone, *rornt;
1863: PetscInt Nct, n;
1865: PetscCall(DMPlexGetCellType(dm, p, &ct));
1866: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1867: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
1868: }
1869: /* Send new root point numbers
1870: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1871: */
1872: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1873: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
1874: PetscCall(PetscSectionSetChart(s, pStart, pEnd));
1875: for (p = pStart; p < pEnd; ++p) {
1876: DMPolytopeType ct;
1877: DMPolytopeType *rct;
1878: PetscInt *rsize, *rcone, *rornt;
1879: PetscInt Nct, n;
1881: PetscCall(DMPlexGetCellType(dm, p, &ct));
1882: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1883: for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
1884: }
1885: PetscCall(PetscSectionSetUp(s));
1886: PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
1887: PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
1888: PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
1889: PetscCall(PetscFree(remoteOffsets));
1890: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1891: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1892: PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
1893: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1894: for (p = pStart; p < pEnd; ++p) {
1895: DMPolytopeType ct;
1896: DMPolytopeType *rct;
1897: PetscInt *rsize, *rcone, *rornt;
1898: PetscInt Nct, n, r, off;
1900: if (!rootdegree[p - pStart]) continue;
1901: PetscCall(PetscSectionGetOffset(s, p, &off));
1902: PetscCall(DMPlexGetCellType(dm, p, &ct));
1903: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1904: for (n = 0, m = 0; n < Nct; ++n) {
1905: for (r = 0; r < rsize[n]; ++r, ++m) {
1906: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1907: rootPointsNew[off + m] = pNew;
1908: }
1909: }
1910: }
1911: PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1912: PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1913: PetscCall(PetscSFDestroy(&rsf));
1914: PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
1915: PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
1916: for (l = 0, m = 0; l < numLeaves; ++l) {
1917: const PetscInt p = localPoints[l];
1918: DMPolytopeType ct;
1919: DMPolytopeType *rct;
1920: PetscInt *rsize, *rcone, *rornt;
1921: PetscInt Nct, n, r, q, off;
1923: PetscCall(PetscSectionGetOffset(s, p, &off));
1924: PetscCall(DMPlexGetCellType(dm, p, &ct));
1925: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1926: for (n = 0, q = 0; n < Nct; ++n) {
1927: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1928: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1929: localPointsNew[m] = pNew;
1930: remotePointsNew[m].index = rootPointsNew[off + q];
1931: remotePointsNew[m].rank = remotePoints[l].rank;
1932: }
1933: }
1934: }
1935: PetscCall(PetscSectionDestroy(&s));
1936: PetscCall(PetscFree(rootPointsNew));
1937: /* SF needs sorted leaves to correctly calculate Gather */
1938: {
1939: PetscSFNode *rp, *rtmp;
1940: PetscInt *lp, *idx, *ltmp, i;
1942: PetscCall(PetscMalloc1(numLeavesNew, &idx));
1943: PetscCall(PetscMalloc1(numLeavesNew, &lp));
1944: PetscCall(PetscMalloc1(numLeavesNew, &rp));
1945: for (i = 0; i < numLeavesNew; ++i) {
1946: 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);
1947: idx[i] = i;
1948: }
1949: PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
1950: for (i = 0; i < numLeavesNew; ++i) {
1951: lp[i] = localPointsNew[idx[i]];
1952: rp[i] = remotePointsNew[idx[i]];
1953: }
1954: ltmp = localPointsNew;
1955: localPointsNew = lp;
1956: rtmp = remotePointsNew;
1957: remotePointsNew = rp;
1958: PetscCall(PetscFree(idx));
1959: PetscCall(PetscFree(ltmp));
1960: PetscCall(PetscFree(rtmp));
1961: }
1962: PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
1963: PetscFunctionReturn(PETSC_SUCCESS);
1964: }
1966: /*
1967: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`.
1969: Not Collective
1971: Input Parameters:
1972: + tr - The `DMPlexTransform`
1973: . ct - The type of the parent cell
1974: . rct - The type of the produced cell
1975: . r - The index of the produced cell
1976: - x - The localized coordinates for the parent cell
1978: Output Parameter:
1979: . xr - The localized coordinates for the produced cell
1981: Level: developer
1983: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
1984: */
1985: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1986: {
1987: PetscFE fe = NULL;
1988: PetscInt cdim, v, *subcellV;
1990: PetscFunctionBegin;
1991: PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
1992: PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
1993: PetscCall(PetscFEGetNumComponents(fe, &cdim));
1994: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
1995: PetscFunctionReturn(PETSC_SUCCESS);
1996: }
1998: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1999: {
2000: DM dm, cdm, cdmCell, cdmNew, cdmCellNew;
2001: PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
2002: Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
2003: const PetscScalar *coords;
2004: PetscScalar *coordsNew;
2005: const PetscReal *maxCell, *Lstart, *L;
2006: PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
2007: PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;
2009: PetscFunctionBegin;
2010: PetscCall(DMPlexTransformGetDM(tr, &dm));
2011: PetscCall(DMGetCoordinateDM(dm, &cdm));
2012: PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
2013: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
2014: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
2015: if (localized) {
2016: /* Localize coordinates of new vertices */
2017: localizeVertices = PETSC_TRUE;
2018: /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
2019: if (!maxCell) localizeCells = PETSC_TRUE;
2020: }
2021: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2022: PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
2023: if (maxCell) {
2024: PetscReal maxCellNew[3];
2026: for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
2027: PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L));
2028: }
2029: PetscCall(DMGetCoordinateDim(rdm, &dE));
2030: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
2031: PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
2032: PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
2033: PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
2034: PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
2035: /* Localization should be inherited */
2036: /* Stefano calculates parent cells for each new cell for localization */
2037: /* Localized cells need coordinates of closure */
2038: for (v = vStartNew; v < vEndNew; ++v) {
2039: PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
2040: PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
2041: }
2042: PetscCall(PetscSectionSetUp(coordSectionNew));
2043: PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));
2045: if (localizeCells) {
2046: PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
2047: PetscCall(DMClone(cdmNew, &cdmCellNew));
2048: PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
2049: PetscCall(DMDestroy(&cdmCellNew));
2051: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
2052: PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
2053: PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
2054: PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
2055: PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));
2057: PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
2058: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2059: for (c = cStart; c < cEnd; ++c) {
2060: PetscInt dof;
2062: PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
2063: if (dof) {
2064: DMPolytopeType ct;
2065: DMPolytopeType *rct;
2066: PetscInt *rsize, *rcone, *rornt;
2067: PetscInt dim, cNew, Nct, n, r;
2069: PetscCall(DMPlexGetCellType(dm, c, &ct));
2070: dim = DMPolytopeTypeGetDim(ct);
2071: PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2072: /* This allows for different cell types */
2073: for (n = 0; n < Nct; ++n) {
2074: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2075: for (r = 0; r < rsize[n]; ++r) {
2076: PetscInt *closure = NULL;
2077: PetscInt clSize, cl, Nv = 0;
2079: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
2080: PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2081: for (cl = 0; cl < clSize * 2; cl += 2) {
2082: if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
2083: }
2084: PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2085: PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
2086: PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
2087: }
2088: }
2089: }
2090: }
2091: PetscCall(PetscSectionSetUp(coordSectionCellNew));
2092: PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
2093: }
2094: PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
2095: {
2096: VecType vtype;
2097: PetscInt coordSizeNew, bs;
2098: const char *name;
2100: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
2101: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
2102: PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
2103: PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
2104: PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
2105: PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
2106: PetscCall(VecGetBlockSize(coordsLocal, &bs));
2107: PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
2108: PetscCall(VecGetType(coordsLocal, &vtype));
2109: PetscCall(VecSetType(coordsLocalNew, vtype));
2110: }
2111: PetscCall(VecGetArrayRead(coordsLocal, &coords));
2112: PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
2113: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2114: /* First set coordinates for vertices */
2115: for (p = pStart; p < pEnd; ++p) {
2116: DMPolytopeType ct;
2117: DMPolytopeType *rct;
2118: PetscInt *rsize, *rcone, *rornt;
2119: PetscInt Nct, n, r;
2120: PetscBool hasVertex = PETSC_FALSE;
2122: PetscCall(DMPlexGetCellType(dm, p, &ct));
2123: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2124: for (n = 0; n < Nct; ++n) {
2125: if (rct[n] == DM_POLYTOPE_POINT) {
2126: hasVertex = PETSC_TRUE;
2127: break;
2128: }
2129: }
2130: if (hasVertex) {
2131: const PetscScalar *icoords = NULL;
2132: const PetscScalar *array = NULL;
2133: PetscScalar *pcoords = NULL;
2134: PetscBool isDG;
2135: PetscInt Nc, Nv, v, d;
2137: PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2139: icoords = pcoords;
2140: Nv = Nc / dEo;
2141: if (ct != DM_POLYTOPE_POINT) {
2142: if (localizeVertices && maxCell) {
2143: PetscScalar anchor[3];
2145: for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
2146: for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
2147: }
2148: }
2149: for (n = 0; n < Nct; ++n) {
2150: if (rct[n] != DM_POLYTOPE_POINT) continue;
2151: for (r = 0; r < rsize[n]; ++r) {
2152: PetscScalar vcoords[3];
2153: PetscInt vNew, off;
2155: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
2156: PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
2157: PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
2158: PetscCall(DMSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
2159: }
2160: }
2161: PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2162: }
2163: }
2164: PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
2165: PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
2166: PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
2167: PetscCall(VecDestroy(&coordsLocalNew));
2168: PetscCall(PetscSectionDestroy(&coordSectionNew));
2169: /* Then set coordinates for cells by localizing */
2170: if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
2171: else {
2172: VecType vtype;
2173: PetscInt coordSizeNew, bs;
2174: const char *name;
2176: PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
2177: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
2178: PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
2179: PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
2180: PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
2181: PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
2182: PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
2183: PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
2184: PetscCall(VecGetType(coordsLocalCell, &vtype));
2185: PetscCall(VecSetType(coordsLocalCellNew, vtype));
2186: PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
2187: PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));
2189: for (p = pStart; p < pEnd; ++p) {
2190: DMPolytopeType ct;
2191: DMPolytopeType *rct;
2192: PetscInt *rsize, *rcone, *rornt;
2193: PetscInt dof = 0, Nct, n, r;
2195: PetscCall(DMPlexGetCellType(dm, p, &ct));
2196: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2197: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
2198: if (dof) {
2199: const PetscScalar *pcoords;
2201: PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords));
2202: for (n = 0; n < Nct; ++n) {
2203: const PetscInt Nr = rsize[n];
2205: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2206: for (r = 0; r < Nr; ++r) {
2207: PetscInt pNew, offNew;
2209: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2210: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2211: cell to the ones it produces. */
2212: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2213: PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
2214: PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
2215: }
2216: }
2217: }
2218: }
2219: PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
2220: PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
2221: PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
2222: PetscCall(VecDestroy(&coordsLocalCellNew));
2223: PetscCall(PetscSectionDestroy(&coordSectionCellNew));
2224: }
2225: PetscFunctionReturn(PETSC_SUCCESS);
2226: }
2228: /*@
2229: DMPlexTransformApply - Execute the transformation, producing another `DM`
2231: Collective
2233: Input Parameters:
2234: + tr - The `DMPlexTransform` object
2235: - dm - The original `DM`
2237: Output Parameter:
2238: . tdm - The transformed `DM`
2240: Level: intermediate
2242: Options Database Keys:
2243: + -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point
2244: . -dm_plex_transform_label_replica_inc <num> - Increment for the label value to be multiplied by the replica number
2245: - -dm_plex_transform_active <name> - Name for active mesh label
2247: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()`
2248: @*/
2249: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
2250: {
2251: DM rdm;
2252: DMPlexInterpolatedFlag interp;
2253: PetscInt pStart, pEnd;
2255: PetscFunctionBegin;
2258: PetscAssertPointer(tdm, 3);
2259: PetscCall(PetscLogEventBegin(DMPLEX_Transform, tr, dm, 0, 0));
2260: PetscCall(DMPlexTransformSetDM(tr, dm));
2262: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
2263: PetscCall(DMSetType(rdm, DMPLEX));
2264: PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
2265: /* Calculate number of new points of each depth */
2266: PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
2267: PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
2268: /* Step 1: Set chart */
2269: PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
2270: PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
2271: /* Step 2: Set cone/support sizes (automatically stratifies) */
2272: PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
2273: /* Step 3: Setup refined DM */
2274: PetscCall(DMSetUp(rdm));
2275: /* Step 4: Set cones and supports (automatically symmetrizes) */
2276: PetscCall(DMPlexTransformSetCones(tr, rdm));
2277: /* Step 5: Create pointSF */
2278: PetscCall(DMPlexTransformCreateSF(tr, rdm));
2279: /* Step 6: Create labels */
2280: PetscCall(DMPlexTransformCreateLabels(tr, rdm));
2281: /* Step 7: Set coordinates */
2282: PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
2283: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm));
2284: // If the original DM was configured from options, the transformed DM should be as well
2285: rdm->setfromoptionscalled = dm->setfromoptionscalled;
2286: PetscCall(PetscLogEventEnd(DMPLEX_Transform, tr, dm, 0, 0));
2287: *tdm = rdm;
2288: PetscFunctionReturn(PETSC_SUCCESS);
2289: }
2291: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2292: {
2293: DMPlexTransform tr;
2294: DM cdm, rcdm;
2295: const char *prefix;
2297: PetscFunctionBegin;
2298: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
2299: PetscCall(PetscObjectSetName((PetscObject)tr, "Adapt Label Transform"));
2300: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2301: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
2302: PetscCall(DMPlexTransformSetDM(tr, dm));
2303: PetscCall(DMPlexTransformSetFromOptions(tr));
2304: PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
2305: PetscCall(DMPlexTransformSetUp(tr));
2306: PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
2307: PetscCall(DMPlexTransformApply(tr, dm, rdm));
2308: PetscCall(DMCopyDisc(dm, *rdm));
2309: PetscCall(DMGetCoordinateDM(dm, &cdm));
2310: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
2311: PetscCall(DMCopyDisc(cdm, rcdm));
2312: PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
2313: PetscCall(DMCopyDisc(dm, *rdm));
2314: PetscCall(DMPlexTransformDestroy(&tr));
2315: ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2316: PetscFunctionReturn(PETSC_SUCCESS);
2317: }