Merge pull request #3 from tboudreaux/feat/multi_block

Vacuum coordinate and multi block

Doxyfile

@@ -48,7 +48,7 @@ PROJECT_NAME           = stroid
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = v0.2.1
+PROJECT_NUMBER         = v0.5.0
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewers a

build-config/meson.build

@@ -2,3 +2,19 @@ subdir('mfem')
 subdir('libconfig')
 subdir('CLI11')
 subdir('magic_enum')
+
+if get_option('build_python')
+    subdir('python')
+    subdir('pybind')
+endif
+
+if get_option('build_python')
+    stroid_pkg_dir    = py_installation.get_install_dir() / 'stroid'
+    stroid_includedir = stroid_pkg_dir / 'include'
+    stroid_libdir     = stroid_pkg_dir / 'lib'
+    stroid_pcdir      = stroid_libdir / 'pkgconfig'
+else
+    stroid_includedir = get_option('includedir')
+    stroid_libdir     = get_option('libdir')
+    stroid_pcdir      = get_option('libdir') / 'pkgconfig'
+endif

build-config/mfem/meson.build

@@ -1,3 +1,6 @@
+mfem_dep = dependency('mfem', required : false)
+
+if not mfem_dep.found()
     cmake = import('cmake')
     mfem_cmake_options = cmake.subproject_options()
     mfem_cmake_options.add_cmake_defines({
@@ -14,3 +17,6 @@ mfem_sp = cmake.subproject(
         'mfem',
         options: mfem_cmake_options)
     mfem_dep = mfem_sp.dependency('mfem')
+else
+    message('Using system-installed MFEM library')
+endif

build-config/pybind/meson.build

@@ -0,0 +1,3 @@
+pybind11_proj = subproject('pybind11')
+pybind11_dep = pybind11_proj.get_variable('pybind11_dep')
+python3_dep = dependency('python3')

build-config/python/meson.build

@@ -0,0 +1,5 @@
+py_installation = import('python').find_installation('python3', pure: false)
+
+py_dep = py_installation.dependency()
+py_module_prefix = ''
+py_module_suffix = 'so'

build-python/meson.build

@@ -0,0 +1,43 @@
+if get_option('build_python')
+    message('Building Python bindings...')
+
+    stroid_py_deps = [
+        py_dep,
+        pybind11_dep,
+        stroid_dep
+    ]
+
+    if host_machine.system() == 'darwin'
+        stroid_ext_rpath = '@loader_path/lib'
+    else
+        stroid_ext_rpath = '$ORIGIN/lib'
+    endif
+
+    py_sources =  [
+        meson.project_source_root() + '/src/python/bindings.cpp',
+        meson.project_source_root() + '/src/python/config/bindings.cpp',
+        meson.project_source_root() + '/src/python/exceptions/bindings.cpp',
+        meson.project_source_root() + '/src/python/IO/bindings.cpp',
+        meson.project_source_root() + '/src/python/refinement/bindings.cpp',
+        meson.project_source_root() + '/src/python/utils/bindings.cpp',
+    ]
+
+    py_mod = py_installation.extension_module(
+        '_stroid',
+        sources: py_sources,
+        dependencies: stroid_py_deps,
+        install: true,
+        link_args: stroid_ext_rpath_args,
+        build_rpath: stroid_ext_rpath,
+        install_rpath: stroid_ext_rpath,
+        subdir: 'stroid',
+    )
+
+    py_installation.install_sources(
+        meson.project_source_root() + '/src/python/stroid/__init__.py',
+        subdir: 'stroid',
+    )
+
+else
+    message('Python bindings disabled')
+endif

configs/conditioned_core.toml

@@ -0,0 +1,21 @@
+[main]
+core_mapping = "multi_block"
+refinement_levels = 2
+order = 3
+include_external_domain = true
+r_core = 0.25
+r_star = 1.0
+r_infinity = 6.0
+flattening = 0.0
+r_instability = 1e-14
+core_steepness = 1.0
+continuity_order = 2
+surface_bdr_id = 1
+inf_bdr_id = 2
+core_id = 1
+envelope_id = 2
+vacuum_id = 3
+
+[main.optimization_methods]
+tmop = false
+smoothstep = true

configs/default_config.toml

@@ -13,3 +13,6 @@ surface_bdr_id = 1
 core_id = 1
 envelope_id = 2
 vacuum_id = 3
+
+[main.optimization_methods]
+smoothstep = true

configs/test_external_domain.toml

@@ -13,3 +13,5 @@ surface_bdr_id = 1
 core_id = 1
 envelope_id = 2
 vacuum_id = 3
+[main.optimization_methods]
+smoothstep = true

configs/test_external_domain_refinement_l1.toml

@@ -14,3 +14,5 @@ core_id = 1
 envelope_id = 2
 vacuum_id = 3
 
+[main.optimization_methods]
+smoothstep = true

configs/test_flattening.toml

@@ -13,3 +13,5 @@ surface_bdr_id = 1
 core_id = 1
 envelope_id = 2
 vacuum_id = 3
+[main.optimization_methods]
+smoothstep = true

configs/test_polynomial_projection.toml

@@ -0,0 +1,17 @@
+[main]
+core_steepness = 1.0
+flattening = 0.2
+include_external_domain = false
+inf_bdr_id = 2
+order = 3
+r_core = 1.5
+r_infinity = 6.0
+r_instability = 1e-14
+r_star = 5.0
+refinement_levels = 1
+surface_bdr_id = 1
+core_id = 1
+envelope_id = 2
+vacuum_id = 3
+[main.optimization_methods]
+smoothstep = true

configs/test_refinement_l2.toml

@@ -14,3 +14,5 @@ core_id = 1
 envelope_id = 2
 vacuum_id = 3
 
+[main.optimization_methods]
+smoothstep = true

configs/test_volume_no_external.toml

@@ -14,3 +14,5 @@ core_id = 1
 envelope_id = 2
 vacuum_id = 3
 
+[main.optimization_methods]
+smoothstep = true

configs/test_volume_spherical_no_external.toml

@@ -14,3 +14,5 @@ core_id = 1
 envelope_id = 2
 vacuum_id = 3
 
+[main.optimization_methods]
+smoothstep = true

configs/test_volume_spherical_with_external.toml

@@ -14,3 +14,5 @@ core_id = 1
 envelope_id = 2
 vacuum_id = 3
 
+[main.optimization_methods]
+smoothstep = true

configs/test_volume_with_external.toml

@@ -14,3 +14,5 @@ core_id = 1
 envelope_id = 2
 vacuum_id = 3
 
+[main.optimization_methods]
+smoothstep = true

meson.build

@@ -1,4 +1,4 @@
-project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.2.1', default_options : ['cpp_std=c++23'])
+project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.5.0', default_options : ['cpp_std=c++23'])
 
 subdir('build-check')
 
@@ -13,6 +13,10 @@ if get_option('build_tools')
     subdir('tools')
 endif
 
+if get_option('build_python')
+    subdir('build-python')
+endif
+
 if get_option('pkg_config')
     pkg = import('pkgconfig')
     pkg.generate(
@@ -20,10 +24,11 @@ if get_option('pkg_config')
         description: 'Stroid multi-block curvilinear mesh generation library',
         version: meson.project_version(),
         libraries: [
-            stroid_lib
+            libstroid
         ],
         subdirs: ['stroid'],
         filebase: 'stroid',
-        install_dir: join_paths(get_option('libdir'), 'pkgconfig')
+        install_dir: join_paths(get_option('libdir'), 'pkgconfig'),
+        requires: ['fourdst_config']
     )
 endif

meson_options.txt

@@ -1,3 +1,4 @@
 option('pkg_config', type: 'boolean', value: false, description: 'generate pkg-config file for stroid')
 option('build_tests', type: 'boolean', value: true, description: 'compile subproject tests')
 option('build_tools', type: 'boolean', value: true, description: 'compile stroid command line tools')
+option('build_python', type: 'boolean', value: true, description: 'compile stroid python bindings')

pyproject.toml

@@ -0,0 +1,25 @@
+[build-system]
+requires = ["meson-python>=0.19.0", "meson>=1.9.1", "pybind11==3.0.0", "fourdst==0.10.6"]
+build-backend = "mesonpy"
+
+[project]
+name = "stroid"
+dynamic = ["version"]
+description = "O-grid mesh generation with multiple domains"
+readme = "README.md"
+license = { file = "LICENSE.txt" }
+
+authors = [
+    {name = "Emily M. Boudreaux", email = "emily@boudreauxmail.com"},
+]
+maintainers = [
+    {name = "Emily M. Boudreaux", email = "emily@boudreauxmail.com"}
+]
+
+[tool.meson-python.args]
+setup = [
+    '-Dbuild_tools=false',
+    '-Dbuild_tests=false',
+    '-Dpkg_config=false'
+]
+install = ['--skip-subprojects']

readme.md

@@ -92,11 +92,17 @@ inf_bdr_id = 2
 core_id = 1
 envelope_id = 2
 vacuum_id = 3
+core_mapping = "multi_block"
+
+
+[main.optimization_methods]
+tmop = false
+smoothstep = true
 ```
 
 <!-- Table of what these parameters do -->
 | Parameter                       | Description                                                                                                                                                                                                                                        | Default       |
-|-------------------------|-----------------------------------------------------------------------------------------------------|---------|
+|---------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|---------------|
 | refinement_levels               | Number of uniform refinement levels to apply to the mesh after generation                                                                                                                                                                          | 4             |
 | order                           | The polynomial order of the finite elements in the mesh                                                                                                                                                                                            | 3             |
 | include_external_domain         | Whether to include an external domain extending to r_infinity                                                                                                                                                                                      | true          |
@@ -111,11 +117,37 @@ vacuum_id = 3
 | core_id                         | The material ID to assign to the core region of the star                                                                                                                                                                                           | 1             |
 | envelope_id                     | The material ID to assign to the envelope region of the star                                                                                                                                                                                       | 2             |
 | vacuum_id                       | The material ID to assign to the vacuum region of the star (if included)                                                                                                                                                                           | 3             |
+| optimization_methods.tmop       | The tmop flag enables or disables the use of TMOP ideal shape unit size metric optimization during mesh generation. This can help improve the quality of the generated mesh, but will dramatically increase the time required for mesh generation. | false         |
+| optimization_methods.smoothstep | The smoothstep flag enables or disables the use of a smoothstep function to transition between the core and envelope regions of the star. This can help improve the quality of the generated mesh                                                  | true          |
+ | core_mapping                    | The core mapping strategy to use for the mesh generation. Options are "spherified" (legacy) or "multi_block" (conditioned). The multi_block strategy is strongly preferred for its improved condition number.                                      | "multi_block" |
 
 
 If no configuration file is provided, stroid will use the default parameters listed above. Further, configuration files
 need only include parameters that differ from the defaults, any parameters not specified will use the default values.
 
+### Conditioned core mapping
+
+There are two core mapping strategies, spherified and multi_block. Generally multi_block should be strongly preferred. The
+`core_mapping = "multi_block"` strategy avoids the radial rank loss at the eight corners of the spherified core
+block. It uses a Cartesian center plus six transition blocks inside the core. The inner cube has circumscribed radius
+`r_core / 2`; its six faces connect linearly to the existing spherical `r_core` interface. If enabled, spheroidal flattening is
+applied afterwards. 
+
+```python
+cfg = stroid.config.MeshConfig(core_mapping="multi_block", refinement_levels=2)
+cfg.optimization_methods = stroid.config.OptimizationMethods(tmop=False)
+mesh = stroid.GenerateMesh(cfg)
+```
+
+The optional, non-installed `geometry_quality_experiment` target may be used to measure the actual high-order geometry
+at quadrature points, vertices, edges, and near-corner probes. You may build and run it explicitly:
+
+```bash
+meson compile -C build geometry_quality_experiment
+build/tools/geometry_quality_experiment --orders 4 --refinements 2 \
+    --contraction-probe --probe-order 3 --output core_comparison.csv
+```
+
 ### C++ Interface
 Stroid can be used as a library in C++ projects. After installation, include the stroid header and link against the stroid library.
 
@@ -138,6 +170,7 @@ int main() {
     stroid::topology::Finalize(*mesh, cfg);
     stroid::topology::PromoteToHighOrder(*mesh, cfg);
     stroid::topology::ProjectMesh(*mesh, cfg);
+    stroid::topology::OptimizeMesh(*mesh, cfg);
 
 
     stroid::IO::ViewMesh(*mesh, "Spheroidal Mesh", stroid::IO::VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID);
@@ -146,7 +179,14 @@ int main() {
 
 ## Example Meshes
 An example mesh with the default configuration parameters is shown below (coloration indicates attribute IDs of different regions):
-![Example Mesh](assets/imgs/ExampleMesh.png)
+![Example Mesh](assets/imgs/ExampleMesh_multi-block.png)
+
+The legacy spherified core mapping strategy is shown below as well
+![Example Spheried Mesh](assets/imgs/ExampleMesh_spherified.png)
+
+Note that both of these meshes are shown with 3 levels of refinement and polynomial order 3. Blue shows the stellar
+domain while purple shows the vacuum domain.
+
 
 ## Funding
 Stroid is developed as part of the 4D-STAR project.

src/include/stroid/IO/mesh.h

@@ -1,7 +1,12 @@
 #pragma once
 #include <string>
+#include <expected>
+#include <istream>
+
 #include "mfem.hpp"
 
+#include "stroid/utils/types.h"
+
 namespace stroid::IO {
     /**
      * @brief Visualization modes for GLVis display.
@@ -15,18 +20,43 @@ namespace stroid::IO {
         BOUNDARY_ELEMENT_ID
     };
 
+    void SaveStroidMesh(const StroidMesh& mesh, const std::string& filename, const std::string& comment="");
+
     /**
      * @brief Save a mesh to MFEM's native `.mesh` format.
      * @param mesh Mesh to serialize.
      * @param filename Output path (including extension).
      */
     void SaveMesh(const mfem::Mesh& mesh, const std::string& filename);
+
+    /**
+     * @brief Overload of SaveMesh which accepts a StroidMesh type and will internally unpack it
+     * @param mesh StroidMesh to serialize.
+     * @param filename Path to save to
+     *
+     * @note This function is a utility wrapper to save a StroidMesh object in MFEM's native .mesh format. Data other than the mesh pointer
+     * in StroidMesh **will not be saved** (e.g. the reference mesh, the number of refinement levels, etc..). If you need to serialize an
+     * entire StroidMesh then please use the stroid::IO::SaveStroidMesh function
+     */
+    void SaveMesh(const stroid::StroidMesh& mesh, const std::string& filename);
     /**
      * @brief Save a mesh as a ParaView VTU dataset.
      * @param mesh Mesh to export.
      * @param exportName Output base name (ParaView will add extensions).
      */
     void SaveVTU(mfem::Mesh& mesh, const std::string& exportName);
+
+    /**
+     * @brief Overload of SaveVTU which accepts a StroidMesh type and will internally unpack it
+     * @param mesh StroidMesh to serialize.
+     * @param filename Path to save to
+     *
+     * @note This function is a utility wrapper to save a StroidMesh object in MFEM's native .mesh format. Data other than the mesh pointer
+     * in StroidMesh **will not be saved** (e.g. the reference mesh, the number of refinement levels, etc..). If you need to serialize an
+     * entire StroidMesh then please use the stroid::IO::SaveStroidVTU function
+     */
+    void SaveVTU(const stroid::StroidMesh& mesh, const std::string& exportName);
+
     /**
      * @brief Stream a mesh to a running GLVis server for interactive viewing.
      * @param mesh Mesh to display.
@@ -36,9 +66,22 @@ namespace stroid::IO {
      * @param visport GLVis server port.
      */
     void ViewMesh(mfem::Mesh &mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport);
+
+    void ViewMesh(const stroid::StroidMesh& mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport);
+
     /**
      * @brief Visualize boundary face valence (1=surface, 2=internal).
      * @param mesh Mesh whose boundary faces are inspected.
      */
-    void VisualizeFaceValence(mfem::Mesh& mesh);
+    void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport);
+
+    void VisualizeFaceValence(const stroid::StroidMesh& mesh, const std::string &vishost, int visport);
+
+    std::expected<StroidMesh, std::string> ParseStroidMesh(std::istream& is);
+    std::expected<StroidMesh, std::string> LoadStroidMesh(const std::string& filename);
+
+#ifdef MFEM_USE_MPI
+    std::expected<StroidMesh, std::string> ParseStroidMesh(std::istream& is, MPI_Comm comm);
+    std::expected<StroidMesh, std::string> LoadStroidMesh(const std::string& filename, MPI_Comm comm);
+#endif
 }

src/include/stroid/config/config.h

@@ -1,6 +1,18 @@
 #pragma once
 
+#include <cstdint>
+#include <optional>
+#include <string>
+#include <format>
+#include <sstream>
+
 namespace stroid::config {
+
+    struct OptimizationMethods {
+        std::optional<bool> tmop{false};
+        std::optional<bool> smoothstep{true};
+    };
+
     /**
      * @brief Configuration parameters for stroid mesh generation.
      *
@@ -15,92 +27,157 @@ namespace stroid::config {
          * @section toml
          * - [main].refinement_levels
          */
-        int refinement_levels = 4;
+        std::optional<int> refinement_levels = 4;
         /**
          * @brief Polynomial order for high-order elements.
          * @section toml
          * - [main].order
          */
-        int order = 3;
+        std::optional<int> order = 3;
         /**
          * @brief Whether to include an external domain extending to `r_infinity`.
          * @section toml
          * - [main].include_external_domain
          */
-        bool include_external_domain = true;
+        std::optional<bool> include_external_domain = true;
 
         /**
          * @brief Radius of the stellar core region.
          * @section toml
          * - [main].r_core
          */
-        double r_core = 1.5;
+        std::optional<double> r_core = 0.25;
         /**
          * @brief Radius of the stellar surface.
          * @section toml
          * - [main].r_star
          */
-        double r_star = 5.0;
+        std::optional<double> r_star = 1.0;
         /**
          * @brief Flattening factor for spheroidal shaping (0 = spherical, >0 = oblate).
          * @section toml
          * - [main].flattening
          */
-        double flattening = 0;
+        std::optional<double> flattening = 0;
 
         /**
          * @brief Outer radius of the external domain when enabled.
          * @section toml
          * - [main].r_infinity
          */
-        double r_infinity = 6.0;
+        std::optional<double> r_infinity = 6.0;
 
         /**
          * @brief Radius inside which transformations are skipped to avoid singularities.
          * @section toml
          * - [main].r_instability
          */
-        double r_instability = 1e-14;
+        std::optional<double> r_instability = 1e-14;
         /**
          * @brief Controls the smoothness/steepness of the core-to-envelope transition.
          * @section toml
          * - [main].core_steepness
          */
-        double core_steepness = 1.0;
+        std::optional<double> core_steepness = 1.0;
+
+        /**
+         * @brief Continuity order for the core-envelope transition (0 = discontinuous, 1 = C1, 2 = C2).
+         * @section toml
+         * - [main].continuity_order
+         */
+        std::optional<size_t> continuity_order = 2;
 
         /**
          * @brief Boundary attribute id for stellar surface
          * @section toml
          * - [main].surface_bdr_id
          */
-        size_t surface_bdr_id = 1;
+        std::optional<size_t> surface_bdr_id = 1;
 
         /**
          * @brief Boundary attribute id for infinity in kelvin mapping
          * @section toml
          * - [main].inf_bdr_id
          */
-        size_t inf_bdr_id = 2;
+        std::optional<size_t> inf_bdr_id = 2;
 
         /**
          * @brief Material attribute id for the core region
          * @section toml
          * - [main].core_id
          */
-        size_t core_id = 1;
+        std::optional<size_t> core_id = 1;
 
         /**
          * @brief Material attribute id for the envelope region
          * @section toml
          * - [main].envelope_id
          */
-        size_t envelope_id = 2;
+        std::optional<size_t> envelope_id = 2;
 
         /**
          * @brief Material attribute id for the external domain (if enabled)
          * @section toml
          * - [main].vacuum_id
          */
-        size_t vacuum_id = 3;
+        std::optional<size_t> vacuum_id = 3;
+
+        std::optional<OptimizationMethods> optimization_methods = OptimizationMethods{true, true};
+
+        /**
+         * @brief Core mapping strategy: legacy "spherified" or conditioned "multi_block".
+         *
+         * spherified generates a either two or three inscribed cubes then projects them into spheres.
+         * multi_block generates a multi-block topology with a single core block and six envelope blocks, then projects the core block into a sphere and the envelope blocks into a spheroid.
+         *
+         * multi_block is strongly preferred for its ~1000x improved condition number, Spherified is only provided for legacy compatibility.
+         *
+         * @section toml
+         * - [main].core_mapping
+         */
+        std::optional<std::string> core_mapping = "multi_block";
+
     };
+
+    inline std::string to_string(const MeshConfig &mesh_config) {
+        auto opt_2_string = [](const OptimizationMethods& opt) {
+            std::stringstream ss;
+            ss << "<OptimizationMethods:";
+            if (*opt.tmop) {
+                ss << " tmop";
+            }
+            if (*opt.smoothstep) {
+                ss << " smoothstep";
+            }
+            ss << ">";
+            return ss.str();
+        };
+
+        std::stringstream ss;
+
+        OptimizationMethods opt = mesh_config.optimization_methods.value_or(OptimizationMethods{false, true});
+        std::string opt_string = opt_2_string(opt);
+
+        ss << "MeshConfig:\n";
+        ss << std::format("  refinement_levels: {}\n", mesh_config.refinement_levels.value_or(4));
+        ss << std::format("  order: {}\n", mesh_config.order.value_or(3));
+        ss << std::format("  include_external_domain: {}\n", mesh_config.include_external_domain.value_or(true));
+        ss << std::format("  r_core: {}\n", mesh_config.r_core.value_or(0.25));
+        ss << std::format("  r_star: {}\n", mesh_config.r_star.value_or(1.0));
+        ss << std::format("  flattening: {}\n", mesh_config.flattening.value_or(0.0));
+        ss << std::format("  r_infinity: {}\n", mesh_config.r_infinity.value_or(6.0));
+        ss << std::format("  r_instability: {}\n", mesh_config.r_instability.value_or(1e-14));
+        ss << std::format("  core_steepness: {}\n", mesh_config.core_steepness.value_or(1.0));
+        ss << std::format("  continuity_order: {}\n", mesh_config.continuity_order.value_or(2));
+        ss << std::format("  surface_bdr_id: {}\n", mesh_config.surface_bdr_id.value_or(1));
+        ss << std::format("  inf_bdr_id: {}\n", mesh_config.inf_bdr_id.value_or(2));
+        ss << std::format("  core_id: {}\n", mesh_config.core_id.value_or(1));
+        ss << std::format("  envelope_id: {}\n", mesh_config.envelope_id.value_or(2));
+        ss << std::format("  vacuum_id: {}\n", mesh_config.vacuum_id.value_or(3));
+        ss << std::format("  optimization_methods: {}\n", opt_string);
+        ss << std::format("  core_mapping: {}\n", mesh_config.core_mapping.value_or("spherified"));
+
+        return ss.str();
+
+    }
 }

src/include/stroid/exceptions/exceptions.h

@@ -0,0 +1,3 @@
+#pragma once
+
+#include "stroid/exceptions/stroid_error.h"

src/include/stroid/exceptions/stroid_error.h

@@ -0,0 +1,25 @@
+#pragma once
+#include <exception>
+#include <string>
+
+namespace stroid::exceptions {
+    class StroidError : public std::exception {
+    public:
+        explicit StroidError(std::string message) : m_msg(std::move(message)) {}
+        const char* what() const noexcept override { return m_msg.c_str(); }
+    private:
+        std::string m_msg;
+    };
+
+    class StroidContinuityError : public StroidError {
+        using StroidError::StroidError;
+    };
+
+    class StroidMeshError : public StroidError {
+        using StroidError::StroidError;
+    };
+
+    class StroidMissingReferenceMesh : public StroidMeshError {
+        using StroidMeshError::StroidMeshError;
+    };
+}

src/include/stroid/meson.build

@@ -21,8 +21,8 @@ config.set('STROID_VERSION_PATCH', ver_parts[2])
 config.set('STROID_VERSION_TAG',   ver_parts[3])
 
 configure_file(
-    input : 'stroid.h.in',
-    output : 'stroid.h',
+    input : 'version.h.in',
+    output : 'version.h',
     configuration : config ,
     install: true,
     install_dir: get_option('includedir') / 'stroid'

src/include/stroid/refinement/uniform.h

@@ -0,0 +1,7 @@
+#pragma once
+
+#include "stroid/utils/types.h"
+
+namespace stroid::refinement {
+    void UniformRefinement(StroidMesh& mesh, size_t levels);
+}

src/include/stroid/stroid.h

@@ -4,8 +4,13 @@
 #include "stroid/topology/topology.h"
 #include "stroid/topology/mapping.h"
 #include "stroid/topology/curvilinear.h"
+#include "stroid/topology/optimize.h"
 #include "stroid/utils/mesh_utils.h"
 #include "stroid/IO/mesh.h"
+#include "stroid/utils/types.h"
+#include "stroid/refinement/uniform.h"
+#include "stroid/utils/mesh_stats.h"
+#include "stroid/version.h"
 
 /**
  * @namespace stroid
@@ -44,46 +49,38 @@
  * @endcode
  */
 namespace stroid {
-    /**
-     * @brief Version helpers for the stroid library.
-     */
-    struct version {
-        static constexpr int major = @STROID_VERSION_MAJOR@;
-        static constexpr int minor = @STROID_VERSION_MINOR@;
-        static constexpr int patch = @STROID_VERSION_PATCH@;
-        static constexpr const char* tag = "@STROID_VERSION_TAG@";
+    inline StroidMesh GenerateMesh(const fourdst::config::Config<stroid::config::MeshConfig>& cfg) {
+        StroidMesh sm;
+        sm.type = MFEM_MESH_TYPE::SERIAL;
+        sm.config = *cfg;
+        auto reference = stroid::topology::BuildSkeleton(cfg);
+        stroid::topology::Finalize(*reference, cfg);
+        sm.refinement_levels = cfg->refinement_levels.value_or(0);
 
-        static std::string toString() {
-            std::string versionStr = std::to_string(major) + "." +
-                                     std::to_string(minor) + "." +
-                                     std::to_string(patch);
-            if (std::string(tag) != "") {
-                versionStr += "-" + std::string(tag);
+        sm.reference_mesh = std::move(reference);
+        sm.mesh = utils::BuildProjected(*sm.reference_mesh, cfg);
+        if (cfg->optimization_methods.has_value() && cfg->optimization_methods.value().tmop.has_value() && cfg->optimization_methods.value().tmop.value()) {
+            stroid::topology::ApplyTMOP(*sm.mesh, cfg);
         }
-            return versionStr;
-        }
-
-        friend std::ostream& operator<<(std::ostream& os, const version&) {
-            os << toString();
-            return os;
+        sm.exterior_coordinate = stroid::topology::BuildExteriorCoordinate(*sm.mesh, *sm.reference_mesh, cfg);
+        return sm;
     }
+    inline StroidMesh GenerateMesh(const stroid::config::MeshConfig& config) {
+        fourdst::config::Config<config::MeshConfig> cfg;
+        auto Mutator = [&config](config::MeshConfig& orig) {
+            orig = config;
         };
+
+        cfg.mutate(Mutator);
+        return GenerateMesh(cfg);
+    }
+    inline StroidMesh GenerateMesh(const std::string& filename) {
+        fourdst::config::Config<stroid::config::MeshConfig> config;
+        config.load(filename);
+        return GenerateMesh(config);
+    }
 }
 
-/**
- * @namespace std
- * @brief Standard library extensions used by stroid.
- *
- * Provides a `std::formatter` specialization for `stroid::version` so it can
- * be used with `std::format` and related APIs.
- */
-// Overload format struct
-template <>
-struct std::formatter<stroid::version> : std::formatter<std::string> {
-    auto format(const stroid::version& v, auto& ctx) {
-        return std::formatter<std::string>::format(stroid::version::toString(), ctx);
-    }
-};
 
 /**
  * @namespace stroid::config

src/include/stroid/topology/curvilinear.h

@@ -3,6 +3,7 @@
 #include "mfem.hpp"
 #include "stroid/config/config.h"
 #include "fourdst/config/config.h"
+#include "stroid/utils/types.h"
 
 namespace stroid::topology {
     /**
@@ -18,4 +19,17 @@ namespace stroid::topology {
      * @param config Mesh configuration (uses radii, flattening, and mapping parameters).
      */
     void ProjectMesh(mfem::Mesh& mesh, const fourdst::config::Config<config::MeshConfig> &config);
+
+    /**
+     * @brief Build a scalar grid function representing the compactification coordinate for a mesh. This ranges from 0-1 with 0 at the stellar surface and 1 at the compactified infinity.
+     * @param mesh Reference to the underlying serial MFEM mesh which has been promoted to high-order and projected into the curvilinear domain.
+     * @param reference_mesh reference to the underlying serial which has not been promoted to high-order or projected into the curvilinear domain. This is used to compute the compactification coordinate.
+     * @param config Config file
+     * @return Unique pointer to a scalar mesh field representing the compactification coordinate.
+     */
+    std::unique_ptr<ScalarMeshField> BuildExteriorCoordinate(
+        mfem::Mesh& mesh,
+        mfem::Mesh& reference_mesh,
+        const fourdst::config::Config<config::MeshConfig>& config
+    );
 }

src/include/stroid/topology/mapping.h

@@ -28,8 +28,23 @@ namespace stroid::topology {
     /**
      * @brief Map a point from the initial block topology to the curvilinear domain.
      * @param pos Position vector updated in-place.
-     * @param config Mesh configuration (uses radii, flattening, instability radius, and core steepness).
+     * @param config Mesh configuration (uses radii, flattening, and `core_mapping`).
+     * The `multi_block` strategy requires the matching skeleton from BuildSkeleton;
+     * changing only the mapping on a legacy core element is not supported.
      * @param attribute_id Element attribute ID (currently unused).
      */
     void TransformPoint(mfem::Vector& pos, const fourdst::config::Config<config::MeshConfig> &config, int attribute_id);
+
+    /**
+     * @brief Compute the compactification coordinate for a point in the curvilinear domain. This ranges from 0-1 with 0 at the stellar surface and 1 at the compactified infinity.
+     * @param logical_position Logical position of the point in the curvilinear domain.
+     * @param attribute Element attribute ID (used to determine the exterior coordinate).
+     * @param config Mesh configuration (uses radii and flattening).
+     * @return Compactification coordinate ranging from 0 (stellar surface) to 1 (compactified infinity).
+     */
+    double ComputeExteriorCoordinate(
+        const mfem::Vector& logical_position,
+        int attribute,
+        const fourdst::config::Config<config::MeshConfig>& config
+    );
 }

src/include/stroid/topology/optimize.h

@@ -0,0 +1,19 @@
+#pragma once
+
+#include "mfem.hpp"
+#include "fourdst/config/base.h"
+#include "stroid/config/config.h"
+
+
+namespace stroid::topology {
+
+    /**
+     * @breif Apply target matrix optimization to improve conditioning of the mesh
+     */
+    void ApplyTMOP(mfem::Mesh& mesh, const fourdst::config::Config<config::MeshConfig> &config);
+
+    /**
+     *@breif Helper to call TMOP if the correct flags are set
+     */
+    void OptimizeMesh(mfem::Mesh& mesh, const fourdst::config::Config<config::MeshConfig> &cfg);
+}

src/include/stroid/topology/topology.h

@@ -8,7 +8,9 @@
 namespace stroid::topology {
     /**
      * @brief Build the initial multi-block mesh topology for the star model.
-     * @param config Mesh configuration (uses radii and domain flags).
+     * @param config Mesh configuration (uses radii, domain flags, and `core_mapping`).
+     * The legacy `spherified` core uses one block; `multi_block` uses an
+     * inner Cartesian block and six core transition blocks.
      * @return Newly allocated mesh skeleton (not yet refined or curved).
      */
     std::unique_ptr<mfem::Mesh> BuildSkeleton(const fourdst::config::Config<config::MeshConfig> & config);

src/include/stroid/utils/mesh_stats.h

@@ -0,0 +1,174 @@
+#pragma once
+
+#include "mfem.hpp"
+#include "stroid/utils/types.h"
+#include "stroid/config/config.h"
+
+#include <cstdint>
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace stroid::stats {
+    enum class MeshStatFeatures : uint32_t {
+        NONE                 = 0u,
+        RADIUS               = 1u << 0,
+        AXES                 = 1u << 1,
+        ELLIPTICITY          = 1u << 2,
+        BOWING               = 1u << 3,
+        CONFORMITY           = 1u << 4,
+        JACOBIAN             = 1u << 5,
+        VOLUME_AREA          = 1u << 6,
+        ELEMENT_COUNT        = 1u << 7,
+        MESH_SIZE            = 1u << 8,
+        OUTER_BOUNDS         = 1u << 9,
+        CENTROID             = 1u << 10,
+        CONFIG_META          = 1u << 11,
+        BOUNDING_BOX         = 1u << 12,
+    };
+
+    constexpr MeshStatFeatures operator|(MeshStatFeatures lhs, MeshStatFeatures rhs) {
+        return static_cast<MeshStatFeatures>(static_cast<uint32_t>(lhs) | static_cast<uint32_t>(rhs));
+    }
+
+    constexpr MeshStatFeatures operator&(MeshStatFeatures lhs, MeshStatFeatures rhs) {
+        return static_cast<MeshStatFeatures>(static_cast<uint32_t>(lhs) & static_cast<uint32_t>(rhs));
+    }
+
+    constexpr bool has_feature(MeshStatFeatures feature, MeshStatFeatures set) {
+        return (static_cast<uint32_t>(set) & static_cast<uint32_t>(feature)) != 0u;
+    }
+
+    inline constexpr MeshStatFeatures MESH_STAT_DEFAULT =
+        MeshStatFeatures::RADIUS | MeshStatFeatures::AXES | MeshStatFeatures::ELLIPTICITY |
+        MeshStatFeatures::CONFORMITY | MeshStatFeatures::CONFIG_META;
+
+    inline constexpr auto MESH_STAT_ALL = static_cast<MeshStatFeatures>(0xFFFFFFFFu);
+
+    struct RadiusStats {
+        double min = 0, max = 0, mean = 0, stddev = 0;
+        long n_samples = 0;
+    };
+
+    struct AxisStats {
+        double semi_major = 0;
+        double semi_minor = 0;
+    };
+
+    struct EllipticityStats {
+        double flattening = 0;
+        double polar_equatorial = 1;
+        double radius_uniformity = 1;
+    };
+
+    struct BowingStats {
+        double max_inward = 0;
+        double max_outward = 0;
+        double rms = 0;
+        double worst_at_radius = 0;
+    };
+
+    struct ConformityStats {
+        bool conforming = true;
+        long n_nonconforming_faces = 0;
+    };
+
+    struct JacobianStats {
+        double detJ_min;
+        double detJ_max;
+        long   n_flipped;
+        double min_detJ_ratio;
+        double worst_ratio_at_radius;
+        double detJ_min_at_radius;
+        long   n_elements;
+    };
+    struct VolumeAreaStats {
+        double stellar_volume = 0, surface_area = 0;
+        double analytic_volume = 0, analytic_area = 0;
+    };
+
+    struct ElementCounts {
+        long total = 0, core = 0, envelope = 0, vacuum = 0, other = 0;
+        long n_vertices = 0;
+    };
+
+    struct MeshSizeStats {
+        double h_min = 0, h_max = 0, h_mean = 0, h_stddev = 0;
+    };
+
+    struct OuterBoundsStats {
+        double min = 0, max = 0, mean = 0;
+        long n_samples = 0;
+    };
+
+    struct CentroidStats {
+        double x = 0, y = 0, z = 0, offset = 0;
+    };
+
+    struct ConfigMeta {
+        double r_core = 0, r_star = 0, flattening = 0, r_infinity = 0;
+        int geom_order = 0;
+        size_t refinement_levels = 0;
+        bool has_external_domain = true;
+    };
+
+    struct BoundingBox {
+        double xMin = 0, xMax = 0, yMin = 0, yMax = 0, zMin = 0, zMax = 0;
+        bool valid = false;
+
+        [[nodiscard]] double dx() const {return xMax - xMin;}
+        [[nodiscard]] double dy() const {return yMax - yMin;}
+        [[nodiscard]] double dz() const {return zMax - zMin;}
+        [[nodiscard]] double diag() const {
+            const double a = dx(), b = dy(), c = dz();
+            return std::sqrt(a*a + b*b + c*c);
+        }
+    };
+
+    struct BoundingBoxStats {
+        BoundingBox core;
+        BoundingBox star;
+        BoundingBox vacuum;
+    };
+
+    struct MeshStats {
+        MeshStatFeatures computed = MeshStatFeatures::NONE;
+        std::optional<RadiusStats>                      radius;
+        std::optional<AxisStats>                        axes;
+        std::optional<EllipticityStats>                 ellipticity;
+        std::optional<BowingStats>                      bowing;
+        std::optional<ConformityStats>                  conformity;
+        std::optional<JacobianStats>                    jacobian;
+        std::optional<JacobianStats>                    jacobian_stellar;
+        std::optional<JacobianStats>                    jacobian_vacuum;
+        std::optional<VolumeAreaStats>                  volume;
+        std::optional<ElementCounts>                    element_counts;
+        std::optional<MeshSizeStats>                    mesh_size;
+        std::optional<OuterBoundsStats>                 outer_bounds;
+        std::optional<CentroidStats>                    centroid;
+        std::optional<ConfigMeta>                       config_meta;
+        std::optional<BoundingBoxStats>                 bounding_box;
+
+        std::vector<std::string>                        warnings;
+        std::vector<std::string>                        errors;
+    };
+
+    MeshStats ComputeMeshStats(const StroidMesh& sm, MeshStatFeatures features = MESH_STAT_DEFAULT, int sample_order = -1);
+
+    std::string to_string(const MeshStats& s);
+
+    inline std::ostream& operator<<(std::ostream& os, const MeshStats& s) {
+        return os << to_string(s);
+    }
+}
+
+template <>
+struct std::formatter<stroid::stats::MeshStats, char> {
+    static constexpr auto parse(const std::format_parse_context& ctx) {
+        return ctx.begin();
+    }
+
+    static auto format(const stroid::stats::MeshStats &s, std::format_context& ctx) {
+        return std::format_to(ctx.out(), "{}", stroid::stats::to_string(s));
+    }
+};

src/include/stroid/utils/mesh_utils.h

@@ -2,6 +2,9 @@
 
 #include "mfem.hpp"
 
+#include "stroid/config/config.h"
+#include "fourdst/config/config.h"
+
 namespace stroid::utils {
     /**
      * @brief Mark elements with negative Jacobian determinant.
@@ -15,4 +18,9 @@ namespace stroid::utils {
      * @param mesh Mesh to scan and update in-place.
      */
     void MarkFlippedBoundaryElements(mfem::Mesh& mesh);
+
+    void ExportJacobianRadialProfile(mfem::Mesh& mesh, const std::string& filename);
+
+    std::unique_ptr<mfem::Mesh> BuildProjected(const mfem::Mesh& reference, const fourdst::config::Config<config::MeshConfig>& cfg);
+
 }

src/include/stroid/utils/types.h

@@ -0,0 +1,83 @@
+#pragma once
+
+#include "mfem.hpp"
+
+#include "stroid/config/config.h"
+
+#include <memory>
+#include <expected>
+#include <string>
+#include <unordered_map>
+#include <variant>
+
+namespace stroid {
+    enum class MFEM_MESH_TYPE {
+        SERIAL,
+        PARALLEL
+    };
+
+    struct ScalarMeshField {
+        std::unique_ptr<mfem::FiniteElementSpace> space;
+        std::unique_ptr<mfem::GridFunction> values;
+    };
+
+    struct StroidMesh {
+        MFEM_MESH_TYPE type;
+        std::unique_ptr<mfem::Mesh> mesh;
+        std::unique_ptr<mfem::Mesh> reference_mesh;
+        std::unique_ptr<ScalarMeshField> exterior_coordinate;
+        config::MeshConfig config;
+        size_t refinement_levels;
+
+        [[nodiscard]] std::expected<mfem::Mesh*, std::string> as_mesh() const {
+            if (type == MFEM_MESH_TYPE::SERIAL) {
+                return mesh.get();
+            }
+            return std::unexpected{"Mesh is not serial. Try calling as_par_mesh()"};
+        }
+
+        [[nodiscard]] std::expected<mfem::Mesh*, std::string> ref_as_mesh() const {
+            if (type == MFEM_MESH_TYPE::SERIAL) {
+                return reference_mesh.get();
+            }
+            return std::unexpected{"Reference mesh is not serial. Try calling as_par_mesh()"};
+        }
+
+        [[nodiscard]] std::expected<std::unordered_map<std::string, std::variant<int, double, std::string, bool>>, std::string> mesh_stats(bool use_ref_mesh = false) const {
+            if (type != MFEM_MESH_TYPE::SERIAL) {
+                return std::unexpected{"Mesh is not serial. Mesh stats currently only supports serial meshes."};
+            }
+
+            mfem::Mesh* umesh;
+            if (use_ref_mesh) {
+                umesh = reference_mesh.get();
+            } else {
+                umesh = mesh.get();
+            }
+            std::unordered_map<std::string, std::variant<int, double, std::string, bool>> mesh_stats;
+            mesh_stats.emplace("num_elements", umesh->GetNE());
+            mesh_stats.emplace("num_vertices", umesh->GetNV());
+            mesh_stats.emplace("num_edges", umesh->GetNEdges());
+            mesh_stats.emplace("num_faces", umesh->GetNFaces());
+            mesh_stats.emplace("num_boundary_elements", umesh->GetNBE());
+            mesh_stats.emplace("max_bdr_attribute_id", umesh->bdr_attributes.Max());
+            mesh_stats.emplace("min_bdr_attribute_id", umesh->bdr_attributes.Min());
+            mesh_stats.emplace("max_element_attribute_id", umesh->attributes.Max());
+            mesh_stats.emplace("min_element_attribute_id", umesh->attributes.Min());
+
+            return mesh_stats;
+        }
+
+        std::unique_ptr<StroidMesh> clone() const {
+            std::unique_ptr<StroidMesh> new_mesh;
+            new_mesh->type = type;
+            new_mesh->mesh = std::make_unique<mfem::Mesh>(*mesh);
+            new_mesh->reference_mesh = std::make_unique<mfem::Mesh>(*reference_mesh);
+            new_mesh->config = config;
+            new_mesh->refinement_levels = refinement_levels;
+
+            return new_mesh;
+        }
+
+    };
+}

src/include/stroid/version.h.in

@@ -0,0 +1,46 @@
+#pragma once
+
+#include <string>
+#include <ostream>
+
+namespace stroid {
+    /**
+     * @brief Version helpers for the stroid library.
+     */
+    struct version {
+        static constexpr int major = @STROID_VERSION_MAJOR@;
+        static constexpr int minor = @STROID_VERSION_MINOR@;
+        static constexpr int patch = @STROID_VERSION_PATCH@;
+        static constexpr const char* tag = "@STROID_VERSION_TAG@";
+
+        static std::string toString() {
+            std::string versionStr = std::to_string(major) + "." +
+                                     std::to_string(minor) + "." +
+                                     std::to_string(patch);
+            if (std::string(tag) != "") {
+                versionStr += "-" + std::string(tag);
+            }
+            return versionStr;
+        }
+
+        friend std::ostream& operator<<(std::ostream& os, const version&) {
+            os << toString();
+            return os;
+        }
+    };
+}
+
+/**
+ * @namespace std
+ * @brief Standard library extensions used by stroid.
+ *
+ * Provides a `std::formatter` specialization for `stroid::version` so it can
+ * be used with `std::format` and related APIs.
+ */
+// Overload format struct
+template <>
+struct std::formatter<stroid::version> : std::formatter<std::string> {
+    auto format(const stroid::version& v, auto& ctx) {
+        return std::formatter<std::string>::format(stroid::version::toString(), ctx);
+    }
+};

src/lib/IO/mesh.cpp

@@ -1,20 +1,653 @@
 #include "mfem.hpp"
 #include "stroid/config/config.h"
 #include "stroid/IO/mesh.h"
+#include "stroid/topology/curvilinear.h"
+
+#include <algorithm>
+#include <charconv>
+#include <cmath>
+
+#include "stroid/version.h"
 
 #include <fstream>
+#include <iomanip>
 #include <iostream>
 #include <cstdint>
+#include <format>
+#include <chrono>
+#include <string>
+#include <string_view>
+#include <expected>
+#include <stdexcept>
+#include <concepts>
+#include <limits>
+#include <vector>
 
 namespace stroid::IO {
 
+    namespace {
+        std::string format_header(const StroidMesh& mesh, const std::string& comment) {
+            auto now = std::chrono::system_clock::now();
+            version v;
+
+            std::stringstream vs;
+            vs << v;
+
+            std::string header = std::format(R"(# STROID MESH
+# NOTE: STROID MESH IS A THIN WRAPPER AROUND MFEM's NATIVE MESH FORMAT
+# STRUCTURE:
+#           - Type                : Serial or Parallel (S for Serial, P for Parallel)
+#           - mesh                : the primary computational domain which can be of n order and be h-refined
+#           - reference mesh      : a reference, linear order mesh, used to ensure that the primary mesh remains well formed
+#           - exterior coordinate : a scalar material coordinate which is zero at the stellar surface and one at infinity
+#           - config              : The configuration options initially used to generate the mesh
+#           - refinement-levels   : the total number of refinement levels the primary mesh has been subjected too
+# NOTE: EACH BLOCK OF DATA IS STORED BETWEEN "BEGIN BLOCK <NAME>\n ... \nEND BLOCK <NAME>
+#       PARSING THE UNDERLYING MFEM NATIVE MESH FORMAT CAN BE DONE WITH MFEM'S STREAM READER
+#       IF YOU EXTRACT THE RAW CONTENTS BETWEEN THOSE LINES
+BEGIN BLOCK HEADER
+    MESH_TYPE:{}
+    REFINEMENT_LEVELS:{}
+    DATE_CREATED:{:%Y-%m-%d}
+    COMMENT:{}
+    STROID_VERSION:{}
+END BLOCK HEADER)",
+            mesh.type == MFEM_MESH_TYPE::PARALLEL ? "P" : "S",
+            mesh.refinement_levels,
+            now,
+            comment,
+            vs.str(),
+            mesh.refinement_levels
+            );
+        return header;
+        }
+
+        std::string format_primary_mesh(const StroidMesh& mesh) {
+            std::stringstream ss;
+            ss.precision(std::numeric_limits<double>::max_digits10);
+            mesh.mesh->Print(ss);
+
+            std::string pmesh = std::format("BEGIN BLOCK PMESH\n{}END BLOCK PMESH", ss.str());
+
+            return pmesh;
+        }
+
+        template <typename T>
+        std::string format_opt(const std::optional<T> opt, T default_val) {
+            if (opt.has_value()) {
+                return std::format("{}", opt.value());
+            }
+            return std::format("{}", default_val);
+
+        }
+
+        std::string format_reference_mesh(const StroidMesh& mesh) {
+            std::stringstream ss;
+            ss.precision(std::numeric_limits<double>::max_digits10);
+            mesh.reference_mesh->Print(ss);
+
+            std::string rmesh = std::format("BEGIN BLOCK RMESH\n{}END BLOCK RMESH", ss.str());
+            return rmesh;
+        }
+
+        std::string format_config(const StroidMesh& mesh) {
+            config::MeshConfig d;
+
+            config::OptimizationMethods d_opt = d.optimization_methods.value_or(config::OptimizationMethods{false, true});
+            config::OptimizationMethods m_opt = mesh.config.optimization_methods.value_or(d_opt);
+
+            std::string config_str = std::format(R"(BEGIN BLOCK CONFIG
+# refiniment_levels: Initial number of levels of refinmenet, note the value in the header may be more up to date
+# std::optional<int>
+# default: 4
+refinement_levels:{}
+
+# order: Polynomial / geometric order to use when constructing the mesh
+# std::optional<int>
+# default: 3
+order:{}
+
+# include_external_domain: Whether or not to include the external domain in the mesh generally used for applying boundary conditions at infinity
+# std::optional<bool>
+# default: true
+include_external_domain:{}
+
+# r_core: the radius of the stellar core region (in reference space)
+# std::optional<double>
+# default: 0.25
+r_core:{}
+
+# r_star: the radius of the stellar surface (in reference space)
+# std::optional<double>
+# default: 1.0
+r_star:{}
+
+# flattening: the flattening of the star (in reference space) where 0 is spherical and >0 is oblate. Note that this parameter is not equivalent to solving for the structure of a rotating model
+# std::optional<float>
+# default: 0.0
+flattening:{}
+
+# r_infinity: the radius of the outer boundary of the mesh (in reference space)
+# std::optional<double>
+# default: 6.0
+r_infinity:{}
+
+# r_instability: the radius inside which computations of geometry are skipped to avoid a core singularity
+# std::optional<double>
+# default: 1e-14
+r_instability:{}
+
+# core_steepness: Controls the rate of transition of the core-to-envelope transition
+# std::optional<double>
+# default: 1.0
+core_steepness:{}
+
+# continuity_order: order of continuity to force from teh core-envelope transition (0 = discontinuous, 1=C1 continuity, etc...)
+# std::optional<double>
+# default: 2
+continuity_order:{}
+
+# surface_bdr_id: the boundary id to tag the stellar surface boundary elements as
+# std::optional<size_t>
+# default: 1
+surface_bdr_id:{}
+
+# inf_bdr_id: the boundary id to tag the outer boundary elements as
+# std::optional<size_t>
+# default: 2
+inf_bdr_id:{}
+
+# core_id: the material attribute to tag elements in the core region as
+# std::optional<size_t>
+# default 1
+core_id:{}
+
+# envelope_id: the material attribute to tag elements in the envelope as
+# std::optional<size_t>
+# default 2
+envelope_id:{}
+
+# vacuum_id: the material attribute to tag elements in the vacuum region as
+# std::optional<size_t>
+# default 3
+vacuum_id:{}
+
+# optimization_method: struct for storing which optimization methods are being used
+# includes tmop and smoothstep booleans
+optimization_methods-tmop:{}
+optimization_methods-smoothstep:{}
+
+# core_mapping: Core mapping strategy, either spherified or multi_block
+# std::optional<std::string>
+# default: spherified
+core_mapping:{}
+END BLOCK CONFIG)",
+            format_opt(mesh.config.refinement_levels, d.refinement_levels.value()),
+            format_opt(mesh.config.order, d.order.value()),
+            format_opt(mesh.config.include_external_domain, d.include_external_domain.value()),
+            format_opt(mesh.config.r_core, d.r_core.value()),
+            format_opt(mesh.config.r_star, d.r_star.value()),
+            format_opt(mesh.config.flattening, d.flattening.value()),
+            format_opt(mesh.config.r_infinity, d.r_infinity.value()),
+            format_opt(mesh.config.r_instability, d.r_instability.value()),
+            format_opt(mesh.config.core_steepness, d.core_steepness.value()),
+            format_opt(mesh.config.continuity_order, d.continuity_order.value()),
+            format_opt(mesh.config.surface_bdr_id, d.surface_bdr_id.value()),
+            format_opt(mesh.config.inf_bdr_id, d.inf_bdr_id.value()),
+            format_opt(mesh.config.core_id, d.core_id.value()),
+            format_opt(mesh.config.envelope_id, d.envelope_id.value()),
+            format_opt(mesh.config.vacuum_id, d.vacuum_id.value()),
+            m_opt.tmop.value_or(false),
+            m_opt.smoothstep.value_or(true),
+            format_opt(mesh.config.core_mapping, d.core_mapping.value()));
+
+            return config_str;
+        }
+
+        std::string format_exterior_coordinate(const StroidMesh& mesh) {
+            const bool include_external_domain = mesh.config.include_external_domain.value_or(true);
+
+            if (!include_external_domain) {
+                if (mesh.exterior_coordinate) throw std::runtime_error("A mesh without an external domain cannot contain an exterior-coordinate field.");
+                return "BEGIN BLOCK EXTERIOR_COORDINATE\nPRESENT:false\nEND BLOCK EXTERIOR_COORDINATE";
+            }
+
+            if (!mesh.exterior_coordinate || !mesh.exterior_coordinate->space || !mesh.exterior_coordinate->values) {
+                throw std::runtime_error("A mesh with an external domain must contain a complete exterior-coordinate field before it can be saved.");
+            }
+            if (mesh.exterior_coordinate->space->GetMesh() != mesh.mesh.get()) {
+                throw std::runtime_error("The exterior-coordinate finite-element space is attached to the wrong mesh.");
+            }
+            if (mesh.exterior_coordinate->values->FESpace() != mesh.exterior_coordinate->space.get()) {
+                throw std::runtime_error("The exterior-coordinate grid function is attached to the wrong finite-element space.");
+            }
+
+            const int scalar_dofs = mesh.exterior_coordinate->space->GetNDofs();
+            if (mesh.exterior_coordinate->values->Size() != scalar_dofs) {
+                throw std::runtime_error("The exterior-coordinate grid function has an invalid size.");
+            }
+
+            std::stringstream ss;
+            ss << std::setprecision(std::numeric_limits<double>::max_digits10);
+            ss << "BEGIN BLOCK EXTERIOR_COORDINATE\n";
+            ss << "PRESENT:true\n";
+            ss << "NDOFS:" << scalar_dofs << '\n';
+            ss << "VALUES:\n";
+
+            for (int dof = 0; dof < scalar_dofs; ++dof) {
+                const double coordinate = (*mesh.exterior_coordinate->values)(dof);
+                if (!std::isfinite(coordinate) || coordinate < 0.0 || coordinate > 1.0) {
+                    throw std::runtime_error(std::format("Exterior-coordinate DOF {} has invalid value {}.", dof, coordinate));
+                }
+                ss << coordinate << '\n';
+            }
+
+            ss << "END BLOCK EXTERIOR_COORDINATE";
+            return ss.str();
+        }
+    }
+
+    namespace {
+        constexpr std::string_view BEGIN_PREFIX = "BEGIN BLOCK ";
+        constexpr std::string_view END_PREFIX   = "END BLOCK ";
+
+        std::string_view trim(std::string_view s) {
+            const auto b = s.find_first_not_of(" \t\r\n");
+            if (b == std::string_view::npos) return {};
+            const auto e = s.find_last_not_of(" \t\r\n");
+            return s.substr(b, e - b + 1);
+        }
+
+
+        std::expected<bool, std::string> parse_bool(std::string_view v) {
+            std::string s(trim(v));
+            std::ranges::transform(s, s.begin(),
+                                   [](const unsigned char c) { return static_cast<char>(std::tolower(c)); });
+            if (s == "true"  || s == "1") return true;
+            if (s == "false" || s == "0") return false;
+            return std::unexpected(std::format("invalid bool value '{}'", v));
+        }
+
+        template <std::integral T>
+        std::expected<T, std::string> parse_int(std::string_view v) {
+            const std::string_view s = trim(v);
+
+            std::string temp(s);
+
+            try {
+                size_t pos = 0;
+
+                if constexpr (std::is_signed_v<T>) {
+                    long long val = std::stoll(temp, &pos);
+
+                    if (pos != temp.size() || val < std::numeric_limits<T>::min() || val > std::numeric_limits<T>::max()) {
+                        return std::unexpected(std::format("invalid integer value '{}'", v));
+                    }
+                    return static_cast<T>(val);
+
+                } else {
+                    unsigned long long val = std::stoull(temp, &pos);
+
+                    if (pos != temp.size() || val > std::numeric_limits<T>::max()) {
+                        return std::unexpected(std::format("invalid integer value '{}'", v));
+                    }
+                    return static_cast<T>(val);
+                }
+            } catch (const std::exception&) {
+                return std::unexpected(std::format("invalid integer value '{}'", v));
+            }
+        }
+
+        std::expected<double, std::string> parse_double(std::string_view v) {
+            const std::string_view s = trim(v);
+            std::string temp(s);
+
+            try {
+                size_t pos = 0;
+                double out = std::stod(temp, &pos);
+
+                if (pos != temp.size()) {
+                    return std::unexpected(std::format("invalid floating-point value '{}'", v));
+                }
+                return out;
+
+            } catch (const std::exception&) {
+                return std::unexpected(std::format("invalid floating-point value '{}'", v));
+            }
+        }
+
+        std::expected<std::map<std::string, std::string>, std::string> extract_blocks(std::istream& is) {
+            std::map<std::string, std::string> blocks;
+            std::string line;
+            std::string current;
+            std::string buffer;
+            bool in_block = false;
+
+            while (std::getline(is, line)) {
+                const std::string_view t = trim(line);
+
+                if (!in_block) {
+                    if (t.starts_with(BEGIN_PREFIX)) {
+                        current = std::string(trim(t.substr(BEGIN_PREFIX.size())));
+                        if (current.empty())
+                            return std::unexpected("found 'BEGIN BLOCK' with no block name");
+                        if (blocks.contains(current))
+                            return std::unexpected(std::format("duplicate block '{}'", current));
+                        buffer.clear();
+                        in_block = true;
+                    }
+                } else {
+                    if (t.starts_with(END_PREFIX)) {
+                        if (const std::string end_name(trim(t.substr(END_PREFIX.size()))); end_name != current)
+                            return std::unexpected(std::format(
+                                "mismatched block markers: opened '{}' but closed '{}'",
+                                current, end_name));
+                        blocks.emplace(std::move(current), std::move(buffer));
+                        current.clear();
+                        buffer.clear();
+                        in_block = false;
+                    } else {
+                        std::string_view raw = line;
+                        if (!raw.empty() && raw.back() == '\r') raw.remove_suffix(1);
+                        buffer.append(raw);
+                        buffer.push_back('\n');
+                    }
+                }
+            }
+
+            if (in_block)
+                return std::unexpected(std::format("unterminated block '{}' (missing END BLOCK)", current));
+            return blocks;
+        }
+
+        std::expected<void, std::string> parse_header(const std::string& content, StroidMesh& out) {
+            std::istringstream iss(content);
+            std::string line;
+            std::optional<MFEM_MESH_TYPE> type;
+            std::optional<size_t> ref_levels;
+
+            while (std::getline(iss, line)) {
+                const std::string_view t = trim(line);
+                if (t.empty() || t.starts_with('#')) continue;
+
+                const auto colon = t.find(':');
+                if (colon == std::string_view::npos) continue;
+
+                const std::string_view key = trim(t.substr(0, colon));
+                const std::string_view val = trim(t.substr(colon + 1));
+
+                if (key == "MESH_TYPE") {
+                    if      (val == "P") type = MFEM_MESH_TYPE::PARALLEL;
+                    else if (val == "S") type = MFEM_MESH_TYPE::SERIAL;
+                    else return std::unexpected(std::format("unknown MESH_TYPE '{}'", val));
+                } else if (key == "REFINEMENT_LEVELS") {
+                    auto r = parse_int<size_t>(val);
+                    if (!r) return std::unexpected("REFINEMENT_LEVELS: " + r.error());
+                    ref_levels = *r;
+                }
+            }
+
+            if (!type) return std::unexpected("HEADER block missing MESH_TYPE");
+            out.type = *type;
+            out.refinement_levels = ref_levels.value_or(0);
+            return {};
+        }
+
+        std::expected<config::MeshConfig, std::string> parse_config(const std::string& content) {
+            config::MeshConfig cfg;
+            config::OptimizationMethods opt =
+                cfg.optimization_methods.value_or(config::OptimizationMethods{});
+
+            using Handler = std::function<std::expected<void, std::string>(std::string_view)>;
+
+            auto as_int    = [](std::optional<int>* f)    { return [f](const std::string_view v) -> std::expected<void, std::string> { auto r = parse_int<int>(v);    if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; };
+            auto as_size   = [](std::optional<size_t>* f) { return [f](const std::string_view v) -> std::expected<void, std::string> { auto r = parse_int<size_t>(v); if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; };
+            auto as_double = [](std::optional<double>* f) { return [f](const std::string_view v) -> std::expected<void, std::string> { auto r = parse_double(v);     if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; };
+            auto as_bool   = [](std::optional<bool>* f)   { return [f](const std::string_view v) -> std::expected<void, std::string> { auto r = parse_bool(v);       if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; };
+            auto as_string = [](std::optional<std::string>* f) { return [f](const std::string_view v) -> std::expected<void, std::string> { *f = std::string(v); return {}; }; };
+
+            const std::unordered_map<std::string_view, Handler> handlers = {
+                {"refinement_levels",               as_int(&cfg.refinement_levels)},
+                {"order",                           as_int(&cfg.order)},
+                {"include_external_domain",         as_bool(&cfg.include_external_domain)},
+                {"r_core",                          as_double(&cfg.r_core)},
+                {"r_star",                          as_double(&cfg.r_star)},
+                {"flattening",                      as_double(&cfg.flattening)},
+                {"r_infinity",                      as_double(&cfg.r_infinity)},
+                {"r_instability",                   as_double(&cfg.r_instability)},
+                {"core_steepness",                  as_double(&cfg.core_steepness)},
+                {"continuity_order",                as_size(&cfg.continuity_order)},
+                {"surface_bdr_id",                  as_size(&cfg.surface_bdr_id)},
+                {"inf_bdr_id",                      as_size(&cfg.inf_bdr_id)},
+                {"core_id",                         as_size(&cfg.core_id)},
+                {"envelope_id",                     as_size(&cfg.envelope_id)},
+                {"vacuum_id",                       as_size(&cfg.vacuum_id)},
+                {"optimization_methods-tmop",       as_bool(&opt.tmop)},
+                {"optimization_methods-smoothstep", as_bool(&opt.smoothstep)},
+                {"core_mapping",                    as_string(&cfg.core_mapping)},
+            };
+
+            std::istringstream iss(content);
+            std::string line;
+            while (std::getline(iss, line)) {
+                const std::string_view t = trim(line);
+                if (t.empty() || t.starts_with('#')) continue;
+
+                const auto colon = t.find(':');
+                if (colon == std::string_view::npos) continue;
+
+                const std::string_view key = trim(t.substr(0, colon));
+                const std::string_view val = trim(t.substr(colon + 1));
+
+                const auto it = handlers.find(key);
+                if (it == handlers.end()) continue;
+                if (auto r = it->second(val); !r)
+                    return std::unexpected(std::format("{}: {}", key, r.error()));
+            }
+
+            cfg.optimization_methods = opt;
+            return cfg;
+        }
+
+        std::expected<std::unique_ptr<mfem::Mesh>, std::string> load_serial_mesh(const std::string& raw) {
+            if (trim(raw).empty()) return std::unexpected("empty mesh block");
+            std::istringstream iss(raw);
+            try {
+                return std::make_unique<mfem::Mesh>(iss);
+            } catch (const std::exception& e) {
+                return std::unexpected(std::string("MFEM failed to parse mesh: ") + e.what());
+            }
+        }
+
+        struct ParsedMeta {
+            StroidMesh  mesh;
+            std::string pmesh_raw;
+            std::string rmesh_raw;
+            std::optional<std::string> exterior_coordinate_raw;
+        };
+
+        struct ParsedExteriorCoordinate {
+            bool present{false};
+            int scalar_dofs{0};
+            std::vector<double> values;
+        };
+
+        std::expected<ParsedExteriorCoordinate, std::string> parse_exterior_coordinate(const std::string& content) {
+            ParsedExteriorCoordinate parsed;
+            std::optional<bool> present;
+            std::optional<int> scalar_dofs;
+            bool reading_values = false;
+            std::istringstream iss(content);
+            std::string line;
+
+            while (std::getline(iss, line)) {
+                const std::string_view value = trim(line);
+                if (value.empty() || value.starts_with('#')) continue;
+
+                if (reading_values) {
+                    auto coordinate = parse_double(value);
+                    if (!coordinate) return std::unexpected("EXTERIOR_COORDINATE value -> " + coordinate.error());
+                    parsed.values.push_back(*coordinate);
+                    continue;
+                }
+
+                const auto colon = value.find(':');
+                if (colon == std::string_view::npos) return std::unexpected(std::format("invalid EXTERIOR_COORDINATE line '{}'.", value));
+
+                const std::string_view key = trim(value.substr(0, colon));
+                const std::string_view field_value = trim(value.substr(colon + 1));
+
+                if (key == "PRESENT") {
+                    auto result = parse_bool(field_value);
+                    if (!result) return std::unexpected("EXTERIOR_COORDINATE PRESENT -> " + result.error());
+                    present = *result;
+                } else if (key == "NDOFS") {
+                    auto result = parse_int<int>(field_value);
+                    if (!result) return std::unexpected("EXTERIOR_COORDINATE NDOFS -> " + result.error());
+                    scalar_dofs = *result;
+                } else if (key == "VALUES") {
+                    if (!field_value.empty()) return std::unexpected("EXTERIOR_COORDINATE VALUES must not contain an inline value.");
+                    reading_values = true;
+                } else {
+                    return std::unexpected(std::format("unknown EXTERIOR_COORDINATE key '{}'.", key));
+                }
+            }
+
+            if (!present.has_value()) return std::unexpected("EXTERIOR_COORDINATE block is missing PRESENT.");
+            parsed.present = *present;
+
+            if (!parsed.present) {
+                if (scalar_dofs.has_value() || !parsed.values.empty()) return std::unexpected("An absent exterior coordinate cannot contain NDOFS or VALUES.");
+                return parsed;
+            }
+
+            if (!scalar_dofs.has_value() || *scalar_dofs < 0) return std::unexpected("EXTERIOR_COORDINATE block has an invalid or missing NDOFS.");
+            if (static_cast<int>(parsed.values.size()) != *scalar_dofs) {
+                return std::unexpected(std::format("EXTERIOR_COORDINATE expected {} values but found {}.", *scalar_dofs, parsed.values.size()));
+            }
+
+            parsed.scalar_dofs = *scalar_dofs;
+            return parsed;
+        }
+
+        std::expected<void, std::string> restore_exterior_coordinate(StroidMesh& mesh, const std::optional<std::string>& raw) {
+            fourdst::config::Config<config::MeshConfig> config;
+            config.mutate([&mesh](config::MeshConfig& value) { value = mesh.config; });
+
+            try {
+                mesh.exterior_coordinate = topology::BuildExteriorCoordinate(*mesh.mesh, *mesh.reference_mesh, config);
+            } catch (const std::exception& exception) {
+                return std::unexpected(std::string("failed to reconstruct exterior coordinate: ") + exception.what());
+            }
+
+            if (!raw.has_value()) return {};
+
+            auto parsed = parse_exterior_coordinate(*raw);
+            if (!parsed) return std::unexpected(parsed.error());
+
+            const bool include_external_domain = mesh.config.include_external_domain.value_or(true);
+            if (!parsed->present) {
+                if (include_external_domain) return std::unexpected("EXTERIOR_COORDINATE is absent even though the mesh includes an external domain.");
+                if (mesh.exterior_coordinate) return std::unexpected("An exterior-coordinate field was reconstructed for a mesh without an external domain.");
+                return {};
+            }
+
+            if (!include_external_domain) return std::unexpected("EXTERIOR_COORDINATE is present for a mesh without an external domain.");
+            if (!mesh.exterior_coordinate || !mesh.exterior_coordinate->space || !mesh.exterior_coordinate->values) {
+                return std::unexpected("Unable to allocate the exterior-coordinate field while loading the mesh.");
+            }
+            if (parsed->scalar_dofs != mesh.exterior_coordinate->space->GetNDofs()) {
+                return std::unexpected(std::format("EXTERIOR_COORDINATE contains {} DOFs but the reconstructed space has {}.", parsed->scalar_dofs, mesh.exterior_coordinate->space->GetNDofs()));
+            }
+
+            constexpr double consistency_tolerance = 1.0e-12;
+
+            for (int dof = 0; dof < parsed->scalar_dofs; ++dof) {
+                const double stored_coordinate = parsed->values[static_cast<size_t>(dof)];
+                const double reconstructed_coordinate = (*mesh.exterior_coordinate->values)(dof);
+
+                if (!std::isfinite(stored_coordinate) || stored_coordinate < 0.0 || stored_coordinate > 1.0) {
+                    return std::unexpected(std::format("EXTERIOR_COORDINATE DOF {} has invalid stored value {}.", dof, stored_coordinate));
+                }
+                if (std::abs(stored_coordinate - reconstructed_coordinate) > consistency_tolerance) {
+                    return std::unexpected(std::format("EXTERIOR_COORDINATE DOF {} is inconsistent with the reference mesh: stored value {}, reconstructed value {}.", dof, stored_coordinate, reconstructed_coordinate));
+                }
+
+                (*mesh.exterior_coordinate->values)(dof) = stored_coordinate;
+            }
+
+            return {};
+        }
+
+        std::expected<ParsedMeta, std::string> parse_metadata(std::istream& is) {
+            auto blocks = extract_blocks(is);
+            if (!blocks) return std::unexpected(blocks.error());
+
+            auto need = [&](std::string_view name) -> std::expected<std::string, std::string> {
+                const auto it = blocks->find(std::string(name));
+                if (it == blocks->end())
+                    return std::unexpected(std::format("missing required block '{}'", name));
+                return it->second;
+            };
+
+            ParsedMeta pm{};
+
+            const auto header = need("HEADER");
+            if (!header) return std::unexpected(header.error());
+            if (auto r = parse_header(*header, pm.mesh); !r) return std::unexpected(r.error());
+
+            const auto config = need("CONFIG");
+            if (!config) return std::unexpected(config.error());
+            auto cfg = parse_config(*config);
+            if (!cfg) return std::unexpected("CONFIG block -> " + cfg.error());
+            pm.mesh.config = std::move(*cfg);
+
+            const auto pmesh = need("PMESH");
+            if (!pmesh) return std::unexpected(pmesh.error());
+            pm.pmesh_raw = *pmesh;
+
+            const auto rmesh = need("RMESH");
+            if (!rmesh) return std::unexpected(rmesh.error());
+            pm.rmesh_raw = *rmesh;
+
+            if (const auto exterior_coordinate = blocks->find("EXTERIOR_COORDINATE"); exterior_coordinate != blocks->end()) {
+                pm.exterior_coordinate_raw = exterior_coordinate->second;
+            }
+
+            return pm;
+        }
+
+    }
+
+    void SaveStroidMesh(const StroidMesh &mesh, const std::string &filename, const std::string &comment) {
+        std::ofstream ofs(filename);
+
+        // First Write a header with some information
+        std::string header = format_header(mesh, comment);
+
+        std::string pmesh = format_primary_mesh(mesh);
+        std::string rmesh = format_reference_mesh(mesh);
+
+        std::string config = format_config(mesh);
+        std::string exterior_coordinate = format_exterior_coordinate(mesh);
+
+        ofs << header << "\n";
+        ofs << pmesh << "\n";
+        ofs << rmesh << "\n";
+        ofs << config << "\n";
+        ofs << exterior_coordinate << "\n";
+    }
 
     void SaveMesh(const mfem::Mesh& mesh, const std::string& filename) {
         std::ofstream ofs(filename);
-        ofs.precision(8);
+        ofs.precision(std::numeric_limits<double>::max_digits10);
         mesh.Print(ofs);
     }
 
+    void SaveMesh(const stroid::StroidMesh &mesh, const std::string &filename) {
+        SaveMesh(*mesh.mesh, filename);
+    }
+
     void SaveVTU(mfem::Mesh &mesh, const std::string &exportName) {
         mfem::ParaViewDataCollection pd(exportName, &mesh);
         pd.SetDataFormat(mfem::VTKFormat::BINARY);
@@ -22,6 +655,10 @@ namespace stroid::IO {
         pd.Save();
     }
 
+    void SaveVTU(const stroid::StroidMesh &mesh, const std::string &exportName) {
+        SaveVTU(*mesh.mesh, exportName);
+    }
+
     void ViewMesh(mfem::Mesh &mesh, const std::string& title, const VISUALIZATION_MODE mode, const std::string &vishost, int visport) {
         mfem::socketstream sol_sock(vishost.c_str(), visport);
         if (!sol_sock.is_open()) {
@@ -61,7 +698,12 @@ namespace stroid::IO {
         sol_sock << "keys iMj\n";
         sol_sock << std::flush;
     }
-    void VisualizeFaceValence(mfem::Mesh& mesh) {
+
+    void ViewMesh(const stroid::StroidMesh &mesh, const std::string &title, VISUALIZATION_MODE mode, const std::string &vishost, int visport) {
+        ViewMesh(*mesh.mesh, title, mode, vishost, visport);
+    }
+
+    void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport) {
         mfem::L2_FECollection fec(0, 3);
         mfem::FiniteElementSpace fes(&mesh, &fec);
         mfem::GridFunction valence_gf(&fes);
@@ -78,13 +720,83 @@ namespace stroid::IO {
         }
 
         // View in GLVis
-        char vishost[] = "localhost";
-        int  visport   = 19916;
-        mfem::socketstream sol_sock(vishost, visport);
+        mfem::socketstream sol_sock(vishost.c_str(), visport);
         if (sol_sock.is_open()) {
             sol_sock << "solution\n" << mesh << valence_gf;
             sol_sock << "window_title 'Boundary Valence: 1=Surface, 2=Internal'\n";
             sol_sock << "keys am\n" << std::flush;
         }
     }
+
+    void VisualizeFaceValence(const stroid::StroidMesh &mesh, const std::string &vishost, int visport) {
+        VisualizeFaceValence(*mesh.mesh, vishost, visport);
+    }
+
+    std::expected<StroidMesh, std::string> ParseStroidMesh(std::istream& is) {
+        auto pm = parse_metadata(is);
+        if (!pm) return std::unexpected(pm.error());
+
+        if (pm->mesh.type != MFEM_MESH_TYPE::SERIAL) {
+            return std::unexpected(
+                "parsed a PARALLEL StroidMesh, but ParseStroidMesh(std::istream&) can only "
+                "reconstruct serial meshes; use the MPI-aware overload "
+                "ParseStroidMesh(std::istream&, MPI_Comm) (requires MFEM_USE_MPI)");
+        }
+
+        auto m = load_serial_mesh(pm->pmesh_raw);
+        if (!m) return std::unexpected("PMESH -> " + m.error());
+        auto rm = load_serial_mesh(pm->rmesh_raw);
+        if (!rm) return std::unexpected("RMESH -> " + rm.error());
+
+        pm->mesh.mesh           = std::move(*m);
+        pm->mesh.reference_mesh = std::move(*rm);
+        if (auto result = restore_exterior_coordinate(pm->mesh, pm->exterior_coordinate_raw); !result) return std::unexpected(result.error());
+        return std::move(pm->mesh);
+    }
+
+    std::expected<StroidMesh, std::string> LoadStroidMesh(const std::string& filename) {
+        std::ifstream ifs(filename);
+        if (!ifs.is_open())
+            return std::unexpected(std::format("could not open file '{}'", filename));
+        return ParseStroidMesh(ifs);
+    }
+
+#ifdef MFEM_USE_MPI
+    std::expected<StroidMesh, std::string> ParseStroidMesh(std::istream& is, MPI_Comm comm) {
+        auto pm = parse_metadata(is);
+        if (!pm) return std::unexpected(pm.error());
+
+        auto build = [&](const std::string& raw)
+            -> std::expected<std::unique_ptr<mfem::Mesh>, std::string> {
+            if (trim(raw).empty()) return std::unexpected("empty mesh block");
+            std::istringstream iss(raw);
+            try {
+                if (pm->mesh.type == MFEM_MESH_TYPE::PARALLEL)
+                    return std::unique_ptr<mfem::Mesh>(new mfem::ParMesh(comm, iss));
+                return std::make_unique<mfem::Mesh>(iss);
+            } catch (const std::exception& e) {
+                return std::unexpected(std::string("MFEM failed to parse mesh: ") + e.what());
+            }
+        };
+
+        auto m = build(pm->pmesh_raw);
+        if (!m) return std::unexpected("PMESH -> " + m.error());
+        auto rm = build(pm->rmesh_raw);
+        if (!rm) return std::unexpected("RMESH -> " + rm.error());
+
+        pm->mesh.mesh           = std::move(*m);
+        pm->mesh.reference_mesh = std::move(*rm);
+        if (auto result = restore_exterior_coordinate(pm->mesh, pm->exterior_coordinate_raw); !result) return std::unexpected(result.error());
+        return std::move(pm->mesh);
+    }
+
+    std::expected<StroidMesh, std::string> LoadStroidMesh(const std::string& filename, MPI_Comm comm) {
+        std::ifstream ifs(filename);
+        if (!ifs.is_open())
+            return std::unexpected(std::format("could not open file '{}'", filename));
+        return ParseStroidMesh(ifs, comm);
+    }
+#endif // MFEM_USE_MPI
+
+
 }

src/lib/refinement/uniform.cpp

@@ -0,0 +1,42 @@
+#include "mfem.hpp"
+
+#include "stroid/refinement/uniform.h"
+#include "stroid/utils/types.h"
+#include "stroid/utils/mesh_utils.h"
+#include "stroid/exceptions/exceptions.h"
+#include "stroid/topology/curvilinear.h"
+#include "stroid/topology/topology.h"
+#include "stroid/topology/optimize.h"
+
+namespace stroid::refinement {
+    void UniformRefinement(StroidMesh &mesh, const size_t levels) {
+        if (!mesh.reference_mesh) {
+            throw exceptions::StroidMissingReferenceMesh("UniformRefinement requires a reference mesh to be present in the StroidMesh object. This should be present by construction and the fact that is is missing represents a bug. Please report this to the stroid developers on GitHub or by email at emily.boudreaux@dartmouth.edu");
+        }
+
+        if (levels == 0) {
+            return;
+        }
+
+        if (!mesh.mesh) {
+            throw exceptions::StroidMissingReferenceMesh("UniformRefinement requires a primary mesh to be present in the StroidMesh object. This should be present by construction and the fact that it is missing represents a bug. Please report this to the stroid developers on GitHub or by email at emily.boudreaux@dartmouth.edu");
+        }
+        mesh.exterior_coordinate.reset();
+        for (size_t i = 0; i < levels; i++) {
+            mesh.reference_mesh->UniformRefinement();
+        }
+
+        mesh.refinement_levels += levels;
+
+        fourdst::config::Config<config::MeshConfig> cfg;
+        auto Mutator = [&mesh](config::MeshConfig& orig) {
+            orig = mesh.config;
+        };
+
+        cfg.mutate(Mutator);
+
+        mesh.mesh = utils::BuildProjected(*mesh.reference_mesh, cfg);
+        topology::OptimizeMesh(*mesh.mesh, cfg);
+        mesh.exterior_coordinate = topology::BuildExteriorCoordinate(*mesh.mesh, *mesh.reference_mesh, cfg);
+    }
+}

src/lib/topology/curvilinear.cpp

@@ -3,9 +3,44 @@
 
 #include <iostream>
 
+namespace {
+    double compute_exterior_coordinate(
+        const mfem::Vector& logical_position,
+        const int attribute,
+        const fourdst::config::Config<stroid::config::MeshConfig>& config
+    ) {
+        if (!config->include_external_domain.value_or(true) || attribute != static_cast<int>(config->vacuum_id.value_or(3))) return 0.0;
+
+        const double r_star = config->r_star.value_or(1.0);
+        const double r_infinity = config->r_infinity.value_or(6.0);
+        const double radial_extent = r_infinity - r_star;
+
+        if (!std::isfinite(r_star) || !std::isfinite(r_infinity) || r_star <= 0.0 || radial_extent <= 0.0) {
+            throw std::invalid_argument("Exterior-coordinate construction requires finite radii with 0 < r_star < r_infinity.");
+        }
+
+        double logical_radius = 0.0;
+        for (int d = 0; d < logical_position.Size(); ++d) {
+            if (!std::isfinite(logical_position(d))) throw std::runtime_error("Reference mesh produced a non-finite logical position.");
+            logical_radius = std::max(logical_radius, std::abs(logical_position(d)));
+        }
+
+        double coordinate = (logical_radius - r_star) / radial_extent;
+        const double tolerance = 1024.0 * std::numeric_limits<double>::epsilon() * std::max({1.0, std::abs(r_star), std::abs(r_infinity)}) / radial_extent;
+
+        if (coordinate < -tolerance || coordinate > 1.0 + tolerance) {
+            throw std::runtime_error(std::format("Logical exterior coordinate {} lies outside [0, 1].", coordinate));
+        }
+
+        if (std::abs(coordinate) <= tolerance) coordinate = 0.0;
+        if (std::abs(coordinate - 1.0) <= tolerance) coordinate = 1.0;
+        return coordinate;
+    }
+}
+
 namespace stroid::topology {
     void PromoteToHighOrder(mfem::Mesh &mesh, const fourdst::config::Config<config::MeshConfig> &config) {
-        const auto* fec = new mfem::H1_FECollection(config->order, mesh.Dimension());
+        const auto* fec = new mfem::H1_FECollection(config->order.value(), mesh.Dimension());
         auto* fes = new mfem::FiniteElementSpace(&mesh, fec, mesh.SpaceDimension());
         mesh.SetNodalFESpace(fes);
     }
@@ -53,16 +88,90 @@ namespace stroid::topology {
             }
         }
 
-        // for (int i = 0; i < nDofs; ++i) {
-        //     for (int d = 0; d < vDim; ++d) {
-        //         pos(d) = nodes(fes->DofToVDof(i, d));
-        //     }
-        //
-        //     TransformPoint(pos, config, 0);
-        //
-        //     for (int d = 0; d < vDim; ++d) {
-        //         nodes(fes->DofToVDof(i, d)) = pos(d);
-        //     }
-        // }
+    }
+
+    std::unique_ptr<ScalarMeshField> BuildExteriorCoordinate(
+        mfem::Mesh& mesh,
+        mfem::Mesh& reference_mesh,
+        const fourdst::config::Config<config::MeshConfig>& config
+    ) {
+        if (!config->include_external_domain.value_or(true)) return nullptr;
+        if (mesh.Dimension() != reference_mesh.Dimension() || mesh.SpaceDimension() != reference_mesh.SpaceDimension()) {
+            throw std::invalid_argument("Primary and reference meshes must have matching dimensions when constructing the exterior coordinate.");
+        }
+        if (mesh.GetNE() != reference_mesh.GetNE()) {
+            throw std::invalid_argument("Primary and reference meshes must have the same number of elements when constructing the exterior coordinate.");
+        }
+        if (mesh.GetNodalFESpace() == nullptr) {
+            throw std::invalid_argument("Exterior-coordinate construction requires a primary mesh with a nodal finite-element space.");
+        }
+
+        for (int element_id = 0; element_id < mesh.GetNE(); ++element_id) {
+            if (mesh.GetElementGeometry(element_id) != reference_mesh.GetElementGeometry(element_id)) {
+                throw std::invalid_argument(std::format("Primary and reference element {} have different geometries.", element_id));
+            }
+            if (mesh.GetAttribute(element_id) != reference_mesh.GetAttribute(element_id)) {
+                throw std::invalid_argument(std::format("Primary and reference element {} have different attributes.", element_id));
+            }
+        }
+
+        auto field = std::make_unique<ScalarMeshField>();
+        const mfem::FiniteElementCollection* collection = mesh.GetNodalFESpace()->FEColl();
+        field->space = std::make_unique<mfem::FiniteElementSpace>(&mesh, collection);
+        field->values = std::make_unique<mfem::GridFunction>(field->space.get());
+        *field->values = 0.0;
+
+        const int scalar_dofs = field->space->GetNDofs();
+        std::vector<bool> processed(static_cast<size_t>(scalar_dofs), false);
+        mfem::Array<int> element_dofs;
+        mfem::Vector logical_position(reference_mesh.SpaceDimension());
+
+        const double consistency_tolerance = 4096.0 * std::numeric_limits<double>::epsilon();
+
+        for (int element_id = 0; element_id < mesh.GetNE(); ++element_id) {
+            const mfem::FiniteElement& element = *field->space->GetFE(element_id);
+            const mfem::IntegrationRule& nodes = element.GetNodes();
+            mfem::ElementTransformation* reference_transformation = reference_mesh.GetElementTransformation(element_id);
+
+            if (reference_transformation == nullptr) throw std::runtime_error(std::format("Reference element {} has no element transformation.", element_id));
+
+            field->space->GetElementDofs(element_id, element_dofs);
+            if (nodes.GetNPoints() != element_dofs.Size()) {
+                throw std::runtime_error(std::format("Element {} has {} nodal points but {} scalar DOFs.", element_id, nodes.GetNPoints(), element_dofs.Size()));
+            }
+
+            for (int local_dof = 0; local_dof < element_dofs.Size(); ++local_dof) {
+                const int encoded_dof = element_dofs[local_dof];
+                const int global_dof = encoded_dof >= 0 ? encoded_dof : -1 - encoded_dof;
+
+                if (global_dof < 0 || global_dof >= scalar_dofs) {
+                    throw std::runtime_error(std::format("Element {} references invalid scalar DOF {}.", element_id, global_dof));
+                }
+
+                reference_transformation->Transform(nodes.IntPoint(local_dof), logical_position);
+                const double coordinate = compute_exterior_coordinate(logical_position, mesh.GetAttribute(element_id), config);
+
+                if (processed[static_cast<size_t>(global_dof)]) {
+                    const double existing_coordinate = (*field->values)(global_dof);
+                    if (std::abs(existing_coordinate - coordinate) > consistency_tolerance) {
+                        throw std::runtime_error(std::format("Exterior coordinate is inconsistent at shared scalar DOF {}: existing value {}, new value {} from element {}.", global_dof, existing_coordinate, coordinate, element_id));
+                    }
+                    continue;
+                }
+
+                (*field->values)(global_dof) = coordinate;
+                processed[static_cast<size_t>(global_dof)] = true;
+            }
+        }
+
+        for (int dof = 0; dof < scalar_dofs; ++dof) {
+            if (!processed[static_cast<size_t>(dof)]) throw std::runtime_error(std::format("Exterior-coordinate scalar DOF {} was not assigned.", dof));
+            const double coordinate = (*field->values)(dof);
+            if (!std::isfinite(coordinate) || coordinate < 0.0 || coordinate > 1.0) {
+                throw std::runtime_error(std::format("Exterior-coordinate scalar DOF {} has invalid value {}.", dof, coordinate));
+            }
+        }
+
+        return field;
     }
 }

src/lib/topology/mapping.cpp

@@ -1,6 +1,64 @@
 #include "stroid/topology/mapping.h"
+#include "stroid/exceptions/exceptions.h"
 #include <cmath>
 #include <algorithm>
+#include <array>
+#include <utility>
+#include <format>
+#include <string>
+
+namespace {
+    template<int n, int k>
+    consteval int nCr() {
+        if constexpr (k > n) {
+            return 0;
+        } else {
+            if constexpr (constexpr int kk = (k * 2 > n) ? (n - k) : k; kk == 0) {
+                return 1;
+            } else {
+                int result = n;
+
+                for (int i = 2; i <= kk; ++i) {
+                    result *= (n - i + 1);
+                    result /= i;
+                }
+
+                return result;
+            }
+        }
+    }
+
+    template <int n>
+    double GeneralizedSmoothstep(const double x) {
+        if (x <= 0.0) return 0.0;
+        if (x >= 1.0) return 1.0;
+
+        double sum = 0.0;
+
+        auto compute_term = [&]<std::size_t k>(std::integral_constant<std::size_t, k>) {
+            return nCr<n + k, k>() * std::pow(1.0 - x, k);
+        };
+
+        auto unroller = [&]<std::size_t... ks>(std::index_sequence<ks...>) {
+            return (compute_term(std::integral_constant<std::size_t, ks>{}) + ...);
+        };
+
+        sum = unroller(std::make_index_sequence<n + 1>{});
+        return sum * std::pow(x, n + 1);
+    }
+
+    template <std::size_t... Is>
+    constexpr auto make_smoothstep_dispatch_table(std::index_sequence<Is...>) {
+        return std::array<double(*)(double), sizeof...(Is)>{
+            &GeneralizedSmoothstep<Is + 1>...
+        };
+    }
+
+    constexpr int MAX_SMOOTHSTEP_ORDER = 10;
+    constexpr auto smoothstep_dispatch = make_smoothstep_dispatch_table(
+        std::make_index_sequence<MAX_SMOOTHSTEP_ORDER>{}
+    );
+}
 
 namespace stroid::topology {
     void ApplyEquiangular(mfem::Vector &pos) {
@@ -29,48 +87,64 @@ namespace stroid::topology {
     }
 
     void ApplySpheroidal(mfem::Vector &pos, const fourdst::config::Config<config::MeshConfig> &config) {
-        pos(2) *= (1.0 - config->flattening);
+        pos(2) *= (1.0 - config->flattening.value());
     }
 
     void TransformPoint(mfem::Vector &pos, const fourdst::config::Config<config::MeshConfig> &config, int attribute_id) {
-        double l_inf = 0.0;
-        for (int i = 0; i < pos.Size(); ++i) {
-            l_inf = std::max(l_inf, std::abs(pos(i)));
+        double X = pos(0);
+        double Y = pos(1);
+        double Z = pos(2);
+
+        double maxAbs = std::max({std::abs(X), std::abs(Y), std::abs(Z)});
+        const bool multi_block = config->core_mapping.value_or("spherified") == "multi_block";
+        const double inner_radius = config->r_core.value() / 2.0;
+        if (multi_block && maxAbs <= inner_radius) {
+            pos /= std::sqrt(3.0);
+            ApplySpheroidal(pos, config);
+            return;
         }
+        if (!multi_block && maxAbs < 1e-14) return;
 
-        if (l_inf < config->r_instability) return;
+        double cx = X / maxAbs;
+        double cy = Y / maxAbs;
+        double cz = Z / maxAbs;
 
-        // Gnomonic projection
-        const double r_log = pos.Norml2();
-        mfem::Vector unit_dir = pos;
-        unit_dir /= r_log;
+        double sx = cx * std::sqrt(1.0 - cy*cy/2.0 - cz*cz/2.0 + cy*cy*cz*cz/3.0);
+        double sy = cy * std::sqrt(1.0 - cx*cx/2.0 - cz*cz/2.0 + cx*cx*cz*cz/3.0);
+        double sz = cz * std::sqrt(1.0 - cx*cx/2.0 - cy*cy/2.0 + cx*cx*cy*cy/3.0);
 
-        ApplyEquiangular(unit_dir);
-        unit_dir /= unit_dir.Norml2(); // Re-normalize
+        mfem::Vector unit_dir(3);
+        unit_dir(0) = sx;
+        unit_dir(1) = sy;
+        unit_dir(2) = sz;
 
-        if (l_inf <= config->r_core) {
-            const double t = l_inf / config->r_core;
-            double alpha = std::pow(t, config->core_steepness);
-
-            // Smoothstep function to apply C1 continuity
-            alpha = alpha * alpha * (3.0 - 2.0 * alpha);
-
-            mfem::Vector pos_cartesian = pos;
-            mfem::Vector pos_spherical = unit_dir;
-
-            pos_spherical *= l_inf;
-
-
-            for (int d = 0; d < pos.Size(); ++d) {
-                pos(d) = (1.0 - alpha) * pos_cartesian(d) + alpha * pos_spherical(d);
+        if (maxAbs <= config->r_core.value()) {
+            if (multi_block) {
+                const double t = (maxAbs - inner_radius) / inner_radius;
+                const double inner_scale = inner_radius / std::sqrt(3.0);
+                pos(0) = (1.0 - t) * inner_scale * cx + t * config->r_core.value() * sx;
+                pos(1) = (1.0 - t) * inner_scale * cy + t * config->r_core.value() * sy;
+                pos(2) = (1.0 - t) * inner_scale * cz + t * config->r_core.value() * sz;
+                ApplySpheroidal(pos, config);
+                return;
             }
+            double nx = X / config->r_core.value();
+            double ny = Y / config->r_core.value();
+            double nz = Z / config->r_core.value();
+
+            pos(0) = nx * std::sqrt(1.0 - ny*ny/2.0 - nz*nz/2.0 + ny*ny*nz*nz/3.0);
+            pos(1) = ny * std::sqrt(1.0 - nx*nx/2.0 - nz*nz/2.0 + nx*nx*nz*nz/3.0);
+            pos(2) = nz * std::sqrt(1.0 - nx*nx/2.0 - ny*ny/2.0 + nx*nx*ny*ny/3.0);
+
+            pos *= config->r_core.value();
 
             ApplySpheroidal(pos, config);
             return;
         }
-        if (l_inf <= config->r_star) {
-            const double xi = (l_inf - config->r_core) / (config->r_star - config->r_core);
-            const double r_phys = config->r_core + xi * (config->r_star - config->r_core);
+
+        if (maxAbs <= config->r_star.value()) {
+            const double xi = (maxAbs - config->r_core.value()) / (config->r_star.value() - config->r_core.value());
+            const double r_phys = config->r_core.value() + xi * (config->r_star.value() - config->r_core.value());
 
             pos = unit_dir;
             pos *= r_phys;
@@ -78,9 +152,37 @@ namespace stroid::topology {
             ApplySpheroidal(pos, config);
         } else {
             pos = unit_dir;
-            pos *= l_inf;
+            pos *= maxAbs;
 
             ApplySpheroidal(pos, config);
         }
     }
+
+    double ComputeExteriorCoordinate(
+        const mfem::Vector& logical_position,
+        const int attribute,
+        const fourdst::config::Config<config::MeshConfig>& config
+    ) {
+        if (!config->include_external_domain.value() || attribute != static_cast<int>(config->vacuum_id.value())) return 0.0;
+
+        const double logical_radius = std::max({
+            std::abs(logical_position(0)),
+            std::abs(logical_position(1)),
+            std::abs(logical_position(2))
+        });
+
+        const double r_star = config->r_star.value();
+        const double r_infinity = config->r_infinity.value();
+        const double coordinate = (logical_radius - r_star) / (r_infinity - r_star);
+        constexpr double tolerance = 64.0 * std::numeric_limits<double>::epsilon();
+
+        if (coordinate < -tolerance || coordinate > 1.0 + tolerance) {
+            throw std::runtime_error("Logical exterior coordinate lies outside [0, 1].");
+        }
+
+        if (std::abs(coordinate) <= tolerance) return 0.0;
+        if (std::abs(coordinate - 1.0) <= tolerance) return 1.0;
+        return coordinate;
+    }
+
 }

src/lib/topology/optimize.cpp

@@ -0,0 +1,238 @@
+#include "mfem.hpp"
+
+#include <thread>
+#include <atomic>
+#include <chrono>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <algorithm>
+
+#include "stroid/topology/optimize.h"
+
+
+#include <clocale>
+#include <cstdlib>
+#include <cstring>
+#include <string>
+
+namespace stroid::utils::term_support {
+    inline bool locale_name_looks_utf8(const char* localeName) {
+        if (localeName == nullptr) return false;
+
+        const std::string localeString(localeName);
+
+        return localeString.find("UTF-8") != std::string::npos ||
+               localeString.find("utf-8") != std::string::npos ||
+               localeString.find("utf8")  != std::string::npos ||
+               localeString.find("UTF8")  != std::string::npos;
+    }
+
+    inline bool unicode_output_is_usable() {
+        const char* ctypeLocale = std::setlocale(LC_CTYPE, "");
+        if (locale_name_looks_utf8(ctypeLocale)) {
+            return true;
+        }
+
+        const char* lcAllEnv = std::getenv("LC_ALL");
+        if (locale_name_looks_utf8(lcAllEnv)) {
+            return true;
+        }
+
+        const char* lcCtypeEnv = std::getenv("LC_CTYPE");
+        if (locale_name_looks_utf8(lcCtypeEnv)) {
+            return true;
+        }
+
+        const char* langEnv = std::getenv("LANG");
+        if (locale_name_looks_utf8(langEnv)) {
+            return true;
+        }
+
+        return false;
+    }
+}
+
+namespace stroid::topology {
+class TMOPProgressBar : public mfem::IterativeSolverMonitor {
+    private:
+        double r0_ = -1.0;
+        double rtol_;
+        int bar_width_;
+
+        std::atomic<bool> done_{false};
+        std::atomic<double> progress_{0.0};
+        std::atomic<int> iter_{0};
+        std::atomic<double> res_{0.0};
+
+        std::thread spinner_thread_;
+
+        void Spin() {
+            std::vector<std::string> spin_chars;
+            if (!utils::term_support::unicode_output_is_usable()) {
+                spin_chars= {"|", "/", "-", "\\"};
+            } else {
+                spin_chars = {"▉", "▊", "▋", "▌", "▍", "▎", "▏", "▎", "▍", "▌", "▋", "▊", "▉"};
+            }
+            int spin_idx = 0;
+
+            while (!done_.load()) {
+                Draw(spin_chars[spin_idx]);
+                spin_idx = (spin_idx + 1) % spin_chars.size();
+                std::this_thread::sleep_for(std::chrono::milliseconds(100));
+            }
+        }
+
+        void Draw(const std::string& spinner) {
+            const double p = progress_.load();
+            const int pos = static_cast<int>(bar_width_ * p);
+
+            std::cout << "\r[" << spinner << "] TMOP Relaxation [";
+            for (int i = 0; i < bar_width_; ++i) {
+                if (i < pos) std::cout << "=";
+                else if (i == pos) std::cout << ">";
+                else std::cout << " ";
+            }
+            std::cout << "] " << std::setw(3) << static_cast<int>(p * 100.0) << "% "
+                      << "(Iter: " << std::setw(2) << iter_.load()
+                      << ", Res: " << std::scientific << std::setprecision(2) << res_.load() << ") " << std::flush;
+        }
+
+    public:
+        TMOPProgressBar(double rel_tol, int width = 50)
+            : rtol_(rel_tol), bar_width_(width) {
+            spinner_thread_ = std::thread(&TMOPProgressBar::Spin, this);
+        }
+
+        ~TMOPProgressBar() override {
+            if (!done_.load()) {
+                done_ = true;
+                if (spinner_thread_.joinable()) {
+                    spinner_thread_.join();
+                }
+            }
+        }
+
+        void MonitorResidual(int it, double norm, const mfem::Vector &r, bool final) override {
+            if (it == 0 || r0_ < 0.0) {
+                r0_ = norm;
+            }
+
+            iter_ = it;
+            res_ = norm;
+
+            double p = 0.0;
+            const double target_norm = r0_ * rtol_;
+
+            if (norm <= target_norm || final) {
+                p = 1.0;
+            } else if (norm < r0_ && r0_ > 0.0 && target_norm > 0.0) {
+                const double log_start = std::log10(r0_);
+                const double log_current = std::log10(norm);
+                const double log_target = std::log10(target_norm);
+                p = (log_start - log_current) / (log_start - log_target);
+                p = std::clamp(p, 0.0, 1.0);
+            }
+
+            progress_ = p;
+
+            if (final) {
+                done_ = true;
+                if (spinner_thread_.joinable()) {
+                    spinner_thread_.join();
+                }
+
+                Draw("*");
+                std::cout << std::endl;
+            }
+        }
+    };
+    void ApplyTMOP(mfem::Mesh &mesh, const fourdst::config::Config<config::MeshConfig> &config) {
+        const mfem::FiniteElementSpace* cfes = mesh.GetNodalFESpace();
+        mfem::FiniteElementSpace* fes = const_cast<mfem::FiniteElementSpace*>(cfes);
+
+        if (!fes) {
+            std::cerr << "Error: Mesh has no nodal finite element space. Call PromoteToHighOrder first." << std::endl;
+            return;
+        }
+
+        const int max_bdr_attr = mesh.bdr_attributes.Size() > 0 ? mesh.bdr_attributes.Max() : 0;
+        mfem::Array<int> ess_bdr(max_bdr_attr);
+        ess_bdr = 0.0;
+
+        if (max_bdr_attr >= config->surface_bdr_id.value()) {
+            ess_bdr[config->surface_bdr_id.value() - 1] = 1;
+        }
+
+        if (config->include_external_domain.value() && max_bdr_attr >= config->inf_bdr_id.value()) {
+            ess_bdr[config->inf_bdr_id.value() - 1] = 1;
+        }
+
+        mfem::Array<int> ess_tdof_list;
+        fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+        mfem::TMOP_QualityMetric* metric = new mfem::TMOP_Metric_302();
+        mfem::TargetConstructor* target_c = new mfem::TargetConstructor(mfem::TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE);
+        target_c->SetNodes(*mesh.GetNodes());
+        mfem::TMOP_Integrator* tmop_integrator = new mfem::TMOP_Integrator(metric, target_c);
+
+        mfem::NonlinearForm a(fes);
+        a.AddDomainIntegrator(tmop_integrator);
+        a.SetEssentialTrueDofs(ess_tdof_list);
+
+        mfem::GridFunction* nodes = mesh.GetNodes();
+        mfem::Vector x(*nodes);
+        mfem::Vector b(a.Height());
+        b = 0.0;
+
+        mfem::MINRESSolver minres;
+        minres.SetMaxIter(750);
+        minres.SetRelTol(1e-5);
+        minres.SetAbsTol(0.0);
+        minres.SetPrintLevel(-1);
+
+        mfem::DSmoother jacobi(1, 1.0, 1);
+        jacobi.SetPositiveDiagonal(true);
+        minres.SetPreconditioner(jacobi);
+
+        const int quad_order = 2 * fes->GetMaxElementOrder() + 3;
+        const mfem::IntegrationRule &ir = mfem::IntRules.Get(mesh.GetTypicalElementGeometry(), quad_order);
+
+        double min_detJ = std::numeric_limits<double>::infinity();
+        for (int i = 0; i < mesh.GetNE(); i++) {
+            mfem::ElementTransformation *T = mesh.GetElementTransformation(i);
+            for (int j = 0; j < ir.GetNPoints(); j++) {
+                T->SetIntPoint(&ir.IntPoint(j));
+                min_detJ = std::min(min_detJ, T->Jacobian().Det());
+            }
+        }
+
+        constexpr double newton_rtol = 1e-8;
+        mfem::TMOPNewtonSolver newton(ir, 0);
+        newton.SetPreconditioner(minres);
+        newton.SetOperator(a);
+        newton.SetMaxIter(50);
+        newton.SetRelTol(newton_rtol);
+        newton.SetAbsTol(0.0);
+        newton.SetMinDetPtr(&min_detJ);
+        newton.SetPrintLevel(0);
+
+        TMOPProgressBar progress_bar(newton_rtol);
+        newton.SetMonitor(progress_bar);
+
+        std::cout << "Applying TMOP optimization to mesh. Note this may take a long time. Depending on your mesh resolution expect to wait up to the order of 10s of minutes..." << std::endl;
+        newton.Mult(b, x);
+        *nodes = x;
+
+        mesh.NodesUpdated();
+
+        delete metric;
+        delete target_c;
+    }
+
+    void OptimizeMesh(mfem::Mesh& mesh, const fourdst::config::Config<config::MeshConfig> &cfg) {
+        if (cfg->optimization_methods.has_value() && cfg->optimization_methods.value().tmop.has_value() && cfg->optimization_methods.value().tmop.value()) {
+            ApplyTMOP(mesh, cfg);
+        }
+    }
+}

src/lib/topology/topology.cpp

@@ -1,6 +1,8 @@
 #include "mfem.hpp"
 #include <vector>
 #include <memory>
+#include <cmath>
+#include <stdexcept>
 
 #include "stroid/config/config.h"
 #include "fourdst/config/config.h"
@@ -8,9 +10,31 @@
 namespace stroid::topology {
 
     std::unique_ptr<mfem::Mesh> BuildSkeleton(const fourdst::config::Config<config::MeshConfig> & config) {
-        int nVert = config->include_external_domain ? 24 : 16;
-        int nElem = config->include_external_domain ? 13 : 7;
-        int nBev  = config->include_external_domain ? 12 : 6;
+        const std::string core_mapping = config->core_mapping.value_or("spherified");
+        if (core_mapping != "spherified" && core_mapping != "multi_block") {
+            throw std::invalid_argument("Unknown core_mapping: " + core_mapping);
+        }
+
+        const bool multi_block = core_mapping == "multi_block";
+        const bool include_external_domain = config->include_external_domain.value_or(true);
+        if (multi_block) {
+            const double r_core = config->r_core.value();
+            const double r_star = config->r_star.value();
+            const double r_infinity = config->r_infinity.value_or(6.0);
+            const double flattening = config->flattening.value();
+            if (!std::isfinite(r_core) || !std::isfinite(r_star) || r_core <= 0.0 || r_star <= r_core ||
+                (include_external_domain && (!std::isfinite(r_infinity) || r_infinity <= r_star))) {
+                throw std::invalid_argument("multi_block requires 0 < r_core < r_star < r_infinity (when external).");
+            }
+            if (!std::isfinite(flattening) || flattening >= 1.0) {
+                throw std::invalid_argument("multi_block requires finite flattening < 1.");
+            }
+        }
+
+        const int offset = multi_block ? 8 : 0;
+        int nVert = (include_external_domain ? 24 : 16) + offset;
+        int nElem = (include_external_domain ? 13 : 7) + (multi_block ? 6 : 0);
+        int nBev  = include_external_domain ? 12 : 6;
 
         auto mesh = std::make_unique<mfem::Mesh>(3, nVert, nElem, nBev, 3);
 
@@ -21,14 +45,17 @@ namespace stroid::topology {
                         mesh->AddVertex(x, y, z);
         };
 
-        add_box(config->r_core);
-        add_box(config->r_star);
-        if (config->include_external_domain) {
-            add_box(config->r_infinity);
+        if (multi_block) {
+            add_box(config->r_core.value() / 2.0);
+        }
+        add_box(config->r_core.value());
+        add_box(config->r_star.value());
+        if (include_external_domain) {
+            add_box(config->r_infinity.value());
         }
 
         const int core_v[8] = {0, 1, 3, 2, 4, 5, 7, 6};
-        mesh->AddHex(core_v, config->core_id);
+        mesh->AddHex(core_v, config->core_id.value());
 
         std::vector<std::array<int, 8>> stellar_shells = {
             {8, 9, 11, 10, 0, 1, 3, 2},
@@ -38,11 +65,18 @@ namespace stroid::topology {
             {1, 3, 7, 5, 9, 11, 15, 13},  // +X face
             {0, 4, 6, 2, 8, 12, 14, 10}   // -X face
         };
+        if (multi_block) {
             for (const auto & shell : stellar_shells) {
-            mesh->AddHex(shell.data(), config->envelope_id);
+                mesh->AddHex(shell.data(), config->core_id.value());
+            }
+        }
+        for (const auto & shell : stellar_shells) {
+            auto vertices = shell;
+            for (auto & vertex : vertices) vertex += offset;
+            mesh->AddHex(vertices.data(), config->envelope_id.value());
         }
 
-        if (config->include_external_domain) {
+        if (include_external_domain) {
             std::vector<std::array<int, 8>> vacuum_shells;
             vacuum_shells.push_back({8, 9, 13, 12, 16, 17, 21, 20});
             vacuum_shells.push_back({9, 11, 15, 13, 17, 19, 23, 21});
@@ -51,7 +85,9 @@ namespace stroid::topology {
             vacuum_shells.push_back({12, 13, 15, 14, 20, 21, 23, 22});
             vacuum_shells.push_back({10, 11, 9, 8, 18, 19, 17, 16});
             for (const auto & shell : vacuum_shells) {
-                mesh->AddHex(shell.data(), config->vacuum_id);
+                auto vertices = shell;
+                for (auto & vertex : vertices) vertex += offset;
+                mesh->AddHex(vertices.data(), config->vacuum_id.value());
             }
         }
 
@@ -66,10 +102,12 @@ namespace stroid::topology {
         };
 
         for (const auto& bdr: surface_bdr_quads) {
-            mesh->AddBdrQuad(bdr, config->surface_bdr_id);
+            int vertices[4];
+            for (int i = 0; i < 4; ++i) vertices[i] = bdr[i] + offset;
+            mesh->AddBdrQuad(vertices, config->surface_bdr_id.value());
         }
 
-        if (config->include_external_domain) {
+        if (include_external_domain) {
             const int inf_bdr_quads[6][4] = {
                 {16, 17, 21, 20},
                 {17, 19, 23, 21},
@@ -80,13 +118,16 @@ namespace stroid::topology {
             };
 
             for (const auto& bdr: inf_bdr_quads) {
-                mesh->AddBdrQuad(bdr, config->inf_bdr_id);
+                int vertices[4];
+                for (int i = 0; i < 4; ++i) vertices[i] = bdr[i] + offset;
+                mesh->AddBdrQuad(vertices, config->inf_bdr_id.value());
             }
         }
 
         return mesh;
     }
 
+    // ReSharper disable once CppUseInternalLinkage
     void Finalize(mfem::Mesh& mesh, const fourdst::config::Config<config::MeshConfig> &config) {
         mesh.FinalizeTopology();
         mesh.Finalize();

src/lib/utils/mesh_stats.cpp

@@ -0,0 +1,485 @@
+#include "stroid/utils/mesh_stats.h"
+
+#include <cmath>
+#include <limits>
+#include <format>
+#include <string>
+#include <algorithm>
+
+namespace stroid::stats {
+    namespace {
+        double SpheroidRadius(const double ux, const double uy, const double uz,
+                              const double r_star, const double flattening) {
+            const double a = r_star;
+            const double c = r_star * (1.0 - flattening);
+            const double inv = (ux*ux + uy*uy) / (a*a) + (uz*uz) / (c*c);
+            return (inv > 0.0) ? 1.0 / std::sqrt(inv) : 0.0;
+        }
+    }
+
+    MeshStats ComputeMeshStats(const StroidMesh& sm, MeshStatFeatures features, int sample_order) {
+        MeshStats out;
+        out.computed = features;
+
+        auto mesh_or = sm.as_mesh();
+        if (!mesh_or) {
+            out.errors.push_back(mesh_or.error());
+            return out;
+        }
+        mfem::Mesh* mesh = *mesh_or;
+        if (!mesh) {
+            out.errors.push_back("StroidMesh has no stored computational mesh to compute stats against");
+            return out;   // BUGFIX: was falling through to a null deref
+        }
+
+        const auto& cfg = sm.config;
+        const double r_star     = cfg.r_star.value_or(1.0);
+        const double flattening = cfg.flattening.value_or(0.0);
+        const int    surf_bdr   = static_cast<int>(cfg.surface_bdr_id.value_or(-99));
+        const int    inf_bdr    = static_cast<int>(cfg.inf_bdr_id.value_or(-99));
+        const int    core_id    = static_cast<int>(cfg.core_id.value_or(-99));
+        const int    env_id     = static_cast<int>(cfg.envelope_id.value_or(-99));
+        const int    vac_id     = static_cast<int>(cfg.vacuum_id.value_or(-99));
+
+        int geom_order = -99;
+        if (mesh->GetNodes()) {
+            geom_order = mesh->GetNodes()->FESpace()->GetMaxElementOrder();
+        }
+        const int sorder = (sample_order > 0) ? sample_order : (2 * geom_order + 4);
+
+        if (has_feature(features, MeshStatFeatures::CONFIG_META)) {
+            ConfigMeta meta;
+            meta.r_core             = cfg.r_core.value_or(-99.99);
+            meta.r_star             = cfg.r_star.value_or(-99.99);
+            meta.r_infinity         = cfg.r_infinity.value_or(-99.99);
+            meta.flattening         = flattening;
+            meta.geom_order         = geom_order;
+            meta.refinement_levels  = sm.refinement_levels;
+            meta.has_external_domain = cfg.include_external_domain.value_or(false);
+            out.config_meta = meta;
+        }
+
+        if (has_feature(features, MeshStatFeatures::CONFORMITY)) {
+            ConformityStats conformity;
+            conformity.conforming = mesh->Conforming();
+            conformity.n_nonconforming_faces = conformity.conforming ? 0 : -99; // TODO: count
+            out.conformity = conformity;
+        }
+
+        // ============================ SURFACE PASS ============================
+        const bool needs_surface =
+            has_feature(features, MeshStatFeatures::RADIUS) ||
+            has_feature(features, MeshStatFeatures::AXES) ||
+            has_feature(features, MeshStatFeatures::ELLIPTICITY) ||
+            has_feature(features, MeshStatFeatures::BOWING) ||
+            has_feature(features, MeshStatFeatures::VOLUME_AREA);
+
+        if (needs_surface) {
+            double r_min = std::numeric_limits<double>::max();
+            double r_max = std::numeric_limits<double>::lowest();
+            double r_sum = 0.0, r_sum_sq = 0.0;
+            double a_eq = 0.0, c_pol = 0.0;
+            double bow_in = 0.0, bow_out = 0.0, bow_sum_sq = 0.0, bow_worst_r = 0.0;
+            double bow_worst_mag = 0.0;
+            double area = 0.0;
+            long   n_samples = 0;
+            mfem::Vector phys;
+
+            for (int b = 0; b < mesh->GetNBE(); ++b) {
+                if (mesh->GetBdrAttribute(b) != surf_bdr) continue;
+                mfem::ElementTransformation* T = mesh->GetBdrElementTransformation(b);
+                const mfem::IntegrationRule& ir = mfem::IntRules.Get(T->GetGeometryType(), sorder);
+
+                for (int q = 0; q < ir.GetNPoints(); ++q) {
+                    const mfem::IntegrationPoint& ip = ir.IntPoint(q);
+                    T->SetIntPoint(&ip);
+                    T->Transform(ip, phys);
+
+                    const double x = phys(0), y = phys(1), z = phys(2);
+                    const double r = phys.Norml2();
+                    const double rho = std::sqrt(x*x + y*y);
+                    const double w = T->Weight() * ip.weight;
+
+                    r_min = std::min(r_min, r);
+                    r_max = std::max(r_max, r);
+                    r_sum += r; r_sum_sq += r*r;
+                    a_eq = std::max(a_eq, rho);
+                    c_pol = std::max(c_pol, std::abs(z));
+                    area += w;
+                    ++n_samples;
+
+                    if (has_feature(features, MeshStatFeatures::BOWING) && r > 1e-14) {
+                        const double rt = SpheroidRadius(x/r, y/r, z/r, r_star, flattening);
+                        const double dev = r - rt;
+                        bow_in  = std::min(bow_in, dev);
+                        bow_out = std::max(bow_out, dev);
+                        bow_sum_sq += dev * dev;
+                        if (std::abs(dev) > bow_worst_mag) {
+                            bow_worst_mag = std::abs(dev);
+                            bow_worst_r = r;
+                        }
+                    }
+                }
+            }
+
+            if (n_samples == 0) {
+                out.warnings.push_back("No samples were collected from the surface boundary. Check that the "
+                                       "surface boundary ID is correct. Lacking a surface pass prevents the "
+                                       "following from reporting accurate results: "
+                                       "[RADIUS, AXES, ELLIPTICITY, BOWING, VOLUME_AREA]");
+            } else {
+                if (has_feature(features, MeshStatFeatures::RADIUS)) {
+                    const double mean = r_sum / n_samples;
+                    const double var  = std::max(0.0, r_sum_sq / n_samples - mean * mean);
+                    out.radius = RadiusStats{
+                        .min = r_min, .max = r_max, .mean = mean,
+                        .stddev = std::sqrt(var), .n_samples = n_samples
+                    };
+                }
+                if (has_feature(features, MeshStatFeatures::AXES)) {
+                    out.axes = AxisStats{ .semi_major = a_eq, .semi_minor = c_pol };
+                }
+                if (has_feature(features, MeshStatFeatures::ELLIPTICITY)) {
+                    out.ellipticity = EllipticityStats{
+                        .flattening        = (a_eq > 0) ? (a_eq - c_pol) / a_eq : 0.0,
+                        .polar_equatorial  = (a_eq > 0) ? c_pol / a_eq : 1.0,
+                        .radius_uniformity = (r_max > 0) ? r_min / r_max : 1.0,
+                    };
+                }
+                if (has_feature(features, MeshStatFeatures::BOWING)) {
+                    out.bowing = BowingStats{
+                        .max_inward = bow_in, .max_outward = bow_out,
+                        .rms = std::sqrt(bow_sum_sq / n_samples), .worst_at_radius = bow_worst_r,
+                    };
+                }
+                if (has_feature(features, MeshStatFeatures::VOLUME_AREA)) {
+                    const double a = r_star, c = r_star * (1.0 - flattening);
+                    out.volume = VolumeAreaStats{
+                        .surface_area = area,
+                        .analytic_area = (flattening == 0) ? 4.0 * M_PI * a * a : -99.99
+                    };
+                }
+            }
+        }
+
+        // ============================ VOLUME PASS ============================
+        const bool needs_volume =
+            has_feature(features, MeshStatFeatures::VOLUME_AREA) ||
+            has_feature(features, MeshStatFeatures::JACOBIAN) ||
+            has_feature(features, MeshStatFeatures::ELEMENT_COUNT) ||
+            has_feature(features, MeshStatFeatures::MESH_SIZE) ||
+            has_feature(features, MeshStatFeatures::CENTROID) ||
+            has_feature(features, MeshStatFeatures::BOUNDING_BOX);   // BUGFIX: was a dangling ';'
+
+        if (needs_volume) {
+            const bool need_bbox = has_feature(features, MeshStatFeatures::BOUNDING_BOX);
+            const bool need_jac  = has_feature(features, MeshStatFeatures::JACOBIAN);
+
+            // Bounding-box accumulators (seeded inverted so empty regions stay invalid).
+            double cxmin=+std::numeric_limits<double>::max(), cxmax=-std::numeric_limits<double>::max();
+            double cymin=cxmin, cymax=cxmax, czmin=cxmin, czmax=cxmax;                            // core
+            double sxmin=cxmin, sxmax=cxmax, symin=cxmin, symax=cxmax, szmin=cxmin, szmax=cxmax;  // stellar
+            double vxmin=cxmin, vxmax=cxmax, vymin=cxmin, vymax=cxmax, vzmin=cxmin, vzmax=cxmax;  // vacuum
+            long n_core_box=0, n_stel_box=0, n_vac_box=0;
+            mfem::Vector bphys;
+
+            // Per-region Jacobian accumulators.
+            struct JacAccum {
+                double detJ_min = std::numeric_limits<double>::max();
+                double detJ_max = -std::numeric_limits<double>::max();
+                double min_ratio = 1.0;
+                long   n_flipped = 0;
+                long   n_elem = 0;
+                double worst_ratio_r = -1.0;
+                double min_detJ_r = -1.0;
+            };
+            JacAccum all_acc, stel_acc, vac_acc;
+
+            auto jac_update = [](JacAccum& a, double dmin, double dmax, bool flip, double r) {
+                ++a.n_elem;
+                if (dmin < a.detJ_min) { a.detJ_min = dmin; a.min_detJ_r = r; }
+                if (dmax > a.detJ_max)   a.detJ_max = dmax;
+                if (dmax > 1e-30) {
+                    const double ratio = dmin / dmax;
+                    if (ratio < a.min_ratio) { a.min_ratio = ratio; a.worst_ratio_r = r; }
+                }
+                if (flip) ++a.n_flipped;
+            };
+
+            double vol = 0.0, cx = 0.0, cy = 0.0, cz = 0.0;
+            double h_min = std::numeric_limits<double>::max();
+            double h_max = -std::numeric_limits<double>::max();
+            double h_sum = 0.0, h_sum_sq = 0.0;
+            long   n_core = 0, n_env = 0, n_vac = 0, n_other = 0;
+            mfem::Vector phys;
+
+            for (int e = 0; e < mesh->GetNE(); ++e) {
+                const int attr = mesh->GetAttribute(e);
+                if      (attr == core_id) ++n_core;
+                else if (attr == env_id)  ++n_env;
+                else if (attr == vac_id)  ++n_vac;
+                else                      ++n_other;
+
+                const bool stellar = (attr == core_id || attr == env_id);
+
+                if (has_feature(features, MeshStatFeatures::MESH_SIZE) && stellar) {
+                    const double h = mesh->GetElementSize(e);
+                    h_min = std::min(h_min, h);
+                    h_max = std::max(h_max, h);
+                    h_sum += h; h_sum_sq += h*h;
+                }
+
+                mfem::ElementTransformation* T = mesh->GetElementTransformation(e);
+                const mfem::IntegrationRule& ir = mfem::IntRules.Get(T->GetGeometryType(), sorder);
+
+                double e_detmin = std::numeric_limits<double>::max();
+                double e_detmax = -std::numeric_limits<double>::max();
+                bool   e_flip = false;
+
+                for (int q = 0; q < ir.GetNPoints(); ++q) {
+                    const mfem::IntegrationPoint& ip = ir.IntPoint(q);
+                    T->SetIntPoint(&ip);
+
+                    if (need_bbox) {
+                        T->Transform(ip, bphys);
+                        const double X = bphys(0), Y = bphys(1), Z = bphys(2);
+                        if (attr == core_id) {
+                            cxmin=std::min(cxmin,X); cxmax=std::max(cxmax,X);
+                            cymin=std::min(cymin,Y); cymax=std::max(cymax,Y);
+                            czmin=std::min(czmin,Z); czmax=std::max(czmax,Z);
+                            ++n_core_box;
+                        }
+                        if (stellar) {   // stellar = core U envelope
+                            sxmin=std::min(sxmin,X); sxmax=std::max(sxmax,X);
+                            symin=std::min(symin,Y); symax=std::max(symax,Y);
+                            szmin=std::min(szmin,Z); szmax=std::max(szmax,Z);
+                            ++n_stel_box;
+                        }
+                        if (attr == vac_id) {
+                            vxmin=std::min(vxmin,X); vxmax=std::max(vxmax,X);
+                            vymin=std::min(vymin,Y); vymax=std::max(vymax,Y);
+                            vzmin=std::min(vzmin,Z); vzmax=std::max(vzmax,Z);
+                            ++n_vac_box;
+                        }
+                    }
+
+                    const double dJ = T->Jacobian().Det();
+                    e_detmin = std::min(e_detmin, dJ);
+                    e_detmax = std::max(e_detmax, dJ);
+                    if (dJ < 0.0) e_flip = true;
+
+                    if (stellar && (has_feature(features, MeshStatFeatures::VOLUME_AREA) ||
+                                    has_feature(features, MeshStatFeatures::CENTROID))) {
+                        const double w = std::abs(dJ) * ip.weight;
+                        vol += w;
+                        if (has_feature(features, MeshStatFeatures::CENTROID)) {
+                            T->Transform(ip, phys);
+                            cx += w*phys(0); cy += w*phys(1); cz += w*phys(2);
+                        }
+                    }
+                }
+
+                if (need_jac) {
+                    // Representative element radius (center) for locating the worst element.
+                    const mfem::IntegrationPoint& cip = mfem::Geometries.GetCenter(T->GetGeometryType());
+                    T->SetIntPoint(&cip);
+                    mfem::Vector cpt;
+                    T->Transform(cip, cpt);
+                    const double er = cpt.Norml2();
+
+                    jac_update(all_acc, e_detmin, e_detmax, e_flip, er);
+                    if (stellar)             jac_update(stel_acc, e_detmin, e_detmax, e_flip, er);
+                    else if (attr == vac_id) jac_update(vac_acc,  e_detmin, e_detmax, e_flip, er);
+                }
+            }
+
+            if (has_feature(features, MeshStatFeatures::ELEMENT_COUNT)) {
+                out.element_counts = ElementCounts{
+                    .total = n_core + n_env + n_vac + n_other,
+                    .core = n_core, .envelope = n_env, .vacuum = n_vac, .other = n_other,
+                    .n_vertices = mesh->GetNV(),
+                };
+            }
+
+            if (need_jac) {
+                auto finalize = [](const JacAccum& a) {
+                    JacobianStats j;
+                    j.n_elements           = a.n_elem;
+                    j.detJ_min             = (a.n_elem > 0) ? a.detJ_min : 0.0;
+                    j.detJ_max             = (a.n_elem > 0) ? a.detJ_max : 0.0;
+                    j.min_detJ_ratio       = a.min_ratio;
+                    j.n_flipped            = a.n_flipped;
+                    j.worst_ratio_at_radius = a.worst_ratio_r;
+                    j.detJ_min_at_radius    = a.min_detJ_r;
+                    return j;
+                };
+                out.jacobian = finalize(all_acc);
+                if (stel_acc.n_elem > 0) out.jacobian_stellar = finalize(stel_acc);
+                if (vac_acc.n_elem  > 0) out.jacobian_vacuum  = finalize(vac_acc);
+            }
+
+            if (has_feature(features, MeshStatFeatures::MESH_SIZE)) {
+                const long ns = n_core + n_env;
+                const double mean = (ns > 0) ? h_sum / ns : 0.0;
+                const double var  = (ns > 0) ? std::max(0.0, h_sum_sq / ns - mean * mean) : 0.0;
+                out.mesh_size = MeshSizeStats{
+                    .h_min = h_min, .h_max = h_max, .h_mean = mean, .h_stddev = std::sqrt(var),
+                };
+            }
+
+            if (has_feature(features, MeshStatFeatures::VOLUME_AREA)) {
+                if (!out.volume) out.volume.emplace();   // BUGFIX: surface pass may not have created it
+                out.volume->stellar_volume = vol;
+                const double a = r_star, c = r_star * (1.0 - flattening);
+                out.volume->analytic_volume = (4.0 / 3.0) * M_PI * a * a * c;
+            }
+
+            if (has_feature(features, MeshStatFeatures::CENTROID)) {
+                if (vol <= 0) {
+                    out.warnings.push_back("Stellar volume is zero or negative, cannot compute centroid.");
+                } else {
+                    const double ccx = cx / vol, ccy = cy / vol, ccz = cz / vol;  // BUGFIX: normalize
+                    out.centroid = CentroidStats{
+                        .x = ccx, .y = ccy, .z = ccz,
+                        .offset = std::sqrt(ccx*ccx + ccy*ccy + ccz*ccz)
+                    };
+                }
+            }
+
+            if (need_bbox) {
+                BoundingBoxStats bb;
+                auto fill = [](BoundingBox& box, long n,
+                               double xmn,double xmx,double ymn,double ymx,double zmn,double zmx) {
+                    if (n > 0) {
+                        box.valid = true;
+                        box.xMin=xmn; box.xMax=xmx;
+                        box.yMin=ymn; box.yMax=ymx;
+                        box.zMin=zmn; box.zMax=zmx;
+                    }
+                };
+                fill(bb.core,   n_core_box, cxmin,cxmax,cymin,cymax,czmin,czmax);
+                fill(bb.star,   n_stel_box, sxmin,sxmax,symin,symax,szmin,szmax);
+                fill(bb.vacuum, n_vac_box,  vxmin,vxmax,vymin,vymax,vzmin,vzmax);
+                out.bounding_box = bb;
+            }
+        }
+
+        // ============================ OUTER BOUND PASS ============================
+        if (has_feature(features, MeshStatFeatures::OUTER_BOUNDS)) {
+            double r_min = std::numeric_limits<double>::max();
+            double r_max = std::numeric_limits<double>::lowest();
+            double r_sum = 0.0;
+            long   n_samples = 0;
+            mfem::Vector phys;
+
+            for (int b = 0; b < mesh->GetNBE(); ++b) {
+                if (mesh->GetBdrAttribute(b) != inf_bdr) continue;
+                mfem::ElementTransformation* T = mesh->GetBdrElementTransformation(b);
+                const mfem::IntegrationRule& ir = mfem::IntRules.Get(T->GetGeometryType(), sorder);
+                for (int q = 0; q < ir.GetNPoints(); ++q) {
+                    T->SetIntPoint(&ir.IntPoint(q));
+                    T->Transform(ir.IntPoint(q), phys);
+                    const double r = phys.Norml2();
+                    r_min = std::min(r_min, r); r_max = std::max(r_max, r);
+                    r_sum += r; ++n_samples;
+                }
+            }
+
+            if (n_samples == 0) {
+                out.warnings.push_back("No samples found on the outer boundary, cannot compute outer bounds.");
+            } else {
+                out.outer_bounds = OuterBoundsStats{
+                    .min = r_min, .max = r_max, .mean = r_sum / n_samples, .n_samples = n_samples
+                };
+            }
+        }
+
+        return out;
+    }
+
+    std::string to_string(const MeshStats& s) {
+        std::string o = "MeshStats:\n";
+        auto line = [&](const std::string& l){ o += " =>" + l + "\n"; };
+
+        if (s.config_meta) {
+            const auto& m = *s.config_meta;
+            line(std::format(
+                "config: r_core={:0.4f}, r_star={:0.4f}, r_inf={:0.4f}, flattening={:0.4f}, "
+                "geometric order={}, refinement levels={}",
+                m.r_core, m.r_star, m.r_infinity, m.flattening, m.geom_order, m.refinement_levels));
+        }
+        if (s.radius) {
+            const auto& r = *s.radius;
+            line(std::format("radius: min={:.6f} max={:.6f} mean={:.6f} std={:.3E} (n={})",
+                 r.min, r.max, r.mean, r.stddev, r.n_samples));
+        }
+        if (s.axes) {
+            line(std::format("axes: semi_major(eq)={:.6f} semi_minor(pol)={:.6f}",
+                 s.axes->semi_major, s.axes->semi_minor));
+        }
+        if (s.ellipticity) {
+            const auto& e = *s.ellipticity;
+            line(std::format("ellipticity: flattening={:.5f} c/a={:.5f} r_min/r_max={:.5f}",
+                 e.flattening, e.polar_equatorial, e.radius_uniformity));
+        }
+        if (s.bowing) {
+            const auto& b = *s.bowing;
+            line(std::format("bowing: max_inward={:.3E} max_outward={:.3E} rms={:.3E}",
+                 b.max_inward, b.max_outward, b.rms));
+        }
+        if (s.conformity) {
+            line(std::format("conforming: {}", s.conformity->conforming));
+        }
+
+        auto jac_line = [&](const std::string& label, const JacobianStats& j) {
+            line(std::format(
+                "jacobian[{}]: detJ=[{:.3E},{:.3E}] min_ratio={:.3E} (@r={:.4f}) "
+                "min_detJ@r={:.4f} flipped={} n={}",
+                label, j.detJ_min, j.detJ_max, j.min_detJ_ratio, j.worst_ratio_at_radius,
+                j.detJ_min_at_radius, j.n_flipped, j.n_elements));
+        };
+        if (s.jacobian)         jac_line("all",     *s.jacobian);
+        if (s.jacobian_stellar) jac_line("stellar", *s.jacobian_stellar);
+        if (s.jacobian_vacuum)  jac_line("vacuum",  *s.jacobian_vacuum);
+
+        if (s.volume) {
+            const auto& v = *s.volume;
+            line(std::format("volume={:.6f} (analytic {:.6f})  area={:.6f}",
+                 v.stellar_volume, v.analytic_volume, v.surface_area));
+        }
+        if (s.element_counts) {
+            const auto& c = *s.element_counts;
+            line(std::format("elements: total={} core={} env={} vac={} other={} NV={}",
+                 c.total, c.core, c.envelope, c.vacuum, c.other, c.n_vertices));
+        }
+        if (s.mesh_size) {
+            line(std::format("h: min={:.4E} max={:.4E} mean={:.4E} std={:.4E}",
+                 s.mesh_size->h_min, s.mesh_size->h_max, s.mesh_size->h_mean, s.mesh_size->h_stddev));
+        }
+        if (s.bounding_box) {
+            const auto& bb = *s.bounding_box;
+            auto bline = [&](const char* nm, const BoundingBox& b){
+                if (b.valid)
+                    line(std::format("bbox[{}]: x[{:.4f},{:.4f}] y[{:.4f},{:.4f}] z[{:.4f},{:.4f}]",
+                         nm, b.xMin,b.xMax, b.yMin,b.yMax, b.zMin,b.zMax));
+                else
+                    line(std::format("bbox[{}]: <absent>", nm));
+            };
+            bline("core", bb.core);
+            bline("star", bb.star);
+            bline("vacuum", bb.vacuum);
+        }
+        if (s.outer_bounds) {
+            line(std::format("outer: min={:.4f} max={:.4f} mean={:.4f}",
+                 s.outer_bounds->min, s.outer_bounds->max, s.outer_bounds->mean));
+        }
+        if (s.centroid) {
+            line(std::format("centroid: x={:.6f} y={:.6f} z={:.6f} offset={:.6f}",
+                 s.centroid->x, s.centroid->y, s.centroid->z, s.centroid->offset));
+        }
+
+        for (const auto& w : s.warnings) line("WARNING: " + w);
+        for (const auto& e : s.errors)   line(std::format("ERROR: {}", e));
+        return o;
+    }
+}

src/lib/utils/mesh_utils.cpp

@@ -2,6 +2,8 @@
 #include "mfem.hpp"
 #include <print>
 
+#include "stroid/topology/curvilinear.h"
+
 namespace stroid::utils {
     void MarkFlippedElements(mfem::Mesh& mesh) {
         for (int i = 0; i < mesh.GetNE(); i++) {
@@ -51,4 +53,44 @@ namespace stroid::utils {
             }
         }
     }
+
+    void ExportJacobianRadialProfile(mfem::Mesh& mesh, const std::string& filename) {
+        std::ofstream ofs(filename);
+
+        if (!ofs.good()) {
+            throw std::runtime_error(std::format("Stroid: Could not open file {} for writing Jacobian radial profile", filename));
+        }
+
+        ofs << "Radius,DetJ,Attribute,ElementID\n";
+        ofs.precision(10);
+
+        const int sample_order = 2 * mesh.GetNodes()->FESpace()->GetMaxElementOrder() + 2;
+        for (int i = 0; i < mesh.GetNE(); ++i) {
+            mfem::ElementTransformation *T = mesh.GetElementTransformation(i);
+            const int attr = mesh.GetAttribute(i);
+
+            const mfem::IntegrationRule &ir = mfem::IntRules.Get(T->GetGeometryType(), sample_order);
+
+            for (int j = 0; j < ir.GetNPoints(); ++j) {
+                T->SetIntPoint(&ir.IntPoint(j));
+
+                mfem::Vector pos;
+                T->Transform(ir.IntPoint(j), pos);
+
+                const double r = pos.Norml2();
+                const double detJ = T->Jacobian().Det();
+
+                ofs << r << "," << detJ << "," << attr << ',' << i << "\n";
+            }
+        }
+        ofs.close();
+        std::println("Jacobian radial profile exported to {}", filename);
+    }
+
+    std::unique_ptr<mfem::Mesh> BuildProjected(const mfem::Mesh& reference, const fourdst::config::Config<config::MeshConfig>& cfg) {
+        auto projected = std::make_unique<mfem::Mesh>(reference);
+        topology::PromoteToHighOrder(*projected, cfg);
+        topology::ProjectMesh(*projected, cfg);
+        return projected;
+    }
 }

src/meson.build

@@ -10,25 +10,80 @@ stroid_sources = files(
     'lib/topology/curvilinear.cpp',
     'lib/topology/mapping.cpp',
     'lib/topology/topology.cpp',
+    'lib/topology/optimize.cpp',
     'lib/IO/mesh.cpp',
     'lib/utils/mesh_utils.cpp',
+    'lib/utils/mesh_stats.cpp',
+    'lib/refinement/uniform.cpp',
 )
 
-stroid_lib = static_library(
-    'libstroid',
+if get_option('build_python')
+    if host_machine.system() == 'darwin'
+        stroid_lib_rpath_args = [
+            '-Wl,-rpath,@loader_path',
+            '-Wl,-rpath,@loader_path/../../fourdst/lib',
+            '-Wl,-rpath,@loader_path/../../fourdst/lib/vendor',
+        ]
+        stroid_lib_rpath = ''
+
+        stroid_ext_rpath_args = [
+            '-Wl,-rpath,@loader_path/lib',
+            '-Wl,-rpath,@loader_path/../fourdst/lib',
+            '-Wl,-rpath,@loader_path/../fourdst/lib/vendor',
+        ]
+        stroid_ext_rpath = ''
+    else
+        stroid_lib_rpath_args = []
+        stroid_lib_rpath = '$ORIGIN:' + '$ORIGIN/../../fourdst/lib:' + '$ORIGIN/../../fourdst/lib/vendor'
+
+        stroid_ext_rpath_args = []
+        stroid_ext_rpath = '$ORIGIN/lib:' + '$ORIGIN/../fourdst/lib:' + '$ORIGIN/../fourdst/lib/vendor'
+    endif
+
+    libstroid = static_library(
+        'stroid',
+        stroid_sources,
+        include_directories: stroid_include_files,
+        dependencies: dependencies,
+        install_dir: stroid_libdir,
+        link_args: stroid_lib_rpath_args,
+        build_rpath: stroid_lib_rpath,
+        install_rpath: stroid_lib_rpath,
+    )
+else
+    libstroid = static_library(
+        'stroid',
         stroid_sources,
         include_directories: stroid_include_files,
         dependencies: dependencies,
         install: true
     )
+endif
+
+
+if get_option('build_python')
+    stroid_iface_dep = declare_dependency(
+        dependencies: dependencies
+    ).partial_dependency(compile_args: true, includes: true)
 
     stroid_dep = declare_dependency(
-    link_with: stroid_lib,
+        link_with: libstroid,
+        include_directories: stroid_include_files,
+        dependencies: [stroid_iface_dep]
+    )
+else
+    stroid_dep = declare_dependency(
+        link_with: libstroid,
         include_directories: stroid_include_files,
         dependencies: dependencies
     )
+endif
 
+meson.override_dependency('stroid', stroid_dep)
+
+message('stroid include dir: ' + stroid_includedir)
 install_subdir(
     'include/stroid',
-    install_dir: get_option('includedir') / 'stroid'
+    install_dir: stroid_includedir,
+    exclude_files: ['version.h.in']
 )

src/python/IO/bindings.cpp

@@ -0,0 +1,77 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include "bindings.h"
+
+#include "stroid/IO/mesh.h"
+
+namespace py = pybind11;
+
+void register_io_bindings(pybind11::module_ &m) {
+    py::enum_<stroid::IO::VISUALIZATION_MODE>(m, "VISUALIZATION_MODE")
+        .value("NONE", stroid::IO::VISUALIZATION_MODE::NONE)
+        .value("ELEMENT_ID", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID)
+        .value("BOUNDARY_ELEMENT_ID", stroid::IO::VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID)
+        .export_values();
+
+    m.def(
+        "SaveStroidMesh",
+        &stroid::IO::SaveStroidMesh,
+        py::arg("mesh"),
+        py::arg("filename"),
+        py::arg("comment")="",
+        "Save a Stroid mesh to a file."
+    );
+    m.def(
+        "SaveMesh",
+        py::overload_cast<const stroid::StroidMesh&, const std::string&>(&stroid::IO::SaveMesh),
+        py::arg("mesh"),
+        py::arg("filename")
+    );
+    m.def(
+        "SaveVTU",
+        py::overload_cast<const stroid::StroidMesh&, const std::string&>(&stroid::IO::SaveVTU),
+        py::arg("mesh"),
+        py::arg("filename")
+    );
+    m.def(
+        "ViewMesh",
+        py::overload_cast<const stroid::StroidMesh&, const std::string&, stroid::IO::VISUALIZATION_MODE, const std::string&, int>(&stroid::IO::ViewMesh),
+        py::arg("mesh"),
+        py::arg("title")="",
+        py::arg("mode")=stroid::IO::VISUALIZATION_MODE::ELEMENT_ID,
+        py::arg("host")="localhost",
+        py::arg("port")=19916
+    );
+
+    m.def(
+        "VisualizeFaceValence",
+        py::overload_cast<const stroid::StroidMesh&, const std::string&, int>(&stroid::IO::VisualizeFaceValence),
+        py::arg("mesh"),
+        py::arg("host")="localhost",
+        py::arg("port")=19916
+    );
+
+    m.def(
+        "ParseStroidMesh",
+        [](const std::string& buf) {
+            std::stringstream ss;
+            ss << buf;
+            auto r = stroid::IO::ParseStroidMesh(ss);
+            if (!r.has_value()) {
+                throw std::runtime_error("Parsing failed: " + r.error());
+            }
+            return std::move(r.value());
+        }
+    );
+
+    m.def(
+        "LoadStroidMesh",
+        [](const std::string& filename) {
+            auto r = stroid::IO::LoadStroidMesh(filename);
+            if (!r.has_value()) {
+                throw std::runtime_error("Loading " + filename + " failed: " + r.error());
+            }
+            return std::move(r.value());
+        }
+    );
+}

src/python/IO/bindings.h

@@ -0,0 +1,5 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+
+void register_io_bindings(pybind11::module_& m);

src/python/bindings.cpp

@@ -0,0 +1,36 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include "config/bindings.h"
+#include "exceptions/bindings.h"
+#include "IO/bindings.h"
+#include "refinement/bindings.h"
+#include "utils/bindings.h"
+
+#include "stroid/exceptions/stroid_error.h"
+
+#include "stroid/stroid.h"
+#include "stroid/version.h"
+
+PYBIND11_MODULE(_stroid, m) {
+    m.doc() = "Python bindings for stroid library.";
+
+    register_utils_bindings(m);
+
+    auto exceptionsMod = m.def_submodule("exceptions", "Exceptions Bindings");
+    register_exceptions_bindings(exceptionsMod);
+
+    auto configMod = m.def_submodule("config", "Config Bindings");
+    register_config_bindings(configMod);
+
+    auto IOMod = m.def_submodule("IO", "IO Bindings");
+    register_io_bindings(IOMod);
+
+    auto refinementMod = m.def_submodule("refinement", "Refinement Bindings");
+    register_refinement_bindings(refinementMod);
+
+    m.def("GenerateMesh", pybind11::overload_cast<const stroid::config::MeshConfig&>(&stroid::GenerateMesh), "Generate a mesh from a MeshConfig object.");
+    m.def("GenerateMesh", pybind11::overload_cast<const std::string&>(&stroid::GenerateMesh), "Generate a mesh from a config file path.");
+
+
+}

src/python/config/bindings.cpp

@@ -0,0 +1,218 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include "bindings.h"
+
+#include "stroid/config/config.h"
+
+namespace py = pybind11;
+
+void register_config_bindings(pybind11::module_& m) {
+    py::class_<stroid::config::OptimizationMethods>(m, "OptimizationMethods")
+        .def(py::init([](bool tmop, bool smoothstep) {
+            return stroid::config::OptimizationMethods{tmop, smoothstep};
+        }), py::arg("tmop") = false, py::arg("smoothstep") = true)
+    .def_property("tmop",
+        [](const stroid::config::OptimizationMethods& self) {
+            return self.tmop;
+        },
+        [](stroid::config::OptimizationMethods& self, bool value) {
+            self.tmop = value;
+        }
+    )
+    .def_property("smoothstep",
+        [](const stroid::config::OptimizationMethods& self) {
+            return self.smoothstep;
+        },
+        [](stroid::config::OptimizationMethods& self, bool value) {
+            self.smoothstep = value;
+        }
+    );
+
+    py::class_<stroid::config::MeshConfig>(m, "MeshConfig")
+        .def(py::init([](py::kwargs kwargs) {
+            int ref_level = 4, order = 3;
+            size_t continuity_order = 2, surface_bdr_id = 1, inf_bdr_id = 2, core_id = 1, envelope_id = 2, vacuum_id=3;
+            bool include_external_domain = true;
+            double r_core = 0.25, r_star = 1.0, flattening = 0.0, r_inf = 6.0, r_instability = 1e-14, core_steepness = 1.0;
+            stroid::config::OptimizationMethods opt_method{.tmop = false, .smoothstep = true};
+
+            return stroid::config::MeshConfig{
+                .refinement_levels = kwargs.contains("refinement_levels") ? kwargs["refinement_levels"].cast<int>() : ref_level,
+                .order = kwargs.contains("order") ? kwargs["order"].cast<int>() : order,
+                .include_external_domain = kwargs.contains("include_external_domain") ? kwargs["include_external_domain"].cast<bool>() : include_external_domain,
+                .r_core = kwargs.contains("r_core") ? kwargs["r_core"].cast<double>() : r_core,
+                .r_star = kwargs.contains("r_star") ? kwargs["r_star"].cast<double>() : r_star,
+                .flattening = kwargs.contains("flattening") ? kwargs["flattening"].cast<double>() : flattening,
+                .r_infinity = kwargs.contains("r_infinity") ? kwargs["r_infinity"].cast<double>() : r_inf,
+                .r_instability = kwargs.contains("r_instability") ? kwargs["r_instability"].cast<double>() : r_instability,
+                .core_steepness = kwargs.contains("core_steepness") ? kwargs["core_steepness"].cast<double>() : core_steepness,
+                .continuity_order = kwargs.contains("continuity_order") ? kwargs["continuity_order"].cast<size_t>() : continuity_order,
+                .surface_bdr_id = kwargs.contains("surface_bdr_id") ? kwargs["surface_bdr_id"].cast<size_t>() : surface_bdr_id,
+                .inf_bdr_id = kwargs.contains("inf_bdr_id") ? kwargs["inf_bdr_id"].cast<size_t>() : inf_bdr_id,
+                .core_id = kwargs.contains("core_id") ? kwargs["core_id"].cast<size_t>() : core_id,
+                .envelope_id = kwargs.contains("envelope_id") ? kwargs["envelope_id"].cast<size_t>() : envelope_id,
+                .vacuum_id = kwargs.contains("vacuum_id") ? kwargs["vacuum_id"].cast<size_t>() : vacuum_id,
+                .optimization_methods = kwargs.contains("optimization_methods") ? kwargs["optimization_methods"].cast<stroid::config::OptimizationMethods>() : opt_method,
+                .core_mapping = kwargs.contains("core_mapping") ? kwargs["core_mapping"].cast<std::string>() : "spherified"
+            };
+        }))
+        .def_property(
+            "refinement_levels",
+            [](const stroid::config::MeshConfig& self) {
+                return self.refinement_levels;
+            },
+            [](stroid::config::MeshConfig& self, int value) {
+                self.refinement_levels = value;
+            }
+        )
+        .def_property(
+            "order",
+            [](const stroid::config::MeshConfig& self) {
+                return self.order;
+            },
+            [](stroid::config::MeshConfig& self, int value) {
+                self.order = value;
+            }
+        )
+        .def_property(
+            "include_external_domain",
+            [](const stroid::config::MeshConfig& self) {
+                return self.include_external_domain;
+            },
+            [](stroid::config::MeshConfig& self, bool value) {
+                self.include_external_domain = value;
+            }
+        )
+        .def_property(
+            "r_core",
+            [](const stroid::config::MeshConfig& self) {
+                return self.r_core;
+            },
+            [](stroid::config::MeshConfig& self, int value) {
+                self.order = value;
+            }
+        )
+        .def_property(
+            "r_star",
+            [](const stroid::config::MeshConfig& self) {
+                return self.r_star;
+            },
+            [](stroid::config::MeshConfig& self, double value) {
+                self.r_star = value;
+            }
+        )
+        .def_property(
+            "flattening",
+            [](const stroid::config::MeshConfig& self) {
+                return self.flattening;
+            },
+            [](stroid::config::MeshConfig& self, double value) {
+                self.flattening = value;
+            }
+        )
+        .def_property(
+            "r_infinity",
+            [](const stroid::config::MeshConfig& self) {
+                return self.r_infinity;
+            },
+            [](stroid::config::MeshConfig& self, double value) {
+                self.r_infinity = value;
+            }
+        )
+        .def_property(
+            "r_instability",
+            [](const stroid::config::MeshConfig& self) {
+                return self.r_instability;
+            },
+            [](stroid::config::MeshConfig& self, double value) {
+                self.r_instability = value;
+            }
+        )
+        .def_property(
+            "core_steepness",
+            [](const stroid::config::MeshConfig& self) {
+                return self.core_steepness;
+            },
+            [](stroid::config::MeshConfig& self, double value) {
+                self.core_steepness = value;
+            }
+        )
+        .def_property(
+            "continuity_order",
+            [](const stroid::config::MeshConfig& self) {
+                return self.continuity_order;
+            },
+            [](stroid::config::MeshConfig& self, size_t value) {
+                self.continuity_order = value;
+            }
+        )
+        .def_property(
+            "surface_bdr_id",
+            [](const stroid::config::MeshConfig& self) {
+                return self.surface_bdr_id;
+            },
+            [](stroid::config::MeshConfig& self, size_t value) {
+                self.surface_bdr_id = value;
+            }
+        )
+        .def_property(
+            "inf_bdr_id",
+            [](const stroid::config::MeshConfig& self) {
+                return self.inf_bdr_id;
+            },
+            [](stroid::config::MeshConfig& self, size_t value) {
+                self.inf_bdr_id = value;
+            }
+        )
+        .def_property(
+            "core_id",
+            [](const stroid::config::MeshConfig& self) {
+                return self.core_id;
+            },
+            [](stroid::config::MeshConfig& self, size_t value) {
+                self.core_id = value;
+            }
+        )
+        .def_property(
+            "envelope_id",
+            [](const stroid::config::MeshConfig& self) {
+                return self.envelope_id;
+            },
+            [](stroid::config::MeshConfig& self, size_t value) {
+                self.envelope_id = value;
+            }
+        )
+        .def_property(
+            "vacuum_id",
+            [](const stroid::config::MeshConfig& self) {
+                return self.vacuum_id;
+            },
+            [](stroid::config::MeshConfig& self, size_t value) {
+                self.vacuum_id = value;
+            }
+        )
+        .def_property(
+            "optimization_methods",
+            [](const stroid::config::MeshConfig& self) {
+                return self.optimization_methods;
+            },
+            [](stroid::config::MeshConfig& self, stroid::config::OptimizationMethods value) {
+                self.optimization_methods = value;
+            }
+        )
+        .def_property(
+            "core_mapping",
+            [](const stroid::config::MeshConfig& self) {
+                return self.core_mapping;
+            },
+            [](stroid::config::MeshConfig& self, const std::string& value) {
+                if (value != "spherified" && value != "multi_block") {
+                    throw std::invalid_argument("Invalid core_mapping value. Must be 'spherified' or 'multi_block'.");
+                }
+                self.core_mapping = value;
+            }
+        )
+        .def("__repr__", [](const stroid::config::MeshConfig& self) {
+            return stroid::config::to_string(self);
+        });
+}

src/python/config/bindings.h

@@ -0,0 +1,5 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+
+void register_config_bindings(pybind11::module_& m);

src/python/exceptions/bindings.cpp

@@ -0,0 +1,14 @@
+#include <pybind11/pybind11.h>
+
+#include "bindings.h"
+#include "stroid/exceptions/exceptions.h"
+
+namespace py = pybind11;
+
+
+void register_exceptions_bindings(py::module_& m) {
+    py::register_exception<stroid::exceptions::StroidError>(m, "StroidError");
+    py::register_exception<stroid::exceptions::StroidContinuityError>(m, "StroidContinuityError", m.attr("StroidError"));
+    py::register_exception<stroid::exceptions::StroidMeshError>(m, "StroidMeshError", m.attr("StroidError"));
+    py::register_exception<stroid::exceptions::StroidMissingReferenceMesh>(m, "StroidMissingReferenceMesh", m.attr("StroidMeshError"));
+}

src/python/exceptions/bindings.h

@@ -0,0 +1,5 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+
+void register_exceptions_bindings(pybind11::module_& m);

src/python/refinement/bindings.cpp

@@ -0,0 +1,11 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include "bindings.h"
+
+#include "stroid/refinement/uniform.h"
+
+namespace py = pybind11;
+
+void register_refinement_bindings(pybind11::module_ &m) {
+    m.def("UniformRefinement", &stroid::refinement::UniformRefinement, py::arg("mesh"), py::arg("levels"), "Perform uniform refinement without breaking the higher order structure");
+}

src/python/refinement/bindings.h

@@ -0,0 +1,5 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+
+void register_refinement_bindings(pybind11::module_& m);

src/python/stroid/__init__.py

@@ -0,0 +1,45 @@
+import io
+import sys
+
+from ._stroid import *
+
+from ._stroid import config
+from ._stroid import exceptions
+from ._stroid import IO
+from ._stroid import refinement
+from ._stroid import stats
+from ._stroid import GenerateMesh
+from ._stroid import StroidMesh
+
+sys.modules['stroid.config'] = config
+sys.modules['stroid.exceptions'] = exceptions
+sys.modules['stroid.IO'] = IO
+sys.modules["stroid.refinement"] = refinement
+sys.modules["stroid.stats"] = stats
+
+
+__all__ = ['config', 'exceptions', 'IO', 'refinement', 'stats', 'GenerateMesh', 'StroidMesh']
+
+import importlib.metadata
+
+
+try:
+    _meta = importlib.metadata.metadata('stroid')
+    __version__ = _meta['Version']
+    __license__ = _meta['License']
+    __description__ = _meta['Summary']
+    __author__ = 'Emily M. Boudreaux'
+    __url__ = 'https://github.com/4D-STAR/stroid'
+except importlib.metadata.PackageNotFoundError :
+    __version__ = 'unknown - Package not installed'
+    __license__ = 'GNU General Public License v3.0'
+    __email__ = 'emily.boudreaux@dartmouth.edu'
+    __url__ = 'https://github.com/4D-STAR/stroid'
+
+
+import os
+from pathlib import Path
+from typing import List
+
+_PACKAGE_DIR = Path(__file__).resolve().parent
+

src/python/utils/bindings.cpp

@@ -0,0 +1,187 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include "bindings.h"
+
+#include "stroid/utils/types.h"
+#include "stroid/utils/mesh_stats.h"
+#include "stroid/utils/mesh_utils.h"
+
+namespace py = pybind11;
+
+void register_stats_bindings(pybind11::module_ &m) {
+    auto statsMod = m.def_submodule("stats", "Stats Bindings");
+    py::enum_<stroid::stats::MeshStatFeatures>(statsMod, "MeshStatFeatures", py::arithmetic())
+     .value("NONE", stroid::stats::MeshStatFeatures::NONE)
+     .value("RADIUS", stroid::stats::MeshStatFeatures::RADIUS)
+     .value("AXES", stroid::stats::MeshStatFeatures::AXES)
+     .value("ELLIPTICITY", stroid::stats::MeshStatFeatures::ELLIPTICITY)
+     .value("BOWING", stroid::stats::MeshStatFeatures::BOWING)
+     .value("CONFORMITY", stroid::stats::MeshStatFeatures::CONFORMITY)
+     .value("JACOBIAN", stroid::stats::MeshStatFeatures::JACOBIAN)
+     .value("VOLUME_AREA", stroid::stats::MeshStatFeatures::VOLUME_AREA)
+     .value("ELEMENT_COUNT", stroid::stats::MeshStatFeatures::ELEMENT_COUNT)
+     .value("MESH_SIZE", stroid::stats::MeshStatFeatures::MESH_SIZE)
+     .value("OUTER_BOUNDS", stroid::stats::MeshStatFeatures::OUTER_BOUNDS)
+     .value("CENTROID", stroid::stats::MeshStatFeatures::CENTROID)
+     .value("CONFIG_META", stroid::stats::MeshStatFeatures::CONFIG_META)
+     .value("BOUNDING_BOX", stroid::stats::MeshStatFeatures::BOUNDING_BOX)
+     .export_values();
+
+    py::class_<stroid::stats::RadiusStats>(statsMod, "RadiusStats")
+    .def_readonly("min", &stroid::stats::RadiusStats::min)
+    .def_readonly("max", &stroid::stats::RadiusStats::max)
+    .def_readonly("mean", &stroid::stats::RadiusStats::mean)
+    .def_readonly("stddev", &stroid::stats::RadiusStats::stddev)
+    .def_readonly("n_samples", &stroid::stats::RadiusStats::n_samples);
+
+    py::class_<stroid::stats::AxisStats>(statsMod, "AxisStats")
+    .def_readonly("semi_major", &stroid::stats::AxisStats::semi_major)
+    .def_readonly("semi_minor", &stroid::stats::AxisStats::semi_minor);
+
+    py::class_<stroid::stats::EllipticityStats>(statsMod, "EllipticityStats")
+    .def_readonly("flattening", &stroid::stats::EllipticityStats::flattening)
+    .def_readonly("polar_equatorial", &stroid::stats::EllipticityStats::polar_equatorial)
+    .def_readonly("radius_uniformity", &stroid::stats::EllipticityStats::radius_uniformity);
+
+    py::class_<stroid::stats::BowingStats>(statsMod, "BowingStats")
+    .def_readonly("max_inward", &stroid::stats::BowingStats::max_inward)
+    .def_readonly("max_outward", &stroid::stats::BowingStats::max_outward)
+    .def_readonly("rms", &stroid::stats::BowingStats::rms)
+    .def_readonly("worst_at_radius", &stroid::stats::BowingStats::worst_at_radius);
+
+    py::class_<stroid::stats::ConformityStats>(statsMod, "ConformityStats")
+    .def_readonly("conforming", &stroid::stats::ConformityStats::conforming)
+    .def_readonly("n_nonconforming_faces", &stroid::stats::ConformityStats::n_nonconforming_faces);
+
+    py::class_<stroid::stats::JacobianStats>(statsMod, "JacobianStats")
+    .def_readonly("detJ_min", &stroid::stats::JacobianStats::detJ_min)
+    .def_readonly("detJ_max", &stroid::stats::JacobianStats::detJ_max)
+    .def_readonly("n_flipped", &stroid::stats::JacobianStats::n_flipped)
+    .def_readonly("min_detJ_ratio", &stroid::stats::JacobianStats::min_detJ_ratio)
+    .def_readonly("worst_ratio_at_radius", &stroid::stats::JacobianStats::worst_ratio_at_radius)
+    .def_readonly("detJ_min_at_radius", &stroid::stats::JacobianStats::detJ_min_at_radius)
+    .def_readonly("n_elements", &stroid::stats::JacobianStats::n_elements);
+
+    py::class_<stroid::stats::VolumeAreaStats>(statsMod, "VolumeAreaStats")
+    .def_readonly("stellar_volume", &stroid::stats::VolumeAreaStats::stellar_volume)
+    .def_readonly("surface_area", &stroid::stats::VolumeAreaStats::surface_area)
+    .def_readonly("analytic_volume", &stroid::stats::VolumeAreaStats::analytic_volume)
+    .def_readonly("analytic_area", &stroid::stats::VolumeAreaStats::analytic_area);
+
+    py::class_<stroid::stats::ElementCounts>(statsMod, "ElementCounts")
+    .def_readonly("total", &stroid::stats::ElementCounts::total)
+    .def_readonly("core", &stroid::stats::ElementCounts::core)
+    .def_readonly("envelope", &stroid::stats::ElementCounts::envelope)
+    .def_readonly("vacuum", &stroid::stats::ElementCounts::vacuum)
+    .def_readonly("other", &stroid::stats::ElementCounts::other)
+    .def_readonly("n_vertices", &stroid::stats::ElementCounts::n_vertices);
+
+    py::class_<stroid::stats::MeshSizeStats>(statsMod, "MeshSizeStats")
+    .def_readonly("h_min", &stroid::stats::MeshSizeStats::h_min)
+    .def_readonly("h_max", &stroid::stats::MeshSizeStats::h_max)
+    .def_readonly("h_mean", &stroid::stats::MeshSizeStats::h_mean)
+    .def_readonly("h_stddev", &stroid::stats::MeshSizeStats::h_stddev);
+
+    py::class_<stroid::stats::OuterBoundsStats>(statsMod, "OuterBoundsStats")
+    .def_readonly("min", &stroid::stats::OuterBoundsStats::min)
+    .def_readonly("max", &stroid::stats::OuterBoundsStats::max)
+    .def_readonly("mean", &stroid::stats::OuterBoundsStats::mean)
+    .def_readonly("n_samples", &stroid::stats::OuterBoundsStats::n_samples);
+
+    py::class_<stroid::stats::CentroidStats>(statsMod, "CentroidStats")
+    .def_readonly("x", &stroid::stats::CentroidStats::x)
+    .def_readonly("y", &stroid::stats::CentroidStats::y)
+    .def_readonly("z", &stroid::stats::CentroidStats::z)
+    .def_readonly("offset", &stroid::stats::CentroidStats::offset);
+
+    py::class_<stroid::stats::ConfigMeta>(statsMod, "ConfigMeta")
+    .def_readonly("r_core", &stroid::stats::ConfigMeta::r_core)
+    .def_readonly("r_star", &stroid::stats::ConfigMeta::r_star)
+    .def_readonly("flattening", &stroid::stats::ConfigMeta::flattening)
+    .def_readonly("r_infinity", &stroid::stats::ConfigMeta::r_infinity)
+    .def_readonly("geom_order", &stroid::stats::ConfigMeta::geom_order)
+    .def_readonly("refinement_levels", &stroid::stats::ConfigMeta::refinement_levels)
+    .def_readonly("has_external_domain", &stroid::stats::ConfigMeta::has_external_domain);
+
+    py::class_<stroid::stats::BoundingBox>(statsMod, "BoundingBox")
+    .def_readonly("xMin", &stroid::stats::BoundingBox::xMin)
+    .def_readonly("xMax", &stroid::stats::BoundingBox::xMax)
+    .def_readonly("yMin", &stroid::stats::BoundingBox::yMin)
+    .def_readonly("yMax", &stroid::stats::BoundingBox::yMax)
+    .def_readonly("zMin", &stroid::stats::BoundingBox::zMin)
+    .def_readonly("zMax", &stroid::stats::BoundingBox::zMax)
+    .def_readonly("valid", &stroid::stats::BoundingBox::valid)
+    .def("dx", &stroid::stats::BoundingBox::dx)
+    .def("dy", &stroid::stats::BoundingBox::dy)
+    .def("dz", &stroid::stats::BoundingBox::dz)
+    .def("diag", &stroid::stats::BoundingBox::diag);
+
+    py::class_<stroid::stats::BoundingBoxStats>(statsMod, "BoundingBoxStats")
+    .def_readonly("core", &stroid::stats::BoundingBoxStats::core)
+    .def_readonly("star", &stroid::stats::BoundingBoxStats::star)
+    .def_readonly("vacuum", &stroid::stats::BoundingBoxStats::vacuum);
+
+    py::class_<stroid::stats::MeshStats>(statsMod, "MeshStats")
+    .def_readonly("computed", &stroid::stats::MeshStats::computed)
+    .def_readonly("radius", &stroid::stats::MeshStats::radius)
+    .def_readonly("axes", &stroid::stats::MeshStats::axes)
+    .def_readonly("ellipticity", &stroid::stats::MeshStats::ellipticity)
+    .def_readonly("bowing", &stroid::stats::MeshStats::bowing)
+    .def_readonly("conformity", &stroid::stats::MeshStats::conformity)
+    .def_readonly("jacobian", &stroid::stats::MeshStats::jacobian)
+    .def_readonly("jacobian_stellar", &stroid::stats::MeshStats::jacobian_stellar)
+    .def_readonly("jacobian_vacuum", &stroid::stats::MeshStats::jacobian_vacuum)
+    .def_readonly("volume", &stroid::stats::MeshStats::volume)
+    .def_readonly("element_counts", &stroid::stats::MeshStats::element_counts)
+    .def_readonly("mesh_size", &stroid::stats::MeshStats::mesh_size)
+    .def_readonly("outer_bounds", &stroid::stats::MeshStats::outer_bounds)
+    .def_readonly("centroid", &stroid::stats::MeshStats::centroid)
+    .def_readonly("config_meta", &stroid::stats::MeshStats::config_meta)
+    .def_readonly("bounding_box", &stroid::stats::MeshStats::bounding_box)
+    .def_readonly("warnings", &stroid::stats::MeshStats::warnings)
+    .def_readonly("errors", &stroid::stats::MeshStats::errors)
+    .def("__repr__", [](const stroid::stats::MeshStats& self) {
+       return stroid::stats::to_string(self);
+    });
+
+    statsMod.attr("MESH_STAT_DEFAULT") = stroid::stats::MESH_STAT_DEFAULT;
+    statsMod.attr("MESH_STAT_ALL") = stroid::stats::MESH_STAT_ALL;
+
+    statsMod.def(
+        "ComputeMeshStats",
+        &stroid::stats::ComputeMeshStats,
+        py::arg("mesh"),
+        py::arg("features") = stroid::stats::MESH_STAT_DEFAULT,
+        py::arg("sample_order")=-1
+    );
+}
+
+void register_type_bindings(py::module_ &m) {
+    py::enum_<stroid::MFEM_MESH_TYPE>(m, "MFEM_MESH_TYPE")
+    .value("SERIAL", stroid::MFEM_MESH_TYPE::SERIAL)
+    .value("PARALLEL", stroid::MFEM_MESH_TYPE::PARALLEL)
+    .export_values();
+
+    py::class_<stroid::StroidMesh>(m, "StroidMesh")
+    .def_property_readonly("type", [](const stroid::StroidMesh& self) {
+        return (self.type == stroid::MFEM_MESH_TYPE::SERIAL) ? "SERIAL" : "PARALLEL";
+    })
+    .def_readonly("config", &stroid::StroidMesh::config)
+    .def_readonly("refinement_levels", &stroid::StroidMesh::refinement_levels)
+    .def("has_mesh", [](const stroid::StroidMesh& self) {
+        return self.mesh != nullptr;
+    })
+    .def("has_rmesh", [](const stroid::StroidMesh& self) {
+        return self.reference_mesh != nullptr;
+    })
+    .def("mesh_stats", &stroid::StroidMesh::mesh_stats)
+    .def("__repr__", [](const stroid::StroidMesh& self) {
+       return std::format("<StroidMesh [{}]: NE: {}, NV: {}>", (self.type == stroid::MFEM_MESH_TYPE::SERIAL) ? "SERIAL" : "PARALLEL", self.mesh->GetNE(), self.mesh->GetNV());
+    });
+}
+
+void register_utils_bindings(pybind11::module_ &m) {
+    register_type_bindings(m);
+    register_stats_bindings(m);
+}
+

src/python/utils/bindings.h

@@ -0,0 +1,5 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+
+void register_utils_bindings(pybind11::module_& m);

subprojects/libconfig.wrap

@@ -1,4 +1,4 @@
 [wrap-git]
 url = https://github.com/4D-STAR/libconfig.git
-revision = v2.0.5
+revision = v2.2.2
 depth = 1

subprojects/packagefiles/pybind11/LICENSE.build

@@ -0,0 +1,19 @@
+Copyright (c) 2021 The Meson development team
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

subprojects/packagefiles/pybind11/meson.build

@@ -0,0 +1,8 @@
+project('pybind11', 'cpp',
+        version : 'v3.0.0',
+        license : 'BSD-3-Clause')
+
+pybind11_incdir = include_directories('include')
+
+pybind11_dep = declare_dependency(
+    include_directories : pybind11_incdir)

subprojects/pybind11.wrap

@@ -0,0 +1,8 @@
+[wrap-git]
+url = https://github.com/pybind/pybind11.git
+revision = v3.0.0
+depth = 1
+patch_directory = pybind11
+
+[provide]
+pybind11 = pybind11_dep

tests/sandbox/sandbox_test.cpp

@@ -1,14 +1,12 @@
-#include "fourdst/config/config.h"
-#include "stroid/config/config.h"
-#include "stroid/IO/mesh.h"
-#include "stroid/topology/curvilinear.h"
-#include "stroid/topology/mapping.h"
-#include "stroid/topology/topology.h"
+#include "stroid/stroid.h"
 
 #include "mfem.hpp"
 
 #include <print>
 
+#include "stroid/topology/optimize.h"
+#include "stroid/utils/mesh_utils.h"
+
 struct SandboxConfig {
     std::string host = "localhost";
     int port = 19916;
@@ -22,16 +20,26 @@ int main() {
     MeshConfig mesh_cfg;
     mesh_cfg.load("default.toml");
 
-    UserConfig user_cfg;
+    // const UserConfig user_cfg;
+    //
+    //
+    // std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(mesh_cfg);
+    // stroid::topology::Finalize(*mesh, mesh_cfg);
+    // stroid::topology::PromoteToHighOrder(*mesh, mesh_cfg);
+    // stroid::topology::ProjectMesh(*mesh, mesh_cfg);
+    //
+    // if (mesh_cfg->optimization_methods.has_value() && mesh_cfg->optimization_methods.value().tmop.has_value() && mesh_cfg->optimization_methods.value().tmop.value()) {
+    //     stroid::topology::ApplyTMOP(*mesh, mesh_cfg);
+    // }
+    //
+    // stroid::IO::ViewMesh(*mesh, "Sandbox Mesh", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, user_cfg->host, user_cfg->port);
+    // stroid::IO::SaveMesh(*mesh, "sandbox.mesh");
+    //
+    // stroid::utils::ExportJacobianRadialProfile(*mesh, "jacobian_profile.csv");
 
 
-    std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(mesh_cfg);
-    stroid::topology::Finalize(*mesh, mesh_cfg);
-    stroid::topology::PromoteToHighOrder(*mesh, mesh_cfg);
-    stroid::topology::ProjectMesh(*mesh, mesh_cfg);
-    stroid::IO::ViewMesh(*mesh, "Sandbox Mesh", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, user_cfg->host, user_cfg->port);
+    stroid::StroidMesh mesh = stroid::GenerateMesh(mesh_cfg);
+    stroid::IO::SaveStroidMesh(mesh, "sandbox.mesh");
 
     return 0;
-
-
 }

tools/geometry_quality_experiment.cpp

@@ -0,0 +1,348 @@
+#include "stroid/stroid.h"
+#include "CLI/CLI.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <map>
+#include <memory>
+#include <string>
+#include <vector>
+
+namespace {
+    constexpr double infinity = std::numeric_limits<double>::infinity();
+
+    struct SampleLocation {
+        int element = -1;
+        mfem::IntegrationPoint point;
+        std::string source = "none";
+    };
+
+    struct ConditioningStats {
+        int elements = 0;
+        size_t samples = 0;
+        size_t nonpositive_samples = 0;
+        double min_det = infinity;
+        double min_sigma = infinity;
+        double max_condition = 0.0;
+        double min_scaled_jacobian = infinity;
+        double contraction_boundary = infinity;
+        SampleLocation det_location;
+        SampleLocation condition_location;
+        SampleLocation contraction_location;
+    };
+
+    class ScopedOutputRedirect {
+        std::streambuf* original;
+    public:
+        ScopedOutputRedirect() : original(std::cout.rdbuf(std::cerr.rdbuf())) {}
+        ~ScopedOutputRedirect() { std::cout.rdbuf(original); }
+    };
+
+    std::vector<mfem::IntegrationPoint> ClosedSamples(int grid_points) {
+        std::vector<mfem::IntegrationPoint> points;
+        for (int i = 0; i < grid_points; ++i) {
+            for (int j = 0; j < grid_points; ++j) {
+                for (int k = 0; k < grid_points; ++k) {
+                    mfem::IntegrationPoint point;
+                    point.Set3(static_cast<double>(i) / (grid_points - 1),
+                               static_cast<double>(j) / (grid_points - 1),
+                               static_cast<double>(k) / (grid_points - 1));
+                    points.push_back(point);
+                }
+            }
+        }
+        for (const double offset : {0.005, 0.010885670927, 0.02}) {
+            for (int corner = 0; corner < 8; ++corner) {
+                mfem::IntegrationPoint point;
+                point.Set3((corner & 1) ? 1.0 - offset : offset,
+                           (corner & 2) ? 1.0 - offset : offset,
+                           (corner & 4) ? 1.0 - offset : offset);
+                points.push_back(point);
+            }
+        }
+        return points;
+    }
+
+    double ColumnNorm(const mfem::DenseMatrix& matrix, int column) {
+        double norm_squared = 0.0;
+        for (int row = 0; row < 3; ++row) {
+            norm_squared += matrix(row, column) * matrix(row, column);
+        }
+        return std::sqrt(norm_squared);
+    }
+
+    double ContractionBoundary(const mfem::DenseMatrix& jacobian,
+                               const mfem::DenseMatrix& direction) {
+        std::array<double, 4> coefficients{};
+        mfem::DenseMatrix mixed(3);
+        for (int mask = 0; mask < 8; ++mask) {
+            int degree = 0;
+            for (int column = 0; column < 3; ++column) {
+                const bool use_direction = (mask & (1 << column)) != 0;
+                degree += use_direction;
+                for (int row = 0; row < 3; ++row) {
+                    mixed(row, column) = use_direction ? direction(row, column) : jacobian(row, column);
+                }
+            }
+            coefficients[degree] += mixed.Det();
+        }
+        if (!(coefficients[0] > 0.0)) return 0.0;
+
+        auto polynomial = [&](double x) {
+            return ((coefficients[3] * x + coefficients[2]) * x + coefficients[1]) * x + coefficients[0];
+        };
+        std::vector<double> breaks{0.0, 1.0};
+        auto add_break = [&](double x) {
+            if (std::isfinite(x) && x > 0.0 && x < 1.0) breaks.push_back(x);
+        };
+        const double a = 3.0 * coefficients[3];
+        const double b = 2.0 * coefficients[2];
+        const double c = coefficients[1];
+        if (a == 0.0) {
+            if (b != 0.0) add_break(-c / b);
+        } else {
+            const double discriminant = b * b - 4.0 * a * c;
+            if (discriminant >= 0.0) {
+                const double q = -0.5 * (b + std::copysign(std::sqrt(discriminant), b));
+                if (q == 0.0) {
+                    add_break(-b / (2.0 * a));
+                } else {
+                    add_break(q / a);
+                    add_break(c / q);
+                }
+            }
+        }
+        std::sort(breaks.begin(), breaks.end());
+        const double scale = std::abs(coefficients[0]) + std::abs(coefficients[1])
+                           + std::abs(coefficients[2]) + std::abs(coefficients[3]);
+        const double tolerance = 64.0 * std::numeric_limits<double>::epsilon() * scale;
+        for (size_t i = 1; i < breaks.size(); ++i) {
+            double right = breaks[i];
+            const double value = polynomial(right);
+            if (value > tolerance) continue;
+            if (std::abs(value) <= tolerance) return right;
+            double left = breaks[i - 1];
+            for (int iteration = 0; iteration < 64; ++iteration) {
+                const double middle = 0.5 * (left + right);
+                if (polynomial(middle) > 0.0) left = middle;
+                else right = middle;
+            }
+            return right;
+        }
+        return infinity;
+    }
+
+    std::unique_ptr<mfem::GridFunction> BuildContractionProbe(
+        stroid::StroidMesh& mesh, mfem::FiniteElementSpace& space, double stellar_radius
+    ) {
+        auto values = std::make_unique<mfem::GridFunction>(&space);
+        *values = 0.0;
+        std::vector<bool> processed(space.GetNDofs(), false);
+        mfem::Array<int> dofs;
+        mfem::Vector physical(3), logical(3);
+        for (int element = 0; element < mesh.mesh->GetNE(); ++element) {
+            const auto& nodes = space.GetFE(element)->GetNodes();
+            space.GetElementDofs(element, dofs);
+            auto* physical_transform = mesh.mesh->GetElementTransformation(element);
+            auto* logical_transform = mesh.reference_mesh->GetElementTransformation(element);
+            for (int local = 0; local < dofs.Size(); ++local) {
+                const int dof = dofs[local] >= 0 ? dofs[local] : -1 - dofs[local];
+                if (processed[dof]) continue;
+                physical_transform->Transform(nodes.IntPoint(local), physical);
+                logical_transform->Transform(nodes.IntPoint(local), logical);
+                const double radius = physical.Norml2();
+                const double logical_radius = std::max({std::abs(logical(0)), std::abs(logical(1)), std::abs(logical(2))});
+                const double fraction = std::min(logical_radius / stellar_radius, 1.0);
+                for (int component = 0; component < 3; ++component) {
+                    (*values)(space.DofToVDof(dof, component)) = radius > 0.0
+                        ? -stellar_radius * fraction * fraction * physical(component) / radius : 0.0;
+                }
+                processed[dof] = true;
+            }
+        }
+        return values;
+    }
+
+    std::map<int, ConditioningStats> InspectMesh(stroid::StroidMesh& mesh, int order,
+        int grid_points, bool contraction_probe, int probe_order) {
+        std::map<int, ConditioningStats> result;
+        const auto closed_samples = ClosedSamples(grid_points);
+        mfem::H1_FECollection probe_collection(probe_order, 3);
+        mfem::FiniteElementSpace probe_space(mesh.mesh.get(), &probe_collection, 3);
+        std::unique_ptr<mfem::GridFunction> probe;
+        if (contraction_probe) probe = BuildContractionProbe(mesh, probe_space, mesh.config.r_star.value());
+        mfem::DenseMatrix probe_values, probe_shape, direction(3);
+        mfem::Array<int> probe_dofs;
+
+        for (int element = 0; element < mesh.mesh->GetNE(); ++element) {
+            const int attribute = mesh.mesh->GetAttribute(element);
+            ++result[attribute].elements;
+            ++result[0].elements;
+            auto* transform = mesh.mesh->GetElementTransformation(element);
+            const bool inspect_probe = contraction_probe && attribute != static_cast<int>(mesh.config.vacuum_id.value());
+            if (inspect_probe) {
+                probe_space.GetElementDofs(element, probe_dofs);
+                probe_values.SetSize(3, probe_dofs.Size());
+                probe_shape.SetSize(probe_dofs.Size(), 3);
+                for (int local = 0; local < probe_dofs.Size(); ++local) {
+                    const int dof = probe_dofs[local] >= 0 ? probe_dofs[local] : -1 - probe_dofs[local];
+                    for (int component = 0; component < 3; ++component) {
+                        probe_values(component, local) = (*probe)(probe_space.DofToVDof(dof, component));
+                    }
+                }
+            }
+            auto inspect_point = [&](const mfem::IntegrationPoint& point, const std::string& source) {
+                transform->SetIntPoint(&point);
+                const mfem::DenseMatrix& jacobian = transform->Jacobian();
+                const double determinant = jacobian.Det();
+                const double sigma_min = jacobian.CalcSingularvalue(2);
+                const double sigma_max = jacobian.CalcSingularvalue(0);
+                const double condition = sigma_min > 0.0 ? sigma_max / sigma_min : infinity;
+                const double denominator = ColumnNorm(jacobian, 0) * ColumnNorm(jacobian, 1) * ColumnNorm(jacobian, 2);
+                const double scaled_jacobian = denominator > 0.0 ? determinant / denominator : 0.0;
+                double boundary = infinity;
+                if (inspect_probe) {
+                    probe_space.GetFE(element)->CalcDShape(point, probe_shape);
+                    mfem::Mult(probe_values, probe_shape, direction);
+                    boundary = ContractionBoundary(jacobian, direction);
+                }
+                for (const int region : {0, attribute}) {
+                    auto& stats = result[region];
+                    ++stats.samples;
+                    if (!(determinant > 0.0)) ++stats.nonpositive_samples;
+                    if (determinant < stats.min_det) {
+                        stats.min_det = determinant;
+                        stats.det_location = {element, point, source};
+                    }
+                    stats.min_sigma = std::min(stats.min_sigma, sigma_min);
+                    stats.min_scaled_jacobian = std::min(stats.min_scaled_jacobian, scaled_jacobian);
+                    if (condition > stats.max_condition) {
+                        stats.max_condition = condition;
+                        stats.condition_location = {element, point, source};
+                    }
+                    if (boundary < stats.contraction_boundary) {
+                        stats.contraction_boundary = boundary;
+                        stats.contraction_location = {element, point, source};
+                    }
+                }
+            };
+            const auto& quadrature = mfem::IntRules.Get(transform->GetGeometryType(), 2 * order + 4);
+            for (int point = 0; point < quadrature.GetNPoints(); ++point) {
+                inspect_point(quadrature.IntPoint(point), "quadrature");
+            }
+            for (const auto& point : closed_samples) inspect_point(point, "closed_grid_and_corner_probes");
+        }
+        return result;
+    }
+
+    void WriteLocation(std::ostream& output, const SampleLocation& location, stroid::StroidMesh& mesh) {
+        output << ',' << location.element << ',' << location.source;
+        if (location.element < 0) {
+            output << ",nan,nan,nan,nan,nan,nan,nan,nan,nan";
+            return;
+        }
+        mfem::Vector physical(3), logical(3);
+        mesh.mesh->GetElementTransformation(location.element)->Transform(location.point, physical);
+        mesh.reference_mesh->GetElementTransformation(location.element)->Transform(location.point, logical);
+        output << ',' << location.point.x << ',' << location.point.y << ',' << location.point.z;
+        for (int component = 0; component < 3; ++component) output << ',' << physical(component);
+        for (int component = 0; component < 3; ++component) output << ',' << logical(component);
+    }
+}
+
+int main(int argc, char** argv) {
+    std::vector<int> orders{1, 2, 3, 4, 5, 6};
+    std::vector<int> refinements{0, 1, 2};
+    std::vector<std::string> mappings{"spherified", "multi_block"};
+    std::string output_path;
+    int grid_points = 5;
+    int probe_order = 3;
+    double core_radius = 0.25;
+    double infinity_radius = 5.0;
+    double flattening = 0.0;
+    bool no_external = false;
+    bool contraction_probe = false;
+    CLI::App app{"Compare signed Jacobians and conditioning of the actual high-order STROID mesh; TMOP is disabled."};
+    app.add_option("--orders", orders, "Geometry orders, comma separated")->delimiter(',')->check(CLI::Range(1, 8));
+    app.add_option("--refinements", refinements, "Uniform refinement levels, comma separated")->delimiter(',')->check(CLI::Range(0, 3));
+    app.add_option("--mappings", mappings, "Core mappings, comma separated")->delimiter(',')->check(CLI::IsMember({"spherified", "multi_block"}));
+    app.add_option("--grid-points", grid_points, "Closed tensor grid points per coordinate, plus near-corner probes")->check(CLI::Range(2, 15));
+    app.add_option("--core-radius", core_radius, "Core radius; stellar radius is one");
+    app.add_option("--infinity-radius", infinity_radius, "Outer reference radius");
+    app.add_option("--flattening", flattening, "Spheroidal flattening");
+    app.add_option("--output", output_path, "New CSV output file; defaults to stdout");
+    app.add_flag("--no-external", no_external, "Omit exterior domain");
+    app.add_flag("--contraction-probe", contraction_probe, "Inspect an interpolated unit logical-radius-squared radial contraction in stellar elements");
+    app.add_option("--probe-order", probe_order, "H1 displacement order for the optional contraction probe")->check(CLI::Range(1, 8));
+    try {
+        app.parse(argc, argv);
+    } catch (const CLI::ParseError& error) {
+        return app.exit(error);
+    }
+
+    try {
+        if (!std::isfinite(core_radius) || core_radius <= 0.0 || core_radius >= 1.0
+            || !std::isfinite(infinity_radius) || infinity_radius <= 1.0
+            || !std::isfinite(flattening) || flattening < 0.0 || flattening >= 1.0) {
+            throw std::invalid_argument("Require 0 < core-radius < 1 < infinity-radius and 0 <= flattening < 1.");
+        }
+        std::ofstream file;
+        if (!output_path.empty()) {
+            if (std::filesystem::exists(output_path)) throw std::runtime_error("Refusing to overwrite existing output: " + output_path);
+            file.open(output_path);
+            if (!file) throw std::runtime_error("Could not open output: " + output_path);
+        }
+        std::ostream& output = output_path.empty() ? std::cout : file;
+        output << std::setprecision(17);
+        output << "mapping,order,refinement,r_core,r_star,r_infinity,flattening,external,grid_points,quadrature_order,probe_order,attribute,elements,samples,nonpositive_samples,min_signed_det,min_sigma,max_condition,min_scaled_jacobian,contraction_boundary_up_to_one";
+        for (const std::string prefix : {"det", "condition", "contraction"}) {
+            output << ',' << prefix << "_element," << prefix << "_source," << prefix << "_xi," << prefix << "_eta," << prefix << "_zeta," << prefix << "_x," << prefix << "_y," << prefix << "_z," << prefix << "_logical_x," << prefix << "_logical_y," << prefix << "_logical_z";
+        }
+        output << '\n';
+        std::cerr << "Sampling actual FE geometry, not the analytical map. Attribute 0 aggregates all regions.\n"
+                     "Signed determinants and scaled Jacobians retain orientation; inf boundary means no sampled root through alpha=1.\n"
+                     "The optional contraction probe is a diagnostic field, not a Newton correction or a production exterior extension.\n";
+        for (const auto& mapping : mappings) {
+            for (const int order : orders) {
+                for (const int refinement : refinements) {
+                    stroid::config::MeshConfig config;
+                    config.core_mapping = mapping;
+                    config.order = order;
+                    config.refinement_levels = refinement;
+                    config.r_core = core_radius;
+                    config.r_star = 1.0;
+                    config.r_infinity = infinity_radius;
+                    config.flattening = flattening;
+                    config.include_external_domain = !no_external;
+                    config.optimization_methods = stroid::config::OptimizationMethods{false, true};
+                    std::cerr << "Inspecting " << mapping << ", order " << order << ", refinement " << refinement << '\n';
+                    stroid::StroidMesh mesh;
+                    {
+                        ScopedOutputRedirect redirect;
+                        mesh = stroid::GenerateMesh(config);
+                    }
+                    const auto regions = InspectMesh(mesh, order, grid_points, contraction_probe, probe_order);
+                    for (const auto& [attribute, stats] : regions) {
+                        output << mapping << ',' << order << ',' << refinement << ',' << core_radius << ",1," << infinity_radius << ',' << flattening << ',' << !no_external << ',' << grid_points << ',' << 2 * order + 4 << ',' << (contraction_probe ? probe_order : 0) << ',' << attribute << ',' << stats.elements << ',' << stats.samples << ',' << stats.nonpositive_samples << ',' << stats.min_det << ',' << stats.min_sigma << ',' << stats.max_condition << ',' << stats.min_scaled_jacobian << ',' << stats.contraction_boundary;
+                        WriteLocation(output, stats.det_location, mesh);
+                        WriteLocation(output, stats.condition_location, mesh);
+                        WriteLocation(output, stats.contraction_location, mesh);
+                        output << '\n';
+                    }
+                    output.flush();
+                    if (!output) throw std::runtime_error("Failed to write experiment output.");
+                }
+            }
+        }
+    } catch (const std::exception& error) {
+        std::cerr << "Geometry quality experiment failed: " << error.what() << '\n';
+        return 1;
+    }
+    return 0;
+}

tools/meson.build

@@ -1 +1,4 @@
 executable('stroid', 'stroid.cpp', dependencies: [stroid_dep, cli11_dep, magic_enum_dep], install: true)
+
+# Opt-in diagnostic driver; deliberately not part of the installed API/tools.
+executable('geometry_quality_experiment', 'geometry_quality_experiment.cpp', dependencies: [stroid_dep, cli11_dep], build_by_default: false, install: false)

tools/python/visualize_winding.py

@@ -1,46 +0,0 @@
-import matplotlib.pyplot as plt
-
-class Box:
-    def __init__(self, scale, vc="red", ec="blue", offset=0):
-        self.scale = scale
-        self.offset = offset
-        self.verticies = [[scale, -scale, -scale], [scale, scale, -scale], [-scale, -scale, -scale], [-scale, scale, -scale], [scale, -scale, scale], [scale, scale, scale], [-scale, -scale, scale], [-scale, scale, scale]]
-        self.edges = [[0, 1], [0, 2], [0, 4], [1, 3], [1, 5], [3, 7], [3, 2], [2, 6], [4, 5], [4, 6], [5,  7], [7, 6]] 
-        self.edge_color = ec
-        self.vertex_color = vc
-         
-    def plot(self, ax, vlabel=True, elabel=False):
-        for vertID, vert in enumerate(self.verticies):
-            ax.scatter(vert[0], vert[1], vert[2], c=self.vertex_color)
-            if vlabel:
-                ax.text(vert[0], vert[1], vert[2], f"{self.offset + vertID}", fontsize=25)
-        for edge in self.edges:
-            ax.plot([self.verticies[edge[0]][0], self.verticies[edge[1]][0]], [self.verticies[edge[0]][1], self.verticies[edge[1]][1]], [self.verticies[edge[0]][2], self.verticies[edge[1]][2]], color=self.edge_color)
-
-class Wedge:
-    def __init__(self, A, B, ec="green"):
-        self.A = A
-        self.B = B
-        self.edge_color = ec
-    def plot(self, ax):
-        for vA, vB in zip(self.A.verticies, self.B.verticies):
-            ax.plot([vA[0], vB[0]], [vA[1], vB[1]], [vA[2], vB[2]], color=self.edge_color)
-
-def main():
-    core = Box(0.5)
-    envelope = Box(2, offset=8)
-    star = Wedge(core, envelope)
-    infinity = Box(5, offset=16)
-    vacuum = Wedge(envelope, infinity)
-
-    fig, ax = plt.subplots(1, 1, figsize=(10, 10), subplot_kw={"projection": "3d"})
-    core.plot(ax)
-    envelope.plot(ax)
-    star.plot(ax)
-    infinity.plot(ax)
-    vacuum.plot(ax)
-    ax.view_init(30, 30)
-    plt.show()
-
-if __name__ == "__main__":
-    main()

tools/stroid.cpp

@@ -66,8 +66,10 @@ int main(int argc, char** argv) {
 
     auto* generate = app.add_subcommand("generate", "Generate a multi-block mesh");
     auto* info = app.add_subcommand("info", "Access information about stroid");
+    auto* view = app.add_subcommand("view", "Display a mesh with glvis");
 
     std::optional<std::string> config_filename;
+    std::optional<std::string> mesh_file;
     std::string output_filename = "stroid.mesh";
     bool view_mesh = false;
     bool no_save = false;
@@ -81,6 +83,37 @@ int main(int argc, char** argv) {
     generate->add_option("--glvis-port", glvis_port, "GLVis server port")->capture_default_str();
     generate->add_option("-o,--output", output_filename, "Output filename base")->capture_default_str();
 
+    view->add_option("--host", glvis_host, "GLVis server host")->capture_default_str();
+    view->add_option("--port", glvis_port, "GLVis server port")->capture_default_str();
+    view->add_option("-f,--file", mesh_file, "Path to .mesh file")->check(CLI::ExistingFile);
+
+    auto to_lower = [](std::string s) {
+        std::string out;
+        out.reserve(s.size());
+        std::ranges::transform(s, std::back_inserter(out), [](unsigned char c) {
+            return static_cast<char>(std::tolower(c));
+        });
+        return out;
+    };
+
+    std::map<std::string, stroid::IO::VISUALIZATION_MODE> mode_map;
+    for (auto [value, name] : magic_enum::enum_entries<stroid::IO::VISUALIZATION_MODE>()) {
+        mode_map[to_lower(std::string(name))] = value;
+    }
+
+    stroid::IO::VISUALIZATION_MODE selected_mode;
+
+    view->add_option("-v,--vis-mode", selected_mode, "Select Visualization mode")
+        ->transform(CLI::CheckedTransformer(mode_map, CLI::ignore_case))
+        ->default_val(stroid::IO::VISUALIZATION_MODE::ELEMENT_ID);
+    view->add_flag_callback("-l,--list", [&]() {
+        std::println("Available Visualization Modes:");
+        for (const auto &name: mode_map | std::views::keys) {
+            std::println("\t - {}", name);
+        }
+        exit(0);
+    });
+
     for (auto [value, name_view] : magic_enum::enum_entries<MESH_FORMATS>()) {
         std::string name{name_view};
         std::ranges::transform(name, name.begin(), ::tolower);
@@ -108,8 +141,6 @@ int main(int argc, char** argv) {
         });
     }
 
-    // generate->require_subcommand(1);
-
     info->add_flag_callback("-v,--version", []() {
         std::println("Stroid Version {}", stroid::version::toString());
         exit(0);
@@ -127,15 +158,31 @@ int main(int argc, char** argv) {
         return app.exit(e);
     }
 
+    if (*view) {
+        if (!mesh_file.has_value()) {
+            throw std::runtime_error("Mesh file must be specified");
+        }
+        mfem::Mesh mesh(mesh_file.value().c_str());
+        stroid::IO::ViewMesh(mesh,
+                             "Mesh Viewer - Colored by Element ID",
+                             selected_mode,
+                             glvis_host,
+                             glvis_port);
+        exit(0);
+
+    }
+
     if (*generate) {
         if (config_filename.has_value()) {
             cfg.load(config_filename.value());
         }
 
+
         const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
         stroid::topology::Finalize(*mesh, cfg);
         stroid::topology::PromoteToHighOrder(*mesh, cfg);
         stroid::topology::ProjectMesh(*mesh, cfg);
+        stroid::topology::OptimizeMesh(*mesh, cfg);
 
         if (!no_save) {
             const std::string& final_path = output_filename;
@@ -195,7 +242,7 @@ int main(int argc, char** argv) {
                                  glvis_port);
         }
     } else if (!*info) {
-        std::println("Usage: {} [generate|info] --help", argv[0]);
+        std::println("Usage: {} [generate|info|view] --help", argv[0]);
     }
 
     return 0;

utils/build-wheels-linux_aarch64.sh

@@ -0,0 +1,144 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+if [[ $# -lt 1 ]]; then
+  echo "Usage: $0 <git-repo-url> [fourdst-wheels-dir]"
+  echo "  fourdst-wheels-dir: optional local directory of fourdst wheels to"
+  echo "  install from instead of PyPI (for bootstrapping a new stroid wheel)"
+  exit 1
+fi
+
+REPO_URL="$1"
+LOCAL_FOURDST_WHEELS="${2:-}"
+WORK_DIR="$(pwd)"
+WHEEL_DIR="${WORK_DIR}/wheels_linux_aarch64"
+
+echo "➤ Creating wheel output directory at ${WHEEL_DIR}"
+mkdir -p "${WHEEL_DIR}"
+
+TMPDIR="$(mktemp -d)"
+echo "➤ Cloning ${REPO_URL} → ${TMPDIR}/project"
+git clone "${REPO_URL}" "${TMPDIR}/project"
+
+DOCKER_MOUNTS=(-v "${WHEEL_DIR}":/io/wheels -v "${TMPDIR}/project":/io/project)
+if [[ -n "${LOCAL_FOURDST_WHEELS}" ]]; then
+  DOCKER_MOUNTS+=(-v "${LOCAL_FOURDST_WHEELS}":/io/fourdst-wheels)
+fi
+
+for IMAGE in \
+  tboudreaux/manylinux_2_28_aarch64_boost_1_88_0:latest
+do
+  docker run --rm \
+    "${DOCKER_MOUNTS[@]}" \
+    "${IMAGE}" \
+    /bin/bash -uxo pipefail -c '
+      cd /io/project
+
+      PKG="$(sed -n "s/^name *= *\"\(.*\)\"/\1/p" pyproject.toml | head -n1)"
+      PKG="${PKG//-/_}"   # wheel filename normalization
+
+      BOOT_PY=/opt/python/cp312-cp312/bin/python
+      "$BOOT_PY" -m pip install --quiet meson
+      VERSION="$("$BOOT_PY" -c "
+import json, subprocess, sys
+out = subprocess.check_output(
+    [sys.executable, \"-m\", \"mesonbuild.mesonmain\", \"introspect\",
+     \"meson.build\", \"--projectinfo\"])
+print(json.loads(out)[\"version\"])
+" 2>/dev/null || true)"
+      if [ -z "$VERSION" ]; then
+        VERSION="$(grep -oE "version *: *.[0-9][0-9a-zA-Z.+-]*" meson.build | head -n1 | grep -oE "[0-9][0-9a-zA-Z.+-]*" || true)"
+      fi
+      if [ -z "$VERSION" ]; then
+        echo "ERROR: could not determine project version; refusing to guess for skip logic"
+        exit 1
+      fi
+      echo "➤ Building ${PKG} ${VERSION}"
+
+      FOURDST_PIN="$(grep -oE "fourdst==[0-9][0-9a-zA-Z.]*" pyproject.toml | head -n1 || true)"
+
+      if [ -d /io/fourdst-wheels ]; then
+        export PIP_FIND_LINKS=/io/fourdst-wheels
+      fi
+
+      build_one() {
+        set -e
+        local PY="$1" PYTAG="$2"
+
+        "$PY" -m pip install --upgrade pip setuptools wheel meson meson-python
+
+        local BUILD_WHEEL_DIR
+        BUILD_WHEEL_DIR="$(mktemp -d)"
+        CC=clang CXX=clang++ "$PY" -m pip wheel . --no-deps \
+          -w "$BUILD_WHEEL_DIR" -vv
+
+        local CURRENT_WHEEL
+        CURRENT_WHEEL="$(find "$BUILD_WHEEL_DIR" -name "*.whl" | head -n1)"
+
+        if [ -n "$FOURDST_PIN" ]; then
+          "$PY" -m pip install --force-reinstall "$FOURDST_PIN"
+          local FOURDST_LIB_PATH
+          FOURDST_LIB_PATH="$("$PY" -c "import fourdst, os; print(os.pathsep.join(fourdst.get_lib_dirs()))")"
+          LD_LIBRARY_PATH="$FOURDST_LIB_PATH" auditwheel repair \
+            --exclude "libcomposition.so*" \
+            --exclude "liblogging.so*" \
+            --exclude "libconst.so*" \
+            --exclude "libreflect_cpp.so*" \
+            -w /io/wheels "$CURRENT_WHEEL"
+
+          local REPAIRED
+          REPAIRED="$(find /io/wheels -name "${PKG}-${VERSION}-${PYTAG}-*manylinux*.whl" | head -n1)"
+          if [ -z "$REPAIRED" ]; then
+            echo "ERROR: repaired wheel for ${PYTAG} not found after auditwheel"
+            return 1
+          fi
+          if unzip -l "$REPAIRED" | grep -E "libcomposition|liblogging|libconst[^a-z]|libreflect_cpp"; then
+            echo "ERROR: repaired wheel contains vendored fourdst libraries"
+            rm -f "$REPAIRED"
+            return 1
+          fi
+        else
+          auditwheel repair -w /io/wheels "$CURRENT_WHEEL"
+        fi
+
+        rm -rf "$BUILD_WHEEL_DIR"
+      }
+
+      FAILED_TAGS=""
+      SKIPPED_TAGS=""
+      BUILT_TAGS=""
+
+      for PY in /opt/python/*/bin/python; do
+        PYTAG="$(basename "$(dirname "$(dirname "$PY")")")"
+
+        if compgen -G "/io/wheels/${PKG}-${VERSION}-${PYTAG}-*manylinux*.whl" > /dev/null; then
+          echo "➤ ${PYTAG}: wheel for ${PKG} ${VERSION} already present — skipping"
+          SKIPPED_TAGS="${SKIPPED_TAGS} ${PYTAG}"
+          continue
+        fi
+
+        echo "================================================================"
+        echo "➤ ${PYTAG}: building ${PKG} ${VERSION}"
+        echo "================================================================"
+        if ( build_one "$PY" "$PYTAG" ); then
+          BUILT_TAGS="${BUILT_TAGS} ${PYTAG}"
+        else
+          echo "✗ ${PYTAG}: BUILD FAILED — continuing with remaining versions"
+          FAILED_TAGS="${FAILED_TAGS} ${PYTAG}"
+        fi
+      done
+
+      echo "================================================================"
+      echo "Summary for ${PKG} ${VERSION}:"
+      echo "  built:  ${BUILT_TAGS:- none}"
+      echo "  skipped:${SKIPPED_TAGS:- none}"
+      echo "  failed: ${FAILED_TAGS:- none}"
+      echo "================================================================"
+
+      if [ -n "$FAILED_TAGS" ]; then
+        echo "✗ Some builds failed:${FAILED_TAGS}"
+        exit 1
+      fi
+      echo "Linux wheels ready in /io/wheels"
+    '
+done

utils/build-wheels-linux_x86_64.sh

@@ -0,0 +1,148 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+if [[ $# -lt 1 ]]; then
+  echo "Usage: $0 <git-repo-url> [fourdst-wheels-dir]"
+  echo "  fourdst-wheels-dir: optional local directory of fourdst wheels to"
+  echo "  install from instead of PyPI (for bootstrapping a new stroid wheel)"
+  exit 1
+fi
+
+REPO_URL="$1"
+LOCAL_FOURDST_WHEELS="${2:-}"
+WORK_DIR="$(pwd)"
+WHEEL_DIR="${WORK_DIR}/wheels_linux_x86_64"
+
+echo "➤ Creating wheel output directory at ${WHEEL_DIR}"
+mkdir -p "${WHEEL_DIR}"
+
+TMPDIR="$(mktemp -d)"
+echo "➤ Cloning ${REPO_URL} → ${TMPDIR}/project"
+git clone "${REPO_URL}" "${TMPDIR}/project"
+
+DOCKER_MOUNTS=(-v "${WHEEL_DIR}":/io/wheels -v "${TMPDIR}/project":/io/project)
+if [[ -n "${LOCAL_FOURDST_WHEELS}" ]]; then
+  DOCKER_MOUNTS+=(-v "${LOCAL_FOURDST_WHEELS}":/io/fourdst-wheels)
+fi
+
+for IMAGE in \
+  tboudreaux/manylinux_2_28_x86_64_boost_1_88_0:latest
+do
+  docker run --rm \
+    "${DOCKER_MOUNTS[@]}" \
+    "${IMAGE}" \
+    /bin/bash -uxo pipefail -c '
+      cd /io/project
+
+      PKG="$(sed -n "s/^name *= *\"\(.*\)\"/\1/p" pyproject.toml | head -n1)"
+      PKG="${PKG//-/_}"   # wheel filename normalization
+
+      BOOT_PY=/opt/python/cp312-cp312/bin/python
+      "$BOOT_PY" -m pip install --quiet meson
+      VERSION="$("$BOOT_PY" -c "
+import json, subprocess, sys
+out = subprocess.check_output(
+    [sys.executable, \"-m\", \"mesonbuild.mesonmain\", \"introspect\",
+     \"meson.build\", \"--projectinfo\"])
+print(json.loads(out)[\"version\"])
+" 2>/dev/null || true)"
+      if [ -z "$VERSION" ]; then
+        # fallback: literal version in project()
+        VERSION="$(grep -oE "version *: *.[0-9][0-9a-zA-Z.+-]*" meson.build | head -n1 | grep -oE "[0-9][0-9a-zA-Z.+-]*" || true)"
+      fi
+      if [ -z "$VERSION" ]; then
+        echo "ERROR: could not determine project version; refusing to guess for skip logic"
+        exit 1
+      fi
+      echo "➤ Building ${PKG} ${VERSION}"
+
+      FOURDST_PIN="$(grep -oE "fourdst==[0-9][0-9a-zA-Z.]*" pyproject.toml | head -n1 || true)"
+
+      if [ -d /io/fourdst-wheels ]; then
+        export PIP_FIND_LINKS=/io/fourdst-wheels
+      fi
+
+      build_one() {
+        set -e
+        local PY="$1" PYTAG="$2"
+
+        "$PY" -m pip install --upgrade pip setuptools wheel meson meson-python
+
+        # Build into a per-iteration temp dir so we repair exactly the
+        # wheel we just built.
+        local BUILD_WHEEL_DIR
+        BUILD_WHEEL_DIR="$(mktemp -d)"
+        CC=clang CXX=clang++ "$PY" -m pip wheel . --no-deps \
+          -w "$BUILD_WHEEL_DIR" -vv
+
+        local CURRENT_WHEEL
+        CURRENT_WHEEL="$(find "$BUILD_WHEEL_DIR" -name "*.whl" | head -n1)"
+
+        if [ -n "$FOURDST_PIN" ]; then
+          "$PY" -m pip install --force-reinstall "$FOURDST_PIN"
+          local FOURDST_LIB_PATH
+          FOURDST_LIB_PATH="$("$PY" -c "import fourdst, os; print(os.pathsep.join(fourdst.get_lib_dirs()))")"
+          LD_LIBRARY_PATH="$FOURDST_LIB_PATH" auditwheel repair \
+            --exclude "libcomposition.so*" \
+            --exclude "liblogging.so*" \
+            --exclude "libconst.so*" \
+            --exclude "libreflect_cpp.so*" \
+            -w /io/wheels "$CURRENT_WHEEL"
+
+          # Post-repair sanity check on the wheel we just produced
+          local REPAIRED
+          REPAIRED="$(find /io/wheels -name "${PKG}-${VERSION}-${PYTAG}-*manylinux*.whl" | head -n1)"
+          if [ -z "$REPAIRED" ]; then
+            echo "ERROR: repaired wheel for ${PYTAG} not found after auditwheel"
+            return 1
+          fi
+          if unzip -l "$REPAIRED" | grep -E "libcomposition|liblogging|libconst[^a-z]|libreflect_cpp"; then
+            echo "ERROR: repaired wheel contains vendored fourdst libraries"
+            rm -f "$REPAIRED"
+            return 1
+          fi
+        else
+          auditwheel repair -w /io/wheels "$CURRENT_WHEEL"
+        fi
+
+        rm -rf "$BUILD_WHEEL_DIR"
+      }
+
+      FAILED_TAGS=""
+      SKIPPED_TAGS=""
+      BUILT_TAGS=""
+
+      for PY in /opt/python/*/bin/python; do
+        PYTAG="$(basename "$(dirname "$(dirname "$PY")")")"
+
+        if compgen -G "/io/wheels/${PKG}-${VERSION}-${PYTAG}-*manylinux*.whl" > /dev/null; then
+          echo "➤ ${PYTAG}: wheel for ${PKG} ${VERSION} already present — skipping"
+          SKIPPED_TAGS="${SKIPPED_TAGS} ${PYTAG}"
+          continue
+        fi
+
+        echo "================================================================"
+        echo "➤ ${PYTAG}: building ${PKG} ${VERSION}"
+        echo "================================================================"
+        if ( build_one "$PY" "$PYTAG" ); then
+          BUILT_TAGS="${BUILT_TAGS} ${PYTAG}"
+        else
+          echo "✗ ${PYTAG}: BUILD FAILED — continuing with remaining versions"
+          FAILED_TAGS="${FAILED_TAGS} ${PYTAG}"
+        fi
+      done
+
+      echo "================================================================"
+      echo "Summary for ${PKG} ${VERSION}:"
+      echo "  built:  ${BUILT_TAGS:- none}"
+      echo "  skipped:${SKIPPED_TAGS:- none}"
+      echo "  failed: ${FAILED_TAGS:- none}"
+      echo "================================================================"
+
+      if [ -n "$FAILED_TAGS" ]; then
+        echo "✗ Some builds failed:${FAILED_TAGS}"
+        exit 1
+      fi
+      echo "Linux wheels ready in /io/wheels"
+    '
+done

utils/build-wheels-macos_aarch64.sh

@@ -0,0 +1,108 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+if [[ $(uname -m) != "arm64" ]]; then
+  echo "Error: This script is intended to run on an Apple Silicon (arm64) Mac."
+  exit 1
+fi
+
+if [[ $# -lt 1 ]]; then
+  echo "Usage: $0 <git-repo-url> [fourdst-wheels-dir]"
+  echo "  fourdst-wheels-dir: optional local directory of fourdst wheels to"
+  echo "  install from instead of PyPI (for bootstrapping a new stroid wheel)"
+  exit 1
+fi
+
+for TOOL in pyenv cmake git; do
+  if ! command -v "$TOOL" &> /dev/null; then
+    echo "Error: ${TOOL} not found."
+    exit 1
+  fi
+done
+
+REPO_URL="$1"
+LOCAL_FOURDST_WHEELS="${2:-}"
+WORK_DIR="$(pwd)"
+WHEEL_DIR="${WORK_DIR}/wheels_macos_aarch64_tmp"
+FINAL_WHEEL_DIR="${WORK_DIR}/wheels_macos_aarch64"
+
+echo "➤ Creating wheel output directories"
+mkdir -p "${WHEEL_DIR}"
+mkdir -p "${FINAL_WHEEL_DIR}"
+
+export MACOSX_DEPLOYMENT_TARGET=15.0
+
+
+TMPDIR="$(mktemp -d)"
+echo "➤ Cloning ${REPO_URL} → ${TMPDIR}/project"
+git clone --depth 1 "${REPO_URL}" "${TMPDIR}/project"
+cd "${TMPDIR}/project"
+
+FOURDST_PIN="$(grep -oE 'fourdst==[0-9][0-9a-zA-Z.]*' pyproject.toml | head -n1 || true)"
+if [[ -n "${FOURDST_PIN}" ]]; then
+  echo "➤ Project depends on ${FOURDST_PIN}; wheel repair will exclude fourdst libraries"
+fi
+
+PYTHON_VERSIONS=("3.9.23" "3.10.18" "3.11.13" "3.12.11" "3.13.5" "3.14.0rc1" "3.14.0rc1t")
+
+eval "$(pyenv init -)"
+
+for PY_VERSION in "${PYTHON_VERSIONS[@]}"; do
+  (
+    set -e
+
+    pyenv shell "${PY_VERSION}"
+    PY="$(pyenv which python)"
+
+    echo "----------------------------------------------------------------"
+    echo "➤ Building for $($PY --version) on macOS arm64"
+    echo "----------------------------------------------------------------"
+
+    "$PY" -m pip install --upgrade pip setuptools wheel meson-python delocate
+    "$PY" -m pip install meson==1.9.1
+
+    if [[ -n "${FOURDST_PIN}" ]]; then
+      if [[ -n "${LOCAL_FOURDST_WHEELS}" ]]; then
+        "$PY" -m pip install --force-reinstall \
+          --find-links "${LOCAL_FOURDST_WHEELS}" "${FOURDST_PIN}"
+      else
+        "$PY" -m pip install --force-reinstall "${FOURDST_PIN}"
+      fi
+    fi
+
+    echo "➤ Building wheel with ccache enabled"
+    echo "➤ Found meson version $(meson --version)"
+
+    CC="ccache clang" CXX="ccache clang++" \
+      "$PY" -m pip wheel . --no-deps --no-build-isolation \
+        -w "${WHEEL_DIR}" -v
+    CURRENT_WHEEL=$(find "${WHEEL_DIR}" -name "*.whl" | head -n 1)
+
+    echo "➤ Repairing wheel with delocate"
+    if [[ -n "${FOURDST_PIN}" ]]; then
+        delocate-wheel --require-archs arm64 \
+          -e composition -e logging -e const -e reflect_cpp \
+          -w "${FINAL_WHEEL_DIR}" -v "$CURRENT_WHEEL"
+    else
+        delocate-wheel --require-archs arm64 \
+          -w "${FINAL_WHEEL_DIR}" -v "$CURRENT_WHEEL"
+    fi
+
+    REPAIRED_WHEEL="${FINAL_WHEEL_DIR}/$(basename "$CURRENT_WHEEL")"
+
+    if [[ -n "${FOURDST_PIN}" ]]; then
+      if unzip -l "${REPAIRED_WHEEL}" | grep -E 'libcomposition|liblogging|libconst|libreflect_cpp'; then
+        echo "ERROR: repaired wheel contains vendored fourdst libraries"
+        exit 1
+      fi
+    fi
+
+    rm "$CURRENT_WHEEL"
+  )
+done
+
+# Cleanup
+rm -rf "${TMPDIR}"
+rm -rf "${WHEEL_DIR}"
+
+echo "All builds complete. Artifacts in ${FINAL_WHEEL_DIR}"

utils/installPyEnvVersions.sh

@@ -0,0 +1,14 @@
+#!/bin/bash
+
+pyenv install 3.8.20
+pyenv install 3.9.23
+pyenv install 3.10.18
+pyenv install 3.11.13
+pyenv install 3.12.11
+pyenv install 3.13.5
+pyenv install 3.13.5t
+pyenv install 3.14.0rc1
+pyenv install 3.14.0rc1t
+pyenv install pypy3.10-7.3.19
+pyenv install pypy3.11-7.3.20
+

utils/readme.md

@@ -0,0 +1,28 @@
+# Wheel Generation
+This directory contains scripts to generate precompiled python wheels for GridFire
+
+# Notes
+- MacOS wheels can only be generated on macos
+- aarch64 wheels can only be generated on aarch64 machines
+- x86_64 wheels can only be generated on x86_64 machines
+- linux wheels can be generated on any linux machine, but the target architecture must match the machine architecture
+- Running each script will take **a very long time** (could be upwards of half of a day depending on your system) and will require roughly 2GB of disk space
+- When generating MacOS wheels, you must have all the correct versions of python installed with `pyenv`. Run the script `utils/wheels/installPyEnvVersions.sh` to install the correct versions of python.
+
+# Usage
+Once you know you are on the correct machine, run the script for your desired architecture and operating system. For example, to generate a macos x86_64 wheel, run:
+
+```bash
+./build-wheels-macos-aarch64.sh https://github.com/4D-STAR/GridFire
+```
+
+Once you have all the wheels generated (which will likely require multiple systems), copy all the wheels into a single
+directory (lets assume its called `wheels` and in the root of the directory) and then run (from the root of the repository):
+
+```bash
+python -m pip install --upgrade build
+python -m build --sdist --outdir wheels
+twine upload wheels/*
+```
+
+Thie will also take a while (it needs to upload all the wheels to PyPI) but will result in all the wheels being uploaded to PyPI.