feat(poly): started work on penalty term in variational form

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)

src/poly/solver/private/polySolver.cpp

@@ -21,6 +21,7 @@
 #include "mfem.hpp"
 
 #include <string>
+#include <stdexcept>
 
 #include "polySolver.h"
 #include "polyMFEMUtils.h"
@@ -75,6 +76,12 @@ PolySolver::PolySolver(double n, double order, mfem::Mesh& mesh_)
   nonlinearForm(std::make_unique<mfem::NonlinearForm>(feSpace.get())),
   u(std::make_unique<mfem::GridFunction>(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<mfem::VectorFunctionCoefficient>(mesh.SpaceDimension(), polycoeff::diffusionCoeff);
     nonlinearSourceCoeff = std::make_unique<mfem::FunctionCoefficient>(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<bool>("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;

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

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<int> m_originElementIDs;
+    mfem::Array<mfem::IntegrationPoint> 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<double> 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<int> zeroSlopeConnectedElements;
-      std::vector<mfem::IntegrationPoint> zeroSlopeIPs;
-      std::vector<mfem::Array<int>> zeroSlopeDofs;
-
-      std::unique_ptr<mfem::GridFunction> 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<double> 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