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: */
10: static PetscErrorCode DMLabelInitialize_Ephemeral(DMLabel);
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;
21: PetscFunctionBegin;
22: PetscCall(DMLabelEphemeralGetTransform(label, &tr));
23: PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
24: PetscCall(DMPlexTransformGetDM(tr, &dm));
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;
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: }
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;
57: PetscFunctionBegin;
58: PetscCall(DMLabelEphemeralGetTransform(label, &tr));
59: PetscCall(DMLabelEphemeralGetLabel(label, &olabel));
60: PetscCall(DMPlexTransformGetDM(tr, &dm));
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;
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: }
90: static PetscErrorCode DMLabelSetUp_Ephemeral(DMLabel label)
91: {
92: DMLabel olabel;
93: IS valueIS;
94: const PetscInt *values;
95: PetscInt defValue, Nv;
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: }
112: static PetscErrorCode DMLabelView_Ephemeral_Ascii(DMLabel label, PetscViewer viewer)
113: {
114: DMLabel olabel;
115: PetscMPIInt rank;
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;
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;
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: }
153: static PetscErrorCode DMLabelView_Ephemeral(DMLabel label, PetscViewer viewer)
154: {
155: PetscBool iascii;
157: PetscFunctionBegin;
158: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
159: if (iascii) PetscCall(DMLabelView_Ephemeral_Ascii(label, viewer));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: static PetscErrorCode DMLabelDuplicate_Ephemeral(DMLabel label, DMLabel *labelnew)
164: {
165: PetscObject obj;
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: }
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: }
186: /*MC
187: DMLABELEPHEMERAL = "ephemeral" - This type of `DMLabel` is used to label ephemeral meshes.
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.
194: Level: intermediate
196: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelCreate()`, `DMLabelSetType()`
197: M*/
199: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel label)
200: {
201: PetscFunctionBegin;
203: PetscCall(DMLabelInitialize_Ephemeral(label));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: /*@
208: DMLabelEphemeralGetLabel - Get the base label for this ephemeral label
210: Not Collective
212: Input Parameter:
213: . label - the `DMLabel`
215: Output Parameter:
216: . olabel - the base label for this ephemeral label
218: Level: intermediate
220: Note:
221: Ephemeral labels are produced automatically from a base label and a `DMPlexTransform`.
223: .seealso: `DMLabelEphemeralSetLabel()`, `DMLabelEphemeralGetTransform()`, `DMLabelSetType()`
224: @*/
225: PetscErrorCode DMLabelEphemeralGetLabel(DMLabel label, DMLabel *olabel)
226: {
227: PetscFunctionBegin;
228: PetscCall(PetscObjectQuery((PetscObject)label, "__original_label__", (PetscObject *)olabel));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*@
233: DMLabelEphemeralSetLabel - Set the base label for this ephemeral label
235: Not Collective
237: Input Parameters:
238: + label - the `DMLabel`
239: - olabel - the base label for this ephemeral label
241: Level: intermediate
243: Note:
244: Ephemeral labels are produced automatically from a base label and a `DMPlexTransform`.
246: .seealso: `DMLabelEphemeralGetLabel()`, `DMLabelEphemeralSetTransform()`, `DMLabelSetType()`
247: @*/
248: PetscErrorCode DMLabelEphemeralSetLabel(DMLabel label, DMLabel olabel)
249: {
250: PetscFunctionBegin;
251: PetscCall(PetscObjectCompose((PetscObject)label, "__original_label__", (PetscObject)olabel));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }