From e3afe90f37d93cb30cbbc2e25b9b418d40bee43b Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 2 Apr 2025 14:57:37 -0400 Subject: [PATCH] feat(poly): moved to a block form for poly essential dofs can be applied to both theta and phi (grad theta) if we move to a block form. I have done this derivation and made that change so that we can properly apply the central boundary condition to the slope --- src/poly/solver/meson.build | 16 +- src/poly/solver/private/polySolver.cpp | 300 ++++++++++++++------- src/poly/solver/public/polySolver.h | 71 ++--- src/poly/utils/meson.build | 19 +- src/poly/utils/private/integrators.cpp | 118 ++++++++ src/poly/utils/private/operator.cpp | 122 +++++++++ src/poly/utils/private/polyIO.cpp | 43 --- src/poly/utils/private/polyMFEMUtils.cpp | 326 ----------------------- src/poly/utils/public/integrators.h | 85 ++++++ src/poly/utils/public/operator.h | 45 ++++ src/poly/utils/public/polyIO.h | 36 --- src/poly/utils/public/polyMFEMUtils.h | 190 ------------- 12 files changed, 640 insertions(+), 731 deletions(-) create mode 100644 src/poly/utils/private/integrators.cpp create mode 100644 src/poly/utils/private/operator.cpp delete mode 100644 src/poly/utils/private/polyIO.cpp delete mode 100644 src/poly/utils/private/polyMFEMUtils.cpp create mode 100644 src/poly/utils/public/integrators.h create mode 100644 src/poly/utils/public/operator.h delete mode 100644 src/poly/utils/public/polyIO.h delete mode 100644 src/poly/utils/public/polyMFEMUtils.h diff --git a/src/poly/solver/meson.build b/src/poly/solver/meson.build index 246d414..ea20457 100644 --- a/src/poly/solver/meson.build +++ b/src/poly/solver/meson.build @@ -26,11 +26,23 @@ polySolver_headers = files( 'public/polySolver.h' ) +dependencies = [ + mfem_dep, + meshio_dep, + polycoeff_dep, + polyutils_dep, + macros_dep, + probe_dep, + quill_dep, + config_dep, + resourceManager_dep +] + libPolySolver = static_library('polySolver', polySolver_sources, include_directories : include_directories('./public'), cpp_args: ['-fvisibility=default'], - dependencies: [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep, macros_dep, probe_dep, quill_dep, config_dep], + dependencies: dependencies, install: true ) @@ -38,5 +50,5 @@ polysolver_dep = declare_dependency( include_directories : include_directories('./public'), link_with : libPolySolver, sources : polySolver_sources, - dependencies : [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep, macros_dep, probe_dep, quill_dep, config_dep] + dependencies : dependencies ) \ No newline at end of file diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 3e5f95a..4b24351 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -23,12 +23,16 @@ #include #include #include +#include #include "polySolver.h" -#include "polyMFEMUtils.h" +#include "integrators.h" #include "polyCoeff.h" -#include "probe.h" #include "config.h" +#include "probe.h" +#include "resourceManager.h" +#include "resourceManagerTypes.h" +#include "operator.h" #include "quill/LogMacros.h" @@ -63,116 +67,152 @@ namespace laneEmden { } } -PolySolver::PolySolver(double n, double order, mfem::Mesh& mesh_) -: logger(logManager.getLogger("log")), - n(n), - order(order), - mesh(mesh_), - feCollection(std::make_unique(order, mesh.SpaceDimension())), - feSpace(std::make_unique(&mesh, feCollection.get())), - compositeIntegrator(std::make_unique()), - nonlinearForm(std::make_unique(feSpace.get())), - u(std::make_unique(feSpace.get())) { +PolySolver::PolySolver(double n, double order) { // --- 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)); + LOG_ERROR(m_logger, "The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is {}", m_polytropicIndex); + 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(m_polytropicIndex)); } - diffusionCoeff = std::make_unique(mesh.SpaceDimension(), polycoeff::diffusionCoeff); - nonlinearSourceCoeff = std::make_unique(polycoeff::nonlinearSourceCoeff); + m_polytropicIndex = n; + m_feOrder = order; - assembleNonlinearForm(); + ResourceManager& rm = ResourceManager::getInstance(); + const Resource& resource = rm.getResource("mesh:polySphere"); + const auto &meshIO = std::get>(resource); + meshIO->LinearRescale(polycoeff::x1(m_polytropicIndex)); + m_mesh = std::make_unique(meshIO->GetMesh()); + + // Use feOrder - 1 for the RT space to satisfy Brezzi-Babuska condition + // for the H1 and RT spaces + m_fecH1 = std::make_unique(m_feOrder, m_mesh->SpaceDimension()); + m_fecRT = std::make_unique(m_feOrder - 1, m_mesh->SpaceDimension()); + + m_feTheta = std::make_unique(m_mesh.get(), m_fecH1.get()); + m_fePhi = std::make_unique(m_mesh.get(), m_fecRT.get()); + + m_theta = std::make_unique(m_feTheta.get()); + m_phi = std::make_unique(m_fePhi.get()); + + assembleBlockSystem(); } PolySolver::~PolySolver() {} -void PolySolver::assembleNonlinearForm() { - // Add the \int_{\Omega}\nabla v\cdot\nabla\theta d\Omegaterm - auto wrappedDiffusionIntegrator = std::make_unique( - new mfem::DiffusionIntegrator(*diffusionCoeff) - ); - compositeIntegrator->add_integrator(wrappedDiffusionIntegrator.release()); +void PolySolver::assembleBlockSystem() { - // Add the \int_{\Omega}v\theta^{n} d\Omega term - auto nonlinearIntegrator = std::make_unique(*nonlinearSourceCoeff, n); - compositeIntegrator->add_integrator(nonlinearIntegrator.release()); + // Start by defining the block structure of the system + // Block 0: Theta (scalar space, uses m_feTheta) + // Block 1: Phi (vector space, uses m_fePhi) + mfem::Array feSpaces; + feSpaces.Append(m_feTheta.get()); + feSpaces.Append(m_fePhi.get()); - // Add the contraint term \gamma(\nabla \theta(0)\cdot\nabla v(0))^{2} - double gamma = config.get("Poly:Solver:Constraint:Gamma", 1e4); - auto constraintIntegrator = std::make_unique( - new polyMFEMUtils::ConstraintIntegrator(gamma, &mesh) - ); - compositeIntegrator->add_integrator(constraintIntegrator.release()); + // Create the block offsets. These define the start of each block in the combined vector. + // Block offsets will be [0, thetaDofs, thetaDofs + phiDofs] + mfem::Array blockOffsets; + blockOffsets.SetSize(3); + blockOffsets[0] = 0; + blockOffsets[1] = feSpaces[0]->GetVSize(); + blockOffsets[2] = feSpaces[1]->GetVSize(); + blockOffsets.PartialSum(); + + // Coefficients + mfem::ConstantCoefficient negOneCoeff(-1.0); + mfem::ConstantCoefficient oneCoeff(1.0); + + mfem::Vector negOneVec(mfem::Vector(3)); + mfem::Vector oneVec(mfem::Vector(3)); + negOneVec = -1.0; + oneVec = 1.0; + + mfem::VectorConstantCoefficient negOneVCoeff(negOneVec); + mfem::VectorConstantCoefficient oneVCoeff(oneVec); + + // Add integrators to block form one by one + // We add integrators cooresponding to each term in the weak form + // The block form of the residual matrix + // ⎡ 0 -M ⎤ ⎡ θ ⎤ + ⎡f(θ)⎤ = ⎡ 0 ⎤ = R(X) + // ⎣ -Q D ⎦ ⎣ Φ ⎦ ⎣ 0 ⎦ ⎣ 0 ⎦ + // This then simplifies to + // ⎡f(θ) - MΘ ⎤ = ⎡ 0 ⎤ = R + // ⎣ -Qɸ Dθ ⎦ ⎣ 0 ⎦ + // Here M, Q, and D are + // M = ∫∇ψᶿ·Nᵠ dV (MixedVectorWeakDivergenceIntegrator) + // D = ∫ψᵠ·Nᵠ dV (VectorFEMassIntegrator) + // Q = ∫ψᵠ·∇Nᶿ dV (MixedVectorGradientIntegrator) + // f(θ) = ∫ψᶿ·θⁿ dV (NonlinearPowerIntegrator) + // Here ψᶿ and ψᵠ are the test functions for the theta and phi spaces, respectively + // Nᵠ and Nᶿ are the basis functions for the theta and phi spaces, respectively + // A full derivation of the weak form can be found in the 4DSSE documentation + + // --- Assemble the MixedBilinear and Bilinear forms (M, D, and Q) --- + auto Mform = std::make_unique(m_feTheta.get(), m_fePhi.get()); + auto Qform = std::make_unique(m_fePhi.get(), m_feTheta.get()); + auto Dform = std::make_unique(m_fePhi.get()); + + // TODO: Check the sign on all of the integrators + Mform->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator(negOneCoeff)); + Qform->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator(negOneVCoeff)); + Dform->AddDomainIntegrator(new mfem::VectorFEMassIntegrator(oneCoeff)); + + Mform->Assemble(); + Mform->Finalize(); + + Qform->Assemble(); + Qform->Finalize(); + + Dform->Assemble(); + Dform->Finalize(); + + // --- Assemble the NonlinearForm (f) --- + auto fform = std::make_unique(m_feTheta.get()); + fform->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(oneCoeff, m_polytropicIndex)); + // TODO: Add essential boundary conditions to the nonlinear form + + // -- Build the BlockOperator -- + m_polytropOperator = std::make_unique( + std::move(Mform), + std::move(Qform), + std::move(Dform), + std::move(fform), + blockOffsets + ); - nonlinearForm->AddDomainIntegrator(compositeIntegrator.release()); } void PolySolver::solve(){ // --- Set the initial guess for the solution --- - mfem::FunctionCoefficient initCoeff ( - [this](const mfem::Vector &x) { - double r = x.Norml2(); - // double theta = laneEmden::thetaSerieseExpansion(r, n, 10); - // return theta; - double radius = Probe::getMeshRadius(mesh); - double u = 1/radius; + setInitialGuess(); - return -std::pow((u*r), 2)+1.0; - } - ); - u->ProjectCoefficient(initCoeff); - if (config.get("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; - centerPoint(1, 0) = 0.0; - centerPoint(2, 0) = 0.0; + // --- Set the essential true dofs for the operator --- + mfem::Array theta_ess_tdof_list, phi_ess_tdof_list; + std::tie(theta_ess_tdof_list, phi_ess_tdof_list) = getEssentialTrueDof(); + m_polytropOperator->SetEssentialTrueDofs(theta_ess_tdof_list, phi_ess_tdof_list); - mfem::Array elementIDs; - mfem::Array ips; - mesh.FindPoints(centerPoint, elementIDs, ips); - mfem::Array centerDofs; - mfem::Array tempDofs; - for (int i = 0; i < elementIDs.Size(); i++) { - feSpace->GetElementDofs(elementIDs[i], tempDofs); - centerDofs.Append(tempDofs); - } - mfem::Array ess_tdof_list; - mfem::Array ess_brd(mesh.bdr_attributes.Max()); - ess_brd = 1; - feSpace->GetEssentialTrueDofs(ess_brd, ess_tdof_list); - // combine the essential dofs with the center dofs - ess_tdof_list.Append(centerDofs); - nonlinearForm->SetEssentialTrueDofs(ess_tdof_list); - // Set the center elemID to be the Dirichlet boundary + // --- Load configuration parameters --- + double newtonRelTol = m_config.get("Poly:Solver:Newton:RelTol", 1e-7); + double newtonAbsTol = m_config.get("Poly:Solver:Newton:AbsTol", 1e-7); + int newtonMaxIter = m_config.get("Poly:Solver:Newton:MaxIter", 200); + int newtonPrintLevel = m_config.get("Poly:Solver:Newton:PrintLevel", 1); - // double alpha = config.get("Poly:Solver:Newton:Alpha", 1e2); - double newtonRelTol = config.get("Poly:Solver:Newton:RelTol", 1e-7); - double newtonAbsTol = config.get("Poly:Solver:Newton:AbsTol", 1e-7); - int newtonMaxIter = config.get("Poly:Solver:Newton:MaxIter", 200); - int newtonPrintLevel = config.get("Poly:Solver:Newton:PrintLevel", 1); + double gmresRelTol = m_config.get("Poly:Solver:GMRES:RelTol", 1e-10); + double gmresAbsTol = m_config.get("Poly:Solver:GMRES:AbsTol", 1e-12); + int gmresMaxIter = m_config.get("Poly:Solver:GMRES:MaxIter", 2000); + int gmresPrintLevel = m_config.get("Poly:Solver:GMRES:PrintLevel", 0); - double gmresRelTol = config.get("Poly:Solver:GMRES:RelTol", 1e-10); - double gmresAbsTol = config.get("Poly:Solver:GMRES:AbsTol", 1e-12); - int gmresMaxIter = config.get("Poly:Solver:GMRES:MaxIter", 2000); - int gmresPrintLevel = config.get("Poly:Solver:GMRES:PrintLevel", 0); + LOG_DEBUG(m_logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel); + LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); - LOG_INFO(logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel); - LOG_INFO(logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); - - std::vector zeroSlopeCoordinate = {0.0, 0.0, 0.0}; - // polyMFEMUtils::ZeroSlopeNewtonSolver newtonSolver(alpha, zeroSlopeCoordinate); + // --- Set up the Newton solver --- mfem::NewtonSolver newtonSolver; newtonSolver.SetRelTol(newtonRelTol); newtonSolver.SetAbsTol(newtonAbsTol); newtonSolver.SetMaxIter(newtonMaxIter); newtonSolver.SetPrintLevel(newtonPrintLevel); - newtonSolver.SetOperator(*nonlinearForm); + newtonSolver.SetOperator(*m_polytropOperator); mfem::GMRESSolver gmresSolver; gmresSolver.SetRelTol(gmresRelTol); gmresSolver.SetAbsTol(gmresAbsTol); @@ -181,24 +221,96 @@ void PolySolver::solve(){ newtonSolver.SetSolver(gmresSolver); // newtonSolver.SetAdaptiveLinRtol(); - mfem::Vector B(feSpace->GetTrueVSize()); + mfem::Vector B(m_feTheta->GetTrueVSize()); B = 0.0; - newtonSolver.Mult(B, *u); + newtonSolver.Mult(B, *m_theta); - Probe::glVisView(*u, mesh, "solution"); + // --- Save and view the solution --- + saveAndViewSolution(); + + +} + +std::pair, mfem::Array> PolySolver::getEssentialTrueDof() { + mfem::Array theta_ess_tdof_list; + mfem::Array phi_ess_tdof_list; + + mfem::Array centerDofs = findCenterElement(); + + phi_ess_tdof_list.Append(centerDofs); + + mfem::Array ess_brd(m_mesh->bdr_attributes.Max()); + ess_brd = 1; + m_feTheta->GetEssentialTrueDofs(ess_brd, theta_ess_tdof_list); + // combine the essential dofs with the center dofs + theta_ess_tdof_list.Append(centerDofs); + + return std::make_pair(theta_ess_tdof_list, phi_ess_tdof_list); +} + +mfem::Array PolySolver::findCenterElement() { + mfem::Array centerDofs; + mfem::DenseMatrix centerPoint(m_mesh->SpaceDimension(), 1); + centerPoint(0, 0) = 0.0; + centerPoint(1, 0) = 0.0; + centerPoint(2, 0) = 0.0; + + mfem::Array elementIDs; + mfem::Array ips; + m_mesh->FindPoints(centerPoint, elementIDs, ips); + mfem::Array tempDofs; + for (int i = 0; i < elementIDs.Size(); i++) { + m_feTheta->GetElementDofs(elementIDs[i], tempDofs); + centerDofs.Append(tempDofs); + } + return centerDofs; +} + +void PolySolver::setInitialGuess() { + // --- Set the initial guess for the solution --- + mfem::FunctionCoefficient thetaInitGuess ( + [this](const mfem::Vector &x) { + double r = x.Norml2(); + double radius = Probe::getMeshRadius(*m_mesh); + double u = 1/radius; + + return -std::pow((u*r), 2)+1.0; + } + ); + mfem::VectorFunctionCoefficient phiInitGuess (m_mesh->SpaceDimension(), + [this](const mfem::Vector &x, mfem::Vector &v) { + double radius = Probe::getMeshRadius(*m_mesh); + double u = -1/std::pow(radius,2); + v(0) = 2*std::abs(x(0))*u; + v(1) = 2*std::abs(x(1))*u; + v(2) = 2*std::abs(x(2))*u; + } + ); + m_theta->ProjectCoefficient(thetaInitGuess); + m_phi->ProjectCoefficient(phiInitGuess); + if (m_config.get("Poly:Solver:ViewInitialGuess", false)) { + Probe::glVisView(*m_theta, *m_mesh, "initialGuess"); + } + +} + +void PolySolver::saveAndViewSolution() { + bool doView = m_config.get("Poly:Output:View", false); + if (doView) { + Probe::glVisView(*m_theta, *m_mesh, "solution"); + } // --- Extract the Solution --- - bool write11DSolution = config.get("Poly:Output:1D:Save", true); + bool write11DSolution = m_config.get("Poly:Output:1D:Save", true); if (write11DSolution) { - std::string solutionPath = config.get("Poly:Output:1D:Path", "polytropeSolution_1D.csv"); - double rayCoLatitude = config.get("Poly:Output:1D:RayCoLatitude", 0.0); - double rayLongitude = config.get("Poly:Output:1D:RayLongitude", 0.0); - int raySamples = config.get("Poly:Output:1D:RaySamples", 100); + std::string solutionPath = m_config.get("Poly:Output:1D:Path", "polytropeSolution_1D.csv"); + double rayCoLatitude = m_config.get("Poly:Output:1D:RayCoLatitude", 0.0); + double rayLongitude = m_config.get("Poly:Output:1D:RayLongitude", 0.0); + int raySamples = m_config.get("Poly:Output:1D:RaySamples", 100); std::vector rayDirection = {rayCoLatitude, rayLongitude}; - Probe::getRaySolution(*u, *feSpace, rayDirection, raySamples, solutionPath); + Probe::getRaySolution(*m_theta, *m_feTheta, rayDirection, raySamples, solutionPath); } - } \ No newline at end of file diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 0e9a75c..92ceec9 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -2,13 +2,11 @@ #define POLYSOLVER_H #include "mfem.hpp" -#include -#include #include +#include -#include "meshIO.h" -#include "polyCoeff.h" -#include "polyMFEMUtils.h" +#include "integrators.h" +#include "operator.h" #include "config.h" #include "probe.h" #include "quill/Logger.h" @@ -20,35 +18,44 @@ namespace laneEmden { } class PolySolver { -private: - Config& config = Config::getInstance(); - Probe::LogManager& logManager = Probe::LogManager::getInstance(); - quill::Logger* logger; - double n, order; - mfem::Mesh& mesh; - - std::unique_ptr feCollection; - std::unique_ptr feSpace; - - std::unique_ptr compositeIntegrator; - std::unique_ptr nonlinearForm; - - std::unique_ptr u; - - std::unique_ptr diffusionCoeff; - std::unique_ptr nonlinearSourceCoeff; - - void assembleNonlinearForm(); - - -public: - PolySolver(double n, double order, mfem::Mesh& mesh_); +public: // Public methods + PolySolver(double n, double order); ~PolySolver(); + void solve(); - mfem::Mesh& getMesh() { return mesh; } - mfem::GridFunction& getSolution() { return *u; } - double getN() { return n; } - double getOrder() { return order; } + + double getN() { return m_polytropicIndex; } + double getOrder() { return m_feOrder; } + mfem::Mesh* getMesh() { return m_mesh.get(); } + mfem::GridFunction& getSolution() { return *m_theta; } + +private: // Private Attributes + Config& m_config = Config::getInstance(); + Probe::LogManager& m_logManager = Probe::LogManager::getInstance(); + quill::Logger* m_logger = m_logManager.getLogger("log"); + double m_polytropicIndex, m_feOrder; + + std::unique_ptr m_mesh; + std::unique_ptr m_fecH1; + std::unique_ptr m_fecRT; + + std::unique_ptr m_feTheta; + std::unique_ptr m_fePhi; + + std::unique_ptr m_theta; + std::unique_ptr m_phi; + + std::unique_ptr m_polytropOperator; + + +private: // Private methods + void assembleBlockSystem(); + std::pair, mfem::Array> getEssentialTrueDof(); + mfem::Array findCenterElement(); + void setInitialGuess(); + void saveAndViewSolution(); + + }; #endif // POLYSOLVER_H \ No newline at end of file diff --git a/src/poly/utils/meson.build b/src/poly/utils/meson.build index 1f3f150..e97ca94 100644 --- a/src/poly/utils/meson.build +++ b/src/poly/utils/meson.build @@ -19,20 +19,23 @@ # # *********************************************************************** # polyutils_sources = files( - 'private/polyIO.cpp', - 'private/polyMFEMUtils.cpp' + 'private/integrators.cpp', + 'private/operator.cpp' ) -polyutils_headers = files( - 'public/polyIO.h', - 'public/polyMFEMUtils.h' -) +dependencies = [ + mfem_dep, + macros_dep, + probe_dep, + quill_dep, + config_dep, +] libpolyutils = static_library('polyutils', polyutils_sources, include_directories : include_directories('./public'), cpp_args: ['-fvisibility=default'], - dependencies: [mfem_dep, macros_dep, probe_dep, quill_dep, config_dep], + dependencies: dependencies, install: true ) @@ -40,5 +43,5 @@ polyutils_dep = declare_dependency( include_directories : include_directories('./public'), link_with : libpolyutils, sources : polyutils_sources, - dependencies : [mfem_dep, macros_dep, probe_dep, quill_dep, config_dep] + dependencies : dependencies ) diff --git a/src/poly/utils/private/integrators.cpp b/src/poly/utils/private/integrators.cpp new file mode 100644 index 0000000..d678c25 --- /dev/null +++ b/src/poly/utils/private/integrators.cpp @@ -0,0 +1,118 @@ +/* *********************************************************************** +// +// Copyright (C) 2025 -- The 4D-STAR Collaboration +// File Author: Emily Boudreaux +// Last Modified: March 19, 2025 +// +// 4DSSE is free software; you can use it and/or modify +// it under the terms and restrictions the GNU General Library Public +// License version 3 (GPLv3) as published by the Free Software Foundation. +// +// 4DSSE is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU Library General Public License for more details. +// +// You should have received a copy of the GNU Library General Public License +// along with this software; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// *********************************************************************** */ +#include "mfem.hpp" +#include +#include +#include +#include + +#include "quill/LogMacros.h" + +#include "integrators.h" +#include "debug.h" + + +namespace polyMFEMUtils { + NonlinearPowerIntegrator::NonlinearPowerIntegrator( + mfem::Coefficient &coeff, + double n) : + + m_coeff(coeff), + m_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, m_polytropicIndex); + + double coeff_val = m_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++) { + 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 = m_coeff.Eval(Trans, ip); + + + // Calculate the Jacobian + double u_safe = std::max(u_val, 0.0); + double d_u_nl = coeff_val * m_polytropicIndex * std::pow(u_safe, m_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; + } + } + + } + } +} // namespace polyMFEMUtils \ No newline at end of file diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp new file mode 100644 index 0000000..4de7639 --- /dev/null +++ b/src/poly/utils/private/operator.cpp @@ -0,0 +1,122 @@ +#include "operator.h" +#include "mfem.hpp" +#include "linalg/vector.hpp" +#include + +PolytropeOperator::PolytropeOperator( + + std::unique_ptr M, + std::unique_ptr Q, + std::unique_ptr D, + std::unique_ptr f, + const mfem::Array &blockOffsets) : + + mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector + m_blockOffsets(blockOffsets), + m_jacobian(nullptr) { + + m_M = std::move(M); + m_Q = std::move(Q); + m_D = std::move(D); + m_f = std::move(f); + + + m_Mmat = std::make_unique(m_M->SpMat()); + m_Qmat = std::make_unique(m_Q->SpMat()); + m_Dmat = std::make_unique(m_D->SpMat()); + + m_negM_op = std::make_unique(m_Mmat.get(), -1.0); + m_negQ_op = std::make_unique(m_Qmat.get(), -1.0); + + MFEM_ASSERT(m_Mmat.get() != nullptr, "Matrix m_Mmat is null in PolytropeOperator constructor"); + MFEM_ASSERT(m_Qmat.get() != nullptr, "Matrix m_Qmat is null in PolytropeOperator constructor"); + MFEM_ASSERT(m_Dmat.get() != nullptr, "Matrix m_Dmat is null in PolytropeOperator constructor"); + MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator constructor"); +} + + +void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { + // -- Create BlockVector views for input x and output y + mfem::BlockVector x_block(const_cast(x), m_blockOffsets); + mfem::BlockVector y_block(y, m_blockOffsets); + + // -- Get Vector views for individual blocks + const mfem::Vector &x_theta = x_block.GetBlock(0); + const mfem::Vector &x_phi = x_block.GetBlock(1); + mfem::Vector &y_R0 = y_block.GetBlock(0); // Residual Block 0 (theta) + mfem::Vector &y_R1 = y_block.GetBlock(1); // Residual Block 1 (phi) + + int theta_size = m_blockOffsets[1] - m_blockOffsets[0]; + int phi_size = m_blockOffsets[2] - m_blockOffsets[1]; + + mfem::Vector f_term(theta_size); + mfem::Vector Mphi_term(theta_size); + mfem::Vector Dphi_term(phi_size); + mfem::Vector Qtheta_term(phi_size); + + // Caucluate R0 and R1 terms + // R0 = f(θ) - Mɸ + // R1 = Dɸ - Qθ + + MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult"); + MFEM_ASSERT(m_Mmat.get() != nullptr, "SparseMatrix m_Mmat is null in PolytropeOperator::Mult"); + MFEM_ASSERT(m_Dmat.get() != nullptr, "SparseMatrix m_Dmat is null in PolytropeOperator::Mult"); + MFEM_ASSERT(m_Qmat.get() != nullptr, "SparseMatrix m_Qmat is null in PolytropeOperator::Mult"); + + m_f->Mult(x_theta, f_term); + m_Mmat->Mult(x_phi, Mphi_term); + m_Dmat->Mult(x_phi, Dphi_term); + m_Qmat->Mult(x_theta, Qtheta_term); + + subtract(f_term, Mphi_term, y_R0); + subtract(Dphi_term, Qtheta_term, y_R1); + + + // -- Apply essential boundary conditions -- + for (int i = 0; i < m_theta_ess_tofs.Size(); i++) { + int idx = m_theta_ess_tofs[i]; + if (idx >= 0 && idx < y_R0.Size()) { + y_block.GetBlock(0)[idx] = 0.0; // Zero out the essential theta dofs in the bilinear form + } + } + + for (int i = 0; i < m_phi_ess_tofs.Size(); i++) { + int idx = m_phi_ess_tofs[i]; + if (idx >= 0 && idx < y_R1.Size()) { + y_block.GetBlock(1)[idx] = 0.0; // Zero out the essential phi dofs in the bilinear form + } + } + +} +mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const { + // -- Get the gradient of f -- + mfem::BlockVector x_block(const_cast(x), m_blockOffsets); + const mfem::Vector& x_theta = x_block.GetBlock(0); + + mfem::Operator& J00 = m_f->GetGradient(x_theta); + if (m_jacobian == nullptr) { + m_jacobian = std::make_unique(m_blockOffsets); + m_jacobian->SetBlock(0, 0, &J00); // df/dθ (state-dependent) + m_jacobian->SetBlock(0, 1, m_negM_op.get()); // -M (constant) + m_jacobian->SetBlock(1, 0, m_negQ_op.get()); // -Q (constant) + m_jacobian->SetBlock(1, 1, m_Dmat.get()); // D (constant) + } else { + // The Jacobian already exists, we only need to update the first block + // since the other blocks have a constant derivitive (they are linear) + m_jacobian->SetBlock(0, 0, &J00); + } + + return *m_jacobian; +} + +void PolytropeOperator::SetEssentialTrueDofs(const mfem::Array &theta_ess_tofs, + const mfem::Array &phi_ess_tofs) { + m_theta_ess_tofs = theta_ess_tofs; + m_phi_ess_tofs = phi_ess_tofs; + + if (m_f) { + m_f->SetEssentialTrueDofs(theta_ess_tofs); + } else { + MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs"); + } +} \ No newline at end of file diff --git a/src/poly/utils/private/polyIO.cpp b/src/poly/utils/private/polyIO.cpp deleted file mode 100644 index 96f0cd9..0000000 --- a/src/poly/utils/private/polyIO.cpp +++ /dev/null @@ -1,43 +0,0 @@ -/* *********************************************************************** -// -// Copyright (C) 2025 -- The 4D-STAR Collaboration -// File Author: Emily Boudreaux -// Last Modified: February 14, 2025 -// -// 4DSSE is free software; you can use it and/or modify -// it under the terms and restrictions the GNU General Library Public -// License version 3 (GPLv3) as published by the Free Software Foundation. -// -// 4DSSE is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -// See the GNU Library General Public License for more details. -// -// You should have received a copy of the GNU Library General Public License -// along with this software; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// *********************************************************************** */ -#include "mfem.hpp" -#include -#include - -#include "polyIO.h" - -void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename) { - std::ofstream file(filename); - if (!file.is_open()) { - std::cerr << "Error: Could not open " << filename << " for writing." << std::endl; - return; - } - - file << "xi,u\n"; // CSV header - - for (int i = 0; i < u.Size(); i++) { - double xi = mesh.GetVertex(i)[0]; // Get spatial coordinate - file << xi << "," << u[i] << "\n"; // Write to CSV - } - - file.close(); - std::cout << "Solution written to " << filename << std::endl; -} \ No newline at end of file diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp deleted file mode 100644 index 33f18b2..0000000 --- a/src/poly/utils/private/polyMFEMUtils.cpp +++ /dev/null @@ -1,326 +0,0 @@ -/* *********************************************************************** -// -// Copyright (C) 2025 -- The 4D-STAR Collaboration -// File Author: Emily Boudreaux -// Last Modified: March 19, 2025 -// -// 4DSSE is free software; you can use it and/or modify -// it under the terms and restrictions the GNU General Library Public -// License version 3 (GPLv3) as published by the Free Software Foundation. -// -// 4DSSE is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -// See the GNU Library General Public License for more details. -// -// You should have received a copy of the GNU Library General Public License -// along with this software; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// *********************************************************************** */ -#include "mfem.hpp" -#include -#include -#include -#include -#include -#include - -#include "quill/LogMacros.h" - -#include "polyMFEMUtils.h" -#include "debug.h" - - -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++) { - 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; - } - } - - } - } - - BilinearIntegratorWrapper::BilinearIntegratorWrapper( - mfem::BilinearFormIntegrator *integratorInput - ) : integrator(integratorInput) { } - - BilinearIntegratorWrapper::~BilinearIntegratorWrapper() { - delete integrator; - } - - 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 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); - } - - CompositeNonlinearIntegrator::CompositeNonlinearIntegrator() { } - - - CompositeNonlinearIntegrator::~CompositeNonlinearIntegrator() { } - - 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); - } - } - - // TODO: break this up into smaller functions - // TODO: think of a more efficient way to find connected elements - ConstraintIntegrator::ConstraintIntegrator(const double gamma, mfem::Mesh* mesh): m_gamma(gamma), m_mesh(mesh) { - 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; - - m_originCoordinateVector.SetSize(3); - m_originCoordinateVector = 0.0; - - m_mesh->FindPoints(m_originCoordinateMatrix, m_originElementIDs, m_originIntegrationPoints); - - if (m_originElementIDs.Size() == 0) { - LOG_ERROR(m_logger, "The origin point is not found in the mesh."); - throw std::runtime_error("The origin point is not found in the mesh."); - } - - LOG_INFO(m_logger, "The origin point is found in the mesh."); - - // NOTE (EMB, March 2025): This function as it is currently written will break if the mesh is refined after being passed to the constructor - // This may or may not be an issue (it does seem unlikley that the mesh would be refined after being passed to the constructor) - // But if something mysteriously breaks in the future this is may be a good place to start looking - mfem::Table* VETable = m_mesh->GetVertexToElementTable(); - const int nVertices = VETable->Size(); - LOG_INFO(m_logger, "The number of vertices in the mesh is {}", nVertices); - std::vector originVertexIds; - mfem::Array connectedElements; - - // -- Get all vertices connected to the origin element -- - for (int vertexID = 0; vertexID < nVertices; vertexID++) { - VETable->GetRow(vertexID, connectedElements); - for (int j = 0; j < connectedElements.Size(); j++) { - if (connectedElements[j] == m_originElementIDs[0]) { - originVertexIds.push_back(vertexID); - } - } - - } - double minDistanceToOrigin = std::numeric_limits::max(); - int minDistanceVertexId = -1; - - // -- Get the vertex closest to the origin ID -- - for (const auto &vertexId : originVertexIds) { - mfem::Vector vertex; - const double* vcoord = m_mesh->GetVertex(vertexId); - // --- Note if this is run with a 2D or 1D mesh this may lead to a segfault --- - // TODO: Add a check for the dimension of the mesh - double distance = vcoord[0]*vcoord[0] + vcoord[1]*vcoord[1] + vcoord[2]*vcoord[2]; - if (distance < minDistanceToOrigin) { - minDistanceToOrigin = distance; - minDistanceVertexId = vertexId; - } - - - } - if (minDistanceVertexId == -1 || minDistanceToOrigin > 1e-10) { - LOG_ERROR(m_logger, "The origin vertex is not found in the mesh."); - throw std::runtime_error("The origin vertex is not found in the mesh."); - } - // -- Find all elements connected to the origin vertex by looping through the VE table - VETable->GetRow(minDistanceVertexId, m_originConnectedElementIds); - - LOG_INFO(m_logger, "Found {} elements connected to the origin vertex.", m_originConnectedElementIds.Size()); - - } - - - void ConstraintIntegrator::AssembleElementMatrix(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, mfem::DenseMatrix &elmat) { - int elemID = Trans.ElementNo; - - // -- Check if the element is connected to the origin vertex -- - bool isConnected = m_originConnectedElementIds.Find(elemID) != -1 ? true : false; - - if (!isConnected) { - elmat = 0.0; - return; - } - - - // -- Compute the derivitives using MFEM's build in shape routines -- - - int numDoF = el.GetDof(); - int dim = m_mesh->Dimension(); - - // -- Map the origin in physical space to the reference space of the element -- - mfem::Vector originReferenceCoordinate(dim); - mfem::IntegrationPoint originIntegrationPoint; - Trans.TransformBack(m_originCoordinateVector, originIntegrationPoint); - - // -- Compute the derivitives of the shape functions at the origin -- - mfem::DenseMatrix dshape(numDoF, dim); - el.CalcDShape(originIntegrationPoint, dshape); - - // -- Transform derivitives from reference space to physical space using the inverse of the Jacobian -- - mfem::DenseMatrix invJac(dim, dim); - invJac = Trans.InverseJacobian(); - - mfem::DenseMatrix dshapePhysical(numDoF, dim); - dshapePhysical = 0.0; - - for (int dofID = 0; dofID < numDoF; dofID++) { - for (int dimID = 0; dimID < dim; dimID++) { - for (int i = 0; i < dim; i++) { - dshapePhysical(dofID, dimID) += dshape(dofID, i) * invJac(i, dimID); - } - } - } - - // -- Assemble the element matrix contribution = gamma * (grad(phi_i) dot grad(phi_j)) -- - elmat.SetSize(numDoF); - elmat = 0.0; - for (int i = 0; i < numDoF; i++) { - for (int j = 0; j < numDoF; j++) { - double dotProduct = 0.0; - for (int dimID = 0; dimID < dim; dimID++) { - dotProduct += dshapePhysical(i, dimID) * dshapePhysical(j, dimID); - } - elmat(i, j) += m_gamma * dotProduct; - } - } - } - - -} // namespace polyMFEMUtils diff --git a/src/poly/utils/public/integrators.h b/src/poly/utils/public/integrators.h new file mode 100644 index 0000000..960a7a5 --- /dev/null +++ b/src/poly/utils/public/integrators.h @@ -0,0 +1,85 @@ +/* *********************************************************************** +// +// Copyright (C) 2025 -- The 4D-STAR Collaboration +// File Author: Emily Boudreaux +// Last Modified: March 19, 2025 +// +// 4DSSE is free software; you can use it and/or modify +// it under the terms and restrictions the GNU General Library Public +// License version 3 (GPLv3) as published by the Free Software Foundation. +// +// 4DSSE is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU Library General Public License for more details. +// +// You should have received a copy of the GNU Library General Public License +// along with this software; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// *********************************************************************** */ +#ifndef POLYMFEMUTILS_H +#define POLYMFEMUTILS_H + +#include "mfem.hpp" +#include +#include +#include "config.h" +#include "probe.h" + +#include "quill/LogMacros.h" + + + +/** + * @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: + Config& m_config = Config::getInstance(); + Probe::LogManager& m_logManager = Probe::LogManager::getInstance(); + quill::Logger* m_logger = m_logManager.getLogger("log"); + mfem::Coefficient &m_coeff; + double m_polytropicIndex; + public: + /** + * @brief Constructor for NonlinearPowerIntegrator. + * + * @param coeff The function coefficient. + * @param n The polytropic index. + */ + NonlinearPowerIntegrator(mfem::Coefficient &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; + }; +} // namespace polyMFEMUtils + +#endif // POLYMFEMUTILS_H \ No newline at end of file diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h new file mode 100644 index 0000000..00d02e9 --- /dev/null +++ b/src/poly/utils/public/operator.h @@ -0,0 +1,45 @@ +#ifndef POLY_UTILS_OPERATOR_H +#define POLY_UTILS_OPERATOR_H + +#include "mfem.hpp" +#include + +class PolytropeOperator : public mfem::Operator { +public: + PolytropeOperator( + std::unique_ptr M, + std::unique_ptr Q, + std::unique_ptr D, + std::unique_ptr f, + const mfem::Array &blockOffsets); + ~PolytropeOperator() override = default; + + void Mult(const mfem::Vector &x, mfem::Vector &y) const override; + mfem::Operator& GetGradient(const mfem::Vector &x) const override; + + void SetEssentialTrueDofs(const mfem::Array &theta_ess_tofs, + const mfem::Array &phi_ess_tofs); + +private: + std::unique_ptr m_M; + std::unique_ptr m_Q; + std::unique_ptr m_D; + std::unique_ptr m_f; + + const mfem::Array &m_blockOffsets; + + mfem::Array m_theta_ess_tofs; + mfem::Array m_phi_ess_tofs; + + std::unique_ptr m_Mmat; + std::unique_ptr m_Qmat; + std::unique_ptr m_Dmat; + + + std::unique_ptr m_negM_op; + std::unique_ptr m_negQ_op; + mutable std::unique_ptr m_jacobian; +}; + + +#endif // POLY_UTILS_OPERATOR_H \ No newline at end of file diff --git a/src/poly/utils/public/polyIO.h b/src/poly/utils/public/polyIO.h deleted file mode 100644 index acf2237..0000000 --- a/src/poly/utils/public/polyIO.h +++ /dev/null @@ -1,36 +0,0 @@ -/* *********************************************************************** -// -// Copyright (C) 2025 -- The 4D-STAR Collaboration -// File Author: Emily Boudreaux -// Last Modified: February 14, 2025 -// -// 4DSSE is free software; you can use it and/or modify -// it under the terms and restrictions the GNU General Library Public -// License version 3 (GPLv3) as published by the Free Software Foundation. -// -// 4DSSE is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -// See the GNU Library General Public License for more details. -// -// You should have received a copy of the GNU Library General Public License -// along with this software; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// *********************************************************************** */ -#ifndef POLY_IO_H -#define POLY_IO_H - -#include "mfem.hpp" -#include - -/** - * @brief Writes the solution to a CSV file. - * - * @param u The GridFunction containing the solution. - * @param mesh The mesh associated with the solution. - * @param filename The name of the CSV file to write to. - */ -void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename); - -#endif // POLY_IO_H \ No newline at end of file diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h deleted file mode 100644 index d194d8a..0000000 --- a/src/poly/utils/public/polyMFEMUtils.h +++ /dev/null @@ -1,190 +0,0 @@ -/* *********************************************************************** -// -// Copyright (C) 2025 -- The 4D-STAR Collaboration -// File Author: Emily Boudreaux -// Last Modified: March 19, 2025 -// -// 4DSSE is free software; you can use it and/or modify -// it under the terms and restrictions the GNU General Library Public -// License version 3 (GPLv3) as published by the Free Software Foundation. -// -// 4DSSE is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -// See the GNU Library General Public License for more details. -// -// You should have received a copy of the GNU Library General Public License -// along with this software; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// *********************************************************************** */ -#ifndef POLYMFEMUTILS_H -#define POLYMFEMUTILS_H - -#include "mfem.hpp" -#include -#include -#include "config.h" -#include "probe.h" -#include - - - -/** - * @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: - Config& config = Config::getInstance(); - mfem::Coefficient &coeff_; - double polytropicIndex; - public: - /** - * @brief Constructor for NonlinearPowerIntegrator. - * - * @param coeff The function coefficient. - * @param n The polytropic index. - */ - NonlinearPowerIntegrator(mfem::Coefficient &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; - 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; - }; - - 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 m_originElementIDs; - mfem::Array m_originIntegrationPoints; - mfem::DenseMatrix m_originCoordinateMatrix; - mfem::Vector m_originCoordinateVector; - mfem::Mesh* m_mesh; - mfem::Array m_originConnectedElementIds; - public: - ConstraintIntegrator(const double gamma, mfem::Mesh* mesh); - ~ConstraintIntegrator() = default; - - void AssembleElementMatrix(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, mfem::DenseMatrix &elmat) override; - - }; - -} // namespace polyMFEMUtils - -#endif // POLYMFEMUTILS_H \ No newline at end of file