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: /*@C
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];
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: PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
341: /* process any options handlers added with PetscObjectAddOptionsHandler() */
342: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject));
343: PetscOptionsEnd();
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: /*@
348: DMPlexTransformDestroy - Destroys a `DMPlexTransform`
350: Collective
352: Input Parameter:
353: . tr - the transform object to destroy
355: Level: beginner
357: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
358: @*/
359: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
360: {
361: PetscInt c;
363: PetscFunctionBegin;
364: if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
366: if (--((PetscObject)*tr)->refct > 0) {
367: *tr = NULL;
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: PetscTryTypeMethod(*tr, destroy);
372: PetscCall(DMDestroy(&(*tr)->dm));
373: PetscCall(DMLabelDestroy(&(*tr)->active));
374: PetscCall(DMLabelDestroy(&(*tr)->trType));
375: PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld));
376: PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew));
377: PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew));
378: PetscCall(PetscFree((*tr)->offset));
379: PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd));
380: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
381: PetscCall(PetscFEDestroy(&(*tr)->coordFE[c]));
382: PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c]));
383: }
384: if ((*tr)->trVerts) {
385: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
386: DMPolytopeType *rct;
387: PetscInt *rsize, *rcone, *rornt, Nct, n, r;
389: if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) {
390: PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
391: for (n = 0; n < Nct; ++n) {
392: if (rct[n] == DM_POLYTOPE_POINT) continue;
393: for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r]));
394: PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]]));
395: }
396: }
397: PetscCall(PetscFree((*tr)->trSubVerts[c]));
398: PetscCall(PetscFree((*tr)->trVerts[c]));
399: }
400: }
401: PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts));
402: PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom));
403: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
404: PetscCall(PetscHeaderDestroy(tr));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
409: {
410: DMLabel trType = tr->trType;
411: PetscInt c, cN, *off;
413: PetscFunctionBegin;
414: if (trType) {
415: DM dm;
416: IS rtIS;
417: const PetscInt *reftypes;
418: PetscInt Nrt, r;
420: PetscCall(DMPlexTransformGetDM(tr, &dm));
421: PetscCall(DMLabelGetNumValues(trType, &Nrt));
422: PetscCall(DMLabelGetValueIS(trType, &rtIS));
423: PetscCall(ISGetIndices(rtIS, &reftypes));
424: PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off));
425: for (r = 0; r < Nrt; ++r) {
426: const PetscInt rt = reftypes[r];
427: IS rtIS;
428: const PetscInt *points;
429: DMPolytopeType ct;
430: PetscInt np, p;
432: PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
433: PetscCall(ISGetLocalSize(rtIS, &np));
434: PetscCall(ISGetIndices(rtIS, &points));
435: if (!np) continue;
436: p = points[0];
437: PetscCall(ISRestoreIndices(rtIS, &points));
438: PetscCall(ISDestroy(&rtIS));
439: PetscCall(DMPlexGetCellType(dm, p, &ct));
440: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
441: const DMPolytopeType ctNew = (DMPolytopeType)cN;
442: DMPolytopeType *rct;
443: PetscInt *rsize, *cone, *ornt;
444: PetscInt Nct, n, s;
446: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
447: off[r * DM_NUM_POLYTOPES + ctNew] = -1;
448: break;
449: }
450: off[r * DM_NUM_POLYTOPES + ctNew] = 0;
451: for (s = 0; s <= r; ++s) {
452: const PetscInt st = reftypes[s];
453: DMPolytopeType sct;
454: PetscInt q, qrt;
456: PetscCall(DMLabelGetStratumIS(trType, st, &rtIS));
457: PetscCall(ISGetLocalSize(rtIS, &np));
458: PetscCall(ISGetIndices(rtIS, &points));
459: if (!np) continue;
460: q = points[0];
461: PetscCall(ISRestoreIndices(rtIS, &points));
462: PetscCall(ISDestroy(&rtIS));
463: PetscCall(DMPlexGetCellType(dm, q, &sct));
464: PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt));
465: 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);
466: if (st == rt) {
467: for (n = 0; n < Nct; ++n)
468: if (rct[n] == ctNew) break;
469: if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
470: break;
471: }
472: for (n = 0; n < Nct; ++n) {
473: if (rct[n] == ctNew) {
474: PetscInt sn;
476: PetscCall(DMLabelGetStratumSize(trType, st, &sn));
477: off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
478: }
479: }
480: }
481: }
482: }
483: PetscCall(ISRestoreIndices(rtIS, &reftypes));
484: PetscCall(ISDestroy(&rtIS));
485: } else {
486: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off));
487: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
488: const DMPolytopeType ct = (DMPolytopeType)c;
489: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
490: const DMPolytopeType ctNew = (DMPolytopeType)cN;
491: DMPolytopeType *rct;
492: PetscInt *rsize, *cone, *ornt;
493: PetscInt Nct, n, i;
495: 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) {
496: off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
497: continue;
498: }
499: off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
500: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
501: const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i];
502: const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];
504: PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt));
505: if (ict == ct) {
506: for (n = 0; n < Nct; ++n)
507: if (rct[n] == ctNew) break;
508: if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
509: break;
510: }
511: for (n = 0; n < Nct; ++n)
512: if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
513: }
514: }
515: }
516: }
517: *offset = off;
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: /*@
522: DMPlexTransformSetUp - Create the tables that drive the transform
524: Input Parameter:
525: . tr - The `DMPlexTransform` object
527: Level: intermediate
529: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
530: @*/
531: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
532: {
533: DM dm;
534: DMPolytopeType ctCell;
535: PetscInt pStart, pEnd, p, c, celldim = 0;
537: PetscFunctionBegin;
539: if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
540: PetscTryTypeMethod(tr, setup);
541: PetscCall(DMPlexTransformGetDM(tr, &dm));
542: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
543: if (pEnd > pStart) {
544: PetscCall(DMPlexGetCellType(dm, 0, &ctCell));
545: } else {
546: PetscInt dim;
548: PetscCall(DMGetDimension(dm, &dim));
549: switch (dim) {
550: case 0:
551: ctCell = DM_POLYTOPE_POINT;
552: break;
553: case 1:
554: ctCell = DM_POLYTOPE_SEGMENT;
555: break;
556: case 2:
557: ctCell = DM_POLYTOPE_TRIANGLE;
558: break;
559: case 3:
560: ctCell = DM_POLYTOPE_TETRAHEDRON;
561: break;
562: default:
563: ctCell = DM_POLYTOPE_UNKNOWN;
564: }
565: }
566: PetscCall(DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
567: for (p = pStart; p < pEnd; ++p) {
568: DMPolytopeType ct;
569: DMPolytopeType *rct;
570: PetscInt *rsize, *cone, *ornt;
571: PetscInt Nct, n;
573: PetscCall(DMPlexGetCellType(dm, p, &ct));
574: 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);
575: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
576: for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
577: }
578: PetscCall(DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
579: /* Construct sizes and offsets for each cell type */
580: if (!tr->ctStart) {
581: PetscInt *ctS, *ctSN, *ctC, *ctCN;
583: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
584: PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
585: for (p = pStart; p < pEnd; ++p) {
586: DMPolytopeType ct;
587: DMPolytopeType *rct;
588: PetscInt *rsize, *cone, *ornt;
589: PetscInt Nct, n;
591: PetscCall(DMPlexGetCellType(dm, p, &ct));
592: 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);
593: ++ctC[ct];
594: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
595: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
596: }
597: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
598: const PetscInt cto = tr->ctOrderOld[c];
599: const PetscInt cton = tr->ctOrderOld[c + 1];
600: const PetscInt ctn = tr->ctOrderNew[c];
601: const PetscInt ctnn = tr->ctOrderNew[c + 1];
603: ctS[cton] = ctS[cto] + ctC[cto];
604: ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
605: }
606: PetscCall(PetscFree2(ctC, ctCN));
607: tr->ctStart = ctS;
608: tr->ctStartNew = ctSN;
609: }
610: PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
611: // Compute depth information
612: tr->depth = -1;
613: for (c = 0; c < DM_NUM_POLYTOPES; ++c)
614: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
615: PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
616: for (PetscInt d = 0; d <= tr->depth; ++d) {
617: tr->depthStart[d] = PETSC_MAX_INT;
618: tr->depthEnd[d] = -1;
619: }
620: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
621: const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);
623: if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
624: tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
625: tr->depthEnd[dep] = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
626: }
627: tr->setupcalled = PETSC_TRUE;
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: /*@
632: DMPlexTransformGetDM - Get the base `DM` for the transform
634: Input Parameter:
635: . tr - The `DMPlexTransform` object
637: Output Parameter:
638: . dm - The original `DM` which will be transformed
640: Level: intermediate
642: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
643: @*/
644: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
645: {
646: PetscFunctionBegin;
648: PetscAssertPointer(dm, 2);
649: *dm = tr->dm;
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: /*@
654: DMPlexTransformSetDM - Set the base `DM` for the transform
656: Input Parameters:
657: + tr - The `DMPlexTransform` object
658: - dm - The original `DM` which will be transformed
660: Level: intermediate
662: Note:
663: The user does not typically call this, as it is called by `DMPlexTransformApply()`.
665: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
666: @*/
667: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
668: {
669: PetscFunctionBegin;
672: PetscCall(PetscObjectReference((PetscObject)dm));
673: PetscCall(DMDestroy(&tr->dm));
674: tr->dm = dm;
675: PetscFunctionReturn(PETSC_SUCCESS);
676: }
678: /*@
679: DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform
681: Input Parameter:
682: . tr - The `DMPlexTransform` object
684: Output Parameter:
685: . active - The `DMLabel` indicating which points will be transformed
687: Level: intermediate
689: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
690: @*/
691: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
692: {
693: PetscFunctionBegin;
695: PetscAssertPointer(active, 2);
696: *active = tr->active;
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: /*@
701: DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform
703: Input Parameters:
704: + tr - The `DMPlexTransform` object
705: - active - The original `DM` which will be transformed
707: Level: intermediate
709: Note:
710: This only applies to transforms that can operator on a subset of the mesh, listed in [](plex_transform_table).
712: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
713: @*/
714: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
715: {
716: PetscFunctionBegin;
719: PetscCall(PetscObjectReference((PetscObject)active));
720: PetscCall(DMLabelDestroy(&tr->active));
721: tr->active = active;
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
726: {
727: PetscFunctionBegin;
728: if (!tr->coordFE[ct]) {
729: PetscInt dim, cdim;
731: dim = DMPolytopeTypeGetDim(ct);
732: PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
733: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
734: {
735: PetscDualSpace dsp;
736: PetscQuadrature quad;
737: DM K;
738: PetscFEGeom *cg;
739: PetscScalar *Xq;
740: PetscReal *xq, *wq;
741: PetscInt Nq, q;
743: PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
744: PetscCall(PetscMalloc1(Nq * cdim, &xq));
745: for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
746: PetscCall(PetscMalloc1(Nq, &wq));
747: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
748: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
749: PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
750: PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));
752: PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
753: PetscCall(PetscDualSpaceGetDM(dsp, &K));
754: PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]));
755: cg = tr->refGeom[ct];
756: PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
757: PetscCall(PetscQuadratureDestroy(&quad));
758: }
759: }
760: *fe = tr->coordFE[ct];
761: PetscFunctionReturn(PETSC_SUCCESS);
762: }
764: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
765: {
766: PetscInt dim, cdim;
768: PetscFunctionBegin;
769: PetscCall(DMGetDimension(dm, &dim));
770: PetscCall(DMSetDimension(tdm, dim));
771: PetscCall(DMGetCoordinateDim(dm, &cdim));
772: PetscCall(DMSetCoordinateDim(tdm, cdim));
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: /*@
777: DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`
779: Input Parameters:
780: + tr - The `DMPlexTransform` object
781: - dm - The original `DM`
783: Output Parameter:
784: . tdm - The transformed `DM`
786: Level: advanced
788: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
789: @*/
790: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
791: {
792: PetscFunctionBegin;
793: PetscUseTypeMethod(tr, setdimensions, dm, tdm);
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
798: {
799: PetscFunctionBegin;
800: if (pStart) *pStart = 0;
801: if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
802: PetscFunctionReturn(PETSC_SUCCESS);
803: }
805: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
806: {
807: PetscInt ctNew;
809: PetscFunctionBegin;
811: PetscAssertPointer(celltype, 3);
812: /* TODO Can do bisection since everything is sorted */
813: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
814: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
816: if (cell >= ctSN && cell < ctEN) break;
817: }
818: PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
819: *celltype = (DMPolytopeType)ctNew;
820: PetscFunctionReturn(PETSC_SUCCESS);
821: }
823: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
824: {
825: PetscFunctionBegin;
827: if (start) *start = tr->ctStartNew[celltype];
828: if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
833: {
834: PetscFunctionBegin;
836: *depth = tr->depth;
837: PetscFunctionReturn(PETSC_SUCCESS);
838: }
840: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
841: {
842: PetscFunctionBegin;
844: if (start) *start = tr->depthStart[depth];
845: if (end) *end = tr->depthEnd[depth];
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: /*@
850: DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.
852: Not Collective
854: Input Parameters:
855: + tr - The `DMPlexTransform`
856: . ct - The type of the original point which produces the new point
857: . ctNew - The type of the new point
858: . p - The original point which produces the new point
859: - r - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`
861: Output Parameter:
862: . pNew - The new point number
864: Level: developer
866: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
867: @*/
868: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
869: {
870: DMPolytopeType *rct;
871: PetscInt *rsize, *cone, *ornt;
872: PetscInt rt, Nct, n, off, rp;
873: DMLabel trType = tr->trType;
874: PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
875: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
876: PetscInt newp = ctSN, cind;
878: PetscFunctionBeginHot;
879: 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);
880: PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
881: if (trType) {
882: PetscCall(DMLabelGetValueIndex(trType, rt, &cind));
883: PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
884: 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);
885: } else {
886: cind = ct;
887: rp = p - ctS;
888: }
889: off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
890: 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);
891: newp += off;
892: for (n = 0; n < Nct; ++n) {
893: if (rct[n] == ctNew) {
894: if (rsize[n] && r >= rsize[n])
895: 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]);
896: newp += rp * rsize[n] + r;
897: break;
898: }
899: }
901: 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);
902: *pNew = newp;
903: PetscFunctionReturn(PETSC_SUCCESS);
904: }
906: /*@
907: DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.
909: Not Collective
911: Input Parameters:
912: + tr - The `DMPlexTransform`
913: - pNew - The new point number
915: Output Parameters:
916: + ct - The type of the original point which produces the new point
917: . ctNew - The type of the new point
918: . p - The original point which produces the new point
919: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
921: Level: developer
923: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
924: @*/
925: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
926: {
927: DMLabel trType = tr->trType;
928: DMPolytopeType *rct, ctN;
929: PetscInt *rsize, *cone, *ornt;
930: PetscInt rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
931: PetscInt offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;
933: PetscFunctionBegin;
934: PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
935: if (trType) {
936: DM dm;
937: IS rtIS;
938: const PetscInt *reftypes;
939: PetscInt Nrt, r, rtStart;
941: PetscCall(DMPlexTransformGetDM(tr, &dm));
942: PetscCall(DMLabelGetNumValues(trType, &Nrt));
943: PetscCall(DMLabelGetValueIS(trType, &rtIS));
944: PetscCall(ISGetIndices(rtIS, &reftypes));
945: for (r = 0; r < Nrt; ++r) {
946: const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];
948: if (tr->ctStartNew[ctN] + off > pNew) continue;
949: /* Check that any of this refinement type exist */
950: /* TODO Actually keep track of the number produced here instead */
951: if (off > offset) {
952: rt = reftypes[r];
953: offset = off;
954: }
955: }
956: PetscCall(ISRestoreIndices(rtIS, &reftypes));
957: PetscCall(ISDestroy(&rtIS));
958: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
959: /* TODO Map refinement types to cell types */
960: PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL));
961: PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
962: for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
963: PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
965: if ((rtStart >= ctS) && (rtStart < ctE)) break;
966: }
967: PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt);
968: } else {
969: for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
970: const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];
972: if (tr->ctStartNew[ctN] + off > pNew) continue;
973: if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
974: /* TODO Actually keep track of the number produced here instead */
975: if (off > offset) {
976: ctO = ctTmp;
977: offset = off;
978: }
979: }
980: rt = -1;
981: PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
982: }
983: ctS = tr->ctStart[ctO];
984: ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
985: if (trType) {
986: for (rtS = ctS; rtS < ctE; ++rtS) {
987: PetscInt val;
988: PetscCall(DMLabelGetValue(trType, rtS, &val));
989: if (val == rt) break;
990: }
991: 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);
992: } else rtS = ctS;
993: PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
994: 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);
995: for (n = 0; n < Nct; ++n) {
996: if (rct[n] == ctN) {
997: PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c;
999: if (trType) {
1000: for (c = ctS; c < ctE; ++c) {
1001: PetscCall(DMLabelGetValue(trType, c, &val));
1002: if (val == rt) {
1003: if (tmp < rsize[n]) break;
1004: tmp -= rsize[n];
1005: }
1006: }
1007: PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew);
1008: rp = c - ctS;
1009: rO = tmp;
1010: } else {
1011: // This assumes that all points of type ctO transform the same way
1012: rp = tmp / rsize[n];
1013: rO = tmp % rsize[n];
1014: }
1015: break;
1016: }
1017: }
1018: PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
1019: pO = rp + ctS;
1020: 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);
1021: if (ct) *ct = (DMPolytopeType)ctO;
1022: if (ctNew) *ctNew = (DMPolytopeType)ctN;
1023: if (p) *p = pO;
1024: if (r) *r = rO;
1025: PetscFunctionReturn(PETSC_SUCCESS);
1026: }
1028: /*@
1029: DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.
1031: Input Parameters:
1032: + tr - The `DMPlexTransform` object
1033: . source - The source cell type
1034: - p - The source point, which can also determine the refine type
1036: Output Parameters:
1037: + rt - The refine type for this point
1038: . Nt - The number of types produced by this point
1039: . target - An array of length `Nt` giving the types produced
1040: . size - An array of length `Nt` giving the number of cells of each type produced
1041: . cone - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
1042: - ornt - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point
1044: Level: advanced
1046: Notes:
1047: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
1048: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
1049: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
1050: .vb
1051: the number of cones to be taken, or 0 for the current cell
1052: the cell cone point number at each level from which it is subdivided
1053: the replica number r of the subdivision.
1054: .ve
1055: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
1056: .vb
1057: Nt = 2
1058: target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
1059: size = {1, 2}
1060: cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
1061: ornt = { 0, 0, 0, 0}
1062: .ve
1064: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
1065: @*/
1066: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1067: {
1068: PetscFunctionBegin;
1069: PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
1070: PetscFunctionReturn(PETSC_SUCCESS);
1071: }
1073: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1074: {
1075: PetscFunctionBegin;
1076: *rnew = r;
1077: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
1078: PetscFunctionReturn(PETSC_SUCCESS);
1079: }
1081: /* Returns the same thing */
1082: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1083: {
1084: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
1085: static PetscInt vertexS[] = {1};
1086: static PetscInt vertexC[] = {0};
1087: static PetscInt vertexO[] = {0};
1088: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
1089: static PetscInt edgeS[] = {1};
1090: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1091: static PetscInt edgeO[] = {0, 0};
1092: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
1093: static PetscInt tedgeS[] = {1};
1094: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1095: static PetscInt tedgeO[] = {0, 0};
1096: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
1097: static PetscInt triS[] = {1};
1098: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1099: static PetscInt triO[] = {0, 0, 0};
1100: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
1101: static PetscInt quadS[] = {1};
1102: 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};
1103: static PetscInt quadO[] = {0, 0, 0, 0};
1104: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
1105: static PetscInt tquadS[] = {1};
1106: 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};
1107: static PetscInt tquadO[] = {0, 0, 0, 0};
1108: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
1109: static PetscInt tetS[] = {1};
1110: 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};
1111: static PetscInt tetO[] = {0, 0, 0, 0};
1112: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
1113: static PetscInt hexS[] = {1};
1114: 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};
1115: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
1116: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
1117: static PetscInt tripS[] = {1};
1118: 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};
1119: static PetscInt tripO[] = {0, 0, 0, 0, 0};
1120: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
1121: static PetscInt ttripS[] = {1};
1122: 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};
1123: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
1124: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
1125: static PetscInt tquadpS[] = {1};
1126: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
1127: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
1128: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
1129: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
1130: static PetscInt pyrS[] = {1};
1131: 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};
1132: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
1134: PetscFunctionBegin;
1135: if (rt) *rt = 0;
1136: switch (source) {
1137: case DM_POLYTOPE_POINT:
1138: *Nt = 1;
1139: *target = vertexT;
1140: *size = vertexS;
1141: *cone = vertexC;
1142: *ornt = vertexO;
1143: break;
1144: case DM_POLYTOPE_SEGMENT:
1145: *Nt = 1;
1146: *target = edgeT;
1147: *size = edgeS;
1148: *cone = edgeC;
1149: *ornt = edgeO;
1150: break;
1151: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1152: *Nt = 1;
1153: *target = tedgeT;
1154: *size = tedgeS;
1155: *cone = tedgeC;
1156: *ornt = tedgeO;
1157: break;
1158: case DM_POLYTOPE_TRIANGLE:
1159: *Nt = 1;
1160: *target = triT;
1161: *size = triS;
1162: *cone = triC;
1163: *ornt = triO;
1164: break;
1165: case DM_POLYTOPE_QUADRILATERAL:
1166: *Nt = 1;
1167: *target = quadT;
1168: *size = quadS;
1169: *cone = quadC;
1170: *ornt = quadO;
1171: break;
1172: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1173: *Nt = 1;
1174: *target = tquadT;
1175: *size = tquadS;
1176: *cone = tquadC;
1177: *ornt = tquadO;
1178: break;
1179: case DM_POLYTOPE_TETRAHEDRON:
1180: *Nt = 1;
1181: *target = tetT;
1182: *size = tetS;
1183: *cone = tetC;
1184: *ornt = tetO;
1185: break;
1186: case DM_POLYTOPE_HEXAHEDRON:
1187: *Nt = 1;
1188: *target = hexT;
1189: *size = hexS;
1190: *cone = hexC;
1191: *ornt = hexO;
1192: break;
1193: case DM_POLYTOPE_TRI_PRISM:
1194: *Nt = 1;
1195: *target = tripT;
1196: *size = tripS;
1197: *cone = tripC;
1198: *ornt = tripO;
1199: break;
1200: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1201: *Nt = 1;
1202: *target = ttripT;
1203: *size = ttripS;
1204: *cone = ttripC;
1205: *ornt = ttripO;
1206: break;
1207: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1208: *Nt = 1;
1209: *target = tquadpT;
1210: *size = tquadpS;
1211: *cone = tquadpC;
1212: *ornt = tquadpO;
1213: break;
1214: case DM_POLYTOPE_PYRAMID:
1215: *Nt = 1;
1216: *target = pyrT;
1217: *size = pyrS;
1218: *cone = pyrC;
1219: *ornt = pyrO;
1220: break;
1221: default:
1222: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1223: }
1224: PetscFunctionReturn(PETSC_SUCCESS);
1225: }
1227: /*@
1228: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
1230: Not Collective
1232: Input Parameters:
1233: + tr - The `DMPlexTransform`
1234: . sct - The source point cell type, from whom the new cell is being produced
1235: . sp - The source point
1236: . so - The orientation of the source point in its enclosing parent
1237: . tct - The target point cell type
1238: . r - The replica number requested for the produced cell type
1239: - o - The orientation of the replica
1241: Output Parameters:
1242: + rnew - The replica number, given the orientation of the parent
1243: - onew - The replica orientation, given the orientation of the parent
1245: Level: advanced
1247: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1248: @*/
1249: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1250: {
1251: PetscFunctionBeginHot;
1252: PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1257: {
1258: DM dm;
1259: PetscInt pStart, pEnd, p, pNew;
1261: PetscFunctionBegin;
1262: PetscCall(DMPlexTransformGetDM(tr, &dm));
1263: /* Must create the celltype label here so that we do not automatically try to compute the types */
1264: PetscCall(DMCreateLabel(rdm, "celltype"));
1265: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1266: for (p = pStart; p < pEnd; ++p) {
1267: DMPolytopeType ct;
1268: DMPolytopeType *rct;
1269: PetscInt *rsize, *rcone, *rornt;
1270: PetscInt Nct, n, r;
1272: PetscCall(DMPlexGetCellType(dm, p, &ct));
1273: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1274: for (n = 0; n < Nct; ++n) {
1275: for (r = 0; r < rsize[n]; ++r) {
1276: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1277: PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
1278: PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
1279: }
1280: }
1281: }
1282: /* Let the DM know we have set all the cell types */
1283: {
1284: DMLabel ctLabel;
1285: DM_Plex *plex = (DM_Plex *)rdm->data;
1287: PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
1288: PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
1289: }
1290: PetscFunctionReturn(PETSC_SUCCESS);
1291: }
1293: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1294: {
1295: DMPolytopeType ctNew;
1297: PetscFunctionBegin;
1299: PetscAssertPointer(coneSize, 3);
1300: PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
1301: *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType)ctNew);
1302: PetscFunctionReturn(PETSC_SUCCESS);
1303: }
1305: /* The orientation o is for the interior of the cell p */
1306: 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[])
1307: {
1308: DM dm;
1309: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1310: const PetscInt *cone;
1311: PetscInt c, coff = *coneoff, ooff = *orntoff;
1313: PetscFunctionBegin;
1314: PetscCall(DMPlexTransformGetDM(tr, &dm));
1315: PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
1316: for (c = 0; c < csizeNew; ++c) {
1317: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1318: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1319: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1320: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1321: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1322: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1323: const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1324: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1325: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1326: PetscInt lc;
1328: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1329: for (lc = 0; lc < fn; ++lc) {
1330: const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po);
1331: const PetscInt acp = rcone[coff++];
1332: const PetscInt pcp = parr[acp * 2];
1333: const PetscInt pco = parr[acp * 2 + 1];
1334: const PetscInt *ppornt;
1336: ppp = pp;
1337: pp = pcone[pcp];
1338: PetscCall(DMPlexGetCellType(dm, pp, &pct));
1339: // Restore the parent cone from the last iterate
1340: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
1341: PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
1342: PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
1343: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1344: PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
1345: }
1346: if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
1347: pr = rcone[coff++];
1348: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1349: PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
1350: PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
1351: orntNew[c] = fo;
1352: }
1353: PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
1354: *coneoff = coff;
1355: *orntoff = ooff;
1356: PetscFunctionReturn(PETSC_SUCCESS);
1357: }
1359: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1360: {
1361: DM dm;
1362: DMPolytopeType ct;
1363: PetscInt *coneNew, *orntNew;
1364: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1366: PetscFunctionBegin;
1367: PetscCall(DMPlexTransformGetDM(tr, &dm));
1368: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1369: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1370: PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1371: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1372: for (p = pStart; p < pEnd; ++p) {
1373: PetscInt coff, ooff;
1374: DMPolytopeType *rct;
1375: PetscInt *rsize, *rcone, *rornt;
1376: PetscInt Nct, n, r;
1378: PetscCall(DMPlexGetCellType(dm, p, &ct));
1379: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1380: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1381: const DMPolytopeType ctNew = rct[n];
1383: for (r = 0; r < rsize[n]; ++r) {
1384: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1385: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
1386: PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
1387: PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
1388: }
1389: }
1390: }
1391: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
1392: PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
1393: PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
1394: PetscCall(DMPlexSymmetrize(rdm));
1395: PetscCall(DMPlexStratify(rdm));
1396: PetscFunctionReturn(PETSC_SUCCESS);
1397: }
1399: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1400: {
1401: DM dm;
1402: DMPolytopeType ct, qct;
1403: DMPolytopeType *rct;
1404: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1405: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1407: PetscFunctionBegin;
1409: PetscAssertPointer(cone, 4);
1410: PetscAssertPointer(ornt, 5);
1411: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1412: PetscCall(DMPlexTransformGetDM(tr, &dm));
1413: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1414: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1415: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1416: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1417: for (n = 0; n < Nct; ++n) {
1418: const DMPolytopeType ctNew = rct[n];
1419: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1420: PetscInt Nr = rsize[n], fn, c;
1422: if (ctNew == qct) Nr = r;
1423: for (nr = 0; nr < Nr; ++nr) {
1424: for (c = 0; c < csizeNew; ++c) {
1425: ++coff; /* Cell type of new cone point */
1426: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1427: coff += fn;
1428: ++coff; /* Replica number of new cone point */
1429: ++ooff; /* Orientation of new cone point */
1430: }
1431: }
1432: if (ctNew == qct) break;
1433: }
1434: PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1435: *cone = qcone;
1436: *ornt = qornt;
1437: PetscFunctionReturn(PETSC_SUCCESS);
1438: }
1440: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1441: {
1442: DM dm;
1443: DMPolytopeType ct, qct;
1444: DMPolytopeType *rct;
1445: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1446: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1448: PetscFunctionBegin;
1450: if (cone) PetscAssertPointer(cone, 3);
1451: if (ornt) PetscAssertPointer(ornt, 4);
1452: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1453: PetscCall(DMPlexTransformGetDM(tr, &dm));
1454: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1455: PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1456: PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
1457: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1458: for (n = 0; n < Nct; ++n) {
1459: const DMPolytopeType ctNew = rct[n];
1460: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1461: PetscInt Nr = rsize[n], fn, c;
1463: if (ctNew == qct) Nr = r;
1464: for (nr = 0; nr < Nr; ++nr) {
1465: for (c = 0; c < csizeNew; ++c) {
1466: ++coff; /* Cell type of new cone point */
1467: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1468: coff += fn;
1469: ++coff; /* Replica number of new cone point */
1470: ++ooff; /* Orientation of new cone point */
1471: }
1472: }
1473: if (ctNew == qct) break;
1474: }
1475: PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
1476: if (cone) *cone = qcone;
1477: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
1478: if (ornt) *ornt = qornt;
1479: else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
1480: PetscFunctionReturn(PETSC_SUCCESS);
1481: }
1483: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1484: {
1485: DM dm;
1487: PetscFunctionBegin;
1489: PetscCall(DMPlexTransformGetDM(tr, &dm));
1490: if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
1491: if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
1492: PetscFunctionReturn(PETSC_SUCCESS);
1493: }
1495: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1496: {
1497: PetscInt ict;
1499: PetscFunctionBegin;
1500: PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
1501: for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1502: const DMPolytopeType ct = (DMPolytopeType)ict;
1503: DMPlexTransform reftr;
1504: DM refdm, trdm;
1505: Vec coordinates;
1506: const PetscScalar *coords;
1507: DMPolytopeType *rct;
1508: PetscInt *rsize, *rcone, *rornt;
1509: PetscInt Nct, n, r, pNew;
1510: PetscInt trdim, vStart, vEnd, Nc;
1511: const PetscInt debug = 0;
1512: const char *typeName;
1514: /* Since points are 0-dimensional, coordinates make no sense */
1515: if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue;
1516: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
1517: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
1518: PetscCall(DMPlexTransformSetDM(reftr, refdm));
1519: PetscCall(DMPlexTransformGetType(tr, &typeName));
1520: PetscCall(DMPlexTransformSetType(reftr, typeName));
1521: PetscCall(DMPlexTransformSetUp(reftr));
1522: PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));
1524: PetscCall(DMGetDimension(trdm, &trdim));
1525: PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
1526: tr->trNv[ct] = vEnd - vStart;
1527: PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
1528: PetscCall(VecGetLocalSize(coordinates, &Nc));
1529: 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);
1530: PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
1531: PetscCall(VecGetArrayRead(coordinates, &coords));
1532: PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
1533: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1535: PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
1536: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1537: for (n = 0; n < Nct; ++n) {
1538: /* Since points are 0-dimensional, coordinates make no sense */
1539: if (rct[n] == DM_POLYTOPE_POINT) continue;
1540: PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
1541: for (r = 0; r < rsize[n]; ++r) {
1542: PetscInt *closure = NULL;
1543: PetscInt clSize, cl, Nv = 0;
1545: PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
1546: PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
1547: PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1548: for (cl = 0; cl < clSize * 2; cl += 2) {
1549: const PetscInt sv = closure[cl];
1551: if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1552: }
1553: PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
1554: 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]);
1555: }
1556: }
1557: if (debug) {
1558: DMPolytopeType *rct;
1559: PetscInt *rsize, *rcone, *rornt;
1560: PetscInt v, dE = trdim, d, off = 0;
1562: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
1563: for (v = 0; v < tr->trNv[ct]; ++v) {
1564: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1565: for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
1566: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1567: }
1569: PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1570: for (n = 0; n < Nct; ++n) {
1571: if (rct[n] == DM_POLYTOPE_POINT) continue;
1572: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
1573: for (r = 0; r < rsize[n]; ++r) {
1574: PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
1575: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
1576: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1577: }
1578: }
1579: }
1580: PetscCall(DMDestroy(&refdm));
1581: PetscCall(DMDestroy(&trdm));
1582: PetscCall(DMPlexTransformDestroy(&reftr));
1583: }
1584: PetscFunctionReturn(PETSC_SUCCESS);
1585: }
1587: /*@C
1588: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
1590: Input Parameters:
1591: + tr - The `DMPlexTransform` object
1592: - ct - The cell type
1594: Output Parameters:
1595: + Nv - The number of transformed vertices in the closure of the reference cell of given type
1596: - trVerts - The coordinates of these vertices in the reference cell
1598: Level: developer
1600: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
1601: @*/
1602: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1603: {
1604: PetscFunctionBegin;
1605: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1606: if (Nv) *Nv = tr->trNv[ct];
1607: if (trVerts) *trVerts = tr->trVerts[ct];
1608: PetscFunctionReturn(PETSC_SUCCESS);
1609: }
1611: /*@C
1612: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
1614: Input Parameters:
1615: + tr - The `DMPlexTransform` object
1616: . ct - The cell type
1617: . rct - The subcell type
1618: - r - The subcell index
1620: Output Parameter:
1621: . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()`
1623: Level: developer
1625: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
1626: @*/
1627: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1628: {
1629: PetscFunctionBegin;
1630: if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
1631: PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
1632: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1633: PetscFunctionReturn(PETSC_SUCCESS);
1634: }
1636: /* Computes new vertex as the barycenter, or centroid */
1637: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1638: {
1639: PetscInt v, d;
1641: PetscFunctionBeginHot;
1642: PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
1643: for (d = 0; d < dE; ++d) out[d] = 0.0;
1644: for (v = 0; v < Nv; ++v)
1645: for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1646: for (d = 0; d < dE; ++d) out[d] /= Nv;
1647: PetscFunctionReturn(PETSC_SUCCESS);
1648: }
1650: /*@
1651: DMPlexTransformMapCoordinates - Calculate new coordinates for produced points
1653: Not collective
1655: Input Parameters:
1656: + tr - The `DMPlexTransform`
1657: . pct - The cell type of the parent, from whom the new cell is being produced
1658: . ct - The type being produced
1659: . p - The original point
1660: . r - The replica number requested for the produced cell type
1661: . Nv - Number of vertices in the closure of the parent cell
1662: . dE - Spatial dimension
1663: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
1665: Output Parameter:
1666: . out - The coordinates of the new vertices
1668: Level: intermediate
1670: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
1671: @*/
1672: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1673: {
1674: PetscFunctionBeginHot;
1675: if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1676: PetscFunctionReturn(PETSC_SUCCESS);
1677: }
1679: /*
1680: DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label
1682: Not Collective
1684: Input Parameters:
1685: + tr - The `DMPlexTransform`
1686: . label - The label in the transformed mesh
1687: . pp - The parent point in the original mesh
1688: . pct - The cell type of the parent point
1689: . p - The point in the transformed mesh
1690: . ct - The cell type of the point
1691: . r - The replica number of the point
1692: - val - The label value of the parent point
1694: Level: developer
1696: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
1697: */
1698: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
1699: {
1700: PetscFunctionBeginHot;
1701: if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
1702: PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
1703: PetscFunctionReturn(PETSC_SUCCESS);
1704: }
1706: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1707: {
1708: DM dm;
1709: IS valueIS;
1710: const PetscInt *values;
1711: PetscInt defVal, Nv, val;
1713: PetscFunctionBegin;
1714: PetscCall(DMPlexTransformGetDM(tr, &dm));
1715: PetscCall(DMLabelGetDefaultValue(label, &defVal));
1716: PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
1717: PetscCall(DMLabelGetValueIS(label, &valueIS));
1718: PetscCall(ISGetLocalSize(valueIS, &Nv));
1719: PetscCall(ISGetIndices(valueIS, &values));
1720: for (val = 0; val < Nv; ++val) {
1721: IS pointIS;
1722: const PetscInt *points;
1723: PetscInt numPoints, p;
1725: /* Ensure refined label is created with same number of strata as
1726: * original (even if no entries here). */
1727: PetscCall(DMLabelAddStratum(labelNew, values[val]));
1728: PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
1729: PetscCall(ISGetLocalSize(pointIS, &numPoints));
1730: PetscCall(ISGetIndices(pointIS, &points));
1731: for (p = 0; p < numPoints; ++p) {
1732: const PetscInt point = points[p];
1733: DMPolytopeType ct;
1734: DMPolytopeType *rct;
1735: PetscInt *rsize, *rcone, *rornt;
1736: PetscInt Nct, n, r, pNew = 0;
1738: PetscCall(DMPlexGetCellType(dm, point, &ct));
1739: PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1740: for (n = 0; n < Nct; ++n) {
1741: for (r = 0; r < rsize[n]; ++r) {
1742: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
1743: PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
1744: }
1745: }
1746: }
1747: PetscCall(ISRestoreIndices(pointIS, &points));
1748: PetscCall(ISDestroy(&pointIS));
1749: }
1750: PetscCall(ISRestoreIndices(valueIS, &values));
1751: PetscCall(ISDestroy(&valueIS));
1752: PetscFunctionReturn(PETSC_SUCCESS);
1753: }
1755: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1756: {
1757: DM dm;
1758: PetscInt numLabels, l;
1760: PetscFunctionBegin;
1761: PetscCall(DMPlexTransformGetDM(tr, &dm));
1762: PetscCall(DMGetNumLabels(dm, &numLabels));
1763: for (l = 0; l < numLabels; ++l) {
1764: DMLabel label, labelNew;
1765: const char *lname;
1766: PetscBool isDepth, isCellType;
1768: PetscCall(DMGetLabelName(dm, l, &lname));
1769: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
1770: if (isDepth) continue;
1771: PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
1772: if (isCellType) continue;
1773: PetscCall(DMCreateLabel(rdm, lname));
1774: PetscCall(DMGetLabel(dm, lname, &label));
1775: PetscCall(DMGetLabel(rdm, lname, &labelNew));
1776: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1777: }
1778: PetscFunctionReturn(PETSC_SUCCESS);
1779: }
1781: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1782: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1783: {
1784: DM dm;
1785: PetscInt Nf, f, Nds, s;
1787: PetscFunctionBegin;
1788: PetscCall(DMPlexTransformGetDM(tr, &dm));
1789: PetscCall(DMGetNumFields(dm, &Nf));
1790: for (f = 0; f < Nf; ++f) {
1791: DMLabel label, labelNew;
1792: PetscObject obj;
1793: const char *lname;
1795: PetscCall(DMGetField(rdm, f, &label, &obj));
1796: if (!label) continue;
1797: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1798: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1799: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1800: PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
1801: PetscCall(DMLabelDestroy(&labelNew));
1802: }
1803: PetscCall(DMGetNumDS(dm, &Nds));
1804: for (s = 0; s < Nds; ++s) {
1805: DMLabel label, labelNew;
1806: const char *lname;
1808: PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
1809: if (!label) continue;
1810: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1811: PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
1812: PetscCall(RefineLabel_Internal(tr, label, labelNew));
1813: PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
1814: PetscCall(DMLabelDestroy(&labelNew));
1815: }
1816: PetscFunctionReturn(PETSC_SUCCESS);
1817: }
1819: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1820: {
1821: DM dm;
1822: PetscSF sf, sfNew;
1823: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
1824: const PetscInt *localPoints;
1825: const PetscSFNode *remotePoints;
1826: PetscInt *localPointsNew;
1827: PetscSFNode *remotePointsNew;
1828: PetscInt pStartNew, pEndNew, pNew;
1829: /* Brute force algorithm */
1830: PetscSF rsf;
1831: PetscSection s;
1832: const PetscInt *rootdegree;
1833: PetscInt *rootPointsNew, *remoteOffsets;
1834: PetscInt numPointsNew, pStart, pEnd, p;
1836: PetscFunctionBegin;
1837: PetscCall(DMPlexTransformGetDM(tr, &dm));
1838: PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
1839: PetscCall(DMGetPointSF(dm, &sf));
1840: PetscCall(DMGetPointSF(rdm, &sfNew));
1841: /* Calculate size of new SF */
1842: PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
1843: if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
1844: for (l = 0; l < numLeaves; ++l) {
1845: const PetscInt p = localPoints[l];
1846: DMPolytopeType ct;
1847: DMPolytopeType *rct;
1848: PetscInt *rsize, *rcone, *rornt;
1849: PetscInt Nct, n;
1851: PetscCall(DMPlexGetCellType(dm, p, &ct));
1852: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1853: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
1854: }
1855: /* Send new root point numbers
1856: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1857: */
1858: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1859: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
1860: PetscCall(PetscSectionSetChart(s, pStart, pEnd));
1861: for (p = pStart; p < pEnd; ++p) {
1862: DMPolytopeType ct;
1863: DMPolytopeType *rct;
1864: PetscInt *rsize, *rcone, *rornt;
1865: PetscInt Nct, n;
1867: PetscCall(DMPlexGetCellType(dm, p, &ct));
1868: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1869: for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
1870: }
1871: PetscCall(PetscSectionSetUp(s));
1872: PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
1873: PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
1874: PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
1875: PetscCall(PetscFree(remoteOffsets));
1876: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1877: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1878: PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
1879: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1880: for (p = pStart; p < pEnd; ++p) {
1881: DMPolytopeType ct;
1882: DMPolytopeType *rct;
1883: PetscInt *rsize, *rcone, *rornt;
1884: PetscInt Nct, n, r, off;
1886: if (!rootdegree[p - pStart]) continue;
1887: PetscCall(PetscSectionGetOffset(s, p, &off));
1888: PetscCall(DMPlexGetCellType(dm, p, &ct));
1889: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1890: for (n = 0, m = 0; n < Nct; ++n) {
1891: for (r = 0; r < rsize[n]; ++r, ++m) {
1892: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1893: rootPointsNew[off + m] = pNew;
1894: }
1895: }
1896: }
1897: PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1898: PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
1899: PetscCall(PetscSFDestroy(&rsf));
1900: PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
1901: PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
1902: for (l = 0, m = 0; l < numLeaves; ++l) {
1903: const PetscInt p = localPoints[l];
1904: DMPolytopeType ct;
1905: DMPolytopeType *rct;
1906: PetscInt *rsize, *rcone, *rornt;
1907: PetscInt Nct, n, r, q, off;
1909: PetscCall(PetscSectionGetOffset(s, p, &off));
1910: PetscCall(DMPlexGetCellType(dm, p, &ct));
1911: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
1912: for (n = 0, q = 0; n < Nct; ++n) {
1913: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1914: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
1915: localPointsNew[m] = pNew;
1916: remotePointsNew[m].index = rootPointsNew[off + q];
1917: remotePointsNew[m].rank = remotePoints[l].rank;
1918: }
1919: }
1920: }
1921: PetscCall(PetscSectionDestroy(&s));
1922: PetscCall(PetscFree(rootPointsNew));
1923: /* SF needs sorted leaves to correctly calculate Gather */
1924: {
1925: PetscSFNode *rp, *rtmp;
1926: PetscInt *lp, *idx, *ltmp, i;
1928: PetscCall(PetscMalloc1(numLeavesNew, &idx));
1929: PetscCall(PetscMalloc1(numLeavesNew, &lp));
1930: PetscCall(PetscMalloc1(numLeavesNew, &rp));
1931: for (i = 0; i < numLeavesNew; ++i) {
1932: 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);
1933: idx[i] = i;
1934: }
1935: PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
1936: for (i = 0; i < numLeavesNew; ++i) {
1937: lp[i] = localPointsNew[idx[i]];
1938: rp[i] = remotePointsNew[idx[i]];
1939: }
1940: ltmp = localPointsNew;
1941: localPointsNew = lp;
1942: rtmp = remotePointsNew;
1943: remotePointsNew = rp;
1944: PetscCall(PetscFree(idx));
1945: PetscCall(PetscFree(ltmp));
1946: PetscCall(PetscFree(rtmp));
1947: }
1948: PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
1949: PetscFunctionReturn(PETSC_SUCCESS);
1950: }
1952: /*
1953: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`.
1955: Not Collective
1957: Input Parameters:
1958: + tr - The `DMPlexTransform`
1959: . ct - The type of the parent cell
1960: . rct - The type of the produced cell
1961: . r - The index of the produced cell
1962: - x - The localized coordinates for the parent cell
1964: Output Parameter:
1965: . xr - The localized coordinates for the produced cell
1967: Level: developer
1969: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
1970: */
1971: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1972: {
1973: PetscFE fe = NULL;
1974: PetscInt cdim, v, *subcellV;
1976: PetscFunctionBegin;
1977: PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
1978: PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
1979: PetscCall(PetscFEGetNumComponents(fe, &cdim));
1980: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
1981: PetscFunctionReturn(PETSC_SUCCESS);
1982: }
1984: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1985: {
1986: DM dm, cdm, cdmCell, cdmNew, cdmCellNew;
1987: PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
1988: Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
1989: const PetscScalar *coords;
1990: PetscScalar *coordsNew;
1991: const PetscReal *maxCell, *Lstart, *L;
1992: PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1993: PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;
1995: PetscFunctionBegin;
1996: PetscCall(DMPlexTransformGetDM(tr, &dm));
1997: PetscCall(DMGetCoordinateDM(dm, &cdm));
1998: PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
1999: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
2000: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
2001: if (localized) {
2002: /* Localize coordinates of new vertices */
2003: localizeVertices = PETSC_TRUE;
2004: /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
2005: if (!maxCell) localizeCells = PETSC_TRUE;
2006: }
2007: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2008: PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
2009: if (maxCell) {
2010: PetscReal maxCellNew[3];
2012: for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
2013: PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L));
2014: }
2015: PetscCall(DMGetCoordinateDim(rdm, &dE));
2016: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
2017: PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
2018: PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
2019: PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
2020: PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
2021: /* Localization should be inherited */
2022: /* Stefano calculates parent cells for each new cell for localization */
2023: /* Localized cells need coordinates of closure */
2024: for (v = vStartNew; v < vEndNew; ++v) {
2025: PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
2026: PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
2027: }
2028: PetscCall(PetscSectionSetUp(coordSectionNew));
2029: PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));
2031: if (localizeCells) {
2032: PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
2033: PetscCall(DMClone(cdmNew, &cdmCellNew));
2034: PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
2035: PetscCall(DMDestroy(&cdmCellNew));
2037: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
2038: PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
2039: PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
2040: PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
2041: PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));
2043: PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
2044: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2045: for (c = cStart; c < cEnd; ++c) {
2046: PetscInt dof;
2048: PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
2049: if (dof) {
2050: DMPolytopeType ct;
2051: DMPolytopeType *rct;
2052: PetscInt *rsize, *rcone, *rornt;
2053: PetscInt dim, cNew, Nct, n, r;
2055: PetscCall(DMPlexGetCellType(dm, c, &ct));
2056: dim = DMPolytopeTypeGetDim(ct);
2057: PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2058: /* This allows for different cell types */
2059: for (n = 0; n < Nct; ++n) {
2060: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2061: for (r = 0; r < rsize[n]; ++r) {
2062: PetscInt *closure = NULL;
2063: PetscInt clSize, cl, Nv = 0;
2065: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
2066: PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2067: for (cl = 0; cl < clSize * 2; cl += 2) {
2068: if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
2069: }
2070: PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
2071: PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
2072: PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
2073: }
2074: }
2075: }
2076: }
2077: PetscCall(PetscSectionSetUp(coordSectionCellNew));
2078: PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
2079: }
2080: PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
2081: {
2082: VecType vtype;
2083: PetscInt coordSizeNew, bs;
2084: const char *name;
2086: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
2087: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
2088: PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
2089: PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
2090: PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
2091: PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
2092: PetscCall(VecGetBlockSize(coordsLocal, &bs));
2093: PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
2094: PetscCall(VecGetType(coordsLocal, &vtype));
2095: PetscCall(VecSetType(coordsLocalNew, vtype));
2096: }
2097: PetscCall(VecGetArrayRead(coordsLocal, &coords));
2098: PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
2099: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2100: /* First set coordinates for vertices */
2101: for (p = pStart; p < pEnd; ++p) {
2102: DMPolytopeType ct;
2103: DMPolytopeType *rct;
2104: PetscInt *rsize, *rcone, *rornt;
2105: PetscInt Nct, n, r;
2106: PetscBool hasVertex = PETSC_FALSE;
2108: PetscCall(DMPlexGetCellType(dm, p, &ct));
2109: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2110: for (n = 0; n < Nct; ++n) {
2111: if (rct[n] == DM_POLYTOPE_POINT) {
2112: hasVertex = PETSC_TRUE;
2113: break;
2114: }
2115: }
2116: if (hasVertex) {
2117: const PetscScalar *icoords = NULL;
2118: const PetscScalar *array = NULL;
2119: PetscScalar *pcoords = NULL;
2120: PetscBool isDG;
2121: PetscInt Nc, Nv, v, d;
2123: PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2125: icoords = pcoords;
2126: Nv = Nc / dEo;
2127: if (ct != DM_POLYTOPE_POINT) {
2128: if (localizeVertices && maxCell) {
2129: PetscScalar anchor[3];
2131: for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
2132: for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
2133: }
2134: }
2135: for (n = 0; n < Nct; ++n) {
2136: if (rct[n] != DM_POLYTOPE_POINT) continue;
2137: for (r = 0; r < rsize[n]; ++r) {
2138: PetscScalar vcoords[3];
2139: PetscInt vNew, off;
2141: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
2142: PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
2143: PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
2144: PetscCall(DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
2145: }
2146: }
2147: PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
2148: }
2149: }
2150: PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
2151: PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
2152: PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
2153: PetscCall(VecDestroy(&coordsLocalNew));
2154: PetscCall(PetscSectionDestroy(&coordSectionNew));
2155: /* Then set coordinates for cells by localizing */
2156: if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
2157: else {
2158: VecType vtype;
2159: PetscInt coordSizeNew, bs;
2160: const char *name;
2162: PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
2163: PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
2164: PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
2165: PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
2166: PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
2167: PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
2168: PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
2169: PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
2170: PetscCall(VecGetType(coordsLocalCell, &vtype));
2171: PetscCall(VecSetType(coordsLocalCellNew, vtype));
2172: PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
2173: PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));
2175: for (p = pStart; p < pEnd; ++p) {
2176: DMPolytopeType ct;
2177: DMPolytopeType *rct;
2178: PetscInt *rsize, *rcone, *rornt;
2179: PetscInt dof = 0, Nct, n, r;
2181: PetscCall(DMPlexGetCellType(dm, p, &ct));
2182: PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
2183: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
2184: if (dof) {
2185: const PetscScalar *pcoords;
2187: PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords));
2188: for (n = 0; n < Nct; ++n) {
2189: const PetscInt Nr = rsize[n];
2191: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2192: for (r = 0; r < Nr; ++r) {
2193: PetscInt pNew, offNew;
2195: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2196: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2197: cell to the ones it produces. */
2198: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
2199: PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
2200: PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
2201: }
2202: }
2203: }
2204: }
2205: PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
2206: PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
2207: PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
2208: PetscCall(VecDestroy(&coordsLocalCellNew));
2209: PetscCall(PetscSectionDestroy(&coordSectionCellNew));
2210: }
2211: PetscFunctionReturn(PETSC_SUCCESS);
2212: }
2214: /*@
2215: DMPlexTransformApply - Execute the transformation, producing another `DM`
2217: Collective
2219: Input Parameters:
2220: + tr - The `DMPlexTransform` object
2221: - dm - The original `DM`
2223: Output Parameter:
2224: . tdm - The transformed `DM`
2226: Level: intermediate
2228: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()`
2229: @*/
2230: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
2231: {
2232: DM rdm;
2233: DMPlexInterpolatedFlag interp;
2234: PetscInt pStart, pEnd;
2236: PetscFunctionBegin;
2239: PetscAssertPointer(tdm, 3);
2240: PetscCall(PetscLogEventBegin(DMPLEX_Transform, tr, dm, 0, 0));
2241: PetscCall(DMPlexTransformSetDM(tr, dm));
2243: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
2244: PetscCall(DMSetType(rdm, DMPLEX));
2245: PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
2246: /* Calculate number of new points of each depth */
2247: PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
2248: PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
2249: /* Step 1: Set chart */
2250: PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
2251: PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
2252: /* Step 2: Set cone/support sizes (automatically stratifies) */
2253: PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
2254: /* Step 3: Setup refined DM */
2255: PetscCall(DMSetUp(rdm));
2256: /* Step 4: Set cones and supports (automatically symmetrizes) */
2257: PetscCall(DMPlexTransformSetCones(tr, rdm));
2258: /* Step 5: Create pointSF */
2259: PetscCall(DMPlexTransformCreateSF(tr, rdm));
2260: /* Step 6: Create labels */
2261: PetscCall(DMPlexTransformCreateLabels(tr, rdm));
2262: /* Step 7: Set coordinates */
2263: PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
2264: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm));
2265: // If the original DM was configured from options, the transformed DM should be as well
2266: rdm->setfromoptionscalled = dm->setfromoptionscalled;
2267: PetscCall(PetscLogEventEnd(DMPLEX_Transform, tr, dm, 0, 0));
2268: *tdm = rdm;
2269: PetscFunctionReturn(PETSC_SUCCESS);
2270: }
2272: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2273: {
2274: DMPlexTransform tr;
2275: DM cdm, rcdm;
2276: const char *prefix;
2278: PetscFunctionBegin;
2279: PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
2280: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
2281: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
2282: PetscCall(DMPlexTransformSetDM(tr, dm));
2283: PetscCall(DMPlexTransformSetFromOptions(tr));
2284: PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
2285: PetscCall(DMPlexTransformSetUp(tr));
2286: PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
2287: PetscCall(DMPlexTransformApply(tr, dm, rdm));
2288: PetscCall(DMCopyDisc(dm, *rdm));
2289: PetscCall(DMGetCoordinateDM(dm, &cdm));
2290: PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
2291: PetscCall(DMCopyDisc(cdm, rcdm));
2292: PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
2293: PetscCall(DMCopyDisc(dm, *rdm));
2294: PetscCall(DMPlexTransformDestroy(&tr));
2295: ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2296: PetscFunctionReturn(PETSC_SUCCESS);
2297: }