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