147923291SMatthew G. Knepley #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 247923291SMatthew G. Knepley 38c6c5593SMatthew G. Knepley #include <petsc/private/petscfeimpl.h> 48c6c5593SMatthew G. Knepley 5496733e2SMatthew G. Knepley static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscReal time, PetscFECellGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[], 68c6c5593SMatthew G. Knepley PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, 78c6c5593SMatthew G. Knepley PetscScalar values[]) 847923291SMatthew G. Knepley { 98c6c5593SMatthew G. Knepley PetscDS prob; 10496733e2SMatthew G. Knepley PetscInt Nf, *Nc, f, spDim, d, c, v; 118c6c5593SMatthew G. Knepley PetscErrorCode ierr; 128c6c5593SMatthew G. Knepley 138c6c5593SMatthew G. Knepley PetscFunctionBeginHot; 148c6c5593SMatthew G. Knepley ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 158c6c5593SMatthew G. Knepley ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 16496733e2SMatthew G. Knepley ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 178c6c5593SMatthew G. Knepley /* Get values for closure */ 188c6c5593SMatthew G. Knepley for (f = 0, v = 0; f < Nf; ++f) { 198c6c5593SMatthew G. Knepley void * const ctx = ctxs ? ctxs[f] : NULL; 208c6c5593SMatthew G. Knepley 218c6c5593SMatthew G. Knepley if (!sp[f]) continue; 228c6c5593SMatthew G. Knepley ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 238c6c5593SMatthew G. Knepley for (d = 0; d < spDim; ++d) { 248c6c5593SMatthew G. Knepley if (funcs[f]) { 258c6c5593SMatthew G. Knepley if (isFE[f]) {ierr = PetscDualSpaceApply(sp[f], d, time, fegeom, Nc[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 268c6c5593SMatthew G. Knepley else {ierr = PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 278c6c5593SMatthew G. Knepley } else { 288c6c5593SMatthew G. Knepley for (c = 0; c < Nc[f]; ++c) values[v+c] = 0.0; 298c6c5593SMatthew G. Knepley } 308c6c5593SMatthew G. Knepley v += Nc[f]; 318c6c5593SMatthew G. Knepley } 328c6c5593SMatthew G. Knepley } 338c6c5593SMatthew G. Knepley PetscFunctionReturn(0); 348c6c5593SMatthew G. Knepley } 358c6c5593SMatthew G. Knepley 36496733e2SMatthew G. Knepley static PetscErrorCode DMProjectPoint_Field_Private(DM dm, DM dmAux, PetscReal time, Vec localU, Vec localA, PetscFECellGeom *fegeom, PetscDualSpace sp[], PetscInt p, 370c898c18SMatthew G. Knepley PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 388c6c5593SMatthew G. Knepley void (**funcs)(PetscInt, PetscInt, PetscInt, 398c6c5593SMatthew G. Knepley const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 408c6c5593SMatthew G. Knepley const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 418c6c5593SMatthew G. Knepley PetscReal, const PetscReal[], PetscScalar[]), void **ctxs, 428c6c5593SMatthew G. Knepley PetscScalar values[]) 438c6c5593SMatthew G. Knepley { 448c6c5593SMatthew G. Knepley PetscDS prob, probAux = NULL; 458c6c5593SMatthew G. Knepley PetscSection section, sectionAux = NULL; 46496733e2SMatthew G. Knepley PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *refSpaceDer, *refSpaceDerAux = NULL; 478c6c5593SMatthew G. Knepley PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 488c6c5593SMatthew G. Knepley PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL; 498c6c5593SMatthew G. Knepley PetscReal *x; 50496733e2SMatthew G. Knepley PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 511ba34526SMatthew G. Knepley PetscInt dim, dimAux = 0, Nf, NfAux = 0, f, spDim, d, c, v, tp = 0; 528c6c5593SMatthew G. Knepley PetscErrorCode ierr; 538c6c5593SMatthew G. Knepley 548c6c5593SMatthew G. Knepley PetscFunctionBeginHot; 558c6c5593SMatthew G. Knepley ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 568c6c5593SMatthew G. Knepley ierr = PetscDSGetSpatialDimension(prob, &dim);CHKERRQ(ierr); 578c6c5593SMatthew G. Knepley ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 58496733e2SMatthew G. Knepley ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 59496733e2SMatthew G. Knepley ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 608c6c5593SMatthew G. Knepley ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 618c6c5593SMatthew G. Knepley ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 6211d189daSMatthew G. Knepley ierr = PetscDSGetEvaluationArrays(prob, &u, NULL /*&u_t*/, &u_x);CHKERRQ(ierr); 63496733e2SMatthew G. Knepley ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 648c6c5593SMatthew G. Knepley ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 65496733e2SMatthew G. Knepley ierr = DMPlexVecGetClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 668c6c5593SMatthew G. Knepley if (dmAux) { 671ba34526SMatthew G. Knepley DMLabel spmap; 681ba34526SMatthew G. Knepley PetscInt subp; 691ba34526SMatthew G. Knepley 701ba34526SMatthew G. Knepley ierr = DMPlexGetSubpointMap(dmAux, &spmap);CHKERRQ(ierr); 711ba34526SMatthew G. Knepley if (spmap) { 721ba34526SMatthew G. Knepley IS subpointIS; 731ba34526SMatthew G. Knepley const PetscInt *subpoints; 741ba34526SMatthew G. Knepley PetscInt numSubpoints; 751ba34526SMatthew G. Knepley 761ba34526SMatthew G. Knepley ierr = DMPlexCreateSubpointIS(dmAux, &subpointIS);CHKERRQ(ierr); 771ba34526SMatthew G. Knepley ierr = ISGetLocalSize(subpointIS, &numSubpoints);CHKERRQ(ierr); 781ba34526SMatthew G. Knepley ierr = ISGetIndices(subpointIS, &subpoints);CHKERRQ(ierr); 791ba34526SMatthew G. Knepley ierr = PetscFindInt(p, numSubpoints, subpoints, &subp);CHKERRQ(ierr); 801ba34526SMatthew G. Knepley if (subp < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", p); 811ba34526SMatthew G. Knepley ierr = ISRestoreIndices(subpointIS, &subpoints);CHKERRQ(ierr); 821ba34526SMatthew G. Knepley ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 831ba34526SMatthew G. Knepley } else subp = p; 848c6c5593SMatthew G. Knepley ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 851ba34526SMatthew G. Knepley ierr = PetscDSGetSpatialDimension(probAux, &dimAux);CHKERRQ(ierr); 868c6c5593SMatthew G. Knepley ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 87496733e2SMatthew G. Knepley ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 88496733e2SMatthew G. Knepley ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 898c6c5593SMatthew G. Knepley ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 908c6c5593SMatthew G. Knepley ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 918c6c5593SMatthew G. Knepley ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 9211d189daSMatthew G. Knepley ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);CHKERRQ(ierr); 93496733e2SMatthew G. Knepley ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 941ba34526SMatthew G. Knepley ierr = DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);CHKERRQ(ierr); 958c6c5593SMatthew G. Knepley } 968c6c5593SMatthew G. Knepley /* Get values for closure */ 978c6c5593SMatthew G. Knepley for (f = 0, v = 0; f < Nf; ++f) { 988c6c5593SMatthew G. Knepley if (!sp[f]) continue; 998c6c5593SMatthew G. Knepley ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 1008c6c5593SMatthew G. Knepley for (d = 0; d < spDim; ++d) { 1018c6c5593SMatthew G. Knepley if (funcs[f]) { 1028c6c5593SMatthew G. Knepley PetscQuadrature quad; 1030c898c18SMatthew G. Knepley const PetscReal *points; 1048c6c5593SMatthew G. Knepley PetscInt numPoints, q; 1058c6c5593SMatthew G. Knepley 1068c6c5593SMatthew G. Knepley ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 1070c898c18SMatthew G. Knepley ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, NULL);CHKERRQ(ierr); 1080c898c18SMatthew G. Knepley for (q = 0; q < numPoints; ++q, ++tp) { 1098c6c5593SMatthew G. Knepley CoordinatesRefToReal(dim, dim, fegeom->v0, fegeom->J, &points[q*dim], x); 110496733e2SMatthew G. Knepley EvaluateFieldJets(dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, refSpaceDer, fegeom->invJ, coefficients, coefficients_t, u, u_x, u_t); 1111ba34526SMatthew G. Knepley if (probAux) {EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, tp, basisTabAux, basisDerTabAux, refSpaceDerAux, fegeom->invJ, coefficientsAux, coefficientsAux_t, a, a_x, a_t);} 1128c6c5593SMatthew G. Knepley (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, x, &values[v]); 1138c6c5593SMatthew G. Knepley } 1148c6c5593SMatthew G. Knepley } else { 1158c6c5593SMatthew G. Knepley for (c = 0; c < Nc[f]; ++c) values[v+c] = 0.0; 1168c6c5593SMatthew G. Knepley } 1178c6c5593SMatthew G. Knepley v += Nc[f]; 1188c6c5593SMatthew G. Knepley } 1198c6c5593SMatthew G. Knepley } 1208c6c5593SMatthew G. Knepley ierr = DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 1218c6c5593SMatthew G. Knepley if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);CHKERRQ(ierr);} 1228c6c5593SMatthew G. Knepley PetscFunctionReturn(0); 1238c6c5593SMatthew G. Knepley } 1248c6c5593SMatthew G. Knepley 1250c898c18SMatthew G. Knepley static PetscErrorCode DMProjectPoint_Private(DM dm, DM dmAux, PetscInt h, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[], 126496733e2SMatthew G. Knepley PetscDualSpace sp[], PetscInt p, PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 1278c6c5593SMatthew G. Knepley DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[]) 1288c6c5593SMatthew G. Knepley { 1298c6c5593SMatthew G. Knepley PetscFECellGeom fegeom; 1308c6c5593SMatthew G. Knepley PetscFVCellGeom fvgeom; 1318c6c5593SMatthew G. Knepley PetscInt dim, dimEmbed; 1328c6c5593SMatthew G. Knepley PetscErrorCode ierr; 1338c6c5593SMatthew G. Knepley 1348c6c5593SMatthew G. Knepley PetscFunctionBeginHot; 1358c6c5593SMatthew G. Knepley ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1368c6c5593SMatthew G. Knepley ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 1378c6c5593SMatthew G. Knepley if (hasFE) { 13843048b12SMatthew G. Knepley ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, fegeom.invJ, &fegeom.detJ);CHKERRQ(ierr); 1398c6c5593SMatthew G. Knepley fegeom.dim = dim - h; 1408c6c5593SMatthew G. Knepley fegeom.dimEmbed = dimEmbed; 1418c6c5593SMatthew G. Knepley } 1428c6c5593SMatthew G. Knepley if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 1438c6c5593SMatthew G. Knepley switch (type) { 1448c6c5593SMatthew G. Knepley case DM_BC_ESSENTIAL: 1458c6c5593SMatthew G. Knepley case DM_BC_NATURAL: 146496733e2SMatthew G. Knepley ierr = DMProjectPoint_Func_Private(dm, time, &fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);CHKERRQ(ierr);break; 1478c6c5593SMatthew G. Knepley case DM_BC_ESSENTIAL_FIELD: 1488c6c5593SMatthew G. Knepley case DM_BC_NATURAL_FIELD: 149496733e2SMatthew G. Knepley ierr = DMProjectPoint_Field_Private(dm, dmAux, time, localU, localA, &fegeom, sp, p, 1500c898c18SMatthew G. Knepley basisTab, basisDerTab, basisTabAux, basisDerTabAux, 1510c898c18SMatthew G. Knepley (void (**)(PetscInt, PetscInt, PetscInt, 1520c898c18SMatthew G. Knepley const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1530c898c18SMatthew G. Knepley const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1540c898c18SMatthew G. Knepley PetscReal, const PetscReal[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 1558c6c5593SMatthew G. Knepley default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type); 1568c6c5593SMatthew G. Knepley } 1578c6c5593SMatthew G. Knepley PetscFunctionReturn(0); 1588c6c5593SMatthew G. Knepley } 1598c6c5593SMatthew G. Knepley 1608c6c5593SMatthew G. Knepley PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, 1618c6c5593SMatthew G. Knepley DMLabel label, PetscInt numIds, const PetscInt ids[], 1628c6c5593SMatthew G. Knepley DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, 1638c6c5593SMatthew G. Knepley InsertMode mode, Vec localX) 1648c6c5593SMatthew G. Knepley { 1658c6c5593SMatthew G. Knepley DM dmAux = NULL; 1660c898c18SMatthew G. Knepley PetscDS prob, probAux = NULL; 1678c6c5593SMatthew G. Knepley Vec localA = NULL; 16847923291SMatthew G. Knepley PetscSection section; 1698c6c5593SMatthew G. Knepley PetscDualSpace *sp, *cellsp; 1700c898c18SMatthew G. Knepley PetscReal **basisTab = NULL, **basisDerTab = NULL, **basisTabAux = NULL, **basisDerTabAux = NULL; 1718c6c5593SMatthew G. Knepley PetscInt *Nc; 1722af58022SMatthew G. Knepley PetscInt dim, dimEmbed, minHeight, maxHeight, h, Nf, NfAux = 0, f; 17347923291SMatthew G. Knepley PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE; 17447923291SMatthew G. Knepley PetscErrorCode ierr; 17547923291SMatthew G. Knepley 17647923291SMatthew G. Knepley PetscFunctionBegin; 1778c6c5593SMatthew G. Knepley ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1782af58022SMatthew G. Knepley ierr = DMPlexGetVTKCellHeight(dm, &minHeight);CHKERRQ(ierr); 1798c6c5593SMatthew G. Knepley ierr = DMPlexGetMaxProjectionHeight(dm, &maxHeight);CHKERRQ(ierr); 1802af58022SMatthew G. Knepley maxHeight = PetscMax(maxHeight, minHeight); 1818c6c5593SMatthew G. Knepley if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %d not in [0, %d)\n", maxHeight, dim);} 1820c898c18SMatthew G. Knepley ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1838c6c5593SMatthew G. Knepley ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 18447923291SMatthew G. Knepley ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 18547923291SMatthew G. Knepley ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1868c6c5593SMatthew G. Knepley ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1878c6c5593SMatthew G. Knepley ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);CHKERRQ(ierr); 1880c898c18SMatthew G. Knepley if (dmAux) { 1890c898c18SMatthew G. Knepley ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1900c898c18SMatthew G. Knepley ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1910c898c18SMatthew G. Knepley } 192496733e2SMatthew G. Knepley ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 193496733e2SMatthew G. Knepley ierr = PetscMalloc2(Nf, &isFE, Nf, &sp);CHKERRQ(ierr); 1948c6c5593SMatthew G. Knepley if (maxHeight > 0) {ierr = PetscMalloc1(Nf, &cellsp);CHKERRQ(ierr);} 1958c6c5593SMatthew G. Knepley else {cellsp = sp;} 1968c6c5593SMatthew G. Knepley if (localU && localU != localX) {ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);CHKERRQ(ierr);} 1978c6c5593SMatthew G. Knepley /* Get cell dual spaces */ 19847923291SMatthew G. Knepley for (f = 0; f < Nf; ++f) { 19947923291SMatthew G. Knepley PetscObject obj; 20047923291SMatthew G. Knepley PetscClassId id; 20147923291SMatthew G. Knepley 20247923291SMatthew G. Knepley ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 20347923291SMatthew G. Knepley ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 20447923291SMatthew G. Knepley if (id == PETSCFE_CLASSID) { 20547923291SMatthew G. Knepley PetscFE fe = (PetscFE) obj; 20647923291SMatthew G. Knepley 20747923291SMatthew G. Knepley hasFE = PETSC_TRUE; 20847923291SMatthew G. Knepley isFE[f] = PETSC_TRUE; 20947923291SMatthew G. Knepley ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 21047923291SMatthew G. Knepley } else if (id == PETSCFV_CLASSID) { 21147923291SMatthew G. Knepley PetscFV fv = (PetscFV) obj; 21247923291SMatthew G. Knepley 21347923291SMatthew G. Knepley hasFV = PETSC_TRUE; 21447923291SMatthew G. Knepley isFE[f] = PETSC_FALSE; 21547923291SMatthew G. Knepley ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 21647923291SMatthew G. Knepley } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 21747923291SMatthew G. Knepley } 2180c898c18SMatthew G. Knepley if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 2190c898c18SMatthew G. Knepley PetscFE fem; 2200c898c18SMatthew G. Knepley PetscReal *points; 2210c898c18SMatthew G. Knepley PetscInt numPoints, spDim, d; 2220c898c18SMatthew G. Knepley 2232af58022SMatthew G. Knepley if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation"); 2240c898c18SMatthew G. Knepley numPoints = 0; 2250c898c18SMatthew G. Knepley for (f = 0; f < Nf; ++f) { 2260c898c18SMatthew G. Knepley ierr = PetscDualSpaceGetDimension(cellsp[f], &spDim);CHKERRQ(ierr); 2270c898c18SMatthew G. Knepley for (d = 0; d < spDim; ++d) { 2280c898c18SMatthew G. Knepley if (funcs[f]) { 2290c898c18SMatthew G. Knepley PetscQuadrature quad; 2300c898c18SMatthew G. Knepley PetscInt Nq; 2310c898c18SMatthew G. Knepley 2320c898c18SMatthew G. Knepley ierr = PetscDualSpaceGetFunctional(cellsp[f], d, &quad);CHKERRQ(ierr); 2330c898c18SMatthew G. Knepley ierr = PetscQuadratureGetData(quad, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2340c898c18SMatthew G. Knepley numPoints += Nq; 2350c898c18SMatthew G. Knepley } 2360c898c18SMatthew G. Knepley } 2370c898c18SMatthew G. Knepley } 2380c898c18SMatthew G. Knepley ierr = PetscMalloc1(numPoints*dim, &points);CHKERRQ(ierr); 2390c898c18SMatthew G. Knepley numPoints = 0; 2400c898c18SMatthew G. Knepley for (f = 0; f < Nf; ++f) { 2410c898c18SMatthew G. Knepley ierr = PetscDualSpaceGetDimension(cellsp[f], &spDim);CHKERRQ(ierr); 2420c898c18SMatthew G. Knepley for (d = 0; d < spDim; ++d) { 2430c898c18SMatthew G. Knepley if (funcs[f]) { 2440c898c18SMatthew G. Knepley PetscQuadrature quad; 2450c898c18SMatthew G. Knepley const PetscReal *qpoints; 2460c898c18SMatthew G. Knepley PetscInt Nq, q; 2470c898c18SMatthew G. Knepley 2480c898c18SMatthew G. Knepley ierr = PetscDualSpaceGetFunctional(cellsp[f], d, &quad);CHKERRQ(ierr); 2490c898c18SMatthew G. Knepley ierr = PetscQuadratureGetData(quad, NULL, &Nq, &qpoints, NULL);CHKERRQ(ierr); 2500c898c18SMatthew G. Knepley for (q = 0; q < Nq*dim; ++q) points[numPoints*dim+q] = qpoints[q]; 2510c898c18SMatthew G. Knepley numPoints += Nq; 2520c898c18SMatthew G. Knepley } 2530c898c18SMatthew G. Knepley } 2540c898c18SMatthew G. Knepley } 2550c898c18SMatthew G. Knepley ierr = PetscMalloc4(Nf, &basisTab, Nf, &basisDerTab, NfAux, &basisTabAux, NfAux, &basisDerTabAux);CHKERRQ(ierr); 2560c898c18SMatthew G. Knepley for (f = 0; f < Nf; ++f) { 2570c898c18SMatthew G. Knepley if (!isFE[f]) continue; 2580c898c18SMatthew G. Knepley ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 2590c898c18SMatthew G. Knepley ierr = PetscFEGetTabulation(fem, numPoints, points, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 2600c898c18SMatthew G. Knepley } 2610c898c18SMatthew G. Knepley for (f = 0; f < NfAux; ++f) { 2620c898c18SMatthew G. Knepley ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 2630c898c18SMatthew G. Knepley ierr = PetscFEGetTabulation(fem, numPoints, points, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 2640c898c18SMatthew G. Knepley } 2650c898c18SMatthew G. Knepley ierr = PetscFree(points);CHKERRQ(ierr); 2660c898c18SMatthew G. Knepley } 26747923291SMatthew G. Knepley /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */ 2682af58022SMatthew G. Knepley for (h = minHeight; h <= maxHeight; h++) { 269*bd1e2261SMatthew G. Knepley PetscInt effectiveHeight = h - minHeight; 2708c6c5593SMatthew G. Knepley PetscScalar *values; 2718c6c5593SMatthew G. Knepley PetscBool *fieldActive; 27284f7fa7aSMatthew G. Knepley PetscInt pStart, pEnd, p, spDim, totDim, numValues; 2738c6c5593SMatthew G. Knepley 27447923291SMatthew G. Knepley ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 2758c6c5593SMatthew G. Knepley if (!h) { 2768c6c5593SMatthew G. Knepley PetscInt cEndInterior; 2778c6c5593SMatthew G. Knepley 2788c6c5593SMatthew G. Knepley ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2798c6c5593SMatthew G. Knepley pEnd = cEndInterior < 0 ? pEnd : cEndInterior; 2808c6c5593SMatthew G. Knepley } 28147923291SMatthew G. Knepley if (pEnd <= pStart) continue; 28247923291SMatthew G. Knepley /* Compute totDim, the number of dofs in the closure of a point at this height */ 28347923291SMatthew G. Knepley totDim = 0; 28447923291SMatthew G. Knepley for (f = 0; f < Nf; ++f) { 285bec1f4aaSMatthew G. Knepley if (!effectiveHeight) { 28647923291SMatthew G. Knepley sp[f] = cellsp[f]; 28747923291SMatthew G. Knepley } else { 288bec1f4aaSMatthew G. Knepley ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr); 28947923291SMatthew G. Knepley if (!sp[f]) continue; 29047923291SMatthew G. Knepley } 29147923291SMatthew G. Knepley ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 2928c6c5593SMatthew G. Knepley totDim += spDim*Nc[f]; 29347923291SMatthew G. Knepley } 29447923291SMatthew G. Knepley ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 2958c6c5593SMatthew G. Knepley if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim); 29647923291SMatthew G. Knepley if (!totDim) continue; 29747923291SMatthew G. Knepley /* Loop over points at this height */ 29847923291SMatthew G. Knepley ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 2998c6c5593SMatthew G. Knepley ierr = DMGetWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 3008c6c5593SMatthew G. Knepley for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 3018c6c5593SMatthew G. Knepley if (label) { 3028c6c5593SMatthew G. Knepley PetscInt i; 30347923291SMatthew G. Knepley 30447923291SMatthew G. Knepley for (i = 0; i < numIds; ++i) { 30547923291SMatthew G. Knepley IS pointIS; 30647923291SMatthew G. Knepley const PetscInt *points; 3078c6c5593SMatthew G. Knepley PetscInt n; 30847923291SMatthew G. Knepley 30947923291SMatthew G. Knepley ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 31047923291SMatthew G. Knepley if (!pointIS) continue; /* No points with that id on this process */ 31147923291SMatthew G. Knepley ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 31247923291SMatthew G. Knepley ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 31347923291SMatthew G. Knepley for (p = 0; p < n; ++p) { 31447923291SMatthew G. Knepley const PetscInt point = points[p]; 31547923291SMatthew G. Knepley 31647923291SMatthew G. Knepley if ((point < pStart) || (point >= pEnd)) continue; 317*bd1e2261SMatthew G. Knepley ierr = DMProjectPoint_Private(dm, dmAux, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, point, basisTab, basisDerTab, basisTabAux, basisDerTabAux, type, funcs, ctxs, fieldActive, values); 31847923291SMatthew G. Knepley if (ierr) { 31947923291SMatthew G. Knepley PetscErrorCode ierr2; 32047923291SMatthew G. Knepley ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2); 3218c6c5593SMatthew G. Knepley ierr2 = DMRestoreWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr2); 32247923291SMatthew G. Knepley CHKERRQ(ierr); 32347923291SMatthew G. Knepley } 32447923291SMatthew G. Knepley ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr); 32547923291SMatthew G. Knepley } 32647923291SMatthew G. Knepley ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 32747923291SMatthew G. Knepley ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 32847923291SMatthew G. Knepley } 3298c6c5593SMatthew G. Knepley } else { 3308c6c5593SMatthew G. Knepley for (p = pStart; p < pEnd; ++p) { 331*bd1e2261SMatthew G. Knepley ierr = DMProjectPoint_Private(dm, dmAux, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, p, basisTab, basisDerTab, basisTabAux, basisDerTabAux, type, funcs, ctxs, fieldActive, values); 3328c6c5593SMatthew G. Knepley if (ierr) { 3338c6c5593SMatthew G. Knepley PetscErrorCode ierr2; 3348c6c5593SMatthew G. Knepley ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2); 3358c6c5593SMatthew G. Knepley ierr2 = DMRestoreWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr2); 3368c6c5593SMatthew G. Knepley CHKERRQ(ierr); 3378c6c5593SMatthew G. Knepley } 3388c6c5593SMatthew G. Knepley ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, values, mode);CHKERRQ(ierr); 3398c6c5593SMatthew G. Knepley } 3408c6c5593SMatthew G. Knepley } 34147923291SMatthew G. Knepley ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 3428c6c5593SMatthew G. Knepley ierr = DMRestoreWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 34347923291SMatthew G. Knepley } 3448c6c5593SMatthew G. Knepley /* Cleanup */ 3450c898c18SMatthew G. Knepley if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 3460c898c18SMatthew G. Knepley PetscFE fem; 3470c898c18SMatthew G. Knepley 3480c898c18SMatthew G. Knepley for (f = 0; f < Nf; ++f) { 3490c898c18SMatthew G. Knepley if (!isFE[f]) continue; 3500c898c18SMatthew G. Knepley ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 3510c898c18SMatthew G. Knepley ierr = PetscFERestoreTabulation(fem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 3520c898c18SMatthew G. Knepley } 3530c898c18SMatthew G. Knepley for (f = 0; f < NfAux; ++f) { 3540c898c18SMatthew G. Knepley ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 3550c898c18SMatthew G. Knepley ierr = PetscFERestoreTabulation(fem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 3560c898c18SMatthew G. Knepley } 3570c898c18SMatthew G. Knepley ierr = PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);CHKERRQ(ierr); 3580c898c18SMatthew G. Knepley } 359496733e2SMatthew G. Knepley ierr = PetscFree2(isFE, sp);CHKERRQ(ierr); 3608c6c5593SMatthew G. Knepley if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 3618c6c5593SMatthew G. Knepley PetscFunctionReturn(0); 36247923291SMatthew G. Knepley } 3638c6c5593SMatthew G. Knepley 3648c6c5593SMatthew G. Knepley PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 3658c6c5593SMatthew G. Knepley { 3668c6c5593SMatthew G. Knepley PetscErrorCode ierr; 3678c6c5593SMatthew G. Knepley 3688c6c5593SMatthew G. Knepley PetscFunctionBegin; 3698c6c5593SMatthew G. Knepley ierr = DMProjectLocal_Generic_Plex(dm, time, localX, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 3708c6c5593SMatthew G. Knepley PetscFunctionReturn(0); 3718c6c5593SMatthew G. Knepley } 3728c6c5593SMatthew G. Knepley 3738c6c5593SMatthew G. Knepley PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 3748c6c5593SMatthew G. Knepley { 3758c6c5593SMatthew G. Knepley PetscErrorCode ierr; 3768c6c5593SMatthew G. Knepley 3778c6c5593SMatthew G. Knepley PetscFunctionBegin; 3788c6c5593SMatthew G. Knepley ierr = DMProjectLocal_Generic_Plex(dm, time, localX, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 37947923291SMatthew G. Knepley PetscFunctionReturn(0); 38047923291SMatthew G. Knepley } 38147923291SMatthew G. Knepley 3828c6c5593SMatthew G. Knepley PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 38347923291SMatthew G. Knepley void (**funcs)(PetscInt, PetscInt, PetscInt, 38447923291SMatthew G. Knepley const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 38547923291SMatthew G. Knepley const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 38647923291SMatthew G. Knepley PetscReal, const PetscReal[], PetscScalar[]), 38747923291SMatthew G. Knepley InsertMode mode, Vec localX) 38847923291SMatthew G. Knepley { 38947923291SMatthew G. Knepley PetscErrorCode ierr; 39047923291SMatthew G. Knepley 39147923291SMatthew G. Knepley PetscFunctionBegin; 3928c6c5593SMatthew G. Knepley ierr = DMProjectLocal_Generic_Plex(dm, time, localU, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 3938c6c5593SMatthew G. Knepley PetscFunctionReturn(0); 39447923291SMatthew G. Knepley } 39547923291SMatthew G. Knepley 3968c6c5593SMatthew G. Knepley PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], Vec localU, 3978c6c5593SMatthew G. Knepley void (**funcs)(PetscInt, PetscInt, PetscInt, 3988c6c5593SMatthew G. Knepley const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 3998c6c5593SMatthew G. Knepley const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 4008c6c5593SMatthew G. Knepley PetscReal, const PetscReal[], PetscScalar[]), 4018c6c5593SMatthew G. Knepley InsertMode mode, Vec localX) 4028c6c5593SMatthew G. Knepley { 4038c6c5593SMatthew G. Knepley PetscErrorCode ierr; 40447923291SMatthew G. Knepley 4058c6c5593SMatthew G. Knepley PetscFunctionBegin; 4068c6c5593SMatthew G. Knepley ierr = DMProjectLocal_Generic_Plex(dm, time, localU, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 40747923291SMatthew G. Knepley PetscFunctionReturn(0); 40847923291SMatthew G. Knepley } 409