feat(poly): lagrangian constrained weak form of 3D lane-Emden

added a basic implimentation of the 3D lane emden equation using a lagrangian multiplier to constrain the value at the center of a spherical domain
2025-02-20 15:28:00 -05:00
parent deab5be0c1
commit 1fd1e624f2
4 changed files with 89 additions and 85 deletions
--- a/src/poly/solver/private/polySolver.cpp
+++ b/src/poly/solver/private/polySolver.cpp
@@ -13,67 +13,56 @@
 // 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) 
+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());
-
+  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())),
+  diffusionCoeff(std::make_unique<mfem::VectorConstantCoefficient>([&](){
+        mfem::Vector diffusionCoeffVec(mesh.SpaceDimension());
+        diffusionCoeffVec = 1.0;
+        return diffusionCoeffVec;
+    }())),
+  nonLinearSourceCoeff(std::make_unique<mfem::ConstantCoefficient>(1.0)),
+  gaussianCoeff(std::make_unique<polyMFEMUtils::GaussianCoefficient>(0.1)) {
    assembleNonlinearForm();
    assembleConstraintForm();
 }

-PolySolver::assembleNonlinearForm() {
+PolySolver::~PolySolver() {}
+
+void PolySolver::assembleNonlinearForm() {
    // Add the \int_{\Omega}\nabla v\cdot\nabla\theta d\Omegaterm
-    compositeIntegrator->add_integrator(
-        new polyMFEMUtils::BilinearIntegratorWrapper(
-            new mfem::DiffusionIntegrator(diffusionCoeff.get()),
-        )
+    auto wrappedDiffusionIntegrator = std::make_unique<polyMFEMUtils::BilinearIntegratorWrapper>(
+        new mfem::DiffusionIntegrator(*diffusionCoeff)
    );
+    compositeIntegrator->add_integrator(wrappedDiffusionIntegrator.release());

    // Add the \int_{\Omega}v\theta^{n} d\Omega term
-    compositeIntegrator->add_integrator(
-        new polyMFEMUtils::NonlinearPowerIntegrator(
-            nonLinearSourceCoeff.get(),
-            n
-        )
-    );
+    auto nonLinearIntegrator = std::make_unique<polyMFEMUtils::NonlinearPowerIntegrator>(*nonLinearSourceCoeff, n);
+    compositeIntegrator->add_integrator(nonLinearIntegrator.release());

-    compositeIntegrator->add_integrator(
-        new polyMFEMUtils::ConstraintIntegrator(
-            *gaussianCoeff
-        )
-    );
+    // Add the \int_{\Omega}v\eta(r) d\Omega term
+    auto constraintIntegrator = std::make_unique<polyMFEMUtils::ConstraintIntegrator>(*gaussianCoeff);
+    compositeIntegrator->add_integrator(constraintIntegrator.release());

-    nonlinearForm->AddDomainIntegrator(compositeIntegrator.get());
+    nonlinearForm->AddDomainIntegrator(compositeIntegrator.release());
 }

-PolySolver::assembleConstraintForm() {
-    C->AddDomainIntegrator(
-        new mfem::DomainLFIntegrator(
-            *gaussianCoeff
-        )
-    );
+void PolySolver::assembleConstraintForm() {
+    auto constraintIntegrator = std::make_unique<mfem::DomainLFIntegrator>(*gaussianCoeff);
+    C->AddDomainIntegrator(constraintIntegrator.release());
    C->Assemble();
 }

-PolySolver::solve(){
+void PolySolver::solve(){
    // --- Set the initial guess for the solution ---
    mfem::FunctionCoefficient initCoeff (
        [this](const mfem::Vector &x) {
@@ -86,6 +75,10 @@ PolySolver::solve(){
    int lambdaDofOffset  = feSpace->GetTrueVSize(); // Get the size of θ space
    int totalTrueDofs = lambdaDofOffset + lambdaFeSpace->GetTrueVSize();

+    std::cout << "Total True Dofs: " << totalTrueDofs << std::endl;
+    std::cout << "Lambda Dof Offset: " << lambdaDofOffset << std::endl;
+    std::cout << "Lambda True Dofs: " << lambdaFeSpace->GetTrueVSize() << std::endl;
+    std::cout << "U True Dofs: " << feSpace->GetTrueVSize() << std::endl;
    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");
    }
@@ -93,7 +86,8 @@ PolySolver::solve(){
    mfem::Vector U(totalTrueDofs);
    U = 0.0;

-    u->GetTrueDofs(U.GetBlock(0));
+    mfem::Vector u_view(U.GetData(), lambdaDofOffset);
+    u->GetTrueDofs(u_view);

    // --- Setup the Newton Solver ---
    mfem::NewtonSolver newtonSolver;
@@ -113,7 +107,7 @@ PolySolver::solve(){
    // TODO: Change numeric tolerance to grab from the tol module

    // --- Setup the Augmented Operator ---
-    polyMFEMUtils::AugmentedOperator aug_op(nonlinearForm.get(), C.get(), lambdaDofOffset);
+    polyMFEMUtils::AugmentedOperator aug_op(*nonlinearForm, *C, lambdaDofOffset);
    newtonSolver.SetOperator(aug_op);

    // --- Create the RHS of the augmented system ---
@@ -136,7 +130,8 @@ PolySolver::solve(){
    newtonSolver.Mult(B, U);

    // --- Extract the Solution ---
-    u->Distribute(U.GetBlock(0));
+    mfem::Vector u_sol_view(U.GetData(), lambdaDofOffset);
+    u->SetData(u_sol_view);
    double lambda = U[lambdaDofOffset];

    std::cout << "λ = " << lambda << std::endl;
--- a/src/poly/solver/public/polySolver.h
+++ b/src/poly/solver/public/polySolver.h
@@ -8,6 +8,7 @@

 #include "meshIO.h"
 #include "polyCoeff.h"
+#include "polyMFEMUtils.h"


 class PolySolver {
@@ -28,8 +29,6 @@ private:

    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;
@@ -40,6 +39,7 @@ private:

 public:
    PolySolver(double n, double order);
+    ~PolySolver();
    void solve();
 };

--- a/src/poly/utils/private/polyMFEMUtils.cpp
+++ b/src/poly/utils/private/polyMFEMUtils.cpp
@@ -132,11 +132,7 @@ namespace polyMFEMUtils {
  CompositeNonlinearIntegrator::CompositeNonlinearIntegrator() { }


-  CompositeNonlinearIntegrator::~CompositeNonlinearIntegrator() {
-    for (size_t i = 0; i < integrators.size(); i++) {
-      delete integrators[i];
-    }
-  }
+  CompositeNonlinearIntegrator::~CompositeNonlinearIntegrator() { }

  void CompositeNonlinearIntegrator::add_integrator(mfem::NonlinearFormIntegrator *integrator) {
    integrators.push_back(integrator);
@@ -186,7 +182,7 @@ namespace polyMFEMUtils {
    mfem::Vector shape(dof);
    const int intOrder = 2 * el.GetOrder() + 3;

-    mfem::IntegrationRule ir(el.GetGeomType(), intOrder);
+    const mfem::IntegrationRule &ir = mfem::IntRules.Get(el.GetGeomType(), intOrder);

    for (int i = 0; i < ir.GetNPoints(); i++) {
      const mfem::IntegrationPoint &ip = ir.IntPoint(i);
@@ -210,7 +206,7 @@ namespace polyMFEMUtils {
   : 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) {
+  double GaussianCoefficient::Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) {
    mfem::Vector r;
    T.Transform(ip, r);
    double r2 = r * r;
@@ -221,11 +217,11 @@ namespace polyMFEMUtils {
  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_.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_), 
+    nfl(nfl_),
+    C(C_), 
    lambdaDofOffset(lambdaDofOffset_),
    lastJacobian(nullptr) {}

@@ -236,51 +232,59 @@ namespace polyMFEMUtils {

    // 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
+    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Ω
+    for (int i = 0; i < lambdaDofOffset; i++) {
+      y[i] = F[i];
+    }
+    mfem::GridFunction u_gf(C.FESpace());
+    u_gf.SetData(u);
+    y[lambdaDofOffset] = C.operator()(u_gf);

    // 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Ω)
+    mfem::GridFunction constraint_gf(C.FESpace());
+    constraint_gf = 0.0;

-    y.AddBlockVector(0, lambda_C); // Add the term to the residual
+    mfem::Vector constraint_tdofs;
+    C.Assemble();
+    lambda_C = C.GetData();
+    lambda_C *= -lambda;  // Multiply by -λ (this is now the term −λ ∫ vη(r)dΩ)
+
+    for (int i = 0; i < lambdaDofOffset; i++) {
+      y[i] += lambda_C[i];
+    }
  }

  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::SparseMatrix *Jnfl_sparse = dynamic_cast<mfem::SparseMatrix*>(&nfl.GetGradient(u));
+    if (!Jnfl_sparse) {
      MFEM_ABORT("GetGradient did not return a SparseMatrix");
    }

+    mfem::SparseMatrix *J_aug = new mfem::SparseMatrix(height, width);
+    // Copy the original Jacobian into the augmented Jacobian
+    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->Elem(i, j));
+      }
+    }
+
    // Bottom-left C (the constraint matrix)
    mfem::Vector CVec(lambdaDofOffset);
-    C->GetVector(CVec);
+    mfem::GridFunction tempCGrid(C.FESpace());
+    C.Assemble();
+    CVec = C.GetData();
    for (int i = 0; i < CVec.Size(); i++) {
      J_aug->Set(lambdaDofOffset, i, CVec[i]);
    }
@@ -297,7 +301,7 @@ namespace polyMFEMUtils {
    return *lastJacobian;
  }

-  AugmentedOperator::~AugAugmentedOperator() {
+  AugmentedOperator::~AugmentedOperator() {
    delete lastJacobian;
  }
 } // namespace polyMFEMUtils
--- a/src/poly/utils/public/polyMFEMUtils.h
+++ b/src/poly/utils/public/polyMFEMUtils.h
@@ -1,3 +1,6 @@
+#ifndef POLYMFEMUTILS_H
+#define POLYMFEMUTILS_H
+
 #include "mfem.hpp"
 #include <string>

@@ -183,10 +186,10 @@ namespace polyMFEMUtils {

  class AugmentedOperator : public mfem::Operator {
    private:
-      std::unique_ptr<mfem::NonlinearForm> nfl;
-      std::unique_ptr<mfem::LinearForm> C;
+      mfem::NonlinearForm &nfl;
+      mfem::LinearForm &C;
      int lambdaDofOffset;
-      mfem::SparseMatrix *lastJacobian = nullptr;
+      mutable mfem::SparseMatrix *lastJacobian = nullptr;

    public:
      AugmentedOperator(mfem::NonlinearForm &nfl_, mfem::LinearForm &C_, int lambdaDofOffset_);
@@ -196,4 +199,6 @@ namespace polyMFEMUtils {

      virtual mfem::Operator &GetGradient(const mfem::Vector &x) const override;
  };
-} // namespace polyMFEMUtils
+} // namespace polyMFEMUtils
+
+#endif // POLYMFEMUTILS_H