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