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