Actual source code: dmlabeleph.c
1: #include <petsc/private/dmlabelimpl.h>
2: #include <petscdmlabelephemeral.h>
4: /*
5: An emphemeral label is read-only
7: This initial implementation makes the assumption that the produced points have unique parents. If this is
8: not satisfied, hash tables can be introduced to ensure produced points are unique.
9: */
11: static PetscErrorCode DMLabelEphemeralComputeStratumSize_Private(DMLabel label, PetscInt value)
12: {
13: DMPlexTransform tr;
14: DM dm;
15: DMLabel olabel;
16: IS opointIS;
17: const PetscInt *opoints;
18: PetscInt Np = 0, Nop, op, v;
20: PetscFunctionBegin;
21: PetscCall(DMLabelEphemeralGetTransform(label, &tr));
22: PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
23: PetscCall(DMPlexTransformGetDM(tr, &dm));
25: PetscCall(DMLabelLookupStratum(olabel, value, &v));
26: PetscCall(DMLabelGetStratumIS(olabel, value, &opointIS));
27: PetscCall(ISGetLocalSize(opointIS, &Nop));
28: PetscCall(ISGetIndices(opointIS, &opoints));
29: for (op = 0; op < Nop; ++op) {
30: const PetscInt point = opoints[op];
31: DMPolytopeType ct;
32: DMPolytopeType *rct;
33: PetscInt *rsize, *rcone, *rornt;
34: PetscInt Nct;
36: PetscCall(DMPlexGetCellType(dm, point, &ct));
37: PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
38: for (PetscInt n = 0; n < Nct; ++n) Np += rsize[n];
39: }
40: PetscCall(ISRestoreIndices(opointIS, &opoints));
41: PetscCall(ISDestroy(&opointIS));
42: label->stratumSizes[v] = Np;
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: static PetscErrorCode DMLabelGetStratumIS_Ephemeral(DMLabel label, PetscInt v, IS *stratum)
47: {
48: DMPlexTransform tr;
49: DM dm;
50: DMLabel olabel;
51: IS opointIS;
52: const PetscInt *opoints;
53: PetscInt *points;
54: PetscInt Np, p, Nop, op;
56: PetscFunctionBegin;
57: PetscCall(DMLabelEphemeralGetTransform(label, &tr));
58: PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
59: PetscCall(DMPlexTransformGetDM(tr, &dm));
61: PetscCall(DMLabelGetStratumSize_Private(label, v, &Np));
62: PetscCall(PetscMalloc1(Np, &points));
63: PetscUseTypeMethod(olabel, getstratumis, v, &opointIS);
64: PetscCall(ISGetLocalSize(opointIS, &Nop));
65: PetscCall(ISGetIndices(opointIS, &opoints));
66: for (op = 0, p = 0; op < Nop; ++op) {
67: const PetscInt point = opoints[op];
68: DMPolytopeType ct;
69: DMPolytopeType *rct;
70: PetscInt *rsize, *rcone, *rornt;
71: PetscInt Nct, n, r, pNew = 0;
73: PetscCall(DMPlexGetCellType(dm, point, &ct));
74: PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
75: for (n = 0; n < Nct; ++n) {
76: for (r = 0; r < rsize[n]; ++r) {
77: PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
78: points[p++] = pNew;
79: }
80: }
81: }
82: PetscCheck(p == Np, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of stratum points %" PetscInt_FMT " != %" PetscInt_FMT " precomputed size", p, Np);
83: PetscCall(ISRestoreIndices(opointIS, &opoints));
84: PetscCall(ISDestroy(&opointIS));
85: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Np, points, PETSC_OWN_POINTER, stratum));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: static PetscErrorCode DMLabelSetUp_Ephemeral(DMLabel label)
90: {
91: DMLabel olabel;
92: IS valueIS;
93: const PetscInt *values;
94: PetscInt defValue, Nv;
96: PetscFunctionBegin;
97: PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
98: PetscCall(DMLabelGetDefaultValue(olabel, &defValue));
99: PetscCall(DMLabelSetDefaultValue(label, defValue));
100: PetscCall(DMLabelGetNumValues(olabel, &Nv));
101: PetscCall(DMLabelGetValueIS(olabel, &valueIS));
102: PetscCall(ISGetIndices(valueIS, &values));
103: PetscCall(DMLabelAddStrataIS(label, valueIS));
104: for (PetscInt v = 0; v < Nv; ++v) PetscCall(DMLabelEphemeralComputeStratumSize_Private(label, values[v]));
105: PetscCall(ISRestoreIndices(valueIS, &values));
106: PetscCall(ISDestroy(&valueIS));
107: label->readonly = PETSC_TRUE;
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: static PetscErrorCode DMLabelView_Ephemeral_Ascii(DMLabel label, PetscViewer viewer)
112: {
113: DMLabel olabel;
114: PetscMPIInt rank;
116: PetscFunctionBegin;
117: PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
118: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
119: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
120: if (olabel) {
121: IS valueIS;
122: const PetscInt *values;
123: PetscInt Nv, v;
124: const char *name;
126: PetscCall(PetscObjectGetName((PetscObject)label, &name));
127: PetscCall(PetscViewerASCIIPrintf(viewer, "Ephemeral Label '%s':\n", name));
128: PetscCall(DMLabelGetNumValues(olabel, &Nv));
129: PetscCall(DMLabelGetValueIS(olabel, &valueIS));
130: PetscCall(ISGetIndices(valueIS, &values));
131: for (v = 0; v < Nv; ++v) {
132: IS pointIS;
133: const PetscInt value = values[v];
134: const PetscInt *points;
135: PetscInt n, p;
137: PetscCall(DMLabelGetStratumIS(olabel, value, &pointIS));
138: PetscCall(ISGetIndices(pointIS, &points));
139: PetscCall(ISGetLocalSize(pointIS, &n));
140: for (p = 0; p < n; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, points[p], value));
141: PetscCall(ISRestoreIndices(pointIS, &points));
142: PetscCall(ISDestroy(&pointIS));
143: }
144: PetscCall(ISRestoreIndices(valueIS, &values));
145: PetscCall(ISDestroy(&valueIS));
146: }
147: PetscCall(PetscViewerFlush(viewer));
148: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode DMLabelView_Ephemeral(DMLabel label, PetscViewer viewer)
153: {
154: PetscBool iascii;
156: PetscFunctionBegin;
157: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
158: if (iascii) PetscCall(DMLabelView_Ephemeral_Ascii(label, viewer));
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: static PetscErrorCode DMLabelDuplicate_Ephemeral(DMLabel label, DMLabel *labelnew)
163: {
164: PetscObject olabel;
166: PetscFunctionBegin;
167: PetscCall(PetscObjectQuery((PetscObject)label, "__original_label__", &olabel));
168: PetscCall(PetscObjectCompose((PetscObject)*labelnew, "__original_label__", olabel));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: static PetscErrorCode DMLabelInitialize_Ephemeral(DMLabel label)
173: {
174: PetscFunctionBegin;
175: label->ops->view = DMLabelView_Ephemeral;
176: label->ops->setup = DMLabelSetUp_Ephemeral;
177: label->ops->duplicate = DMLabelDuplicate_Ephemeral;
178: label->ops->getstratumis = DMLabelGetStratumIS_Ephemeral;
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: /*MC
183: DMLABELEPHEMERAL = "ephemeral" - This type of `DMLabel` is used to label ephemeral meshes.
185: Ephemeral meshes are never concretely instantiated, but rather the answers to queries are created on the fly from a base mesh and a `DMPlexTransform` object.
186: For example, we could integrate over a refined mesh without ever storing the entire thing, just producing each cell closure one at a time. The label consists
187: of a base label and the same transform. Stratum are produced on demand, so that the time is in $O(max(M, N))$ where $M$ is the size of the original stratum,
188: and $N$ is the size of the output stratum. Queries take the same time, since we cannot sort the output.
190: Level: intermediate
192: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelCreate()`, `DMLabelSetType()`
193: M*/
195: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel label)
196: {
197: PetscFunctionBegin;
199: PetscCall(DMLabelInitialize_Ephemeral(label));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: DMLabelEphemeralGetLabel - Get the base label for this ephemeral label
206: Not Collective
208: Input Parameter:
209: . label - the `DMLabel`
211: Output Parameter:
212: . olabel - the base label for this ephemeral label
214: Level: intermediate
216: Note:
217: Ephemeral labels are produced automatically from a base label and a `DMPlexTransform`.
219: .seealso: `DMLabelEphemeralSetLabel()`, `DMLabelEphemeralGetTransform()`, `DMLabelSetType()`
220: @*/
221: PetscErrorCode DMLabelEphemeralGetLabel(DMLabel label, DMLabel *olabel)
222: {
223: PetscFunctionBegin;
224: PetscCall(PetscObjectQuery((PetscObject)label, "__original_label__", (PetscObject *)olabel));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /*@
229: DMLabelEphemeralSetLabel - Set the base label for this ephemeral label
231: Not Collective
233: Input Parameters:
234: + label - the `DMLabel`
235: - olabel - the base label for this ephemeral label
237: Level: intermediate
239: Note:
240: Ephemeral labels are produced automatically from a base label and a `DMPlexTransform`.
242: .seealso: `DMLabelEphemeralGetLabel()`, `DMLabelEphemeralSetTransform()`, `DMLabelSetType()`
243: @*/
244: PetscErrorCode DMLabelEphemeralSetLabel(DMLabel label, DMLabel olabel)
245: {
246: PetscFunctionBegin;
247: PetscCall(PetscObjectCompose((PetscObject)label, "__original_label__", (PetscObject)olabel));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }