feat(poly): added first pass implimentation of 3D constrained lane-emden solver

This has not currently been tested and this commit should not be viewed as scientifically complete

src/poly/coeff/meson.build

@@ -1,22 +1,22 @@
 polyCoeff_sources = files(
-  'private/coeff.cpp'
+  'private/polyCoeff.cpp'
 )
 
 polyCoeff_headers = files(
-  'public/coeff.h'
+  'public/polyCoeff.h'
 )
 
 libPolyCoeff = static_library('polyCoeff',
   polyCoeff_sources,
-  include_directories : include_directories('.'),
+  include_directories : include_directories('./public'),
   cpp_args: ['-fvisibility=default'],
   dependencies: [mfem_dep],
   install: true
 )
 
 
-polyCoeff_dep = declare_dependency(
+polycoeff_dep = declare_dependency(
-  include_directories : include_directories('.'),
+  include_directories : include_directories('./public'),
   link_with : libPolyCoeff,
   sources : polyCoeff_sources,
   dependencies : [mfem_dep]

src/poly/coeff/private/polyCoeff.cpp

@@ -1,40 +1,23 @@
 #include "mfem.hpp"
 #include <cmath>
 
-#include "coeff.h"
+#include "polyCoeff.h"
 
-/**
+namespace polycoeff{
- * @brief Computes the xi coefficient function.
- * 
- * @param x Input vector.
- * @return double The computed xi coefficient.
- */
   double xi_coeff_func(const mfem::Vector &x)
   {
     return std::pow(x(0), 2);
   }
 
-/**
- * @brief Computes the vector xi coefficient function.
- * 
- * @param x Input vector.
- * @param v Output vector to store the computed xi coefficient.
- */
   void vec_xi_coeff_func(const mfem::Vector &x, mfem::Vector &v)
   {
     v.SetSize(1);
     v[0] = -std::pow(x(0), 2);
   }
 
-/**
- * @brief Computes the initial guess for theta.
- * 
- * @param x Input vector.
- * @param root Root value used in the computation.
- * @return double The initial guess for theta.
- */
   double theta_initial_guess(const mfem::Vector &x, double root)
   {
     double xi = x[0];
     return 1 - std::pow(xi / root, 2);
   }
+}

src/poly/coeff/public/polyCoeff.h

@@ -1,8 +1,35 @@
+#ifndef POLYCOEFF_H
+#define POLYCOEFF_H
+
 #include "mfem.hpp"
 #include <cmath>
 
+namespace polycoeff
+{
+  /**
+   * @brief Computes the xi coefficient function.
+   * 
+   * @param x Input vector.
+   * @return double The computed xi coefficient.
+   */
     double xi_coeff_func(const mfem::Vector &x);
 
+  /**
+   * @brief Computes the vector xi coefficient function.
+   * 
+   * @param x Input vector.
+   * @param v Output vector to store the computed xi coefficient.
+   */
     void vec_xi_coeff_func(const mfem::Vector &x, mfem::Vector &v);
 
+  /**
+   * @brief Computes the initial guess for theta.
+   * 
+   * @param x Input vector.
+   * @param root Root value used in the computation.
+   * @return double The initial guess for theta.
+   */
     double theta_initial_guess(const mfem::Vector &x, double root);
+} // namespace polyCoeff
+
+#endif // POLYCOEFF_H

src/poly/solver/meson.build

@@ -0,0 +1,22 @@
+polySolver_sources = files(
+  'private/polySolver.cpp'
+)
+
+polySolver_headers = files(
+  'public/polySolver.h'
+)
+
+libPolySolver = static_library('polySolver',
+  polySolver_sources,
+  include_directories : include_directories('./public'),
+  cpp_args: ['-fvisibility=default'],
+  dependencies: [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep],
+  install: true
+)
+
+polysolver_dep = declare_dependency(
+  include_directories : include_directories('./public'),
+  link_with : libPolySolver,
+  sources : polySolver_sources,
+  dependencies : [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep]
+)

src/poly/solver/private/polySolver.cpp

@@ -0,0 +1,144 @@
+#include "mfem.hpp"
+
+#include <string>
+#include <iostream>
+#include <memory>
+
+#include "meshIO.h"
+#include "polySolver.h"
+#include "polyMFEMUtils.h"
+#include "polyCoeff.h"
+
+
+// TODO: Come back to this and think of a better way to get the mesh file
+const std::string SPHERICAL_MESH = std::string(getenv("MESON_SOURCE_ROOT")) + "/src/resources/mesh/sphere.msh";
+
+PolySolver::PolySolver(double n, double order) 
+: n(n),
+  order(order),
+  meshIO(SPHERICAL_MESH),
+  mesh(meshIO.GetMesh()),
+  gaussianCoeff(std::make_unique<polyMFEMUtils::GaussianCoefficient>(0.1)),
+  diffusionCoeff(std::make_unique<mfem::VectorConstantCoefficient>(mfem::Vector(mesh.SpaceDimension()))),
+  nonLinearSourceCoeff(std::make_unique<mfem::ConstantCoefficient>(1.0))
+  {
+    (*diffusionCoeff).GetVec() = 1.0;
+    feCollection = std::make_unique<mfem::H1_FECollection>(order, mesh.SpaceDimension());
+
+    feSpace = std::make_unique<mfem::FiniteElementSpace>(&mesh, feCollection.get());
+    lambdaFeSpace = std::make_unique<mfem::FiniteElementSpace>(&mesh, feCollection.get(), 1); // Scalar space for lambda
+
+    compositeIntegrator = std::make_unique<polyMFEMUtils::CompositeNonlinearIntegrator>();
+    nonlinearForm = std::make_unique<mfem::NonlinearForm>(feSpace.get());
+
+    C = std::make_unique<mfem::LinearForm>(feSpace.get());
+
+    u = std::make_unique<mfem::GridFunction>(feSpace.get());
+
+    assembleNonlinearForm();
+    assembleConstraintForm();
+}
+
+PolySolver::assembleNonlinearForm() {
+    // Add the \int_{\Omega}\nabla v\cdot\nabla\theta d\Omegaterm
+    compositeIntegrator->add_integrator(
+        new polyMFEMUtils::BilinearIntegratorWrapper(
+            new mfem::DiffusionIntegrator(diffusionCoeff.get()),
+        )
+    );
+
+    // Add the \int_{\Omega}v\theta^{n} d\Omega term
+    compositeIntegrator->add_integrator(
+        new polyMFEMUtils::NonlinearPowerIntegrator(
+            nonLinearSourceCoeff.get(),
+            n
+        )
+    );
+
+    compositeIntegrator->add_integrator(
+        new polyMFEMUtils::ConstraintIntegrator(
+            *gaussianCoeff
+        )
+    );
+
+    nonlinearForm->AddDomainIntegrator(compositeIntegrator.get());
+}
+
+PolySolver::assembleConstraintForm() {
+    C->AddDomainIntegrator(
+        new mfem::DomainLFIntegrator(
+            *gaussianCoeff
+        )
+    );
+    C->Assemble();
+}
+
+PolySolver::solve(){
+    // --- Set the initial guess for the solution ---
+    mfem::FunctionCoefficient initCoeff (
+        [this](const mfem::Vector &x) {
+            return 1.0; // Update this to be a better init guess
+        }
+    );
+    u->ProjectCoefficient(initCoeff);
+
+    // --- Combine DOFs (u and λ) into a single vector ---
+    int lambdaDofOffset  = feSpace->GetTrueVSize(); // Get the size of θ space
+    int totalTrueDofs = lambdaDofOffset + lambdaFeSpace->GetTrueVSize();
+
+    if (totalTrueDofs != lambdaDofOffset + 1) {
+        throw std::runtime_error("The total number of true dofs is not equal to the sum of the lambda offset and the lambda space size");
+    }
+
+    mfem::Vector U(totalTrueDofs);
+    U = 0.0;
+
+    u->GetTrueDofs(U.GetBlock(0));
+
+    // --- Setup the Newton Solver ---
+    mfem::NewtonSolver newtonSolver;
+    newtonSolver.SetRelTol(1e-8);
+    newtonSolver.SetAbsTol(1e-10);
+    newtonSolver.SetMaxIter(200);
+    newtonSolver.SetPrintLevel(1);
+
+    // --- Setup the GMRES Solver ---
+    // --- GMRES is good for indefinite systems ---
+    mfem::GMRESSolver gmresSolver;
+    gmresSolver.SetRelTol(1e-10);
+    gmresSolver.SetAbsTol(1e-12);
+    gmresSolver.SetMaxIter(2000);
+    gmresSolver.SetPrintLevel(0);
+    newtonSolver.SetSolver(gmresSolver);
+    // TODO: Change numeric tolerance to grab from the tol module
+
+    // --- Setup the Augmented Operator ---
+    polyMFEMUtils::AugmentedOperator aug_op(nonlinearForm.get(), C.get(), lambdaDofOffset);
+    newtonSolver.SetOperator(aug_op);
+
+    // --- Create the RHS of the augmented system ---
+    mfem::Vector B(totalTrueDofs);
+    B = 0.0;
+
+    // Set the constraint value (∫η(r) dΩ) in the last entry of B
+    // This sets the last entry to 1.0, this mighht be a problem later on...
+    mfem::ConstantCoefficient one(1.0);
+    mfem::LinearForm constraint_rhs(lambdaFeSpace.get());
+    constraint_rhs.AddDomainIntegrator(
+        new mfem::DomainLFIntegrator(*gaussianCoeff)
+    );
+
+    constraint_rhs.Assemble();
+    B[lambdaDofOffset] = constraint_rhs(0); // Get that single value for the rhs. Only one value because it's a scalar space
+
+
+    // --- Solve the augmented system ---
+    newtonSolver.Mult(B, U);
+
+    // --- Extract the Solution ---
+    u->Distribute(U.GetBlock(0));
+    double lambda = U[lambdaDofOffset];
+
+    std::cout << "λ = " << lambda << std::endl;
+    // TODO : Add a way to get the solution out of the solver
+}

src/poly/solver/public/polySolver.h

@@ -0,0 +1,46 @@
+#ifndef POLYSOLVER_H
+#define POLYSOLVER_H
+
+#include "mfem.hpp"
+#include <iostream>
+#include <string>
+#include <memory>
+
+#include "meshIO.h"
+#include "polyCoeff.h"
+
+
+class PolySolver {
+private:
+    double n, order;
+    MeshIO meshIO;
+    mfem::Mesh& mesh;
+    
+    std::unique_ptr<mfem::H1_FECollection> feCollection;
+
+    std::unique_ptr<mfem::FiniteElementSpace> feSpace;
+    std::unique_ptr<mfem::FiniteElementSpace> lambdaFeSpace;
+
+    std::unique_ptr<polyMFEMUtils::CompositeNonlinearIntegrator> compositeIntegrator;
+
+    std::unique_ptr<mfem::NonlinearForm> nonlinearForm;
+    std::unique_ptr<mfem::LinearForm> C; // For the constraint equation
+
+    std::unique_ptr<mfem::GridFunction> u;
+
+    std::unique_ptr<mfem::Vector> oneVec;
+
+    std::unique_ptr<mfem::VectorConstantCoefficient> diffusionCoeff;
+    std::unique_ptr<mfem::ConstantCoefficient> nonLinearSourceCoeff;
+    std::unique_ptr<polyMFEMUtils::GaussianCoefficient> gaussianCoeff;
+
+    void assembleNonlinearForm();
+    void assembleConstraintForm();
+
+
+public:
+    PolySolver(double n, double order);
+    void solve();
+};
+
+#endif // POLYSOLVER_H

src/poly/utils/meson.build

@@ -10,14 +10,14 @@ polyutils_headers = files(
 
 libpolyutils = static_library('polyutils',
   polyutils_sources,
-  include_directories : include_directories('.'),
+  include_directories : include_directories('./public'),
   cpp_args: ['-fvisibility=default'],
   dependencies: [mfem_dep],
   install: true
 )
 
-libpolyutils_dep = declare_dependency(
+polyutils_dep = declare_dependency(
-  include_directories : include_directories('.'),
+  include_directories : include_directories('./public'),
   link_with : libpolyutils,
   sources : polyutils_sources,
   dependencies : [mfem_dep]

src/poly/utils/private/polyMFEMUtils.cpp

@@ -2,11 +2,13 @@
 #include <string>
 #include <iostream>
 #include <cmath>
+#include <numbers>
 
 #include "polyMFEMUtils.h"
 
+namespace polyMFEMUtils {
   NonlinearPowerIntegrator::NonlinearPowerIntegrator(
- mfem::FunctionCoefficient &coeff,
+  mfem::Coefficient &coeff,
   double n) : coeff_(coeff), polytropicIndex(n) { 
 
   }
@@ -60,7 +62,7 @@ void NonlinearPowerIntegrator::AssembleElementGrad (
     mfem::Vector shape(dof);
 
     for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
-    mfem::IntegrationPoint ip = ir->IntPoint(iqp);
+      const mfem::IntegrationPoint &ip = ir->IntPoint(iqp);
       Trans.SetIntPoint(&ip);
       double weight = ip.weight * Trans.Weight();
 
@@ -173,3 +175,129 @@ void CompositeNonlinearIntegrator::AssembleElementGrad(
       elmat.Add(1.0, temp);
     }
   }
+
+  ConstraintIntegrator::ConstraintIntegrator(mfem::Coefficient &eta_) : eta(eta_) {}
+
+  void ConstraintIntegrator::AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) {
+    int dof = el.GetDof();
+    elvect.SetSize(dof);
+    elvect = 0.0;
+
+    mfem::Vector shape(dof);
+    const int intOrder = 2 * el.GetOrder() + 3;
+
+    mfem::IntegrationRule ir(el.GetGeomType(), intOrder);
+
+    for (int i = 0; i < ir.GetNPoints(); i++) {
+      const mfem::IntegrationPoint &ip = ir.IntPoint(i);
+      Trans.SetIntPoint(&ip);
+      el.CalcShape(ip, shape);
+
+      double eta_val = eta.Eval(Trans, ip);
+
+      double weight = ip.weight * Trans.Weight();
+      elvect.Add(eta_val * weight, shape);
+    }
+  }
+
+  void ConstraintIntegrator::AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) {
+    int dof = el.GetDof();
+    elmat.SetSize(dof);
+    elmat = 0.0;
+  }
+
+  GaussianCoefficient::GaussianCoefficient(double stdDev_)
+   : stdDev(stdDev_),
+     norm_coeff(1.0/(std::pow(std::sqrt(2*std::numbers::pi)*stdDev_,3))) {}
+
+  GaussianCoefficient::Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) {
+    mfem::Vector r;
+    T.Transform(ip, r);
+    double r2 = r * r;
+    return norm_coeff * std::exp(-r2/(2*std::pow(stdDev, 2)));
+  }
+
+
+  AugmentedOperator::AugmentedOperator(mfem::NonlinearForm &nfl_, mfem::LinearForm &C_, int lambdaDofOffset_)
+    : 
+    mfem::Operator( // Call the base class constructor
+      nfl_->FESpace()->GetTrueVSize()+1, // Sets the height attribute (+1 for the lambda component)
+      nfl_->FESpace()->GetTrueVSize()+1 // Sets the width attribute (+1 for the lambda component)
+    ),
+    nfl(&nfl_),
+    C(&C_), 
+    lambdaDofOffset(lambdaDofOffset_),
+    lastJacobian(nullptr) {}
+
+  void AugmentedOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
+    // Select the lambda component of the input vector and seperate it from the θ component
+    mfem::Vector u(x.GetData(), lambdaDofOffset);
+    double lambda = x[lambdaDofOffset];
+
+    // Compute the residual of the nonlinear form (F(u) - λC)
+    mfem::Vector F(lambdaDofOffset);
+    nfl->Mult(u, F); // This now includes the -λ∫vη(r) dΩ term
+
+    // Compute the transpose of C for the off diagonal terms of the augmented operator
+    y.SetSize(height);
+    y = 0.0;
+
+    y.SetBlock(0, F);
+    y[lambdaDofOffset] = (*C)(u); // Cᵀ×u. This is equivalent to ∫ η(r)θ dΩ
+
+    // add -lambda * C to the residual
+    mfem::Vector lambda_C(lambdaDofOffset);
+    C->GetVector(lambda_C) // Get the coefficient vector for C. I.e. ∫ vη(r) dΩ
+    lambda_C *= -lambda; // Multiply by -λ (this is now the term −λ ∫ vη(r)dΩ)
+
+    y.AddBlockVector(0, lambda_C); // Add the term to the residual
+  }
+
+  mfem::Operator &AugmentedOperator::GetGradient(const mfem::Vector &x) const {
+    // Select the lambda component of the input vector and seperate it from the θ component
+    mfem::Vector u(x.GetData(), lambdaDofOffset);
+    double lambda = x[lambdaDofOffset];
+
+    // --- Create the augmented Jacobian matrix ---
+    mfem::Operator *Jnfl = nfl->GetGradient(u); // Get the gradient of the nonlinear form
+
+    mfem::SparseMatrix *J_aug = new mfem::SparseMatrix(height, width);
+
+    // Fill in the blocks of the augmented Jacobian matrix
+    // Top-Left: Jacobian of the nonlinear form
+    mfem::SparseMatrix *Jnfl_sparse = dynamic_cast<mfem::SparseMatrix*>(Jnfl);
+
+    // Copy the original Jacobian into the augmented Jacobian
+    if (Jnfl_sparse) { // Checks if the dynamic cast was successful
+      for (int i = 0; i < Jnfl_sparse->Height(); i++) {
+        for (int j = 0; j < Jnfl_sparse->Width(); j++) {
+          J_aug->Set(i, j, Jnfl_sparse->Get(i, j));
+        }
+      }
+    } else {
+      MFEM_ABORT("GetGradient did not return a SparseMatrix");
+    }
+
+    // Bottom-left C (the constraint matrix)
+    mfem::Vector CVec(lambdaDofOffset);
+    C->GetVector(CVec);
+    for (int i = 0; i < CVec.Size(); i++) {
+      J_aug->Set(lambdaDofOffset, i, CVec[i]);
+    }
+
+    // Top-right -Cᵀ (the negative transpose of the constraint matrix)
+    for (int i = 0; i < CVec.Size(); i++) {
+      J_aug->Set(i, lambdaDofOffset, -CVec[i]);
+    }
+
+    J_aug->Finalize();
+
+    delete lastJacobian;
+    lastJacobian = J_aug;
+    return *lastJacobian;
+  }
+
+  AugmentedOperator::~AugAugmentedOperator() {
+    delete lastJacobian;
+  }
+} // namespace polyMFEMUtils

src/poly/utils/public/polyMFEMUtils.h

@@ -2,14 +2,23 @@
 #include <string>
 
 
-void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename);
+/**
+ * @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: 
-  mfem::FunctionCoefficient coeff_;
+    mfem::Coefficient &coeff_;
     double polytropicIndex;
   public:
     /**
@@ -18,7 +27,7 @@ public:
      * @param coeff The function coefficient.
      * @param n The polytropic index.
      */
-  NonlinearPowerIntegrator(mfem::FunctionCoefficient &coeff, double n);
+    NonlinearPowerIntegrator(mfem::Coefficient &coeff, double n);
 
     /**
      * @brief Assembles the element vector.
@@ -126,3 +135,65 @@ class CompositeNonlinearIntegrator: public mfem::NonlinearFormIntegrator {
        */
       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 &eta;
+
+    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:
+      std::unique_ptr<mfem::NonlinearForm> nfl;
+      std::unique_ptr<mfem::LinearForm> C;
+      int lambdaDofOffset;
+      mfem::SparseMatrix *lastJacobian = nullptr;
+
+    public:
+      AugmentedOperator(mfem::NonlinearForm &nfl_, mfem::LinearForm &C_, int lambdaDofOffset_);
+      ~AugmentedOperator();
+
+      virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+      virtual mfem::Operator &GetGradient(const mfem::Vector &x) const override;
+  };
+} // namespace polyMFEMUtils