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 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; 19 20 PetscFunctionBegin; 21 PetscCall(DMLabelEphemeralGetTransform(label, &tr)); 22 PetscCall(DMLabelEphemeralGetLabel(label, &olabel)); 23 PetscCall(DMPlexTransformGetDM(tr, &dm)); 24 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; 35 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 } 45 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; 55 56 PetscFunctionBegin; 57 PetscCall(DMLabelEphemeralGetTransform(label, &tr)); 58 PetscCall(DMLabelEphemeralGetLabel(label, &olabel)); 59 PetscCall(DMPlexTransformGetDM(tr, &dm)); 60 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; 72 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 } 88 89 static PetscErrorCode DMLabelSetUp_Ephemeral(DMLabel label) 90 { 91 DMLabel olabel; 92 IS valueIS; 93 const PetscInt *values; 94 PetscInt defValue, Nv; 95 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 } 110 111 static PetscErrorCode DMLabelView_Ephemeral_Ascii(DMLabel label, PetscViewer viewer) 112 { 113 DMLabel olabel; 114 PetscMPIInt rank; 115 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; 125 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; 136 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 } 151 152 static PetscErrorCode DMLabelView_Ephemeral(DMLabel label, PetscViewer viewer) 153 { 154 PetscBool iascii; 155 156 PetscFunctionBegin; 157 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 158 if (iascii) PetscCall(DMLabelView_Ephemeral_Ascii(label, viewer)); 159 PetscFunctionReturn(PETSC_SUCCESS); 160 } 161 162 static PetscErrorCode DMLabelDuplicate_Ephemeral(DMLabel label, DMLabel *labelnew) 163 { 164 PetscObject olabel; 165 166 PetscFunctionBegin; 167 PetscCall(PetscObjectQuery((PetscObject)label, "__original_label__", &olabel)); 168 PetscCall(PetscObjectCompose((PetscObject)*labelnew, "__original_label__", olabel)); 169 PetscFunctionReturn(PETSC_SUCCESS); 170 } 171 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 } 181 182 /*MC 183 DMLABELEPHEMERAL = "ephemeral" - This type of `DMLabel` is used to label ephemeral meshes. 184 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. 189 190 Level: intermediate 191 192 .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelCreate()`, `DMLabelSetType()` 193 M*/ 194 195 PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel label) 196 { 197 PetscFunctionBegin; 198 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 199 PetscCall(DMLabelInitialize_Ephemeral(label)); 200 PetscFunctionReturn(PETSC_SUCCESS); 201 } 202 203 /*@ 204 DMLabelEphemeralGetLabel - Get the base label for this ephemeral label 205 206 Not collective 207 208 Input Parameter: 209 . label - the DMLabel 210 211 Output Paramater: 212 . olabel - the base label for this ephemeral label 213 214 Note: 215 Ephemeral labels are produced automatically from a base label and a `DMPlexTransform`. 216 217 Level: intermediate 218 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 } 227 228 /*@ 229 DMLabelEphemeralSetLabel - Set the base label for this ephemeral label 230 231 Not collective 232 233 Input Parameters: 234 + label - the DMLabel 235 - olabel - the base label for this ephemeral label 236 237 Note: 238 Ephemeral labels are produced automatically from a base label and a `DMPlexTransform`. 239 240 Level: intermediate 241 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 } 250