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 = EvaluateFieldJets(prob, dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, &fegeom, coefficients, coefficients_t, u, u_x, u_t);CHKERRQ(ierr); 153 if (probAux) {ierr = EvaluateFieldJets(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.invJ = &fgeom->invJ[tp*dE*dE]; 249 fegeom.detJ = &fgeom->detJ[tp]; 250 fegeom.n = &fgeom->n[tp*dE]; 251 } 252 ierr = EvaluateFieldJets(prob, dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, &fegeom, coefficients, coefficients_t, u, u_x, u_t);CHKERRQ(ierr); 253 if (probAux) {ierr = EvaluateFieldJets(probAux, dimAux, NfAux, NbAux, NcAux, tp, basisTabAux, basisDerTabAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);CHKERRQ(ierr);} 254 (*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]); 255 } 256 ierr = PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);CHKERRQ(ierr); 257 ierr = DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 258 v += spDim; 259 } 260 ierr = DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 261 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);CHKERRQ(ierr);} 262 PetscFunctionReturn(0); 263 } 264 265 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[], 266 PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[], 267 PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 268 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[]) 269 { 270 PetscFVCellGeom fvgeom; 271 PetscInt dim, dimEmbed; 272 PetscErrorCode ierr; 273 274 PetscFunctionBeginHot; 275 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 276 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 277 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 278 switch (type) { 279 case DM_BC_ESSENTIAL: 280 case DM_BC_NATURAL: 281 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; 282 case DM_BC_ESSENTIAL_FIELD: 283 case DM_BC_NATURAL_FIELD: 284 ierr = DMProjectPoint_Field_Private(dm, prob, dmAux, probAux, time, localU, localA, fegeom, sp, p, Ncc, comps, 285 basisTab, basisDerTab, basisTabAux, basisDerTabAux, 286 (void (**)(PetscInt, PetscInt, PetscInt, 287 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 288 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 289 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 290 case DM_BC_ESSENTIAL_BD_FIELD: 291 ierr = DMProjectPoint_BdField_Private(dm, prob, dmAux, probAux, time, localU, localA, fegeom, sp, p, Ncc, comps, 292 basisTab, basisDerTab, basisTabAux, basisDerTabAux, 293 (void (**)(PetscInt, PetscInt, PetscInt, 294 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 295 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 296 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 297 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type); 298 } 299 PetscFunctionReturn(0); 300 } 301 302 static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints) 303 { 304 PetscReal *points; 305 PetscInt f, numPoints; 306 PetscErrorCode ierr; 307 308 PetscFunctionBegin; 309 numPoints = 0; 310 for (f = 0; f < Nf; ++f) { 311 if (funcs[f]) { 312 PetscQuadrature fAllPoints; 313 PetscInt fNumPoints; 314 315 ierr = PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);CHKERRQ(ierr); 316 ierr = PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);CHKERRQ(ierr); 317 numPoints += fNumPoints; 318 } 319 } 320 ierr = PetscMalloc1(dim*numPoints,&points);CHKERRQ(ierr); 321 numPoints = 0; 322 for (f = 0; f < Nf; ++f) { 323 PetscInt spDim; 324 325 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 326 if (funcs[f]) { 327 PetscQuadrature fAllPoints; 328 PetscInt qdim, fNumPoints, q; 329 const PetscReal *fPoints; 330 331 ierr = PetscDualSpaceGetAllPoints(sp[f],&fAllPoints);CHKERRQ(ierr); 332 ierr = PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);CHKERRQ(ierr); 333 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); 334 for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q]; 335 numPoints += fNumPoints; 336 } 337 } 338 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);CHKERRQ(ierr); 339 ierr = PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);CHKERRQ(ierr); 340 PetscFunctionReturn(0); 341 } 342 343 static PetscErrorCode DMGetFirstLabelEntry_Private(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *lStart, PetscDS *prob) 344 { 345 DMLabel depthLabel; 346 PetscInt dim, cdepth, ls = -1, i; 347 PetscErrorCode ierr; 348 349 PetscFunctionBegin; 350 if (lStart) *lStart = -1; 351 if (!label) PetscFunctionReturn(0); 352 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 353 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 354 cdepth = dim - height; 355 for (i = 0; i < numIds; ++i) { 356 IS pointIS; 357 const PetscInt *points; 358 PetscInt pdepth; 359 360 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 361 if (!pointIS) continue; /* No points with that id on this process */ 362 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 363 ierr = DMLabelGetValue(depthLabel, points[0], &pdepth);CHKERRQ(ierr); 364 if (pdepth == cdepth) { 365 ls = points[0]; 366 if (prob) {ierr = DMGetCellDS(dm, ls, prob);CHKERRQ(ierr);} 367 } 368 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 369 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 370 if (ls >= 0) break; 371 } 372 if (lStart) *lStart = ls; 373 PetscFunctionReturn(0); 374 } 375 376 /* 377 This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM 378 379 There are several different scenarios: 380 381 1) Volumetric mesh with volumetric auxiliary data 382 383 Here minHeight=0 since we loop over cells. 384 385 2) Boundary mesh with boundary auxiliary data 386 387 Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation. 388 389 3) Volumetric mesh with boundary auxiliary data 390 391 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. 392 393 The maxHeight is used to support enforcement of constraints in DMForest. 394 395 If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it. 396 397 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. 398 - 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 399 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 400 dual spaces for the boundary from our input spaces. 401 - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them. 402 403 We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration 404 405 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. 406 */ 407 static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, 408 PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[], 409 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, 410 InsertMode mode, Vec localX) 411 { 412 DM dmAux = NULL, tdm; 413 PetscDS prob = NULL, probAux = NULL; 414 Vec localA = NULL, tv; 415 PetscSection section; 416 PetscDualSpace *sp, *cellsp; 417 PetscReal **basisTab = NULL, **basisDerTab = NULL, **basisTabAux = NULL, **basisDerTabAux = NULL; 418 PetscInt *Nc; 419 PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, Nf, NfAux = 0, f; 420 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE, transform; 421 DMField coordField; 422 DMLabel depthLabel; 423 PetscQuadrature allPoints = NULL; 424 PetscErrorCode ierr; 425 426 PetscFunctionBegin; 427 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 428 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);CHKERRQ(ierr); 429 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 430 ierr = DMPlexGetVTKCellHeight(dm, &minHeight);CHKERRQ(ierr); 431 ierr = DMGetBasisTransformDM_Internal(dm, &tdm);CHKERRQ(ierr); 432 ierr = DMGetBasisTransformVec_Internal(dm, &tv);CHKERRQ(ierr); 433 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 434 /* Auxiliary information can only be used with interpolation of field functions */ 435 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 436 if (dmAux && !localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector"); 437 if (!minHeight && dmAux) { 438 DMLabel spmap; 439 440 /* If dmAux is a surface, then force the projection to take place over a surface */ 441 ierr = DMPlexGetSubpointMap(dmAux, &spmap);CHKERRQ(ierr); 442 if (spmap) { 443 ierr = DMPlexGetVTKCellHeight(dmAux, &minHeight);CHKERRQ(ierr); 444 auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE; 445 } 446 } 447 } 448 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 449 ierr = DMPlexGetDepthLabel(dm,&depthLabel);CHKERRQ(ierr); 450 ierr = DMPlexGetMaxProjectionHeight(dm, &maxHeight);CHKERRQ(ierr); 451 maxHeight = PetscMax(maxHeight, minHeight); 452 if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim); 453 ierr = DMGetFirstLabelEntry_Private(dm, label, numIds, ids, 0, NULL, &prob);CHKERRQ(ierr); 454 if (!prob) {ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);} 455 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 456 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 457 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 458 if (dmAux) { 459 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 460 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 461 } 462 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 463 ierr = PetscMalloc2(Nf, &isFE, Nf, &sp);CHKERRQ(ierr); 464 if (maxHeight > 0) {ierr = PetscMalloc1(Nf, &cellsp);CHKERRQ(ierr);} 465 else {cellsp = sp;} 466 if (localU && localU != localX) {ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);CHKERRQ(ierr);} 467 /* Get cell dual spaces */ 468 for (f = 0; f < Nf; ++f) { 469 PetscObject obj; 470 PetscClassId id; 471 472 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 473 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 474 if (id == PETSCFE_CLASSID) { 475 PetscFE fe = (PetscFE) obj; 476 477 hasFE = PETSC_TRUE; 478 isFE[f] = PETSC_TRUE; 479 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 480 } else if (id == PETSCFV_CLASSID) { 481 PetscFV fv = (PetscFV) obj; 482 483 hasFV = PETSC_TRUE; 484 isFE[f] = PETSC_FALSE; 485 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 486 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 487 } 488 ierr = DMGetCoordinateField(dm,&coordField);CHKERRQ(ierr); 489 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 490 PetscInt effectiveHeight = auxBd ? minHeight : 0; 491 PetscFE fem, subfem; 492 const PetscReal *points; 493 PetscInt numPoints; 494 495 if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation"); 496 for (f = 0; f < Nf; ++f) { 497 if (!effectiveHeight) {sp[f] = cellsp[f];} 498 else {ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr);} 499 } 500 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);CHKERRQ(ierr); 501 ierr = PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);CHKERRQ(ierr); 502 ierr = PetscMalloc4(Nf, &basisTab, Nf, &basisDerTab, NfAux, &basisTabAux, NfAux, &basisDerTabAux);CHKERRQ(ierr); 503 for (f = 0; f < Nf; ++f) { 504 if (!isFE[f]) continue; 505 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 506 if (!effectiveHeight) {subfem = fem;} 507 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 508 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 509 } 510 for (f = 0; f < NfAux; ++f) { 511 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 512 if (!effectiveHeight || auxBd) {subfem = fem;} 513 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 514 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 515 } 516 } 517 /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */ 518 for (h = minHeight; h <= maxHeight; h++) { 519 PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight); 520 PetscDS probEff = prob; 521 PetscScalar *values; 522 PetscBool *fieldActive; 523 PetscInt maxDegree; 524 PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues; 525 IS heightIS; 526 527 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 528 ierr = DMGetFirstLabelEntry_Private(dm, label, numIds, ids, h, &lStart, NULL);CHKERRQ(ierr); 529 ierr = DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);CHKERRQ(ierr); 530 if (!h) { 531 PetscInt cEndInterior; 532 533 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 534 pEnd = cEndInterior < 0 ? pEnd : cEndInterior; 535 } 536 if (pEnd <= pStart) { 537 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 538 continue; 539 } 540 /* Compute totDim, the number of dofs in the closure of a point at this height */ 541 totDim = 0; 542 for (f = 0; f < Nf; ++f) { 543 if (!effectiveHeight) { 544 sp[f] = cellsp[f]; 545 } else { 546 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr); 547 if (!sp[f]) continue; 548 } 549 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 550 totDim += spDim; 551 } 552 ierr = DMPlexVecGetClosure(dm, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);CHKERRQ(ierr); 553 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim); 554 if (!totDim) { 555 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 556 continue; 557 } 558 if (effectiveHeight) {ierr = PetscDSGetHeightSubspace(prob, effectiveHeight, &probEff);CHKERRQ(ierr);} 559 /* Loop over points at this height */ 560 ierr = DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 561 ierr = DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 562 for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 563 if (label) { 564 PetscInt i; 565 566 for (i = 0; i < numIds; ++i) { 567 IS pointIS, isectIS; 568 const PetscInt *points; 569 PetscInt n; 570 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 571 PetscQuadrature quad = NULL; 572 573 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 574 if (!pointIS) continue; /* No points with that id on this process */ 575 ierr = ISIntersect(pointIS,heightIS,&isectIS);CHKERRQ(ierr); 576 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 577 if (!isectIS) continue; 578 ierr = ISGetLocalSize(isectIS, &n);CHKERRQ(ierr); 579 ierr = ISGetIndices(isectIS, &points);CHKERRQ(ierr); 580 ierr = DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);CHKERRQ(ierr); 581 if (maxDegree <= 1) { 582 ierr = DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);CHKERRQ(ierr); 583 } 584 if (!quad) { 585 if (!h && allPoints) { 586 quad = allPoints; 587 allPoints = NULL; 588 } else { 589 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);CHKERRQ(ierr); 590 } 591 } 592 ierr = DMFieldCreateFEGeom(coordField,isectIS,quad,PETSC_FALSE,&fegeom);CHKERRQ(ierr); 593 for (p = 0; p < n; ++p) { 594 const PetscInt point = points[p]; 595 596 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 597 ierr = PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 598 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); 599 if (ierr) { 600 PetscErrorCode ierr2; 601 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 602 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 603 CHKERRQ(ierr); 604 } 605 if (transform) {ierr = DMPlexBasisTransformPoint_Internal(dm, tdm, tv, point, fieldActive, PETSC_FALSE, values);CHKERRQ(ierr);} 606 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, Ncc, comps, values, mode);CHKERRQ(ierr); 607 } 608 ierr = PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 609 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 610 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 611 ierr = ISRestoreIndices(isectIS, &points);CHKERRQ(ierr); 612 ierr = ISDestroy(&isectIS);CHKERRQ(ierr); 613 } 614 } else { 615 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 616 PetscQuadrature quad = NULL; 617 IS pointIS; 618 619 ierr = ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);CHKERRQ(ierr); 620 ierr = DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);CHKERRQ(ierr); 621 if (maxDegree <= 1) { 622 ierr = DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);CHKERRQ(ierr); 623 } 624 if (!quad) { 625 if (!h && allPoints) { 626 quad = allPoints; 627 allPoints = NULL; 628 } else { 629 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);CHKERRQ(ierr); 630 } 631 } 632 ierr = DMFieldCreateFEGeom(coordField,pointIS,quad,PETSC_FALSE,&fegeom);CHKERRQ(ierr); 633 for (p = pStart; p < pEnd; ++p) { 634 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 635 ierr = PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);CHKERRQ(ierr); 636 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); 637 if (ierr) { 638 PetscErrorCode ierr2; 639 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 640 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 641 CHKERRQ(ierr); 642 } 643 if (transform) {ierr = DMPlexBasisTransformPoint_Internal(dm, tdm, tv, p, fieldActive, PETSC_FALSE, values);CHKERRQ(ierr);} 644 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, Ncc, comps, values, mode);CHKERRQ(ierr); 645 } 646 ierr = PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);CHKERRQ(ierr); 647 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 648 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 649 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 650 } 651 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 652 ierr = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 653 ierr = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 654 } 655 /* Cleanup */ 656 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 657 PetscInt effectiveHeight = auxBd ? minHeight : 0; 658 PetscFE fem, subfem; 659 660 for (f = 0; f < Nf; ++f) { 661 if (!isFE[f]) continue; 662 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 663 if (!effectiveHeight) {subfem = fem;} 664 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 665 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 666 } 667 for (f = 0; f < NfAux; ++f) { 668 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 669 if (!effectiveHeight || auxBd) {subfem = fem;} 670 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 671 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 672 } 673 ierr = PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);CHKERRQ(ierr); 674 } 675 ierr = PetscQuadratureDestroy(&allPoints);CHKERRQ(ierr); 676 ierr = PetscFree2(isFE, sp);CHKERRQ(ierr); 677 if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 678 PetscFunctionReturn(0); 679 } 680 681 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, 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, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 687 PetscFunctionReturn(0); 688 } 689 690 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) 691 { 692 PetscErrorCode ierr; 693 694 PetscFunctionBegin; 695 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 696 PetscFunctionReturn(0); 697 } 698 699 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 700 void (**funcs)(PetscInt, PetscInt, PetscInt, 701 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 702 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 703 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 704 InsertMode mode, Vec localX) 705 { 706 PetscErrorCode ierr; 707 708 PetscFunctionBegin; 709 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 710 PetscFunctionReturn(0); 711 } 712 713 PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 714 void (**funcs)(PetscInt, PetscInt, PetscInt, 715 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 716 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 717 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 718 InsertMode mode, Vec localX) 719 { 720 PetscErrorCode ierr; 721 722 PetscFunctionBegin; 723 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 724 PetscFunctionReturn(0); 725 } 726