1 #include <petsc/private/dmlabelimpl.h> /*I "petscdmlabelephemeral.h" I*/
2 #include <petscdmlabelephemeral.h>
3
4 /*
5 An emphemeral label is read-only
6
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 */
10 static PetscErrorCode DMLabelInitialize_Ephemeral(DMLabel);
11
DMLabelEphemeralComputeStratumSize_Private(DMLabel label,PetscInt value)12 static PetscErrorCode DMLabelEphemeralComputeStratumSize_Private(DMLabel label, PetscInt value)
13 {
14 DMPlexTransform tr;
15 DM dm;
16 DMLabel olabel;
17 IS opointIS;
18 const PetscInt *opoints;
19 PetscInt Np = 0, Nop, op, v;
20
21 PetscFunctionBegin;
22 PetscCall(DMLabelEphemeralGetTransform(label, &tr));
23 PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
24 PetscCall(DMPlexTransformGetDM(tr, &dm));
25
26 PetscCall(DMLabelLookupStratum(olabel, value, &v));
27 PetscCall(DMLabelGetStratumIS(olabel, value, &opointIS));
28 PetscCall(ISGetLocalSize(opointIS, &Nop));
29 PetscCall(ISGetIndices(opointIS, &opoints));
30 for (op = 0; op < Nop; ++op) {
31 const PetscInt point = opoints[op];
32 DMPolytopeType ct;
33 DMPolytopeType *rct;
34 PetscInt *rsize, *rcone, *rornt;
35 PetscInt Nct;
36
37 PetscCall(DMPlexGetCellType(dm, point, &ct));
38 PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
39 for (PetscInt n = 0; n < Nct; ++n) Np += rsize[n];
40 }
41 PetscCall(ISRestoreIndices(opointIS, &opoints));
42 PetscCall(ISDestroy(&opointIS));
43 label->stratumSizes[v] = Np;
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46
DMLabelGetStratumIS_Ephemeral(DMLabel label,PetscInt v,IS * stratum)47 static PetscErrorCode DMLabelGetStratumIS_Ephemeral(DMLabel label, PetscInt v, IS *stratum)
48 {
49 DMPlexTransform tr;
50 DM dm;
51 DMLabel olabel;
52 IS opointIS;
53 const PetscInt *opoints;
54 PetscInt *points;
55 PetscInt Np, p, Nop, op;
56
57 PetscFunctionBegin;
58 PetscCall(DMLabelEphemeralGetTransform(label, &tr));
59 PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
60 PetscCall(DMPlexTransformGetDM(tr, &dm));
61
62 PetscCall(DMLabelGetStratumSize_Private(label, v, &Np));
63 PetscCall(PetscMalloc1(Np, &points));
64 PetscUseTypeMethod(olabel, getstratumis, v, &opointIS);
65 PetscCall(ISGetLocalSize(opointIS, &Nop));
66 PetscCall(ISGetIndices(opointIS, &opoints));
67 for (op = 0, p = 0; op < Nop; ++op) {
68 const PetscInt point = opoints[op];
69 DMPolytopeType ct;
70 DMPolytopeType *rct;
71 PetscInt *rsize, *rcone, *rornt;
72 PetscInt Nct, n, r, pNew = 0;
73
74 PetscCall(DMPlexGetCellType(dm, point, &ct));
75 PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
76 for (n = 0; n < Nct; ++n) {
77 for (r = 0; r < rsize[n]; ++r) {
78 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
79 points[p++] = pNew;
80 }
81 }
82 }
83 PetscCheck(p == Np, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of stratum points %" PetscInt_FMT " != %" PetscInt_FMT " precomputed size", p, Np);
84 PetscCall(ISRestoreIndices(opointIS, &opoints));
85 PetscCall(ISDestroy(&opointIS));
86 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Np, points, PETSC_OWN_POINTER, stratum));
87 PetscFunctionReturn(PETSC_SUCCESS);
88 }
89
DMLabelSetUp_Ephemeral(DMLabel label)90 static PetscErrorCode DMLabelSetUp_Ephemeral(DMLabel label)
91 {
92 DMLabel olabel;
93 IS valueIS;
94 const PetscInt *values;
95 PetscInt defValue, Nv;
96
97 PetscFunctionBegin;
98 PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
99 PetscCall(DMLabelGetDefaultValue(olabel, &defValue));
100 PetscCall(DMLabelSetDefaultValue(label, defValue));
101 PetscCall(DMLabelGetNumValues(olabel, &Nv));
102 PetscCall(DMLabelGetValueIS(olabel, &valueIS));
103 PetscCall(ISGetIndices(valueIS, &values));
104 PetscCall(DMLabelAddStrataIS(label, valueIS));
105 for (PetscInt v = 0; v < Nv; ++v) PetscCall(DMLabelEphemeralComputeStratumSize_Private(label, values[v]));
106 PetscCall(ISRestoreIndices(valueIS, &values));
107 PetscCall(ISDestroy(&valueIS));
108 label->readonly = PETSC_TRUE;
109 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111
DMLabelView_Ephemeral_Ascii(DMLabel label,PetscViewer viewer)112 static PetscErrorCode DMLabelView_Ephemeral_Ascii(DMLabel label, PetscViewer viewer)
113 {
114 DMLabel olabel;
115 PetscMPIInt rank;
116
117 PetscFunctionBegin;
118 PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
119 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
120 PetscCall(PetscViewerASCIIPushSynchronized(viewer));
121 if (olabel) {
122 IS valueIS;
123 const PetscInt *values;
124 PetscInt Nv, v;
125 const char *name;
126
127 PetscCall(PetscObjectGetName((PetscObject)label, &name));
128 PetscCall(PetscViewerASCIIPrintf(viewer, "Ephemeral Label '%s':\n", name));
129 PetscCall(DMLabelGetNumValues(olabel, &Nv));
130 PetscCall(DMLabelGetValueIS(olabel, &valueIS));
131 PetscCall(ISGetIndices(valueIS, &values));
132 for (v = 0; v < Nv; ++v) {
133 IS pointIS;
134 const PetscInt value = values[v];
135 const PetscInt *points;
136 PetscInt n, p;
137
138 PetscCall(DMLabelGetStratumIS(olabel, value, &pointIS));
139 PetscCall(ISGetIndices(pointIS, &points));
140 PetscCall(ISGetLocalSize(pointIS, &n));
141 for (p = 0; p < n; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, points[p], value));
142 PetscCall(ISRestoreIndices(pointIS, &points));
143 PetscCall(ISDestroy(&pointIS));
144 }
145 PetscCall(ISRestoreIndices(valueIS, &values));
146 PetscCall(ISDestroy(&valueIS));
147 }
148 PetscCall(PetscViewerFlush(viewer));
149 PetscCall(PetscViewerASCIIPopSynchronized(viewer));
150 PetscFunctionReturn(PETSC_SUCCESS);
151 }
152
DMLabelView_Ephemeral(DMLabel label,PetscViewer viewer)153 static PetscErrorCode DMLabelView_Ephemeral(DMLabel label, PetscViewer viewer)
154 {
155 PetscBool isascii;
156
157 PetscFunctionBegin;
158 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
159 if (isascii) PetscCall(DMLabelView_Ephemeral_Ascii(label, viewer));
160 PetscFunctionReturn(PETSC_SUCCESS);
161 }
162
DMLabelDuplicate_Ephemeral(DMLabel label,DMLabel * labelnew)163 static PetscErrorCode DMLabelDuplicate_Ephemeral(DMLabel label, DMLabel *labelnew)
164 {
165 PetscObject obj;
166
167 PetscFunctionBegin;
168 PetscCall(PetscObjectQuery((PetscObject)label, "__original_label__", &obj));
169 PetscCall(PetscObjectCompose((PetscObject)*labelnew, "__original_label__", obj));
170 PetscCall(PetscObjectQuery((PetscObject)label, "__transform__", &obj));
171 PetscCall(PetscObjectCompose((PetscObject)*labelnew, "__transform__", obj));
172 PetscCall(DMLabelInitialize_Ephemeral(*labelnew));
173 PetscFunctionReturn(PETSC_SUCCESS);
174 }
175
DMLabelInitialize_Ephemeral(DMLabel label)176 static PetscErrorCode DMLabelInitialize_Ephemeral(DMLabel label)
177 {
178 PetscFunctionBegin;
179 label->ops->view = DMLabelView_Ephemeral;
180 label->ops->setup = DMLabelSetUp_Ephemeral;
181 label->ops->duplicate = DMLabelDuplicate_Ephemeral;
182 label->ops->getstratumis = DMLabelGetStratumIS_Ephemeral;
183 PetscFunctionReturn(PETSC_SUCCESS);
184 }
185
186 /*MC
187 DMLABELEPHEMERAL = "ephemeral" - This type of `DMLabel` is used to label ephemeral meshes.
188
189 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.
190 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
191 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,
192 and $N$ is the size of the output stratum. Queries take the same time, since we cannot sort the output.
193
194 Level: intermediate
195
196 .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelCreate()`, `DMLabelSetType()`
197 M*/
198
DMLabelCreate_Ephemeral(DMLabel label)199 PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel label)
200 {
201 PetscFunctionBegin;
202 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1);
203 PetscCall(DMLabelInitialize_Ephemeral(label));
204 PetscFunctionReturn(PETSC_SUCCESS);
205 }
206
207 /*@
208 DMLabelEphemeralGetLabel - Get the base label for this ephemeral label
209
210 Not Collective
211
212 Input Parameter:
213 . label - the `DMLabel`
214
215 Output Parameter:
216 . olabel - the base label for this ephemeral label
217
218 Level: intermediate
219
220 Note:
221 Ephemeral labels are produced automatically from a base label and a `DMPlexTransform`.
222
223 .seealso: `DMLabelEphemeralSetLabel()`, `DMLabelEphemeralGetTransform()`, `DMLabelSetType()`
224 @*/
DMLabelEphemeralGetLabel(DMLabel label,DMLabel * olabel)225 PetscErrorCode DMLabelEphemeralGetLabel(DMLabel label, DMLabel *olabel)
226 {
227 PetscFunctionBegin;
228 PetscCall(PetscObjectQuery((PetscObject)label, "__original_label__", (PetscObject *)olabel));
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
232 /*@
233 DMLabelEphemeralSetLabel - Set the base label for this ephemeral label
234
235 Not Collective
236
237 Input Parameters:
238 + label - the `DMLabel`
239 - olabel - the base label for this ephemeral label
240
241 Level: intermediate
242
243 Note:
244 Ephemeral labels are produced automatically from a base label and a `DMPlexTransform`.
245
246 .seealso: `DMLabelEphemeralGetLabel()`, `DMLabelEphemeralSetTransform()`, `DMLabelSetType()`
247 @*/
DMLabelEphemeralSetLabel(DMLabel label,DMLabel olabel)248 PetscErrorCode DMLabelEphemeralSetLabel(DMLabel label, DMLabel olabel)
249 {
250 PetscFunctionBegin;
251 PetscCall(PetscObjectCompose((PetscObject)label, "__original_label__", (PetscObject)olabel));
252 PetscFunctionReturn(PETSC_SUCCESS);
253 }
254