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