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