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