From 1fd1e624f205e2a23cc39972968f92808f49e7c5 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 20 Feb 2025 15:28:00 -0500 Subject: [PATCH] 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 --- src/poly/solver/private/polySolver.cpp | 85 +++++++++++------------- src/poly/solver/public/polySolver.h | 4 +- src/poly/utils/private/polyMFEMUtils.cpp | 72 ++++++++++---------- src/poly/utils/public/polyMFEMUtils.h | 13 ++-- 4 files changed, 89 insertions(+), 85 deletions(-) diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 86d0a25..84365f4 100644 --- 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(0.1)), - diffusionCoeff(std::make_unique(mfem::Vector(mesh.SpaceDimension()))), - nonLinearSourceCoeff(std::make_unique(1.0)) - { - (*diffusionCoeff).GetVec() = 1.0; - feCollection = std::make_unique(order, mesh.SpaceDimension()); - - feSpace = std::make_unique(&mesh, feCollection.get()); - lambdaFeSpace = std::make_unique(&mesh, feCollection.get(), 1); // Scalar space for lambda - - compositeIntegrator = std::make_unique(); - nonlinearForm = std::make_unique(feSpace.get()); - - C = std::make_unique(feSpace.get()); - - u = std::make_unique(feSpace.get()); - + feCollection(std::make_unique(order, mesh.SpaceDimension())), + feSpace(std::make_unique(&mesh, feCollection.get())), + lambdaFeSpace(std::make_unique(&mesh, feCollection.get(), 1)), // Scalar space for lambda + compositeIntegrator(std::make_unique()), + nonlinearForm(std::make_unique(feSpace.get())), + C(std::make_unique(feSpace.get())), + u(std::make_unique(feSpace.get())), + diffusionCoeff(std::make_unique([&](){ + mfem::Vector diffusionCoeffVec(mesh.SpaceDimension()); + diffusionCoeffVec = 1.0; + return diffusionCoeffVec; + }())), + nonLinearSourceCoeff(std::make_unique(1.0)), + gaussianCoeff(std::make_unique(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( + 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(*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(*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(*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; diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 5ec4d0e..3ef60c3 100644 --- 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 u; - std::unique_ptr oneVec; - std::unique_ptr diffusionCoeff; std::unique_ptr nonLinearSourceCoeff; std::unique_ptr gaussianCoeff; @@ -40,6 +39,7 @@ private: public: PolySolver(double n, double order); + ~PolySolver(); void solve(); }; diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp index 4158dbf..6a3b569 100644 --- 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(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(&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 \ No newline at end of file diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h index 3b7097f..8bddb01 100644 --- 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 @@ -183,10 +186,10 @@ namespace polyMFEMUtils { class AugmentedOperator : public mfem::Operator { private: - std::unique_ptr nfl; - std::unique_ptr 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 \ No newline at end of file +} // namespace polyMFEMUtils + +#endif // POLYMFEMUTILS_H \ No newline at end of file