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, transform; 11 PetscErrorCode ierr; 12 13 PetscFunctionBeginHot; 14 ierr = DMGetCoordinateDim(dm,&coordDim);CHKERRQ(ierr); 15 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 16 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 17 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 18 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 19 /* Get values for closure */ 20 isAffine = fegeom->isAffine; 21 for (f = 0, v = 0, tp = 0; f < Nf; ++f) { 22 void * const ctx = ctxs ? ctxs[f] : NULL; 23 24 if (!sp[f]) continue; 25 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 26 if (funcs[f]) { 27 if (isFE[f]) { 28 PetscQuadrature allPoints; 29 PetscInt q, dim, numPoints; 30 const PetscReal *points; 31 PetscScalar *pointEval; 32 PetscReal *x; 33 DM rdm; 34 35 ierr = PetscDualSpaceGetDM(sp[f],&rdm);CHKERRQ(ierr); 36 ierr = PetscDualSpaceGetAllPoints(sp[f], &allPoints);CHKERRQ(ierr); 37 ierr = PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);CHKERRQ(ierr); 38 ierr = DMGetWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 39 ierr = DMGetWorkArray(rdm,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 if (transform) {ierr = DMPlexBasisTransformApplyReal_Internal(dm, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx);CHKERRQ(ierr); v0 = x;} 50 ierr = (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f]*q], ctx);CHKERRQ(ierr); 51 } 52 ierr = PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);CHKERRQ(ierr); 53 ierr = DMRestoreWorkArray(rdm,coordDim,MPIU_REAL,&x);CHKERRQ(ierr); 54 ierr = DMRestoreWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 55 v += spDim; 56 } else { 57 for (d = 0; d < spDim; ++d, ++v) { 58 ierr = PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr); 59 } 60 } 61 } else { 62 for (d = 0; d < spDim; d++, v++) {values[v] = 0.0;} 63 } 64 } 65 PetscFunctionReturn(0); 66 } 67 68 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[], 69 PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 70 void (**funcs)(PetscInt, PetscInt, PetscInt, 71 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 72 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 73 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, 74 PetscScalar values[]) 75 { 76 PetscSection section, sectionAux = NULL; 77 PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *refSpaceDer, *refSpaceDerAux = NULL, *bc; 78 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 79 PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL; 80 const PetscScalar *constants; 81 PetscReal *x; 82 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 83 const PetscInt dE = fegeom->dimEmbed; 84 PetscInt dimAux = 0, numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0; 85 PetscBool isAffine, transform; 86 PetscErrorCode ierr; 87 88 PetscFunctionBeginHot; 89 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 90 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 91 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 92 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 93 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 94 ierr = PetscDSGetEvaluationArrays(prob, &u, &bc /*&u_t*/, &u_x);CHKERRQ(ierr); 95 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 96 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 97 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 98 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 99 ierr = DMPlexVecGetClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 100 if (dmAux) { 101 PetscInt subp; 102 103 ierr = DMPlexGetAuxiliaryPoint(dm, dmAux, p, &subp);CHKERRQ(ierr); 104 ierr = PetscDSGetSpatialDimension(probAux, &dimAux);CHKERRQ(ierr); 105 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 106 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 107 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 108 ierr = DMGetSection(dmAux, §ionAux);CHKERRQ(ierr); 109 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 110 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 111 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);CHKERRQ(ierr); 112 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 113 ierr = DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);CHKERRQ(ierr); 114 } 115 /* Get values for closure */ 116 isAffine = fegeom->isAffine; 117 for (f = 0, v = 0; f < Nf; ++f) { 118 PetscQuadrature allPoints; 119 PetscInt q, dim, numPoints; 120 const PetscReal *points; 121 PetscScalar *pointEval; 122 DM dm; 123 124 if (!sp[f]) continue; 125 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 126 if (!funcs[f]) { 127 for (d = 0; d < spDim; d++, v++) values[v] = 0.; 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(dE, 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 if (transform) {ierr = DMPlexBasisTransformApplyReal_Internal(dm, v0, PETSC_TRUE, dE, v0, x, dm->transformCtx);CHKERRQ(ierr); v0 = x;} 149 (*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]); 150 } 151 ierr = PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);CHKERRQ(ierr); 152 ierr = DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 153 v += spDim; 154 } 155 ierr = DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 156 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);CHKERRQ(ierr);} 157 PetscFunctionReturn(0); 158 } 159 160 static PetscErrorCode DMProjectPoint_BdField_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[], 161 PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 162 void (**funcs)(PetscInt, PetscInt, PetscInt, 163 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 164 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 165 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, 166 PetscScalar values[]) 167 { 168 PetscSection section, sectionAux = NULL; 169 PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *refSpaceDer, *refSpaceDerAux = NULL, *bc; 170 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 171 PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL; 172 const PetscScalar *constants; 173 PetscReal *x; 174 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 175 const PetscInt dE = fegeom->dimEmbed; 176 PetscInt dimAux = 0, numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0; 177 PetscBool isAffine; 178 PetscErrorCode ierr; 179 180 PetscFunctionBeginHot; 181 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 182 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 183 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 184 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 185 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 186 ierr = PetscDSGetEvaluationArrays(prob, &u, &bc /*&u_t*/, &u_x);CHKERRQ(ierr); 187 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 188 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 189 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 190 ierr = DMPlexVecGetClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 191 if (dmAux) { 192 DMLabel spmap; 193 PetscInt subp = p; 194 195 /* If dm is a submesh, do not get subpoint */ 196 ierr = DMPlexGetSubpointMap(dm, &spmap);CHKERRQ(ierr); 197 if (!spmap) {ierr = DMPlexGetSubpoint(dmAux, p, &subp);CHKERRQ(ierr);} 198 ierr = PetscDSGetSpatialDimension(probAux, &dimAux);CHKERRQ(ierr); 199 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 200 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 201 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 202 ierr = DMGetSection(dmAux, §ionAux);CHKERRQ(ierr); 203 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 204 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 205 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);CHKERRQ(ierr); 206 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 207 ierr = DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);CHKERRQ(ierr); 208 } 209 /* Get values for closure */ 210 isAffine = fegeom->isAffine; 211 for (f = 0, v = 0; f < Nf; ++f) { 212 PetscQuadrature allPoints; 213 PetscInt q, dim, numPoints; 214 const PetscReal *points; 215 PetscScalar *pointEval; 216 DM dm; 217 218 if (!sp[f]) continue; 219 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 220 if (!funcs[f]) { 221 for (d = 0; d < spDim; d++, v++) values[v] = 0.; 222 continue; 223 } 224 ierr = PetscDualSpaceGetDM(sp[f],&dm);CHKERRQ(ierr); 225 ierr = PetscDualSpaceGetAllPoints(sp[f], &allPoints);CHKERRQ(ierr); 226 ierr = PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);CHKERRQ(ierr); 227 ierr = DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 228 for (q = 0; q < numPoints; ++q, ++tp) { 229 const PetscReal *v0; 230 const PetscReal *invJ; 231 const PetscReal *n; 232 233 if (isAffine) { 234 CoordinatesRefToReal(dE, dim, fegeom->xi, fegeom->v, fegeom->J, &points[q*dim], x); 235 v0 = x; 236 invJ = fegeom->invJ; 237 n = fegeom->n; 238 } else { 239 v0 = &fegeom->v[tp*dE]; 240 invJ = &fegeom->invJ[tp*dE*dE]; 241 n = &fegeom->n[tp*dE]; 242 } 243 EvaluateFieldJets(dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, refSpaceDer, invJ, coefficients, coefficients_t, u, u_x, u_t); 244 if (probAux) {EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, tp, basisTabAux, basisDerTabAux, refSpaceDerAux, invJ, coefficientsAux, coefficientsAux_t, a, a_x, a_t);} 245 (*funcs[f])(dE, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, v0, n, numConstants, constants, &pointEval[Nc[f]*q]); 246 } 247 ierr = PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);CHKERRQ(ierr); 248 ierr = DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 249 v += spDim; 250 } 251 ierr = DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 252 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);CHKERRQ(ierr);} 253 PetscFunctionReturn(0); 254 } 255 256 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[], 257 PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[], 258 PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 259 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[]) 260 { 261 PetscFVCellGeom fvgeom; 262 PetscInt dim, dimEmbed; 263 PetscErrorCode ierr; 264 265 PetscFunctionBeginHot; 266 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 267 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 268 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 269 switch (type) { 270 case DM_BC_ESSENTIAL: 271 case DM_BC_NATURAL: 272 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; 273 case DM_BC_ESSENTIAL_FIELD: 274 case DM_BC_NATURAL_FIELD: 275 ierr = DMProjectPoint_Field_Private(dm, prob, dmAux, probAux, time, localU, localA, fegeom, sp, p, Ncc, comps, 276 basisTab, basisDerTab, basisTabAux, basisDerTabAux, 277 (void (**)(PetscInt, PetscInt, PetscInt, 278 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 279 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 280 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 281 case DM_BC_ESSENTIAL_BD_FIELD: 282 ierr = DMProjectPoint_BdField_Private(dm, prob, dmAux, probAux, time, localU, localA, fegeom, sp, p, Ncc, comps, 283 basisTab, basisDerTab, basisTabAux, basisDerTabAux, 284 (void (**)(PetscInt, PetscInt, PetscInt, 285 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 286 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 287 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 288 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type); 289 } 290 PetscFunctionReturn(0); 291 } 292 293 static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints) 294 { 295 PetscReal *points; 296 PetscInt f, numPoints; 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 numPoints = 0; 301 for (f = 0; f < Nf; ++f) { 302 if (funcs[f]) { 303 PetscQuadrature fAllPoints; 304 PetscInt fNumPoints; 305 306 ierr = PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);CHKERRQ(ierr); 307 ierr = PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);CHKERRQ(ierr); 308 numPoints += fNumPoints; 309 } 310 } 311 ierr = PetscMalloc1(dim*numPoints,&points);CHKERRQ(ierr); 312 numPoints = 0; 313 for (f = 0; f < Nf; ++f) { 314 PetscInt spDim; 315 316 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 317 if (funcs[f]) { 318 PetscQuadrature fAllPoints; 319 PetscInt qdim, fNumPoints, q; 320 const PetscReal *fPoints; 321 322 ierr = PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);CHKERRQ(ierr); 323 ierr = PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);CHKERRQ(ierr); 324 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); 325 for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q]; 326 numPoints += fNumPoints; 327 } 328 } 329 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);CHKERRQ(ierr); 330 ierr = PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);CHKERRQ(ierr); 331 PetscFunctionReturn(0); 332 } 333 334 static PetscErrorCode DMGetFirstLabelEntry_Private(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *lStart, PetscDS *prob) 335 { 336 DMLabel depthLabel; 337 PetscInt dim, cdepth, ls = -1, i; 338 PetscErrorCode ierr; 339 340 PetscFunctionBegin; 341 if (lStart) *lStart = -1; 342 if (!label) PetscFunctionReturn(0); 343 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 344 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 345 cdepth = dim - height; 346 for (i = 0; i < numIds; ++i) { 347 IS pointIS; 348 const PetscInt *points; 349 PetscInt pdepth; 350 351 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 352 if (!pointIS) continue; /* No points with that id on this process */ 353 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 354 ierr = DMLabelGetValue(depthLabel, points[0], &pdepth);CHKERRQ(ierr); 355 if (pdepth == cdepth) { 356 ls = points[0]; 357 if (prob) {ierr = DMGetCellDS(dm, ls, prob);CHKERRQ(ierr);} 358 } 359 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 360 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 361 if (ls >= 0) break; 362 } 363 if (lStart) *lStart = ls; 364 PetscFunctionReturn(0); 365 } 366 367 /* 368 This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM 369 370 There are several different scenarios: 371 372 1) Volumetric mesh with volumetric auxiliary data 373 374 Here minHeight=0 since we loop over cells. 375 376 2) Boundary mesh with boundary auxiliary data 377 378 Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation. 379 380 3) Volumetric mesh with boundary auxiliary data 381 382 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. 383 384 The maxHeight is used to support enforcement of constraints in DMForest. 385 386 If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it. 387 388 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. 389 - 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 390 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 391 dual spaces for the boundary from our input spaces. 392 - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them. 393 394 We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration 395 396 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. 397 */ 398 static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, 399 PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[], 400 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, 401 InsertMode mode, Vec localX) 402 { 403 DM dmAux = NULL, tdm; 404 PetscDS prob = NULL, probAux = NULL; 405 Vec localA = NULL, tv; 406 PetscSection section; 407 PetscDualSpace *sp, *cellsp; 408 PetscReal **basisTab = NULL, **basisDerTab = NULL, **basisTabAux = NULL, **basisDerTabAux = NULL; 409 PetscInt *Nc; 410 PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, Nf, NfAux = 0, f; 411 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE, transform; 412 DMField coordField; 413 DMLabel depthLabel; 414 PetscQuadrature allPoints = NULL; 415 PetscErrorCode ierr; 416 417 PetscFunctionBegin; 418 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 419 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);CHKERRQ(ierr); 420 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 421 ierr = DMPlexGetVTKCellHeight(dm, &minHeight);CHKERRQ(ierr); 422 ierr = DMGetBasisTransformDM_Internal(dm, &tdm);CHKERRQ(ierr); 423 ierr = DMGetBasisTransformVec_Internal(dm, &tv);CHKERRQ(ierr); 424 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 425 /* Auxiliary information can only be used with interpolation of field functions */ 426 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 427 if (dmAux && !localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector"); 428 if (!minHeight && dmAux) { 429 DMLabel spmap; 430 431 /* If dmAux is a surface, then force the projection to take place over a surface */ 432 ierr = DMPlexGetSubpointMap(dmAux, &spmap);CHKERRQ(ierr); 433 if (spmap) { 434 ierr = DMPlexGetVTKCellHeight(dmAux, &minHeight);CHKERRQ(ierr); 435 auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE; 436 } 437 } 438 } 439 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 440 ierr = DMPlexGetDepthLabel(dm,&depthLabel);CHKERRQ(ierr); 441 ierr = DMPlexGetMaxProjectionHeight(dm, &maxHeight);CHKERRQ(ierr); 442 maxHeight = PetscMax(maxHeight, minHeight); 443 if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim); 444 ierr = DMGetFirstLabelEntry_Private(dm, label, numIds, ids, 0, NULL, &prob);CHKERRQ(ierr); 445 if (!prob) {ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);} 446 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 447 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 448 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 449 if (dmAux) { 450 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 451 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 452 } 453 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 454 ierr = PetscMalloc2(Nf, &isFE, Nf, &sp);CHKERRQ(ierr); 455 if (maxHeight > 0) {ierr = PetscMalloc1(Nf, &cellsp);CHKERRQ(ierr);} 456 else {cellsp = sp;} 457 if (localU && localU != localX) {ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);CHKERRQ(ierr);} 458 /* Get cell dual spaces */ 459 for (f = 0; f < Nf; ++f) { 460 PetscObject obj; 461 PetscClassId id; 462 463 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 464 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 465 if (id == PETSCFE_CLASSID) { 466 PetscFE fe = (PetscFE) obj; 467 468 hasFE = PETSC_TRUE; 469 isFE[f] = PETSC_TRUE; 470 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 471 } else if (id == PETSCFV_CLASSID) { 472 PetscFV fv = (PetscFV) obj; 473 474 hasFV = PETSC_TRUE; 475 isFE[f] = PETSC_FALSE; 476 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 477 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 478 } 479 ierr = DMGetCoordinateField(dm,&coordField);CHKERRQ(ierr); 480 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 481 PetscInt effectiveHeight = auxBd ? minHeight : 0; 482 PetscFE fem, subfem; 483 const PetscReal *points; 484 PetscInt numPoints; 485 486 if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation"); 487 for (f = 0; f < Nf; ++f) { 488 if (!effectiveHeight) {sp[f] = cellsp[f];} 489 else {ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr);} 490 } 491 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);CHKERRQ(ierr); 492 ierr = PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);CHKERRQ(ierr); 493 ierr = PetscMalloc4(Nf, &basisTab, Nf, &basisDerTab, NfAux, &basisTabAux, NfAux, &basisDerTabAux);CHKERRQ(ierr); 494 for (f = 0; f < Nf; ++f) { 495 if (!isFE[f]) continue; 496 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 497 if (!effectiveHeight) {subfem = fem;} 498 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 499 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 500 } 501 for (f = 0; f < NfAux; ++f) { 502 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 503 if (!effectiveHeight || auxBd) {subfem = fem;} 504 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 505 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 506 } 507 } 508 /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */ 509 for (h = minHeight; h <= maxHeight; h++) { 510 PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight); 511 PetscDS probEff = prob; 512 PetscScalar *values; 513 PetscBool *fieldActive; 514 PetscInt maxDegree; 515 PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues; 516 IS heightIS; 517 518 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 519 ierr = DMGetFirstLabelEntry_Private(dm, label, numIds, ids, h, &lStart, NULL);CHKERRQ(ierr); 520 ierr = DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);CHKERRQ(ierr); 521 if (!h) { 522 PetscInt cEndInterior; 523 524 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 525 pEnd = cEndInterior < 0 ? pEnd : cEndInterior; 526 } 527 if (pEnd <= pStart) { 528 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 529 continue; 530 } 531 /* Compute totDim, the number of dofs in the closure of a point at this height */ 532 totDim = 0; 533 for (f = 0; f < Nf; ++f) { 534 if (!effectiveHeight) { 535 sp[f] = cellsp[f]; 536 } else { 537 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr); 538 if (!sp[f]) continue; 539 } 540 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 541 totDim += spDim; 542 } 543 ierr = DMPlexVecGetClosure(dm, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);CHKERRQ(ierr); 544 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim); 545 if (!totDim) { 546 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 547 continue; 548 } 549 if (effectiveHeight) {ierr = PetscDSGetHeightSubspace(prob, effectiveHeight, &probEff);CHKERRQ(ierr);} 550 /* Loop over points at this height */ 551 ierr = DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 552 ierr = DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 553 for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 554 if (label) { 555 PetscInt i; 556 557 for (i = 0; i < numIds; ++i) { 558 IS pointIS, isectIS; 559 const PetscInt *points; 560 PetscInt n; 561 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 562 PetscQuadrature quad = NULL; 563 564 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 565 if (!pointIS) continue; /* No points with that id on this process */ 566 ierr = ISIntersect(pointIS,heightIS,&isectIS);CHKERRQ(ierr); 567 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 568 if (!isectIS) continue; 569 ierr = ISGetLocalSize(isectIS, &n);CHKERRQ(ierr); 570 ierr = ISGetIndices(isectIS, &points);CHKERRQ(ierr); 571 ierr = DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);CHKERRQ(ierr); 572 if (maxDegree <= 1) { 573 ierr = DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);CHKERRQ(ierr); 574 } 575 if (!quad) { 576 if (!h && allPoints) { 577 quad = allPoints; 578 allPoints = NULL; 579 } else { 580 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);CHKERRQ(ierr); 581 } 582 } 583 ierr = DMFieldCreateFEGeom(coordField,isectIS,quad,PETSC_FALSE,&fegeom);CHKERRQ(ierr); 584 for (p = 0; p < n; ++p) { 585 const PetscInt point = points[p]; 586 587 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 588 ierr = PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 589 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); 590 if (ierr) { 591 PetscErrorCode ierr2; 592 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 593 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 594 CHKERRQ(ierr); 595 } 596 if (transform) {ierr = DMPlexBasisTransformPoint_Internal(dm, tdm, tv, point, fieldActive, PETSC_FALSE, values);CHKERRQ(ierr);} 597 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, Ncc, comps, values, mode);CHKERRQ(ierr); 598 } 599 ierr = PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 600 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 601 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 602 ierr = ISRestoreIndices(isectIS, &points);CHKERRQ(ierr); 603 ierr = ISDestroy(&isectIS);CHKERRQ(ierr); 604 } 605 } else { 606 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 607 PetscQuadrature quad = NULL; 608 IS pointIS; 609 610 ierr = ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);CHKERRQ(ierr); 611 ierr = DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);CHKERRQ(ierr); 612 if (maxDegree <= 1) { 613 ierr = DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);CHKERRQ(ierr); 614 } 615 if (!quad) { 616 if (!h && allPoints) { 617 quad = allPoints; 618 allPoints = NULL; 619 } else { 620 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);CHKERRQ(ierr); 621 } 622 } 623 ierr = DMFieldCreateFEGeom(coordField,pointIS,quad,PETSC_FALSE,&fegeom);CHKERRQ(ierr); 624 for (p = pStart; p < pEnd; ++p) { 625 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 626 ierr = PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);CHKERRQ(ierr); 627 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); 628 if (ierr) { 629 PetscErrorCode ierr2; 630 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 631 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 632 CHKERRQ(ierr); 633 } 634 if (transform) {ierr = DMPlexBasisTransformPoint_Internal(dm, tdm, tv, p, fieldActive, PETSC_FALSE, values);CHKERRQ(ierr);} 635 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, Ncc, comps, values, mode);CHKERRQ(ierr); 636 } 637 ierr = PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);CHKERRQ(ierr); 638 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 639 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 640 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 641 } 642 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 643 ierr = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 644 ierr = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 645 } 646 /* Cleanup */ 647 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 648 PetscInt effectiveHeight = auxBd ? minHeight : 0; 649 PetscFE fem, subfem; 650 651 for (f = 0; f < Nf; ++f) { 652 if (!isFE[f]) continue; 653 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 654 if (!effectiveHeight) {subfem = fem;} 655 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 656 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 657 } 658 for (f = 0; f < NfAux; ++f) { 659 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 660 if (!effectiveHeight || auxBd) {subfem = fem;} 661 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 662 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 663 } 664 ierr = PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);CHKERRQ(ierr); 665 } 666 ierr = PetscQuadratureDestroy(&allPoints);CHKERRQ(ierr); 667 ierr = PetscFree2(isFE, sp);CHKERRQ(ierr); 668 if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 669 PetscFunctionReturn(0); 670 } 671 672 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 673 { 674 PetscErrorCode ierr; 675 676 PetscFunctionBegin; 677 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 678 PetscFunctionReturn(0); 679 } 680 681 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) 682 { 683 PetscErrorCode ierr; 684 685 PetscFunctionBegin; 686 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 687 PetscFunctionReturn(0); 688 } 689 690 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 691 void (**funcs)(PetscInt, PetscInt, PetscInt, 692 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 693 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 694 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 695 InsertMode mode, Vec localX) 696 { 697 PetscErrorCode ierr; 698 699 PetscFunctionBegin; 700 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 701 PetscFunctionReturn(0); 702 } 703 704 PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 705 void (**funcs)(PetscInt, PetscInt, PetscInt, 706 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 707 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 708 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 709 InsertMode mode, Vec localX) 710 { 711 PetscErrorCode ierr; 712 713 PetscFunctionBegin; 714 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 715 PetscFunctionReturn(0); 716 } 717