feat(poly): added skeleton of polytrope model

the polytrope module will be used as an initial guess to the solver. A skeleton of this has been imported from https://github.com/tboudreaux/FEMPolytrope

This module will need major updates still to handle 3D, proper boundary conditions, and to incorporate it with the rest of our meshing scheme

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]
+)

src/poly/coeff/private/polyCoeff.cpp

@@ -0,0 +1,18 @@
+#include "mfem.hpp"
+#include <cmath>
+
+#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);
+}

src/poly/coeff/public/polyCoeff.h

@@ -0,0 +1,8 @@
+#include "mfem.hpp"
+#include <cmath>
+
+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);

src/poly/meson.build

@@ -0,0 +1,3 @@
+subdir('coeff')
+subdir('utils')
+subdir('solver')

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]
+)

src/poly/utils/private/polyIO.cpp

@@ -0,0 +1,23 @@
+#include "mfem.hpp"
+#include <string>
+#include<fstream>
+
+#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;
+}

src/poly/utils/private/polyMFEMUtils.cpp

@@ -0,0 +1,175 @@
+#include "mfem.hpp"
+#include <string>
+#include <iostream>
+#include <cmath>
+
+#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);
+  }
+}

src/poly/utils/public/polyMFEMUtils.h

@@ -0,0 +1,43 @@
+#include "mfem.hpp"
+#include <string>
+
+
+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<mfem::NonlinearFormIntegrator*> 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;
+};