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