1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #include <petsc/private/petscfeimpl.h> 4 5 /* 6 DMProjectPoint_Func_Private - Interpolate the given function in the output basis on the given point 7 8 Input Parameters: 9 + dm - The output DM 10 . ds - The output DS 11 . dmIn - The input DM 12 . dsIn - The input DS 13 . time - The time for this evaluation 14 . fegeom - The FE geometry for this point 15 . fvgeom - The FV geometry for this point 16 . isFE - Flag indicating whether each output field has an FE discretization 17 . sp - The output PetscDualSpace for each field 18 . funcs - The evaluation function for each field 19 - ctxs - The user context for each field 20 21 Output Parameter: 22 . values - The value for each dual basis vector in the output dual space 23 24 Level: developer 25 26 .seealso: DMProjectPoint_Field_Private() 27 */ 28 static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, PetscReal time, PetscFEGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[], 29 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, 30 PetscScalar values[]) 31 { 32 PetscInt coordDim, Nf, *Nc, f, spDim, d, v, tp; 33 PetscBool isAffine, isHybrid, transform; 34 PetscErrorCode ierr; 35 36 PetscFunctionBeginHot; 37 ierr = DMGetCoordinateDim(dmIn, &coordDim);CHKERRQ(ierr); 38 ierr = DMHasBasisTransform(dmIn, &transform);CHKERRQ(ierr); 39 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 40 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 41 ierr = PetscDSGetHybrid(ds, &isHybrid);CHKERRQ(ierr); 42 /* Get values for closure */ 43 isAffine = fegeom->isAffine; 44 for (f = 0, v = 0, tp = 0; f < Nf; ++f) { 45 void * const ctx = ctxs ? ctxs[f] : NULL; 46 47 if (!sp[f]) continue; 48 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 49 if (funcs[f]) { 50 if (isFE[f]) { 51 PetscQuadrature allPoints; 52 PetscInt q, dim, numPoints; 53 const PetscReal *points; 54 PetscScalar *pointEval; 55 PetscReal *x; 56 DM rdm; 57 58 ierr = PetscDualSpaceGetDM(sp[f],&rdm);CHKERRQ(ierr); 59 ierr = PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);CHKERRQ(ierr); 60 ierr = PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);CHKERRQ(ierr); 61 ierr = DMGetWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 62 ierr = DMGetWorkArray(rdm,coordDim,MPIU_REAL,&x);CHKERRQ(ierr); 63 for (q = 0; q < numPoints; q++, tp++) { 64 const PetscReal *v0; 65 66 if (isAffine) { 67 const PetscReal *refpoint = &points[q*dim]; 68 PetscReal injpoint[3] = {0., 0., 0.}; 69 70 if (dim != fegeom->dim) { 71 if (isHybrid) { 72 /* We just need to inject into the higher dimensional space assuming the last dimension is collapsed */ 73 for (d = 0; d < dim; ++d) injpoint[d] = refpoint[d]; 74 refpoint = injpoint; 75 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reference spatial dimension %D != %D dual basis spatial dimension", fegeom->dim, dim); 76 } 77 CoordinatesRefToReal(coordDim, fegeom->dim, fegeom->xi, fegeom->v, fegeom->J, refpoint, x); 78 v0 = x; 79 } else { 80 v0 = &fegeom->v[tp*coordDim]; 81 } 82 if (transform) {ierr = DMPlexBasisTransformApplyReal_Internal(dmIn, v0, PETSC_TRUE, coordDim, v0, x, dm->transformCtx);CHKERRQ(ierr); v0 = x;} 83 ierr = (*funcs[f])(coordDim, time, v0, Nc[f], &pointEval[Nc[f]*q], ctx);CHKERRQ(ierr); 84 } 85 /* Transform point evaluations pointEval[q,c] */ 86 ierr = PetscDualSpacePullback(sp[f], fegeom, numPoints, Nc[f], pointEval);CHKERRQ(ierr); 87 ierr = PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);CHKERRQ(ierr); 88 ierr = DMRestoreWorkArray(rdm,coordDim,MPIU_REAL,&x);CHKERRQ(ierr); 89 ierr = DMRestoreWorkArray(rdm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 90 v += spDim; 91 if (isHybrid && (f < Nf-1)) { 92 for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim]; 93 } 94 } else { 95 for (d = 0; d < spDim; ++d, ++v) { 96 ierr = PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr); 97 } 98 } 99 } else { 100 for (d = 0; d < spDim; d++, v++) values[v] = 0.; 101 if (isHybrid && (f < Nf-1)) { 102 for (d = 0; d < spDim; d++, v++) values[v] = 0.; 103 } 104 } 105 } 106 PetscFunctionReturn(0); 107 } 108 109 /* 110 DMProjectPoint_Field_Private - Interpolate a function of the given field, in the input basis, using the output basis on the given point 111 112 Input Parameters: 113 + dm - The output DM 114 . ds - The output DS 115 . dmIn - The input DM 116 . dsIn - The input DS 117 . dmAux - The auxiliary DM, which is always for the input space 118 . dsAux - The auxiliary DS, which is always for the input space 119 . time - The time for this evaluation 120 . localU - The local solution 121 . localA - The local auziliary fields 122 . cgeom - The FE geometry for this point 123 . sp - The output PetscDualSpace for each field 124 . p - The point in the output DM 125 . T - Input basis and derviatives for each field tabulated on the quadrature points 126 . TAux - Auxiliary basis and derivatives for each aux field tabulated on the quadrature points 127 . funcs - The evaluation function for each field 128 - ctxs - The user context for each field 129 130 Output Parameter: 131 . values - The value for each dual basis vector in the output dual space 132 133 Note: Not supported for FV 134 135 Level: developer 136 137 .seealso: DMProjectPoint_Field_Private() 138 */ 139 static PetscErrorCode DMProjectPoint_Field_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *cgeom, PetscDualSpace sp[], PetscInt p, 140 PetscTabulation *T, PetscTabulation *TAux, 141 void (**funcs)(PetscInt, PetscInt, PetscInt, 142 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 143 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 144 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, 145 PetscScalar values[]) 146 { 147 PetscSection section, sectionAux = NULL; 148 PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc; 149 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 150 PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL; 151 const PetscScalar *constants; 152 PetscReal *x; 153 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc; 154 PetscFEGeom fegeom; 155 const PetscInt dE = cgeom->dimEmbed; 156 PetscInt numConstants, Nf, NfIn, NfAux = 0, f, spDim, d, v, inp, tp = 0; 157 PetscBool isAffine, isHybrid, transform; 158 PetscErrorCode ierr; 159 160 PetscFunctionBeginHot; 161 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 162 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 163 ierr = PetscDSGetHybrid(ds, &isHybrid);CHKERRQ(ierr); 164 ierr = PetscDSGetNumFields(dsIn, &NfIn);CHKERRQ(ierr); 165 ierr = PetscDSGetComponentOffsets(dsIn, &uOff);CHKERRQ(ierr); 166 ierr = PetscDSGetComponentDerivativeOffsets(dsIn, &uOff_x);CHKERRQ(ierr); 167 ierr = PetscDSGetEvaluationArrays(dsIn, &u, &bc /*&u_t*/, &u_x);CHKERRQ(ierr); 168 ierr = PetscDSGetWorkspace(dsIn, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 169 ierr = PetscDSGetConstants(dsIn, &numConstants, &constants);CHKERRQ(ierr); 170 ierr = DMHasBasisTransform(dmIn, &transform);CHKERRQ(ierr); 171 ierr = DMGetLocalSection(dmIn, §ion);CHKERRQ(ierr); 172 ierr = DMGetEnclosurePoint(dmIn, dm, encIn, p, &inp);CHKERRQ(ierr); 173 ierr = DMPlexVecGetClosure(dmIn, section, localU, inp, NULL, &coefficients);CHKERRQ(ierr); 174 if (dmAux) { 175 PetscInt subp; 176 177 ierr = DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);CHKERRQ(ierr); 178 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 179 ierr = DMGetLocalSection(dmAux, §ionAux);CHKERRQ(ierr); 180 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 181 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 182 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);CHKERRQ(ierr); 183 ierr = DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);CHKERRQ(ierr); 184 } 185 /* Get values for closure */ 186 isAffine = cgeom->isAffine; 187 fegeom.dim = cgeom->dim; 188 fegeom.dimEmbed = cgeom->dimEmbed; 189 if (isAffine) { 190 fegeom.v = x; 191 fegeom.xi = cgeom->xi; 192 fegeom.J = cgeom->J; 193 fegeom.invJ = cgeom->invJ; 194 fegeom.detJ = cgeom->detJ; 195 } 196 for (f = 0, v = 0; f < Nf; ++f) { 197 PetscQuadrature allPoints; 198 PetscInt q, dim, numPoints; 199 const PetscReal *points; 200 PetscScalar *pointEval; 201 DM dm; 202 203 if (!sp[f]) continue; 204 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 205 if (!funcs[f]) { 206 for (d = 0; d < spDim; d++, v++) values[v] = 0.; 207 if (isHybrid && (f < Nf-1)) { 208 for (d = 0; d < spDim; d++, v++) values[v] = 0.; 209 } 210 continue; 211 } 212 ierr = PetscDualSpaceGetDM(sp[f],&dm);CHKERRQ(ierr); 213 ierr = PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);CHKERRQ(ierr); 214 ierr = PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);CHKERRQ(ierr); 215 ierr = DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 216 for (q = 0; q < numPoints; ++q, ++tp) { 217 if (isAffine) { 218 CoordinatesRefToReal(dE, cgeom->dim, fegeom.xi, cgeom->v, fegeom.J, &points[q*dim], x); 219 } else { 220 fegeom.v = &cgeom->v[tp*dE]; 221 fegeom.J = &cgeom->J[tp*dE*dE]; 222 fegeom.invJ = &cgeom->invJ[tp*dE*dE]; 223 fegeom.detJ = &cgeom->detJ[tp]; 224 } 225 ierr = PetscFEEvaluateFieldJets_Internal(dsIn, NfIn, 0, tp, T, &fegeom, coefficients, coefficients_t, u, u_x, u_t);CHKERRQ(ierr); 226 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &fegeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);CHKERRQ(ierr);} 227 if (transform) {ierr = DMPlexBasisTransformApplyReal_Internal(dmIn, fegeom.v, PETSC_TRUE, dE, fegeom.v, fegeom.v, dm->transformCtx);CHKERRQ(ierr);} 228 (*funcs[f])(dE, NfIn, 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]); 229 } 230 ierr = PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);CHKERRQ(ierr); 231 ierr = DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 232 v += spDim; 233 /* TODO: For now, set both sides equal, but this should use info from other support cell */ 234 if (isHybrid && (f < Nf-1)) { 235 for (d = 0; d < spDim; d++, v++) values[v] = values[v - spDim]; 236 } 237 } 238 ierr = DMPlexVecRestoreClosure(dmIn, section, localU, inp, NULL, &coefficients);CHKERRQ(ierr); 239 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);CHKERRQ(ierr);} 240 PetscFunctionReturn(0); 241 } 242 243 static PetscErrorCode DMProjectPoint_BdField_Private(DM dm, PetscDS ds, DM dmIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscReal time, Vec localU, Vec localA, PetscFEGeom *fgeom, PetscDualSpace sp[], PetscInt p, 244 PetscTabulation *T, PetscTabulation *TAux, 245 void (**funcs)(PetscInt, PetscInt, PetscInt, 246 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 247 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 248 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, 249 PetscScalar values[]) 250 { 251 PetscSection section, sectionAux = NULL; 252 PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *bc; 253 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 254 PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL; 255 const PetscScalar *constants; 256 PetscReal *x; 257 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nc; 258 PetscFEGeom fegeom, cgeom; 259 const PetscInt dE = fgeom->dimEmbed; 260 PetscInt numConstants, Nf, NfAux = 0, f, spDim, d, v, tp = 0; 261 PetscBool isAffine; 262 PetscErrorCode ierr; 263 264 PetscFunctionBeginHot; 265 if (dm != dmIn) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not yet upgraded to use different input DM"); 266 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 267 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 268 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 269 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 270 ierr = PetscDSGetEvaluationArrays(ds, &u, &bc /*&u_t*/, &u_x);CHKERRQ(ierr); 271 ierr = PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 272 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 273 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 274 ierr = DMPlexVecGetClosure(dmIn, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 275 if (dmAux) { 276 PetscInt subp; 277 278 ierr = DMGetEnclosurePoint(dmAux, dm, encAux, p, &subp);CHKERRQ(ierr); 279 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 280 ierr = DMGetLocalSection(dmAux, §ionAux);CHKERRQ(ierr); 281 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 282 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 283 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL /*&a_t*/, &a_x);CHKERRQ(ierr); 284 ierr = DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);CHKERRQ(ierr); 285 } 286 /* Get values for closure */ 287 isAffine = fgeom->isAffine; 288 fegeom.n = 0; 289 fegeom.J = 0; 290 fegeom.v = 0; 291 fegeom.xi = 0; 292 cgeom.dim = fgeom->dim; 293 cgeom.dimEmbed = fgeom->dimEmbed; 294 if (isAffine) { 295 fegeom.v = x; 296 fegeom.xi = fgeom->xi; 297 fegeom.J = fgeom->J; 298 fegeom.invJ = fgeom->invJ; 299 fegeom.detJ = fgeom->detJ; 300 fegeom.n = fgeom->n; 301 302 cgeom.J = fgeom->suppJ[0]; 303 cgeom.invJ = fgeom->suppInvJ[0]; 304 cgeom.detJ = fgeom->suppDetJ[0]; 305 } 306 for (f = 0, v = 0; f < Nf; ++f) { 307 PetscQuadrature allPoints; 308 PetscInt q, dim, numPoints; 309 const PetscReal *points; 310 PetscScalar *pointEval; 311 DM dm; 312 313 if (!sp[f]) continue; 314 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 315 if (!funcs[f]) { 316 for (d = 0; d < spDim; d++, v++) values[v] = 0.; 317 continue; 318 } 319 ierr = PetscDualSpaceGetDM(sp[f],&dm);CHKERRQ(ierr); 320 ierr = PetscDualSpaceGetAllData(sp[f], &allPoints, NULL);CHKERRQ(ierr); 321 ierr = PetscQuadratureGetData(allPoints,&dim,NULL,&numPoints,&points,NULL);CHKERRQ(ierr); 322 ierr = DMGetWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 323 for (q = 0; q < numPoints; ++q, ++tp) { 324 if (isAffine) { 325 CoordinatesRefToReal(dE, fgeom->dim, fegeom.xi, fgeom->v, fegeom.J, &points[q*dim], x); 326 } else { 327 fegeom.v = &fgeom->v[tp*dE]; 328 fegeom.J = &fgeom->J[tp*dE*dE]; 329 fegeom.invJ = &fgeom->invJ[tp*dE*dE]; 330 fegeom.detJ = &fgeom->detJ[tp]; 331 fegeom.n = &fgeom->n[tp*dE]; 332 333 cgeom.J = &fgeom->suppJ[0][tp*dE*dE]; 334 cgeom.invJ = &fgeom->suppInvJ[0][tp*dE*dE]; 335 cgeom.detJ = &fgeom->suppDetJ[0][tp]; 336 } 337 /* TODO We should use cgeom here, instead of fegeom, however the geometry coming in through fgeom does not have the support cell geometry */ 338 ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, tp, T, &cgeom, coefficients, coefficients_t, u, u_x, u_t);CHKERRQ(ierr); 339 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, tp, TAux, &cgeom, coefficientsAux, coefficientsAux_t, a, a_x, a_t);CHKERRQ(ierr);} 340 (*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]); 341 } 342 ierr = PetscDualSpaceApplyAll(sp[f], pointEval, &values[v]);CHKERRQ(ierr); 343 ierr = DMRestoreWorkArray(dm,numPoints*Nc[f],MPIU_SCALAR,&pointEval);CHKERRQ(ierr); 344 v += spDim; 345 } 346 ierr = DMPlexVecRestoreClosure(dmIn, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 347 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);CHKERRQ(ierr);} 348 PetscFunctionReturn(0); 349 } 350 351 static PetscErrorCode DMProjectPoint_Private(DM dm, PetscDS ds, DM dmIn, DMEnclosureType encIn, PetscDS dsIn, DM dmAux, DMEnclosureType encAux, PetscDS dsAux, PetscFEGeom *fegeom, PetscInt effectiveHeight, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[], 352 PetscDualSpace sp[], PetscInt p, PetscTabulation *T, PetscTabulation *TAux, 353 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[]) 354 { 355 PetscFVCellGeom fvgeom; 356 PetscInt dim, dimEmbed; 357 PetscErrorCode ierr; 358 359 PetscFunctionBeginHot; 360 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 361 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 362 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 363 switch (type) { 364 case DM_BC_ESSENTIAL: 365 case DM_BC_NATURAL: 366 ierr = DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);CHKERRQ(ierr);break; 367 case DM_BC_ESSENTIAL_FIELD: 368 case DM_BC_NATURAL_FIELD: 369 ierr = DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux, 370 (void (**)(PetscInt, PetscInt, PetscInt, 371 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 372 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 373 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 374 case DM_BC_ESSENTIAL_BD_FIELD: 375 ierr = DMProjectPoint_BdField_Private(dm, ds, dmIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, T, TAux, 376 (void (**)(PetscInt, PetscInt, PetscInt, 377 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 378 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 379 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 380 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type); 381 } 382 PetscFunctionReturn(0); 383 } 384 385 static PetscErrorCode PetscDualSpaceGetAllPointsUnion(PetscInt Nf, PetscDualSpace *sp, PetscInt dim, void (**funcs)(void), PetscQuadrature *allPoints) 386 { 387 PetscReal *points; 388 PetscInt f, numPoints; 389 PetscErrorCode ierr; 390 391 PetscFunctionBegin; 392 numPoints = 0; 393 for (f = 0; f < Nf; ++f) { 394 if (funcs[f]) { 395 PetscQuadrature fAllPoints; 396 PetscInt fNumPoints; 397 398 ierr = PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);CHKERRQ(ierr); 399 ierr = PetscQuadratureGetData(fAllPoints, NULL, NULL, &fNumPoints, NULL, NULL);CHKERRQ(ierr); 400 numPoints += fNumPoints; 401 } 402 } 403 ierr = PetscMalloc1(dim*numPoints,&points);CHKERRQ(ierr); 404 numPoints = 0; 405 for (f = 0; f < Nf; ++f) { 406 if (funcs[f]) { 407 PetscQuadrature fAllPoints; 408 PetscInt qdim, fNumPoints, q; 409 const PetscReal *fPoints; 410 411 ierr = PetscDualSpaceGetAllData(sp[f],&fAllPoints, NULL);CHKERRQ(ierr); 412 ierr = PetscQuadratureGetData(fAllPoints, &qdim, NULL, &fNumPoints, &fPoints, NULL);CHKERRQ(ierr); 413 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); 414 for (q = 0; q < fNumPoints*dim; ++q) points[numPoints*dim+q] = fPoints[q]; 415 numPoints += fNumPoints; 416 } 417 } 418 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);CHKERRQ(ierr); 419 ierr = PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);CHKERRQ(ierr); 420 PetscFunctionReturn(0); 421 } 422 423 static PetscErrorCode DMGetFirstLabelEntry_Private(DM dm, DM odm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt height, PetscInt *lStart, PetscDS *ds) 424 { 425 DM plex; 426 DMEnclosureType enc; 427 DMLabel depthLabel; 428 PetscInt dim, cdepth, ls = -1, i; 429 PetscErrorCode ierr; 430 431 PetscFunctionBegin; 432 if (lStart) *lStart = -1; 433 if (!label) PetscFunctionReturn(0); 434 ierr = DMGetEnclosureRelation(dm, odm, &enc);CHKERRQ(ierr); 435 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 436 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 437 ierr = DMPlexGetDepthLabel(plex, &depthLabel);CHKERRQ(ierr); 438 cdepth = dim - height; 439 for (i = 0; i < numIds; ++i) { 440 IS pointIS; 441 const PetscInt *points; 442 PetscInt pdepth, point; 443 444 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 445 if (!pointIS) continue; /* No points with that id on this process */ 446 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 447 ierr = DMGetEnclosurePoint(dm, odm, enc, points[0], &point);CHKERRQ(ierr); 448 ierr = DMLabelGetValue(depthLabel, point, &pdepth);CHKERRQ(ierr); 449 if (pdepth == cdepth) { 450 ls = point; 451 if (ds) {ierr = DMGetCellDS(dm, ls, ds);CHKERRQ(ierr);} 452 } 453 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 454 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 455 if (ls >= 0) break; 456 } 457 if (lStart) *lStart = ls; 458 ierr = DMDestroy(&plex);CHKERRQ(ierr); 459 PetscFunctionReturn(0); 460 } 461 462 /* 463 This function iterates over a manifold, and interpolates the input function/field using the basis provided by the DS in our DM 464 465 There are several different scenarios: 466 467 1) Volumetric mesh with volumetric auxiliary data 468 469 Here minHeight=0 since we loop over cells. 470 471 2) Boundary mesh with boundary auxiliary data 472 473 Here minHeight=1 since we loop over faces. This normally happens since we hang cells off of our boundary meshes to facilitate computation. 474 475 3) Volumetric mesh with boundary auxiliary data 476 477 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. 478 479 The maxHeight is used to support enforcement of constraints in DMForest. 480 481 If localU is given and not equal to localX, we call DMPlexInsertBoundaryValues() to complete it. 482 483 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. 484 - 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 485 is for the volumetric mesh, but we are iterating over a surface, then the effective height is nonzero. When the effective height is nonzero, we need to extract 486 dual spaces for the boundary from our input spaces. 487 - After extracting all quadrature points, we tabulate the input fields and auxiliary fields on them. 488 489 We check that the #dof(closure(p)) == #dual basis functionals(p) for a representative p in the iteration 490 491 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. 492 */ 493 static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, 494 PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[], 495 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, 496 InsertMode mode, Vec localX) 497 { 498 DM plex, dmIn, plexIn, dmAux = NULL, plexAux = NULL, tdm; 499 DMEnclosureType encIn, encAux; 500 PetscDS ds = NULL, dsIn = NULL, dsAux = NULL; 501 Vec localA = NULL, tv; 502 PetscSection section; 503 PetscDualSpace *sp, *cellsp; 504 PetscTabulation *T = NULL, *TAux = NULL; 505 PetscInt *Nc; 506 PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, Nf, NfIn, NfAux = 0, f; 507 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE, isHybrid = PETSC_FALSE, transform; 508 DMField coordField; 509 DMLabel depthLabel; 510 PetscQuadrature allPoints = NULL; 511 PetscErrorCode ierr; 512 513 PetscFunctionBegin; 514 if (localU) {ierr = VecGetDM(localU, &dmIn);CHKERRQ(ierr);} 515 else {dmIn = dm;} 516 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 517 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);CHKERRQ(ierr); 518 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 519 ierr = DMConvert(dmIn, DMPLEX, &plexIn);CHKERRQ(ierr); 520 ierr = DMGetEnclosureRelation(dmIn, dm, &encIn);CHKERRQ(ierr); 521 ierr = DMGetEnclosureRelation(dmAux, dm, &encAux);CHKERRQ(ierr); 522 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 523 ierr = DMPlexGetVTKCellHeight(plex, &minHeight);CHKERRQ(ierr); 524 ierr = DMGetBasisTransformDM_Internal(dm, &tdm);CHKERRQ(ierr); 525 ierr = DMGetBasisTransformVec_Internal(dm, &tv);CHKERRQ(ierr); 526 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 527 /* Auxiliary information can only be used with interpolation of field functions */ 528 if (dmAux) { 529 ierr = DMConvert(dmAux, DMPLEX, &plexAux);CHKERRQ(ierr); 530 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) { 531 if (!localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector"); 532 if (!minHeight) { 533 DMLabel spmap; 534 535 /* If dmAux is a surface, then force the projection to take place over a surface */ 536 ierr = DMPlexGetSubpointMap(plexAux, &spmap);CHKERRQ(ierr); 537 if (spmap) { 538 ierr = DMPlexGetVTKCellHeight(plexAux, &minHeight);CHKERRQ(ierr); 539 auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE; 540 } 541 } 542 } 543 } 544 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr); 545 ierr = DMPlexGetDepthLabel(plex, &depthLabel);CHKERRQ(ierr); 546 ierr = DMPlexGetMaxProjectionHeight(plex, &maxHeight);CHKERRQ(ierr); 547 maxHeight = PetscMax(maxHeight, minHeight); 548 if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim); 549 ierr = DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, 0, NULL, &ds);CHKERRQ(ierr); 550 if (!ds) {ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);} 551 ierr = DMGetFirstLabelEntry_Private(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);CHKERRQ(ierr); 552 if (!dsIn) {ierr = DMGetDS(dmIn, &dsIn);CHKERRQ(ierr);} 553 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 554 ierr = PetscDSGetNumFields(dsIn, &NfIn);CHKERRQ(ierr); 555 ierr = PetscDSGetHybrid(ds, &isHybrid);CHKERRQ(ierr); 556 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 557 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 558 if (dmAux) { 559 ierr = DMGetDS(dmAux, &dsAux);CHKERRQ(ierr); 560 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 561 } 562 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 563 ierr = PetscMalloc2(Nf, &isFE, Nf, &sp);CHKERRQ(ierr); 564 if (maxHeight > 0) {ierr = PetscMalloc1(Nf, &cellsp);CHKERRQ(ierr);} 565 else {cellsp = sp;} 566 if (localU && localU != localX) {ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);CHKERRQ(ierr);} 567 /* Get cell dual spaces */ 568 for (f = 0; f < Nf; ++f) { 569 PetscDiscType disctype; 570 571 ierr = PetscDSGetDiscType_Internal(ds, f, &disctype);CHKERRQ(ierr); 572 if (disctype == PETSC_DISC_FE) { 573 PetscFE fe; 574 575 isFE[f] = PETSC_TRUE; 576 hasFE = PETSC_TRUE; 577 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); 578 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 579 } else if (disctype == PETSC_DISC_FV) { 580 PetscFV fv; 581 582 isFE[f] = PETSC_FALSE; 583 hasFV = PETSC_TRUE; 584 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fv);CHKERRQ(ierr); 585 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 586 } else { 587 isFE[f] = PETSC_FALSE; 588 cellsp[f] = NULL; 589 } 590 } 591 ierr = DMGetCoordinateField(dm,&coordField);CHKERRQ(ierr); 592 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) { 593 PetscInt effectiveHeight = auxBd ? minHeight : 0; 594 PetscFE fem, subfem; 595 PetscDiscType disctype; 596 const PetscReal *points; 597 PetscInt numPoints; 598 599 if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation"); 600 for (f = 0; f < Nf; ++f) { 601 if (!effectiveHeight) {sp[f] = cellsp[f];} 602 else {ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr);} 603 } 604 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);CHKERRQ(ierr); 605 ierr = PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);CHKERRQ(ierr); 606 ierr = PetscMalloc2(NfIn, &T, NfAux, &TAux);CHKERRQ(ierr); 607 for (f = 0; f < NfIn; ++f) { 608 ierr = PetscDSGetDiscType_Internal(dsIn, f, &disctype);CHKERRQ(ierr); 609 if (disctype != PETSC_DISC_FE) continue; 610 ierr = PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);CHKERRQ(ierr); 611 if (!effectiveHeight) {subfem = fem;} 612 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 613 ierr = PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);CHKERRQ(ierr); 614 } 615 for (f = 0; f < NfAux; ++f) { 616 ierr = PetscDSGetDiscType_Internal(dsAux, f, &disctype);CHKERRQ(ierr); 617 if (disctype != PETSC_DISC_FE) continue; 618 ierr = PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 619 if (!effectiveHeight || auxBd) {subfem = fem;} 620 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 621 ierr = PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);CHKERRQ(ierr); 622 } 623 } 624 /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */ 625 for (h = minHeight; h <= maxHeight; h++) { 626 PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight); 627 PetscDS dsEff = ds; 628 PetscScalar *values; 629 PetscBool *fieldActive; 630 PetscInt maxDegree; 631 PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues; 632 IS heightIS; 633 634 /* Note we assume that dm and dmIn share the same topology */ 635 ierr = DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);CHKERRQ(ierr); 636 ierr = DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, h, &lStart, NULL);CHKERRQ(ierr); 637 ierr = DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);CHKERRQ(ierr); 638 if (pEnd <= pStart) { 639 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 640 continue; 641 } 642 /* Compute totDim, the number of dofs in the closure of a point at this height */ 643 totDim = 0; 644 for (f = 0; f < Nf; ++f) { 645 if (!effectiveHeight) { 646 sp[f] = cellsp[f]; 647 } else { 648 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr); 649 } 650 if (!sp[f]) continue; 651 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 652 totDim += spDim; 653 if (isHybrid && (f < Nf-1)) totDim += spDim; 654 } 655 ierr = DMPlexVecGetClosure(plex, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);CHKERRQ(ierr); 656 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim); 657 if (!totDim) { 658 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 659 continue; 660 } 661 if (effectiveHeight) {ierr = PetscDSGetHeightSubspace(ds, effectiveHeight, &dsEff);CHKERRQ(ierr);} 662 /* Loop over points at this height */ 663 ierr = DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 664 ierr = DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 665 for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 666 if (label) { 667 PetscInt i; 668 669 for (i = 0; i < numIds; ++i) { 670 IS pointIS, isectIS; 671 const PetscInt *points; 672 PetscInt n; 673 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 674 PetscQuadrature quad = NULL; 675 676 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 677 if (!pointIS) continue; /* No points with that id on this process */ 678 ierr = ISIntersect(pointIS,heightIS,&isectIS);CHKERRQ(ierr); 679 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 680 if (!isectIS) continue; 681 ierr = ISGetLocalSize(isectIS, &n);CHKERRQ(ierr); 682 ierr = ISGetIndices(isectIS, &points);CHKERRQ(ierr); 683 ierr = DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);CHKERRQ(ierr); 684 if (maxDegree <= 1) { 685 ierr = DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);CHKERRQ(ierr); 686 } 687 if (!quad) { 688 if (!h && allPoints) { 689 quad = allPoints; 690 allPoints = NULL; 691 } else { 692 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,isHybrid ? dim-h-1 : dim-h,funcs,&quad);CHKERRQ(ierr); 693 } 694 } 695 ierr = DMFieldCreateFEGeom(coordField,isectIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);CHKERRQ(ierr); 696 for (p = 0; p < n; ++p) { 697 const PetscInt point = points[p]; 698 699 ierr = PetscArrayzero(values, numValues);CHKERRQ(ierr); 700 ierr = PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 701 ierr = DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsIn, plexAux, encAux, dsAux, chunkgeom, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, point, T, TAux, type, funcs, ctxs, fieldActive, values); 702 if (ierr) { 703 PetscErrorCode ierr2; 704 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 705 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 706 CHKERRQ(ierr); 707 } 708 if (transform) {ierr = DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);CHKERRQ(ierr);} 709 ierr = DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);CHKERRQ(ierr); 710 } 711 ierr = PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 712 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 713 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 714 ierr = ISRestoreIndices(isectIS, &points);CHKERRQ(ierr); 715 ierr = ISDestroy(&isectIS);CHKERRQ(ierr); 716 } 717 } else { 718 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 719 PetscQuadrature quad = NULL; 720 IS pointIS; 721 722 ierr = ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);CHKERRQ(ierr); 723 ierr = DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);CHKERRQ(ierr); 724 if (maxDegree <= 1) { 725 ierr = DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);CHKERRQ(ierr); 726 } 727 if (!quad) { 728 if (!h && allPoints) { 729 quad = allPoints; 730 allPoints = NULL; 731 } else { 732 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);CHKERRQ(ierr); 733 } 734 } 735 ierr = DMFieldCreateFEGeom(coordField,pointIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);CHKERRQ(ierr); 736 for (p = pStart; p < pEnd; ++p) { 737 ierr = PetscArrayzero(values, numValues);CHKERRQ(ierr); 738 ierr = PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);CHKERRQ(ierr); 739 ierr = DMProjectPoint_Private(dm, dsEff, plexIn, encIn, dsIn, plexAux, encAux, dsAux, chunkgeom, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, p, T, TAux, type, funcs, ctxs, fieldActive, values); 740 if (ierr) { 741 PetscErrorCode ierr2; 742 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 743 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 744 CHKERRQ(ierr); 745 } 746 if (transform) {ierr = DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);CHKERRQ(ierr);} 747 ierr = DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);CHKERRQ(ierr); 748 } 749 ierr = PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);CHKERRQ(ierr); 750 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 751 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 752 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 753 } 754 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 755 ierr = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 756 ierr = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 757 } 758 /* Cleanup */ 759 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) { 760 PetscInt effectiveHeight = auxBd ? minHeight : 0; 761 PetscFE fem, subfem; 762 763 for (f = 0; f < NfIn; ++f) { 764 ierr = PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);CHKERRQ(ierr); 765 if (!effectiveHeight) {subfem = fem;} 766 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 767 ierr = PetscTabulationDestroy(&T[f]);CHKERRQ(ierr); 768 } 769 for (f = 0; f < NfAux; ++f) { 770 ierr = PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 771 if (!effectiveHeight || auxBd) {subfem = fem;} 772 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 773 ierr = PetscTabulationDestroy(&TAux[f]);CHKERRQ(ierr); 774 } 775 ierr = PetscFree2(T, TAux);CHKERRQ(ierr); 776 } 777 ierr = PetscQuadratureDestroy(&allPoints);CHKERRQ(ierr); 778 ierr = PetscFree2(isFE, sp);CHKERRQ(ierr); 779 if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 780 ierr = DMDestroy(&plex);CHKERRQ(ierr); 781 ierr = DMDestroy(&plexIn);CHKERRQ(ierr); 782 if (dmAux) {ierr = DMDestroy(&plexAux);CHKERRQ(ierr);} 783 PetscFunctionReturn(0); 784 } 785 786 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 787 { 788 PetscErrorCode ierr; 789 790 PetscFunctionBegin; 791 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 792 PetscFunctionReturn(0); 793 } 794 795 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) 796 { 797 PetscErrorCode ierr; 798 799 PetscFunctionBegin; 800 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 801 PetscFunctionReturn(0); 802 } 803 804 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 805 void (**funcs)(PetscInt, PetscInt, PetscInt, 806 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 807 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 808 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 809 InsertMode mode, Vec localX) 810 { 811 PetscErrorCode ierr; 812 813 PetscFunctionBegin; 814 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 815 PetscFunctionReturn(0); 816 } 817 818 PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 819 void (**funcs)(PetscInt, PetscInt, PetscInt, 820 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 821 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 822 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 823 InsertMode mode, Vec localX) 824 { 825 PetscErrorCode ierr; 826 827 PetscFunctionBegin; 828 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 829 PetscFunctionReturn(0); 830 } 831 832 PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 833 void (**funcs)(PetscInt, PetscInt, PetscInt, 834 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 835 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 836 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 837 InsertMode mode, Vec localX) 838 { 839 PetscErrorCode ierr; 840 841 PetscFunctionBegin; 842 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_BD_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 843 PetscFunctionReturn(0); 844 } 845