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 static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints) 163 { 164 PetscReal *points; 165 PetscInt f, numPoints; 166 PetscErrorCode ierr; 167 168 PetscFunctionBegin; 169 numPoints = 0; 170 for (f = 0; f < Nf; ++f) { 171 if (funcs[f]) { 172 PetscQuadrature fAllPoints; 173 PetscInt fNumPoints; 174 175 ierr = PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);CHKERRQ(ierr); 176 ierr = PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);CHKERRQ(ierr); 177 numPoints += fNumPoints; 178 } 179 } 180 ierr = PetscMalloc1(dim*numPoints,&points);CHKERRQ(ierr); 181 numPoints = 0; 182 for (f = 0; f < Nf; ++f) { 183 PetscInt spDim; 184 185 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 186 if (funcs[f]) { 187 PetscQuadrature fAllPoints; 188 PetscInt fNumPoints, q; 189 const PetscReal *fPoints; 190 191 ierr = PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);CHKERRQ(ierr); 192 ierr = PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, &fPoints, NULL);CHKERRQ(ierr); 193 for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q]; 194 numPoints += fNumPoints; 195 } 196 } 197 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);CHKERRQ(ierr); 198 ierr = PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 202 /* 203 This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM 204 205 There are several different scenarios: 206 207 1) Volumetric mesh with volumetric auxiliary data 208 209 Here minHeight=0 since we loop over cells. 210 211 2) Boundary mesh with boundary auxiliary data 212 213 Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation. 214 215 3) Volumetric mesh with boundary auxiliary data 216 217 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. 218 219 The maxHeight is used to support enforcement of constraints in DMForest. 220 221 If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it. 222 223 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. 224 - 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 225 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 226 dual spaces for the boundary from our input spaces. 227 - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them. 228 229 We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration 230 231 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. 232 */ 233 static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, 234 PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[], 235 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, 236 InsertMode mode, Vec localX) 237 { 238 DM dmAux = NULL; 239 PetscDS prob, probAux = NULL; 240 Vec localA = NULL; 241 PetscSection section; 242 PetscDualSpace *sp, *cellsp; 243 PetscReal **basisTab = NULL, **basisDerTab = NULL, **basisTabAux = NULL, **basisDerTabAux = NULL; 244 PetscInt *Nc; 245 PetscInt dim, dimEmbed, minHeight, maxHeight, h, Nf, NfAux = 0, f; 246 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE; 247 PetscErrorCode ierr; 248 249 PetscFunctionBegin; 250 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 251 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);CHKERRQ(ierr); 252 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 253 ierr = DMPlexGetVTKCellHeight(dm, &minHeight);CHKERRQ(ierr); 254 /* Auxiliary information can only be used with interpolation of field functions */ 255 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 256 if (!minHeight && dmAux) { 257 DMLabel spmap; 258 259 /* If dmAux is a surface, then force the projection to take place over a surface */ 260 ierr = DMPlexGetSubpointMap(dmAux, &spmap);CHKERRQ(ierr); 261 if (spmap) {minHeight = 1; auxBd = PETSC_TRUE;} 262 } 263 } 264 ierr = DMPlexGetMaxProjectionHeight(dm, &maxHeight);CHKERRQ(ierr); 265 maxHeight = PetscMax(maxHeight, minHeight); 266 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %d not in [0, %d)\n", maxHeight, dim);} 267 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 268 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 269 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 270 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 271 if (dmAux) { 272 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 273 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 274 } 275 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 276 ierr = PetscMalloc2(Nf, &isFE, Nf, &sp);CHKERRQ(ierr); 277 if (maxHeight > 0) {ierr = PetscMalloc1(Nf, &cellsp);CHKERRQ(ierr);} 278 else {cellsp = sp;} 279 if (localU && localU != localX) {ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);CHKERRQ(ierr);} 280 /* Get cell dual spaces */ 281 for (f = 0; f < Nf; ++f) { 282 PetscObject obj; 283 PetscClassId id; 284 285 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 286 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 287 if (id == PETSCFE_CLASSID) { 288 PetscFE fe = (PetscFE) obj; 289 290 hasFE = PETSC_TRUE; 291 isFE[f] = PETSC_TRUE; 292 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 293 } else if (id == PETSCFV_CLASSID) { 294 PetscFV fv = (PetscFV) obj; 295 296 hasFV = PETSC_TRUE; 297 isFE[f] = PETSC_FALSE; 298 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 299 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 300 } 301 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 302 PetscInt effectiveHeight = auxBd ? minHeight : 0; 303 PetscFE fem, subfem; 304 PetscReal *points; 305 PetscInt numPoints, spDim, qdim = 0, d; 306 307 if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation"); 308 numPoints = 0; 309 for (f = 0; f < Nf; ++f) { 310 if (!effectiveHeight) {sp[f] = cellsp[f];} 311 else {ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr);} 312 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 313 for (d = 0; d < spDim; ++d) { 314 if (funcs[f]) { 315 PetscQuadrature quad; 316 PetscInt Nq; 317 318 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 319 ierr = PetscQuadratureGetData(quad, &qdim, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 320 numPoints += Nq; 321 } 322 } 323 } 324 ierr = PetscMalloc1(numPoints*qdim, &points);CHKERRQ(ierr); 325 numPoints = 0; 326 for (f = 0; f < Nf; ++f) { 327 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 328 for (d = 0; d < spDim; ++d) { 329 if (funcs[f]) { 330 PetscQuadrature quad; 331 const PetscReal *qpoints; 332 PetscInt Nq, q; 333 334 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 335 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL);CHKERRQ(ierr); 336 for (q = 0; q < Nq*qdim; ++q) points[numPoints*qdim+q] = qpoints[q]; 337 numPoints += Nq; 338 } 339 } 340 } 341 ierr = PetscMalloc4(Nf, &basisTab, Nf, &basisDerTab, NfAux, &basisTabAux, NfAux, &basisDerTabAux);CHKERRQ(ierr); 342 for (f = 0; f < Nf; ++f) { 343 if (!isFE[f]) continue; 344 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 345 if (!effectiveHeight) {subfem = fem;} 346 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 347 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 348 } 349 for (f = 0; f < NfAux; ++f) { 350 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 351 if (!effectiveHeight || auxBd) {subfem = fem;} 352 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 353 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 354 } 355 ierr = PetscFree(points);CHKERRQ(ierr); 356 } 357 /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */ 358 for (h = minHeight; h <= maxHeight; h++) { 359 PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight); 360 PetscDS probEff = prob; 361 PetscScalar *values; 362 PetscBool *fieldActive; 363 PetscInt pStart, pEnd, p, spDim, totDim, numValues; 364 365 if (effectiveHeight) {ierr = PetscDSGetHeightSubspace(prob, effectiveHeight, &probEff);CHKERRQ(ierr);} 366 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 367 if (!h) { 368 PetscInt cEndInterior; 369 370 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 371 pEnd = cEndInterior < 0 ? pEnd : cEndInterior; 372 } 373 if (pEnd <= pStart) continue; 374 /* Compute totDim, the number of dofs in the closure of a point at this height */ 375 totDim = 0; 376 for (f = 0; f < Nf; ++f) { 377 if (!effectiveHeight) { 378 sp[f] = cellsp[f]; 379 } else { 380 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr); 381 if (!sp[f]) continue; 382 } 383 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 384 totDim += spDim; 385 } 386 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 387 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim); 388 if (!totDim) continue; 389 /* Loop over points at this height */ 390 ierr = DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 391 ierr = DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 392 for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 393 if (label) { 394 PetscInt i; 395 396 for (i = 0; i < numIds; ++i) { 397 IS pointIS; 398 const PetscInt *points; 399 PetscInt n; 400 401 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 402 if (!pointIS) continue; /* No points with that id on this process */ 403 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 404 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 405 for (p = 0; p < n; ++p) { 406 const PetscInt point = points[p]; 407 408 if ((point < pStart) || (point >= pEnd)) continue; 409 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 410 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); 411 if (ierr) { 412 PetscErrorCode ierr2; 413 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 414 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 415 CHKERRQ(ierr); 416 } 417 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, Ncc, comps, values, mode);CHKERRQ(ierr); 418 } 419 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 420 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 421 } 422 } else { 423 for (p = pStart; p < pEnd; ++p) { 424 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 425 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); 426 if (ierr) { 427 PetscErrorCode ierr2; 428 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 429 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 430 CHKERRQ(ierr); 431 } 432 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, Ncc, comps, values, mode);CHKERRQ(ierr); 433 } 434 } 435 ierr = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 436 ierr = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 437 } 438 /* Cleanup */ 439 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 440 PetscInt effectiveHeight = auxBd ? minHeight : 0; 441 PetscFE fem, subfem; 442 443 for (f = 0; f < Nf; ++f) { 444 if (!isFE[f]) continue; 445 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 446 if (!effectiveHeight) {subfem = fem;} 447 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 448 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 449 } 450 for (f = 0; f < NfAux; ++f) { 451 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 452 if (!effectiveHeight || auxBd) {subfem = fem;} 453 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 454 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 455 } 456 ierr = PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);CHKERRQ(ierr); 457 } 458 ierr = PetscFree2(isFE, sp);CHKERRQ(ierr); 459 if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 460 PetscFunctionReturn(0); 461 } 462 463 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 464 { 465 PetscErrorCode ierr; 466 467 PetscFunctionBegin; 468 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 469 PetscFunctionReturn(0); 470 } 471 472 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) 473 { 474 PetscErrorCode ierr; 475 476 PetscFunctionBegin; 477 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 478 PetscFunctionReturn(0); 479 } 480 481 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 482 void (**funcs)(PetscInt, PetscInt, PetscInt, 483 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 484 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 485 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 486 InsertMode mode, Vec localX) 487 { 488 PetscErrorCode ierr; 489 490 PetscFunctionBegin; 491 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 492 PetscFunctionReturn(0); 493 } 494 495 PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 496 void (**funcs)(PetscInt, PetscInt, PetscInt, 497 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 498 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 499 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 500 InsertMode mode, Vec localX) 501 { 502 PetscErrorCode ierr; 503 504 PetscFunctionBegin; 505 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 506 PetscFunctionReturn(0); 507 } 508