feat(discretization): Added SystemLayout and associated helpers to allow for more robust access to parts of the block system

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

@@ -14,6 +14,17 @@ namespace serif::discretization::blocks {
      */
     struct block {};
 
+    /**
+     * @enum RepresentedValueKind
+     * @brief Enum representing the kind of value a block represents, either a physical quantity (value) or a residual (residual).
+     *
+     * Generally this is mostly used in the solver/normalization system.
+     */
+    enum class RepresentedValueKind {
+        value,             ///< A value representing a physical quantity (e.g., temperature, pressure).
+        residual           ///< A residual representing the difference between the current state and the desired state in a numerical method.
+    };
+
     /**
      * @struct residual_block_base
      * @brief Base class for residual blocks, derived from block.

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

@@ -0,0 +1,16 @@
+#pragma once
+
+/**
+ * @file keys.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief Keys for the shared layout of a discretized state and residual in a block-structured system.
+ */
+
+
+namespace serif::discretization::blocks::keys {
+    /**
+     * @brief Share one independently ordered state and residual layout across numerical services.
+     */
+    struct Layout final {};
+}

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

@@ -0,0 +1,47 @@
+#pragma once
+
+#include <limits>
+#include <vector>
+#include "serif/discretization/blocks/system_layout.hpp"
+#include "serif/utils/error/errors.hpp"
+
+namespace serif::discretization::blocks {
+    /**
+     * @brief Select named packed coordinates in their requested order, independently for values and residuals.
+     *
+     * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals.
+     * @tparam Layout The block layout type, e.g. SystemLayout.
+     * @tparam Coordinates The coordinate tags to select, e.g. keys::Density,
+     *
+     * @param layout The block layout from which to select coordinates.
+     * @param keys A list of coordinate tags to select. Note that this is not used at runtime, but rather to select the correct overloads of layout.size and layout.offset.
+     *
+     * @returns A vector of packed indices corresponding to the requested coordinates in the order requested.
+     *
+     * Given some layout (see @ref SystemLayout) and a list of coordinate tags return the corresponding packed indices in the order requested.
+     * The result is suitable for indexing a contiguous vector of values or residuals.
+     */
+    template <RepresentedValueKind Kind, typename Layout, typename... Coordinates>
+    [[nodiscard]] std::vector<int> coordinate_indices(
+        const Layout& layout,
+        utils::dependencies::KeyList<Coordinates...> /*keys*/
+    ) {
+        utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_DOF_COUNT>(layout.template total_size<Kind>() <= static_cast<std::size_t>(std::numeric_limits<int>::max()),"Named coordinate selection exceeds the backend index range; use a layout with representable packed indices.");
+
+        std::vector<int> result;
+        result.reserve((std::size_t{0} + ... + layout.template size<Kind, Coordinates>()));
+
+        const auto append = [&]<typename Coordinate>() {
+            const int offset = static_cast<int>(layout.template offset<Kind, Coordinate>());
+            const int size = static_cast<int>(layout.template size<Kind, Coordinate>());
+
+            for (int index = 0; index < size; ++index) {
+                result.push_back(offset + index);
+            }
+        };
+
+        (append.template operator()<Coordinates>(), ...);
+
+        return result;
+    }
+}

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

@@ -0,0 +1,188 @@
+#pragma once
+
+/**
+ * @file system_layout.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief Describe the independently ordered physical state and residual coordinates of a discrete system. This is an
+ * accounting layer to keep index offsets straight in memory.
+ */
+
+#include <array>
+#include <cstddef>
+#include <span>
+#include <type_traits>
+
+#include "serif/discretization/blocks/base.hpp"
+#include "serif/utils/dependencies/key_list_operations.hpp"
+
+namespace serif::discretization::blocks
+{
+    /**
+     * @brief Sum some memory span with checking for if any element is too large to fit in a size_t
+     *
+     * @param sizes A span of sizes to sum
+     * @return The sum of the sizes
+     *
+     * @throws std::overflow_error If an element of sizes is too large to fit in a size_t
+     *
+     * @todo Really this does not need to be specific to this module, this might more generally be named something like "safe sum" and moved to a more general location.
+     */
+    [[nodiscard]] std::size_t checked_block_size(std::span<const std::size_t> sizes);
+
+    /**
+     * @brief Describe the independently ordered physical state and residual coordinates of a discrete system.
+     *
+     * @tparam ValuesT A KeyList of the physical coordinates represented in the state vector.
+     * @tparam ResidualsT A KeyList of the physical coordinates represented in the residual vector.
+     *
+     */
+    template <utils::dependencies::IsKeyList ValuesT, utils::dependencies::IsKeyList ResidualsT>
+    struct SystemLayout final {
+        using Values = std::remove_cvref_t<ValuesT>;                ///< The physical coordinates represented in the state vector.
+        using Residuals = std::remove_cvref_t<ResidualsT>;          ///< The physical coordinates represented in the residual vector.
+
+        /**
+         * @brief Select the coordinate tag at compile time based on the kind of the represented values selected (e.g. RepresentedValueKind::values or RepresentedValueKind::residuals).
+         *
+         * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals.
+         *
+         * @requires Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual
+         *
+         * @note This requires clause is not strictly necessary right now since we have a template type argument; however, if we ever add
+         * more RepresentedValueKind types, this will ensure that we don't accidentally allow them to be used here.
+         */
+        template <RepresentedValueKind Kind>
+            requires(Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual)
+        using Coordinates = std::conditional_t<Kind == RepresentedValueKind::value, Values, Residuals>;
+
+        std::array<std::size_t, utils::dependencies::key_count_v<Values>> state_sizes;         ///< The coefficient counts for each physical coordinate in the state vector, in the order of Values.
+        std::array<std::size_t, utils::dependencies::key_count_v<Residuals>> residual_sizes;   ///< The coefficient counts for each physical coordinate in the residual vector, in the order of Residuals.
+
+        /**
+         * @brief Get the sizes of the physical coordinates for the requested kind.
+         *
+         * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals.
+         * @returns A span of sizes for the physical coordinates of the requested kind.
+         *
+         * @requires Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual
+         */
+        template <RepresentedValueKind Kind>
+            requires(Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual)
+        [[nodiscard]] std::span<const std::size_t> sizes() const noexcept {
+            if constexpr (Kind == RepresentedValueKind::value) {
+                return state_sizes;
+            }
+            else {
+                return residual_sizes;
+            }
+        }
+
+        /**
+         * @brief Locate a physical identity on the requested side of the system.
+         *
+         * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals.
+         * @tparam Coordinate The coordinate type to locate.
+         * @returns The index of the coordinate in the requested kind.
+         *
+         * @requires Coordinate is contained in the KeyList of the requested kind.
+         */
+        template <RepresentedValueKind Kind, typename Coordinate>
+            requires utils::dependencies::ContainsKey<Coordinate, Coordinates<Kind>>
+        [[nodiscard]] static consteval std::size_t index() {
+            if constexpr (Kind == RepresentedValueKind::value) {
+                return utils::dependencies::key_index<Coordinate, Values>();
+            }
+            else {
+                return utils::dependencies::key_index<Coordinate, Residuals>();
+            }
+        }
+
+        /**
+         * @brief Read the coefficient count for a physical coordinate.
+         *
+         * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals.
+         * @tparam Coordinate The coordinate type to locate.
+         * @returns The coefficient count for the specified coordinate.
+         *
+         * @requires Coordinate is contained in the KeyList of the requested kind.
+         */
+        template <RepresentedValueKind Kind, typename Coordinate>
+            requires utils::dependencies::ContainsKey<Coordinate, Coordinates<Kind>>
+        [[nodiscard]] std::size_t size() const noexcept {
+            return sizes<Kind>()[index<Kind, Coordinate>()]; // A note from EMB on Oct 14 2026 to a future dev. Maybe I am sorry I picked C++? We end up with stuff like the syntax on the left... But perhaps we would one way or the other.
+        }
+
+        /**
+         * @brief Locate a coordinate in packed storage, rejecting size overflow before creating a view. Note
+         * that this specifically locates the first coefficient of a coordinate, or the offset to where that coordinate's coefficients begin in packed storage.
+         *
+         * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals.
+         * @tparam Coordinate The coordinate type to locate.
+         * @returns The offset of the coordinate in the requested kind. (e.g. the index of the first coefficient of that coordinate in packed storage.)
+         *
+         * @requires Coordinate is contained in the KeyList of the requested kind.
+         */
+        template <RepresentedValueKind Kind, typename Coordinate>
+            requires utils::dependencies::ContainsKey<Coordinate, Coordinates<Kind>>
+        [[nodiscard]] std::size_t offset() const {
+            return checked_block_size(sizes<Kind>().first(index<Kind, Coordinate>()));
+        }
+
+        /**
+         * @brief Count all coefficients of one kind
+         *
+         * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals.
+         * @returns The total coefficient count for the requested kind.
+         */
+        template <RepresentedValueKind Kind>
+            requires(Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual)
+        [[nodiscard]] std::size_t total_size() const {
+            return checked_block_size(sizes<Kind>());
+        }
+
+        /**
+         * @brief Default equality comparison
+         */
+        [[nodiscard]] bool operator==(const SystemLayout&) const = default;
+    };
+
+    namespace detail
+    {
+        /**
+         * @brief Default false specialization for is_system_layout, specialized to true for SystemLayout<Values, Residuals>.
+         * @tparam Candidate The type to check for being a SystemLayout.
+         */
+        template <typename Candidate>
+        inline constexpr bool is_system_layout = false;
+
+        /**
+         * @brief True specialization for is_system_layout, specialized to true for SystemLayout<Values, Residuals>.
+         * @tparam Values The KeyList of physical coordinates represented in the state vector.
+         * @tparam Residuals The KeyList of physical coordinates represented in the residual vector.
+         */
+        template <typename Values, typename Residuals>
+        inline constexpr bool is_system_layout<SystemLayout<Values, Residuals>> = true;
+    } // namespace detail
+
+    /**
+     * @brief Compile time check for a SystemLayout type, e.g. SystemLayout<Values, Residuals>.
+     */
+    template <typename Candidate>
+    concept SystemLayoutType = detail::is_system_layout<std::remove_cvref_t<Candidate>>;
+
+    /**
+     * @brief Build both coordinate orderings from the same physical size resolver.
+     * @param size A callable whose templated operator()<Coordinate>() returns that coordinate's coefficient count.
+     */
+    template <SystemLayoutType Layout, typename SizeResolver>
+    [[nodiscard]] constexpr std::remove_cvref_t<Layout> make_system_layout(const SizeResolver& size) {
+        using Result = std::remove_cvref_t<Layout>;
+
+        const auto sizes = [&]<typename... Coordinates>(utils::dependencies::KeyList<Coordinates...>) {
+            return std::array<std::size_t, sizeof...(Coordinates)>{static_cast<std::size_t>(size.template operator()<Coordinates>())...};
+        };
+
+        return {.state_sizes = sizes(typename Result::Values{}), .residual_sizes = sizes(typename Result::Residuals{})};
+    }
+} // namespace serif::discretization::blocks

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

@@ -67,4 +67,7 @@ namespace serif::discretization::blocks {
         static constexpr int order = FunctionalSpace::order; ///< The order of the functional space associated with the term. e.g. 1st order, 2nd order, etc.
         static constexpr std::size_t rank = FunctionalSpace::rank::value; ///< The rank of the functional space associated with the term.
     };
-}
+    /** @brief A unique scalar boundary value is available for continuous scalar fields and scalar globals. */
+    template <typename TermT>
+    concept ScalarTraceTerm = GlobalScalarTerm<TermT> || (SpatialTerm<TermT> && SpatialTermTraits<TermT>::rank == 0 && SpatialTermTraits<TermT>::family == functional_space::Family::h1);
+}

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

@@ -0,0 +1,113 @@
+#pragma once
+
+/**
+ * @file vector_view.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief This file includes a class for viewing some SystemLayout as a span of coefficients without copying the underlying storage.
+ */
+
+#include <span>
+#include <ranges>
+#include <stdexcept>
+#include <type_traits>
+
+#include "serif/discretization/blocks/system_layout.hpp"
+#include "serif/utils/error/errors.hpp"
+
+namespace serif::discretization::blocks
+{
+    /**
+     * @brief Build a non-owning view of a packed layout as a span.
+     *
+     * @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals.
+     * @tparam Layout The block layout type, e.g. SystemLayout.
+     * @tparam Scalar The scalar type of the coefficients, e.g. double.
+     *
+     * @requires Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual
+     *
+     * @note The layout must remain unchanged while the view is used. Probably we want to, in the future, add some way to
+     * validate this; however, I am at the moment at a loss for how to do so.
+     */
+    template <RepresentedValueKind Kind, SystemLayoutType Layout, typename Scalar = double>
+        requires(Kind == RepresentedValueKind::value || Kind == RepresentedValueKind::residual)
+    class BlockVectorView final
+    {
+    public:
+        using LayoutType = std::remove_cvref_t<Layout>;    ///< Alias for the layout type, e.g. SystemLayout<Values, Residuals>.
+
+        /**
+         * @brief Construct one of these views from a layout and a span of coefficients.
+         *
+         * @param layout The layout to use for the view.
+         * @param values The span of coefficients to view.
+         *
+         * @throws std::domain_error If the span of coefficients does not match the total size of the layout.
+         *
+         * @todo Really this could go in many places; however, I think we should move away from throwing std exceptions and introduce our own exception hierarchy.
+         */
+        BlockVectorView(
+            const LayoutType& layout,
+            std::span<Scalar> values
+        ) :
+        m_layout(layout),
+        m_values(values) {
+            const std::size_t required = layout.template total_size<Kind>();
+
+            utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_DOF_COUNT, std::domain_error>(values.size() == required, "Packed storage has {} coefficients but the system layout requires {}. Bind the matching state or residual vector.", values.size(), required);
+        }
+
+        /**
+         * @brief Reject temporary metadata whose offsets would dangle before the coefficients are accessed.
+         */
+        BlockVectorView(const LayoutType&&, std::span<Scalar>) = delete;
+
+        /**
+         * @brief Prevent implicit const-span conversion from retaining a temporary owning container.
+         */
+        template <typename Range>
+            requires(!std::is_lvalue_reference_v<Range> && !std::ranges::borrowed_range<Range>)
+        BlockVectorView(const LayoutType&, Range&&) = delete;
+
+        /**
+         * @brief Access the requested physical coordinate while retaining storage constness.
+         *
+         * @tparam Coordinate The coordinate tag to access
+         *
+         * @requires A layout type with kind Kind must contain the requested coordinate tag.
+         */
+        template <typename Coordinate>
+            requires utils::dependencies::ContainsKey<Coordinate, typename LayoutType::template Coordinates<Kind>>
+        [[nodiscard]] std::span<Scalar> get() const
+        {
+            const std::size_t offset = m_layout.template offset<Kind, Coordinate>();
+            const std::size_t count = m_layout.template size<Kind, Coordinate>();
+
+            utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_DOF_COUNT>(offset <= m_values.size() && count <= m_values.size() - offset,"The system layout changed after this vector view was bound. Rebind the view to matching storage.");
+
+            return m_values.subspan(offset, count);
+        }
+
+    private:
+        const LayoutType& m_layout; ///< The layout to use for the view. The layout must remain unchanged while the view is used.
+        std::span<Scalar> m_values; ///< The span of coefficients to view.
+    };
+
+    /**
+     * @brief A convenient alias for a view of the physical state coefficients.
+     *
+     * @tparam Layout The block layout type, e.g. SystemLayout.
+     * @tparam Scalar The scalar type of the coefficients, e.g. double.
+     */
+    template <SystemLayoutType Layout, typename Scalar = double>
+    using StateView = BlockVectorView<RepresentedValueKind::value, Layout, Scalar>;
+
+    /**
+     * @brief A convenient alias for a view of the residual coefficients.
+     *
+     * @tparam Layout The block layout type, e.g. SystemLayout.
+     * @tparam Scalar The scalar type of the coefficients, e.g. double.
+     */
+    template <SystemLayoutType Layout, typename Scalar = double>
+    using ResidualView = BlockVectorView<RepresentedValueKind::residual, Layout, Scalar>;
+} // namespace serif::discretization::blocks