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