/* *********************************************************************** // // Copyright (C) 2025 -- The 4D-STAR Collaboration // File Author: Emily Boudreaux // Last Modified: February 14, 2025 // // 4DSSE is free software; you can use it and/or modify // it under the terms and restrictions the GNU General Library Public // License version 3 (GPLv3) as published by the Free Software Foundation. // // 4DSSE is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. // See the GNU Library General Public License for more details. // // You should have received a copy of the GNU Library General Public License // along with this software; if not, write to the Free Software // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // // *********************************************************************** */ #ifndef POLYMFEMUTILS_H #define POLYMFEMUTILS_H #include "mfem.hpp" #include #include #include #include "config.h" #include "probe.h" #include "quill/LogMacros.h" /** * @file polyMFEMUtils.h * @brief A collection of utilities for working with MFEM and solving the lane-emden equation. */ /** * @namespace polyMFEMUtils * @brief A namespace for utilities for working with MFEM and solving the lane-emden equation. */ namespace polyMFEMUtils { /** * @brief A class for nonlinear power integrator. */ class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator { private: Config& config = Config::getInstance(); mfem::Coefficient &coeff_; double polytropicIndex; public: /** * @brief Constructor for NonlinearPowerIntegrator. * * @param coeff The function coefficient. * @param n The polytropic index. */ NonlinearPowerIntegrator(mfem::Coefficient &coeff, double n); /** * @brief Assembles the element vector. * * @param el The finite element. * @param Trans The element transformation. * @param elfun The element function. * @param elvect The element vector to be assembled. */ virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; /** * @brief Assembles the element gradient. * * @param el The finite element. * @param Trans The element transformation. * @param elfun The element function. * @param elmat The element matrix to be assembled. */ virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; }; /** * @brief A wrapper class for bilinear integrator. */ class BilinearIntegratorWrapper : public mfem::NonlinearFormIntegrator { private: mfem::BilinearFormIntegrator *integrator; public: /** * @brief Constructor for BilinearIntegratorWrapper. * * @param integratorInput The bilinear form integrator input. */ BilinearIntegratorWrapper(mfem::BilinearFormIntegrator *integratorInput); /** * @brief Destructor for BilinearIntegratorWrapper. */ virtual ~BilinearIntegratorWrapper(); /** * @brief Assembles the element vector. * * @param el The finite element. * @param Trans The element transformation. * @param elfun The element function. * @param elvect The element vector to be assembled. */ virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; /** * @brief Assembles the element gradient. * * @param el The finite element. * @param Trans The element transformation. * @param elfun The element function. * @param elmat The element matrix to be assembled. */ virtual void AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; }; /** * @brief A class for composite nonlinear integrator. */ class CompositeNonlinearIntegrator: public mfem::NonlinearFormIntegrator { private: std::vector integrators; public: /** * @brief Constructor for CompositeNonlinearIntegrator. */ CompositeNonlinearIntegrator(); /** * @brief Destructor for CompositeNonlinearIntegrator. */ virtual ~CompositeNonlinearIntegrator(); /** * @brief Adds an integrator to the composite integrator. * * @param integrator The nonlinear form integrator to add. */ void add_integrator(mfem::NonlinearFormIntegrator *integrator); /** * @brief Assembles the element vector. * * @param el The finite element. * @param Trans The element transformation. * @param elfun The element function. * @param elvect The element vector to be assembled. */ virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; /** * @brief Assembles the element gradient. * * @param el The finite element. * @param Trans The element transformation. * @param elfun The element function. * @param elmat The element matrix to be assembled. */ virtual void AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; }; /** * @brief A class for constraint integrator. */ class ConstraintIntegrator: public mfem::NonlinearFormIntegrator { private: mfem::Coefficient η public: /** * @brief Constructor for ConstraintIntegrator. * * @param eta The coefficient. */ ConstraintIntegrator(mfem::Coefficient &eta_); /** * @brief Assembles the element vector. * * @param el The finite element. * @param Trans The element transformation. * @param elfun The element function. * @param elvect The element vector to be assembled. */ virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; /** * @brief Assembles the element gradient. * * @param el The finite element. * @param Trans The element transformation. * @param elfun The element function. * @param elmat The element matrix to be assembled. */ virtual void AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; }; class GaussianCoefficient : public mfem::Coefficient { private: double stdDev; double norm_coeff; public: GaussianCoefficient(double stdDev); virtual double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override; }; class AugmentedOperator : public mfem::Operator { private: Config& config = Config::getInstance(); mfem::NonlinearForm &nfl; mfem::LinearForm &C; double C_val; int lambdaDofOffset; mutable mfem::SparseMatrix *lastJacobian = nullptr; public: AugmentedOperator(mfem::NonlinearForm &nfl_, mfem::LinearForm &C_, int lambdaDofOffset_, double C_val_); ~AugmentedOperator(); virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const override; virtual mfem::Operator &GetGradient(const mfem::Vector &x) const override; }; /** * @brief Calculates the Gaussian integral. * * @param mesh The mesh. * @param gaussianCoeff The Gaussian coefficient. * @return The Gaussian integral. */ double calculateGaussianIntegral(mfem::Mesh &mesh, polyMFEMUtils::GaussianCoefficient &gaussianCoeff); class ZeroSlopeNewtonSolver : public mfem::NewtonSolver { private: Config& config = Config::getInstance(); Probe::LogManager& logManager = Probe::LogManager::getInstance(); quill::Logger* logger = logManager.getLogger("log"); double alpha; // The penalty term for the flat slope at zero std::vector zeroSlopeCoordinate; // The coordinate of the zero slope point int zeroSlopeElemID = -1; double zeroIPReferenceCoord[4] = {0.0, 0.0, 0.0, 1.0}; mfem::IntegrationPoint zeroIP; mfem::Array zeroSlopeConnectedElements; std::vector zeroSlopeIPs; std::vector> zeroSlopeDofs; std::unique_ptr u_gf; mutable mfem::SparseMatrix *grad = nullptr; void ComputeConstrainedResidual(const mfem::Vector &x, mfem::Vector &r) const; void ComputeConstrainedGradient(const mfem::Vector &x) const; public: ZeroSlopeNewtonSolver(double alpha_, std::vector zeroSlopeCoordinate_); ~ZeroSlopeNewtonSolver(); // virtual void ProcessNewState(const mfem::Vector &x) const; virtual void SetOperator(const mfem::Operator &op) override; void Mult(const mfem::Vector &b, mfem::Vector &x) const override; }; } // namespace polyMFEMUtils #endif // POLYMFEMUTILS_H