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