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; 370 PetscInt maxDegree; 371 PetscInt pStart, pEnd, p, spDim, totDim, numValues; 372 IS heightIS; 373 374 if (effectiveHeight) {ierr = PetscDSGetHeightSubspace(prob, effectiveHeight, &probEff);CHKERRQ(ierr);} 375 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 376 ierr = DMLabelGetStratumIS(depthLabel,depth - h,&heightIS);CHKERRQ(ierr); 377 if (!h) { 378 PetscInt cEndInterior; 379 380 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 381 pEnd = cEndInterior < 0 ? pEnd : cEndInterior; 382 } 383 if (pEnd <= pStart) { 384 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 385 continue; 386 } 387 /* Compute totDim, the number of dofs in the closure of a point at this height */ 388 totDim = 0; 389 for (f = 0; f < Nf; ++f) { 390 if (!effectiveHeight) { 391 sp[f] = cellsp[f]; 392 } else { 393 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr); 394 if (!sp[f]) continue; 395 } 396 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 397 totDim += spDim; 398 } 399 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 400 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim); 401 if (!totDim) { 402 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 403 continue; 404 } 405 /* Loop over points at this height */ 406 ierr = DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 407 ierr = DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 408 for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 409 if (label) { 410 PetscInt i; 411 412 for (i = 0; i < numIds; ++i) { 413 IS pointIS, isectIS; 414 const PetscInt *points; 415 PetscInt n; 416 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 417 PetscQuadrature quad = NULL; 418 419 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 420 if (!pointIS) continue; /* No points with that id on this process */ 421 ierr = ISIntersect(pointIS,heightIS,&isectIS);CHKERRQ(ierr); 422 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 423 if (!isectIS) continue; 424 ierr = ISGetLocalSize(isectIS, &n);CHKERRQ(ierr); 425 ierr = ISGetIndices(isectIS, &points);CHKERRQ(ierr); 426 ierr = DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);CHKERRQ(ierr); 427 if (maxDegree <= 1) { 428 ierr = DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);CHKERRQ(ierr); 429 } 430 if (!quad) { 431 if (!h && allPoints) { 432 quad = allPoints; 433 allPoints = NULL; 434 } else { 435 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);CHKERRQ(ierr); 436 } 437 } 438 ierr = DMFieldCreateFEGeom(coordField,isectIS,quad,PETSC_FALSE,&fegeom);CHKERRQ(ierr); 439 for (p = 0; p < n; ++p) { 440 const PetscInt point = points[p]; 441 442 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 443 ierr = PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 444 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); 445 if (ierr) { 446 PetscErrorCode ierr2; 447 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 448 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 449 CHKERRQ(ierr); 450 } 451 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, Ncc, comps, values, mode);CHKERRQ(ierr); 452 } 453 ierr = PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 454 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 455 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 456 ierr = ISRestoreIndices(isectIS, &points);CHKERRQ(ierr); 457 ierr = ISDestroy(&isectIS);CHKERRQ(ierr); 458 } 459 } else { 460 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 461 PetscQuadrature quad = NULL; 462 IS pointIS; 463 464 ierr = ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);CHKERRQ(ierr); 465 ierr = DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);CHKERRQ(ierr); 466 if (maxDegree <= 1) { 467 ierr = DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);CHKERRQ(ierr); 468 } 469 if (!quad) { 470 if (!h && allPoints) { 471 quad = allPoints; 472 allPoints = NULL; 473 } else { 474 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);CHKERRQ(ierr); 475 } 476 } 477 ierr = DMFieldCreateFEGeom(coordField,pointIS,quad,PETSC_FALSE,&fegeom);CHKERRQ(ierr); 478 for (p = pStart; p < pEnd; ++p) { 479 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 480 ierr = PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);CHKERRQ(ierr); 481 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); 482 if (ierr) { 483 PetscErrorCode ierr2; 484 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 485 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 486 CHKERRQ(ierr); 487 } 488 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, Ncc, comps, values, mode);CHKERRQ(ierr); 489 } 490 ierr = PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);CHKERRQ(ierr); 491 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 492 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 493 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 494 } 495 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 496 ierr = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 497 ierr = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 498 } 499 /* Cleanup */ 500 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 501 PetscInt effectiveHeight = auxBd ? minHeight : 0; 502 PetscFE fem, subfem; 503 504 for (f = 0; f < Nf; ++f) { 505 if (!isFE[f]) continue; 506 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 507 if (!effectiveHeight) {subfem = fem;} 508 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 509 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 510 } 511 for (f = 0; f < NfAux; ++f) { 512 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 513 if (!effectiveHeight || auxBd) {subfem = fem;} 514 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 515 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 516 } 517 ierr = PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);CHKERRQ(ierr); 518 } 519 ierr = PetscQuadratureDestroy(&allPoints);CHKERRQ(ierr); 520 ierr = PetscFree2(isFE, sp);CHKERRQ(ierr); 521 if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 522 PetscFunctionReturn(0); 523 } 524 525 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, 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, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 531 PetscFunctionReturn(0); 532 } 533 534 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) 535 { 536 PetscErrorCode ierr; 537 538 PetscFunctionBegin; 539 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 540 PetscFunctionReturn(0); 541 } 542 543 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 544 void (**funcs)(PetscInt, PetscInt, PetscInt, 545 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 546 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 547 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 548 InsertMode mode, Vec localX) 549 { 550 PetscErrorCode ierr; 551 552 PetscFunctionBegin; 553 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 554 PetscFunctionReturn(0); 555 } 556 557 PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 558 void (**funcs)(PetscInt, PetscInt, PetscInt, 559 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 560 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 561 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 562 InsertMode mode, Vec localX) 563 { 564 PetscErrorCode ierr; 565 566 PetscFunctionBegin; 567 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 568 PetscFunctionReturn(0); 569 } 570