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