Merge branch 'feature/polytrope' into feature/meshing

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,40 @@
+#include "mfem.hpp"
+#include <cmath>
+
+#include "coeff.h"
+
+/**
+ * @brief Computes the xi coefficient function.
+ * 
+ * @param x Input vector.
+ * @return double The computed xi coefficient.
+ */
+double xi_coeff_func(const mfem::Vector &x)
+{
+  return std::pow(x(0), 2);
+}
+
+/**
+ * @brief Computes the vector xi coefficient function.
+ * 
+ * @param x Input vector.
+ * @param v Output vector to store the computed xi coefficient.
+ */
+void vec_xi_coeff_func(const mfem::Vector &x, mfem::Vector &v)
+{
+  v.SetSize(1);
+  v[0] = -std::pow(x(0), 2);
+}
+
+/**
+ * @brief Computes the initial guess for theta.
+ * 
+ * @param x Input vector.
+ * @param root Root value used in the computation.
+ * @return double The initial guess for theta.
+ */
+double theta_initial_guess(const mfem::Vector &x, double root)
+{
+  double xi = x[0];
+  return 1 - std::pow(xi / root, 2);
+}

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 "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;
+}

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/polyIO.h

@@ -0,0 +1,16 @@
+#ifndef POLY_IO_H
+#define POLY_IO_H
+
+#include "mfem.hpp"
+#include <string>
+
+/**
+ * @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

src/poly/utils/public/polyMFEMUtils.h

@@ -0,0 +1,128 @@
+#include "mfem.hpp"
+#include <string>
+
+
+void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename);
+
+/**
+ * @brief A class for nonlinear power integrator.
+ */
+class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator {
+private: 
+  mfem::FunctionCoefficient coeff_;
+  double polytropicIndex;
+public:
+  /**
+   * @brief Constructor for NonlinearPowerIntegrator.
+   * 
+   * @param coeff The function coefficient.
+   * @param n The polytropic index.
+   */
+  NonlinearPowerIntegrator(mfem::FunctionCoefficient &coeff, double n);
+
+  /**
+   * @brief Assembles the element vector.
+   * 
+   * @param el The finite element.
+   * @param Trans The element transformation.
+   * @param elfun The element function.
+   * @param elvect The element vector to be assembled.
+   */
+  virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override;
+
+  /**
+   * @brief Assembles the element gradient.
+   * 
+   * @param el The finite element.
+   * @param Trans The element transformation.
+   * @param elfun The element function.
+   * @param elmat The element matrix to be assembled.
+   */
+  virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override;
+};
+
+/**
+ * @brief A wrapper class for bilinear integrator.
+ */
+class BilinearIntegratorWrapper : public mfem::NonlinearFormIntegrator
+{
+private:
+   mfem::BilinearFormIntegrator *integrator;
+public:
+   /**
+    * @brief Constructor for BilinearIntegratorWrapper.
+    * 
+    * @param integratorInput The bilinear form integrator input.
+    */
+   BilinearIntegratorWrapper(mfem::BilinearFormIntegrator *integratorInput);
+
+   /**
+    * @brief Destructor for BilinearIntegratorWrapper.
+    */
+   virtual ~BilinearIntegratorWrapper();
+
+   /**
+    * @brief Assembles the element vector.
+    * 
+    * @param el The finite element.
+    * @param Trans The element transformation.
+    * @param elfun The element function.
+    * @param elvect The element vector to be assembled.
+    */
+   virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override;
+
+   /**
+    * @brief Assembles the element gradient.
+    * 
+    * @param el The finite element.
+    * @param Trans The element transformation.
+    * @param elfun The element function.
+    * @param elmat The element matrix to be assembled.
+    */
+   virtual void AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override;
+};
+
+/**
+ * @brief A class for composite nonlinear integrator.
+ */
+class CompositeNonlinearIntegrator: public mfem::NonlinearFormIntegrator {
+  private:
+    std::vector<mfem::NonlinearFormIntegrator*> 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;
+};