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 = 0; 329 fegeom.J = 0; 330 fegeom.v = 0; 331 fegeom.xi = 0; 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 PetscSection section; 543 PetscDualSpace *sp, *cellsp; 544 PetscTabulation *T = NULL, *TAux = NULL; 545 PetscInt *Nc; 546 PetscInt dim, dimEmbed, depth, minHeight, maxHeight, h, Nf, NfIn, NfAux = 0, f; 547 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE, isHybrid = PETSC_FALSE, transform; 548 DMField coordField; 549 DMLabel depthLabel; 550 PetscQuadrature allPoints = NULL; 551 PetscErrorCode ierr; 552 553 PetscFunctionBegin; 554 if (localU) {ierr = VecGetDM(localU, &dmIn);CHKERRQ(ierr);} 555 else {dmIn = dm;} 556 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 557 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);CHKERRQ(ierr); 558 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 559 ierr = DMConvert(dmIn, DMPLEX, &plexIn);CHKERRQ(ierr); 560 ierr = DMGetEnclosureRelation(dmIn, dm, &encIn);CHKERRQ(ierr); 561 ierr = DMGetEnclosureRelation(dmAux, dm, &encAux);CHKERRQ(ierr); 562 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 563 ierr = DMPlexGetVTKCellHeight(plex, &minHeight);CHKERRQ(ierr); 564 ierr = DMGetBasisTransformDM_Internal(dm, &tdm);CHKERRQ(ierr); 565 ierr = DMGetBasisTransformVec_Internal(dm, &tv);CHKERRQ(ierr); 566 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 567 /* Auxiliary information can only be used with interpolation of field functions */ 568 if (dmAux) { 569 ierr = DMConvert(dmAux, DMPLEX, &plexAux);CHKERRQ(ierr); 570 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) { 571 if (!localA) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Missing localA vector"); 572 if (!minHeight) { 573 DMLabel spmap; 574 575 /* If dmAux is a surface, then force the projection to take place over a surface */ 576 ierr = DMPlexGetSubpointMap(plexAux, &spmap);CHKERRQ(ierr); 577 if (spmap) { 578 ierr = DMPlexGetVTKCellHeight(plexAux, &minHeight);CHKERRQ(ierr); 579 auxBd = minHeight ? PETSC_TRUE : PETSC_FALSE; 580 } 581 } 582 } 583 } 584 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr); 585 ierr = DMPlexGetDepthLabel(plex, &depthLabel);CHKERRQ(ierr); 586 ierr = DMPlexGetMaxProjectionHeight(plex, &maxHeight);CHKERRQ(ierr); 587 maxHeight = PetscMax(maxHeight, minHeight); 588 if (maxHeight < 0 || maxHeight > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %D not in [0, %D)", maxHeight, dim); 589 ierr = DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, 0, NULL, &ds);CHKERRQ(ierr); 590 if (!ds) {ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);} 591 ierr = DMGetFirstLabelEntry_Private(dmIn, dm, label, numIds, ids, 0, NULL, &dsIn);CHKERRQ(ierr); 592 if (!dsIn) {ierr = DMGetDS(dmIn, &dsIn);CHKERRQ(ierr);} 593 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 594 ierr = PetscDSGetNumFields(dsIn, &NfIn);CHKERRQ(ierr); 595 ierr = PetscDSGetHybrid(ds, &isHybrid);CHKERRQ(ierr); 596 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 597 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 598 if (dmAux) { 599 ierr = DMGetDS(dmAux, &dsAux);CHKERRQ(ierr); 600 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 601 } 602 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 603 ierr = PetscMalloc2(Nf, &isFE, Nf, &sp);CHKERRQ(ierr); 604 if (maxHeight > 0) {ierr = PetscMalloc1(Nf, &cellsp);CHKERRQ(ierr);} 605 else {cellsp = sp;} 606 if (localU && localU != localX) {ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, localU, time, NULL, NULL, NULL);CHKERRQ(ierr);} 607 /* Get cell dual spaces */ 608 for (f = 0; f < Nf; ++f) { 609 PetscDiscType disctype; 610 611 ierr = PetscDSGetDiscType_Internal(ds, f, &disctype);CHKERRQ(ierr); 612 if (disctype == PETSC_DISC_FE) { 613 PetscFE fe; 614 615 isFE[f] = PETSC_TRUE; 616 hasFE = PETSC_TRUE; 617 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); 618 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 619 } else if (disctype == PETSC_DISC_FV) { 620 PetscFV fv; 621 622 isFE[f] = PETSC_FALSE; 623 hasFV = PETSC_TRUE; 624 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fv);CHKERRQ(ierr); 625 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 626 } else { 627 isFE[f] = PETSC_FALSE; 628 cellsp[f] = NULL; 629 } 630 } 631 ierr = DMGetCoordinateField(dm,&coordField);CHKERRQ(ierr); 632 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) { 633 PetscInt effectiveHeight = auxBd ? minHeight : 0; 634 PetscFE fem, subfem; 635 PetscDiscType disctype; 636 const PetscReal *points; 637 PetscInt numPoints; 638 639 if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation"); 640 for (f = 0; f < Nf; ++f) { 641 if (!effectiveHeight) {sp[f] = cellsp[f];} 642 else {ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr);} 643 } 644 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-effectiveHeight,funcs,&allPoints);CHKERRQ(ierr); 645 ierr = PetscQuadratureGetData(allPoints,NULL,NULL,&numPoints,&points,NULL);CHKERRQ(ierr); 646 ierr = PetscMalloc2(NfIn, &T, NfAux, &TAux);CHKERRQ(ierr); 647 for (f = 0; f < NfIn; ++f) { 648 ierr = PetscDSGetDiscType_Internal(dsIn, f, &disctype);CHKERRQ(ierr); 649 if (disctype != PETSC_DISC_FE) continue; 650 ierr = PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);CHKERRQ(ierr); 651 if (!effectiveHeight) {subfem = fem;} 652 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 653 ierr = PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &T[f]);CHKERRQ(ierr); 654 } 655 for (f = 0; f < NfAux; ++f) { 656 ierr = PetscDSGetDiscType_Internal(dsAux, f, &disctype);CHKERRQ(ierr); 657 if (disctype != PETSC_DISC_FE) continue; 658 ierr = PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 659 if (!effectiveHeight || auxBd) {subfem = fem;} 660 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 661 ierr = PetscFECreateTabulation(subfem, 1, numPoints, points, 1, &TAux[f]);CHKERRQ(ierr); 662 } 663 } 664 /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */ 665 for (h = minHeight; h <= maxHeight; h++) { 666 PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight); 667 PetscDS dsEff = ds; 668 PetscScalar *values; 669 PetscBool *fieldActive; 670 PetscInt maxDegree; 671 PetscInt pStart, pEnd, p, lStart, spDim, totDim, numValues; 672 IS heightIS; 673 674 /* Note we assume that dm and dmIn share the same topology */ 675 ierr = DMPlexGetSimplexOrBoxCells(plex, h, &pStart, &pEnd);CHKERRQ(ierr); 676 ierr = DMGetFirstLabelEntry_Private(dm, dm, label, numIds, ids, h, &lStart, NULL);CHKERRQ(ierr); 677 ierr = DMLabelGetStratumIS(depthLabel, depth - h, &heightIS);CHKERRQ(ierr); 678 if (pEnd <= pStart) { 679 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 680 continue; 681 } 682 /* Compute totDim, the number of dofs in the closure of a point at this height */ 683 totDim = 0; 684 for (f = 0; f < Nf; ++f) { 685 if (!effectiveHeight) { 686 sp[f] = cellsp[f]; 687 } else { 688 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr); 689 } 690 if (!sp[f]) continue; 691 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 692 totDim += spDim; 693 if (isHybrid && (f < Nf-1)) totDim += spDim; 694 } 695 ierr = DMPlexVecGetClosure(plex, section, localX, lStart < 0 ? pStart : lStart, &numValues, NULL);CHKERRQ(ierr); 696 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); 697 if (!totDim) { 698 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 699 continue; 700 } 701 if (effectiveHeight) {ierr = PetscDSGetHeightSubspace(ds, effectiveHeight, &dsEff);CHKERRQ(ierr);} 702 /* Loop over points at this height */ 703 ierr = DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 704 ierr = DMGetWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 705 for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 706 if (label) { 707 PetscInt i; 708 709 for (i = 0; i < numIds; ++i) { 710 IS pointIS, isectIS; 711 const PetscInt *points; 712 PetscInt n; 713 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 714 PetscQuadrature quad = NULL; 715 716 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 717 if (!pointIS) continue; /* No points with that id on this process */ 718 ierr = ISIntersect(pointIS,heightIS,&isectIS);CHKERRQ(ierr); 719 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 720 if (!isectIS) continue; 721 ierr = ISGetLocalSize(isectIS, &n);CHKERRQ(ierr); 722 ierr = ISGetIndices(isectIS, &points);CHKERRQ(ierr); 723 ierr = DMFieldGetDegree(coordField,isectIS,NULL,&maxDegree);CHKERRQ(ierr); 724 if (maxDegree <= 1) { 725 ierr = DMFieldCreateDefaultQuadrature(coordField,isectIS,&quad);CHKERRQ(ierr); 726 } 727 if (!quad) { 728 if (!h && allPoints) { 729 quad = allPoints; 730 allPoints = NULL; 731 } else { 732 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,isHybrid ? dim-h-1 : dim-h,funcs,&quad);CHKERRQ(ierr); 733 } 734 } 735 ierr = DMFieldCreateFEGeom(coordField,isectIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);CHKERRQ(ierr); 736 for (p = 0; p < n; ++p) { 737 const PetscInt point = points[p]; 738 739 ierr = PetscArrayzero(values, numValues);CHKERRQ(ierr); 740 ierr = PetscFEGeomGetChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 741 ierr = DMPlexSetActivePoint(dm, point);CHKERRQ(ierr); 742 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); 743 if (ierr) { 744 PetscErrorCode ierr2; 745 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 746 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 747 CHKERRQ(ierr); 748 } 749 if (transform) {ierr = DMPlexBasisTransformPoint_Internal(plex, tdm, tv, point, fieldActive, PETSC_FALSE, values);CHKERRQ(ierr);} 750 ierr = DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, point, Ncc, comps, label, ids[i], values, mode);CHKERRQ(ierr); 751 } 752 ierr = PetscFEGeomRestoreChunk(fegeom,p,p+1,&chunkgeom);CHKERRQ(ierr); 753 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 754 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 755 ierr = ISRestoreIndices(isectIS, &points);CHKERRQ(ierr); 756 ierr = ISDestroy(&isectIS);CHKERRQ(ierr); 757 } 758 } else { 759 PetscFEGeom *fegeom = NULL, *chunkgeom = NULL; 760 PetscQuadrature quad = NULL; 761 IS pointIS; 762 763 ierr = ISCreateStride(PETSC_COMM_SELF,pEnd-pStart,pStart,1,&pointIS);CHKERRQ(ierr); 764 ierr = DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);CHKERRQ(ierr); 765 if (maxDegree <= 1) { 766 ierr = DMFieldCreateDefaultQuadrature(coordField,pointIS,&quad);CHKERRQ(ierr); 767 } 768 if (!quad) { 769 if (!h && allPoints) { 770 quad = allPoints; 771 allPoints = NULL; 772 } else { 773 ierr = PetscDualSpaceGetAllPointsUnion(Nf,sp,dim-h,funcs,&quad);CHKERRQ(ierr); 774 } 775 } 776 ierr = DMFieldCreateFEGeom(coordField,pointIS,quad,(effectiveHeight && h == minHeight)?PETSC_TRUE:PETSC_FALSE,&fegeom);CHKERRQ(ierr); 777 for (p = pStart; p < pEnd; ++p) { 778 ierr = PetscArrayzero(values, numValues);CHKERRQ(ierr); 779 ierr = PetscFEGeomGetChunk(fegeom,p-pStart,p-pStart+1,&chunkgeom);CHKERRQ(ierr); 780 ierr = DMPlexSetActivePoint(dm, p);CHKERRQ(ierr); 781 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); 782 if (ierr) { 783 PetscErrorCode ierr2; 784 ierr2 = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr2); 785 ierr2 = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr2); 786 CHKERRQ(ierr); 787 } 788 if (transform) {ierr = DMPlexBasisTransformPoint_Internal(plex, tdm, tv, p, fieldActive, PETSC_FALSE, values);CHKERRQ(ierr);} 789 ierr = DMPlexVecSetFieldClosure_Internal(plex, section, localX, fieldActive, p, Ncc, comps, NULL, -1, values, mode);CHKERRQ(ierr); 790 } 791 ierr = PetscFEGeomRestoreChunk(fegeom,p-pStart,pStart-p+1,&chunkgeom);CHKERRQ(ierr); 792 ierr = PetscFEGeomDestroy(&fegeom);CHKERRQ(ierr); 793 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 794 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 795 } 796 ierr = ISDestroy(&heightIS);CHKERRQ(ierr); 797 ierr = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 798 ierr = DMRestoreWorkArray(dm, Nf, MPI_INT, &fieldActive);CHKERRQ(ierr); 799 } 800 /* Cleanup */ 801 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_ESSENTIAL_BD_FIELD || type == DM_BC_NATURAL_FIELD) { 802 PetscInt effectiveHeight = auxBd ? minHeight : 0; 803 PetscFE fem, subfem; 804 805 for (f = 0; f < NfIn; ++f) { 806 ierr = PetscDSGetDiscretization(dsIn, f, (PetscObject *) &fem);CHKERRQ(ierr); 807 if (!effectiveHeight) {subfem = fem;} 808 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 809 ierr = PetscTabulationDestroy(&T[f]);CHKERRQ(ierr); 810 } 811 for (f = 0; f < NfAux; ++f) { 812 ierr = PetscDSGetDiscretization(dsAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 813 if (!effectiveHeight || auxBd) {subfem = fem;} 814 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 815 ierr = PetscTabulationDestroy(&TAux[f]);CHKERRQ(ierr); 816 } 817 ierr = PetscFree2(T, TAux);CHKERRQ(ierr); 818 } 819 ierr = PetscQuadratureDestroy(&allPoints);CHKERRQ(ierr); 820 ierr = PetscFree2(isFE, sp);CHKERRQ(ierr); 821 if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 822 ierr = DMDestroy(&plex);CHKERRQ(ierr); 823 ierr = DMDestroy(&plexIn);CHKERRQ(ierr); 824 if (dmAux) {ierr = DMDestroy(&plexAux);CHKERRQ(ierr);} 825 PetscFunctionReturn(0); 826 } 827 828 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 829 { 830 PetscErrorCode ierr; 831 832 PetscFunctionBegin; 833 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 834 PetscFunctionReturn(0); 835 } 836 837 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) 838 { 839 PetscErrorCode ierr; 840 841 PetscFunctionBegin; 842 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 843 PetscFunctionReturn(0); 844 } 845 846 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 847 void (**funcs)(PetscInt, PetscInt, PetscInt, 848 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 849 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 850 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 851 InsertMode mode, Vec localX) 852 { 853 PetscErrorCode ierr; 854 855 PetscFunctionBegin; 856 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 857 PetscFunctionReturn(0); 858 } 859 860 PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 861 void (**funcs)(PetscInt, PetscInt, PetscInt, 862 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 863 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 864 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 865 InsertMode mode, Vec localX) 866 { 867 PetscErrorCode ierr; 868 869 PetscFunctionBegin; 870 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 871 PetscFunctionReturn(0); 872 } 873 874 PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 875 void (**funcs)(PetscInt, PetscInt, PetscInt, 876 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 877 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 878 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 879 InsertMode mode, Vec localX) 880 { 881 PetscErrorCode ierr; 882 883 PetscFunctionBegin; 884 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); 885 PetscFunctionReturn(0); 886 } 887