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