diff --git a/src/poly/coeff/meson.build b/src/poly/coeff/meson.build new file mode 100644 index 0000000..112a1ff --- /dev/null +++ b/src/poly/coeff/meson.build @@ -0,0 +1,23 @@ +polyCoeff_sources = files( + 'private/coeff.cpp' +) + +polyCoeff_headers = files( + 'public/coeff.h' +) + +libPolyCoeff = static_library('polyCoeff', + polyCoeff_sources, + include_directories : include_directories('.'), + cpp_args: ['-fvisibility=default'], + dependencies: [mfem_dep], + install: true +) + + +polyCoeff_dep = declare_dependency( + include_directories : include_directories('.'), + link_with : libPolyCoeff, + sources : polyCoeff_sources, + dependencies : [mfem_dep] +) \ No newline at end of file diff --git a/src/poly/coeff/private/polyCoeff.cpp b/src/poly/coeff/private/polyCoeff.cpp new file mode 100644 index 0000000..ad4d117 --- /dev/null +++ b/src/poly/coeff/private/polyCoeff.cpp @@ -0,0 +1,18 @@ +#include "mfem.hpp" +#include + +#include "coeff.h" + +double xi_coeff_func(const mfem::Vector &x) { + return 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); +} + +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 new file mode 100644 index 0000000..46204f1 --- /dev/null +++ b/src/poly/coeff/public/polyCoeff.h @@ -0,0 +1,8 @@ +#include "mfem.hpp" +#include + +double xi_coeff_func(const mfem::Vector &x); + +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 diff --git a/src/poly/meson.build b/src/poly/meson.build new file mode 100644 index 0000000..d9c272f --- /dev/null +++ b/src/poly/meson.build @@ -0,0 +1,3 @@ +subdir('coeff') +subdir('utils') +subdir('solver') \ No newline at end of file diff --git a/src/poly/solver/meson.build b/src/poly/solver/meson.build new file mode 100644 index 0000000..e69de29 diff --git a/src/poly/utils/meson.build b/src/poly/utils/meson.build new file mode 100644 index 0000000..fc072eb --- /dev/null +++ b/src/poly/utils/meson.build @@ -0,0 +1,24 @@ +polyutils_sources = files( + 'private/polyIO.cpp', + 'private/polyMFEMUtils.cpp' +) + +polyutils_headers = files( + 'public/polyIO.h', + 'public/polyMFEMUtils.h' +) + +libpolyutils = static_library('polyutils', + polyutils_sources, + include_directories : include_directories('.'), + cpp_args: ['-fvisibility=default'], + dependencies: [mfem_dep], + install: true +) + +libpolyutils_dep = declare_dependency( + include_directories : include_directories('.'), + link_with : libpolyutils, + sources : polyutils_sources, + dependencies : [mfem_dep] +) diff --git a/src/poly/utils/private/polyIO.cpp b/src/poly/utils/private/polyIO.cpp new file mode 100644 index 0000000..8eaf25b --- /dev/null +++ b/src/poly/utils/private/polyIO.cpp @@ -0,0 +1,23 @@ +#include "mfem.hpp" +#include +#include + +#include "io.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 new file mode 100644 index 0000000..781efab --- /dev/null +++ b/src/poly/utils/private/polyMFEMUtils.cpp @@ -0,0 +1,175 @@ +#include "mfem.hpp" +#include +#include +#include + +#include "polyMFEMUtils.h" + +NonlinearPowerIntegrator::NonlinearPowerIntegrator( + mfem::FunctionCoefficient &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); + + + // 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() { + 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); + } +} \ No newline at end of file diff --git a/src/poly/utils/public/polyIO.h b/src/poly/utils/public/polyIO.h new file mode 100644 index 0000000..e69de29 diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h new file mode 100644 index 0000000..33985aa --- /dev/null +++ b/src/poly/utils/public/polyMFEMUtils.h @@ -0,0 +1,43 @@ +#include "mfem.hpp" +#include + + +void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename); + +class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator { +private: + mfem::FunctionCoefficient coeff_; + double polytropicIndex; +public: + NonlinearPowerIntegrator(mfem::FunctionCoefficient &coeff, double n); + + virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; + virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; +}; + +class BilinearIntegratorWrapper : public mfem::NonlinearFormIntegrator +{ +private: + mfem::BilinearFormIntegrator *integrator; +public: + BilinearIntegratorWrapper(mfem::BilinearFormIntegrator *integratorInput); + + virtual ~BilinearIntegratorWrapper() ; + + virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; + virtual void AssembleElementGrad(const mfem::FiniteElement &el,mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; +}; + +class CompositeNonlinearIntegrator: public mfem::NonlinearFormIntegrator { + private: + std::vector integrators; + public: + CompositeNonlinearIntegrator(); + + virtual ~CompositeNonlinearIntegrator(); + + void add_integrator(mfem::NonlinearFormIntegrator *integrator); + + virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; + virtual void AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; +}; \ No newline at end of file