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: }