feat(mean_field): added dimensions, discritization, and start of eos

The full rewrite of mean_field into something maintainable is progressing. dimensions is mostly done, discritization (domain, blocks, and fields) is done, and eos is progressing quickly

build-check/meson.build

@@ -8,6 +8,13 @@ if not cpp.has_header('memory') or not cpp.compiles('''
   error('MFEM 4.10 requires a C++17-or-newer compiler.')
 endif
 
+if not cpp.links('''
+  #include <print>
+  int main() { std::println("SERiF"); }
+''', args: platform_cpp_args + platform_link_args, name: 'C++23 std::println support')
+  error('SERiF requires C++23 std::println support. On macOS, use -Dmacos_deployment_target=15.0 or newer with a recent Xcode toolchain.')
+endif
+
 cmake_program = find_program('cmake', required: true)
 # With meson-python, the generated native file pins this to the interpreter
 # building the wheel. For command-line builds, PATH controls the selection.

build-python/meson.build

@@ -2,6 +2,7 @@ python_extension = []
 python_extension_dir = ''
 
 if get_option('build_python')
+  error('Python bindings are not currently supported')
   if is_wasm
     error('The native CPython extension cannot be built for Emscripten. Disable build_python for the WASM profile.')
   endif

meson.build

@@ -29,10 +29,7 @@ endif
 subdir('build-check')
 subdir('build-config')
 subdir('include/meson_mfem_template')
-subdir('examples')
+subdir('src')
-subdir('benchmarks')
-subdir('web')
-subdir('build-python')
 subdir('tests')
 
 summary(
@@ -52,8 +49,6 @@ summary(
     'FMS': mfem_features.get('fms'),
     'Algoim': mfem_features.get('algoim'),
     'Python bindings': get_option('build_python'),
-    'examples': get_option('build_examples'),
-    'benchmarks': get_option('build_benchmarks'),
     'tests': get_option('build_tests'),
   },
   section: 'meson-mfem-template ' + meson.project_version(),

meson_options.txt

@@ -1,8 +1,6 @@
 option('allow_preinstalled', type: 'boolean', value: true, description: 'Allow one complete, feature-compatible system MFEM. Set false for the dependency-hermetic source fallback.')
 option('feature_profile', type: 'combo', choices: ['minimal', 'portable', 'full'], value: 'portable', description: 'Defaults for auto features. full is strict and requires all selected non-bundled TPLs in dependency_prefix.')
 option('build_python', type: 'boolean', value: false, description: 'Build the nanobind Python extension and install native assets inside the Python package.')
-option('build_examples', type: 'boolean', value: true, description: 'Build MFEM capability examples.')
-option('build_benchmarks', type: 'boolean', value: false, description: 'Build the native MPI-capable backend benchmark and register the benchmark suite.')
 option('build_tests', type: 'boolean', value: true, description: 'Register serial, parallel, and Python smoke tests when applicable.')
 option('install_mfem', type: 'boolean', value: true, description: 'Install the source-built MFEM bundle and its runtime dependencies.')
 option('dependency_prefix', type: 'string', value: '', description: 'Explicit prefix containing optional TPLs. Used even when allow_preinstalled=false.')
@@ -11,7 +9,7 @@ option('mfem_cuda_arch', type: 'string', value: 'native', description: 'CUDA arc
 option('mfem_hip_arch', type: 'string', value: '', description: 'HIP GPU target, for example gfx90a. Empty lets the toolchain choose.')
 option('mfem_precision', type: 'combo', choices: ['double', 'single'], value: 'double', description: 'Floating-point precision used consistently by MFEM and Hypre.')
 option('python_install_native', type: 'boolean', value: true, description: 'Install MFEM headers, CMake/pkg-config data, and runtime libraries under the Python package.')
-option('macos_deployment_target', type: 'string', value: '11.0', description: 'Minimum macOS version for native artifacts; ignored on other systems. Empty uses the toolchain default.')
+option('macos_deployment_target', type: 'string', value: '15.0', description: 'Minimum macOS version for native artifacts (15.0 supports C++23 std::print); ignored on other systems. Empty uses the toolchain default.')
 
 option('mfem_mpi', type: 'feature', value: 'auto', description: 'MPI parallel MFEM; implies Hypre.')
 option('mfem_metis', type: 'feature', value: 'auto', description: 'METIS graph partitioning.')

src/include/serif/dimensions/dimensions.hpp

@@ -0,0 +1,155 @@
+#pragma once
+
+#include <string_view>
+
+#include "serif/dimensions/quantities.hpp"
+
+namespace serif::dimensions {
+    struct Dimensionless final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "dimensionless";
+    };
+
+    struct Mass final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "mass";
+    };
+
+    struct Length final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "length";
+    };
+
+    struct Time final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "time";
+    };
+
+    struct Area final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "area";
+    };
+
+    struct Volume final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "volume";
+    };
+
+    struct Density final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "density";
+    };
+
+    struct SurfaceDensity final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "surface_density";
+    };
+
+    struct NumberDensity final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "number_density";
+    };
+
+    struct Pressure final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "pressure";
+    };
+
+    struct Temperature final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "temperature";
+    };
+
+    struct Entropy final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "entropy";
+    };
+
+    struct SpecificEntropy final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "specific_entropy";
+    };
+
+    struct ChemicalPotential final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "chemical_potential";
+    };
+
+    struct Energy final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "energy";
+    };
+
+    struct InternalEnergy final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "internal_energy";
+    };
+
+    struct SpecificEnergy final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "specific_energy";
+    };
+
+    struct SpecificInternalEnergy final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "specific_internal_energy";
+    };
+
+    struct SpecificEnthalpy final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "specific_enthalpy";
+    };
+
+    struct EnergyDensity final : ThermodynamicQuantity {
+        static constexpr std::string_view identifier = "energy_density";
+    };
+
+    struct GravitationalPotential final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "gravitational_potential";
+    };
+
+    struct Velocity final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "velocity";
+    };
+
+    struct Acceleration final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "acceleration";
+    };
+
+    struct Frequency final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "frequency";
+    };
+
+    struct AngularVelocity final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "angular_velocity";
+    };
+
+    struct Momentum final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "momentum";
+    };
+
+    struct AngularMomentum final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "angular_momentum";
+    };
+
+    struct MomentOfInertia final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "moment_of_inertia";
+    };
+
+    struct Force final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "force";
+    };
+
+    struct Torque final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "torque";
+    };
+
+    struct Power final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "power";
+    };
+
+    struct Luminosity final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "luminosity";
+    };
+
+    struct MassFlowRate final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "mass_flow_rate";
+    };
+
+    struct Opacity final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "opacity";
+    };
+
+    struct DynamicViscosity final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "dynamic_viscosity";
+    };
+
+    struct KinematicViscosity final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "kinematic_viscosity";
+    };
+
+    struct MagneticFluxDensity final : PhysicalQuantity {
+        static constexpr std::string_view identifier = "magnetic_flux_density";
+    };
+}

src/include/serif/dimensions/partial.hpp

@@ -0,0 +1,83 @@
+#pragma once
+
+#include "serif/dimensions/quantities.hpp"
+#include "serif/utils/misc/concepts/numeric.hpp"
+
+namespace serif::dimensions {
+    template <ThermodynamicQuantityType Quantity>
+    struct WithRespectTo final {};
+
+    template <ThermodynamicQuantityType OutputQuantity, ThermodynamicQuantityType InputQuantity>
+    class PartialDerivative final {
+    public:
+        explicit constexpr PartialDerivative(const double value) noexcept : m_value(value){}
+
+        [[nodiscard]] constexpr double value() const noexcept {
+            return m_value;
+        }
+
+        friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator+(
+            const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
+            const PartialDerivative<OutputQuantity, InputQuantity>& rhs
+        ) noexcept {
+            return PartialDerivative<OutputQuantity, InputQuantity>{lhs.m_value + rhs.m_value};
+        }
+
+        friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator-(
+            const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
+            const PartialDerivative<OutputQuantity, InputQuantity>& rhs
+        ) noexcept {
+            return PartialDerivative<OutputQuantity, InputQuantity>{lhs.m_value - rhs.m_value};
+        }
+
+        template <utils::misc::concepts::IsNumeric Scalar>
+        friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator*(
+            const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
+            const Scalar& rhs
+        ) noexcept {
+            return PartialDerivative<OutputQuantity, InputQuantity>{lhs.m_value * rhs};
+        }
+
+        template <utils::misc::concepts::IsNumeric Scalar>
+        friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator*(
+            const Scalar& lhs,
+            const PartialDerivative<OutputQuantity, InputQuantity>& rhs
+        ) noexcept {
+            return PartialDerivative<OutputQuantity, InputQuantity>{lhs * rhs.m_value};
+        }
+
+        template <utils::misc::concepts::IsNumeric Scalar>
+        friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator/(
+            const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
+            const Scalar& rhs
+        ) noexcept {
+            return PartialDerivative<OutputQuantity, InputQuantity>{lhs.m_value / rhs};
+        }
+
+        template <utils::misc::concepts::IsNumeric Scalar>
+        friend constexpr std::partial_ordering operator<=>(
+            const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
+            Scalar rhs
+        ) noexcept {
+            return lhs.m_value <=> rhs;
+        }
+
+        template <utils::misc::concepts::IsNumeric Scalar>
+        friend constexpr std::partial_ordering operator<=>(
+            Scalar lhs,
+            const PartialDerivative<OutputQuantity, InputQuantity>& rhs
+        ) noexcept {
+            return lhs <=> rhs.m_value;
+        }
+
+        friend constexpr std::partial_ordering operator<=>(
+            const PartialDerivative<OutputQuantity, InputQuantity> &lhs,
+            const PartialDerivative<OutputQuantity, InputQuantity> &rhs
+        ) noexcept {
+            return lhs.m_value <=> rhs.m_value;
+        }
+
+    private:
+        double m_value;
+    };
+}

src/include/serif/dimensions/quantities.hpp

@@ -0,0 +1,138 @@
+#pragma once
+
+#include <type_traits>
+#include <compare>
+#include <concepts>
+
+#include "serif/utils/misc/std_helper/cleaning.hpp"
+#include "serif/utils/misc/concepts/numeric.hpp"
+
+namespace serif::dimensions {
+    using utils::misc::std_helper::NeitherConstNorVolatile;
+    using utils::misc::std_helper::NeitherConstNorVolatileReference;
+
+    using utils::misc::concepts::IsNumeric;
+
+    struct PhysicalQuantity { };
+    struct ThermodynamicQuantity: PhysicalQuantity { };
+
+    template <typename QuantityCandidate>
+    concept PhysicalQuantityType =
+        std::same_as<QuantityCandidate, NeitherConstNorVolatile<QuantityCandidate>> &&
+        std::derived_from<QuantityCandidate, PhysicalQuantity>;
+
+    template <typename QuantityCandidate>
+    concept ThermodynamicQuantityType =
+        PhysicalQuantityType<QuantityCandidate> &&
+        std::derived_from<QuantityCandidate, ThermodynamicQuantity>;
+
+    template <PhysicalQuantityType Quantity>
+    class QuantityValue final {
+    public:
+        explicit constexpr QuantityValue(const double value) noexcept : m_value(value) {}
+
+        [[nodiscard]] constexpr double value() const noexcept {
+            return m_value;
+        }
+
+        [[nodiscard]] friend constexpr bool operator==(
+            const QuantityValue &,
+            const QuantityValue &
+        ) noexcept = default;
+
+        friend constexpr QuantityValue operator+(
+            const QuantityValue &lhs,
+            const QuantityValue &rhs
+        ) noexcept {
+            return QuantityValue{lhs.m_value + rhs.m_value};
+        }
+
+        friend constexpr QuantityValue operator-(
+            const QuantityValue &lhs,
+            const QuantityValue &rhs
+        ) noexcept {
+            return QuantityValue{lhs.m_value - rhs.m_value};
+        }
+
+        template <IsNumeric Scalar>
+        friend constexpr QuantityValue operator*(
+            const QuantityValue &lhs,
+            const Scalar rhs
+        ) noexcept {
+            return QuantityValue{lhs.m_value * static_cast<double>(rhs)};
+        }
+
+        template <IsNumeric Scalar>
+        friend constexpr QuantityValue operator*(
+            const Scalar lhs,
+            const QuantityValue &rhs
+        ) noexcept {
+            return QuantityValue{static_cast<double>(lhs) * rhs.m_value};
+        }
+
+        template <IsNumeric Scalar>
+        friend constexpr QuantityValue operator/(
+            const QuantityValue &lhs,
+            const Scalar rhs
+        ) noexcept {
+            return QuantityValue{lhs.m_value / static_cast<double>(rhs)};
+        }
+
+        template <IsNumeric Scalar>
+        friend constexpr std::partial_ordering operator<=>(
+            const QuantityValue &lhs,
+            const Scalar rhs
+        ) noexcept {
+            return lhs.m_value <=> static_cast<double>(rhs);
+        }
+
+        template <IsNumeric Scalar>
+        friend constexpr std::partial_ordering operator<=>(
+            const Scalar lhs,
+            const QuantityValue &rhs
+        ) noexcept {
+            return static_cast<double>(lhs) <=> rhs.m_value;
+        }
+
+        friend constexpr std::partial_ordering operator<=>(
+            const QuantityValue &lhs,
+            const QuantityValue &rhs
+        ) noexcept {
+            return lhs.m_value <=> rhs.m_value;
+        }
+
+    private:
+        double m_value;
+    };
+
+    template <typename Candidate>
+    struct IsQuantityValue : std::false_type { };
+
+    template <PhysicalQuantityType Quantity>
+    struct IsQuantityValue<QuantityValue<Quantity>> : std::true_type { };
+
+    template <typename Candidate>
+    concept QuantityValueType = IsQuantityValue<std::remove_cvref_t<Candidate>>::value;
+
+    template <typename Candidate>
+    struct QuantityOf;
+
+    template <PhysicalQuantityType Quantity>
+    struct QuantityOf<QuantityValue<Quantity>> {
+        using Type = Quantity;
+    };
+
+    template <QuantityValueType Value>
+    using QuantityOfT = typename QuantityOf<NeitherConstNorVolatileReference<Value>>::Type;
+
+    template <typename... Quantities>
+    struct QuantityList final {};
+
+    template <std::size_t Index, typename Quantities>
+    struct QuantityAt;
+
+    template <std::size_t Index, typename... Quantities>
+    struct QuantityAt<Index, QuantityList<Quantities...>> {
+        using Type = std::tuple_element_t<Index, std::tuple<Quantities...>>;
+    };
+}

src/include/serif/dimensions/type_alias.hpp

@@ -0,0 +1,44 @@
+#pragma once
+
+#include "serif/dimensions/quantities.hpp"
+#include "serif/dimensions/dimensions.hpp"
+
+namespace serif::dimensions {
+    using DimensionlessValue                             = QuantityValue<Dimensionless>;
+    using MassValue                                      = QuantityValue<Mass>;
+    using LengthValue                                    = QuantityValue<Length>;
+    using TimeValue                                      = QuantityValue<Time>;
+    using AreaValue                                      = QuantityValue<Area>;
+    using VolumeValue                                    = QuantityValue<Volume>;
+    using DensityValue                                   = QuantityValue<Density>;
+    using SurfaceDensityValue                            = QuantityValue<SurfaceDensity>;
+    using NumberDensityValue                             = QuantityValue<NumberDensity>;
+    using PressureValue                                  = QuantityValue<Pressure>;
+    using TemperatureValue                               = QuantityValue<Temperature>;
+    using EntropyValue                                   = QuantityValue<Entropy>;
+    using SpecificEntropyValue                           = QuantityValue<SpecificEntropy>;
+    using ChemicalPotentialValue                         = QuantityValue<ChemicalPotential>;
+    using EnergyValue                                    = QuantityValue<Energy>;
+    using InternalEnergyValue                            = QuantityValue<InternalEnergy>;
+    using SpecificEnergyValue                            = QuantityValue<SpecificEnergy>;
+    using SpecificInternalEnergyValue                    = QuantityValue<SpecificInternalEnergy>;
+    using SpecificEnthalpyValue                          = QuantityValue<SpecificEnthalpy>;
+    using EnergyDensityValue                             = QuantityValue<EnergyDensity>;
+    using GravitationalPotentialValue                    = QuantityValue<GravitationalPotential>;
+    using VelocityValue                                  = QuantityValue<Velocity>;
+    using AccelerationValue                              = QuantityValue<Acceleration>;
+    using FrequencyValue                                 = QuantityValue<Frequency>;
+    using AngularVelocityValue                           = QuantityValue<AngularVelocity>;
+    using MomentumValue                                  = QuantityValue<Momentum>;
+    using AngularMomentumValue                           = QuantityValue<AngularMomentum>;
+    using MomentOfInertiaValue                           = QuantityValue<MomentOfInertia>;
+    using ForceValue                                     = QuantityValue<Force>;
+    using TorqueValue                                    = QuantityValue<Torque>;
+    using PowerValue                                     = QuantityValue<Power>;
+    using LuminosityValue                                = QuantityValue<Luminosity>;
+    using MassFlowRateValue                              = QuantityValue<MassFlowRate>;
+    using OpacityValue                                   = QuantityValue<Opacity>;
+    using DynamicViscosityValue                          = QuantityValue<DynamicViscosity>;
+    using KinematicViscosityValue                        = QuantityValue<KinematicViscosity>;
+    using MagneticFluxDensityValue                       = QuantityValue<MagneticFluxDensity>;
+}

src/include/serif/discretization/blocks/base.hpp

@@ -0,0 +1,48 @@
+#pragma once
+
+namespace serif::discretization::blocks {
+    inline constexpr int DYNAMIC_BLOCK_SIZE = -1;
+
+    struct block {};
+
+    struct residual_block_base : block {
+        static constexpr int static_block_size = DYNAMIC_BLOCK_SIZE;
+    };
+
+    struct value_block_base : block {
+        static constexpr int static_block_size = DYNAMIC_BLOCK_SIZE;
+    };
+
+    template <typename GeneratedValue>
+    struct generated_value_block final : value_block_base {
+        using GeneratedType = GeneratedValue;
+        static constexpr int static_block_size = static_cast<int>(GeneratedValue::scalarArity);
+    };
+
+    template <typename GeneratedResidual>
+    struct generated_residual_block final : residual_block_base {
+        using GeneratedType = GeneratedResidual;
+        static constexpr int static_block_size = static_cast<int>(GeneratedResidual::scalarArity);
+    };
+
+    template <int index_value>
+    struct residual_block final : residual_block_base {
+        static constexpr int index = index_value;
+
+        // ReSharper disable once CppNonExplicitConversionOperator
+        constexpr operator int() const noexcept {
+            return index;
+        }
+    };
+
+    template <int index_value>
+    struct value_block final : value_block_base {
+        static constexpr int index = index_value;
+
+        // ReSharper disable once CppNonExplicitConversionOperator
+        constexpr operator int() const noexcept {
+            return index;
+        }
+    };
+
+}

src/include/serif/discretization/blocks/concepts.hpp

@@ -0,0 +1 @@
+#pragma once

src/include/serif/discretization/blocks/fields.hpp

@@ -0,0 +1,65 @@
+#pragma once
+
+#include "base.hpp"
+#include "serif/discretization/blocks/base.hpp"
+
+namespace serif::discretization::blocks {
+    struct field {};
+    struct term {};
+
+    struct density final : field {
+        struct mass final : term {
+            struct value final : value_block_base {};
+            struct residual final : residual_block_base {};
+        };
+
+        static inline constexpr mass mass_term{};
+    };
+
+    struct displacement final : field {
+        struct geometry final : term {
+            struct value final : value_block_base {};
+            struct residual final : residual_block_base {};
+        };
+
+        static inline constexpr geometry geometry_term{};
+    };
+
+    struct surface_deformation final : field {
+        struct parameters final : term {
+            struct value final : value_block_base {};
+        };
+
+        struct shape_equilibrium final : term {
+            struct residual final : residual_block_base {};
+        };
+
+        static inline constexpr parameters parameters_term{};
+        static inline constexpr shape_equilibrium shape_equilibrium_term{};
+    };
+
+    struct gravity final : field {
+        struct gradient final : term {
+            struct value final : value_block_base {};
+            struct residual final : residual_block_base {};
+        };
+
+        struct potential final : term {
+            struct value final : value_block_base {};
+            struct residual final : residual_block_base {};
+        };
+
+        static inline constexpr potential potential_term{};
+        static inline constexpr gradient gradient_term{};
+    };
+
+    struct enthalpy final : field {
+        struct specific final : term {
+            struct value final : value_block_base { };
+            struct residual final : residual_block_base { };
+        };
+
+        static inline constexpr specific specific_term{};
+    };
+
+}

src/include/serif/discretization/blocks/rows.hpp

@@ -0,0 +1 @@
+#pragma once

src/include/serif/discretization/blocks/traits.hpp

@@ -0,0 +1 @@
+#pragma once

src/include/serif/discretization/domain/concepts.hpp

@@ -0,0 +1,55 @@
+#pragma once
+
+#include <cstddef>
+
+#include "serif/discretization/domain/types.hpp"
+#include "serif/utils/misc/std_helper/cleaning.hpp"
+#include "serif/utils/misc/std_helper/variadic_cleaning.hpp"
+
+namespace serif::discretization::domain {
+    using utils::misc::std_helper::AIsBaseClassOfB;
+    using utils::misc::std_helper::CountVariadicArguments;
+
+    template <typename T>
+    concept IsDomain = AIsBaseClassOfB<T, Domain>;
+
+    template <typename T>
+    concept IsBoundary = AIsBaseClassOfB<T, Boundary>;
+
+    // A DomainSet is several volumes treated as one. StellarDomains
+    // (core + envelope) is the motivating example.
+    template <IsDomain... DomainTs>
+    struct DomainSet {
+        static constexpr std::size_t count = CountVariadicArguments<DomainTs...>();
+    };
+
+    template <IsBoundary... BoundaryTs>
+    struct BoundarySet {
+        static constexpr std::size_t count = CountVariadicArguments<BoundaryTs...>();
+    };
+
+    template <typename T>
+    constexpr bool is_domain_set_v = false;
+
+    template <IsDomain... DomainTs>
+    constexpr bool is_domain_set_v<DomainSet<DomainTs...>> = true;
+
+    template <typename T>
+    concept IsDomainSet = is_domain_set_v<T>;
+
+    template <typename T>
+    concept IsDomainOrSet = IsDomain<T> || IsDomainSet<T>;
+
+    template <IsDomainOrSet... DomainTs>
+    struct DomainOrSetList {
+        static constexpr std::size_t count = CountVariadicArguments<DomainTs...>();
+    };
+
+    // "This volume is the one tagged by that DomainID" / "this surface is the one
+    // tagged by that BoundaryID". Used by the resolvers in resolver.hpp.
+    template <typename DomainT, typename DomainIDT>
+    concept DomainIsIdentifiedBy = utils::misc::std_helper::SameType<DomainT, typename DomainIDT::domain_type>;
+
+    template <typename BoundaryT, typename BoundaryIDT>
+    concept BoundaryIsIdentifiedBy = utils::misc::std_helper::SameType<BoundaryT, typename BoundaryIDT::boundary_type>;
+}

src/include/serif/discretization/domain/ids/boundary.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <string_view>
+
+#include "serif/discretization/domain/concepts.hpp"
+
+// BoundaryID binds a semantic surface (a Boundary) to the integer boundary
+// attribute MFEM stores in the mesh file.
+//
+//     BoundaryID<StellarSurfaceBoundary, 1>
+//         "boundary elements with attribute 1 are the stellar surface"
+namespace serif::discretization::domain::ids {
+    template <IsBoundary BoundaryT, int IdentifierValue>
+    struct BoundaryID {
+        using boundary_type = BoundaryT;
+        static constexpr int ID = IdentifierValue;
+    };
+
+    struct BoundaryDescriptor {
+        std::string_view name;
+        int ID;
+    };
+}

src/include/serif/discretization/domain/ids/concepts.hpp

@@ -0,0 +1,27 @@
+#pragma once
+
+#include "serif/discretization/domain/ids/boundary.hpp"
+#include "serif/discretization/domain/ids/domain.hpp"
+
+namespace serif::discretization::domain::ids {
+    template <typename T>
+    constexpr bool is_domain_id_v = false;
+
+    template <IsDomain DomainT, int IdentifierValue>
+    constexpr bool is_domain_id_v<DomainID<DomainT, IdentifierValue>> = true;
+
+    template <typename T>
+    constexpr bool is_boundary_id_v = false;
+
+    template <IsBoundary BoundaryT, int IdentifierValue>
+    constexpr bool is_boundary_id_v<BoundaryID<BoundaryT, IdentifierValue>> = true;
+
+    template <typename T>
+    concept IsDomainID = is_domain_id_v<T>;
+
+    template <typename T>
+    concept IsBoundaryID = is_boundary_id_v<T>;
+
+    template <typename T>
+    concept IsDomainOrBoundaryID = IsDomainID<T> || IsBoundaryID<T>;
+}

src/include/serif/discretization/domain/ids/domain.hpp

@@ -0,0 +1,24 @@
+#pragma once
+
+#include <string_view>
+
+#include "serif/discretization/domain/concepts.hpp"
+
+// DomainID binds a semantic volume (a Domain) to the integer element attribute
+// MFEM stores in the mesh file. In MFEM's own vocabulary this integer is the
+// element "attribute"; in a mesh generator it is usually the "material id".
+//
+//     DomainID<CoreDomain, 1>   "elements with attribute 1 are the core"
+namespace serif::discretization::domain::ids {
+    template <IsDomain DomainT, int IdentifierValue>
+    struct DomainID {
+        using domain_type = DomainT;
+        static constexpr int ID = IdentifierValue;
+    };
+
+    // Runtime-friendly (name, ID) pair, for logging and error messages.
+    struct DomainDescriptor {
+        std::string_view name;
+        int ID;
+    };
+}

src/include/serif/discretization/domain/ids/lists/concepts.hpp

@@ -0,0 +1,11 @@
+#pragma once
+
+#include "serif/discretization/domain/ids/lists/traits.hpp"
+
+namespace serif::discretization::domain::ids::lists {
+    template <typename T>
+    concept IsDomainIDList = is_domain_id_list_v<T>;
+
+    template <typename T>
+    concept IsBoundaryIDList = is_boundary_id_list_v<T>;
+}

src/include/serif/discretization/domain/ids/lists/lists.hpp

@@ -0,0 +1,31 @@
+#pragma once
+
+#include <array>
+#include <cstddef>
+
+#include "serif/discretization/domain/ids/validation/concepts.hpp"
+#include "serif/utils/misc/std_helper/variadic_cleaning.hpp"
+
+namespace serif::discretization::domain::ids::lists {
+    using utils::misc::std_helper::CountVariadicArguments;
+
+    template <IsDomainID... DomainIDTs>
+        requires validation::AreDistinctDomainIDs<DomainIDTs...>
+    struct DomainIDList {
+        static constexpr std::size_t count = CountVariadicArguments<DomainIDTs...>();
+
+        [[nodiscard]] static constexpr std::array<DomainDescriptor, count> descriptors() noexcept {
+            return {DomainDescriptor{.name = DomainIDTs::domain_type::name, .ID = DomainIDTs::ID}...};
+        }
+    };
+
+    template <IsBoundaryID... BoundaryIDTs>
+        requires validation::AreDistinctBoundaryIDs<BoundaryIDTs...>
+    struct BoundaryIDList {
+        static constexpr std::size_t count = CountVariadicArguments<BoundaryIDTs...>();
+
+        [[nodiscard]] static constexpr std::array<BoundaryDescriptor, count> descriptors() noexcept {
+            return {BoundaryDescriptor{.name = BoundaryIDTs::boundary_type::name, .ID = BoundaryIDTs::ID}...};
+        }
+    };
+}

src/include/serif/discretization/domain/ids/lists/traits.hpp

@@ -0,0 +1,17 @@
+#pragma once
+
+#include "serif/discretization/domain/ids/lists/lists.hpp"
+
+namespace serif::discretization::domain::ids::lists {
+    template <typename T>
+    constexpr bool is_domain_id_list_v = false;
+
+    template <IsDomainID... DomainIDTs>
+    constexpr bool is_domain_id_list_v<DomainIDList<DomainIDTs...>> = true;
+
+    template <typename T>
+    constexpr bool is_boundary_id_list_v = false;
+
+    template <IsBoundaryID... BoundaryIDTs>
+    constexpr bool is_boundary_id_list_v<BoundaryIDList<BoundaryIDTs...>> = true;
+}

src/include/serif/discretization/domain/ids/validation/concepts.hpp

@@ -0,0 +1,42 @@
+#pragma once
+
+#include "serif/discretization/domain/ids/concepts.hpp"
+#include "serif/discretization/domain/ids/validation/unique_ids.hpp"
+#include "serif/utils/misc/concepts/type_uniqueness.hpp"
+
+namespace serif::discretization::domain::ids::validation {
+    using utils::misc::concepts::AllTypesAreUnique;
+
+    // --- ID value validation ---
+    template <typename... DomainIDTs>
+    concept HaveUniqueDomainIDs = all_unique_values(DomainIDTs::ID...);
+
+    template <typename... BoundaryIDTs>
+    concept HaveUniqueBoundaryIDs = all_unique_values(BoundaryIDTs::ID...);
+
+    // --- Type validation (Domains & Boundaries) ---
+    template <typename... DomainIDTs>
+    concept HaveUniqueDomains = AllTypesAreUnique<typename DomainIDTs::domain_type...>;
+
+    template <typename... BoundaryIDTs>
+    concept HaveUniqueBoundaries = AllTypesAreUnique<typename BoundaryIDTs::boundary_type...>;
+
+    // A well formed ID table: every entry is really an ID, no integer is reused,
+    // and no volume (or surface) is tagged twice.
+    //
+    // Note that these are variadic concepts, therefore they must be applied to a pack as a whole
+    // via a requires-clause. Writing them as a template parameter constraint
+    // (e.g. template <AreDistinctDomainIDs... Ts>) will apply them to each element
+    // individually, this checks nothing.
+    template <typename... DomainIDTs>
+    concept AreDistinctDomainIDs =
+        (IsDomainID<DomainIDTs> && ...) &&      // 1. everything in the pack really is a DomainID
+        HaveUniqueDomainIDs<DomainIDTs...> &&   // 2. no repeated integer attribute
+        HaveUniqueDomains<DomainIDTs...>;       // 3. no volume tagged twice
+
+    template <typename... BoundaryIDTs>
+    concept AreDistinctBoundaryIDs =
+        (IsBoundaryID<BoundaryIDTs> && ...) &&      // 1. everything in the pack really is a BoundaryID
+        HaveUniqueBoundaryIDs<BoundaryIDTs...> &&   // 2. no repeated integer attribute
+        HaveUniqueBoundaries<BoundaryIDTs...>;      // 3. no surface tagged twice
+}

src/include/serif/discretization/domain/ids/validation/unique_ids.hpp

@@ -0,0 +1,29 @@
+#pragma once
+
+#include <array>
+#include <concepts>
+#include <type_traits>
+
+#include "serif/utils/misc/std_helper/variadic_cleaning.hpp"
+
+namespace serif::discretization::domain::ids::validation {
+    // Returns true when every value handed to it is distinct.
+    template <typename... Ts>
+        requires (std::equality_comparable<Ts> && ...)
+    [[nodiscard]] consteval bool all_unique_values(Ts... values) noexcept {
+        constexpr std::size_t valueCount = utils::misc::std_helper::CountVariadicArguments<Ts...>();
+
+        if constexpr (valueCount < 2) {
+            return true;
+        } else {
+            const std::array<std::common_type_t<Ts...>, valueCount> vals{values...};
+
+            for (std::size_t i = 0; i < vals.size(); ++i) {
+                for (std::size_t j = i + 1; j < vals.size(); ++j) {
+                    if (vals[i] == vals[j]) return false;
+                }
+            }
+            return true;
+        }
+    }
+}

src/include/serif/discretization/domain/mesh/topology.hpp

@@ -0,0 +1,56 @@
+#pragma once
+
+#include <cstddef>
+#include <vector>
+
+namespace mfem {
+    class Mesh;
+}
+
+// MeshTopology
+namespace serif::discretization::domain::mesh {
+    struct FaceElements {
+        int firstElementID{-1};
+        int secondElementID{-1};
+    };
+
+    class MeshTopology final {
+    public:
+        explicit MeshTopology(const mfem::Mesh &mesh);
+
+        [[nodiscard]] int element_count() const noexcept;
+        [[nodiscard]] int face_count() const noexcept;
+        [[nodiscard]] int boundary_element_count() const noexcept;
+
+        // The MFEM element attribute of a volume element, i.e. which DomainID
+        // the mesh claims this element belongs to.
+        [[nodiscard]] int element_domain_id(int elementID) const;
+
+        // The MFEM boundary attribute of a boundary element, i.e. which
+        // BoundaryID the mesh claims this boundary element belongs to.
+        [[nodiscard]] int boundary_element_boundary_id(int boundaryElementID) const;
+
+        [[nodiscard]] FaceElements face_elements(int faceID) const;
+
+        // Elements sharing an interior face with elementID. Not filtered by
+        // domain: callers that care about one domain filter as they walk.
+        [[nodiscard]] const std::vector<int> &element_neighbors(int elementID) const;
+
+        // The face a boundary element sits on, or -1 if MFEM reports none.
+        [[nodiscard]] int boundary_element_face(int boundaryElementID) const;
+
+        // A face may carry more than one boundary element
+        [[nodiscard]] const std::vector<int> &boundary_elements_on_face(int faceID) const;
+
+    private:
+        std::vector<int> m_elementDomainIDs;
+        std::vector<FaceElements> m_faceElements;
+        std::vector<std::vector<int>> m_elementNeighbors;
+
+        std::vector<int> m_boundaryElementBoundaryIDs;
+        std::vector<int> m_boundaryElementFaceIDs;
+        std::vector<std::vector<int>> m_boundaryElementsByFace;
+
+        static const std::vector<int> s_noElements;
+    };
+}

src/include/serif/discretization/domain/physical_domains.hpp

@@ -0,0 +1,11 @@
+#pragma once
+
+#include "serif/discretization/domain/types.hpp"
+#include "serif/discretization/domain/concepts.hpp"
+
+namespace serif::discretization::domain {
+    // Sets are named in the plural to distinguish them from the single volumes
+    // in types.hpp.
+    using StellarDomains = DomainSet<CoreDomain, EnvelopeDomain>;
+    using AllDomains     = DomainSet<CoreDomain, EnvelopeDomain, VacuumDomain>;
+}

src/include/serif/discretization/domain/relation/concepts.hpp

@@ -0,0 +1,11 @@
+#pragma once
+
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/utils/misc/std_helper/cleaning.hpp"
+
+namespace serif::discretization::domain::relation {
+    using utils::misc::std_helper::AIsBaseClassOfB;
+
+    template <typename R>
+    concept IsRelation = AIsBaseClassOfB<R, DomainRelation>;
+}

src/include/serif/discretization/domain/relation/lists/concepts.hpp

@@ -0,0 +1,8 @@
+#pragma once
+
+#include "serif/discretization/domain/relation/lists/traits.hpp"
+
+namespace serif::discretization::domain::relation::lists {
+    template <typename RL>
+    concept IsRelationList = is_relation_list_v<RL>;
+}

src/include/serif/discretization/domain/relation/lists/relation_list.hpp

@@ -0,0 +1,15 @@
+#pragma once
+
+#include <cstddef>
+
+#include "serif/discretization/domain/relation/concepts.hpp"
+#include "serif/utils/misc/std_helper/variadic_cleaning.hpp"
+
+namespace serif::discretization::domain::relation::lists {
+    using utils::misc::std_helper::CountVariadicArguments;
+
+    template <IsRelation... RelationTs>
+    struct RelationList {
+        static constexpr std::size_t count = CountVariadicArguments<RelationTs...>();
+    };
+}

src/include/serif/discretization/domain/relation/lists/traits.hpp

@@ -0,0 +1,11 @@
+#pragma once
+
+#include "serif/discretization/domain/relation/lists/relation_list.hpp"
+
+namespace serif::discretization::domain::relation::lists {
+    template <typename T>
+    constexpr bool is_relation_list_v = false;
+
+    template <IsRelation... RelationTs>
+    constexpr bool is_relation_list_v<RelationList<RelationTs...>> = true;
+}

src/include/serif/discretization/domain/relation/relations.hpp

@@ -0,0 +1,49 @@
+#pragma once
+
+#include <string_view>
+
+#include "serif/discretization/domain/concepts.hpp"
+#include "serif/utils/misc/concepts/numeric.hpp"
+#include "serif/utils/misc/std_helper/variadic_cleaning.hpp"
+
+// The relations a schema can declare about a mesh. Each one is a statement that
+// must hold for the mesh to be usable, and each has a matching runtime check in
+// schema/validation/.
+namespace serif::discretization::domain::relation {
+    using utils::misc::concepts::SetOfOneOrTwo;
+    using utils::misc::std_helper::CountVariadicArguments;
+
+    struct DomainRelation { };
+
+    // "Every face of A that is not interior to A borders B."
+    template <IsDomainOrSet A, IsDomainOrSet B>
+    struct Inscribed final : DomainRelation {
+        using inner_type = A;
+        using outer_type = B;
+
+        static constexpr std::string_view name = "inscribed";
+    };
+
+    // "A is a single connected blob of elements, not two separate pieces."
+    template <IsDomainOrSet A>
+    struct FullyConnected final : DomainRelation {
+        using domain_type = A;
+
+        static constexpr std::string_view name = "fully_connected";
+    };
+
+    // "BoundaryT is exactly the surface between these one or two volumes."
+    //
+    // One domain  -> the surface between that volume and the outside of the mesh.
+    // Two domains -> the interface between the two volumes.
+    template <IsBoundary BoundaryT, IsDomainOrSet... DomainTs>
+        requires SetOfOneOrTwo<DomainTs...>
+    struct DomainBoundary final : DomainRelation {
+        using boundary_type = BoundaryT;
+        using domain_types  = DomainOrSetList<DomainTs...>;
+
+        static constexpr std::size_t domain_count = CountVariadicArguments<DomainTs...>();
+
+        static constexpr std::string_view name = "domain_boundary";
+    };
+}

src/include/serif/discretization/domain/relation/validation/concepts.hpp

@@ -0,0 +1,8 @@
+#pragma once
+
+#include "serif/discretization/domain/relation/validation/validation.hpp"
+
+namespace serif::discretization::domain::relation::validation {
+    template <typename DomainIDsT, typename BoundaryIDsT, typename RelationsT>
+    concept HaveValidRelationEntities = RelationsUseRegisteredEntities<DomainIDsT, BoundaryIDsT, RelationsT>::value;
+}

src/include/serif/discretization/domain/relation/validation/runtime.hpp

@@ -0,0 +1,73 @@
+#pragma once
+
+#include <cstdint>
+#include <optional>
+
+namespace serif::discretization::domain::relation::validation {
+    enum class RelationValidationFailure : std::uint8_t {
+        None,
+
+        // FullyConnected<DomainT>
+        DomainAbsent,
+        DomainDisconnected,
+
+        // Inscribed<InnerT, OuterT>
+        InnerDomainAbsent,
+        OuterDomainAbsent,
+        InnerDomainHasNoBoundary,
+        InnerDomainTouchesMeshBoundary,
+        InnerDomainTouchesUnexpectedDomain,
+
+        // DomainBoundary<BoundaryT, DomainTs...>
+        DomainBoundaryAbsent,
+        DomainBoundaryTaggedFaceHasWrongTopology,
+        DomainBoundaryTaggedFaceTouchesUnexpectedDomain,
+        DomainBoundaryExpectedFaceIsUntagged,
+        DomainBoundaryExpectedFaceHasWrongID
+    };
+
+    struct RelationValidationResult {
+        RelationValidationFailure failure{RelationValidationFailure::None};
+
+        struct InscribedDiagnostics {
+            int faceID{-1};
+            int innerElementID{-1};
+            int adjacentElementID{-1};
+            int adjacentDomainID{-1};
+        };
+
+        std::optional<InscribedDiagnostics> inscribedDiagnostics = std::nullopt;
+
+        struct ConnectedDiagnostics {
+            int elementID{-1};
+            int domainElementCount{0};
+            int visitedElementCount{0};
+        };
+
+        std::optional<ConnectedDiagnostics> connectedDiagnostics = std::nullopt;
+
+        struct DomainBoundaryDiagnostics {
+            int faceID{-1};
+            int boundaryElementID{-1};
+
+            int expectedBoundaryID{0};
+            std::optional<int> actualBoundaryID = std::nullopt;
+
+            int firstElementID{-1};
+            int secondElementID{-1};
+
+            std::optional<int> firstDomainID  = std::nullopt;
+            std::optional<int> secondDomainID = std::nullopt;
+        };
+
+        std::optional<DomainBoundaryDiagnostics> domainBoundaryDiagnostics = std::nullopt;
+
+        [[nodiscard]] bool valid() const noexcept {
+            return failure == RelationValidationFailure::None;
+        }
+
+        [[nodiscard]] explicit operator bool() const noexcept {
+            return valid();
+        }
+    };
+}

src/include/serif/discretization/domain/relation/validation/validation.hpp

@@ -0,0 +1,47 @@
+#pragma once
+
+#include <type_traits>
+
+#include "serif/discretization/domain/concepts.hpp"
+#include "serif/discretization/domain/ids/lists/concepts.hpp"
+#include "serif/discretization/domain/relation/concepts.hpp"
+#include "serif/discretization/domain/relation/lists/concepts.hpp"
+#include "serif/discretization/domain/relation/lists/relation_list.hpp"
+#include "serif/discretization/domain/resolver.hpp"
+
+// Compile-time check that every volume and surface named by a relation is
+// actually registered in the schema's ID tables.
+namespace serif::discretization::domain::relation::validation {
+    using ids::lists::IsDomainIDList;
+    using ids::lists::IsBoundaryIDList;
+    using lists::IsRelationList;
+
+    // This header is primarily split into a singular and a plural section. The plural section simply folds over a pack to invoke the singular section.
+
+    template <IsRelation RelationT, IsDomainIDList DomainIDsT, IsBoundaryIDList BoundaryIDsT>
+    struct RelationUsesRegisteredEntities;
+
+    template <IsDomainOrSet DomainT, IsDomainIDList DomainIDsT, IsBoundaryIDList BoundaryIDsT>
+    struct RelationUsesRegisteredEntities<FullyConnected<DomainT>, DomainIDsT, BoundaryIDsT> :
+        std::bool_constant<DomainIDResolver<DomainT, DomainIDsT>::registered> { };
+
+    template <IsDomainOrSet InnerDomainT, IsDomainOrSet OuterDomainT, IsDomainIDList DomainIDsT, IsBoundaryIDList BoundaryIDsT>
+    struct RelationUsesRegisteredEntities<Inscribed<InnerDomainT, OuterDomainT>, DomainIDsT, BoundaryIDsT> :
+        std::bool_constant<
+            DomainIDResolver<InnerDomainT, DomainIDsT>::registered &&
+            DomainIDResolver<OuterDomainT, DomainIDsT>::registered> { };
+
+    //Note that the fold below is over DomainTs (the volumes the relation names).
+    template <IsBoundary BoundaryT, IsDomainOrSet... DomainTs, IsDomainIDList DomainIDsT, IsBoundaryIDList BoundaryIDsT>
+    struct RelationUsesRegisteredEntities<DomainBoundary<BoundaryT, DomainTs...>, DomainIDsT, BoundaryIDsT> :
+        std::bool_constant<
+            BoundaryIDResolver<BoundaryT, BoundaryIDsT>::registered &&
+            (DomainIDResolver<DomainTs, DomainIDsT>::registered && ...)> { };
+
+    template <IsDomainIDList DomainIDsT, IsBoundaryIDList BoundaryIDsT, typename RelationListT>
+    struct RelationsUseRegisteredEntities;
+
+    template <IsDomainIDList DomainIDsT, IsBoundaryIDList BoundaryIDsT, IsRelation... RelationTs>
+    struct RelationsUseRegisteredEntities<DomainIDsT, BoundaryIDsT, lists::RelationList<RelationTs...>> :
+        std::bool_constant<(RelationUsesRegisteredEntities<RelationTs, DomainIDsT, BoundaryIDsT>::value && ...)> { };
+}

src/include/serif/discretization/domain/resolver.hpp

@@ -0,0 +1,94 @@
+#pragma once
+
+#include "serif/discretization/domain/concepts.hpp"
+#include "serif/discretization/domain/ids/lists/concepts.hpp"
+#include "serif/discretization/domain/ids/lists/lists.hpp"
+
+namespace serif::discretization::domain {
+    template <typename DomainT, typename DomainIDListT>
+    struct DomainIDResolver;
+
+
+    // In the following two structs there is a lot of folding, again (I am writing this comment after the comment regarding folding in utils/misc/traits/type_uniqueness.hpp was written).
+    // Anyway a brief overview of the folding in these structs follow
+    // # For the single domain specialization:
+    //     1. We make registered using a fold that checks that the domain is tagged by at least one of the IDs in the list.
+    //     2. contains_id then checks that both the domain is tagged by at least one entry and that the entry tagging the domain has the ID value that is passed in.
+    //        note that this latter check is why we cannot simply check registered again.
+    //     3. id first validates that the domain is registered properly (e.g. the domain is tagged by at least one of the IDs in the list) and
+    //        then it uses a fold to return the ID value of the entry that tags the domain using a ternary. Read this as looping over all entries in the
+    //        list, then if the entry tags the domain then accumulate that ID, otherwise accumulate 0. Note that MFEM convention has it that element attributes start
+    //        at one. This is a legitimate safety concern we may want to address at some point given that there is nothing stopping someone from making a mesh with
+    //        an attribute of 0.
+    // # For the multiple domain specialization:
+    //     1. We make a helper alias which lets us query the single domain resolver for each domain in the set.
+    //     2. registered therefore asks if all of the domains in the set are registered with the ID list by folding over all domains and then asking the single domain resolver if that domain is registered with the ID list.
+    //     3. contains_id does basically the same thing as registered, but instead of asking if the domain is registered it asks if the domain contains the ID that is passed in.
+
+    // A specialization to resolve the element attribute for a single Domain
+    template <IsDomain DomainT, ids::IsDomainID... DomainIDTs>
+    struct DomainIDResolver<DomainT, ids::lists::DomainIDList<DomainIDTs...>> {
+        static constexpr bool registered = (DomainIsIdentifiedBy<DomainT, DomainIDTs> || ...);
+
+        [[nodiscard]] static constexpr bool contains_id(const int domainID) noexcept {
+            return ((DomainIsIdentifiedBy<DomainT, DomainIDTs> && DomainIDTs::ID == domainID) || ...);
+        }
+
+        [[nodiscard]] static consteval int id() {
+            static_assert(registered, "Requested domain is not registered in this schema.");
+            return ((DomainIsIdentifiedBy<DomainT, DomainIDTs> ? DomainIDTs::ID : 0) + ...);
+        }
+    };
+
+    // This specialization resolves the element attributes for a set of Domains
+    template <IsDomain... DomainTs, ids::IsDomainID... DomainIDTs>
+    struct DomainIDResolver<DomainSet<DomainTs...>, ids::lists::DomainIDList<DomainIDTs...>> {
+
+        template <typename D>
+        using ResolverFor = DomainIDResolver<D, ids::lists::DomainIDList<DomainIDTs...>>;
+
+        static constexpr bool registered = (ResolverFor<DomainTs>::registered && ...);
+
+        [[nodiscard]] static constexpr bool contains_id(const int domainID) noexcept {
+            return (ResolverFor<DomainTs>::contains_id(domainID) || ...);
+        }
+    };
+
+
+    template <typename BoundaryT, typename BoundaryIDListT>
+    struct BoundaryIDResolver;
+
+
+    // MFEM has no pleasant and clear way to
+    // distinguish the *idea* of a boundary from the integer attribute tagging it.
+    // We elect to use the vocabulary Boundary for the idea and BoundaryID for the MFEM integer (and
+    // Domain / DomainID on the volume side), so the two sides read symmetrically. Generally in MFEM these would be called Boundary and BoundaryAttribute.
+    // Aside from that these follow the same logic as the DomainIDResolver
+    // specializations above, so I will not repeat the comments here.
+    template <IsBoundary BoundaryT, ids::IsBoundaryID... BoundaryIDTs>
+    struct BoundaryIDResolver<BoundaryT, ids::lists::BoundaryIDList<BoundaryIDTs...>> {
+        static constexpr bool registered = (BoundaryIsIdentifiedBy<BoundaryT, BoundaryIDTs> || ...);
+
+        [[nodiscard]] static constexpr bool contains_id(const int boundaryID) noexcept {
+            return ((BoundaryIsIdentifiedBy<BoundaryT, BoundaryIDTs> && BoundaryIDTs::ID == boundaryID) || ...);
+        }
+
+        [[nodiscard]] static consteval int id() {
+            static_assert(registered, "Requested boundary is not registered in this schema.");
+            return ((BoundaryIsIdentifiedBy<BoundaryT, BoundaryIDTs> ? BoundaryIDTs::ID : 0) + ...);
+        }
+    };
+
+    template <IsBoundary... BoundaryTs, ids::IsBoundaryID... BoundaryIDTs>
+    struct BoundaryIDResolver<BoundarySet<BoundaryTs...>, ids::lists::BoundaryIDList<BoundaryIDTs...>> {
+
+        template <typename B>
+        using ResolverFor = BoundaryIDResolver<B, ids::lists::BoundaryIDList<BoundaryIDTs...>>;
+
+        static constexpr bool registered = (ResolverFor<BoundaryTs>::registered && ...);
+
+        [[nodiscard]] static constexpr bool contains_id(const int boundaryID) noexcept {
+            return (ResolverFor<BoundaryTs>::contains_id(boundaryID) || ...);
+        }
+    };
+}

src/include/serif/discretization/domain/schema/concepts.hpp

@@ -0,0 +1,8 @@
+#pragma once
+
+#include "serif/discretization/domain/schema/traits.hpp"
+
+namespace serif::discretization::domain::schema {
+    template <typename T>
+    concept IsSchema = is_schema_v<T>;
+}

src/include/serif/discretization/domain/schema/domain_schema.hpp

@@ -0,0 +1,77 @@
+#pragma once
+
+#include <cstddef>
+
+#include "serif/discretization/domain/concepts.hpp"
+#include "serif/discretization/domain/ids/lists/concepts.hpp"
+#include "serif/discretization/domain/relation/lists/concepts.hpp"
+#include "serif/discretization/domain/relation/validation/concepts.hpp"
+#include "serif/discretization/domain/resolver.hpp"
+
+
+namespace serif::discretization::domain::schema {
+    using ids::lists::IsDomainIDList;
+    using ids::lists::IsBoundaryIDList;
+    using relation::lists::IsRelationList;
+
+    using relation::validation::HaveValidRelationEntities;
+
+    template <IsDomainIDList DomainIDs, IsBoundaryIDList BoundaryIDs, IsRelationList Relations>
+        requires HaveValidRelationEntities<DomainIDs, BoundaryIDs, Relations>
+    struct DomainSchema {
+        using domain_ids_type   = DomainIDs;
+        using boundary_ids_type = BoundaryIDs;
+        using relations_type    = Relations;
+
+        static constexpr std::size_t domain_id_count   = DomainIDs::count;
+        static constexpr std::size_t boundary_id_count = BoundaryIDs::count;
+        static constexpr std::size_t relation_count    = Relations::count;
+
+        [[nodiscard]] static constexpr auto domain_descriptors() noexcept {
+            return DomainIDs::descriptors();
+        }
+
+        [[nodiscard]] static constexpr auto boundary_descriptors() noexcept {
+            return BoundaryIDs::descriptors();
+        }
+
+        template <IsDomainOrSet DomainT>
+        [[nodiscard]] static consteval bool contains_domain() noexcept {
+            return DomainIDResolver<DomainT, DomainIDs>::registered;
+        }
+
+        // Does an element attribute read from a mesh belong to DomainT?
+        template <IsDomainOrSet DomainT>
+        [[nodiscard]] static constexpr bool domain_id_belongs_to(const int domainID) noexcept {
+            static_assert(contains_domain<DomainT>(), "DomainT is not registered in this schema.");
+
+            return DomainIDResolver<DomainT, DomainIDs>::contains_id(domainID);
+        }
+
+        // The single element attribute of a single volume.
+        template <IsDomain DomainT>
+        [[nodiscard]] static consteval int domain_id() noexcept {
+            static_assert(contains_domain<DomainT>(), "DomainT is not registered in this schema.");
+            return DomainIDResolver<DomainT, DomainIDs>::id();
+        }
+
+        template <IsBoundary BoundaryT>
+        [[nodiscard]] static consteval bool contains_boundary() noexcept {
+            return BoundaryIDResolver<BoundaryT, BoundaryIDs>::registered;
+        }
+
+        template <IsBoundary BoundaryT>
+        [[nodiscard]] static consteval int boundary_id() noexcept {
+            static_assert(contains_boundary<BoundaryT>(), "Requested boundary is not registered in this schema.");
+            return BoundaryIDResolver<BoundaryT, BoundaryIDs>::id();
+        }
+
+        // Does a boundary attribute read from a mesh belong to BoundaryT?
+        template <IsBoundary BoundaryT>
+        [[nodiscard]] static constexpr bool boundary_id_matches(const int boundaryID) noexcept {
+            static_assert(contains_boundary<BoundaryT>(), "Requested boundary is not registered in this schema.");
+
+            return BoundaryIDResolver<BoundaryT, BoundaryIDs>::contains_id(boundaryID);
+        }
+    };
+}

src/include/serif/discretization/domain/schema/schemas.hpp

@@ -0,0 +1,50 @@
+#pragma once
+
+#include "serif/discretization/domain/ids/boundary.hpp"
+#include "serif/discretization/domain/ids/domain.hpp"
+#include "serif/discretization/domain/ids/lists/lists.hpp"
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/relation/lists/relation_list.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/schema/domain_schema.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+namespace serif::discretization::domain::schema {
+    using ids::DomainID;
+    using ids::BoundaryID;
+
+    using ids::lists::DomainIDList;
+    using ids::lists::BoundaryIDList;
+
+    using relation::lists::RelationList;
+
+    // A star (core inside envelope) surrounded by a vacuum region, with a
+    // tagged stellar surface and a tagged outer surface.
+    using CoreEnvelopeVacuumDomainSchema = DomainSchema<
+        DomainIDList<
+            DomainID<CoreDomain, 1>,
+            DomainID<EnvelopeDomain, 2>,
+            DomainID<VacuumDomain, 3>
+        >,
+        BoundaryIDList<
+            BoundaryID<StellarSurfaceBoundary, 1>,
+            BoundaryID<InfinitySurfaceBoundary, 2>
+        >,
+        RelationList<
+            // No domain may arrive as two disconnected blobs.
+            relation::FullyConnected<CoreDomain>,
+            relation::FullyConnected<EnvelopeDomain>,
+            relation::FullyConnected<VacuumDomain>,
+
+            // The core sits inside the envelope, and the star as a whole sits
+            // inside the vacuum region.
+            relation::Inscribed<CoreDomain, EnvelopeDomain>,
+            relation::Inscribed<StellarDomains, VacuumDomain>,
+
+            // The stellar surface is the star/vacuum interface; the infinity
+            // surface is the outer edge of the vacuum region.
+            relation::DomainBoundary<StellarSurfaceBoundary, StellarDomains, VacuumDomain>,
+            relation::DomainBoundary<InfinitySurfaceBoundary, VacuumDomain>
+        >
+    >;
+}

src/include/serif/discretization/domain/schema/traits.hpp

@@ -0,0 +1,11 @@
+#pragma once
+
+#include "serif/discretization/domain/schema/domain_schema.hpp"
+
+namespace serif::discretization::domain::schema {
+    template <typename S>
+    constexpr bool is_schema_v = false;
+
+    template <IsDomainIDList DomainIDs, IsBoundaryIDList BoundaryIDs, IsRelationList Relations>
+    constexpr bool is_schema_v<DomainSchema<DomainIDs, BoundaryIDs, Relations>> = true;
+}

src/include/serif/discretization/domain/schema/utils.hpp

@@ -0,0 +1,24 @@
+#pragma once
+
+#include <mfem.hpp>
+
+#include "serif/discretization/domain/concepts.hpp"
+#include "serif/discretization/domain/schema/concepts.hpp"
+
+namespace serif::discretization::domain::schema {
+    // Build the 0/1 marker array MFEM wants when an integrator or an essential
+    // BC should apply only on part of the mesh.
+    template <IsDomainOrSet DomainT, IsSchema SchemaT>
+    [[nodiscard]] mfem::Array<int> make_attribute_marker(const mfem::Mesh &mesh) {
+        static_assert(
+            SchemaT::template contains_domain<DomainT>(),
+            "Requested marker domain is not completely registered in the supplied DomainSchema."
+        );
+
+        mfem::Array<int> marker(mesh.attributes.Max());
+        for (int domainID = 1; domainID <= marker.Size(); ++domainID) {
+            marker[domainID - 1] = SchemaT::template domain_id_belongs_to<DomainT>(domainID) ? 1 : 0;
+        }
+        return marker;
+    }
+}

src/include/serif/discretization/domain/schema/validation/all.hpp

@@ -0,0 +1,33 @@
+#pragma once
+
+#include <mfem.hpp>
+
+#include "serif/discretization/domain/mesh/topology.hpp"
+#include "serif/discretization/domain/schema/concepts.hpp"
+#include "serif/discretization/domain/schema/validation/results.hpp"
+#include "serif/discretization/domain/schema/validation/validator.hpp"
+
+// Every RelationValidator specialization must be visible before
+// SchemaRelationValidator is instantiated, so they are all included here and
+// validate_schema is defined only in this header. Include this file to validate
+// a mesh; including validator.hpp alone is not enough.
+#include "serif/discretization/domain/schema/validation/boundary.hpp"
+#include "serif/discretization/domain/schema/validation/connected.hpp"
+#include "serif/discretization/domain/schema/validation/inscribed.hpp"
+
+namespace serif::discretization::domain::schema::validation {
+    // Validate against connectivity that has already been extracted. Prefer
+    // this overload when validating several schemas against the same mesh.
+    template <IsSchema SchemaT>
+    [[nodiscard]] SchemaValidationResult validate_schema(const mesh::MeshTopology &topology) {
+        using RelationsT = typename SchemaT::relations_type;
+        return SchemaRelationValidator<SchemaT, RelationsT>::validate(topology);
+    }
+
+    // Extract the mesh connectivity once, then check every relation against it.
+    template <IsSchema SchemaT>
+    [[nodiscard]] SchemaValidationResult validate_schema(const mfem::Mesh &mfemMesh) {
+        const mesh::MeshTopology topology{mfemMesh};
+        return validate_schema<SchemaT>(topology);
+    }
+}

src/include/serif/discretization/domain/schema/validation/boundary.hpp

@@ -0,0 +1,229 @@
+#pragma once
+
+#include <optional>
+#include <tuple>
+#include <vector>
+
+#include "serif/discretization/domain/concepts.hpp"
+#include "serif/discretization/domain/mesh/topology.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/relation/validation/runtime.hpp"
+#include "serif/discretization/domain/schema/concepts.hpp"
+#include "serif/discretization/domain/schema/validation/validator.hpp"
+
+namespace serif::discretization::domain::schema::validation {
+    using relation::DomainBoundary;
+
+    using relation::validation::RelationValidationResult;
+    using relation::validation::RelationValidationFailure;
+
+    template <IsBoundary BoundaryT, IsDomainOrSet... DomainTs>
+    struct RelationValidator<DomainBoundary<BoundaryT, DomainTs...>> {
+        template <IsSchema SchemaT>
+        [[nodiscard]] static RelationValidationResult validate(const mesh::MeshTopology &topology) {
+            static_assert(
+                sizeof...(DomainTs) == 1 || sizeof...(DomainTs) == 2,
+                "DomainBoundary requires exactly one or two domains."
+            );
+
+            static_assert(
+                SchemaT::template contains_boundary<BoundaryT>(),
+                "DomainBoundary refers to a boundary which is not registered in the supplied DomainSchema."
+            );
+
+            static_assert(
+                (SchemaT::template contains_domain<DomainTs>() && ...),
+                "DomainBoundary refers to a domain which is not completely registered in the supplied DomainSchema."
+            );
+
+            constexpr int expectedBoundaryID = SchemaT::template boundary_id<BoundaryT>();
+
+            using DomainsTuple = std::tuple<DomainTs...>;
+
+            const auto make_diagnostics = [&](const int faceID, const int boundaryElementID, const std::optional<int> actualBoundaryID) {
+                RelationValidationResult::DomainBoundaryDiagnostics diagnostics{
+                    .faceID             = faceID,
+                    .boundaryElementID  = boundaryElementID,
+                    .expectedBoundaryID = expectedBoundaryID,
+                    .actualBoundaryID   = actualBoundaryID
+                };
+
+                if (faceID < 0 || faceID >= topology.face_count()) {
+                    return diagnostics;
+                }
+
+                const mesh::FaceElements faceElements = topology.face_elements(faceID);
+
+                diagnostics.firstElementID  = faceElements.firstElementID;
+                diagnostics.secondElementID = faceElements.secondElementID;
+
+                if (diagnostics.firstElementID >= 0) {
+                    diagnostics.firstDomainID = topology.element_domain_id(diagnostics.firstElementID);
+                }
+                if (diagnostics.secondElementID >= 0) {
+                    diagnostics.secondDomainID = topology.element_domain_id(diagnostics.secondElementID);
+                }
+
+                return diagnostics;
+            };
+
+            // Check only the cardinality/topological shape required by the relation.
+            //
+            // One-domain form:
+            //
+            //     Domain | computational exterior
+            //
+            // Exactly one adjacent volume element must exist.
+            //
+            // Two-domain form:
+            //
+            //     DomainA | DomainB
+            //
+            // Both adjacent volume elements must exist.
+            const auto has_required_topology = [](const int firstElementID, const int secondElementID) noexcept {
+                if constexpr (sizeof...(DomainTs) == 1) {
+                    const bool firstExists  = firstElementID >= 0;
+                    const bool secondExists = secondElementID >= 0;
+
+                    return firstExists != secondExists;
+                } else {
+                    return firstElementID >= 0 && secondElementID >= 0;
+                }
+            };
+
+            // Determine whether a face is exactly one of the faces described by
+            // DomainBoundary<BoundaryT, DomainTs...>. For two domains, ordering
+            // is intentionally irrelevant.
+            const auto face_matches_domains = [&topology](const int firstElementID, const int secondElementID) noexcept {
+                if constexpr (sizeof...(DomainTs) == 1) {
+                    using DomainT = std::tuple_element_t<0, DomainsTuple>;
+
+                    const bool firstExists  = firstElementID >= 0;
+                    const bool secondExists = secondElementID >= 0;
+
+                    if (firstExists == secondExists) {
+                        return false;
+                    }
+
+                    const int elementID = firstExists ? firstElementID : secondElementID;
+                    const int domainID  = topology.element_domain_id(elementID);
+
+                    return SchemaT::template domain_id_belongs_to<DomainT>(domainID);
+                } else {
+                    using FirstDomainT  = std::tuple_element_t<0, DomainsTuple>;
+                    using SecondDomainT = std::tuple_element_t<1, DomainsTuple>;
+
+                    if (firstElementID < 0 || secondElementID < 0) {
+                        return false;
+                    }
+
+                    const int firstDomainID  = topology.element_domain_id(firstElementID);
+                    const int secondDomainID = topology.element_domain_id(secondElementID);
+
+                    const bool forwardMatch = SchemaT::template domain_id_belongs_to<FirstDomainT>(firstDomainID) &&
+                                              SchemaT::template domain_id_belongs_to<SecondDomainT>(secondDomainID);
+
+                    const bool reverseMatch = SchemaT::template domain_id_belongs_to<SecondDomainT>(firstDomainID) &&
+                                              SchemaT::template domain_id_belongs_to<FirstDomainT>(secondDomainID);
+
+                    return forwardMatch || reverseMatch;
+                }
+            };
+
+            bool foundTaggedBoundary = false;
+
+            // Forward validation: every boundary element carrying BoundaryT must
+            // lie on exactly the interface declared by DomainBoundary.
+            for (int boundaryElementID = 0; boundaryElementID < topology.boundary_element_count(); ++boundaryElementID) {
+                const int boundaryID = topology.boundary_element_boundary_id(boundaryElementID);
+
+                if (boundaryID != expectedBoundaryID) {
+                    continue;
+                }
+
+                foundTaggedBoundary = true;
+
+                const int faceID = topology.boundary_element_face(boundaryElementID);
+
+                const auto [firstElementID, secondElementID] = topology.face_elements(faceID);
+
+                if (!has_required_topology(firstElementID, secondElementID)) {
+                    return {
+                        .failure = RelationValidationFailure::DomainBoundaryTaggedFaceHasWrongTopology,
+                        .domainBoundaryDiagnostics =
+                            std::make_optional<RelationValidationResult::DomainBoundaryDiagnostics>(
+                                make_diagnostics(faceID, boundaryElementID, boundaryID)
+                            )
+                    };
+                }
+
+                if (!face_matches_domains(firstElementID, secondElementID)) {
+                    return {
+                        .failure = RelationValidationFailure::DomainBoundaryTaggedFaceTouchesUnexpectedDomain,
+                        .domainBoundaryDiagnostics =
+                            std::make_optional<RelationValidationResult::DomainBoundaryDiagnostics>(
+                                make_diagnostics(faceID, boundaryElementID, boundaryID)
+                            )
+                    };
+                }
+            }
+
+            bool foundExpectedFace = false;
+
+            // Reverse validation: every face having the declared domain
+            // adjacency must carry BoundaryT.
+            for (int faceID = 0; faceID < topology.face_count(); ++faceID) {
+                const auto [firstElementID, secondElementID] = topology.face_elements(faceID);
+
+                if (!face_matches_domains(firstElementID, secondElementID)) {
+                    continue;
+                }
+
+                foundExpectedFace = true;
+
+                const std::vector<int> &boundaryElementIDs = topology.boundary_elements_on_face(faceID);
+
+                if (boundaryElementIDs.empty()) {
+                    return {
+                        .failure = RelationValidationFailure::DomainBoundaryExpectedFaceIsUntagged,
+                        .domainBoundaryDiagnostics =
+                            std::make_optional<RelationValidationResult::DomainBoundaryDiagnostics>(
+                                make_diagnostics(faceID, -1, std::nullopt)
+                            )
+                    };
+                }
+
+                for (const int boundaryElementID : boundaryElementIDs) {
+                    const int actualBoundaryID = topology.boundary_element_boundary_id(boundaryElementID);
+
+                    if (actualBoundaryID == expectedBoundaryID) {
+                        continue;
+                    }
+
+                    return {
+                        .failure = RelationValidationFailure::DomainBoundaryExpectedFaceHasWrongID,
+                        .domainBoundaryDiagnostics =
+                            std::make_optional<RelationValidationResult::DomainBoundaryDiagnostics>(
+                                make_diagnostics(faceID, boundaryElementID, actualBoundaryID)
+                            )
+                    };
+                }
+            }
+
+            // If neither a correctly tagged boundary nor a face having the
+            // required semantic topology exists, the declared DomainBoundary
+            // simply is not realized by this mesh.
+            if (!foundTaggedBoundary || !foundExpectedFace) {
+                return {
+                    .failure = RelationValidationFailure::DomainBoundaryAbsent,
+                    .domainBoundaryDiagnostics =
+                        std::make_optional<RelationValidationResult::DomainBoundaryDiagnostics>(
+                            make_diagnostics(-1, -1, std::nullopt)
+                        )
+                };
+            }
+
+            return {};
+        }
+    };
+}

src/include/serif/discretization/domain/schema/validation/connected.hpp

@@ -0,0 +1,117 @@
+#pragma once
+
+#include <cstddef>
+#include <optional>
+#include <vector>
+
+#include "serif/discretization/domain/concepts.hpp"
+#include "serif/discretization/domain/mesh/topology.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/relation/validation/runtime.hpp"
+#include "serif/discretization/domain/schema/concepts.hpp"
+#include "serif/discretization/domain/schema/validation/validator.hpp"
+
+// FullyConnected<DomainT>: walk the domain from one of its elements and check we
+// reach all of them.
+namespace serif::discretization::domain::schema::validation {
+    using relation::FullyConnected;
+
+    using relation::validation::RelationValidationResult;
+    using relation::validation::RelationValidationFailure;
+
+    template <IsDomainOrSet DomainT>
+    struct RelationValidator<FullyConnected<DomainT>> {
+        template <IsSchema SchemaT>
+        [[nodiscard]] static RelationValidationResult validate(const mesh::MeshTopology &topology) {
+            static_assert(SchemaT::template contains_domain<DomainT>(), "Domain not present, cannot determine connectedness.");
+
+            const int elementCount = topology.element_count();
+
+            // We should be careful here that I have not accidentally introduced an issue due to the specialized vector bool override. I think this is safe though.
+            std::vector<bool> belongsToDomain(static_cast<std::size_t>(elementCount), false);
+
+            int domainElementCount = 0;
+            int firstDomainElement = -1;
+
+            for (int elementID = 0; elementID < elementCount; ++elementID) {
+                const int domainID = topology.element_domain_id(elementID);
+                const bool belongs = SchemaT::template domain_id_belongs_to<DomainT>(domainID);
+
+                belongsToDomain[static_cast<std::size_t>(elementID)] = belongs;
+
+                if (!belongs) {
+                    continue;
+                }
+
+                ++domainElementCount;
+                if (firstDomainElement < 0) {
+                    firstDomainElement = elementID;
+                }
+            }
+
+            if (domainElementCount == 0) {
+                return {
+                    .failure = RelationValidationFailure::DomainAbsent,
+                    .connectedDiagnostics = std::make_optional<RelationValidationResult::ConnectedDiagnostics>({
+                        .domainElementCount  = 0,
+                        .visitedElementCount = 0
+                    })};
+            }
+
+            std::vector<bool> visited(static_cast<std::size_t>(elementCount), false);
+            std::vector<int> pending;
+
+            pending.reserve(static_cast<std::size_t>(domainElementCount));
+            pending.push_back(firstDomainElement);
+
+            int visitedElementCount = 0;
+
+            while (!pending.empty()) {
+                const int elementID = pending.back();
+                pending.pop_back();
+
+                if (visited[static_cast<std::size_t>(elementID)]) {
+                    continue;
+                }
+
+                visited[static_cast<std::size_t>(elementID)] = true;
+                ++visitedElementCount;
+
+                for (const int neighborElementID : topology.element_neighbors(elementID)) {
+                    const std::size_t neighborIndex = static_cast<std::size_t>(neighborElementID);
+
+                    if (belongsToDomain[neighborIndex] && !visited[neighborIndex]) {
+                        pending.push_back(neighborElementID);
+                    }
+                }
+            }
+
+            if (visitedElementCount == domainElementCount) {
+                return {
+                    .connectedDiagnostics = std::make_optional<RelationValidationResult::ConnectedDiagnostics>(
+                        {.domainElementCount = domainElementCount, .visitedElementCount = visitedElementCount}
+                    )
+                };
+            }
+
+            int disconnectedElementID = -1;
+
+            for (int elementID = 0; elementID < elementCount; ++elementID) {
+                const std::size_t index = static_cast<std::size_t>(elementID);
+
+                if (belongsToDomain[index] && !visited[index]) {
+                    disconnectedElementID = elementID;
+                    break;
+                }
+            }
+
+            return {
+                .failure = RelationValidationFailure::DomainDisconnected,
+                .connectedDiagnostics = std::make_optional<RelationValidationResult::ConnectedDiagnostics>({
+                    .elementID           = disconnectedElementID,
+                    .domainElementCount  = domainElementCount,
+                    .visitedElementCount = visitedElementCount
+                })};
+        }
+    };
+}

src/include/serif/discretization/domain/schema/validation/inscribed.hpp

@@ -0,0 +1,101 @@
+#pragma once
+
+#include <optional>
+
+#include "serif/discretization/domain/concepts.hpp"
+#include "serif/discretization/domain/mesh/topology.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/relation/validation/runtime.hpp"
+#include "serif/discretization/domain/schema/concepts.hpp"
+#include "serif/discretization/domain/schema/validation/validator.hpp"
+
+// Inscribed<InnerT, OuterT>: every face where the inner domain stops must have
+// the outer domain on the far side.
+namespace serif::discretization::domain::schema::validation {
+    using relation::Inscribed;
+
+    using relation::validation::RelationValidationResult;
+    using relation::validation::RelationValidationFailure;
+
+    template <IsDomainOrSet InnerT, IsDomainOrSet OuterT>
+    struct RelationValidator<Inscribed<InnerT, OuterT>> {
+        template <IsSchema SchemaT>
+        [[nodiscard]] static RelationValidationResult validate(const mesh::MeshTopology &topology) {
+            static_assert(
+                SchemaT::template contains_domain<InnerT>(), "The inner domain of an Inscribed relation is not present in the supplied schema. Inscribed cannot be enforced"
+            );
+
+            static_assert(
+                SchemaT::template contains_domain<OuterT>(), "The outer domain of an Inscribed relation is not present in the supplied schema. Inscribed cannot be enforced"
+            );
+
+            bool foundInnerElement  = false;
+            bool foundOuterElement  = false;
+            bool foundInnerBoundary = false;
+
+            for (int elementID = 0; elementID < topology.element_count(); ++elementID) {
+                const int domainID = topology.element_domain_id(elementID);
+
+                foundInnerElement = foundInnerElement || SchemaT::template domain_id_belongs_to<InnerT>(domainID);
+                foundOuterElement = foundOuterElement || SchemaT::template domain_id_belongs_to<OuterT>(domainID);
+            }
+
+            if (!foundInnerElement) {
+                return {.failure = RelationValidationFailure::InnerDomainAbsent};
+            }
+
+            if (!foundOuterElement) {
+                return {.failure = RelationValidationFailure::OuterDomainAbsent};
+            }
+
+            for (int faceID = 0; faceID < topology.face_count(); ++faceID) {
+                const mesh::FaceElements faceElements = topology.face_elements(faceID);
+
+                const int firstElementID  = faceElements.firstElementID;
+                const int secondElementID = faceElements.secondElementID;
+
+                const bool firstIsInner  = firstElementID >= 0 &&
+                    SchemaT::template domain_id_belongs_to<InnerT>(topology.element_domain_id(firstElementID));
+                const bool secondIsInner = secondElementID >= 0 &&
+                    SchemaT::template domain_id_belongs_to<InnerT>(topology.element_domain_id(secondElementID));
+
+                if (firstIsInner == secondIsInner) {
+                    continue;
+                }
+
+                foundInnerBoundary = true;
+
+                const int innerElementID    = firstIsInner ? firstElementID : secondElementID;
+                const int adjacentElementID = firstIsInner ? secondElementID : firstElementID;
+
+                if (adjacentElementID < 0) {
+                    return {
+                        .failure = RelationValidationFailure::InnerDomainTouchesMeshBoundary,
+                        .inscribedDiagnostics = std::make_optional<RelationValidationResult::InscribedDiagnostics>({
+                            .faceID         = faceID,
+                            .innerElementID = innerElementID,
+                        })};
+                }
+
+                const int adjacentDomainID = topology.element_domain_id(adjacentElementID);
+
+                if (!SchemaT::template domain_id_belongs_to<OuterT>(adjacentDomainID)) {
+                    return {
+                        .failure = RelationValidationFailure::InnerDomainTouchesUnexpectedDomain,
+                        .inscribedDiagnostics = std::make_optional<RelationValidationResult::InscribedDiagnostics>({
+                            .faceID            = faceID,
+                            .innerElementID    = innerElementID,
+                            .adjacentElementID = adjacentElementID,
+                            .adjacentDomainID  = adjacentDomainID
+                        })};
+                }
+            }
+
+            if (!foundInnerBoundary) {
+                return {.failure = RelationValidationFailure::InnerDomainHasNoBoundary};
+            }
+
+            return {};
+        }
+    };
+}

src/include/serif/discretization/domain/schema/validation/results.hpp

@@ -0,0 +1,58 @@
+#pragma once
+
+#include <cstddef>
+#include <optional>
+#include <string_view>
+#include <vector>
+#include <format>
+#include <string>
+
+#include "serif/discretization/domain/relation/validation/runtime.hpp"
+
+// The aggregate result of validating a whole schema against a mesh
+namespace serif::discretization::domain::schema::validation {
+    struct SchemaRelationValidationResult {
+        std::size_t relationIndex{0};
+        std::string_view relationName;
+        relation::validation::RelationValidationResult result;
+
+        [[nodiscard]] bool valid() const noexcept;
+        [[nodiscard]] explicit operator bool() const noexcept;
+    };
+
+    struct SchemaValidationResult {
+        std::vector<SchemaRelationValidationResult> relationResults;
+
+        [[nodiscard]] bool valid() const noexcept;
+        [[nodiscard]] explicit operator bool() const noexcept;
+        [[nodiscard]] std::size_t relation_count() const noexcept;
+        [[nodiscard]] std::size_t failed_relation_count() const noexcept;
+        [[nodiscard]] std::size_t passed_relation_count() const noexcept;
+        [[nodiscard]] std::optional<std::size_t> first_failed_relation_index() const noexcept;
+    };
+}
+
+
+template <>
+struct std::formatter<serif::discretization::domain::schema::validation::SchemaValidationResult> {
+    std::formatter<std::string> string_formatter;
+
+    constexpr auto parse(std::format_parse_context& ctx) {
+        return string_formatter.parse(ctx);
+    }
+
+    auto format(const serif::discretization::domain::schema::validation::SchemaValidationResult& result, std::format_context& ctx) const {
+        std::string output;
+        output += "SchemaValidationResult:\n";
+        output += std::format("  Valid: {}\n", result.valid());
+        output += std::format("  Relation count: {}\n", result.relation_count());
+        output += std::format("  Passed relation count: {}\n", result.passed_relation_count());
+        output += std::format("  Failed relation count: {}\n", result.failed_relation_count());
+        if (auto first_failed_index = result.first_failed_relation_index()) {
+            output += std::format("  First failed relation index: {}\n", *first_failed_index);
+        } else {
+            output += "  No failed relations.\n";
+        }
+        return string_formatter.format(output, ctx);
+    }
+};

src/include/serif/discretization/domain/schema/validation/validator.hpp

@@ -0,0 +1,49 @@
+#pragma once
+
+#include <cstddef>
+
+#include "serif/discretization/domain/mesh/topology.hpp"
+#include "serif/discretization/domain/relation/concepts.hpp"
+#include "serif/discretization/domain/relation/lists/relation_list.hpp"
+#include "serif/discretization/domain/relation/validation/runtime.hpp"
+#include "serif/discretization/domain/schema/concepts.hpp"
+#include "serif/discretization/domain/schema/validation/results.hpp"
+
+// The primary RelationValidator template, plus the fold that runs one validator
+// per declared relation.
+//
+// IMPORTANT: every relation kind is validated by an explicit specialization of
+// RelationValidator, and those specializations live in sibling headers
+// (inscribed.hpp, connected.hpp, boundary.hpp). Under headers, instantiating
+// SchemaRelationValidator in a translation unit that has not yet seen all of
+// those specializations is ill-formed with no diagnostic required. That is why
+// validate_schema lives in all.hpp, which includes this header and every
+// specialization. Include all.hpp, not this header, to validate a mesh.
+namespace serif::discretization::domain::schema::validation {
+    using relation::IsRelation;
+
+    template <IsRelation RelationT>
+    struct RelationValidator;
+
+    template <IsSchema SchemaT, typename RelationListT>
+    struct SchemaRelationValidator;
+
+    template <IsSchema SchemaT, IsRelation... RelationTs>
+    struct SchemaRelationValidator<SchemaT, relation::lists::RelationList<RelationTs...>> {
+
+        [[nodiscard]] static SchemaValidationResult validate(const mesh::MeshTopology &topology) {
+            SchemaValidationResult schemaResult;
+
+            schemaResult.relationResults.reserve(sizeof...(RelationTs));
+            std::size_t relationIndex = 0;
+
+            (schemaResult.relationResults.push_back(SchemaRelationValidationResult{
+                     .relationIndex = relationIndex++,
+                     .relationName  = RelationTs::name,
+                     .result        = RelationValidator<RelationTs>::template validate<SchemaT>(topology)
+                 }), ...);
+
+            return schemaResult;
+        }
+    };
+}

src/include/serif/discretization/domain/types.hpp

@@ -0,0 +1,32 @@
+#pragma once
+#include <string_view>
+
+// A domain is a volume while a boundary is
+// a surface. The mapping from these types to
+// the integer attributes stored in a mesh file lives in ids/ (DomainID and
+// BoundaryID) and defined by a schema.
+namespace serif::discretization::domain {
+    struct Domain { };
+
+    struct CoreDomain final : Domain {
+        static constexpr std::string_view name = "core";
+    };
+
+    struct EnvelopeDomain final : Domain {
+        static constexpr std::string_view name = "envelope";
+    };
+
+    struct VacuumDomain final : Domain {
+        static constexpr std::string_view name = "vacuum";
+    };
+
+    struct Boundary { };
+
+    struct StellarSurfaceBoundary final : Boundary {
+        static constexpr std::string_view name = "stellar_surface";
+    };
+
+    struct InfinitySurfaceBoundary final : Boundary {
+        static constexpr std::string_view name = "infinity_surface";
+    };
+}

src/include/serif/eos/evaluation.hpp

@@ -0,0 +1,32 @@
+#pragma once
+#include "serif/eos/models/concepts.hpp"
+
+namespace serif::eos {
+    /*
+     * I find the syntax that follows to be potentially confusing. Specifically the double use of noexcept. For those
+     * unfamiliar a brief explanation follows.
+     *
+     * The first noexcept is used to mark if the overall function (evaluate) is noexcept. However, noexcept can take a boolean. noexcept(true)
+     * tells the compiler that the function will not throw an exception (note if an exception is thrown in a noexcept function std::terminate will be called
+     * and no try-catch block will catch the exception). This allows for the compiler to generate more efficient code. noexcept(false) tells
+     * the compiler that the function may throw an exception (semantically it is equivalent to not including noexcept). The second noexcept
+     * (e.g. noexcept(noexcept()) is used to determine what that boolean is. This is a confusing part of C++ syntax where the same
+     * keyword has different semantic meanings depending on context. In this case the inner noexcept is used
+     * to determine if the expression inside it is noexcept. That is to say if eos.evaluate(...) is noexcept for the particular
+     * templates EOS argument then this function will be noexcept. If eos.evaluate(...) is not noexcept then this function will not be noexcept.
+     */
+    template <dimensions::ThermodynamicQuantityType OutputQ, models::EOSModel EOS, dimensions::QuantityValueType... InputValues>
+        requires models::EOSSupportingRelation<EOS, relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>>
+    [[nodiscard]] constexpr dimensions::QuantityValue<OutputQ> evaluate(const EOS& eos, const InputValues... inputValues)
+    noexcept(noexcept(eos.evaluate(relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>{}, inputValues...))) {
+        return eos.evaluate(relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>{}, inputValues...);
+    };
+
+    template <dimensions::ThermodynamicQuantityType OutputQ, dimensions::ThermodynamicQuantityType InputQ, models::EOSModel EOS, dimensions::QuantityValueType... InputValues>
+        requires models::EOSSupportingPartialDerivative<EOS, relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>, InputQ>
+    [[nodiscard]] constexpr dimensions::PartialDerivative<OutputQ, InputQ> partial_derivative(const EOS& eos, const InputValues... inputValues)
+    noexcept(noexcept(eos.partialDerivative(relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>{}, dimensions::WithRespectTo<InputQ>{}, inputValues...))) {
+        return eos.partial_derivative(relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>{}, dimensions::WithRespectTo<InputQ>{}, inputValues...);
+    }
+
+}

src/include/serif/eos/exceptions.hpp

@@ -0,0 +1,32 @@
+#pragma once
+#include <cstdint>
+#include <stdexcept>
+#include <string>
+
+namespace serif::eos {
+    enum class EOSEvaluationErrorCode : uint8_t {
+        unsupported_relation,
+        unsupported_derivative,
+        wrong_input_count,
+        wrong_input_quantity,
+        nonfinite_input,
+        outside_domain,
+        nonfinite_result,
+        invalid_thermodynamic_input
+    };
+
+    class EOSEvaluationError final : public std::domain_error {
+    public:
+        explicit EOSEvaluationError(
+            const EOSEvaluationErrorCode code,
+            std::string message ):
+        std::domain_error(std::move(message)),
+        m_code(code) {}
+
+        [[nodiscard]] EOSEvaluationErrorCode code() const noexcept {
+            return m_code;
+        }
+    private:
+        EOSEvaluationErrorCode m_code;
+    };
+}

src/include/serif/eos/models/concepts.hpp

@@ -0,0 +1,28 @@
+#pragma once
+
+#include "serif/eos/models/traits.hpp"
+
+namespace serif::eos::models {
+    template <typename EOSCandidate>
+    concept EOSModel = is_equation_of_state_model_v<EOSCandidate>;
+
+    template <typename EOS, typename RelationT>
+    concept EOSSupportingRelation = supports_relation_v<EOS, RelationT>;
+
+    template <typename EOS, typename RelationT, typename InputQuantity>
+    concept EOSSupportingPartialDerivative = supports_partial_derivative_v<EOS, RelationT, InputQuantity>;
+
+    template <typename EOSCandidate>
+    concept StructureSeedEOS = EOSModel<EOSCandidate> &&
+        EOSSupportingRelation<EOSCandidate, relations::SpecificEnthalpyFromDensity>;
+
+    template <typename EOSCandidate>
+    concept BarotropicClosureEOS = EOSModel<EOSCandidate> &&
+        EOSSupportingRelation<EOSCandidate, relations::DensityFromSpecificEnthalpy> &&
+        EOSSupportingPartialDerivative<EOSCandidate, relations::DensityFromSpecificEnthalpy, dimensions::SpecificEnthalpy>;
+
+    template <typename EOSCandidate>
+    concept PressureForceEOS = EOSModel<EOSCandidate> &&
+        EOSSupportingRelation<EOSCandidate, relations::PressureFromSpecificEnthalpy> &&
+        EOSSupportingPartialDerivative<EOSCandidate, relations::PressureFromSpecificEnthalpy, dimensions::SpecificEnthalpy>;
+}

src/include/serif/eos/models/polytropic.hpp

@@ -0,0 +1,67 @@
+#pragma once
+#include <format>
+#include "serif/eos/relations/relations.hpp"
+#include "serif/dimensions/type_alias.hpp"
+#include "serif/dimensions/partial.hpp"
+#include "serif/utils/misc/finite.hpp"
+
+namespace serif::eos::models {
+    struct PolytropeParameters final {
+        double n;
+        double K;
+    };
+
+    class Polytrope final {
+    public:
+        /* The idea here (and with all EOS) is that we force the developer / implementer to enumerate all the thermodynamic relations
+         * that the EOS supports explicitly. This is a small amount of boilerplate; however, it enables effectively a fully automatic
+         * EOS compiler down the line. Much of the EOS system these classes depend on (relations, models, the system in evaluation) are all
+         * there to allow for the compile-time checking of the EOS model and its relations and to allow the developer cross-section
+         * of that to be limited to defining this relation catalog. Every EOS implemented should start (at the very top of the class
+         * with a Relations alias in the public section of the class. Any EOS that does not have this will not compile.
+         */
+        using Relations = relations::RelationCatalog<
+            relations::PressureFromDensity,
+            relations::PressureFromSpecificEnthalpy,
+            relations::SpecificEnthalpyFromDensity,
+            relations::SpecificEnthalpyFromPressure,
+            relations::DensityFromSpecificEnthalpy
+        >;
+
+        explicit Polytrope(const PolytropeParameters& parameters);
+        explicit Polytrope(double n, double K);
+
+        [[nodiscard]] double polytropic_index() const noexcept;
+        [[nodiscard]] double polytropic_constant() const noexcept;
+        [[nodiscard]] double enthalpy_scale() const noexcept;
+
+        [[nodiscard]] dimensions::PressureValue evaluate(relations::PressureFromDensity, dimensions::DensityValue density) const;
+        [[nodiscard]] dimensions::PressureValue evaluate(relations::PressureFromSpecificEnthalpy, dimensions::SpecificEnthalpyValue specific_enthalpy) const;
+        [[nodiscard]] dimensions::SpecificEnthalpyValue evaluate(relations::SpecificEnthalpyFromDensity, dimensions::DensityValue density) const;
+        [[nodiscard]] dimensions::SpecificEnthalpyValue evaluate(relations::SpecificEnthalpyFromPressure, dimensions::PressureValue pressure) const;
+        [[nodiscard]] dimensions::DensityValue evaluate(relations::DensityFromSpecificEnthalpy, dimensions::SpecificEnthalpyValue specific_enthalpy) const;
+
+        [[nodiscard]] dimensions::PartialDerivative<dimensions::Density, dimensions::SpecificEnthalpy> partial_derivative(
+            relations::DensityFromSpecificEnthalpy,
+            dimensions::WithRespectTo<dimensions::SpecificEnthalpy>,
+            dimensions::SpecificEnthalpyValue specific_enthalpy
+        ) const;
+
+        [[nodiscard]] dimensions::PartialDerivative<dimensions::Pressure, dimensions::SpecificEnthalpy> partial_derivative(
+            relations::PressureFromSpecificEnthalpy,
+            dimensions::WithRespectTo<dimensions::SpecificEnthalpy>,
+            dimensions::SpecificEnthalpyValue specific_enthalpy
+        ) const;
+
+        [[nodiscard]] dimensions::PartialDerivative<dimensions::Pressure, dimensions::Density> partial_derivative(
+            relations::PressureFromDensity,
+            dimensions::WithRespectTo<dimensions::Density>,
+            dimensions::DensityValue density
+        ) const;
+
+    private:
+        double m_polytropic_index;
+        double m_polytropic_constant;
+        double m_enthalpy_scale;
+    };
+}

src/include/serif/eos/models/traits.hpp

@@ -0,0 +1,56 @@
+#pragma once
+
+#include "serif/dimensions/quantities.hpp"
+#include "serif/eos/relations/concepts.hpp"
+
+namespace serif::eos::models {
+    template <typename EOS, typename RelationT>
+    struct ImplementsRelation_t : std::false_type{};
+
+    template <typename EOS, typename Output, typename... Inputs>\
+    struct ImplementsRelation_t<EOS, relations::Relation<Output, Inputs...>> :
+        std::bool_constant<
+            requires(const std::remove_cvref_t<EOS> &eos, dimensions::QuantityValue<Inputs>... inputValues) {
+        {
+            eos.evaluate(relations::Relation<Output, Inputs...>{}, inputValues...)
+        } -> std::same_as<dimensions::QuantityValue<Output>>;}>{}; // Read this as a requirement that when evaluating the templates EOS with the given relation and inputs the returned value is of the correct output type.
+
+    template <typename EOS, typename Catalog>
+    struct ImplementsRelationCatalog_t : std::false_type {};
+
+    template <typename EOS, typename... Relations>
+    struct ImplementsRelationCatalog_t<EOS, relations::RelationCatalog<Relations...>> :
+        std::bool_constant<
+            (ImplementsRelation_t<EOS, Relations>::value && ...)>{}; // Fold over all relations in a catalog and ensure the EOS implements all of them.
+
+    template <typename EOSCandidate, typename = void>
+    struct IsEquationOfStateModel_t : std::false_type {};
+
+    template <typename EOSCandidate>
+    struct IsEquationOfStateModel_t<EOSCandidate, std::void_t<typename std::remove_cvref_t<EOSCandidate>::Relations>> :
+        std::bool_constant<
+            relations::ValidRelationCatalog<typename std::remove_cvref_t<EOSCandidate>::Relations> && // First ensure the relation catalog is itself valid
+            ImplementsRelationCatalog_t< // Then ensure that the EOSCandidate implements all relations in the catalog
+                std::remove_cvref_t<EOSCandidate>,
+                typename std::remove_cvref_t<EOSCandidate>::Relations>::value>{ };
+
+    template <typename EOS, typename RelationT, typename InputQuantity>
+    struct ImplementsPartialDerivative_t : std::false_type {};
+
+    template <typename EOS, typename Output, typename... Inputs, typename InputQuantity>
+    struct ImplementsPartialDerivative_t<EOS, relations::Relation<Output, Inputs...>, InputQuantity> :
+        std::bool_constant<
+            (std::same_as<InputQuantity, Inputs> || ...) &&
+            requires (const std::remove_cvref_t<EOS> &eos, dimensions::QuantityValue<Inputs>... inputValues) {
+        {eos.partialDerivative(relations::Relation<Output, Inputs...>{}, dimensions::WithRespectTo<InputQuantity>{}, inputValues...)} -> std::same_as<dimensions::PartialDerivative<Output, Inputs...>>;
+            }>{};
+
+    template <typename EOSCandidate>
+    inline constexpr bool is_equation_of_state_model_v = IsEquationOfStateModel_t<std::remove_cvref_t<EOSCandidate>>::value;
+
+    template <typename EOS, typename RelationT>
+    inline constexpr bool supports_relation_v = is_equation_of_state_model_v<EOS> && ImplementsRelation_t<std::remove_cvref_t<EOS>, std::remove_cvref_t<RelationT>>::value;
+
+    template <typename EOS, typename RelationT, typename InputQuantity>
+    inline constexpr bool supports_partial_derivative_v = is_equation_of_state_model_v<EOS> && ImplementsPartialDerivative_t<std::remove_cvref_t<EOS>, std::remove_cvref_t<RelationT>, InputQuantity>::value;
+}

src/include/serif/eos/relations/catalog.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include "serif/utils/misc/concepts/type_uniqueness.hpp"
+#include "serif/utils/misc/std_helper/variadic_cleaning.hpp"
+#include "serif/dimensions/quantities.hpp"
+
+namespace serif::eos::relations {
+    using utils::misc::std_helper::CountVariadicArguments;
+    using utils::misc::concepts::AllTypesAreUnique;
+
+    template <typename Output, typename... Inputs>
+    struct Relation final {
+        using OutputQuantity = Output;
+        using InputQuantities = dimensions::QuantityList<Inputs...>;
+
+        static constexpr std::size_t inputCount = CountVariadicArguments<Inputs...>();
+    };
+
+    template <typename... Relations>
+    struct RelationCatalog final {
+        static constexpr std::size_t size = CountVariadicArguments<Relations...>();
+    };
+}

src/include/serif/eos/relations/concepts.hpp

@@ -0,0 +1,18 @@
+#pragma once
+#include "serif/eos/relations/traits.hpp"
+
+
+namespace serif::eos::relations {
+    template <typename Catalog, typename RelationT>
+    concept CatalogContainingRelation = relation_catalog_contains_relation_v<Catalog, RelationT>;
+
+    template <typename RelationT, typename Quantity>
+    concept RelationContainingInput = relation_contains_input_v<RelationT, Quantity>;
+
+    template <typename RelationCandidate>
+    concept ValidThermodynamicRelation = is_thermodynamic_relation_v<RelationCandidate>;
+
+    template <typename CandidateCatalog>
+    concept ValidRelationCatalog = is_valid_relation_catalog_v<CandidateCatalog>;
+}
+

src/include/serif/eos/relations/relations.hpp

@@ -0,0 +1,26 @@
+#pragma once
+
+#include "serif/eos/relations/catalog.hpp"
+#include "serif/dimensions/dimensions.hpp"
+
+namespace serif::eos::relations {
+    using PressureFromDensity = Relation<
+        dimensions::Pressure,
+        dimensions::Density>;
+
+    using PressureFromSpecificEnthalpy = Relation<
+        dimensions::Pressure,
+        dimensions::SpecificEnthalpy>;
+
+    using SpecificEnthalpyFromDensity = Relation<
+        dimensions::SpecificEnthalpy,
+        dimensions::Density>;
+
+    using SpecificEnthalpyFromPressure = Relation<
+        dimensions::SpecificEnthalpy,
+        dimensions::Pressure>;
+
+    using DensityFromSpecificEnthalpy = Relation<
+        dimensions::Density,
+        dimensions::SpecificEnthalpy>;
+}

src/include/serif/eos/relations/traits.hpp

@@ -0,0 +1,54 @@
+#pragma once
+
+#include "serif/eos/relations/catalog.hpp"
+
+namespace serif::eos::relations {
+    template <typename Catalog, typename RelationT>
+    struct CatalogContainsRelation : std::false_type {};
+
+    template <typename... Relations, typename RelationT>
+    struct CatalogContainsRelation<RelationCatalog<Relations...>, RelationT> :
+        std::bool_constant<(std::same_as<RelationT, Relations> || ...)> {}; // Fold over all Relations and return true type if any of them match RelationT
+
+    template <typename RelationT, typename Quantity>
+    struct RelationContainsInput : std::false_type {};
+
+    template <typename Output, typename... Inputs, typename Quantity>
+    struct RelationContainsInput<Relation<Output, Inputs...>, Quantity> :
+        std::bool_constant<(std::same_as<Quantity, Inputs> || ...)> {};
+
+    template <typename RelationCandidate>
+    struct IsThermodynamicRelation : std::false_type {};
+
+    template <typename Output, typename... Inputs>
+    struct IsThermodynamicRelation<Relation<Output, Inputs...>> :
+        std::bool_constant<
+            dimensions::ThermodynamicQuantityType<Output> && // Ensure the OutputQuantity is a valid thermodynamic quantity
+            (dimensions::ThermodynamicQuantityType<Inputs> && ...) && // Fold over all InputQuantities to ensure they are valid thermodynamic quantities
+            AllTypesAreUnique<Inputs...>> {}; // Ensure all InputQuantities are unique types
+
+    template <typename CandidateCatalog>
+    struct IsValidRelationCatalog : std::false_type {};
+
+    template <typename... Relations>
+    struct IsValidRelationCatalog<RelationCatalog<Relations...>> :
+        std::bool_constant<
+            (sizeof...(Relations) > 0) && // Ensure the catalog is not empty
+            (IsThermodynamicRelation<Relations>::value && ...) && // Fold over all Relations and ensure they are all valid
+            AllTypesAreUnique<Relations...>> {}; // Ensure all Relations are unique types (e.g. no duplicate relations)
+
+    template <typename Catalog, typename RelationT>
+    inline constexpr bool relation_catalog_contains_relation_v = CatalogContainsRelation<std::remove_cv_t<Catalog>, std::remove_cv_t<RelationT>>::value;
+
+    template <typename RelationT, typename Quantity>
+    inline constexpr bool relation_contains_input_v = RelationContainsInput<std::remove_cv_t<RelationT>, std::remove_cv_t<Quantity>>::value;
+
+    template <typename RelationCandidate>
+    constexpr bool is_thermodynamic_relation_v = IsThermodynamicRelation<RelationCandidate>::value;
+
+    template <typename CatalogCandidate>
+    constexpr bool is_valid_relation_catalog_v = IsValidRelationCatalog<CatalogCandidate>::value;
+
+
+
+}

src/include/serif/utils/misc/concepts/exceptions.hpp

@@ -0,0 +1,15 @@
+#pragma once
+#include <exception>
+#include <type_traits>
+#include <string>
+
+#include "serif/utils/misc/concepts/numeric.hpp"
+
+namespace serif::utils::misc::concepts {
+    template <typename ExceptionCandidate>
+    concept IsException = std::is_base_of_v<std::exception, ExceptionCandidate>;
+
+    template <typename ExceptionCandidate, typename ErrorCode>
+    concept IsExceptionWithErrorCode = IsException<ExceptionCandidate> &&
+                                       std::constructible_from<ExceptionCandidate, ErrorCode, std::string>;
+}

src/include/serif/utils/misc/concepts/numeric.hpp

@@ -0,0 +1,26 @@
+#pragma once
+
+#include <concepts>
+#include <type_traits>
+
+namespace serif::utils::misc::concepts {
+    template <typename T>
+    concept IsNumeric = std::integral<T> || std::floating_point<T>;
+
+    template <typename T>
+    concept IsOrderable = std::totally_ordered<T>;
+
+    template <typename T>
+    concept IsComparableToNumeric = std::totally_ordered_with<T, float> ||
+        std::totally_ordered_with<T, double> ||
+        std::totally_ordered_with<T, long double> ||
+        std::totally_ordered_with<T, int> ||
+        std::totally_ordered_with<T, long> ||
+        std::totally_ordered_with<T, long long>;
+
+    template <typename...Ts>
+    concept SetOfOneOrTwo = (sizeof...(Ts) == 1) || (sizeof...(Ts) == 2);
+
+    template <typename EnumCandidate>
+    concept IsEnum = std::is_enum_v<EnumCandidate>;
+}

src/include/serif/utils/misc/concepts/type_uniqueness.hpp

@@ -0,0 +1,8 @@
+#pragma once
+
+#include "serif/utils/misc/traits/type_uniqueness.hpp"
+
+namespace serif::utils::misc::concepts {
+    template <typename...Ts>
+    concept AllTypesAreUnique = traits::all_unique_types_v<Ts...>;
+}

src/include/serif/utils/misc/finite.hpp

@@ -0,0 +1,44 @@
+#pragma once
+#include <cmath>
+#include <stdexcept>
+#include <format>
+
+#include "serif/utils/misc/concepts/numeric.hpp"
+#include "serif/utils/misc/concepts/exceptions.hpp"
+
+namespace serif::utils::misc {
+    template <typename ExceptionType, concepts::IsNumeric Scalar, concepts::IsEnum ErrorCode>
+    requires concepts::IsExceptionWithErrorCode<ExceptionType, ErrorCode>
+    void validate_finite(const Scalar value, ErrorCode error_code) {
+        if (!std::isfinite(value)) {
+          throw ExceptionType(error_code, std::format("The value must be finite. Instead a value of {} has been provided", value));
+        };
+    }
+
+    template <typename ExceptionType, concepts::IsComparableToNumeric Scalar, concepts::IsEnum ErrorCode>
+    requires concepts::IsExceptionWithErrorCode<ExceptionType, ErrorCode>
+    void validate_greater_than(const Scalar value, const Scalar predicate, ErrorCode error_code) {
+        validate_finite<ExceptionType>(value, error_code);
+        if (value < predicate) {
+            throw ExceptionType(error_code, std::format("The value must be finite and greater than {}. While the value is finite it is also less than {}. A value of {} has been provided", predicate, predicate, value));
+        };
+    }
+
+    template <typename ExceptionType, concepts::IsComparableToNumeric Scalar, concepts::IsEnum ErrorCode>
+    requires concepts::IsExceptionWithErrorCode<ExceptionType, ErrorCode>
+    void validate_less_than(const Scalar value, const Scalar predicate, ErrorCode error_code) {
+        validate_finite<ExceptionType>(value, error_code);
+        if (value > predicate) {
+            throw ExceptionType(error_code, std::format("The value must be finite and less than {}. While the value is finite it is also greater than {}. A value of {} has been provided", predicate, predicate, value));
+        };
+    }
+
+    template <typename ExceptionType, concepts::IsComparableToNumeric Scalar, concepts::IsEnum ErrorCode>
+    requires concepts::IsExceptionWithErrorCode<ExceptionType, ErrorCode>
+    void validate_nonnegativity(const Scalar value, ErrorCode error_code) {
+        validate_finite<ExceptionType>(value, error_code);
+        validate_greater_than<ExceptionType>(value,static_cast<Scalar>(0), error_code);
+    }
+
+
+}

src/include/serif/utils/misc/std_helper/cleaning.hpp

@@ -0,0 +1,17 @@
+#pragma once
+
+#include <type_traits>
+
+namespace serif::utils::misc::std_helper {
+    template< class T >
+    using NeitherConstNorVolatile = std::remove_cv_t<T>;
+
+    template <class T>
+    using NeitherConstNorVolatileReference = std::remove_cvref_t<T>;
+
+    template<class Child, class Parent>
+    constexpr bool AIsBaseClassOfB = std::is_base_of_v<Parent, Child>;
+
+    template <class Ta, class Tb>
+    constexpr bool SameType = std::is_same_v<Ta, Tb>;
+}

src/include/serif/utils/misc/std_helper/variadic_cleaning.hpp

@@ -0,0 +1,9 @@
+#pragma once
+#include <cstddef>
+
+namespace serif::utils::misc::std_helper {
+    template <typename...Ts>
+    [[nodiscard]] consteval std::size_t CountVariadicArguments() noexcept {
+        return sizeof...(Ts);
+    }
+}

src/include/serif/utils/misc/traits/type_uniqueness.hpp

@@ -0,0 +1,19 @@
+#pragma once
+
+#include "serif/utils/misc/std_helper/cleaning.hpp"
+
+namespace serif::utils::misc::traits {
+    template <typename... Ts>
+    constexpr bool all_unique_types_v = true;
+
+
+    // This can be a bit confusing. Whats going on here is that !SameType<First, Rest> && ...
+    // is preforming variadic folding. That is we check if First is the same as all types in
+    // the rest type list. If First is the same as any of those this will eval to true and then be negated to false.
+    // The only way for the first term (!SameType<First, Rest> && ...) to be true is if First is not the same as any of the types in Rest.
+    // Then we recursively call all_unique_types_v<Rest...> to check the rest. This ensures that this is only true if
+    // every single type in the type list is unique.
+    template <typename First, typename... Rest>
+    constexpr bool all_unique_types_v<First, Rest...> = (!std_helper::SameType<First, Rest> && ...) && all_unique_types_v<Rest...>;
+
+}

src/lib/serif/discretization/domain/mesh/topology.cpp

@@ -0,0 +1,118 @@
+#include "serif/discretization/domain/mesh/topology.hpp"
+
+#include <mfem.hpp>
+
+namespace serif::discretization::domain::mesh {
+    const std::vector<int> MeshTopology::s_noElements{};
+
+    namespace {
+        [[nodiscard]] std::size_t as_index(const int id) noexcept {
+            return static_cast<std::size_t>(id);
+        }
+
+        [[nodiscard]] bool in_range(const int id, const std::size_t size) noexcept {
+            return id >= 0 && as_index(id) < size;
+        }
+    }
+
+    MeshTopology::MeshTopology(const mfem::Mesh &mesh)
+        : m_elementDomainIDs(as_index(mesh.GetNE()), 0),
+          m_faceElements(as_index(mesh.GetNumFaces())),
+          m_elementNeighbors(as_index(mesh.GetNE())),
+          m_boundaryElementBoundaryIDs(as_index(mesh.GetNBE()), 0),
+          m_boundaryElementFaceIDs(as_index(mesh.GetNBE()), -1),
+          m_boundaryElementsByFace(as_index(mesh.GetNumFaces())) {
+
+        // Which domain each volume element claims to belong to.
+        for (int elementID = 0; elementID < mesh.GetNE(); ++elementID) {
+            m_elementDomainIDs[as_index(elementID)] = mesh.GetAttribute(elementID);
+        }
+
+        // Face connectivity.
+        for (int faceID = 0; faceID < mesh.GetNumFaces(); ++faceID) {
+            int firstElementID  = -1;
+            int secondElementID = -1;
+
+            mesh.GetFaceElements(faceID, &firstElementID, &secondElementID);
+
+            m_faceElements[as_index(faceID)] = FaceElements{
+                .firstElementID  = firstElementID,
+                .secondElementID = secondElementID
+            };
+
+            if (firstElementID < 0 || secondElementID < 0) {
+                continue;
+            }
+
+            m_elementNeighbors[as_index(firstElementID)].push_back(secondElementID);
+            m_elementNeighbors[as_index(secondElementID)].push_back(firstElementID);
+        }
+
+        // Boundary elements
+        for (int boundaryElementID = 0; boundaryElementID < mesh.GetNBE(); ++boundaryElementID) {
+            m_boundaryElementBoundaryIDs[as_index(boundaryElementID)] = mesh.GetBdrAttribute(boundaryElementID);
+
+            const int faceID = mesh.GetBdrElementFaceIndex(boundaryElementID);
+
+            m_boundaryElementFaceIDs[as_index(boundaryElementID)] = faceID;
+
+            if (in_range(faceID, m_boundaryElementsByFace.size())) {
+                m_boundaryElementsByFace[as_index(faceID)].push_back(boundaryElementID);
+            }
+        }
+    }
+
+    int MeshTopology::element_count() const noexcept {
+        return static_cast<int>(m_elementDomainIDs.size());
+    }
+
+    int MeshTopology::face_count() const noexcept {
+        return static_cast<int>(m_faceElements.size());
+    }
+
+    int MeshTopology::boundary_element_count() const noexcept {
+        return static_cast<int>(m_boundaryElementBoundaryIDs.size());
+    }
+
+    int MeshTopology::element_domain_id(const int elementID) const {
+        if (!in_range(elementID, m_elementDomainIDs.size())) {
+            return 0;
+        }
+        return m_elementDomainIDs[as_index(elementID)];
+    }
+
+    int MeshTopology::boundary_element_boundary_id(const int boundaryElementID) const {
+        if (!in_range(boundaryElementID, m_boundaryElementBoundaryIDs.size())) {
+            return 0;
+        }
+        return m_boundaryElementBoundaryIDs[as_index(boundaryElementID)];
+    }
+
+    FaceElements MeshTopology::face_elements(const int faceID) const {
+        if (!in_range(faceID, m_faceElements.size())) {
+            return FaceElements{};
+        }
+        return m_faceElements[as_index(faceID)];
+    }
+
+    const std::vector<int> &MeshTopology::element_neighbors(const int elementID) const {
+        if (!in_range(elementID, m_elementNeighbors.size())) {
+            return s_noElements;
+        }
+        return m_elementNeighbors[as_index(elementID)];
+    }
+
+    int MeshTopology::boundary_element_face(const int boundaryElementID) const {
+        if (!in_range(boundaryElementID, m_boundaryElementFaceIDs.size())) {
+            return -1;
+        }
+        return m_boundaryElementFaceIDs[as_index(boundaryElementID)];
+    }
+
+    const std::vector<int> &MeshTopology::boundary_elements_on_face(const int faceID) const {
+        if (!in_range(faceID, m_boundaryElementsByFace.size())) {
+            return s_noElements;
+        }
+        return m_boundaryElementsByFace[as_index(faceID)];
+    }
+}

src/lib/serif/discretization/domain/schema/validation/results.cpp

@@ -0,0 +1,52 @@
+#include "serif/discretization/domain/schema/validation/results.hpp"
+
+namespace serif::discretization::domain::schema::validation {
+    bool SchemaRelationValidationResult::valid() const noexcept {
+        return result.valid();
+    }
+
+    SchemaRelationValidationResult::operator bool() const noexcept {
+        return valid();
+    }
+
+    bool SchemaValidationResult::valid() const noexcept {
+        for (const auto &relationResult : relationResults) {
+            if (!relationResult.valid()) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+    SchemaValidationResult::operator bool() const noexcept {
+        return valid();
+    }
+
+    std::size_t SchemaValidationResult::relation_count() const noexcept {
+        return relationResults.size();
+    }
+
+    std::size_t SchemaValidationResult::failed_relation_count() const noexcept {
+        std::size_t failureCount = 0;
+
+        for (const auto &relationResult : relationResults) {
+            if (!relationResult.valid()) {
+                ++failureCount;
+            }
+        }
+        return failureCount;
+    }
+
+    std::size_t SchemaValidationResult::passed_relation_count() const noexcept {
+        return relationResults.size() - failed_relation_count();
+    }
+
+    std::optional<std::size_t> SchemaValidationResult::first_failed_relation_index() const noexcept {
+        for (std::size_t relationIndex = 0; relationIndex < relationResults.size(); ++relationIndex) {
+            if (!relationResults[relationIndex].valid()) {
+                return relationIndex;
+            }
+        }
+        return std::nullopt;
+    }
+}

src/lib/serif/eos/models/polytropic.cpp

@@ -0,0 +1,146 @@
+#include "serif/eos/models/polytropic.hpp"
+#include "serif/utils/misc/finite.hpp"
+#include "serif/eos/exceptions.hpp"
+
+namespace serif::eos::models {
+    using utils::misc::validate_finite;
+    using utils::misc::validate_greater_than;
+    using utils::misc::validate_less_than;
+    using utils::misc::validate_nonnegativity;
+
+    using namespace dimensions;
+    using namespace relations;
+
+    Polytrope::Polytrope(const PolytropeParameters& parameters) : Polytrope(parameters.n, parameters.K) {}
+
+    Polytrope::Polytrope(const double n, const double K):
+    m_polytropic_index(n),
+    m_polytropic_constant(K),
+    m_enthalpy_scale((n + 1.0) * K) {
+        validate_finite<EOSEvaluationError>(m_polytropic_index, EOSEvaluationErrorCode::outside_domain);
+        validate_greater_than<EOSEvaluationError>(m_polytropic_index, 1.0, EOSEvaluationErrorCode::outside_domain);
+        validate_finite<EOSEvaluationError>(m_polytropic_constant, EOSEvaluationErrorCode::outside_domain);
+        validate_greater_than<EOSEvaluationError>(m_polytropic_constant, 0.0, EOSEvaluationErrorCode::outside_domain);
+    }
+
+    [[nodiscard]] double Polytrope::polytropic_index() const noexcept {
+        return m_polytropic_index;
+    }
+
+    [[nodiscard]] double Polytrope::polytropic_constant() const noexcept {
+        return m_polytropic_constant;
+    }
+
+    [[nodiscard]] double Polytrope::enthalpy_scale() const noexcept {
+        return m_enthalpy_scale;
+    }
+
+    [[nodiscard]] PressureValue Polytrope::evaluate(PressureFromDensity, const DensityValue density) const {
+        validate_nonnegativity<EOSEvaluationError>(density.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
+
+        if (density.value() == 0.0) {
+            // P = 0
+            return PressureValue{0.0};
+        }
+
+        // P = K * (ρ)^(1 + 1/n)
+        return PressureValue{
+            m_polytropic_constant * std::pow(density.value(), 1.0 + 1.0 / m_polytropic_index)
+        };
+    }
+
+    [[nodiscard]] SpecificEnthalpyValue Polytrope::evaluate(SpecificEnthalpyFromDensity, const DensityValue density) const {
+        validate_nonnegativity<EOSEvaluationError>(density.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
+        if (density.value() == 0.0) {
+            // h = 0
+            return SpecificEnthalpyValue{0.0};
+        }
+
+        // h = hs * (ρ)^(1/n)
+        return SpecificEnthalpyValue{
+            m_enthalpy_scale * std::pow(density.value(), 1.0 / m_polytropic_index)
+        };
+    }
+
+    [[nodiscard]] DensityValue Polytrope::evaluate(DensityFromSpecificEnthalpy, SpecificEnthalpyValue specific_enthalpy) const {
+        validate_finite<EOSEvaluationError>(specific_enthalpy.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
+
+        if (specific_enthalpy.value() <= 0.0) {
+            // ρ = 0
+            return DensityValue{0.0};
+        }
+
+        // ρ = (h/hs)^n
+        return DensityValue{std::pow(specific_enthalpy.value() / m_enthalpy_scale, m_polytropic_index)};
+    }
+
+    [[nodiscard]] PressureValue Polytrope::evaluate(PressureFromSpecificEnthalpy, SpecificEnthalpyValue specific_enthalpy) const {
+        const DensityValue density = evaluate(DensityFromSpecificEnthalpy{}, specific_enthalpy);
+        if (specific_enthalpy.value() <= 0.0) {
+            // P = 0
+            return PressureValue{0.0};
+        }
+
+        // P = ρ * h / (n + 1)
+        return PressureValue{density.value() * specific_enthalpy.value()/ (m_polytropic_index + 1.0)};
+    }
+
+    [[nodiscard]] SpecificEnthalpyValue Polytrope::evaluate(SpecificEnthalpyFromPressure, PressureValue pressure) const {
+        validate_nonnegativity<EOSEvaluationError>(pressure.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
+        if (pressure.value() == 0.0) {
+            // h = 0
+            return SpecificEnthalpyValue{0.0};
+        }
+
+        const double indexPlusOne = m_polytropic_index + 1.0;
+
+        // h = (n+1) * (K^(n/(n+1))) * (P^(1/(n+1)))
+        return SpecificEnthalpyValue{
+            indexPlusOne * std::pow(m_polytropic_constant, m_polytropic_index / indexPlusOne) *
+            std::pow(pressure.value(), 1.0 / indexPlusOne)
+        };
+    }
+
+    [[nodiscard]] PartialDerivative<Density, SpecificEnthalpy> Polytrope::partial_derivative(DensityFromSpecificEnthalpy, WithRespectTo<SpecificEnthalpy>, const SpecificEnthalpyValue specific_enthalpy) const {
+        validate_finite<EOSEvaluationError>(specific_enthalpy.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
+
+        if (specific_enthalpy.value() < 0.0) {
+            // dρ/dh = 0
+            return PartialDerivative<Density, SpecificEnthalpy>{0.0};
+        }
+
+        if (specific_enthalpy.value() == 0.0) {
+            // dρ/dh = n == 1 -> 1/hs
+            // dρ/dh = else   -> 0
+            return PartialDerivative<Density, SpecificEnthalpy>{
+                m_polytropic_index == 1.0 ? 1.0 / m_enthalpy_scale : 0.0
+            };
+        }
+
+        // dρ/dh = (n/hs) * (h/hs)^(n-1)
+        return PartialDerivative<Density, SpecificEnthalpy>{
+            m_polytropic_index / m_enthalpy_scale * std::pow(specific_enthalpy.value() / m_enthalpy_scale, m_polytropic_index - 1.0)
+        };
+    }
+
+    [[nodiscard]] PartialDerivative<Pressure, SpecificEnthalpy> Polytrope::partial_derivative(PressureFromSpecificEnthalpy, WithRespectTo<SpecificEnthalpy>, const SpecificEnthalpyValue specific_enthalpy) const {
+        const DensityValue density = evaluate(DensityFromSpecificEnthalpy{}, specific_enthalpy);
+
+        // dP/dh = ρ
+        return PartialDerivative<Pressure, SpecificEnthalpy>{density.value()};
+    }
+
+    [[nodiscard]] PartialDerivative<Pressure, Density> Polytrope::partial_derivative(PressureFromDensity, WithRespectTo<Density>, const DensityValue density) const {
+        validate_nonnegativity<EOSEvaluationError>(density.value(), EOSEvaluationErrorCode::invalid_thermodynamic_input);
+
+        if (density.value() == 0.0) {
+            // dP/dρ = 0
+            return PartialDerivative<Pressure, Density>{0.0};
+        }
+
+        // dP/dρ = K * (1 + 1/n) * (ρ)^(1/n)
+        return PartialDerivative<Pressure, Density>{
+            m_polytropic_constant * (1.0 + 1.0 / m_polytropic_index) * std::pow(density.value(), 1.0 / m_polytropic_index)
+        };
+    }
+}

src/meson.build

@@ -0,0 +1,21 @@
+include_dir = include_directories('include')
+
+serif_sources = files(
+  'lib/serif/discretization/domain/mesh/topology.cpp',
+  'lib/serif/discretization/domain/schema/validation/results.cpp',
+  'lib/serif/eos/models/polytropic.cpp',
+)
+
+serif_lib = library(
+  'serif',
+  serif_sources,
+  include_directories: include_dir,
+  dependencies: [mfem_dep],
+  build_rpath: mfem_runtime_prefix == '' ? '' : mfem_runtime_prefix / 'lib',
+)
+
+serif_dep = declare_dependency(
+  include_directories: include_dir,
+  link_with: serif_lib,
+  dependencies: [mfem_dep],
+)

src/python/bindings.cpp

@@ -1,91 +0,0 @@
-#include <meson_mfem_template/config.hpp>
-#include <mfem.hpp>
-#include <nanobind/nanobind.h>
-#include <nanobind/stl/string.h>
-
-#include <cmath>
-#include <stdexcept>
-#include <string>
-
-namespace nb = nanobind;
-
-namespace
-{
-int serial_poisson_dofs(int cells, int order)
-{
-   if (cells < 2 || order < 1)
-   {
-      throw nb::value_error("cells must be >= 2 and order must be >= 1");
-   }
-   mfem::Mesh mesh = mfem::Mesh::MakeCartesian2D(
-      cells, cells, mfem::Element::QUADRILATERAL, true, 1.0, 1.0);
-   mfem::H1_FECollection elements(order, mesh.Dimension());
-   mfem::FiniteElementSpace space(&mesh, &elements);
-   mfem::Array<int> essential_boundary(mesh.bdr_attributes.Max());
-   essential_boundary = 1;
-   mfem::Array<int> essential_dofs;
-   space.GetEssentialTrueDofs(essential_boundary, essential_dofs);
-
-   mfem::ConstantCoefficient one(1.0);
-   mfem::LinearForm rhs(&space);
-   rhs.AddDomainIntegrator(new mfem::DomainLFIntegrator(one));
-   rhs.Assemble();
-   mfem::GridFunction solution(&space);
-   solution = 0.0;
-   mfem::BilinearForm diffusion(&space);
-   diffusion.AddDomainIntegrator(new mfem::DiffusionIntegrator(one));
-   diffusion.Assemble();
-
-   mfem::OperatorPtr matrix;
-   mfem::Vector linear_rhs;
-   mfem::Vector linear_solution;
-   diffusion.FormLinearSystem(
-      essential_dofs, solution, rhs, matrix, linear_solution, linear_rhs);
-   mfem::GSSmoother smoother(static_cast<mfem::SparseMatrix &>(*matrix));
-   mfem::CGSolver solver;
-   solver.SetOperator(*matrix);
-   solver.SetPreconditioner(smoother);
-   solver.SetRelTol(1e-10);
-   solver.SetMaxIter(200);
-   solver.SetPrintLevel(0);
-   solver.Mult(linear_rhs, linear_solution);
-   if (!solver.GetConverged())
-   {
-      throw std::runtime_error("MFEM CG solve did not converge");
-   }
-   return space.GetTrueVSize();
-}
-
-nb::dict capabilities()
-{
-   nb::dict result;
-   result["mpi"] = MESON_MFEM_HAS_MPI != 0;
-   result["hypre"] = MESON_MFEM_HAS_HYPRE != 0;
-   result["metis"] = MESON_MFEM_HAS_METIS != 0;
-   result["cuda"] = MESON_MFEM_HAS_CUDA != 0;
-   result["hip"] = MESON_MFEM_HAS_HIP != 0;
-   result["openmp"] = MESON_MFEM_HAS_OPENMP != 0;
-   result["gslib"] = MESON_MFEM_HAS_GSLIB != 0;
-   result["zlib"] = MESON_MFEM_HAS_ZLIB != 0;
-   result["sundials"] = MESON_MFEM_HAS_SUNDIALS != 0;
-   result["ceed"] = MESON_MFEM_HAS_CEED != 0;
-   result["fms"] = MESON_MFEM_HAS_FMS != 0;
-   result["algoim"] = MESON_MFEM_HAS_ALGOIM != 0;
-   result["simd"] = MESON_MFEM_HAS_SIMD != 0;
-   result["bundled_mfem"] = MESON_MFEM_USING_BUNDLED_MFEM != 0;
-   return result;
-}
-} // namespace
-
-NB_MODULE(_core, module)
-{
-   module.doc() = "Small nanobind validation module for the Meson MFEM bundle";
-   module.def("mfem_version", []() { return std::string(MFEM_VERSION_STRING); });
-   module.def("capabilities", &capabilities);
-   module.def(
-      "serial_poisson_dofs",
-      &serial_poisson_dofs,
-      nb::arg("cells") = 4,
-      nb::arg("order") = 2,
-      "Assemble and solve a small MFEM Poisson problem; return true DOFs.");
-}

tests/discritization/domain/concepts.cpp

@@ -0,0 +1,25 @@
+#include <catch2/catch_test_macros.hpp>
+
+#include "serif/discretization/domain/concepts.hpp"
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+#include "serif/tests/test_tags.hpp"
+
+namespace domain = serif::discretization::domain;
+
+TEST_CASE(
+    "Domain Types And Composite Domains Preserve Their Semantic Categories",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    STATIC_REQUIRE(domain::IsDomain<domain::CoreDomain>);
+    STATIC_REQUIRE(domain::IsDomain<domain::EnvelopeDomain>);
+    STATIC_REQUIRE(domain::IsDomain<domain::VacuumDomain>);
+    STATIC_REQUIRE(domain::IsDomainSet<domain::StellarDomains>);
+    STATIC_REQUIRE(domain::IsDomainSet<domain::AllDomains>);
+    STATIC_REQUIRE_FALSE(domain::IsDomain<domain::StellarDomains>);
+    STATIC_REQUIRE(domain::IsDomainOrSet<domain::StellarDomains>);
+    STATIC_REQUIRE(domain::IsBoundary<domain::StellarSurfaceBoundary>);
+    STATIC_REQUIRE(domain::IsBoundary<domain::InfinitySurfaceBoundary>);
+    CHECK(true);
+}

tests/discritization/domain/ids/lists/lists.cpp

@@ -0,0 +1,79 @@
+#include <catch2/catch_test_macros.hpp>
+
+#include "serif/discretization/domain/ids/boundary.hpp"
+#include "serif/discretization/domain/ids/domain.hpp"
+#include "serif/discretization/domain/ids/lists/lists.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+#include "serif/tests/test_tags.hpp"
+#include "serif/tests/discritization/domain/domain_test_utils.hpp"
+
+namespace domain = serif::discretization::domain;
+namespace ids    = domain::ids;
+
+TEST_CASE(
+    "Material Lists Reject Duplicate Ids And Duplicate Semantic Domains",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    STATIC_REQUIRE(
+        domain_test_utils::CanFormDomainIDList<
+            ids::DomainID<domain::CoreDomain, 1>,
+            ids::DomainID<domain::EnvelopeDomain, 2>>
+    );
+
+    STATIC_REQUIRE_FALSE(
+        domain_test_utils::CanFormDomainIDList<
+            ids::DomainID<domain::CoreDomain, 1>,
+            ids::DomainID<domain::EnvelopeDomain, 1>>
+    );
+
+    STATIC_REQUIRE_FALSE(
+        domain_test_utils::CanFormDomainIDList<
+            ids::DomainID<domain::CoreDomain, 1>,
+            ids::DomainID<domain::CoreDomain, 2>>
+    );
+
+    /*
+     * The schema intentionally imposes no convention on the
+     * numerical range or indexing scheme used by a mesh producer.
+     */
+    STATIC_REQUIRE(
+        domain_test_utils::CanFormDomainIDList<
+            ids::DomainID<domain::CoreDomain, 0>,
+            ids::DomainID<domain::EnvelopeDomain, -7>,
+            ids::DomainID<domain::VacuumDomain, 42>>
+    );
+
+    CHECK(true);
+}
+
+TEST_CASE(
+    "Boundary Lists Reject Duplicate Ids And Duplicate Semantic Boundaries",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    STATIC_REQUIRE(
+        domain_test_utils::CanFormBoundaryIDList<
+            ids::BoundaryID<domain::StellarSurfaceBoundary, 1>,
+            ids::BoundaryID<domain::InfinitySurfaceBoundary, 2>>
+    );
+
+    STATIC_REQUIRE_FALSE(
+        domain_test_utils::CanFormBoundaryIDList<
+            ids::BoundaryID<domain::StellarSurfaceBoundary, 1>,
+            ids::BoundaryID<domain::InfinitySurfaceBoundary, 1>>
+    );
+
+    STATIC_REQUIRE_FALSE(
+        domain_test_utils::CanFormBoundaryIDList<
+            ids::BoundaryID<domain::StellarSurfaceBoundary, 1>,
+            ids::BoundaryID<domain::StellarSurfaceBoundary, 2>>
+    );
+
+    STATIC_REQUIRE(
+        domain_test_utils::CanFormBoundaryIDList<
+            ids::BoundaryID<domain::StellarSurfaceBoundary, 0>,
+            ids::BoundaryID<domain::InfinitySurfaceBoundary, -13>>
+    );
+
+    CHECK(true);
+}

tests/discritization/domain/relation/relations.cpp

@@ -0,0 +1,36 @@
+#include <catch2/catch_test_macros.hpp>
+
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+#include "serif/tests/test_tags.hpp"
+#include "serif/tests/discritization/domain/domain_test_utils.hpp"
+
+namespace domain = serif::discretization::domain;
+
+TEST_CASE(
+    "Domain Boundary Relations Accept Exactly One Or Two Domains",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    STATIC_REQUIRE(
+        domain_test_utils::CanFormDomainBoundary<
+            domain::InfinitySurfaceBoundary, domain::VacuumDomain>
+    );
+
+    STATIC_REQUIRE(
+        domain_test_utils::CanFormDomainBoundary<
+            domain::StellarSurfaceBoundary, domain::StellarDomains,
+            domain::VacuumDomain>
+    );
+
+    STATIC_REQUIRE_FALSE(domain_test_utils::CanFormDomainBoundary<domain::StellarSurfaceBoundary>);
+
+    STATIC_REQUIRE_FALSE(
+        domain_test_utils::CanFormDomainBoundary<
+            domain::StellarSurfaceBoundary, domain::CoreDomain,
+            domain::EnvelopeDomain, domain::VacuumDomain>
+    );
+
+    CHECK(true);
+}

tests/discritization/domain/schema/domain_schema.cpp

@@ -0,0 +1,55 @@
+#include <catch2/catch_test_macros.hpp>
+
+#include "serif/discretization/domain/ids/boundary.hpp"
+#include "serif/discretization/domain/ids/domain.hpp"
+#include "serif/discretization/domain/ids/lists/lists.hpp"
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/relation/lists/relation_list.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/schema/domain_schema.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+#include "serif/tests/test_tags.hpp"
+#include "serif/tests/discritization/domain/domain_test_utils.hpp"
+
+namespace domain   = serif::discretization::domain;
+namespace ids      = domain::ids;
+namespace relation = domain::relation;
+
+TEST_CASE(
+    "Domain Schemas Reject Relations That Reference Unregistered Entities",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    using IncompleteDomainIDs = ids::lists::DomainIDList<
+        ids::DomainID<domain::CoreDomain, 1>,
+        ids::DomainID<domain::VacuumDomain, 3>>;
+
+    using CompleteDomainIDs = ids::lists::DomainIDList<
+        ids::DomainID<domain::CoreDomain, 1>,
+        ids::DomainID<domain::EnvelopeDomain, 2>,
+        ids::DomainID<domain::VacuumDomain, 3>>;
+
+    using CompleteBoundaryIDs = ids::lists::BoundaryIDList<
+        ids::BoundaryID<domain::StellarSurfaceBoundary, 1>,
+        ids::BoundaryID<domain::InfinitySurfaceBoundary, 2>>;
+
+    using InfinityOnlyBoundaryIDs = ids::lists::BoundaryIDList<
+        ids::BoundaryID<domain::InfinitySurfaceBoundary, 2>>;
+
+    using MissingEnvelopeRelation = relation::lists::RelationList<
+        relation::FullyConnected<domain::EnvelopeDomain>>;
+
+    using MissingBoundaryRelation = relation::lists::RelationList<relation::DomainBoundary<
+        domain::StellarSurfaceBoundary, domain::StellarDomains,
+        domain::VacuumDomain>>;
+
+    STATIC_REQUIRE_FALSE(
+        domain_test_utils::CanFormSchema<IncompleteDomainIDs, CompleteBoundaryIDs, MissingEnvelopeRelation>
+    );
+
+    STATIC_REQUIRE_FALSE(
+        domain_test_utils::CanFormSchema<CompleteDomainIDs, InfinityOnlyBoundaryIDs, MissingBoundaryRelation>
+    );
+
+    CHECK(true);
+}

tests/discritization/domain/schema/schemas.cpp

@@ -0,0 +1,64 @@
+#include <string_view>
+
+#include <catch2/catch_test_macros.hpp>
+
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/schema/concepts.hpp"
+#include "serif/discretization/domain/schema/schemas.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+#include "serif/tests/test_tags.hpp"
+
+namespace domain = serif::discretization::domain;
+namespace schema = domain::schema;
+
+TEST_CASE(
+    "Core Envelope Vacuum Schema Exposes Exact Compile Time And Runtime "
+    "Metadata",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    using SchemaT = schema::CoreEnvelopeVacuumDomainSchema;
+
+    STATIC_REQUIRE(schema::IsSchema<SchemaT>);
+    STATIC_REQUIRE(SchemaT::domain_id_count == 3);
+    STATIC_REQUIRE(SchemaT::boundary_id_count == 2);
+    STATIC_REQUIRE(SchemaT::relation_count == 7);
+
+    constexpr auto domainDescriptors   = SchemaT::domain_descriptors();
+    constexpr auto boundaryDescriptors = SchemaT::boundary_descriptors();
+
+    STATIC_REQUIRE(domainDescriptors[0].name == std::string_view{"core"});
+    STATIC_REQUIRE(domainDescriptors[0].ID == 1);
+
+    STATIC_REQUIRE(domainDescriptors[1].name == std::string_view{"envelope"});
+    STATIC_REQUIRE(domainDescriptors[1].ID == 2);
+
+    STATIC_REQUIRE(domainDescriptors[2].name == std::string_view{"vacuum"});
+    STATIC_REQUIRE(domainDescriptors[2].ID == 3);
+
+    STATIC_REQUIRE(boundaryDescriptors[0].name == std::string_view{"stellar_surface"});
+    STATIC_REQUIRE(boundaryDescriptors[0].ID == 1);
+
+    STATIC_REQUIRE(boundaryDescriptors[1].name == std::string_view{"infinity_surface"});
+    STATIC_REQUIRE(boundaryDescriptors[1].ID == 2);
+
+    STATIC_REQUIRE(SchemaT::template contains_domain<domain::CoreDomain>());
+    STATIC_REQUIRE(SchemaT::template contains_domain<domain::StellarDomains>());
+    STATIC_REQUIRE(SchemaT::template contains_domain<domain::AllDomains>());
+    STATIC_REQUIRE(SchemaT::template domain_id_belongs_to<domain::StellarDomains>(1));
+    STATIC_REQUIRE(SchemaT::template domain_id_belongs_to<domain::StellarDomains>(2));
+    STATIC_REQUIRE_FALSE(SchemaT::template domain_id_belongs_to<domain::StellarDomains>(3));
+    STATIC_REQUIRE(SchemaT::template domain_id_belongs_to<domain::AllDomains>(1));
+    STATIC_REQUIRE(SchemaT::template domain_id_belongs_to<domain::AllDomains>(2));
+    STATIC_REQUIRE(SchemaT::template domain_id_belongs_to<domain::AllDomains>(3));
+    STATIC_REQUIRE(SchemaT::template domain_id<domain::CoreDomain>() == 1);
+    STATIC_REQUIRE(SchemaT::template domain_id<domain::EnvelopeDomain>() == 2);
+    STATIC_REQUIRE(SchemaT::template domain_id<domain::VacuumDomain>() == 3);
+
+    STATIC_REQUIRE(SchemaT::template contains_boundary<domain::StellarSurfaceBoundary>());
+    STATIC_REQUIRE(SchemaT::template contains_boundary<domain::InfinitySurfaceBoundary>());
+    STATIC_REQUIRE(SchemaT::template boundary_id<domain::StellarSurfaceBoundary>() == 1);
+    STATIC_REQUIRE(SchemaT::template boundary_id<domain::InfinitySurfaceBoundary>() == 2);
+
+    CHECK(true);
+}

tests/discritization/domain/schema/utils.cpp

@@ -0,0 +1,42 @@
+#include <catch2/catch_test_macros.hpp>
+#include <mfem.hpp>
+
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/schema/schemas.hpp"
+#include "serif/discretization/domain/schema/utils.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+#include "serif/tests/test_tags.hpp"
+#include "serif/tests/discritization/domain/domain_test_utils.hpp"
+
+namespace domain = serif::discretization::domain;
+namespace schema = domain::schema;
+
+TEST_CASE(
+    "Domain Schema Builds Exact MFEM Attribute Markers",
+    tags::domain &tags::utils &tags::unit
+) {
+    using Schema                         = schema::CoreEnvelopeVacuumDomainSchema;
+
+    const mfem::Mesh mesh                = domain_test_utils::make_layered_mesh();
+
+    const mfem::Array<int> stellarMarker = schema::make_attribute_marker<domain::StellarDomains, Schema>(mesh);
+    const mfem::Array<int> vacuumMarker  = schema::make_attribute_marker<domain::VacuumDomain, Schema>(mesh);
+    const mfem::Array<int> allMarker     = schema::make_attribute_marker<domain::AllDomains, Schema>(mesh);
+
+    REQUIRE(stellarMarker.Size() == 3);
+    REQUIRE(vacuumMarker.Size() == 3);
+    REQUIRE(allMarker.Size() == 3);
+
+    CHECK(stellarMarker[0] == 1);
+    CHECK(stellarMarker[1] == 1);
+    CHECK(stellarMarker[2] == 0);
+
+    CHECK(vacuumMarker[0] == 0);
+    CHECK(vacuumMarker[1] == 0);
+    CHECK(vacuumMarker[2] == 1);
+
+    CHECK(allMarker[0] == 1);
+    CHECK(allMarker[1] == 1);
+    CHECK(allMarker[2] == 1);
+}

tests/discritization/domain/schema/validation/all.cpp

@@ -0,0 +1,182 @@
+#include <array>
+#include <cstddef>
+#include <string_view>
+
+#include <catch2/catch_test_macros.hpp>
+#include <mfem.hpp>
+#include <stroid/stroid.h>
+
+#include "serif/discretization/domain/schema/schemas.hpp"
+#include "serif/discretization/domain/schema/validation/all.hpp"
+
+#include "serif/tests/test_tags.hpp"
+#include "serif/tests/discritization/domain/domain_test_utils.hpp"
+
+namespace domain            = serif::discretization::domain;
+namespace schema            = domain::schema;
+namespace schema_validation = schema::validation;
+
+TEST_CASE(
+    "Complete Schema Validation Accepts A Synthetic Core Envelope Vacuum Mesh",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_layered_mesh();
+
+    const auto validation =
+        schema_validation::validate_schema<schema::CoreEnvelopeVacuumDomainSchema>(mesh);
+
+    REQUIRE(validation.valid());
+    REQUIRE(validation.relationResults.size() == 7);
+    CHECK(validation.failed_relation_count() == 0);
+    CHECK(validation.passed_relation_count() == 7);
+    CHECK_FALSE(validation.first_failed_relation_index().has_value());
+
+    constexpr std::array<std::string_view, 7> expectedRelationNames{"fully_connected", "fully_connected",
+                                                                "fully_connected", "inscribed",
+                                                                "inscribed",       "domain_boundary",
+                                                                "domain_boundary"};
+
+    for (std::size_t relationIndex = 0; relationIndex < expectedRelationNames.size(); ++relationIndex) {
+        CHECK(validation.relationResults[relationIndex].relationIndex == relationIndex);
+        CHECK(validation.relationResults[relationIndex].relationName == expectedRelationNames[relationIndex]);
+        CHECK(validation.relationResults[relationIndex].valid());
+    }
+}
+
+TEST_CASE(
+    "Complete Schema Validation Evaluates Every Relation After A Failure",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    /*
+     * All material topology and the outer vacuum boundary are valid.
+     * Only the Stellar/Vacuum boundary tagging is intentionally absent.
+     */
+    const mfem::Mesh mesh = domain_test_utils::make_layered_mesh(false, true);
+    const auto validation =
+        schema_validation::validate_schema<schema::CoreEnvelopeVacuumDomainSchema>(mesh);
+    CHECK_FALSE(validation.valid());
+    REQUIRE(validation.relationResults.size() == 7);
+    CHECK(validation.failed_relation_count() == 1);
+    CHECK(validation.passed_relation_count() == 6);
+    REQUIRE(validation.first_failed_relation_index().has_value());
+    CHECK(*validation.first_failed_relation_index() == 5);
+
+    for (std::size_t relationIndex = 0; relationIndex < 7; ++relationIndex) {
+        CAPTURE(relationIndex);
+
+        if (relationIndex == 5) {
+            CHECK_FALSE(validation.relationResults[relationIndex].valid());
+            CHECK(
+                validation.relationResults[relationIndex].result.failure ==
+                schema_validation::RelationValidationFailure::DomainBoundaryExpectedFaceIsUntagged
+            );
+            continue;
+        }
+
+        CHECK(validation.relationResults[relationIndex].valid());
+    }
+}
+
+TEST_CASE(
+    "STROID Meshes Satisfy The Core Envelope Vacuum Domain Schema",
+    tags::integration &tags::mesh &tags::utils &tags::domain
+) {
+    constexpr std::array<domain_test_utils::StroidCase, 3> testCases{
+        domain_test_utils::StroidCase{
+            .name = "spherical_low_order", .refinementLevels = 0, .order = 1, .flattening = 0.0
+        },
+        domain_test_utils::StroidCase{.name = "oblate", .refinementLevels = 0, .order = 2, .flattening = 0.15},
+        domain_test_utils::StroidCase{.name = "refined_oblate", .refinementLevels = 1, .order = 2, .flattening = 0.10}
+    };
+
+    for (const auto &[name, refinementLevels, order, flattening] : testCases) {
+        INFO("STROID case = " << name);
+        INFO("Refinement levels = " << refinementLevels);
+        INFO("Order = " << order);
+        INFO("Flattening = " << flattening);
+
+        const stroid::config::MeshConfig config = domain_test_utils::make_stroid_config(refinementLevels, order, flattening);
+        stroid::StroidMesh stroidMesh = stroid::GenerateMesh(config);
+
+        REQUIRE(stroidMesh.reference_mesh != nullptr);
+        REQUIRE(stroidMesh.mesh != nullptr);
+
+        /*
+         * Validate both the reference topology and the projected
+         * physical mesh. The mapping/projection must not alter
+         * material or boundary semantics.
+         */
+        domain_test_utils::check_schema_is_valid<schema::CoreEnvelopeVacuumDomainSchema>(
+            *stroidMesh.reference_mesh
+        );
+
+        domain_test_utils::check_schema_is_valid<schema::CoreEnvelopeVacuumDomainSchema>(
+            *stroidMesh.mesh
+        );
+    }
+}
+
+TEST_CASE(
+    "STROID Material And Boundary Id Conventions Are Fully Schema Driven",
+    tags::integration &tags::mesh &tags::utils &tags::domain
+) {
+    stroid::config::MeshConfig config = domain_test_utils::make_stroid_config(0, 1, 0.0);
+
+    config.core_id                    = 11;
+    config.envelope_id                = 17;
+    config.vacuum_id                  = 29;
+    config.surface_bdr_id             = 101;
+    config.inf_bdr_id                 = 203;
+
+    stroid::StroidMesh stroidMesh     = stroid::GenerateMesh(config);
+
+    REQUIRE(stroidMesh.reference_mesh != nullptr);
+    REQUIRE(stroidMesh.mesh != nullptr);
+
+    /*
+     * The same semantic topology must validate when a mesh generator
+     * uses an entirely different attribute numbering convention.
+     */
+    domain_test_utils::check_schema_is_valid<domain_test_utils::AlternateIdSchema>(*stroidMesh.reference_mesh);
+    domain_test_utils::check_schema_is_valid<domain_test_utils::AlternateIdSchema>(*stroidMesh.mesh);
+
+    /*
+     * Conversely, the production 1/2/3 + 1/2 schema must not silently
+     * accept a mesh generated under another numbering convention.
+     */
+    const auto productionValidation =
+        schema_validation::validate_schema<schema::CoreEnvelopeVacuumDomainSchema>(*stroidMesh.mesh);
+
+    CHECK_FALSE(productionValidation.valid());
+    CHECK(productionValidation.failed_relation_count() > 0);
+}
+
+TEST_CASE(
+    "Complete Schema Validation Rejects A Mesh Without Vacuum",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(3, 3, {2, 2, 2, 2, 1, 2, 2, 2, 2}, {});
+
+    const auto validation = schema_validation::validate_schema<schema::CoreEnvelopeVacuumDomainSchema>(mesh);
+
+    CHECK_FALSE(validation.valid());
+    REQUIRE(validation.relationResults.size() == 7);
+
+    /*
+     * FullyConnected<VacuumDomain>
+     */
+    CHECK_FALSE(validation.relationResults[2].valid());
+    CHECK(
+        validation.relationResults[2].result.failure ==
+        schema_validation::RelationValidationFailure::DomainAbsent
+    );
+
+    /*
+     * Inscribed<StellarDomains, VacuumDomain>
+     */
+    CHECK_FALSE(validation.relationResults[4].valid());
+    CHECK(
+        validation.relationResults[4].result.failure ==
+        schema_validation::RelationValidationFailure::OuterDomainAbsent
+    );
+}

tests/discritization/domain/schema/validation/boundary.cpp

@@ -0,0 +1,179 @@
+#include <vector>
+
+#include <catch2/catch_test_macros.hpp>
+#include <mfem.hpp>
+
+#include "serif/discretization/domain/mesh/topology.hpp"
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/relation/validation/runtime.hpp"
+#include "serif/discretization/domain/schema/schemas.hpp"
+#include "serif/discretization/domain/schema/validation/boundary.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+#include "serif/tests/test_tags.hpp"
+#include "serif/tests/discritization/domain/domain_test_utils.hpp"
+
+namespace domain     = serif::discretization::domain;
+namespace relation   = domain::relation;
+namespace schema     = domain::schema;
+namespace validation = schema::validation;
+
+using domain::mesh::MeshTopology;
+
+TEST_CASE(
+    "Domain Boundary Accepts A Complete Internal Stellar Vacuum Interface",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    std::vector<domain_test_utils::BoundaryEdge> boundaries{{.firstVertexId = 1, .secondVertexId = 4, .attribute = 1}};
+
+    const mfem::Mesh mesh     = domain_test_utils::make_grid_mesh(2, 1, {2, 3}, boundaries);
+    const MeshTopology topology{mesh};
+
+    const auto result         = validation::RelationValidator<relation::DomainBoundary<
+        domain::StellarSurfaceBoundary, domain::StellarDomains,
+        domain::VacuumDomain>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(topology);
+
+    CHECK(result);
+
+    /*
+     * Interface ordering is intentionally semantic rather
+     * than oriented.
+     */
+    const auto reversedResult = validation::RelationValidator<relation::DomainBoundary<
+        domain::StellarSurfaceBoundary, domain::VacuumDomain,
+        domain::StellarDomains>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(topology);
+
+    CHECK(reversedResult);
+}
+
+TEST_CASE(
+    "Domain Boundary Accepts A Complete Exterior Vacuum Boundary",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const std::vector<int> attributes{3};
+    std::vector<domain_test_utils::BoundaryEdge> boundaries;
+
+    domain_test_utils::append_exterior_boundaries(
+        boundaries, attributes, 1, 1, [](const int materialId) { return materialId == 3; }, 2
+    );
+
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(1, 1, attributes, boundaries);
+    const auto result     = validation::RelationValidator<relation::DomainBoundary<
+        domain::InfinitySurfaceBoundary, domain::VacuumDomain>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK(result);
+}
+
+TEST_CASE(
+    "Domain Boundary Rejects A Tagged Internal Face For An Exterior Boundary",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(2, 1, {3, 3}, {{.firstVertexId = 1, .secondVertexId = 4, .attribute = 2}});
+
+    const auto result = validation::RelationValidator<relation::DomainBoundary<
+        domain::InfinitySurfaceBoundary, domain::VacuumDomain>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(
+        result.failure == validation::RelationValidationFailure::DomainBoundaryTaggedFaceHasWrongTopology
+    );
+}
+
+TEST_CASE(
+    "Domain Boundary Rejects A Tagged Exterior Face Of The Wrong Material",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(1, 1, {2}, {{.firstVertexId = 0, .secondVertexId = 1, .attribute = 2}});
+
+    const auto result = validation::RelationValidator<relation::DomainBoundary<
+        domain::InfinitySurfaceBoundary, domain::VacuumDomain>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(
+        result.failure ==
+        validation::RelationValidationFailure::DomainBoundaryTaggedFaceTouchesUnexpectedDomain
+    );
+}
+
+TEST_CASE(
+    "Domain Boundary Rejects A Tagged Internal Interface With Unexpected "
+    "Materials",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh =
+        domain_test_utils::make_grid_mesh(2, 1, {1, 2}, {{.firstVertexId = 1, .secondVertexId = 4, .attribute = 1}});
+
+    const auto result = validation::RelationValidator<relation::DomainBoundary<
+        domain::StellarSurfaceBoundary, domain::StellarDomains,
+        domain::VacuumDomain>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+
+    CHECK(
+        result.failure ==
+        validation::RelationValidationFailure::DomainBoundaryTaggedFaceTouchesUnexpectedDomain
+    );
+}
+
+TEST_CASE(
+    "Domain Boundary Rejects An Untagged Expected Interface",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(2, 1, {2, 3}, {});
+
+    const auto result     = validation::RelationValidator<relation::DomainBoundary<
+        domain::StellarSurfaceBoundary, domain::StellarDomains,
+        domain::VacuumDomain>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(result.failure == validation::RelationValidationFailure::DomainBoundaryExpectedFaceIsUntagged);
+    REQUIRE(result.domainBoundaryDiagnostics.has_value());
+    CHECK(result.domainBoundaryDiagnostics->faceID >= 0);
+    CHECK(result.domainBoundaryDiagnostics->boundaryElementID == -1);
+    CHECK_FALSE(result.domainBoundaryDiagnostics->actualBoundaryID.has_value());
+}
+
+TEST_CASE(
+    "Domain Boundary Rejects An Expected Interface With The Wrong Attribute",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(2, 1, {2, 3}, {{.firstVertexId = 1, .secondVertexId = 4, .attribute = 9}});
+
+    const auto result = validation::RelationValidator<relation::DomainBoundary<
+        domain::StellarSurfaceBoundary, domain::StellarDomains,
+        domain::VacuumDomain>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(
+        result.failure ==
+        validation::RelationValidationFailure::DomainBoundaryExpectedFaceHasWrongID
+    );
+    REQUIRE(result.domainBoundaryDiagnostics.has_value());
+    REQUIRE(result.domainBoundaryDiagnostics->actualBoundaryID.has_value());
+    CHECK(*result.domainBoundaryDiagnostics->actualBoundaryID == 9);
+    CHECK(result.domainBoundaryDiagnostics->expectedBoundaryID == 1);
+}
+
+TEST_CASE(
+    "Domain Boundary Rejects A Relation That Is Not Realized Anywhere",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(2, 1, {2, 2}, {});
+
+    const auto result     = validation::RelationValidator<relation::DomainBoundary<
+        domain::StellarSurfaceBoundary, domain::StellarDomains,
+        domain::VacuumDomain>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(result.failure == validation::RelationValidationFailure::DomainBoundaryAbsent);
+}

tests/discritization/domain/schema/validation/connected.cpp

@@ -0,0 +1,73 @@
+#include <catch2/catch_test_macros.hpp>
+#include <mfem.hpp>
+
+#include "serif/discretization/domain/mesh/topology.hpp"
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/relation/validation/runtime.hpp"
+#include "serif/discretization/domain/schema/schemas.hpp"
+#include "serif/discretization/domain/schema/validation/connected.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+#include "serif/tests/test_tags.hpp"
+#include "serif/tests/discritization/domain/domain_test_utils.hpp"
+
+namespace domain     = serif::discretization::domain;
+namespace relation   = domain::relation;
+namespace schema     = domain::schema;
+namespace validation = schema::validation;
+
+using domain::mesh::MeshTopology;
+
+TEST_CASE(
+    "Connected Accepts Face Connected Atomic And Composite Domains",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh    = domain_test_utils::make_layered_mesh();
+    const MeshTopology topology{mesh};
+
+    const auto coreResult = validation::RelationValidator<relation::FullyConnected<domain::CoreDomain>>::template validate<schema::CoreEnvelopeVacuumDomainSchema>(topology);
+
+    REQUIRE(coreResult);
+    REQUIRE(coreResult.connectedDiagnostics.has_value());
+    CHECK(coreResult.connectedDiagnostics->domainElementCount == 1);
+    CHECK(coreResult.connectedDiagnostics->visitedElementCount == 1);
+
+    const auto stellarResult = validation::RelationValidator<relation::FullyConnected<domain::StellarDomains>>::template validate<schema::CoreEnvelopeVacuumDomainSchema>(topology);
+
+    REQUIRE(stellarResult);
+    REQUIRE(stellarResult.connectedDiagnostics.has_value());
+    CHECK(stellarResult.connectedDiagnostics->domainElementCount == 9);
+    CHECK(stellarResult.connectedDiagnostics->visitedElementCount == 9);
+}
+
+TEST_CASE(
+    "Connected Rejects An Absent Domain",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(2, 1, {2, 2}, {});
+
+    const auto result = validation::RelationValidator<relation::FullyConnected<domain::CoreDomain>>::template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(result.failure == validation::RelationValidationFailure::DomainAbsent);
+    REQUIRE(result.connectedDiagnostics.has_value());
+    CHECK(result.connectedDiagnostics->domainElementCount == 0);
+    CHECK(result.connectedDiagnostics->visitedElementCount == 0);
+}
+
+TEST_CASE(
+    "Connected Rejects Multiple Face Disconnected Components",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(3, 1, {1, 2, 1}, {});
+
+    const auto result = validation::RelationValidator<relation::FullyConnected<domain::CoreDomain>>::template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(result.failure == validation::RelationValidationFailure::DomainDisconnected);
+    REQUIRE(result.connectedDiagnostics.has_value());
+    CHECK(result.connectedDiagnostics->domainElementCount == 2);
+    CHECK(result.connectedDiagnostics->visitedElementCount == 1);
+    CHECK(result.connectedDiagnostics->elementID >= 0);
+}

tests/discritization/domain/schema/validation/inscribed.cpp

@@ -0,0 +1,96 @@
+#include <vector>
+
+#include <catch2/catch_test_macros.hpp>
+#include <mfem.hpp>
+
+#include "serif/discretization/domain/mesh/topology.hpp"
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/relation/validation/runtime.hpp"
+#include "serif/discretization/domain/schema/schemas.hpp"
+#include "serif/discretization/domain/schema/validation/inscribed.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+#include "serif/tests/test_tags.hpp"
+#include "serif/tests/discritization/domain/domain_test_utils.hpp"
+
+namespace domain     = serif::discretization::domain;
+namespace relation   = domain::relation;
+namespace schema     = domain::schema;
+namespace validation = schema::validation;
+
+using domain::mesh::MeshTopology;
+
+TEST_CASE(
+    "Inscribed Accepts Nested Atomic And Composite Domains",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh    = domain_test_utils::make_layered_mesh();
+    const MeshTopology topology{mesh};
+
+    const auto coreResult = validation::RelationValidator<relation::Inscribed<domain::CoreDomain, domain::EnvelopeDomain>>::template validate<schema::CoreEnvelopeVacuumDomainSchema>(topology);
+
+    CHECK(coreResult);
+
+    const auto stellarResult = validation::RelationValidator<relation::Inscribed<domain::StellarDomains, domain::VacuumDomain>>::template validate<schema::CoreEnvelopeVacuumDomainSchema>(topology);
+
+    CHECK(stellarResult);
+}
+
+TEST_CASE(
+    "Inscribed Rejects An Absent Inner Domain",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(2, 2, {2, 2, 2, 2}, {});
+
+    const auto result = validation::RelationValidator<relation::Inscribed<domain::CoreDomain, domain::EnvelopeDomain>>::template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(result.failure == validation::RelationValidationFailure::InnerDomainAbsent);
+}
+
+TEST_CASE(
+    "Inscribed Rejects An Absent Outer Domain",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(1, 1, {1}, {});
+
+    const auto result = validation::RelationValidator<relation::Inscribed<domain::CoreDomain, domain::EnvelopeDomain>>::template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(result.failure == validation::RelationValidationFailure::OuterDomainAbsent);
+}
+
+TEST_CASE(
+    "Inscribed Rejects An Inner Domain Touching The Computational Boundary",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(2, 2, {1, 2, 2, 2}, {});
+
+    const auto result = validation::RelationValidator<relation::Inscribed<domain::CoreDomain, domain::EnvelopeDomain>>::template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(result.failure == validation::RelationValidationFailure::InnerDomainTouchesMeshBoundary);
+    REQUIRE(result.inscribedDiagnostics.has_value());
+    CHECK(result.inscribedDiagnostics->faceID >= 0);
+    CHECK(result.inscribedDiagnostics->innerElementID >= 0);
+    CHECK(result.inscribedDiagnostics->adjacentElementID == -1);
+}
+
+TEST_CASE(
+    "Inscribed Rejects An Inner Domain Touching An Unexpected Material",
+    tags::unit &tags::mesh &tags::utils &tags::domain
+) {
+    std::vector<int> attributes{2, 2, 2, 2, 1, 3, 2, 2, 2};
+
+    const mfem::Mesh mesh = domain_test_utils::make_grid_mesh(3, 3, attributes, {});
+
+    const auto result     = validation::RelationValidator<
+        relation::Inscribed<domain::CoreDomain, domain::EnvelopeDomain>>::
+        template validate<schema::CoreEnvelopeVacuumDomainSchema>(MeshTopology{mesh});
+
+    CHECK_FALSE(result);
+    CHECK(result.failure == validation::RelationValidationFailure::InnerDomainTouchesUnexpectedDomain);
+    REQUIRE(result.inscribedDiagnostics.has_value());
+    CHECK(result.inscribedDiagnostics->adjacentDomainID == 3);
+}

tests/discritization/meson.build

@@ -0,0 +1,12 @@
+discretization_test_sources = files(
+  'domain/concepts.cpp',
+#  'domain/ids/lists/lists.cpp',
+#  'domain/relation/relations.cpp',
+#  'domain/schema/domain_schema.cpp',
+#  'domain/schema/schemas.cpp',
+#  'domain/schema/utils.cpp',
+#  'domain/schema/validation/all.cpp',
+#  'domain/schema/validation/boundary.cpp',
+#  'domain/schema/validation/connected.cpp',
+#  'domain/schema/validation/inscribed.cpp',
+)

tests/include/serif/tests/discritization/domain/domain_test_utils.hpp

@@ -0,0 +1,339 @@
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <string_view>
+#include <vector>
+
+#include <catch2/catch_test_macros.hpp>
+#include <mfem.hpp>
+#include <stroid/stroid.h>
+
+#include "serif/discretization/domain/ids/boundary.hpp"
+#include "serif/discretization/domain/ids/domain.hpp"
+#include "serif/discretization/domain/ids/lists/lists.hpp"
+#include "serif/discretization/domain/physical_domains.hpp"
+#include "serif/discretization/domain/relation/lists/relation_list.hpp"
+#include "serif/discretization/domain/relation/relations.hpp"
+#include "serif/discretization/domain/schema/domain_schema.hpp"
+#include "serif/discretization/domain/schema/validation/all.hpp"
+#include "serif/discretization/domain/types.hpp"
+
+namespace domain_test_utils {
+    namespace domain     = serif::discretization::domain;
+    namespace ids        = domain::ids;
+    namespace relation   = domain::relation;
+    namespace schema     = domain::schema;
+    namespace validation = schema::validation;
+
+    struct UnregisteredDomain final : public domain::Domain {
+        static constexpr std::string_view name = "unregistered_domain";
+    };
+
+    struct UnregisteredBoundary final : public domain::Boundary {
+        static constexpr std::string_view name = "unregistered_boundary";
+    };
+
+    struct BoundaryEdge {
+        int firstVertexId{-1};
+        int secondVertexId{-1};
+        int attribute{0};
+    };
+
+    struct StroidCase {
+        std::string_view name;
+        int refinementLevels{0};
+        int order{1};
+        double flattening{0.0};
+    };
+
+    template <typename... DomainIDTs>
+    concept CanFormDomainIDList = requires { typename ids::lists::DomainIDList<DomainIDTs...>; };
+
+    template <typename... BoundaryIDTs>
+    concept CanFormBoundaryIDList = requires { typename ids::lists::BoundaryIDList<BoundaryIDTs...>; };
+
+    template <typename BoundaryT, typename... DomainTs>
+    concept CanFormDomainBoundary = requires { typename relation::DomainBoundary<BoundaryT, DomainTs...>; };
+
+    template <typename DomainIDsT, typename BoundaryIDsT, typename RelationsT>
+    concept CanFormSchema = requires { typename schema::DomainSchema<DomainIDsT, BoundaryIDsT, RelationsT>; };
+
+    [[nodiscard]] inline int vertex_id(const int xElementCount, const int x, const int y) {
+        return y * (xElementCount + 1) + x;
+    }
+
+    [[nodiscard]] inline int cell_index(const int xElementCount, const int x, const int y) {
+        return y * xElementCount + x;
+    }
+
+    [[nodiscard]]inline int cell_attribute(const std::vector<int> &attributes, const int xElementCount, const int x, const int y) {
+        return attributes.at(static_cast<std::size_t>(cell_index(xElementCount, x, y)));
+    }
+
+    template <typename FirstPredicateT, typename SecondPredicateT>
+    void append_interface_boundaries(
+        std::vector<BoundaryEdge> &boundaries,
+        const std::vector<int> &attributes,
+        const int xElementCount,
+        const int yElementCount,
+        FirstPredicateT firstPredicate,
+        SecondPredicateT secondPredicate,
+        const int boundaryAttribute
+    ) {
+        /*
+         * Vertical internal faces.
+         */
+        for (int y = 0; y < yElementCount; ++y) {
+            for (int x = 1; x < xElementCount; ++x) {
+                const int leftAttribute  = cell_attribute(attributes, xElementCount, x - 1, y);
+                const int rightAttribute = cell_attribute(attributes, xElementCount, x, y);
+                const bool matches       = (firstPredicate(leftAttribute) && secondPredicate(rightAttribute)) ||
+                                           (secondPredicate(leftAttribute) && firstPredicate(rightAttribute));
+
+                if (!matches) {
+                    continue;
+                }
+
+                boundaries.push_back(
+                    {.firstVertexId  = vertex_id(xElementCount, x, y),
+                     .secondVertexId = vertex_id(xElementCount, x, y + 1),
+                     .attribute      = boundaryAttribute}
+                );
+            }
+        }
+
+        /*
+         * Horizontal internal faces.
+         */
+        for (int y = 1; y < yElementCount; ++y) {
+            for (int x = 0; x < xElementCount; ++x) {
+                const int lowerAttribute = cell_attribute(attributes, xElementCount, x, y - 1);
+                const int upperAttribute = cell_attribute(attributes, xElementCount, x, y);
+                const bool matches       = (firstPredicate(lowerAttribute) && secondPredicate(upperAttribute)) ||
+                                           (secondPredicate(lowerAttribute) && firstPredicate(upperAttribute));
+
+                if (!matches) {
+                    continue;
+                }
+
+                boundaries.push_back(
+                    {.firstVertexId  = vertex_id(xElementCount, x, y),
+                     .secondVertexId = vertex_id(xElementCount, x + 1, y),
+                     .attribute      = boundaryAttribute}
+                );
+            }
+        }
+    }
+
+    template <typename PredicateT>
+    void append_exterior_boundaries(
+        std::vector<BoundaryEdge> &boundaries,
+        const std::vector<int> &attributes,
+        const int xElementCount,
+        const int yElementCount,
+        PredicateT predicate,
+        const int boundaryAttribute
+    ) {
+        /*
+         * Bottom.
+         */
+        for (int x = 0; x < xElementCount; ++x) {
+            if (predicate(cell_attribute(attributes, xElementCount, x, 0))) {
+                boundaries.push_back(
+                    {.firstVertexId  = vertex_id(xElementCount, x, 0),
+                     .secondVertexId = vertex_id(xElementCount, x + 1, 0),
+                     .attribute      = boundaryAttribute}
+                );
+            }
+        }
+
+        /*
+         * Top.
+         */
+        for (int x = 0; x < xElementCount; ++x) {
+            if (predicate(cell_attribute(attributes, xElementCount, x, yElementCount - 1))) {
+                boundaries.push_back(
+                    {.firstVertexId  = vertex_id(xElementCount, x, yElementCount),
+                     .secondVertexId = vertex_id(xElementCount, x + 1, yElementCount),
+                     .attribute      = boundaryAttribute}
+                );
+            }
+        }
+
+        /*
+         * Left.
+         */
+        for (int y = 0; y < yElementCount; ++y) {
+            if (predicate(cell_attribute(attributes, xElementCount, 0, y))) {
+                boundaries.push_back(
+                    {.firstVertexId  = vertex_id(xElementCount, 0, y),
+                     .secondVertexId = vertex_id(xElementCount, 0, y + 1),
+                     .attribute      = boundaryAttribute}
+                );
+            }
+        }
+
+        /*
+         * Right.
+         */
+        for (int y = 0; y < yElementCount; ++y) {
+            if (predicate(cell_attribute(attributes, xElementCount, xElementCount - 1, y))) {
+                boundaries.push_back(
+                    {.firstVertexId  = vertex_id(xElementCount, xElementCount, y),
+                     .secondVertexId = vertex_id(xElementCount, xElementCount, y + 1),
+                     .attribute      = boundaryAttribute}
+                );
+            }
+        }
+    }
+
+    [[nodiscard]] inline mfem::Mesh make_grid_mesh(
+        const int xElementCount,
+        const int yElementCount,
+        const std::vector<int> &attributes,
+        const std::vector<BoundaryEdge> &boundaryEdges
+    ) {
+        REQUIRE(static_cast<int>(attributes.size()) == xElementCount * yElementCount);
+
+        mfem::Mesh mesh(
+            2, (xElementCount + 1) * (yElementCount + 1), xElementCount * yElementCount,
+            static_cast<int>(boundaryEdges.size()), 2
+        );
+
+        for (int y = 0; y <= yElementCount; ++y) {
+            for (int x = 0; x <= xElementCount; ++x) {
+                mesh.AddVertex(static_cast<double>(x), static_cast<double>(y));
+            }
+        }
+
+        for (int y = 0; y < yElementCount; ++y) {
+            for (int x = 0; x < xElementCount; ++x) {
+                const int lowerLeft  = vertex_id(xElementCount, x, y);
+                const int lowerRight = vertex_id(xElementCount, x + 1, y);
+                const int upperRight = vertex_id(xElementCount, x + 1, y + 1);
+                const int upperLeft  = vertex_id(xElementCount, x, y + 1);
+
+                mesh.AddQuad(
+                    lowerLeft, lowerRight, upperRight, upperLeft, cell_attribute(attributes, xElementCount, x, y)
+                );
+            }
+        }
+
+        for (const BoundaryEdge &boundary : boundaryEdges) {
+            mesh.AddBdrSegment(boundary.firstVertexId, boundary.secondVertexId, boundary.attribute);
+        }
+
+        mesh.FinalizeTopology(false);
+        mesh.Finalize(false, false);
+
+        REQUIRE(mesh.GetNBE() == static_cast<int>(boundaryEdges.size()));
+
+        return mesh;
+    }
+
+    [[nodiscard]] inline std::vector<int> make_layered_attributes() {
+        constexpr int xElementCount = 5;
+        constexpr int yElementCount = 5;
+
+        std::vector<int> attributes(xElementCount * yElementCount, 3);
+
+        for (int y = 1; y <= 3; ++y) {
+            for (int x = 1; x <= 3; ++x) {
+                attributes[static_cast<std::size_t>(cell_index(xElementCount, x, y))] = 2;
+            }
+        }
+
+        attributes[static_cast<std::size_t>(cell_index(xElementCount, 2, 2))] = 1;
+
+        return attributes;
+    }
+
+    [[nodiscard]] inline mfem::Mesh make_layered_mesh(
+        const bool includeStellarSurface   = true,
+        const bool includeInfinitySurface  = true,
+        const int stellarSurfaceAttribute  = 1,
+        const int infinitySurfaceAttribute = 2
+    ) {
+        constexpr int xElementCount       = 5;
+        constexpr int yElementCount       = 5;
+
+        const std::vector<int> attributes = make_layered_attributes();
+        std::vector<BoundaryEdge> boundaries;
+
+        const auto isStellar = [](const int materialId) { return materialId == 1 || materialId == 2; };
+        const auto isVacuum  = [](const int materialId) { return materialId == 3; };
+
+        if (includeStellarSurface) {
+            append_interface_boundaries(
+                boundaries, attributes, xElementCount, yElementCount, isStellar, isVacuum, stellarSurfaceAttribute
+            );
+        }
+
+        if (includeInfinitySurface) {
+            append_exterior_boundaries(
+                boundaries, attributes, xElementCount, yElementCount, isVacuum, infinitySurfaceAttribute
+            );
+        }
+
+        return make_grid_mesh(xElementCount, yElementCount, attributes, boundaries);
+    }
+
+    template <typename SchemaT>
+    void check_schema_is_valid(const mfem::Mesh &mesh) {
+        const auto validation = domain::schema::validation::validate_schema<SchemaT>(mesh);
+
+        CHECK(validation.relationResults.size() == SchemaT::relation_count);
+
+        for (const auto &relationResult : validation.relationResults) {
+            INFO("Relation index = " << relationResult.relationIndex);
+            INFO("Relation name = " << relationResult.relationName);
+            INFO("Failure enum = " << static_cast<int>(relationResult.result.failure));
+            CHECK(relationResult.valid());
+        }
+
+        CHECK(validation.valid());
+    }
+
+    using AlternateIdSchema = schema::DomainSchema<
+        ids::lists::DomainIDList<
+            ids::DomainID<domain::CoreDomain, 11>,
+            ids::DomainID<domain::EnvelopeDomain, 17>,
+            ids::DomainID<domain::VacuumDomain, 29>>,
+        ids::lists::BoundaryIDList<
+            ids::BoundaryID<domain::StellarSurfaceBoundary, 101>,
+            ids::BoundaryID<domain::InfinitySurfaceBoundary, 203>>,
+        relation::lists::RelationList<
+            relation::FullyConnected<domain::CoreDomain>,
+            relation::FullyConnected<domain::EnvelopeDomain>,
+            relation::FullyConnected<domain::VacuumDomain>,
+            relation::Inscribed<domain::CoreDomain, domain::EnvelopeDomain>,
+            relation::Inscribed<domain::StellarDomains, domain::VacuumDomain>,
+            relation::DomainBoundary<
+                domain::StellarSurfaceBoundary,
+                domain::StellarDomains,
+                domain::VacuumDomain>,
+            relation::DomainBoundary<domain::InfinitySurfaceBoundary, domain::VacuumDomain>>>;
+
+    [[nodiscard]] inline stroid::config::MeshConfig make_stroid_config(
+        const int refinementLevels,
+        const int order,
+        const double flattening
+    ) {
+        stroid::config::MeshConfig config;
+        config.refinement_levels       = refinementLevels;
+        config.order                   = order;
+        config.include_external_domain = true;
+        config.r_core                  = 0.25;
+        config.r_star                  = 1.0;
+        config.r_infinity              = 4.0;
+        config.flattening              = flattening;
+        config.core_id                 = 1;
+        config.envelope_id             = 2;
+        config.vacuum_id               = 3;
+        config.surface_bdr_id          = 1;
+        config.inf_bdr_id              = 2;
+        config.optimization_methods    = stroid::config::OptimizationMethods{.tmop = false, .smoothstep = true};
+        return config;
+    }
+} // namespace domain_test_utils

tests/include/serif/tests/test_tags.hpp

@@ -0,0 +1,71 @@
+#pragma once
+
+/*
+ * Header form of the tag machinery that currently lives in
+ * test_helpers.cppm. Only the tags used by the discretization/domain
+ * tests are carried over; once test_helpers is itself de-moduled this
+ * file should be deleted and the project wide header included instead.
+ */
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+
+#include <catch2/internal/catch_stringref.hpp>
+
+template <std::size_t N> struct Tag {
+    std::array<char, N> chars{};
+
+    // ReSharper disable once CppNonExplicitConvertingConstructor
+    consteval Tag(
+        std::array<
+            char,
+            N> arr
+    )
+        : chars(arr) {
+    }
+
+    // ReSharper disable once CppNonExplicitConversionOperator
+    constexpr operator const char *() const {
+        return chars.data();
+    }
+
+    // ReSharper disable once CppNonExplicitConversionOperator
+    constexpr operator Catch::StringRef() const {
+        return Catch::StringRef(chars.data(), N - 1);
+    }
+
+    template <std::size_t M> consteval Tag<N + M - 1> operator&(const Tag<M> &other) const {
+        std::array<char, N + M - 1> res{};
+        std::ranges::copy(chars.begin(), chars.end() - 1, res.begin());
+        std::ranges::copy(other.chars, res.begin() + (N - 1));
+        return {res};
+    }
+};
+
+template <std::size_t N> consteval auto make_tag(const char (&str)[N]) {
+    std::array<char, N + 2> res{};
+    res[0] = '[';
+    std::ranges::copy(str, str + N - 1, res.begin() + 1);
+    res[N]     = ']';
+    res[N + 1] = '\0';
+    return Tag<N + 2>{res};
+}
+
+template <
+    std::size_t N,
+    std::size_t M>
+consteval auto sub_tag(
+    const Tag<N> &parent,
+    const char (&str)[M]
+) {
+    return parent & make_tag(str);
+}
+
+namespace tags {
+    inline constexpr auto unit        = make_tag("unit");
+    inline constexpr auto mesh        = make_tag("mesh");
+    inline constexpr auto integration = make_tag("integration");
+    inline constexpr auto utils       = make_tag("utils");
+    inline constexpr auto domain      = sub_tag(mesh, "domain");
+} // namespace tags

tests/meson.build

@@ -1,169 +1,13 @@
-if get_option('build_tests') and get_option('build_examples')
-  runtime_environment = environment()
-  if (
-    mfem_runtime_prefix != '' and dependency_prefix != '' and
-    not wheel_carries_native_bundle
-  )
-    runtime_environment.prepend('PATH', dependency_prefix / 'bin')
-    if host_machine.system() == 'darwin'
-      runtime_environment.prepend('DYLD_LIBRARY_PATH', dependency_prefix / 'lib64')
-      runtime_environment.prepend('DYLD_LIBRARY_PATH', dependency_prefix / 'lib')
-    elif host_machine.system() != 'windows' and not is_wasm
-      runtime_environment.prepend('LD_LIBRARY_PATH', dependency_prefix / 'lib64')
-      runtime_environment.prepend('LD_LIBRARY_PATH', dependency_prefix / 'lib')
-    endif
-  endif
-  if mfem_runtime_prefix != ''
-    runtime_environment.prepend('PATH', mfem_runtime_prefix / 'bin')
-    if host_machine.system() == 'darwin'
-      runtime_environment.prepend('DYLD_LIBRARY_PATH', mfem_runtime_prefix / 'lib64')
-      runtime_environment.prepend('DYLD_LIBRARY_PATH', mfem_runtime_prefix / 'lib')
-    elif host_machine.system() != 'windows' and not is_wasm
-      runtime_environment.prepend('LD_LIBRARY_PATH', mfem_runtime_prefix / 'lib64')
-      runtime_environment.prepend('LD_LIBRARY_PATH', mfem_runtime_prefix / 'lib')
-    endif
-  endif
-
-  test(
-    'serial-poisson',
-    serial_example,
-    env: runtime_environment,
-    timeout: 120,
-  )
-
-  if mfem_has_cuda
-    test(
-      'serial-poisson-cuda',
-      serial_example,
-      args: ['cuda'],
-      env: runtime_environment,
-      timeout: 120,
-    )
-  endif
-
-  if mfem_has_mpi
-    if mpi_launcher_from_dependency and not wheel_carries_native_bundle
-      mpi_launcher_program = find_program(
-        dependency_prefix / 'bin' / 'mpiexec',
-        dependency_prefix / 'bin' / 'mpirun',
-        required: true,
-      )
-      mpi_launcher = mpi_launcher_program.full_path()
-    elif mfem_runtime_prefix != ''
-      mpi_launcher = mfem_runtime_prefix / 'bin' / 'mpiexec'
-    else
-      mpi_launcher_program = find_program('mpiexec', 'mpirun', required: true)
-      mpi_launcher = mpi_launcher_program.full_path()
-    endif
-    test(
-      'parallel-hypre-boomeramg',
-      python_build,
-      args: [
-        files('../tools/run_mpi_test.py'),
-        '--launcher', mpi_launcher,
-        '--processes', '2',
-        parallel_example,
-      ],
-      env: runtime_environment,
-      timeout: 180,
-    )
-    if mfem_has_cuda
-      test(
-        'parallel-hypre-boomeramg-cuda',
-        python_build,
-        args: [
-          files('../tools/run_mpi_test.py'),
-          '--launcher', mpi_launcher,
-          '--processes', '2',
-          parallel_example,
-          'cuda',
-        ],
-        env: runtime_environment,
-        timeout: 180,
-      )
-    endif
-  endif
-endif
-
-if get_option('build_tests') and get_option('build_benchmarks')
-  benchmark_smoke_environment = environment()
-  benchmark_smoke_environment.set('OMP_NUM_THREADS', '1')
-  benchmark_smoke_environment.set('OMP_DYNAMIC', 'FALSE')
-  if (
-    mfem_runtime_prefix != '' and dependency_prefix != '' and
-    not wheel_carries_native_bundle
-  )
-    benchmark_smoke_environment.prepend('PATH', dependency_prefix / 'bin')
-    if host_machine.system() == 'darwin'
-      benchmark_smoke_environment.prepend('DYLD_LIBRARY_PATH', dependency_prefix / 'lib64')
-      benchmark_smoke_environment.prepend('DYLD_LIBRARY_PATH', dependency_prefix / 'lib')
-    elif host_machine.system() != 'windows' and not is_wasm
-      benchmark_smoke_environment.prepend('LD_LIBRARY_PATH', dependency_prefix / 'lib64')
-      benchmark_smoke_environment.prepend('LD_LIBRARY_PATH', dependency_prefix / 'lib')
-    endif
-  endif
-  if mfem_runtime_prefix != ''
-    benchmark_smoke_environment.prepend('PATH', mfem_runtime_prefix / 'bin')
-    if host_machine.system() == 'darwin'
-      benchmark_smoke_environment.prepend('DYLD_LIBRARY_PATH', mfem_runtime_prefix / 'lib64')
-      benchmark_smoke_environment.prepend('DYLD_LIBRARY_PATH', mfem_runtime_prefix / 'lib')
-    elif host_machine.system() != 'windows'
-      benchmark_smoke_environment.prepend('LD_LIBRARY_PATH', mfem_runtime_prefix / 'lib64')
-      benchmark_smoke_environment.prepend('LD_LIBRARY_PATH', mfem_runtime_prefix / 'lib')
-    endif
-  endif
-  test(
-    'backend-benchmark-smoke',
-    python_build,
-    args: [
-      files('../tools/run_mpi_test.py'),
-      '--launcher', benchmark_mpi_launcher,
-      '--processes', '1',
-      '--',
-      benchmark_executable,
-      '--device', 'cpu',
-      '--mesh-n', '4',
-      '--order', '2',
-      '--applications', '2',
-      '--minimum-apply-seconds', '0',
-      '--warmup-applications', '1',
-      '--relative-tolerance', '1e-6',
-      '--max-iterations', '200',
-      '--solve',
-    ],
-    depends: benchmark_executable,
-    env: benchmark_smoke_environment,
-    timeout: 180,
-  )
-endif
-
-if get_option('build_tests') and get_option('build_python')
-  python_test_environment = environment()
-  python_test_environment.prepend('PYTHONPATH', python_extension_dir)
-  if mfem_runtime_prefix != ''
-    if host_machine.system() == 'darwin'
-      python_test_environment.prepend('DYLD_LIBRARY_PATH', mfem_runtime_prefix / 'lib64')
-      python_test_environment.prepend('DYLD_LIBRARY_PATH', mfem_runtime_prefix / 'lib')
-    elif host_machine.system() != 'windows'
-      python_test_environment.prepend('LD_LIBRARY_PATH', mfem_runtime_prefix / 'lib64')
-      python_test_environment.prepend('LD_LIBRARY_PATH', mfem_runtime_prefix / 'lib')
-    endif
-  endif
-  test(
-    'python-nanobind',
-    python_build,
-    args: [
-      '-c',
-      'import _core; assert _core.serial_poisson_dofs(3, 1) > 0; assert _core.capabilities()["hypre"] == _core.capabilities()["mpi"]',
-    ],
-    depends: python_extension,
-    env: python_test_environment,
-    timeout: 120,
-  )
-endif
-
 if get_option('build_tests')
-  test('cuda-toolchain-diagnostics', python_build,
+#    subdir('discritization')
-    args: files('test_cuda_toolchain.py'),
+#
-  )
+#    test_include_dir = include_directories('include')
+#
+#    catch2_dep = dependency('Catch2', required:true)
+    stroid_dep = dependency('stroid', required:true)
+#
+#    test_sources = discretization_test_sources
+#    executable('serif_tests', test_sources, dependencies: [serif_dep, catch2_dep, stroid_dep], include_directories: test_include_dir)
+
+    subdir('sandbox')
 endif

tests/sandbox/discritization_sandbox.cpp

@@ -0,0 +1,24 @@
+#include <print>
+
+#include "serif/discretization/domain/schema/schemas.hpp"
+#include "mfem.hpp"
+#include "serif/discretization/domain/schema/validation/all.hpp"
+#include "serif/discretization/domain/schema/validation/results.hpp"
+#include "stroid/stroid.h"
+#include "serif/discretization/blocks/base.hpp"
+
+int main() {
+    using serif::discretization::domain::schema::validation::SchemaValidationResult;
+    using serif::discretization::domain::schema::validation::validate_schema;
+    using serif::discretization::domain::schema::CoreEnvelopeVacuumDomainSchema;
+
+    const auto smesh_result = stroid::IO::LoadStroidMesh("sandbox.smesh");
+    if (not smesh_result.has_value()) {
+        throw std::runtime_error("Failed to load sandbox.smesh");
+    }
+
+    const stroid::StroidMesh& smesh = smesh_result.value();
+
+    SchemaValidationResult result = validate_schema<CoreEnvelopeVacuumDomainSchema>(*smesh.mesh);
+    std::println("{}", result);
+}

tests/sandbox/eos_sandbox.cpp

@@ -0,0 +1,11 @@
+#include <print>
+
+#include "serif/eos/models/polytropic.hpp"
+#include "serif/eos/evaluation.hpp"
+#include "serif/dimensions/type_alias.hpp"
+int main() {
+    const serif::eos::models::Polytrope polytrope(1, 0.6);
+    std::println("Density: {}", serif::eos::evaluate<serif::dimensions::Density>(polytrope, serif::dimensions::SpecificEnthalpyValue{1.0}).value());
+
+
+}

tests/sandbox/finite_sandbox.cpp

@@ -0,0 +1,26 @@
+#include "serif/utils/misc/finite.hpp"
+#include <limits>
+
+enum class [[maybe_unused]] TestErrorCode {
+    NotFinite
+};
+
+namespace {
+    class test_exception : public std::runtime_error {
+    public:
+        test_exception(TestErrorCode error_code, const std::string& message)
+            : std::runtime_error(message), m_error_code(error_code) {}
+
+        TestErrorCode error_code() const noexcept {
+            return m_error_code;
+        }
+    private:
+        TestErrorCode m_error_code;
+    };
+}
+
+int main() {
+    double a = std::numeric_limits<double>::infinity();
+
+    serif::utils::misc::validate_finite<test_exception>(a, TestErrorCode::NotFinite);
+}

tests/sandbox/meson.build

@@ -0,0 +1,14 @@
+sandboxes = [
+    'discritization_sandbox',
+    'eos_sandbox',
+    'finite_sandbox',
+]
+
+foreach sandbox: sandboxes
+  executable(
+    sandbox,
+    sandbox + '.cpp',
+    dependencies: [serif_dep, stroid_dep],
+    build_rpath: mfem_runtime_prefix == '' ? '' : mfem_runtime_prefix / 'lib',
+  )
+endforeach

tests/serif_tests.cpp

@@ -0,0 +1,653 @@
+#include <algorithm>
+#include <catch2/catch_session.hpp>
+#include <catch2/catch_test_case_info.hpp>
+#include <catch2/reporters/catch_reporter_registrars.hpp>
+#include <catch2/reporters/catch_reporter_streaming_base.hpp>
+#include <chrono>
+#include <cstdint>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <mfem.hpp>
+#include <regex>
+#include <sstream>
+#include <string>
+#include <string_view>
+#include <unordered_set>
+#include <utility>
+#include <vector>
+
+#include <CLI/CLI.hpp>
+#include <fourdst/config/config.h>
+
+
+std::string escapeHtml(const std::string &data) {
+    std::string buffer;
+    buffer.reserve(data.size());
+    for (size_t pos = 0; pos != data.size(); ++pos) {
+        switch (data[pos]) {
+        case '&':
+            buffer.append("&amp;");
+            break;
+        case '\"':
+            buffer.append("&quot;");
+            break;
+        case '\'':
+            buffer.append("&apos;");
+            break;
+        case '<':
+            buffer.append("&lt;");
+            break;
+        case '>':
+            buffer.append("&gt;");
+            break;
+        default:
+            buffer.append(&data[pos], 1);
+            break;
+        }
+    }
+    return buffer;
+}
+
+std::string ansiToHtml(const std::string &text) {
+    // Convert text to HTML-safe first
+    std::string htmlEscaped = escapeHtml(text);
+
+    std::ostringstream oss;
+    size_t i        = 0;
+    size_t len      = htmlEscaped.length();
+    int openSpans   = 0;
+
+    auto closeSpans = [&oss, &openSpans]() {
+        while (openSpans > 0) {
+            oss << "</span>";
+            --openSpans;
+        }
+    };
+
+    while (i < len) {
+        // Look for ANSI CSI sequence '\033[' or '\x1b['
+        if ((htmlEscaped[i] == '\033' || htmlEscaped[i] == '\x1b') && i + 1 < len && htmlEscaped[i + 1] == '[') {
+            size_t seqStart = i + 2;
+            size_t seqEnd   = htmlEscaped.find('m', seqStart);
+
+            if (seqEnd != std::string::npos) {
+                std::string codeStr = htmlEscaped.substr(seqStart, seqEnd - seqStart);
+                i                   = seqEnd + 1;
+
+                std::istringstream codeStream(codeStr);
+                std::string codeVal;
+
+                // Defaults if sequence is just \033[m (Reset)
+                if (codeStr.empty()) {
+                    closeSpans();
+                    continue;
+                }
+
+                while (std::getline(codeStream, codeVal, ';')) {
+                    int code = 0;
+                    try {
+                        code = std::stoi(codeVal);
+                    } catch (...) {
+                        continue;
+                    }
+
+                    switch (code) {
+                    case 0: // Reset
+                        closeSpans();
+                        break;
+                    case 1: // Bold
+                        oss << "<span style='font-weight:bold;'>";
+                        openSpans++;
+                        break;
+                    case 2: // Dim
+                        oss << "<span style='opacity:0.7;'>";
+                        openSpans++;
+                        break;
+                    // Standard Foreground Colors
+                    case 30:
+                        oss << "<span style='color:#2c3e50;'>";
+                        openSpans++;
+                        break; // Black
+                    case 31:
+                        oss << "<span style='color:#e74c3c;'>";
+                        openSpans++;
+                        break; // Red
+                    case 32:
+                        oss << "<span style='color:#27ae60;'>";
+                        openSpans++;
+                        break; // Green
+                    case 33:
+                        oss << "<span style='color:#f39c12;'>";
+                        openSpans++;
+                        break; // Yellow
+                    case 34:
+                        oss << "<span style='color:#2980b9;'>";
+                        openSpans++;
+                        break; // Blue
+                    case 35:
+                        oss << "<span style='color:#8e44ad;'>";
+                        openSpans++;
+                        break; // Magenta
+                    case 36:
+                        oss << "<span style='color:#16a085;'>";
+                        openSpans++;
+                        break; // Cyan
+                    case 37:
+                        oss << "<span style='color:#bdc3c7;'>";
+                        openSpans++;
+                        break; // Light Gray
+                    // Bright Foreground Colors
+                    case 90:
+                        oss << "<span style='color:#7f8c8d;'>";
+                        openSpans++;
+                        break; // Dark Gray
+                    case 91:
+                        oss << "<span style='color:#ff6b6b;'>";
+                        openSpans++;
+                        break; // Bright Red
+                    case 92:
+                        oss << "<span style='color:#51cf66;'>";
+                        openSpans++;
+                        break; // Bright Green
+                    case 93:
+                        oss << "<span style='color:#fcc419;'>";
+                        openSpans++;
+                        break; // Bright Yellow
+                    case 94:
+                        oss << "<span style='color:#339af0;'>";
+                        openSpans++;
+                        break; // Bright Blue
+                    case 95:
+                        oss << "<span style='color:#cc5de8;'>";
+                        openSpans++;
+                        break; // Bright Magenta
+                    case 96:
+                        oss << "<span style='color:#22b8cf;'>";
+                        openSpans++;
+                        break; // Bright Cyan
+                    case 97:
+                        oss << "<span style='color:#ffffff;'>";
+                        openSpans++;
+                        break; // White
+                    default:
+                        break;
+                    }
+                }
+                continue;
+            }
+        }
+        oss << htmlEscaped[i];
+        ++i;
+    }
+
+    closeSpans();
+    return oss.str();
+}
+
+std::vector<std::string> wrapText(
+    const std::string &text,
+    size_t width
+) {
+    std::vector<std::string> lines;
+    std::istringstream words(text);
+    std::string word, line;
+
+    while (words >> word) {
+        if (line.length() + word.length() + 1 > width) {
+            if (!line.empty()) {
+                lines.push_back(line);
+                line.clear();
+            }
+            if (word.length() > width) {
+                lines.push_back(word.substr(0, width - 3) + "...");
+                continue;
+            }
+        }
+        if (!line.empty())
+            line += " ";
+        line += word;
+    }
+    if (!line.empty())
+        lines.push_back(line);
+    if (lines.empty())
+        lines.push_back("");
+
+    return lines;
+}
+
+class CheckReporter : public Catch::StreamingReporterBase {
+    struct TestCaseData {
+        std::string name;
+        std::string tags;
+        bool passed;
+        std::size_t assertionsPassed;
+        std::size_t assertionsFailed;
+        double durationSeconds;
+        std::vector<std::string> failureMessages;
+        std::vector<std::string> infoMessages;
+    };
+
+    std::vector<std::string> m_currentFailures;
+    std::vector<std::string> m_currentInfos;
+    std::unordered_set<unsigned int> m_currentInfoSequences;
+    std::vector<TestCaseData> m_testRunData;
+    std::chrono::time_point<std::chrono::steady_clock> m_testStartTime;
+
+    static bool isRootProcess() {
+        int initialized = 0;
+        int finalized   = 0;
+        MPI_Initialized(&initialized);
+        if (initialized == 0) {
+            return true;
+        }
+        MPI_Finalized(&finalized);
+        if (finalized != 0) {
+            return true;
+        }
+
+        int rank = 0;
+        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+        return rank == 0;
+    }
+
+    void captureInfoMessages(Catch::AssertionStats const &assertionStats) {
+        for (auto const &message : assertionStats.infoMessages) {
+            if (m_currentInfoSequences.insert(message.sequence).second) {
+                m_currentInfos.push_back(message.message);
+            }
+        }
+    }
+
+public:
+    explicit CheckReporter(Catch::ReporterConfig &&config) : Catch::StreamingReporterBase(std::move(config)) {
+        // INFO messages are delivered through assertionEnded. Request passing
+        // assertions as well so HTML logging does not depend on Catch2's -s
+        // flag.
+        m_preferences.shouldReportAllAssertions      = true;
+
+        // This reporter does not use assertionStarting events. Disabling them
+        // preserves Catch2's successful-assertion fast path where possible.
+        m_preferences.shouldReportAllAssertionStarts = false;
+    }
+
+    static std::string getDescription() {
+        return "Console reporter with wrapping, tags, live test progress, and collapsible HTML "
+               "export with ANSI color rendering.";
+    }
+
+    void testRunStarting(Catch::TestRunInfo const &_testRunInfo) override {
+        StreamingReporterBase::testRunStarting(_testRunInfo);
+
+        if (!isRootProcess()) {
+            return;
+        }
+
+        std::cout << '\n';
+        std::cout << std::left << std::setw(85) << "Test Case Name"
+                  << "Status    " << std::right << std::setw(8) << "Passed" << std::setw(8) << "Failed" << std::setw(12)
+                  << "Time (s)" << '\n';
+        std::cout << std::string(133, '-') << '\n';
+    }
+
+    void testCaseStarting(Catch::TestCaseInfo const &testInfo) override {
+        StreamingReporterBase::testCaseStarting(testInfo);
+
+        m_testStartTime = std::chrono::steady_clock::now();
+        if (!isRootProcess()) {
+            return;
+        }
+
+        std::string name = testInfo.name;
+        auto wrappedName = wrapText(name, 83);
+
+        // Print progress line, \r to overwrite later, \033[K to clear till end of line
+        std::cout << "\r\033[K" << std::left << std::setw(85) << (wrappedName[0] + " ...") << std::flush;
+    }
+
+    void assertionEnded(Catch::AssertionStats const &assertionStats) override {
+        StreamingReporterBase::assertionEnded(assertionStats);
+
+        // Capture every INFO message encountered by either a passing or failing
+        // assertion. Message sequence IDs prevent a scoped INFO from being
+        // repeated once for every assertion that occurs while it remains
+        // active.
+        captureInfoMessages(assertionStats);
+
+        if (!assertionStats.assertionResult.isOk()) {
+            auto const &result = assertionStats.assertionResult;
+            std::ostringstream oss;
+
+            oss << "    \033[31m-> FAILED:\033[0m " << result.getSourceInfo().file << ":" << result.getSourceInfo().line
+                << '\n';
+            oss << "       " << result.getTestMacroName() << "( " << result.getExpression() << " )\n";
+
+            if (result.hasExpandedExpression()) {
+                oss << "       with expansion:\n"
+                    << "       " << result.getExpandedExpression() << '\n';
+            }
+
+            for (auto const &msg : assertionStats.infoMessages) {
+                oss << "       \033[36m[INFO]\033[0m " << msg.message << '\n';
+            }
+
+            m_currentFailures.push_back(oss.str());
+        }
+    }
+
+    void testCaseEnded(Catch::TestCaseStats const &stats) override {
+        StreamingReporterBase::testCaseEnded(stats);
+
+        auto endTime                          = std::chrono::steady_clock::now();
+        std::chrono::duration<double> elapsed = endTime - m_testStartTime;
+        double duration_s                     = elapsed.count();
+
+        bool passed                           = stats.totals.assertions.allPassed();
+        std::string mark                      = passed ? "\033[32m✓\033[0m" : "\033[31m✗\033[0m";
+
+        std::string name                      = stats.testInfo->name;
+        auto wrappedName                      = wrapText(name, 83);
+
+        if (isRootProcess()) {
+            // Overwrite the loading line with the actual result
+            std::cout << "\r\033[K" << std::left << std::setw(85) << wrappedName[0] << mark << "         " << std::right
+                      << std::setw(8) << stats.totals.assertions.passed << std::setw(8)
+                      << stats.totals.assertions.failed << std::setw(11) << std::fixed << std::setprecision(3)
+                      << duration_s << "s\n";
+
+            for (size_t i = 1; i < wrappedName.size(); ++i) {
+                std::cout << "  \033[90m↳ \033[0m" // Dim indent arrow
+                          << std::left << std::setw(81) << wrappedName[i] << '\n';
+            }
+
+            std::string tagsStr = stats.testInfo->tagsAsString();
+            if (!tagsStr.empty()) {
+                auto wrappedTags = wrapText("Tags: " + tagsStr, 83);
+                for (const auto &line : wrappedTags) {
+                    std::cout << "  \033[36m" << line << "\033[0m\n"; // Cyan
+                }
+            }
+
+            if (!m_currentFailures.empty()) {
+                std::cout << '\n';
+                for (auto const &failure : m_currentFailures) {
+                    std::cout << failure << '\n';
+                }
+                std::cout << std::string(133, '-') << '\n';
+            }
+        }
+
+        std::string tagsStr = stats.testInfo->tagsAsString();
+
+        m_testRunData.push_back(
+            {name, tagsStr, passed, stats.totals.assertions.passed, stats.totals.assertions.failed, duration_s,
+             m_currentFailures, m_currentInfos}
+        );
+
+        m_currentFailures.clear();
+        m_currentInfos.clear();
+        m_currentInfoSequences.clear();
+    }
+
+    void testRunEnded(Catch::TestRunStats const &_testRunStats) override {
+        StreamingReporterBase::testRunEnded(_testRunStats);
+
+        if (!isRootProcess()) {
+            return;
+        }
+
+        std::cout << std::string(133, '=') << '\n';
+
+        auto const &tc = _testRunStats.totals.testCases;
+        auto const &as = _testRunStats.totals.assertions;
+
+        std::string tc_passed_str =
+            tc.passed > 0 ? "\033[32m" + std::to_string(tc.passed) + " passed\033[0m" : "0 passed";
+        std::string tc_failed_str =
+            tc.failed > 0 ? "\033[31m" + std::to_string(tc.failed) + " failed\033[0m" : "0 failed";
+
+        std::string as_passed_str =
+            as.passed > 0 ? "\033[32m" + std::to_string(as.passed) + " passed\033[0m" : "0 passed";
+        std::string as_failed_str =
+            as.failed > 0 ? "\033[31m" + std::to_string(as.failed) + " failed\033[0m" : "0 failed";
+
+        std::cout << "Test Cases: " << tc_passed_str << ", " << tc_failed_str << ", " << tc.total() << " total\n";
+        std::cout << "Assertions: " << as_passed_str << ", " << as_failed_str << ", " << as.total() << " total\n\n";
+
+        generateHtmlReport(_testRunStats);
+    }
+
+private:
+    void generateHtmlReport(Catch::TestRunStats const &stats) {
+        std::ofstream html("test_summary.html");
+        if (!html)
+            return;
+
+        html << "<!DOCTYPE html>\n<html lang='en'>\n<head>\n"
+             << "<meta charset='UTF-8'>\n"
+             << "<meta name='viewport' content='width=device-width, "
+                "initial-scale=1.0'>\n"
+             << "<title>Test Run Summary</title>\n"
+             << "<style>\n"
+             << "body { font-family: -apple-system, BlinkMacSystemFont, 'Segoe "
+                "UI', "
+                "Roboto, Helvetica, Arial, sans-serif; "
+                "background: #f4f6f8; color: #333; margin: 0; padding: 2rem; }\n"
+             << "h1 { color: #2c3e50; border-bottom: 2px solid #e0e0e0; "
+                "padding-bottom: 0.5rem; }\n"
+             << ".summary-cards { display: flex; gap: 1rem; margin-bottom: "
+                "2rem; }\n"
+             << ".card { background: white; padding: 1rem 1.5rem; "
+                "border-radius: "
+                "8px; box-shadow: 0 2px 4px "
+                "rgba(0,0,0,0.05); flex: 1; }\n"
+             << ".card h3 { margin-top: 0; font-size: 0.9rem; color: #7f8c8d; "
+                "text-transform: uppercase; }\n"
+             << ".card p { font-size: 1.5rem; font-weight: bold; margin: 0; }\n"
+             << ".text-green { color: #27ae60; }\n"
+             << ".text-red { color: #e74c3c; }\n"
+             << ".test-item { background: white; border-radius: 8px; padding: "
+                "1rem; "
+                "margin-bottom: 1rem; box-shadow: 0 2px "
+                "4px rgba(0,0,0,0.05); border-left: 5px solid #bdc3c7; }\n"
+             << ".test-item.passed { border-left-color: #27ae60; }\n"
+             << ".test-item.failed { border-left-color: #e74c3c; }\n"
+             << ".test-header { display: flex; justify-content: space-between; "
+                "align-items: flex-start; }\n"
+             << ".test-name { font-size: 1.1rem; font-weight: 600; margin: 0 0 "
+                "0.5rem 0; word-break: break-word; }\n"
+             << ".tags { font-size: 0.8rem; color: #2980b9; background: "
+                "#ebf5fb; "
+                "padding: 2px 6px; border-radius: 4px; "
+                "display: inline-block; margin-top: 4px; }\n"
+             << ".stats { font-size: 0.9rem; color: #7f8c8d; }\n"
+             << "details { margin-top: 0.8rem; background: #f8f9fa; border: 1px "
+                "solid #e9ecef; border-radius: 6px; "
+                "padding: 0.5rem 0.8rem; }\n"
+             << "summary { cursor: pointer; font-weight: 600; color: #34495e; "
+                "user-select: none; font-size: 0.9rem; }\n"
+             << "summary:hover { color: #2980b9; }\n"
+             << "pre { background: #1e293b; color: #f8fafc; padding: 1rem; "
+                "border-radius: 4px; overflow-x: auto; "
+                "font-size: 0.85rem; line-height: 1.4; margin-top: 0.5rem; }\n"
+             << "pre.info-block { background: #0f172a; border-left: 4px solid "
+                "#0284c7; }\n"
+             << "</style>\n</head>\n<body>\n";
+
+        html << "<h1>Test Run Summary</h1>\n";
+
+        // Summary Cards
+        html << "<div class='summary-cards'>\n";
+        html << "<div class='card'><h3>Total Cases</h3><p>" << stats.totals.testCases.total() << "</p></div>\n";
+        html << "<div class='card'><h3>Cases Passed</h3><p class='text-green'>" << stats.totals.testCases.passed
+             << "</p></div>\n";
+        html << "<div class='card'><h3>Cases Failed</h3><p class='text-red'>" << stats.totals.testCases.failed
+             << "</p></div>\n";
+        html << "</div>\n";
+
+        for (const auto &test : m_testRunData) {
+            std::string statusClass = test.passed ? "passed" : "failed";
+            html << "<div class='test-item " << statusClass << "'>\n";
+            html << "  <div class='test-header'>\n";
+            html << "    <div>\n";
+            html << "      <h3 class='test-name'>" << escapeHtml(test.name) << "</h3>\n";
+            if (!test.tags.empty()) {
+                html << "      <div class='tags'>" << escapeHtml(test.tags) << "</div>\n";
+            }
+            html << "    </div>\n";
+            html << "    <div class='stats'>\n";
+            html << "      <span class='text-green'>&#10003; " << test.assertionsPassed << "</span> | ";
+            html << "      <span class='text-red'>&#10007; " << test.assertionsFailed << "</span> | ";
+            html << "      <span style='color: #34495e;'>&#8987; " << std::fixed << std::setprecision(3)
+                 << test.durationSeconds << "s</span>\n";
+            html << "    </div>\n";
+            html << "  </div>\n";
+
+            // Collapsible INFO Messages section with ANSI color rendering
+            if (!test.infoMessages.empty()) {
+                html << "  <details>\n";
+                html << "    <summary>Info Logs (" << test.infoMessages.size() << ")</summary>\n";
+                html << "    <pre class='info-block'>";
+                for (const auto &info : test.infoMessages) {
+                    html << "[INFO] " << ansiToHtml(info) << "\n";
+                }
+                html << "</pre>\n";
+                html << "  </details>\n";
+            }
+
+            // Collapsible Failures section with ANSI color rendering
+            if (!test.failureMessages.empty()) {
+                html << "  <details open>\n";
+                html << "    <summary class='text-red'>Failure Details (" << test.failureMessages.size()
+                     << ")</summary>\n";
+                html << "    <pre>";
+                for (const auto &msg : test.failureMessages) {
+                    html << ansiToHtml(msg) << "\n";
+                }
+                html << "    </pre>\n";
+                html << "  </details>\n";
+            }
+
+            html << "</div>\n";
+        }
+
+        html << "</body>\n</html>\n";
+    }
+};
+
+CATCH_REGISTER_REPORTER(
+    "check",
+    CheckReporter
+)
+
+int main(
+    int argc,
+    char *argv[]
+) {
+    fourdst::config::Config<mean_field::utils::Args> cfg;
+    CLI::App app{"Mean Field Tests"};
+
+    app.allow_extras();
+    app.set_help_flag("--config-help", "Show mean-field configuration options");
+    fourdst::config::register_as_cli(cfg, app);
+
+    std::vector<std::string> config_arguments;
+    std::vector<std::string> forced_catch_arguments;
+    config_arguments.emplace_back(argv[0]);
+
+    bool parsing_catch_arguments = false;
+
+    for (int i = 1; i < argc; ++i) {
+        if (std::string_view(argv[i]) == "--catch2") {
+            parsing_catch_arguments = true;
+            continue;
+        }
+
+        if (parsing_catch_arguments) {
+            forced_catch_arguments.emplace_back(argv[i]);
+        } else {
+            config_arguments.emplace_back(argv[i]);
+        }
+    }
+
+    std::vector<const char *> config_argv;
+    config_argv.reserve(config_arguments.size());
+
+    for (const std::string &argument : config_arguments) {
+        config_argv.push_back(argument.c_str());
+    }
+
+    try {
+        app.parse(static_cast<int>(config_argv.size()), config_argv.data());
+    } catch (const CLI::ParseError &error) {
+        return app.exit(error);
+    }
+
+    std::vector<std::string> catch_arguments;
+    catch_arguments.emplace_back(argv[0]);
+
+    for (const std::string &argument : app.remaining()) {
+        catch_arguments.push_back(argument);
+    }
+
+    for (const std::string &argument : forced_catch_arguments) {
+        catch_arguments.push_back(argument);
+    }
+
+    const auto is_reporter_option = [](const std::string &argument) {
+        return argument == "-r" || argument == "--reporter" || argument.starts_with("-r=") ||
+               argument.starts_with("--reporter=");
+    };
+
+    if (const bool has_reporter = std::ranges::any_of(catch_arguments, is_reporter_option); !has_reporter) {
+        catch_arguments.emplace_back("--reporter");
+        catch_arguments.emplace_back("check");
+    }
+
+    std::vector<const char *> catch_argv;
+    catch_argv.reserve(catch_arguments.size());
+
+    for (const std::string &argument : catch_arguments) {
+        catch_argv.push_back(argument.c_str());
+    }
+
+    Catch::Session session;
+
+    if (const int catch_parse_result = session.applyCommandLine(static_cast<int>(catch_argv.size()), catch_argv.data());
+        catch_parse_result != 0) {
+        return catch_parse_result;
+    }
+
+    mfem::Mpi::Init(argc, argv);
+
+    std::uint32_t synchronized_seed = session.configData().rngSeed;
+    MPI_Bcast(&synchronized_seed, 1, MPI_UINT32_T, 0, MPI_COMM_WORLD);
+    session.configData().rngSeed        = synchronized_seed;
+
+    constexpr std::string device_config = "cpu";
+    mfem::Device device(device_config);
+
+    const int hdiv_max_q1d = mfem::DeviceDofQuadLimits::Get().HDIV_MAX_Q1D;
+    if (mfem::Mpi::Root()) {
+        std::cout << "H(div) maximum Q1D = " << hdiv_max_q1d << '\n';
+        std::cout << "Approximate maximum safe integration order = " << 2 * hdiv_max_q1d - 1 << '\n';
+    }
+
+    mean_field::utils::Args test_args = cfg.main();
+
+    if (app.count("--mesh_file") == 0) {
+        test_args.mesh_file = "sandbox.smesh";
+    }
+
+    if (app.count("--p.rtol") == 0) {
+        test_args.p.rtol = 1.0e-12;
+    }
+
+    if (app.count("--p.atol") == 0) {
+        test_args.p.atol = 1.0e-12;
+    }
+
+    test_utils::set_args(std::move(test_args));
+
+    return session.run();
+}