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