Merge pull request #3 from tboudreaux/feat/multi_block

Vacuum coordinate and multi block

CHANGELOG.md

@@ -0,0 +1,31 @@
+## v0.2.1 (2026-03-20)
+
+### Feat
+
+- **tools**: Added winding visualization
+- **src**: Enabled Vacuume Region Generation
+
+### Fix
+
+- **cuviliniear.cpp**: removed unused BSD proc header
+- **stroid**: exit properly after stroid -v
+
+## v0.1.0 (2026-01-31)
+
+### Feat
+
+- **release**: added macos release script
+- **build-check**: added setup checks for stl compatibility
+- **stroid**: improved CLI & fixed gcc bug
+- **stroid**: added command line and tests
+- **stroid**: first working version
+- **stroid**: initial commit
+
+### Fix
+
+- **linux**: fixed linux compilation
+- **build-check**: removed c lang
+
+### Refactor
+
+- **stroid**: changed default mesh format from VTU to MFEM

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.0
+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,16 +1,22 @@
-cmake = import('cmake')
+mfem_dep = dependency('mfem', required : false)
-mfem_cmake_options = cmake.subproject_options()
-mfem_cmake_options.add_cmake_defines({
-    'MFEM_ENABLE_EXAMPLES': 'OFF',
-    'MFEM_ENABLE_TESTING': 'OFF',
-    'MFEM_ENABLE_MINIAPPS': 'OFF',
-    'MFEM_USE_BENCMARK': 'OFF',
-    'BUILD_SHARED_LIBS': 'OFF',
-    'BUILD_STATIC_LIBS': 'ON',
-})
-mfem_cmake_options.set_install(true)
 
-mfem_sp = cmake.subproject(
+if not mfem_dep.found()
-    'mfem',
+    cmake = import('cmake')
-    options: mfem_cmake_options)
+    mfem_cmake_options = cmake.subproject_options()
-mfem_dep = mfem_sp.dependency('mfem')
+    mfem_cmake_options.add_cmake_defines({
+        'MFEM_ENABLE_EXAMPLES': 'OFF',
+        'MFEM_ENABLE_TESTING': 'OFF',
+        'MFEM_ENABLE_MINIAPPS': 'OFF',
+        'MFEM_USE_BENCMARK': 'OFF',
+        'BUILD_SHARED_LIBS': 'OFF',
+        'BUILD_STATIC_LIBS': 'ON',
+    })
+    mfem_cmake_options.set_install(true)
+
+    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

docs/html/_logo_8png.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/annotated.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/classes.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/config_8h.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/curvilinear_8cpp.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>
@@ -110,12 +110,10 @@ $(function(){initNavTree('curvilinear_8cpp.html',''); initResizable(true); });
 <div class="textblock"><code>#include &quot;<a class="el" href="curvilinear_8h.html">stroid/topology/curvilinear.h</a>&quot;</code><br />
 <code>#include &quot;<a class="el" href="mapping_8h.html">stroid/topology/mapping.h</a>&quot;</code><br />
 <code>#include &lt;iostream&gt;</code><br />
-<code>#include &lt;memory&gt;</code><br />
-<code>#include &lt;sys/proc.h&gt;</code><br />
 </div><div class="textblock"><div class="dynheader">
 Include dependency graph for curvilinear.cpp:</div>
 <div class="dyncontent">
-<div class="center"><div class="zoom"><iframe scrolling="no" frameborder="0" src="curvilinear_8cpp__incl.svg" width="100%" height="438"><p><b>This browser is not able to show SVG: try Firefox, Chrome, Safari, or Opera instead.</b></p></iframe></div></div>
+<div class="center"><iframe scrolling="no" frameborder="0" src="curvilinear_8cpp__incl.svg" width="488" height="184"><p><b>This browser is not able to show SVG: try Firefox, Chrome, Safari, or Opera instead.</b></p></iframe></div>
 </div>
 </div><table class="memberdecls">
 <tr class="heading"><td colspan="2"><h2 class="groupheader"><a id="namespaces" name="namespaces"></a>

docs/html/curvilinear_8cpp__incl.map

@@ -1,15 +1,11 @@
 <map id="src/lib/topology/curvilinear.cpp" name="src/lib/topology/curvilinear.cpp">
-<area shape="rect" id="Node000001" title=" " alt="" coords="348,5,544,31"/>
+<area shape="rect" id="Node000001" title=" " alt="" coords="199,5,395,31"/>
 <area shape="rect" id="Node000002" href="$curvilinear_8h.html" title=" " alt="" coords="207,79,387,105"/>
-<area shape="poly" id="edge1_Node000001_Node000002" title=" " alt="" coords="424,31,337,74,335,69,418,32"/>
+<area shape="poly" id="edge1_Node000001_Node000002" title=" " alt="" coords="299,29,299,64,295,62,295,33"/>
 <area shape="rect" id="Node000006" href="$mapping_8h.html" title=" " alt="" coords="13,79,184,105"/>
-<area shape="poly" id="edge5_Node000001_Node000006" title=" " alt="" coords="386,34,173,78,172,73,388,29"/>
+<area shape="poly" id="edge5_Node000001_Node000006" title=" " alt="" coords="261,33,147,76,145,71,265,30"/>
 <area shape="rect" id="Node000007" title=" " alt="" coords="411,79,482,105"/>
-<area shape="poly" id="edge9_Node000001_Node000007" title=" " alt="" coords="448,29,449,65,444,62,445,33"/>
+<area shape="poly" id="edge9_Node000001_Node000007" title=" " alt="" coords="323,29,410,72,405,72,322,34"/>
-<area shape="rect" id="Node000008" title=" " alt="" coords="506,79,574,105"/>
-<area shape="poly" id="edge10_Node000001_Node000008" title=" " alt="" coords="463,29,515,70,509,69,461,34"/>
-<area shape="rect" id="Node000009" title=" " alt="" coords="598,79,681,105"/>
-<area shape="poly" id="edge11_Node000001_Node000009" title=" " alt="" coords="479,29,595,73,590,73,479,34"/>
 <area shape="rect" id="Node000003" title=" " alt="" coords="5,153,85,178"/>
 <area shape="poly" id="edge2_Node000002_Node000003" title=" " alt="" coords="255,108,100,151,99,146,254,103"/>
 <area shape="rect" id="Node000004" href="$config_8h.html" title=" " alt="" coords="108,153,251,178"/>

docs/html/curvilinear_8cpp__incl.md5

@@ -1 +1 @@
-6950d4e2891a4c4f47ed638d27c5bf68
+31c392d827608579a7e0e3bf6524e7f9

docs/html/curvilinear_8cpp__incl.svg

@@ -4,9 +4,10 @@
 <!-- Generated by graphviz version 14.1.2 (20260124.0452)
  -->
 <!-- Title: src/lib/topology/curvilinear.cpp Pages: 1 -->
-<!--zoomable 138 -->
+<svg width="366pt" height="138pt"
+ viewBox="0.00 0.00 366.00 138.00" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
 
-<svg id="main" version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" xml:space="preserve" onload="init(evt)">
+<svg id="main" version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" xml:space="preserve">
 
 <style type="text/css"><![CDATA[
 .node, .edge {opacity: 0.7;}
@@ -14,53 +15,16 @@
 .edge:hover path { stroke: red; }
 .edge:hover polygon { stroke: red; fill: red; }
 ]]></style>
-
-<defs>
-  <circle id="rim" cx="0" cy="0" r="7"/>
-  <circle id="rim2" cx="0" cy="0" r="3.5"/>
-  <g id="zoomPlus">
-    <use xlink:href="#rim" fill="#404040"><set attributeName="fill" to="#808080" begin="zoomplus.mouseover" end="zoomplus.mouseout"/></use>
-    <path d="M-4,0h8M0,-4v8" fill="none" stroke="white" stroke-width="1.5" pointer-events="none"/>
-  </g>
-  <g id="zoomMin">
-    <use xlink:href="#rim" fill="#404040"><set attributeName="fill" to="#808080" begin="zoomminus.mouseover" end="zoomminus.mouseout"/></use>
-    <path d="M-4,0h8" fill="none" stroke="white" stroke-width="1.5" pointer-events="none"/>
-  </g>
-  <g id="arrowUp" transform="translate(30 24)">
-    <use xlink:href="#rim"/>
-    <path pointer-events="none" fill="none" stroke="white" stroke-width="1.5" d="M0,-3.0v7 M-2.5,-0.5L0,-3.0L2.5,-0.5"/>
-  </g>
-  <g id="arrowRight" transform="rotate(90) translate(36 -43)">
-    <use xlink:href="#rim"/>
-    <path pointer-events="none" fill="none" stroke="white" stroke-width="1.5" d="M0,-3.0v7 M-2.5,-0.5L0,-3.0L2.5,-0.5"/>
-  </g>
-  <g id="arrowDown" transform="rotate(180) translate(-30 -48)">
-    <use xlink:href="#rim"/>
-    <path pointer-events="none" fill="none" stroke="white" stroke-width="1.5" d="M0,-3.0v7 M-2.5,-0.5L0,-3.0L2.5,-0.5"/>
-  </g>
-  <g id="arrowLeft" transform="rotate(270) translate(-36 17)">
-    <use xlink:href="#rim"/>
-    <path pointer-events="none" fill="none" stroke="white" stroke-width="1.5" d="M0,-3.0v7 M-2.5,-0.5L0,-3.0L2.5,-0.5"/>
-  </g>
-  <g id="resetDef">
-    <use xlink:href="#rim2" fill="#404040"><set attributeName="fill" to="#808080" begin="reset.mouseover" end="reset.mouseout"/></use>
-  </g>
-</defs>
-<script type="application/ecmascript">
-var viewWidth = 515;
-var viewHeight = 138;
-var sectionId = 'dynsection-0';
-</script>
 <script type="application/ecmascript" xlink:href="svg.min.js"/>
 <svg id="graph" class="graph">
-<g id="viewport">
+<g id="graph0" class="graph" transform="scale(1 1) rotate(0) translate(4 133.75)">
 <title>src/lib/topology/curvilinear.cpp</title>
 <!-- Node1 -->
 <g id="Node000001" class="node">
 <title>Node1</title>
 <g id="a_Node000001"><a xlink:title=" ">
-<polygon fill="#999999" stroke="#666666" points="404.38,-129.75 257.12,-129.75 257.12,-110.5 404.38,-110.5 404.38,-129.75"/>
+<polygon fill="#999999" stroke="#666666" points="292.38,-129.75 145.12,-129.75 145.12,-110.5 292.38,-110.5 292.38,-129.75"/>
-<text xml:space="preserve" text-anchor="middle" x="330.75" y="-116.25" font-family="Helvetica,sans-Serif" font-size="10.00">src/lib/topology/curvilinear.cpp</text>
+<text xml:space="preserve" text-anchor="middle" x="218.75" y="-116.25" font-family="Helvetica,sans-Serif" font-size="10.00">src/lib/topology/curvilinear.cpp</text>
 </a>
 </g>
 </g>
@@ -77,8 +41,8 @@ var sectionId = 'dynsection-0';
 <g id="edge1_Node000001_Node000002" class="edge">
 <title>Node1&#45;&gt;Node2</title>
 <g id="a_edge1_Node000001_Node000002"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M311.75,-110.09C294.29,-101.79 268.15,-89.36 248.03,-79.79"/>
+<path fill="none" stroke="#63b8ff" d="M218.75,-110.33C218.75,-103.82 218.75,-94.67 218.75,-86.37"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="249.74,-76.74 239.21,-75.6 246.74,-83.06 249.74,-76.74"/>
+<polygon fill="#63b8ff" stroke="#63b8ff" points="222.25,-86.37 218.75,-76.37 215.25,-86.37 222.25,-86.37"/>
 </a>
 </g>
 </g>
@@ -95,8 +59,8 @@ var sectionId = 'dynsection-0';
 <g id="edge5_Node000001_Node000006" class="edge">
 <title>Node1&#45;&gt;Node6</title>
 <g id="a_edge5_Node000001_Node000006"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M286.18,-110.03C241.7,-100.96 173.24,-86.99 125.26,-77.2"/>
+<path fill="none" stroke="#63b8ff" d="M193.47,-110.09C169.39,-101.49 132.89,-88.44 105.77,-78.75"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="126.12,-73.8 115.62,-75.23 124.72,-80.66 126.12,-73.8"/>
+<polygon fill="#63b8ff" stroke="#63b8ff" points="107.17,-75.53 96.58,-75.46 104.81,-82.12 107.17,-75.53"/>
 </a>
 </g>
 </g>
@@ -113,44 +77,8 @@ var sectionId = 'dynsection-0';
 <g id="edge9_Node000001_Node000007" class="edge">
 <title>Node1&#45;&gt;Node7</title>
 <g id="a_edge9_Node000001_Node000007"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M330.75,-110.33C330.75,-103.82 330.75,-94.67 330.75,-86.37"/>
+<path fill="none" stroke="#63b8ff" d="M237.75,-110.09C255.21,-101.79 281.35,-89.36 301.47,-79.79"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="334.25,-86.37 330.75,-76.37 327.25,-86.37 334.25,-86.37"/>
+<polygon fill="#63b8ff" stroke="#63b8ff" points="302.76,-83.06 310.29,-75.6 299.76,-76.74 302.76,-83.06"/>
-</a>
-</g>
-</g>
-<!-- Node8 -->
-<g id="Node000008" class="node">
-<title>Node8</title>
-<g id="a_Node000008"><a xlink:title=" ">
-<polygon fill="#e0e0e0" stroke="#999999" points="426.38,-74.5 375.12,-74.5 375.12,-55.25 426.38,-55.25 426.38,-74.5"/>
-<text xml:space="preserve" text-anchor="middle" x="400.75" y="-61" font-family="Helvetica,sans-Serif" font-size="10.00">memory</text>
-</a>
-</g>
-</g>
-<!-- Node1&#45;&gt;Node8 -->
-<g id="edge10_Node000001_Node000008" class="edge">
-<title>Node1&#45;&gt;Node8</title>
-<g id="a_edge10_Node000001_Node000008"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M342.62,-110.09C352.84,-102.32 367.82,-90.92 380.01,-81.65"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="382.01,-84.53 387.85,-75.69 377.77,-78.96 382.01,-84.53"/>
-</a>
-</g>
-</g>
-<!-- Node9 -->
-<g id="Node000009" class="node">
-<title>Node9</title>
-<g id="a_Node000009"><a xlink:title=" ">
-<polygon fill="#e0e0e0" stroke="#999999" points="506.62,-74.5 444.88,-74.5 444.88,-55.25 506.62,-55.25 506.62,-74.5"/>
-<text xml:space="preserve" text-anchor="middle" x="475.75" y="-61" font-family="Helvetica,sans-Serif" font-size="10.00">sys/proc.h</text>
-</a>
-</g>
-</g>
-<!-- Node1&#45;&gt;Node9 -->
-<g id="edge11_Node000001_Node000009" class="edge">
-<title>Node1&#45;&gt;Node9</title>
-<g id="a_edge11_Node000001_Node000009"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M355.35,-110.09C378.68,-101.52 413.99,-88.56 440.35,-78.88"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="441.43,-82.21 449.61,-75.47 439.02,-75.64 441.43,-82.21"/>
 </a>
 </g>
 </g>
@@ -237,27 +165,6 @@ var sectionId = 'dynsection-0';
 </g>
 </g>
 </svg>
-
-<g id="navigator" transform="translate(0 0)" fill="#404254">
-  <rect fill="#f2f5e9" fill-opacity="0.5" stroke="#606060" stroke-width=".5" x="0" y="0" width="60" height="60"/>
-  <use id="zoomplus" xlink:href="#zoomPlus" x="17" y="9" onmousedown="handleZoom(evt,'in')"/>
-  <use id="zoomminus" xlink:href="#zoomMin" x="42" y="9" onmousedown="handleZoom(evt,'out')"/>
-  <use id="reset" xlink:href="#resetDef" x="30" y="36" onmousedown="handleReset()"/>
-   <use id="arrowup" xlink:href="#arrowUp" x="0" y="0" onmousedown="handlePan(0,-1)"/>
-  <use id="arrowright" xlink:href="#arrowRight" x="0" y="0" onmousedown="handlePan(1,0)"/>
-  <use id="arrowdown" xlink:href="#arrowDown" x="0" y="0" onmousedown="handlePan(0,1)"/>
-  <use id="arrowleft" xlink:href="#arrowLeft" x="0" y="0" onmousedown="handlePan(-1,0)"/>
-</g>
-<svg viewBox="0 0 15 15" width="100%" height="30px" preserveAspectRatio="xMaxYMin meet">
- <g id="arrow_out" transform="scale(0.3 0.3)">
-  <a xlink:href="curvilinear_8cpp__incl_org.svg" target="_base">
-   <rect id="button" ry="5" rx="5" y="6" x="6" height="38" width="38"
-        fill="#f2f5e9" fill-opacity="0.5" stroke="#606060" stroke-width="1.0"/>
-   <path id="arrow"
-     d="M 11.500037,31.436501 C 11.940474,20.09759 22.043105,11.32322 32.158766,21.979434 L 37.068811,17.246167 C 37.068811,17.246167 37.088388,32 37.088388,32 L 22.160133,31.978069 C 22.160133,31.978069 26.997745,27.140456 26.997745,27.140456 C 18.528582,18.264221 13.291696,25.230495 11.500037,31.436501 z"
-     style="fill:#404040;"/>
-  </a>
- </g>
 </svg>
 
 <style type='text/css'>

docs/html/curvilinear_8cpp__incl_org.svg

@@ -4,16 +4,16 @@
 <!-- Generated by graphviz version 14.1.2 (20260124.0452)
  -->
 <!-- Title: src/lib/topology/curvilinear.cpp Pages: 1 -->
-<svg width="515pt" height="138pt"
+<svg width="366pt" height="138pt"
- viewBox="0.00 0.00 515.00 138.00" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
+ viewBox="0.00 0.00 366.00 138.00" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
 <g id="graph0" class="graph" transform="scale(1 1) rotate(0) translate(4 133.75)">
 <title>src/lib/topology/curvilinear.cpp</title>
 <!-- Node1 -->
 <g id="Node000001" class="node">
 <title>Node1</title>
 <g id="a_Node000001"><a xlink:title=" ">
-<polygon fill="#999999" stroke="#666666" points="404.38,-129.75 257.12,-129.75 257.12,-110.5 404.38,-110.5 404.38,-129.75"/>
+<polygon fill="#999999" stroke="#666666" points="292.38,-129.75 145.12,-129.75 145.12,-110.5 292.38,-110.5 292.38,-129.75"/>
-<text xml:space="preserve" text-anchor="middle" x="330.75" y="-116.25" font-family="Helvetica,sans-Serif" font-size="10.00">src/lib/topology/curvilinear.cpp</text>
+<text xml:space="preserve" text-anchor="middle" x="218.75" y="-116.25" font-family="Helvetica,sans-Serif" font-size="10.00">src/lib/topology/curvilinear.cpp</text>
 </a>
 </g>
 </g>
@@ -30,8 +30,8 @@
 <g id="edge1_Node000001_Node000002" class="edge">
 <title>Node1&#45;&gt;Node2</title>
 <g id="a_edge1_Node000001_Node000002"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M311.75,-110.09C294.29,-101.79 268.15,-89.36 248.03,-79.79"/>
+<path fill="none" stroke="#63b8ff" d="M218.75,-110.33C218.75,-103.82 218.75,-94.67 218.75,-86.37"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="249.74,-76.74 239.21,-75.6 246.74,-83.06 249.74,-76.74"/>
+<polygon fill="#63b8ff" stroke="#63b8ff" points="222.25,-86.37 218.75,-76.37 215.25,-86.37 222.25,-86.37"/>
 </a>
 </g>
 </g>
@@ -48,8 +48,8 @@
 <g id="edge5_Node000001_Node000006" class="edge">
 <title>Node1&#45;&gt;Node6</title>
 <g id="a_edge5_Node000001_Node000006"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M286.18,-110.03C241.7,-100.96 173.24,-86.99 125.26,-77.2"/>
+<path fill="none" stroke="#63b8ff" d="M193.47,-110.09C169.39,-101.49 132.89,-88.44 105.77,-78.75"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="126.12,-73.8 115.62,-75.23 124.72,-80.66 126.12,-73.8"/>
+<polygon fill="#63b8ff" stroke="#63b8ff" points="107.17,-75.53 96.58,-75.46 104.81,-82.12 107.17,-75.53"/>
 </a>
 </g>
 </g>
@@ -66,44 +66,8 @@
 <g id="edge9_Node000001_Node000007" class="edge">
 <title>Node1&#45;&gt;Node7</title>
 <g id="a_edge9_Node000001_Node000007"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M330.75,-110.33C330.75,-103.82 330.75,-94.67 330.75,-86.37"/>
+<path fill="none" stroke="#63b8ff" d="M237.75,-110.09C255.21,-101.79 281.35,-89.36 301.47,-79.79"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="334.25,-86.37 330.75,-76.37 327.25,-86.37 334.25,-86.37"/>
+<polygon fill="#63b8ff" stroke="#63b8ff" points="302.76,-83.06 310.29,-75.6 299.76,-76.74 302.76,-83.06"/>
-</a>
-</g>
-</g>
-<!-- Node8 -->
-<g id="Node000008" class="node">
-<title>Node8</title>
-<g id="a_Node000008"><a xlink:title=" ">
-<polygon fill="#e0e0e0" stroke="#999999" points="426.38,-74.5 375.12,-74.5 375.12,-55.25 426.38,-55.25 426.38,-74.5"/>
-<text xml:space="preserve" text-anchor="middle" x="400.75" y="-61" font-family="Helvetica,sans-Serif" font-size="10.00">memory</text>
-</a>
-</g>
-</g>
-<!-- Node1&#45;&gt;Node8 -->
-<g id="edge10_Node000001_Node000008" class="edge">
-<title>Node1&#45;&gt;Node8</title>
-<g id="a_edge10_Node000001_Node000008"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M342.62,-110.09C352.84,-102.32 367.82,-90.92 380.01,-81.65"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="382.01,-84.53 387.85,-75.69 377.77,-78.96 382.01,-84.53"/>
-</a>
-</g>
-</g>
-<!-- Node9 -->
-<g id="Node000009" class="node">
-<title>Node9</title>
-<g id="a_Node000009"><a xlink:title=" ">
-<polygon fill="#e0e0e0" stroke="#999999" points="506.62,-74.5 444.88,-74.5 444.88,-55.25 506.62,-55.25 506.62,-74.5"/>
-<text xml:space="preserve" text-anchor="middle" x="475.75" y="-61" font-family="Helvetica,sans-Serif" font-size="10.00">sys/proc.h</text>
-</a>
-</g>
-</g>
-<!-- Node1&#45;&gt;Node9 -->
-<g id="edge11_Node000001_Node000009" class="edge">
-<title>Node1&#45;&gt;Node9</title>
-<g id="a_edge11_Node000001_Node000009"><a xlink:title=" ">
-<path fill="none" stroke="#63b8ff" d="M355.35,-110.09C378.68,-101.52 413.99,-88.56 440.35,-78.88"/>
-<polygon fill="#63b8ff" stroke="#63b8ff" points="441.43,-82.21 449.61,-75.47 439.02,-75.64 441.43,-82.21"/>
 </a>
 </g>
 </g>

docs/html/curvilinear_8h.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_000004_000002.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_000005_000002.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_000009_000000.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_000010_000002.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_000012_000002.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_31a6c70da1404d38d1b51e38b19548a1.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_3f4f28e99864f6e7d54ad0c93ce0a2a5.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_49e56c817e5e54854c35e136979f97ca.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_68267d1309a1af8e8297ef4c3efbcdba.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_82887d4f338e92537fa9d809bec8c614.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_87d18a4dc5174905bfd7d2dc734defe6.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_936721fd6a792737eaef7358bcd49428.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_9d5018b9ddb63a582aa4c6b91bec2f8c.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_9da86e8abab39b749c6020360d8e4403.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_b0856f6b0d80ccb263b2f415c91f9e17.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_c34d5e8363cf0aa3fabc4f3fad3412a4.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_c85d3e3c5052e9ad9ce18c6863244a25.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/dir_dc36fd7496f87f6989d6d6a9445d600c.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/files.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/functions.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/functions_vars.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/graph_legend.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/index.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>
@@ -226,9 +226,9 @@ C++ Interface</h2>
 <div class="ttc" id="amesh_8h_html"><div class="ttname"><a href="mesh_8h.html">mesh.h</a></div></div>
 <div class="ttc" id="anamespacestroid_1_1_i_o_html_a496f5c16eaffda5922a0b96c1f525dab"><div class="ttname"><a href="namespacestroid_1_1_i_o.html#a496f5c16eaffda5922a0b96c1f525dab">stroid::IO::ViewMesh</a></div><div class="ttdeci">void ViewMesh(mfem::Mesh &amp;mesh, const std::string &amp;title, VISUALIZATION_MODE mode, const std::string &amp;vishost, int visport)</div><div class="ttdoc">Stream a mesh to a running GLVis server for interactive viewing.</div><div class="ttdef"><b>Definition</b> mesh.cpp:25</div></div>
 <div class="ttc" id="anamespacestroid_1_1_i_o_html_ad4048304d8a0c7075d2b2a6e465d0b6eaee96e14c2b71bd59252006289ba464cf"><div class="ttname"><a href="namespacestroid_1_1_i_o.html#ad4048304d8a0c7075d2b2a6e465d0b6eaee96e14c2b71bd59252006289ba464cf">stroid::IO::VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID</a></div><div class="ttdeci">@ BOUNDARY_ELEMENT_ID</div><div class="ttdoc">Color boundary-adjacent elements by boundary attribute/ID.</div><div class="ttdef"><b>Definition</b> mesh.h:15</div></div>
-<div class="ttc" id="anamespacestroid_1_1topology_html_a5907aa2e639cda703d48d177abc37caf"><div class="ttname"><a href="namespacestroid_1_1topology.html#a5907aa2e639cda703d48d177abc37caf">stroid::topology::PromoteToHighOrder</a></div><div class="ttdeci">void PromoteToHighOrder(mfem::Mesh &amp;mesh, const fourdst::config::Config&lt; config::MeshConfig &gt; &amp;config)</div><div class="ttdoc">Promote a mesh to high-order by attaching an H1 nodal finite element space.</div><div class="ttdef"><b>Definition</b> curvilinear.cpp:9</div></div>
+<div class="ttc" id="anamespacestroid_1_1topology_html_a5907aa2e639cda703d48d177abc37caf"><div class="ttname"><a href="namespacestroid_1_1topology.html#a5907aa2e639cda703d48d177abc37caf">stroid::topology::PromoteToHighOrder</a></div><div class="ttdeci">void PromoteToHighOrder(mfem::Mesh &amp;mesh, const fourdst::config::Config&lt; config::MeshConfig &gt; &amp;config)</div><div class="ttdoc">Promote a mesh to high-order by attaching an H1 nodal finite element space.</div><div class="ttdef"><b>Definition</b> curvilinear.cpp:7</div></div>
 <div class="ttc" id="anamespacestroid_1_1topology_html_a62774bcba7ea1a485892dcd4bed6425b"><div class="ttname"><a href="namespacestroid_1_1topology.html#a62774bcba7ea1a485892dcd4bed6425b">stroid::topology::Finalize</a></div><div class="ttdeci">void Finalize(mfem::Mesh &amp;mesh, const fourdst::config::Config&lt; config::MeshConfig &gt; &amp;config)</div><div class="ttdoc">Finalize topology, validate orientation, and apply uniform refinement.</div><div class="ttdef"><b>Definition</b> topology.cpp:90</div></div>
-<div class="ttc" id="anamespacestroid_1_1topology_html_a836ed13e5bac63e7952c3ce4e5532e78"><div class="ttname"><a href="namespacestroid_1_1topology.html#a836ed13e5bac63e7952c3ce4e5532e78">stroid::topology::ProjectMesh</a></div><div class="ttdeci">void ProjectMesh(mfem::Mesh &amp;mesh, const fourdst::config::Config&lt; config::MeshConfig &gt; &amp;config)</div><div class="ttdoc">Project high-order mesh nodes using the configured curvilinear mapping.</div><div class="ttdef"><b>Definition</b> curvilinear.cpp:15</div></div>
+<div class="ttc" id="anamespacestroid_1_1topology_html_a836ed13e5bac63e7952c3ce4e5532e78"><div class="ttname"><a href="namespacestroid_1_1topology.html#a836ed13e5bac63e7952c3ce4e5532e78">stroid::topology::ProjectMesh</a></div><div class="ttdeci">void ProjectMesh(mfem::Mesh &amp;mesh, const fourdst::config::Config&lt; config::MeshConfig &gt; &amp;config)</div><div class="ttdoc">Project high-order mesh nodes using the configured curvilinear mapping.</div><div class="ttdef"><b>Definition</b> curvilinear.cpp:13</div></div>
 <div class="ttc" id="anamespacestroid_1_1topology_html_abc0d8a1fb8e9c5ac0e259e4c93db7892"><div class="ttname"><a href="namespacestroid_1_1topology.html#abc0d8a1fb8e9c5ac0e259e4c93db7892">stroid::topology::BuildSkeleton</a></div><div class="ttdeci">std::unique_ptr&lt; mfem::Mesh &gt; BuildSkeleton(const fourdst::config::Config&lt; config::MeshConfig &gt; &amp;config)</div><div class="ttdoc">Build the initial multi-block mesh topology for the star model.</div><div class="ttdef"><b>Definition</b> topology.cpp:10</div></div>
 <div class="ttc" id="atopology_8h_html"><div class="ttname"><a href="topology_8h.html">topology.h</a></div></div>
 </div><!-- fragment --><h1><a class="anchor" id="autotoc_md8"></a>

docs/html/mainpage_8md.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/mapping_8cpp.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/mapping_8h.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/mesh_8cpp.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/mesh_8h.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/mesh__utils_8cpp.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/mesh__utils_8h.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/namespacemembers.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/namespacemembers_enum.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/namespacemembers_func.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/namespaces.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/namespacestroid.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/namespacestroid_1_1_i_o.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/namespacestroid_1_1config.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/namespacestroid_1_1topology.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/namespacestroid_1_1utils.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/search/all_10.js

@@ -0,0 +1,7 @@
+var searchData=
+[
+  ['vacuum_5fid_0',['vacuum_id',['../structstroid_1_1config_1_1_mesh_config.html#a36c7b05405691393d381086a27ea36ad',1,'stroid::config::MeshConfig']]],
+  ['viewmesh_1',['ViewMesh',['../namespacestroid_1_1_i_o.html#a496f5c16eaffda5922a0b96c1f525dab',1,'stroid::IO']]],
+  ['visualization_5fmode_2',['VISUALIZATION_MODE',['../namespacestroid_1_1_i_o.html#ad4048304d8a0c7075d2b2a6e465d0b6e',1,'stroid::IO']]],
+  ['visualizefacevalence_3',['VisualizeFaceValence',['../namespacestroid_1_1_i_o.html#a8100e130b3a49fdee48bc3c4d4e63963',1,'stroid::IO']]]
+];

docs/html/search/variables_5.js

@@ -0,0 +1,8 @@
+var searchData=
+[
+  ['r_5fcore_0',['r_core',['../structstroid_1_1config_1_1_mesh_config.html#a5c68a895f73dc82a38a8daac22a83ad7',1,'stroid::config::MeshConfig']]],
+  ['r_5finfinity_1',['r_infinity',['../structstroid_1_1config_1_1_mesh_config.html#ac7546899ebbfe191ea3a8bf2403b31eb',1,'stroid::config::MeshConfig']]],
+  ['r_5finstability_2',['r_instability',['../structstroid_1_1config_1_1_mesh_config.html#a4da6d99ff7ba24d2f917e1fd98ddd877',1,'stroid::config::MeshConfig']]],
+  ['r_5fstar_3',['r_star',['../structstroid_1_1config_1_1_mesh_config.html#a3fe80a30990d484dcc39b6f9a0befc05',1,'stroid::config::MeshConfig']]],
+  ['refinement_5flevels_4',['refinement_levels',['../structstroid_1_1config_1_1_mesh_config.html#a8cafcbebf64ae251517118eb152de981',1,'stroid::config::MeshConfig']]]
+];

docs/html/search/variables_6.js

@@ -0,0 +1,4 @@
+var searchData=
+[
+  ['surface_5fbdr_5fid_0',['surface_bdr_id',['../structstroid_1_1config_1_1_mesh_config.html#a33f25ff277aa8834065e04ccb9cdbdda',1,'stroid::config::MeshConfig']]]
+];

docs/html/search/variables_7.js

@@ -0,0 +1,4 @@
+var searchData=
+[
+  ['vacuum_5fid_0',['vacuum_id',['../structstroid_1_1config_1_1_mesh_config.html#a36c7b05405691393d381086a27ea36ad',1,'stroid::config::MeshConfig']]]
+];

docs/html/structstroid_1_1config_1_1_mesh_config-members.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/structstroid_1_1config_1_1_mesh_config.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/topology_8cpp.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

docs/html/topology_8h.html

@@ -31,7 +31,7 @@
  <tr id="projectrow">
   <td id="projectlogo"><img alt="Logo" src="Logo.png"/></td>
   <td id="projectalign">
-   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.0</span>
+   <div id="projectname">stroid<span id="projectnumber">&#160;v0.2.1</span>
    </div>
    <div id="projectbrief">Multi-block curvilinear mesh generation</div>
   </td>

meson.build

@@ -1,4 +1,4 @@
-project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.2.0', 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,30 +92,62 @@ 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 |
+| Parameter                       | Description                                                                                                                                                                                                                                        | Default       |
-|-------------------------|-----------------------------------------------------------------------------------------------------|---------|
+|---------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|---------------|
-| refinement_levels       | Number of uniform refinement levels to apply to the mesh after generation                           | 4       |
+| 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       |
+| 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    |
+| include_external_domain         | Whether to include an external domain extending to r_infinity                                                                                                                                                                                      | true          |
-| r_core                  | The radius of the core region of the star                                                           | 1.5     |
+| r_core                          | The radius of the core region of the star                                                                                                                                                                                                          | 1.5           |
-| r_star                  | The radius of the star                                                                              | 5.0     |
+| r_star                          | The radius of the star                                                                                                                                                                                                                             | 5.0           |
-| flattening              | The flattening factor of the star (0 for spherical, >0 for oblate)                                  | 0       |
+| flattening                      | The flattening factor of the star (0 for spherical, >0 for oblate)                                                                                                                                                                                 | 0             |
-| r_infinity              | The outer radius of the external domain (if included)                                               | 6.0     |
+| r_infinity                      | The outer radius of the external domain (if included)                                                                                                                                                                                              | 6.0           |
-| r_instability           | The radius at which no transformations are applied to the initial topology (to avoid singularities) | 1e-14   |
+| r_instability                   | The radius at which no transformations are applied to the initial topology (to avoid singularities)                                                                                                                                                | 1e-14         |
-| core_steepness          | The steepness of the transition between the core and envelope regions of the star                   | 1.0     |
+| core_steepness                  | The steepness of the transition between the core and envelope regions of the star                                                                                                                                                                  | 1.0           |
-| surface_bdr_id          | The boundary ID to assign to the surface of the star                                                | 1       |
+| surface_bdr_id                  | The boundary ID to assign to the surface of the star                                                                                                                                                                                               | 1             |
-| inf_bdr_id              | The boundary ID to assign to the outer boundary of the external domain (if included)                | 2       |
+| inf_bdr_id                      | The boundary ID to assign to the outer boundary of the external domain (if included)                                                                                                                                                               | 2             |
-| core_id                 | The material ID to assign to the core region of the star                                            | 1       |
+| 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       |
+| 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       |
+| 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',
+    input : 'version.h.in',
-    output : 'stroid.h',
+    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 {
-    /**
+    inline StroidMesh GenerateMesh(const fourdst::config::Config<stroid::config::MeshConfig>& cfg) {
-     * @brief Version helpers for the stroid library.
+        StroidMesh sm;
-     */
+        sm.type = MFEM_MESH_TYPE::SERIAL;
-    struct version {
+        sm.config = *cfg;
-        static constexpr int major = @STROID_VERSION_MAJOR@;
+        auto reference = stroid::topology::BuildSkeleton(cfg);
-        static constexpr int minor = @STROID_VERSION_MINOR@;
+        stroid::topology::Finalize(*reference, cfg);
-        static constexpr int patch = @STROID_VERSION_PATCH@;
+        sm.refinement_levels = cfg->refinement_levels.value_or(0);
-        static constexpr const char* tag = "@STROID_VERSION_TAG@";
 
-        static std::string toString() {
+        sm.reference_mesh = std::move(reference);
-            std::string versionStr = std::to_string(major) + "." +
+        sm.mesh = utils::BuildProjected(*sm.reference_mesh, cfg);
-                                     std::to_string(minor) + "." +
+        if (cfg->optimization_methods.has_value() && cfg->optimization_methods.value().tmop.has_value() && cfg->optimization_methods.value().tmop.value()) {
-                                     std::to_string(patch);
+            stroid::topology::ApplyTMOP(*sm.mesh, cfg);
-            if (std::string(tag) != "") {
-                versionStr += "-" + std::string(tag);
-            }
-            return versionStr;
         }
+        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;
+        };
 
-        friend std::ostream& operator<<(std::ostream& os, const version&) {
+        cfg.mutate(Mutator);
-            os << toString();
+        return GenerateMesh(cfg);
-            return os;
+    }
-        }
+    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";
+        mfem::socketstream sol_sock(vishost.c_str(), visport);
-        int  visport   = 19916;
-        mfem::socketstream sol_sock(vishost, 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

@@ -2,12 +2,45 @@
 #include "stroid/topology/mapping.h"
 
 #include <iostream>
-#include <memory>
+
-#include <sys/proc.h>
+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);
     }
@@ -55,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));
+    std::unique_ptr<ScalarMeshField> BuildExteriorCoordinate(
-        //     }
+        mfem::Mesh& mesh,
-        //
+        mfem::Mesh& reference_mesh,
-        //     TransformPoint(pos, config, 0);
+        const fourdst::config::Config<config::MeshConfig>& config
-        //
+    ) {
-        //     for (int d = 0; d < vDim; ++d) {
+        if (!config->include_external_domain.value_or(true)) return nullptr;
-        //         nodes(fes->DofToVDof(i, d)) = pos(d);
+        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;
     }
 }