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 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 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 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 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 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 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 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 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 @*/ 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 @*/ 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