cellInterpolant.hpp

#ifndef ABLATELIBRARY_CELLINTERPOLANT_HPP
#define ABLATELIBRARY_CELLINTERPOLANT_HPP
 
#include <petsc.h>
#include <vector>
#include "domain/range.hpp"
#include "domain/region.hpp"
#include "domain/subDomain.hpp"
namespace ablate::finiteVolume {
 
class CellInterpolant {
   public:
    using DiscontinuousFluxFunction = PetscErrorCode (*)(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt uOff[], const PetscScalar fieldL[], const PetscScalar fieldR[], const PetscInt aOff[],
                                                         const PetscScalar auxL[], const PetscScalar auxR[], PetscScalar flux[], void* ctx);
 
    using PointFunction = PetscErrorCode (*)(PetscInt dim, PetscReal time, const PetscFVCellGeom* cg, const PetscInt uOff[], const PetscScalar u[], const PetscInt aOff[], const PetscScalar a[],
                                             PetscScalar f[], void* ctx);
 
    struct DiscontinuousFluxFunctionDescription {
        DiscontinuousFluxFunction function;
        void* context;
 
        PetscInt field;
        std::vector<PetscInt> inputFields;
        std::vector<PetscInt> auxFields;
    };
 
    struct PointFunctionDescription {
        PointFunction function;
        void* context;
 
        std::vector<PetscInt> fields;
        std::vector<PetscInt> inputFields;
        std::vector<PetscInt> auxFields;
    };
 
   private:
    std::shared_ptr<ablate::domain::SubDomain> subDomain;
 
    std::vector<DM> gradientCellDms;
 
    void ComputeFluxSourceTerms(DM dm, PetscDS ds, PetscInt totDim, const PetscScalar* xArray, DM dmAux, PetscDS dsAux, PetscInt totDimAux, const PetscScalar* auxArray, DM faceDM,
                                const PetscScalar* faceGeomArray, DM cellDM, const PetscScalar* cellGeomArray, std::vector<DM>& dmGrads, std::vector<const PetscScalar*>& locGradArrays,
                                PetscScalar* locFArray, const std::shared_ptr<domain::Region>& solverRegion, std::vector<CellInterpolant::DiscontinuousFluxFunctionDescription>& rhsFunctions,
                                const ablate::domain::Range& faceRange, const ablate::domain::Range& cellRange);
 
    void ProjectToFace(const std::vector<domain::Field>& fields, PetscDS ds, const PetscFVFaceGeom& faceGeom, PetscInt cellId, const PetscFVCellGeom& cellGeom, DM dm, const PetscScalar* xArray,
                       const std::vector<DM>& dmGrads, const std::vector<const PetscScalar*>& gradArrays, PetscScalar* u, PetscScalar* grad, bool projectField = true);
 
    void ComputeFieldGradients(const domain::Field& field, Vec xLocalVec, Vec& gradLocVec, DM& dmGrad, Vec cellGeomVec, Vec faceGeomVec, const ablate::domain::Range& faceRange,
                               const ablate::domain::Range& cellRange);
 
    static PetscErrorCode ComputeGradientFVM(DM dm, DMLabel regionLabel, PetscInt regionValue, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM* dmGrad);
 
   public:
    CellInterpolant(std::shared_ptr<ablate::domain::SubDomain> subDomain, const std::shared_ptr<domain::Region>& solverRegion, Vec faceGeomVec, Vec cellGeomVec);
    ~CellInterpolant();
 
    void ComputeRHS(PetscReal time, Vec locXVec, Vec locAuxVec, Vec locFVec, const std::shared_ptr<domain::Region>& solverRegion,
                    std::vector<CellInterpolant::DiscontinuousFluxFunctionDescription>& rhsFunctions, const ablate::domain::Range& faceRange, const ablate::domain::Range& cellRange, Vec cellGeomVec,
                    Vec faceGeomVec);
 
    void ComputeRHS(PetscReal time, Vec locXVec, Vec locAuxVec, Vec locFVec, const std::shared_ptr<domain::Region>& solverRegion, std::vector<CellInterpolant::PointFunctionDescription>& rhsFunctions,
                    const ablate::domain::Range& cellRange, Vec cellGeomVec);
};
 
}  // namespace ablate::finiteVolume
 
#endif  // ABLATELIBRARY_CELLINTERPOLANT_HPP