1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #include <petsc/private/petscfeimpl.h> 4 5 static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscDS prob, PetscReal time, PetscFECellGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[], 6 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, 7 PetscScalar values[]) 8 { 9 PetscInt Nf, *Nc, f, totDim, spDim, d, v; 10 PetscErrorCode ierr; 11 12 PetscFunctionBeginHot; 13 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 14 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 15 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 16 /* Get values for closure */ 17 for (f = 0, v = 0; f < Nf; ++f) { 18 void * const ctx = ctxs ? ctxs[f] : NULL; 19 20 if (!sp[f]) continue; 21 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 22 for (d = 0; d < spDim; ++d, ++v) { 23 if (funcs[f]) { 24 if (isFE[f]) {ierr = PetscDualSpaceApply(sp[f], d, time, fegeom, Nc[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 25 else {ierr = PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 26 } else { 27 values[v] = 0.0; 28 } 29 } 30 } 31 PetscFunctionReturn(0); 32 } 33 34 static PetscErrorCode DMProjectPoint_Field_Private(DM dm, PetscDS prob, DM dmAux, PetscDS probAux, PetscReal time, Vec localU, Vec localA, PetscFECellGeom *fegeom, PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[], 35 PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 36 void (**funcs)(PetscInt, PetscInt, PetscInt, 37 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 38 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 39 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, 40 PetscScalar values[]) 41 { 42 PetscSection section, sectionAux = NULL; 43 PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *refSpaceDer, *refSpaceDerAux = NULL, *bc; 44 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 45 PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL; 46 const PetscScalar *constants; 47 PetscReal *x; 48 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 49 const PetscInt dim = fegeom->dim, dimEmbed = fegeom->dimEmbed; 50 PetscInt dimAux = 0, numConstants, Nf, NfAux = 0, f, spDim, d, c, v, tp = 0; 51 PetscErrorCode ierr; 52 53 PetscFunctionBeginHot; 54 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 55 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 56 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 57 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 58 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 59 ierr = PetscDSGetEvaluationArrays(prob, &u, &bc /*&u_t*/, &u_x);CHKERRQ(ierr); 60 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 61 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 62 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 63 ierr = DMPlexVecGetClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 64 if (dmAux) { 65 DMLabel spmap; 66 PetscInt subp; 67 68 ierr = DMPlexGetSubpointMap(dmAux, &spmap);CHKERRQ(ierr); 69 if (spmap) { 70 IS subpointIS; 71 const PetscInt *subpoints; 72 PetscInt numSubpoints; 73 74 ierr = DMPlexCreateSubpointIS(dmAux, &subpointIS);CHKERRQ(ierr); 75 ierr = ISGetLocalSize(subpointIS, &numSubpoints);CHKERRQ(ierr); 76 ierr = ISGetIndices(subpointIS, &subpoints);CHKERRQ(ierr); 77 ierr = PetscFindInt(p, numSubpoints, subpoints, &subp);CHKERRQ(ierr); 78 if (subp < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", p); 79 ierr = ISRestoreIndices(subpointIS, &subpoints);CHKERRQ(ierr); 80 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 81 } else subp = p; 82 ierr = PetscDSGetSpatialDimension(probAux, &dimAux);CHKERRQ(ierr); 83 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 84 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 85 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 86 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 87 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 88 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 89 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);CHKERRQ(ierr); 90 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 91 ierr = DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);CHKERRQ(ierr); 92 } 93 /* Get values for closure */ 94 for (f = 0, v = 0; f < Nf; ++f) { 95 if (!sp[f]) continue; 96 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 97 for (d = 0; d < spDim; ++d, ++v) { 98 if (funcs[f]) { 99 PetscQuadrature quad; 100 const PetscReal *points, *weights; 101 PetscInt numPoints, q; 102 103 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 104 ierr = PetscQuadratureGetData(quad, NULL, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 105 for (q = 0; q < numPoints; ++q, ++tp) { 106 CoordinatesRefToReal(dimEmbed, dim, fegeom->v0, fegeom->J, &points[q*dim], x); 107 EvaluateFieldJets(dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, refSpaceDer, fegeom->invJ, coefficients, coefficients_t, u, u_x, u_t); 108 if (probAux) {EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, tp, basisTabAux, basisDerTabAux, refSpaceDerAux, fegeom->invJ, coefficientsAux, coefficientsAux_t, a, a_x, a_t);} 109 (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, x, numConstants, constants, bc); 110 if (Ncc) { 111 for (c = 0; c < Ncc; ++c) values[v] += bc[comps[c]]*weights[comps[c]]; 112 } else { 113 for (c = 0; c < Nc[f]; ++c) values[v] += bc[c]*weights[c]; 114 } 115 } 116 } else { 117 values[v] = 0.0; 118 } 119 } 120 } 121 ierr = DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 122 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);CHKERRQ(ierr);} 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode DMProjectPoint_Private(DM dm, PetscDS prob, DM dmAux, PetscDS probAux, PetscInt effectiveHeight, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[], 127 PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[], 128 PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 129 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[]) 130 { 131 PetscFECellGeom fegeom; 132 PetscFVCellGeom fvgeom; 133 PetscInt dim, dimEmbed; 134 PetscErrorCode ierr; 135 136 PetscFunctionBeginHot; 137 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 138 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 139 if (hasFE) { 140 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, fegeom.invJ, &fegeom.detJ);CHKERRQ(ierr); 141 fegeom.dim = dim - effectiveHeight; 142 fegeom.dimEmbed = dimEmbed; 143 } 144 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 145 switch (type) { 146 case DM_BC_ESSENTIAL: 147 case DM_BC_NATURAL: 148 ierr = DMProjectPoint_Func_Private(dm, prob, time, &fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);CHKERRQ(ierr);break; 149 case DM_BC_ESSENTIAL_FIELD: 150 case DM_BC_NATURAL_FIELD: 151 ierr = DMProjectPoint_Field_Private(dm, prob, dmAux, probAux, time, localU, localA, &fegeom, sp, p, Ncc, comps, 152 basisTab, basisDerTab, basisTabAux, basisDerTabAux, 153 (void (**)(PetscInt, PetscInt, PetscInt, 154 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 155 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 156 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 157 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type); 158 } 159 PetscFunctionReturn(0); 160 } 161 162 /* 163 This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM 164 165 There are several different scenarios: 166 167 1) Volumetric mesh with volumetric auxiliary data 168 169 Here minHeight=0 since we loop over cells. 170 171 2) Boundary mesh with boundary auxiliary data 172 173 Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation. 174 175 3) Volumetric mesh with boundary auxiliary data 176 177 Here minHeight=1 and auxbd=PETSC_TRUE since we loop over faces and use data only supported on those faces. This is common when imposing Dirichlet boundary conditions. 178 179 The maxHeight is used to support enforcement of constraints in DMForest. 180 181 If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it. 182 183 If we are using an input field (DM_BC_ESSENTIAL_FIELD or DM_BC_NATURAL_FIELD), we need to evaluate it at all the quadrature points of the dual basis functionals. 184 - We use effectiveHeight to mean the height above our incoming DS. For example, if the DS is for a submesh then the effective height is zero, whereas if the DS 185 is for the volumetric mesh, but we are iterating over a surface, then the effective height is nonzero. When th effective height is nonzero, we need to extract 186 dual spaces for the boundary from our input spaces. 187 - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them. 188 189 We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration 190 191 If we have a label, we iterate over those points. This will probably break the maxHeight functionality since we do not check the height of those points. 192 */ 193 static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, 194 PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[], 195 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, 196 InsertMode mode, Vec localX) 197 { 198 DM dmAux = NULL; 199 PetscDS prob, probAux = NULL; 200 Vec localA = NULL; 201 PetscSection section; 202 PetscDualSpace *sp, *cellsp; 203 PetscReal **basisTab = NULL, **basisDerTab = NULL, **basisTabAux = NULL, **basisDerTabAux = NULL; 204 PetscInt *Nc; 205 PetscInt dim, dimEmbed, minHeight, maxHeight, h, Nf, NfAux = 0, f; 206 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE; 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 211 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);CHKERRQ(ierr); 212 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 213 ierr = DMPlexGetVTKCellHeight(dm, &minHeight);CHKERRQ(ierr); 214 /* Auxiliary information can only be used with interpolation of field functions */ 215 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 216 if (!minHeight && dmAux) { 217 DMLabel spmap; 218 219 /* If dmAux is a surface, then force the projection to take place over a surface */ 220 ierr = DMPlexGetSubpointMap(dmAux, &spmap);CHKERRQ(ierr); 221 if (spmap) {minHeight = 1; auxBd = PETSC_TRUE;} 222 } 223 } 224 ierr = DMPlexGetMaxProjectionHeight(dm, &maxHeight);CHKERRQ(ierr); 225 maxHeight = PetscMax(maxHeight, minHeight); 226 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %d not in [0, %d)\n", maxHeight, dim);} 227 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 228 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 229 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 230 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 231 if (dmAux) { 232 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 233 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 234 } 235 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 236 ierr = PetscMalloc2(Nf, &isFE, Nf, &sp);CHKERRQ(ierr); 237 if (maxHeight > 0) {ierr = PetscMalloc1(Nf, &cellsp);CHKERRQ(ierr);} 238 else {cellsp = sp;} 239 if (localU && localU != localX) {ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);CHKERRQ(ierr);} 240 /* Get cell dual spaces */ 241 for (f = 0; f < Nf; ++f) { 242 PetscObject obj; 243 PetscClassId id; 244 245 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 246 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 247 if (id == PETSCFE_CLASSID) { 248 PetscFE fe = (PetscFE) obj; 249 250 hasFE = PETSC_TRUE; 251 isFE[f] = PETSC_TRUE; 252 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 253 } else if (id == PETSCFV_CLASSID) { 254 PetscFV fv = (PetscFV) obj; 255 256 hasFV = PETSC_TRUE; 257 isFE[f] = PETSC_FALSE; 258 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 259 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 260 } 261 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 262 PetscInt effectiveHeight = auxBd ? minHeight : 0; 263 PetscFE fem, subfem; 264 PetscReal *points; 265 PetscInt numPoints, spDim, qdim = 0, d; 266 267 if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation"); 268 numPoints = 0; 269 for (f = 0; f < Nf; ++f) { 270 if (!effectiveHeight) {sp[f] = cellsp[f];} 271 else {ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr);} 272 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 273 for (d = 0; d < spDim; ++d) { 274 if (funcs[f]) { 275 PetscQuadrature quad; 276 PetscInt Nq; 277 278 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 279 ierr = PetscQuadratureGetData(quad, &qdim, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 280 numPoints += Nq; 281 } 282 } 283 } 284 ierr = PetscMalloc1(numPoints*qdim, &points);CHKERRQ(ierr); 285 numPoints = 0; 286 for (f = 0; f < Nf; ++f) { 287 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 288 for (d = 0; d < spDim; ++d) { 289 if (funcs[f]) { 290 PetscQuadrature quad; 291 const PetscReal *qpoints; 292 PetscInt Nq, q; 293 294 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 295 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL);CHKERRQ(ierr); 296 for (q = 0; q < Nq*qdim; ++q) points[numPoints*qdim+q] = qpoints[q]; 297 numPoints += Nq; 298 } 299 } 300 } 301 ierr = PetscMalloc4(Nf, &basisTab, Nf, &basisDerTab, NfAux, &basisTabAux, NfAux, &basisDerTabAux);CHKERRQ(ierr); 302 for (f = 0; f < Nf; ++f) { 303 if (!isFE[f]) continue; 304 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 305 if (!effectiveHeight) {subfem = fem;} 306 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 307 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 308 } 309 for (f = 0; f < NfAux; ++f) { 310 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 311 if (!effectiveHeight || auxBd) {subfem = fem;} 312 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 313 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 314 } 315 ierr = PetscFree(points);CHKERRQ(ierr); 316 } 317 /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */ 318 for (h = minHeight; h <= maxHeight; h++) { 319 PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight); 320 PetscDS probEff = prob; 321 PetscScalar *values; 322 PetscBool *fieldActive; 323 PetscInt pStart, pEnd, p, spDim, totDim, numValues; 324 325 if (effectiveHeight) {ierr = PetscDSGetHeightSubspace(prob, effectiveHeight, &probEff);CHKERRQ(ierr);} 326 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 327 if (!h) { 328 PetscInt cEndInterior; 329 330 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 331 pEnd = cEndInterior < 0 ? pEnd : cEndInterior; 332 } 333 if (pEnd <= pStart) continue; 334 /* Compute totDim, the number of dofs in the closure of a point at this height */ 335 totDim = 0; 336 for (f = 0; f < Nf; ++f) { 337 if (!effectiveHeight) { 338 sp[f] = cellsp[f]; 339 } else { 340 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr); 341 if (!sp[f]) continue; 342 } 343 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 344 totDim += spDim; 345 } 346 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 347 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim); 348 if (!totDim) continue; 349 /* Loop over points at this height */ 350 ierr = DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 351 ierr = DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 352 for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 353 if (label) { 354 PetscInt i; 355 356 for (i = 0; i < numIds; ++i) { 357 IS pointIS; 358 const PetscInt *points; 359 PetscInt n; 360 361 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 362 if (!pointIS) continue; /* No points with that id on this process */ 363 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 364 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 365 for (p = 0; p < n; ++p) { 366 const PetscInt point = points[p]; 367 368 if ((point < pStart) || (point >= pEnd)) continue; 369 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 370 ierr = DMProjectPoint_Private(dm, probEff, dmAux, probAux, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, point, Ncc, comps, basisTab, basisDerTab, basisTabAux, basisDerTabAux, type, funcs, ctxs, fieldActive, values); 371 if (ierr) { 372 PetscErrorCode ierr2; 373 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 374 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 375 CHKERRQ(ierr); 376 } 377 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, Ncc, comps, values, mode);CHKERRQ(ierr); 378 } 379 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 380 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 381 } 382 } else { 383 for (p = pStart; p < pEnd; ++p) { 384 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 385 ierr = DMProjectPoint_Private(dm, probEff, dmAux, probAux, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, p, Ncc, comps, basisTab, basisDerTab, basisTabAux, basisDerTabAux, type, funcs, ctxs, fieldActive, values); 386 if (ierr) { 387 PetscErrorCode ierr2; 388 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 389 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 390 CHKERRQ(ierr); 391 } 392 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, Ncc, comps, values, mode);CHKERRQ(ierr); 393 } 394 } 395 ierr = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 396 ierr = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 397 } 398 /* Cleanup */ 399 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 400 PetscInt effectiveHeight = auxBd ? minHeight : 0; 401 PetscFE fem, subfem; 402 403 for (f = 0; f < Nf; ++f) { 404 if (!isFE[f]) continue; 405 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 406 if (!effectiveHeight) {subfem = fem;} 407 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 408 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 409 } 410 for (f = 0; f < NfAux; ++f) { 411 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 412 if (!effectiveHeight || auxBd) {subfem = fem;} 413 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 414 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 415 } 416 ierr = PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);CHKERRQ(ierr); 417 } 418 ierr = PetscFree2(isFE, sp);CHKERRQ(ierr); 419 if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 420 PetscFunctionReturn(0); 421 } 422 423 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 424 { 425 PetscErrorCode ierr; 426 427 PetscFunctionBegin; 428 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 429 PetscFunctionReturn(0); 430 } 431 432 PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 433 { 434 PetscErrorCode ierr; 435 436 PetscFunctionBegin; 437 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 441 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 442 void (**funcs)(PetscInt, PetscInt, PetscInt, 443 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 444 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 445 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 446 InsertMode mode, Vec localX) 447 { 448 PetscErrorCode ierr; 449 450 PetscFunctionBegin; 451 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 452 PetscFunctionReturn(0); 453 } 454 455 PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 456 void (**funcs)(PetscInt, PetscInt, PetscInt, 457 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 458 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 459 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 460 InsertMode mode, Vec localX) 461 { 462 PetscErrorCode ierr; 463 464 PetscFunctionBegin; 465 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 466 PetscFunctionReturn(0); 467 } 468