From b98f6b6ebd8a833b8566a6912352ee6b546e8843 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 19 Mar 2025 10:09:37 -0400 Subject: [PATCH] feat(poly): started work on penalty term in variational form --- src/poly/coeff/private/polyCoeff.cpp | 12 +-- src/poly/solver/private/polySolver.cpp | 16 +++- src/poly/utils/private/polyMFEMUtils.cpp | 13 +++ src/poly/utils/public/polyMFEMUtils.h | 110 +++-------------------- 4 files changed, 47 insertions(+), 104 deletions(-) diff --git a/src/poly/coeff/private/polyCoeff.cpp b/src/poly/coeff/private/polyCoeff.cpp index 43d0e94..de50ffd 100644 --- a/src/poly/coeff/private/polyCoeff.cpp +++ b/src/poly/coeff/private/polyCoeff.cpp @@ -25,17 +25,17 @@ namespace polycoeff{ double nonlinearSourceCoeff(const mfem::Vector &x) { - double r = x.Norml2(); - return std::pow(r, 2); + return 1; + // double r = x(0)*x(0) + x(1)*x(1) + x(2)*x(2); + // return std::pow(r, 2); } void diffusionCoeff(const mfem::Vector &x, mfem::Vector &v) { v.SetSize(3); - double r = x.Norml2(); - for (int i = 0; i < 3; i++) { - v(i) = -std::pow(r, 2); - } + v = -1; + // double r = x(0)*x(0) + x(1)*x(1) + x(2)*x(2); + // v = -std::pow(r, 2); } double x1(const double n) diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 95d688b..75446ef 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -21,6 +21,7 @@ #include "mfem.hpp" #include +#include #include "polySolver.h" #include "polyMFEMUtils.h" @@ -75,6 +76,12 @@ PolySolver::PolySolver(double n, double order, mfem::Mesh& mesh_) nonlinearForm(std::make_unique(feSpace.get())), u(std::make_unique(feSpace.get())) { + // --- Check the polytropic index --- + if (n > 4.99 || n < 0.0) { + LOG_ERROR(logger, "The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is {}", n); + throw std::runtime_error("The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is " + std::to_string(n)); + } + diffusionCoeff = std::make_unique(mesh.SpaceDimension(), polycoeff::diffusionCoeff); nonlinearSourceCoeff = std::make_unique(polycoeff::nonlinearSourceCoeff); @@ -103,11 +110,18 @@ void PolySolver::solve(){ mfem::FunctionCoefficient initCoeff ( [this](const mfem::Vector &x) { double r = x.Norml2(); - double theta = laneEmden::thetaSerieseExpansion(r, n, 10); + double theta = laneEmden::thetaSerieseExpansion(r, n, 12); return theta; + // double radius = Probe::getMeshRadius(mesh); + // double u = 1/radius; + + // return -std::pow((u*r), 2)+1.0; } ); u->ProjectCoefficient(initCoeff); + if (config.get("Poly:Solver:ViewInitialGuess", false)) { + Probe::glVisView(*u, mesh, "initialGuess"); + } // mfem::DenseMatrix centerPoint(mesh.SpaceDimension(), 7); mfem::DenseMatrix centerPoint(mesh.SpaceDimension(), 1); centerPoint(0, 0) = 0.0; diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp index d633139..5314beb 100644 --- a/src/poly/utils/private/polyMFEMUtils.cpp +++ b/src/poly/utils/private/polyMFEMUtils.cpp @@ -191,4 +191,17 @@ namespace polyMFEMUtils { } } + ConstraintIntegrator::ConstraintIntegrator(const double gamma): m_gamma(gamma) { + LOG_INFO(m_logger, "Initializing Constraint Integrator..."); + m_originCoordinateMatrix.SetSize(3, 1); + m_originCoordinateMatrix(0, 0) = 0.0; + m_originCoordinateMatrix(1, 0) = 0.0; + m_originCoordinateMatrix(2, 0) = 0.0; + } + + void ConstraintIntegrator::AssembleElementMatrix(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, mfem::DenseMatrix &elmat) { + + } + + } // namespace polyMFEMUtils diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h index 036854f..b38a342 100644 --- a/src/poly/utils/public/polyMFEMUtils.h +++ b/src/poly/utils/public/polyMFEMUtils.h @@ -166,107 +166,23 @@ namespace polyMFEMUtils { 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 η + class ConstraintIntegrator: public mfem::BilinearFormIntegrator { + private: + Config& m_config = Config::getInstance(); + Probe::LogManager& logManager = Probe::LogManager::getInstance(); + quill::Logger* m_logger = logManager.getLogger("log"); + const double m_gamma; + mfem::Array m_originElementIDs; + mfem::Array m_originIntegrationPoints; + mfem::DenseMatrix m_originCoordinateMatrix; + public: + ConstraintIntegrator(const double gamma); + ~ConstraintIntegrator() = default; - public: - /** - * @brief Constructor for ConstraintIntegrator. - * - * @param eta The coefficient. - */ - ConstraintIntegrator(mfem::Coefficient &eta_); + void AssembleElementMatrix(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, mfem::DenseMatrix &elmat) override; - /** - * @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 \ No newline at end of file