Merge pull request #12 from tboudreaux/feature/meshing

Added base MeshIO class and MFEM in tree building

.gitignore

@@ -56,4 +56,9 @@ tags
 *.log
 *.cache
 *.private
-*.private/
+*.private/
+
+subprojects/mfem/
+subprojects/tetgen/
+
+.vscode/

Readme.md

@@ -1,8 +1,16 @@
 # New implimentation of 3+1D SSE
+New (as yet unnamed) 4DSSE code. 
 We need an exciting name.
 
-# Building
-In order to build run 
+This code is very early in development and should not be used for scientific purposes yet. 
+
+
+## Building
+In order to build you will need meson installed on your system. The easiest way to do this is to use the python package manager (pip)
+```bash
+pip install meson
+```
+You can then either use the mk script or meson commands automatically. When running either the script or meson commands manually `MFEM` will be pulled from github and built. As part of this a small patch will be applied to the MFEM `CMakeLists.txt` file. This process should only need to happen once as future builds will use the cached version of MFEM in `subprojects` and the cached build files of `MFEM` in `build`.
 ```bash
 ./mk
 ```
@@ -10,3 +18,6 @@ if you want to build with no test suite run
 ```bash
 ./mk --noTest
 ```
+
+## Current Status
+Currently we are working on implimenting modules such as opacity, equation of state, polytrope, and meshing. Builds may not work on any branches at any time. 

build-config/meson.build

@@ -0,0 +1 @@
+subdir('mfem')

build-config/mfem/disable_mfem_selfcheck.patch

@@ -0,0 +1,41 @@
+--- subprojects/mfem/CMakeLists.txt	2025-02-12 15:54:52.454728232 -0500
++++ CMakeLists.txt.bak	2025-02-12 16:08:06.654542689 -0500
+@@ -765,7 +765,7 @@
+ if (MFEM_ENABLE_EXAMPLES)
+   add_subdirectory(examples) #install examples if enabled
+ else()
+-  add_subdirectory(examples EXCLUDE_FROM_ALL)
++  # add_subdirectory(examples EXCLUDE_FROM_ALL)
+ endif()
+ 
+ # Create a target for all miniapps and, optionally, enable it.
+@@ -774,7 +774,7 @@
+ if (MFEM_ENABLE_MINIAPPS)
+   add_subdirectory(miniapps) #install miniapps if enabled
+ else()
+-  add_subdirectory(miniapps EXCLUDE_FROM_ALL)
++  # add_subdirectory(miniapps EXCLUDE_FROM_ALL)
+ endif()
+ 
+ # Target to build all executables, i.e. everything.
+@@ -801,19 +801,7 @@
+   add_dependencies(${MFEM_EXEC_PREREQUISITES_TARGET_NAME} copy_data)
+ endif()
+ 
+-# Add 'check' target - quick test
+-set(MFEM_CHECK_TARGET_NAME ${MFEM_CUSTOM_TARGET_PREFIX}check)
+-if (NOT MFEM_USE_MPI)
+-  add_custom_target(${MFEM_CHECK_TARGET_NAME}
+-    ${CMAKE_CTEST_COMMAND} -R \"^ex1_ser\" -C ${CMAKE_CFG_INTDIR}
+-    USES_TERMINAL)
+-  add_dependencies(${MFEM_CHECK_TARGET_NAME} ex1)
+-else()
+-  add_custom_target(${MFEM_CHECK_TARGET_NAME}
+-    ${CMAKE_CTEST_COMMAND} -R \"^ex1p\" -C ${CMAKE_CFG_INTDIR}
+-    USES_TERMINAL)
+-  add_dependencies(${MFEM_CHECK_TARGET_NAME} ex1p)
+-endif()
++message(STATUS "MFEM Miniapps and Examples disabled by patch!")
+ 
+ #-------------------------------------------------------------------------------
+ # Documentation

build-config/mfem/meson.build

@@ -0,0 +1,24 @@
+cmake = import('cmake')
+patchFile = files('disable_mfem_selfcheck.patch')
+
+patch_check = run_command('grep', '-q', 'MFEM_CHECK_TARGET_NAME', 'subprojects/mfem/CMakeLists.txt', check: false)
+if patch_check.returncode() == 0
+    message('Patching MFEM CMakeLists.txt to remove self checks')
+    run_command('patch', '-p4', '-d', '../../subprojects/mfem', '-i', patchFile[0].full_path(), check: true)
+else
+    message('MFEM CMakeLists.txt already patched')
+endif
+
+mfem_cmake_options = cmake.subproject_options()
+mfem_cmake_options.add_cmake_defines({
+    'MFEM_ENABLE_EXAMPLES': 'OFF',
+    'MFEM_ENABLE_TESTING': 'OFF',
+    'MFEM_ENABLE_MINIAPPS': 'OFF',
+    'MFEM_USE_BENCMARK': 'OFF',
+})
+
+mfem_sp = cmake.subproject(
+    'mfem',
+    options: mfem_cmake_options)
+mfem_dep = mfem_sp.dependency('mfem')
+add_project_arguments('-I' + meson.current_build_dir() + '/subprojects/mfem/__CMake_build/config', language: 'cpp')

meson.build

@@ -3,7 +3,10 @@ project('4DSSE', 'cpp', version: '0.0.1a', default_options: ['cpp_std=c++23'], m
 # Add default visibility for all C++ targets
 add_project_arguments('-fvisibility=default', language: 'cpp')
 
+# Build external dependencies first so that all the embedded resources are available to the other targets
+subdir('build-config')
 
+# Build the main project
 subdir('src')
 if get_option('build_tests')
     subdir('tests')

mk

@@ -1,10 +1,5 @@
 #!/bin/bash
 
-# Check if the build directory is present, and remove it
-if [ -d "build" ]; then
-  rm -rf build
-fi
-
 # Check for the --noTest flag
 if [[ "$1" == "--noTest" ]]; then
   meson setup build -Dbuild_tests=false

src/meshIO/meson.build

@@ -0,0 +1,19 @@
+# Define the library
+meshIO_sources = files(
+  'private/meshIO.cpp',
+)
+
+meshIO_headers = files(
+    'public/meshIO.h'
+)
+
+# Define the libmeshIO library so it can be linked against by other parts of the build system
+libmeshIO = static_library('meshIO',
+    meshIO_sources, 
+    include_directories: include_directories('public'),
+    cpp_args: ['-fvisibility=default'],
+    dependencies: [mfem_dep],
+    install : true)
+
+# Make headers accessible
+install_headers(meshIO_headers, subdir : '4DSSE/meshIO')

src/meshIO/private/meshIO.cpp

@@ -0,0 +1,37 @@
+#include "mfem.hpp"
+#include <string>
+#include <iostream>
+#include <fstream>
+
+#include "meshIO.h"
+
+
+MeshIO::MeshIO(const std::string &mesh_file)
+{
+    mesh_file_ = mesh_file;
+    std::ifstream mesh_stream(mesh_file);
+    if (!mesh_stream)
+    {
+        throw std::runtime_error("Mesh file not found: " + mesh_file);
+        loaded_ = false;
+    }
+    else
+    {
+        mesh_ = mfem::Mesh(mesh_stream, 1, 2);
+        loaded_ = true;
+    }
+}
+
+MeshIO::~MeshIO()
+{
+}
+
+bool MeshIO::IsLoaded() const
+{
+    return loaded_;
+}
+
+mfem::Mesh& MeshIO::GetMesh()
+{
+    return mesh_;
+}

src/meshIO/public/meshIO.h

@@ -0,0 +1,42 @@
+#ifndef MESHIO_H
+#define MESHIO_H
+
+#include "mfem.hpp"
+#include <string>
+
+/**
+ * @brief Class for handling mesh input/output operations.
+ */
+class MeshIO
+{
+private:
+    bool loaded_; ///< Flag to indicate if the mesh is loaded
+    std::string mesh_file_; ///< Filename of the mesh file
+    mfem::Mesh mesh_; ///< The mesh object
+
+public:
+    /**
+     * @brief Constructor that initializes the MeshIO object with a mesh file.
+     * @param mesh_file The name of the mesh file.
+     */
+    MeshIO(const std::string &mesh_file);
+
+    /**
+     * @brief Destructor for the MeshIO class.
+     */
+    ~MeshIO();
+
+    /**
+     * @brief Get the mesh object.
+     * @return Reference to the mesh object.
+     */
+    mfem::Mesh& GetMesh();
+
+    /**
+     * @brief Check if the mesh is loaded.
+     * @return True if the mesh is loaded, false otherwise.
+     */
+    bool IsLoaded() const;
+};
+
+#endif // MESHIO_H

src/meson.build

@@ -4,4 +4,5 @@ subdir('resources')
 # Build the main source code
 subdir('dobj')
 subdir('const')
-subdir('opatIO')
+subdir('opatIO')
+subdir('meshIO')

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

subprojects/mfem.wrap

@@ -0,0 +1,5 @@
+[wrap-git]
+url = https://github.com/mfem/mfem.git
+revision = master
+
+[cmake]

tests/meshIO/meshIOTest.cpp

@@ -0,0 +1,38 @@
+#include <gtest/gtest.h>
+#include "meshIO.h"
+#include <iostream>
+#include <string>
+#include "mfem.hpp"
+
+std::string EXAMPLE_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/src/resources/mesh/sphere.msh";
+
+double ComputeMeshVolume(mfem::Mesh &mesh)
+{
+    double totalVolume = 0.0;
+    for (int i = 0; i < mesh.GetNE(); i++) // Loop over all elements
+    {
+        totalVolume += mesh.GetElementVolume(i);
+    }
+    return totalVolume;
+}
+
+
+class meshIOTest : public ::testing::Test {};
+
+TEST_F(meshIOTest, DefaultConstructor) {
+    EXPECT_NO_THROW(MeshIO meshIO(EXAMPLE_FILENAME));
+}
+
+TEST_F(meshIOTest, IsLoaded) {
+    MeshIO meshIO(EXAMPLE_FILENAME);
+    EXPECT_EQ(meshIO.IsLoaded(), true);
+}
+
+TEST_F(meshIOTest, GetMesh) {
+    MeshIO meshIO(EXAMPLE_FILENAME);
+    mfem::Mesh& mesh = meshIO.GetMesh();
+    EXPECT_EQ(mesh.GetNE(), 670);
+    EXPECT_EQ(mesh.GetNV(), 201);
+    double volume = ComputeMeshVolume(mesh);
+    EXPECT_DOUBLE_EQ(volume, 3.9357596288315868);
+}

tests/meshIO/meson.build

@@ -0,0 +1,22 @@
+# Test files for mesh
+test_sources = [
+    'meshIOTest.cpp',
+]
+
+foreach test_file : test_sources
+  exe_name = test_file.split('.')[0]
+  message('Building test: ' + exe_name)
+
+  # Create an executable target for each test
+  test_exe = executable(
+      exe_name,
+      test_file,
+      dependencies: [gtest_dep, mfem_dep],
+      include_directories: include_directories('../../src/meshIO/public'),
+      link_with: libmeshIO,  # Link the dobj library
+      install_rpath: '@loader_path/../../src'  # Ensure runtime library path resolves correctly
+  )
+
+  # Add the executable as a test
+  test(exe_name, test_exe)
+endforeach

tests/meson.build

@@ -6,6 +6,7 @@ gtest_nomain_dep = dependency('gtest', main: false, required : true)
 subdir('dobj')
 subdir('const')
 subdir('opatIO')
+subdir('meshIO')
 
 # Subdirectories for sandbox tests
 subdir('dobj_sandbox')

utils/meshGeneration/generateMesh.py

@@ -0,0 +1,24 @@
+import pygmsh
+import meshio
+
+import argparse
+
+def generate_spherical_mesh(radius=1, meshSize=0.1):
+    with pygmsh.geo.Geometry() as geo:
+        # Create a spherical (ball) geometry centered at (0,0,0)
+        sphere = geo.add_ball([0, 0, 0], radius)
+        # Generate a 2D mesh (i.e. surface mesh) of the sphere
+        myMesh = geo.generate_mesh(dim=3)
+    return myMesh
+
+def write_mfem_mesh(meshData, filename):
+    meshio.write(filename, meshData, file_format='gmsh22')
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Generate a spherical mesh')
+    parser.add_argument('--radius', type=float, default=1.0, help='Radius of the sphere')
+    parser.add_argument('--meshSize', type=float, default=0.1, help='Mesh size')
+    args = parser.parse_args()
+
+    meshData = generate_spherical_mesh(args.radius, args.meshSize)
+    write_mfem_mesh(meshData, 'sphere.msh')

utils/meshGeneration/readme.md

@@ -0,0 +1,8 @@
+# spherical mesh generation
+A simple script to generate the base spherical mesh which 4DSSE uses to solve equations. 
+This will produce a unit sphere mesh centered on (0, 0, 0) which can be rescaled within the
+4DSSE program.
+
+## Dependecies
+- pygmsh
+- meshio

utils/opatio/src/opatio/opat/opat.py

@@ -61,6 +61,46 @@ defaultHeader = Header(
 )
 
 class OpatIO:
+    """
+    @brief Class for handling OPAT file input/output operations.
+    This class provides methods to validate, manipulate, and save OPAT files. It includes functionalities to validate character arrays, 1D arrays, and 2D arrays, compute checksums, set header information, add tables, and save the OPAT file in both ASCII and binary formats.
+    Attributes:
+        header (Header): The header of the OPAT file.
+        tables (List[Tuple[Tuple[float, float], OPATTable]]): A list of tables in the OPAT file.
+    Methods:
+        validate_char_array_size(s: str, nmax: int) -> bool:
+            Validate the size of a character array.
+        validate_logKappa(logKappa):
+            Validate the logKappa array.
+        validate_1D(arr, name: str):
+            Validate a 1D array.
+        compute_checksum(data: bytes) -> bytes:
+            Compute the SHA-256 checksum of the given data.
+        set_version(version: int) -> int:
+            Set the version of the OPAT file.
+        set_source(source: str) -> str:
+            Set the source information of the OPAT file.
+        set_comment(comment: str) -> str:
+            Set the comment of the OPAT file.
+        add_table(X: float, Z: float, logR: Iterable[float], logT: Iterable[float], logKappa: Iterable[Iterable[float]]):
+            Add a table to the OPAT file.
+        _header_bytes() -> bytes:
+            Convert the header to bytes.
+        _table_bytes(table: OPATTable) -> Tuple[bytes, bytes]:
+            Convert a table to bytes.
+        _tableIndex_bytes(tableIndex: TableIndex) -> bytes:
+            Convert a table index to bytes.
+        __repr__() -> str:
+            Get the string representation of the OpatIO object.
+        _format_table_as_string(table: OPATTable, X: float, Z: float) -> str:
+            Format a table as a string.
+        print_table_indexes(table_indexes: List[TableIndex]) -> str:
+            Print the table indexes.
+        save_as_ascii(filename: str) -> str:
+            Save the OPAT file as an ASCII file.
+        save(filename: str) -> str:
+            Save the OPAT file as a binary file.
+    """
     def __init__(self):
         self.header: Header = defaultHeader
         self.tables: List[Tuple[Tuple[float, float], OPATTable]] = []