From ec37fda35b677574cf0f7fb90472226874d8937e Mon Sep 17 00:00:00 2001
From: Emily Boudreaux <thomas@boudreauxmail.com>
Date: Tue, 29 Apr 2025 15:52:56 -0400
Subject: [PATCH] feat(meshIO): framework for handling curvilinear meshes

meshIO now has a framework to set the curvature of an input mesh. This is still an early and incomplete implimentation. This has also changed the call signature a bit for GetMesh as that now returns a reference to the mesh (which is fully owned by the meshIO object)

BREAKING CHANGE: all calls to meshIO::GetMesh now must expect a mfem::Mesh&
---
 src/meshIO/private/meshIO.cpp            | 42 +++++++++++-------------
 src/meshIO/public/meshIO.h               | 18 +++++-----
 utils/meshGeneration/generateMFEMMesh.py |  0
 3 files changed, 27 insertions(+), 33 deletions(-)
 create mode 100644 utils/meshGeneration/generateMFEMMesh.py

diff --git a/src/meshIO/private/meshIO.cpp b/src/meshIO/private/meshIO.cpp
index 66d028f..57e7310 100644
--- a/src/meshIO/private/meshIO.cpp
+++ b/src/meshIO/private/meshIO.cpp
@@ -34,21 +34,21 @@ MeshIO::MeshIO(const std::string &mesh_file, double scale_factor)
     if (!mesh_stream)
     {
         throw std::runtime_error("Mesh file not found: " + mesh_file);
-        loaded_ = false;
-    }
-    else
-    {
-        mesh_ = mfem::Mesh(mesh_stream, 1, 10);
-        loaded_ = true;
-        if (scale_factor != 1.0) { LinearRescale(scale_factor); }
     }
+    m_mesh = std::make_unique<mfem::Mesh>(mesh_stream, 1, 10);
+
+    bool discont = false;
+    int space_dim = m_mesh->SpaceDimension();
+    int ordering = mfem::Ordering::byVDIM;
+    // m_mesh->SetCurvature(3, discont, space_dim, ordering);
+    m_mesh->Finalize(false, true);
+    loaded_ = true;
+    if (scale_factor != 1.0) { LinearRescale(scale_factor); }
 }
 
-MeshIO::~MeshIO()
-{
-}
+MeshIO::~MeshIO() = default;
 
-void MeshIO::LinearRescale(double scale_factor) {
+void MeshIO::LinearRescale(const double scale_factor) const {
     if (!loaded_) {
         throw std::runtime_error("Mesh not loaded before rescaling.");
     }
@@ -58,20 +58,17 @@ void MeshIO::LinearRescale(double scale_factor) {
     }
 
     // Ensure there are nodes even for linear order meshes
-    mesh_.EnsureNodes();
-    const mfem::GridFunction *nodes = mesh_.GetNodes();
-    double *node_data = nodes->GetData();  // THIS IS KEY
-
-    int data_size = nodes->Size();
+    m_mesh->EnsureNodes();
+    const mfem::GridFunction *nodes = m_mesh->GetNodes();
+    double *node_data = nodes->GetData();
+    const int data_size = nodes->Size();
 
     for (int i = 0; i < data_size; ++i) {
         node_data[i] *= scale_factor;
     }
 
     // nodes->Update(); /updates the fespace
-    mesh_.NodesUpdated();
- 
-
+    m_mesh->NodesUpdated();
 }
 
 bool MeshIO::IsLoaded() const
@@ -79,7 +76,6 @@ bool MeshIO::IsLoaded() const
     return loaded_;
 }
 
-mfem::Mesh& MeshIO::GetMesh()
-{
-    return mesh_;
-}
\ No newline at end of file
+mfem::Mesh &MeshIO::GetMesh() const noexcept {
+    return *m_mesh;
+}
diff --git a/src/meshIO/public/meshIO.h b/src/meshIO/public/meshIO.h
index f96cf63..59c4470 100644
--- a/src/meshIO/public/meshIO.h
+++ b/src/meshIO/public/meshIO.h
@@ -18,8 +18,7 @@
 //   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 //
 // *********************************************************************** */
-#ifndef MESHIO_H
-#define MESHIO_H
+#pragma once
 
 #include "mfem.hpp"
 #include <string>
@@ -32,14 +31,15 @@ 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
+    std::unique_ptr<mfem::Mesh> m_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.
+     * @param scale_factor
      */
-    MeshIO(const std::string &mesh_file, double scale_factor = 1.0);
+    explicit MeshIO(const std::string &mesh_file, double scale_factor = 1.0);
 
     /**
      * @brief Destructor for the MeshIO class.
@@ -50,19 +50,17 @@ public:
      * @brief Rescale the mesh by a linear factor.
      * @param scale_factor The factor by which to scale the mesh.
      */
-    void LinearRescale(double scale_factor);
+    void LinearRescale(double scale_factor) const;
 
     /**
      * @brief Get the mesh object.
-     * @return Reference to the mesh object.
+     * @return Pointer to the mesh object.
      */
-    mfem::Mesh& GetMesh();
+    [[nodiscard]] mfem::Mesh & GetMesh() const noexcept;
 
     /**
      * @brief Check if the mesh is loaded.
      * @return True if the mesh is loaded, false otherwise.
      */
-    bool IsLoaded() const;
+    [[nodiscard]] bool IsLoaded() const;
 };
-
-#endif // MESHIO_H
\ No newline at end of file
diff --git a/utils/meshGeneration/generateMFEMMesh.py b/utils/meshGeneration/generateMFEMMesh.py
new file mode 100644
index 0000000..e69de29