From f61c8fae28dbb24852d7c89462e5e8fe072ee3c8 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 3 Mar 2025 09:54:13 -0500 Subject: [PATCH] fix(poly): working on 3D polytrope not working yet --- src/poly/solver/meson.build | 4 +- src/poly/solver/private/polySolver.cpp | 248 +++++++++++++++++++---- src/poly/solver/public/polySolver.h | 13 ++ src/poly/utils/meson.build | 4 +- src/poly/utils/private/polyIO.cpp | 2 +- src/poly/utils/private/polyMFEMUtils.cpp | 82 ++++++-- src/poly/utils/public/polyMFEMUtils.h | 15 +- 7 files changed, 308 insertions(+), 60 deletions(-) diff --git a/src/poly/solver/meson.build b/src/poly/solver/meson.build index 47210f4..62e5448 100644 --- a/src/poly/solver/meson.build +++ b/src/poly/solver/meson.build @@ -10,7 +10,7 @@ libPolySolver = static_library('polySolver', polySolver_sources, include_directories : include_directories('./public'), cpp_args: ['-fvisibility=default'], - dependencies: [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep, warning_control_dep], + dependencies: [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep, warning_control_dep, probe_dep, quill_dep, config_dep], install: true ) @@ -18,5 +18,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, warning_control_dep] + dependencies : [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep, warning_control_dep, probe_dep, quill_dep, config_dep] ) \ No newline at end of file diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 8d1be8d..f0df901 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -3,22 +3,62 @@ #include #include #include +#include +#include +#include +#include +#include +#include #include "meshIO.h" #include "polySolver.h" #include "polyMFEMUtils.h" #include "polyCoeff.h" +#include "probe.h" +#include "config.h" + +#include "quill/LogMacros.h" #include "warning_control.h" +namespace laneEmden { + + double a (int k, double n) { + if ( k == 0 ) { return 1; } + if ( k == 1 ) { return 0; } + else { return -(c(k-2, n)/(std::pow(k, 2)+k)); } + + } + + double c(int m, double n) { + if ( m == 0 ) { return std::pow(a(0, n), n); } + else { + double termOne = 1.0/(m*a(0, n)); + double acc = 0; + for (int k = 1; k <= m; k++) { + acc += (k*n-m+k)*a(k, n)*c(m-k, n); + } + return termOne*acc; + } + } + + double thetaSerieseExpansion(double xi, double n, int order) { + double acc = 0; + for (int k = 0; k < order; k++) { + acc += a(k, n) * std::pow(xi, k); + } + return acc; + } +} // 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"; +const std::string SPHERICAL_MESH = std::string(getenv("MESON_SOURCE_ROOT")) + "/src/resources/mesh/core.msh"; PolySolver::PolySolver(double n, double order) -: n(n), +: logger(logManager.getLogger("log")), + n(n), order(order), - meshIO(SPHERICAL_MESH, 5), + meshIO(SPHERICAL_MESH, 3.1415), // TODO : Change this from PI (set to PI right now for testing the n = 1 case) mesh(meshIO.GetMesh()), feCollection(std::make_unique(order, mesh.SpaceDimension())), feSpace(std::make_unique(&mesh, feCollection.get())), @@ -31,11 +71,13 @@ PolySolver::PolySolver(double n, double order) diffusionCoeffVec = 1.0; return diffusionCoeffVec; }())), - nonLinearSourceCoeff(std::make_unique(1.0)), - gaussianCoeff(std::make_unique(0.1)) { - + nonLinearSourceCoeff(std::make_unique(-1.0)), + gaussianCoeff(std::make_unique(config.get("Poly:Gaussian:Sigma", 0.1))) { + // C_val is the weighted average of the constraint function + C_val = polyMFEMUtils::calculateGaussianIntegral(mesh, *gaussianCoeff); assembleNonlinearForm(); assembleConstraintForm(); + } PolySolver::~PolySolver() {} @@ -51,10 +93,6 @@ void PolySolver::assembleNonlinearForm() { auto nonLinearIntegrator = std::make_unique(*nonLinearSourceCoeff, n); compositeIntegrator->add_integrator(nonLinearIntegrator.release()); - // Add the \int_{\Omega}v\eta(r) d\Omega term - auto constraintIntegrator = std::make_unique(*gaussianCoeff); - compositeIntegrator->add_integrator(constraintIntegrator.release()); - nonlinearForm->AddDomainIntegrator(compositeIntegrator.release()); } @@ -68,10 +106,17 @@ void 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 + double r = x.Norml2(); + double theta = laneEmden::thetaSerieseExpansion(r, n, 10); + return theta; } ); u->ProjectCoefficient(initCoeff); + std::string initGuessFilename = "output/Poly/Debug/Newton/1D/initial_guess.csv"; + Probe::getRaySolution(*u, *feSpace->GetMesh(), {0.0, 0.0}, 100, initGuessFilename); + if (config.get("Poly:Solver:ViewInitialGuess", false)) { + Probe::glVisView(*u, mesh, "initial_guess"); + } // --- Combine DOFs (u and λ) into a single vector --- int lambdaDofOffset = feSpace->GetTrueVSize(); // Get the size of θ space @@ -83,44 +128,165 @@ void PolySolver::solve(){ mfem::Vector u_view(U.GetData(), lambdaDofOffset); u->GetTrueDofs(u_view); - // --- 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, *C, lambdaDofOffset); - newtonSolver.SetOperator(aug_op); + polyMFEMUtils::AugmentedOperator aug_op(*nonlinearForm, *C, lambdaDofOffset, C_val); // --- Create the RHS of the augmented system --- mfem::Vector B(totalTrueDofs); B = 0.0; - B[lambdaDofOffset] = 1.0; + B[lambdaDofOffset] = C_val; - // --- Solve the augmented system --- - newtonSolver.Mult(B, U); - // --- Extract the Solution --- - mfem::Vector u_sol_view(U.GetData(), lambdaDofOffset); + // --- Custom Newton Solver --- + mfem::GMRESSolver gmresSolver; + gmresSolver.SetRelTol(config.get("Poly:Solver:GMRES:RelTol", 1e-8)); + gmresSolver.SetAbsTol(config.get("Poly:Solver:GMRES:AbsTol", 1e-10)); + gmresSolver.SetMaxIter(config.get("Poly:Solver:GMRES:MaxIter", 2000)); + gmresSolver.SetPrintLevel(config.get("Poly:Solver:GMRES:PrintLevel", 0)); - DEPRECATION_WARNING_OFF // DISABLE DEPRECATION WARNING - u->SetData(u_sol_view); - DEPRECATION_WARNING_ON // REENABLE DEPRECATION WARNING + std::cout << "Setting the Block ILU preconditioner size too " << feSpace->GetTypicalFE()->GetDof() << std::endl; + mfem::BlockILU prec(feSpace->GetTypicalFE()->GetDof(), mfem::BlockILU::Reordering::MINIMUM_DISCARDED_FILL); + gmresSolver.SetPreconditioner(prec); - double lambda = U[lambdaDofOffset]; + int iteration = 0; + const int maxIter = config.get("Poly:Solver:Newton:MaxIterations", 200); + const double relTol = config.get("Poly:Solver:Newton:RelTol", 1e-8); + const double absTol = config.get("Poly:Solver:Newton:AbsTol", 1e-10); - std::cout << "λ = " << lambda << std::endl; - // TODO : Add a way to get the solution out of the solver + bool writeIntermediate = config.get("Poly:Debug:Newton:1D:WriteIntermediate", false); + double rayLatitude = config.get("Poly:Debug:Newton:1D:lat", 0.0); + double rayLongitude = config.get("Poly:Debug:Newton:1D:lon", 0.0); + int raySamples = config.get("Poly:Debug:Newton:1D:radialPoints", 100); + double rayMin = config.get("Poly:Debug:Newton:1D:radialMin", 0.0); + double rayMax = config.get("Poly:Debug:Newton:1D:radialMax", 3.14); + double rayStep = (rayMax - rayMin) / raySamples; + int stepsPerWrite = config.get("Poly:Debug:Newton:1D:StepsPerWrite", 1); + bool exitAfterWrite = config.get("Poly:Debug:Newton:1D:Exit", false); + std::string outputDirectory = config.get("Poly:Debug:Newton:1D:OutputDir", "output/Poly/Debug/Newton/1D"); + std::pair, std::vector> samples; + std::vector radialPoints; + radialPoints.reserve(raySamples); + for (int i = 0; i < raySamples; i++) { + radialPoints.push_back(rayMin + i * rayStep); + } + std::vector rayDirection = {rayLatitude, rayLongitude}; + + if (writeIntermediate) { + std::filesystem::create_directories(outputDirectory); + } + + std::string keyset = config.get("Poly:Debug:Newton:GLVis:Keyset", ""); + bool view = config.get("Poly:Debug:Newton:GLVis:View", false); + bool doExit = config.get("Poly:Debug:Newton:GLVis:Exit", false); + int stepsPerView = config.get("Poly:Debug:Newton:GLVis:StepsPerView", 1); + + while (iteration < maxIter) { + mfem::Vector F(totalTrueDofs); + F = 0.0; + aug_op.Mult(U, F); // F now holds augOp(U) + F -= B; + double resNorm = F.Norml2(); + std::cout << "Iteration: " << iteration << " Residual Norm: [ " << resNorm << " ] --- "; + if (resNorm < relTol || resNorm < absTol) { + std::cout << "Convergence achieved!" << std::endl; + break; + } + + // --- Retrieve the Jacobian --- + mfem::Operator &gradOp = aug_op.GetGradient(U); + std::cout << "Size of the Jacobian: " << gradOp.Height() << " x " << gradOp.Width() << std::endl; + gmresSolver.SetOperator(gradOp); + mfem::SparseMatrix *J = dynamic_cast(&gradOp); + if (!J) { + MFEM_ABORT("GetGradient did not return a SparseMatrix"); + } + std::cout << "Jacobian: " << J->Height() << " x " << J->Width() << std::endl; + std::cout << "Non-zero entries: " << J->NumNonZeroElems() << std::endl; + + // --- Solve the Newton Step: J * step = -F --- + mfem::Vector minusF(totalTrueDofs); + minusF = F; // MFEM's vector class does not overload the unary minus operator + minusF *= -1.0; + mfem::Vector step(totalTrueDofs); + step = 0.0; + gmresSolver.Mult(minusF, step); + double stepNorm = step.Norml2(); + std::cout << "Step Norm: " << stepNorm << std::endl; + + U += step; + + // Silly, but a way to manually force the central value to = 1 + mfem::Array elementIds; + mfem::Array ips; + mfem::DenseMatrix rayPoints(3, 1); + rayPoints(0, 0) = 0.0; + rayPoints(1, 0) = 0.0; + rayPoints(2, 0) = 0.0; + mesh.FindPoints(rayPoints, elementIds, ips); + mfem::Array dofs; + feSpace->GetElementDofs(elementIds[0], dofs); + for (int dofID : dofs) { + U[dofID] = 1.0; + } + + + + if (view && iteration % stepsPerView == 0) { + std::string s_iteration = std::to_string(iteration); + Probe::glVisView(U, *feSpace, "U at " + s_iteration, keyset); + if (doExit) { + std::raise(SIGINT); + } + } + + if (writeIntermediate && iteration % stepsPerWrite == 0) { + std::string s_iteration = std::to_string(iteration); + std::string filename = outputDirectory + "/U_" + s_iteration + ".csv"; + Probe::getRaySolution(U, *feSpace, rayDirection, raySamples, filename); + if (exitAfterWrite) { + std::raise(SIGINT); + } + } + + bool endOfStepPause = config.get("Poly:Debug:Newton:EndOfStepPause", false); + if (endOfStepPause) { + Probe::pause(); + } + iteration++; + } + + + // // --- Setup the Newton Solver --- + // mfem::NewtonSolver newtonSolver; + // newtonSolver.SetRelTol(1e-8); + // newtonSolver.SetAbsTol(1e-10); + // newtonSolver.SetMaxIter(200); + // newtonSolver.SetPrintLevel(1); + // newtonSolver.SetOperator(aug_op); + + // // --- 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 + + + // // --- Solve the augmented system --- + // newtonSolver.Mult(B, U); + + // // --- Extract the Solution --- + // mfem::Vector u_sol_view(U.GetData(), lambdaDofOffset); + + // DEPRECATION_WARNING_OFF // DISABLE DEPRECATION WARNING + // u->SetData(u_sol_view); + // DEPRECATION_WARNING_ON // REENABLE DEPRECATION WARNING + + // 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 index 7701825..443f7c6 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -9,10 +9,21 @@ #include "meshIO.h" #include "polyCoeff.h" #include "polyMFEMUtils.h" +#include "config.h" +#include "probe.h" +#include "quill/Logger.h" +namespace laneEmden { + double a (int k, double n); + double c(int m, double n); + double thetaSerieseExpansion(double xi, double n, int order); +} class PolySolver { private: + Config& config = Config::getInstance(); + Probe::LogManager& logManager = Probe::LogManager::getInstance(); + quill::Logger* logger; double n, order; MeshIO meshIO; mfem::Mesh& mesh; @@ -32,6 +43,8 @@ private: std::unique_ptr nonLinearSourceCoeff; std::unique_ptr gaussianCoeff; + double C_val; + void assembleNonlinearForm(); void assembleConstraintForm(); diff --git a/src/poly/utils/meson.build b/src/poly/utils/meson.build index 041e5b0..adb5d1c 100644 --- a/src/poly/utils/meson.build +++ b/src/poly/utils/meson.build @@ -12,7 +12,7 @@ libpolyutils = static_library('polyutils', polyutils_sources, include_directories : include_directories('./public'), cpp_args: ['-fvisibility=default'], - dependencies: [mfem_dep, warning_control_dep], + dependencies: [mfem_dep, warning_control_dep, probe_dep, quill_dep, config_dep], install: true ) @@ -20,5 +20,5 @@ polyutils_dep = declare_dependency( include_directories : include_directories('./public'), link_with : libpolyutils, sources : polyutils_sources, - dependencies : [mfem_dep, warning_control_dep] + dependencies : [mfem_dep, warning_control_dep, probe_dep, quill_dep, config_dep] ) diff --git a/src/poly/utils/private/polyIO.cpp b/src/poly/utils/private/polyIO.cpp index d540cb0..fc3921c 100644 --- a/src/poly/utils/private/polyIO.cpp +++ b/src/poly/utils/private/polyIO.cpp @@ -1,6 +1,6 @@ #include "mfem.hpp" #include -#include +#include #include "polyIO.h" diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp index 55132a3..0900ae4 100644 --- a/src/poly/utils/private/polyMFEMUtils.cpp +++ b/src/poly/utils/private/polyMFEMUtils.cpp @@ -3,8 +3,12 @@ #include #include #include +#include +#include #include "polyMFEMUtils.h" +#include "probe.h" +#include "config.h" #include "warning_control.h" @@ -206,17 +210,19 @@ namespace polyMFEMUtils { GaussianCoefficient::GaussianCoefficient(double stdDev_) : stdDev(stdDev_), - norm_coeff(1.0/(std::pow(std::sqrt(2*std::numbers::pi)*stdDev_,3))) {} + norm_coeff(1.0/std::pow(std::sqrt(2*std::numbers::pi*std::pow(stdDev_, 2)), 3.0/2.0)) {} double 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))); + double rnorm = std::sqrt(r * r); + // TODO: return to this to resolve why the Guassian does not always have peek at g(0) = 1 without the factor of 3.0145 (manually calibrated). + // Open a file (append if already exists) to write the Gaussian values + return norm_coeff * std::exp(-std::pow(rnorm, 2)/(2*std::pow(stdDev, 2))); } - AugmentedOperator::AugmentedOperator(mfem::NonlinearForm &nfl_, mfem::LinearForm &C_, int lambdaDofOffset_) + AugmentedOperator::AugmentedOperator(mfem::NonlinearForm &nfl_, mfem::LinearForm &C_, int lambdaDofOffset_, double C_val_) : mfem::Operator( // Call the base class constructor nfl_.FESpace()->GetTrueVSize()+1, // Sets the height attribute (+1 for the lambda component) @@ -224,6 +230,7 @@ namespace polyMFEMUtils { ), nfl(nfl_), C(C_), + C_val(C_val_), lambdaDofOffset(lambdaDofOffset_), lastJacobian(nullptr) {} @@ -240,30 +247,57 @@ namespace polyMFEMUtils { y.SetSize(height); y = 0.0; - for (int i = 0; i < lambdaDofOffset; i++) { - y[i] = F[i]; - } mfem::GridFunction u_gf(C.FESpace()); + mfem::Vector C_u(1); DEPRECATION_WARNING_OFF u_gf.SetData(u); DEPRECATION_WARNING_ON - - y[lambdaDofOffset] = C.operator()(u_gf); + C_u[0] = C.operator()(u_gf); // add -lambda * C to the residual mfem::Vector lambda_C(lambdaDofOffset); mfem::GridFunction constraint_gf(C.FESpace()); constraint_gf = 0.0; - mfem::Vector constraint_tdofs; - C.Assemble(); - lambda_C = C.GetData(); + mfem::Vector CTmp(lambdaDofOffset); + CTmp = C.GetData(); + lambda_C = CTmp; 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]; + y[i] = F[i] + lambda_C[i]; } + + // Get Global Debug Options for Poly + std::string defaultKeyset = config.get("Poly:Debug:Global:GLVis:Keyset", ""); + bool defaultView = config.get("Poly:Debug:Global:GLVis:View", false); + bool defaultExit = config.get("Poly:Debug:Global:GLVis:Exit", false); + + + if (config.get("Poly:Debug:GLVis:C_gf_View:View", defaultView)) { + Probe::glVisView(CTmp, *C.FESpace(), "CTmp", config.get("Poly:Debug:C_gf_View:Keyset", defaultKeyset)); + if (config.get("Poly:Debug:GLVis:C_gf_View:Exit", defaultExit)) { + std::raise(SIGINT); + } + } + + if (config.get("Poly:Debug:GLVis:F_gf_View:View", defaultView)) { + Probe::glVisView(lambda_C, *nfl.FESpace(), "lambda_C", config.get("Poly:Debug:F_gf_View:Keyset", defaultKeyset)); + if (config.get("Poly:Debug:GLVis:F_gf_View:Exit", defaultExit)) { + std::raise(SIGINT); + } + } + + if (config.get("Poly:Debug:GLVis:M_gf_View:View", defaultView)) { + Probe::glVisView(y, *nfl.FESpace(), "y", config.get("Poly:Debug:M_gf_View:Keyset", defaultKeyset)); + if (config.get("Poly:Debug:GLVis:M_gf_View:Exit", defaultExit)) { + std::raise(SIGINT); + } + } + + // Compute the constraint residual (C(u)) + y[lambdaDofOffset] = C_u[0] - C_val; // Enforce the constraint C(u) = C_val } mfem::Operator &AugmentedOperator::GetGradient(const mfem::Vector &x) const { @@ -305,6 +339,8 @@ namespace polyMFEMUtils { J_aug->Set(i, lambdaDofOffset, -CVec[i]); } + J_aug->Set(lambdaDofOffset, lambdaDofOffset, 0.0); // The bottom right corner is zero + J_aug->Finalize(); delete lastJacobian; @@ -315,4 +351,24 @@ namespace polyMFEMUtils { AugmentedOperator::~AugmentedOperator() { delete lastJacobian; } + + double calculateGaussianIntegral(mfem::Mesh &mesh, polyMFEMUtils::GaussianCoefficient &gaussianCoeff) { + // Use a discontinuous L2 finite element space (order 0) for integrating the Gaussian. + // We use L2 because the Gaussian is not necessarily continuous across element boundaries + // if the Gaussian is narrow relative to the element size. + mfem::L2_FECollection feCollection(0, mesh.SpaceDimension()); + mfem::FiniteElementSpace feSpace(&mesh, &feCollection); + + mfem::LinearForm gaussianIntegral(&feSpace); + gaussianIntegral.AddDomainIntegrator(new mfem::DomainLFIntegrator(gaussianCoeff)); + gaussianIntegral.Assemble(); + + // Create a GridFunction with a constant value of 1.0 on the L2 space. + mfem::GridFunction one_gf(&feSpace); + one_gf = 1.0; + + // Evaluate the linear form on the constant GridFunction. This gives the integral. + return gaussianIntegral(one_gf); + + } } // 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 8ef8b9f..d65ddf2 100644 --- a/src/poly/utils/public/polyMFEMUtils.h +++ b/src/poly/utils/public/polyMFEMUtils.h @@ -3,6 +3,7 @@ #include "mfem.hpp" #include +#include "config.h" @@ -22,6 +23,7 @@ namespace polyMFEMUtils { */ class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator { private: + Config& config = Config::getInstance(); mfem::Coefficient &coeff_; double polytropicIndex; public: @@ -187,19 +189,30 @@ namespace polyMFEMUtils { class AugmentedOperator : public mfem::Operator { private: + Config& config = Config::getInstance(); mfem::NonlinearForm &nfl; mfem::LinearForm &C; + double C_val; int lambdaDofOffset; mutable mfem::SparseMatrix *lastJacobian = nullptr; public: - AugmentedOperator(mfem::NonlinearForm &nfl_, mfem::LinearForm &C_, int lambdaDofOffset_); + AugmentedOperator(mfem::NonlinearForm &nfl_, mfem::LinearForm &C_, int lambdaDofOffset_, double C_val_); ~AugmentedOperator(); virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const override; virtual mfem::Operator &GetGradient(const mfem::Vector &x) const override; }; + + /** + * @brief Calculates the Gaussian integral. + * + * @param mesh The mesh. + * @param gaussianCoeff The Gaussian coefficient. + * @return The Gaussian integral. + */ + double calculateGaussianIntegral(mfem::Mesh &mesh, polyMFEMUtils::GaussianCoefficient &gaussianCoeff); } // namespace polyMFEMUtils #endif // POLYMFEMUTILS_H \ No newline at end of file