From b939fd68faa9b9be7b9d433185e8bc68610b0f29 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 19 Feb 2025 14:35:15 -0500 Subject: [PATCH] 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 | 10 +- src/poly/coeff/private/polyCoeff.cpp | 49 +-- src/poly/coeff/public/polyCoeff.h | 33 +- src/poly/solver/meson.build | 22 ++ src/poly/solver/private/polySolver.cpp | 144 ++++++++ src/poly/solver/public/polySolver.h | 46 +++ src/poly/utils/meson.build | 6 +- src/poly/utils/private/polyMFEMUtils.cpp | 422 +++++++++++++++-------- src/poly/utils/public/polyMFEMUtils.h | 261 +++++++++----- 9 files changed, 707 insertions(+), 286 deletions(-) create mode 100644 src/poly/solver/private/polySolver.cpp create mode 100644 src/poly/solver/public/polySolver.h diff --git a/src/poly/coeff/meson.build b/src/poly/coeff/meson.build index 112a1ff..4c82cff 100644 --- a/src/poly/coeff/meson.build +++ b/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( - include_directories : include_directories('.'), +polycoeff_dep = declare_dependency( + include_directories : include_directories('./public'), link_with : libPolyCoeff, sources : polyCoeff_sources, dependencies : [mfem_dep] diff --git a/src/poly/coeff/private/polyCoeff.cpp b/src/poly/coeff/private/polyCoeff.cpp index c4076a0..870e817 100644 --- a/src/poly/coeff/private/polyCoeff.cpp +++ b/src/poly/coeff/private/polyCoeff.cpp @@ -1,40 +1,23 @@ #include "mfem.hpp" #include -#include "coeff.h" +#include "polyCoeff.h" -/** - * @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); -} +namespace polycoeff{ + 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); -} + 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); + double theta_initial_guess(const mfem::Vector &x, double root) + { + double xi = x[0]; + return 1 - std::pow(xi / root, 2); + } } \ No newline at end of file diff --git a/src/poly/coeff/public/polyCoeff.h b/src/poly/coeff/public/polyCoeff.h index 46204f1..4a1e81d 100644 --- a/src/poly/coeff/public/polyCoeff.h +++ b/src/poly/coeff/public/polyCoeff.h @@ -1,8 +1,35 @@ +#ifndef POLYCOEFF_H +#define POLYCOEFF_H + #include "mfem.hpp" #include -double xi_coeff_func(const mfem::Vector &x); +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); -void vec_xi_coeff_func(const mfem::Vector &x, mfem::Vector &v); + /** + * @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); -double theta_initial_guess(const mfem::Vector &x, double root); \ No newline at end of file + /** + * @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 \ No newline at end of file diff --git a/src/poly/solver/meson.build b/src/poly/solver/meson.build index e69de29..de06e2c 100644 --- a/src/poly/solver/meson.build +++ b/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] +) \ No newline at end of file diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp new file mode 100644 index 0000000..86d0a25 --- /dev/null +++ b/src/poly/solver/private/polySolver.cpp @@ -0,0 +1,144 @@ +#include "mfem.hpp" + +#include +#include +#include + +#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(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()); + + 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 +} \ No newline at end of file diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h new file mode 100644 index 0000000..5ec4d0e --- /dev/null +++ b/src/poly/solver/public/polySolver.h @@ -0,0 +1,46 @@ +#ifndef POLYSOLVER_H +#define POLYSOLVER_H + +#include "mfem.hpp" +#include +#include +#include + +#include "meshIO.h" +#include "polyCoeff.h" + + +class PolySolver { +private: + double n, order; + MeshIO meshIO; + mfem::Mesh& mesh; + + std::unique_ptr feCollection; + + std::unique_ptr feSpace; + std::unique_ptr lambdaFeSpace; + + std::unique_ptr compositeIntegrator; + + std::unique_ptr nonlinearForm; + std::unique_ptr C; // For the constraint equation + + std::unique_ptr u; + + std::unique_ptr oneVec; + + std::unique_ptr diffusionCoeff; + std::unique_ptr nonLinearSourceCoeff; + std::unique_ptr gaussianCoeff; + + void assembleNonlinearForm(); + void assembleConstraintForm(); + + +public: + PolySolver(double n, double order); + void solve(); +}; + +#endif // POLYSOLVER_H \ No newline at end of file diff --git a/src/poly/utils/meson.build b/src/poly/utils/meson.build index fc072eb..80f4c95 100644 --- a/src/poly/utils/meson.build +++ b/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( - include_directories : include_directories('.'), +polyutils_dep = declare_dependency( + include_directories : include_directories('./public'), link_with : libpolyutils, sources : polyutils_sources, dependencies : [mfem_dep] diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp index 781efab..4158dbf 100644 --- a/src/poly/utils/private/polyMFEMUtils.cpp +++ b/src/poly/utils/private/polyMFEMUtils.cpp @@ -2,174 +2,302 @@ #include #include #include +#include #include "polyMFEMUtils.h" -NonlinearPowerIntegrator::NonlinearPowerIntegrator( - mfem::FunctionCoefficient &coeff, - double n) : coeff_(coeff), polytropicIndex(n) { +namespace polyMFEMUtils { + NonlinearPowerIntegrator::NonlinearPowerIntegrator( + mfem::Coefficient &coeff, + double n) : coeff_(coeff), polytropicIndex(n) { -} - -void NonlinearPowerIntegrator::AssembleElementVector( - const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::Vector &elvect) { - - const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); - int dof = el.GetDof(); - elvect.SetSize(dof); - elvect = 0.0; - - mfem::Vector shape(dof); - - for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { - mfem::IntegrationPoint ip = ir->IntPoint(iqp); - Trans.SetIntPoint(&ip); - double weight = ip.weight * Trans.Weight(); - - el.CalcShape(ip, shape); - - double u_val = 0.0; - for (int j = 0; j < dof; j++) { - u_val += elfun(j) * shape(j); - } - double u_safe = std::max(u_val, 0.0); - double u_nl = std::pow(u_safe, polytropicIndex); - - double coeff_val = coeff_.Eval(Trans, ip); - double x2_u_nl = coeff_val * u_nl; - - for (int i = 0; i < dof; i++){ - elvect(i) += shape(i) * x2_u_nl * weight; - } } -} -void NonlinearPowerIntegrator::AssembleElementGrad ( - const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::DenseMatrix &elmat) { - - const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); - int dof = el.GetDof(); - elmat.SetSize(dof); - elmat = 0.0; - mfem::Vector shape(dof); - - for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { - mfem::IntegrationPoint ip = ir->IntPoint(iqp); - Trans.SetIntPoint(&ip); - double weight = ip.weight * Trans.Weight(); - - el.CalcShape(ip, shape); - - double u_val = 0.0; - - for (int j = 0; j < dof; j++) { - u_val += elfun(j) * shape(j); - } - double coeff_val = coeff_.Eval(Trans, ip); + void NonlinearPowerIntegrator::AssembleElementVector( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::Vector &elvect) { + const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); + int dof = el.GetDof(); + elvect.SetSize(dof); + elvect = 0.0; - // Calculate the Jacobian - double u_safe = std::max(u_val, 0.0); - double d_u_nl = coeff_val * polytropicIndex * std::pow(u_safe, polytropicIndex - 1); - double x2_d_u_nl = d_u_nl; + mfem::Vector shape(dof); - for (int i = 0; i < dof; i++) { + for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { + mfem::IntegrationPoint ip = ir->IntPoint(iqp); + Trans.SetIntPoint(&ip); + double weight = ip.weight * Trans.Weight(); + + el.CalcShape(ip, shape); + + double u_val = 0.0; for (int j = 0; j < dof; j++) { - elmat(i, j) += shape(i) * x2_d_u_nl * shape(j) * weight; + u_val += elfun(j) * shape(j); + } + double u_safe = std::max(u_val, 0.0); + double u_nl = std::pow(u_safe, polytropicIndex); + + double coeff_val = coeff_.Eval(Trans, ip); + double x2_u_nl = coeff_val * u_nl; + + for (int i = 0; i < dof; i++){ + elvect(i) += shape(i) * x2_u_nl * weight; } } - } -} -BilinearIntegratorWrapper::BilinearIntegratorWrapper( - mfem::BilinearFormIntegrator *integratorInput - ) : integrator(integratorInput) { } + void NonlinearPowerIntegrator::AssembleElementGrad ( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::DenseMatrix &elmat) { -BilinearIntegratorWrapper::~BilinearIntegratorWrapper() { - delete integrator; -} + const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); + int dof = el.GetDof(); + elmat.SetSize(dof); + elmat = 0.0; + mfem::Vector shape(dof); -void BilinearIntegratorWrapper::AssembleElementVector( - const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::Vector &elvect) { - int dof = el.GetDof(); - mfem::DenseMatrix elMat(dof); - integrator->AssembleElementMatrix(el, Trans, elMat); - elvect.SetSize(dof); - elvect = 0.0; - for (int i = 0; i < dof; i++) - { - double sum = 0.0; - for (int j = 0; j < dof; j++) - { - sum += elMat(i, j) * elfun(j); - } - elvect(i) = sum; + for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { + const mfem::IntegrationPoint &ip = ir->IntPoint(iqp); + Trans.SetIntPoint(&ip); + double weight = ip.weight * Trans.Weight(); + + el.CalcShape(ip, shape); + + double u_val = 0.0; + + for (int j = 0; j < dof; j++) { + u_val += elfun(j) * shape(j); + } + double coeff_val = coeff_.Eval(Trans, ip); + + + // Calculate the Jacobian + double u_safe = std::max(u_val, 0.0); + double d_u_nl = coeff_val * polytropicIndex * std::pow(u_safe, polytropicIndex - 1); + double x2_d_u_nl = d_u_nl; + + for (int i = 0; i < dof; i++) { + for (int j = 0; j < dof; j++) { + elmat(i, j) += shape(i) * x2_d_u_nl * shape(j) * weight; + } + } + + } } -} -void BilinearIntegratorWrapper::AssembleElementGrad(const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::DenseMatrix &elmat) { - int dof = el.GetDof(); - elmat.SetSize(dof, dof); - elmat = 0.0; - integrator->AssembleElementMatrix(el, Trans, elmat); -} + BilinearIntegratorWrapper::BilinearIntegratorWrapper( + mfem::BilinearFormIntegrator *integratorInput + ) : integrator(integratorInput) { } -CompositeNonlinearIntegrator::CompositeNonlinearIntegrator() { } - - -CompositeNonlinearIntegrator::~CompositeNonlinearIntegrator() { - for (size_t i = 0; i < integrators.size(); i++) { - delete integrators[i]; + BilinearIntegratorWrapper::~BilinearIntegratorWrapper() { + delete integrator; } -} -void CompositeNonlinearIntegrator::add_integrator(mfem::NonlinearFormIntegrator *integrator) { - integrators.push_back(integrator); -} - -void CompositeNonlinearIntegrator::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 temp(dof); - - for (size_t i = 0; i < integrators.size(); i++) { - temp= 0.0; - integrators[i]->AssembleElementVector(el, Trans, elfun, temp); - elvect.Add(1.0, temp); + void BilinearIntegratorWrapper::AssembleElementVector( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::Vector &elvect) { + int dof = el.GetDof(); + mfem::DenseMatrix elMat(dof); + integrator->AssembleElementMatrix(el, Trans, elMat); + elvect.SetSize(dof); + elvect = 0.0; + for (int i = 0; i < dof; i++) + { + double sum = 0.0; + for (int j = 0; j < dof; j++) + { + sum += elMat(i, j) * elfun(j); + } + elvect(i) = sum; + } } -} -void CompositeNonlinearIntegrator::AssembleElementGrad( - const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::DenseMatrix &elmat) { - int dof = el.GetDof(); - elmat.SetSize(dof, dof); - elmat = 0.0; - mfem::DenseMatrix temp(dof); - temp.SetSize(dof, dof); - for (size_t i = 0; i < integrators.size(); i++) { - temp = 0.0; - integrators[i] -> AssembleElementGrad(el, Trans, elfun, temp); - elmat.Add(1.0, temp); + void BilinearIntegratorWrapper::AssembleElementGrad(const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::DenseMatrix &elmat) { + int dof = el.GetDof(); + elmat.SetSize(dof, dof); + elmat = 0.0; + integrator->AssembleElementMatrix(el, Trans, elmat); } -} \ No newline at end of file + + CompositeNonlinearIntegrator::CompositeNonlinearIntegrator() { } + + + CompositeNonlinearIntegrator::~CompositeNonlinearIntegrator() { + for (size_t i = 0; i < integrators.size(); i++) { + delete integrators[i]; + } + } + + void CompositeNonlinearIntegrator::add_integrator(mfem::NonlinearFormIntegrator *integrator) { + integrators.push_back(integrator); + } + + void CompositeNonlinearIntegrator::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 temp(dof); + + for (size_t i = 0; i < integrators.size(); i++) { + temp= 0.0; + integrators[i]->AssembleElementVector(el, Trans, elfun, temp); + elvect.Add(1.0, temp); + } + } + + void CompositeNonlinearIntegrator::AssembleElementGrad( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::DenseMatrix &elmat) { + int dof = el.GetDof(); + elmat.SetSize(dof, dof); + elmat = 0.0; + mfem::DenseMatrix temp(dof); + temp.SetSize(dof, dof); + for (size_t i = 0; i < integrators.size(); i++) { + temp = 0.0; + integrators[i] -> AssembleElementGrad(el, Trans, elfun, temp); + 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(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 \ No newline at end of file diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h index ac07264..3b7097f 100644 --- a/src/poly/utils/public/polyMFEMUtils.h +++ b/src/poly/utils/public/polyMFEMUtils.h @@ -2,109 +2,32 @@ #include -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. + */ + /** - * @brief A class for nonlinear power integrator. + * @namespace polyMFEMUtils + * @brief A namespace for utilities for working with MFEM and solving the lane-emden equation. */ -class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator { -private: - mfem::FunctionCoefficient coeff_; - double polytropicIndex; -public: +namespace polyMFEMUtils { /** - * @brief Constructor for NonlinearPowerIntegrator. - * - * @param coeff The function coefficient. - * @param n The polytropic index. + * @brief A class for nonlinear power integrator. */ - NonlinearPowerIntegrator(mfem::FunctionCoefficient &coeff, double n); - - /** - * @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; -}; - -/** - * @brief A wrapper class for bilinear integrator. - */ -class BilinearIntegratorWrapper : public mfem::NonlinearFormIntegrator -{ -private: - mfem::BilinearFormIntegrator *integrator; -public: - /** - * @brief Constructor for BilinearIntegratorWrapper. - * - * @param integratorInput The bilinear form integrator input. - */ - BilinearIntegratorWrapper(mfem::BilinearFormIntegrator *integratorInput); - - /** - * @brief Destructor for BilinearIntegratorWrapper. - */ - virtual ~BilinearIntegratorWrapper(); - - /** - * @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; -}; - -/** - * @brief A class for composite nonlinear integrator. - */ -class CompositeNonlinearIntegrator: public mfem::NonlinearFormIntegrator { - private: - std::vector integrators; + class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator { + private: + mfem::Coefficient &coeff_; + double polytropicIndex; public: /** - * @brief Constructor for CompositeNonlinearIntegrator. - */ - CompositeNonlinearIntegrator(); - - /** - * @brief Destructor for CompositeNonlinearIntegrator. - */ - virtual ~CompositeNonlinearIntegrator(); - - /** - * @brief Adds an integrator to the composite integrator. + * @brief Constructor for NonlinearPowerIntegrator. * - * @param integrator The nonlinear form integrator to add. + * @param coeff The function coefficient. + * @param n The polytropic index. */ - void add_integrator(mfem::NonlinearFormIntegrator *integrator); + NonlinearPowerIntegrator(mfem::Coefficient &coeff, double n); /** * @brief Assembles the element vector. @@ -124,5 +47,153 @@ class CompositeNonlinearIntegrator: public mfem::NonlinearFormIntegrator { * @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; + }; + + /** + * @brief A wrapper class for bilinear integrator. + */ + class BilinearIntegratorWrapper : public mfem::NonlinearFormIntegrator + { + private: + mfem::BilinearFormIntegrator *integrator; + public: + /** + * @brief Constructor for BilinearIntegratorWrapper. + * + * @param integratorInput The bilinear form integrator input. + */ + BilinearIntegratorWrapper(mfem::BilinearFormIntegrator *integratorInput); + + /** + * @brief Destructor for BilinearIntegratorWrapper. + */ + virtual ~BilinearIntegratorWrapper(); + + /** + * @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; -}; \ No newline at end of file + }; + + /** + * @brief A class for composite nonlinear integrator. + */ + class CompositeNonlinearIntegrator: public mfem::NonlinearFormIntegrator { + private: + std::vector integrators; + public: + /** + * @brief Constructor for CompositeNonlinearIntegrator. + */ + CompositeNonlinearIntegrator(); + + /** + * @brief Destructor for CompositeNonlinearIntegrator. + */ + virtual ~CompositeNonlinearIntegrator(); + + /** + * @brief Adds an integrator to the composite integrator. + * + * @param integrator The nonlinear form integrator to add. + */ + void add_integrator(mfem::NonlinearFormIntegrator *integrator); + + /** + * @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; + }; + + /** + * @brief A class for constraint integrator. + */ + class ConstraintIntegrator: public mfem::NonlinearFormIntegrator { + private: + mfem::Coefficient η + + 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 nfl; + std::unique_ptr 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 \ No newline at end of file