From 87a1026be630c8d179fc0d4332c765cde7b3bf95 Mon Sep 17 00:00:00 2001
From: Emily Boudreaux <thomas@boudreauxmail.com>
Date: Mon, 9 Nov 2026 13:23:08 -0500
Subject: [PATCH] feat(operators): Moved most of operators to backend

Operators defines how SERiF concepts map into MFEM (or generally any backend) concenpts and types. Since we are trying to seperate the specific MFEM implementation out of SERiF most of the code that was here has moved to backend. However, what remains is some general code related to the block layout of a systen and the contributions of different equations to different parts of that block layout
---
 src/include/serif/operators/concepts.hpp      | 100 ++++++++++++++++++
 src/include/serif/operators/contributions.hpp |  85 +++++++++++++++
 .../serif/operators/linearized_block.hpp      |   8 ++
 .../operators/linearized_system_view.hpp      |  64 +++++++++++
 src/include/serif/operators/plan.hpp          |  85 +++++++++++++++
 5 files changed, 342 insertions(+)
 create mode 100644 src/include/serif/operators/concepts.hpp
 create mode 100644 src/include/serif/operators/contributions.hpp
 create mode 100644 src/include/serif/operators/linearized_block.hpp
 create mode 100644 src/include/serif/operators/linearized_system_view.hpp
 create mode 100644 src/include/serif/operators/plan.hpp

diff --git a/src/include/serif/operators/concepts.hpp b/src/include/serif/operators/concepts.hpp
new file mode 100644
index 00000000..a8dba3a3
--- /dev/null
+++ b/src/include/serif/operators/concepts.hpp
@@ -0,0 +1,100 @@
+#pragma once
+
+#include <concepts>
+#include <type_traits>
+
+#include "serif/discretization/blocks/traits.hpp"
+#include "serif/discretization/functional_space/space.hpp"
+#include "serif/equations/concepts.hpp"
+
+namespace serif::operators {
+    // These restrictions describe the existing volume implementation,
+    // not which physical equations may be declared in serif::equations.
+    template <typename TermT>
+    concept ScalarValueTerm =
+        discretization::blocks::GlobalScalarTerm<TermT> ||
+        (discretization::blocks::SpatialTerm<TermT> &&
+        discretization::blocks::SpatialTermTraits<TermT>::rank == 0 &&
+        (
+            discretization::blocks::SpatialTermTraits<TermT>::family == discretization::functional_space::Family::h1 ||
+            discretization::blocks::SpatialTermTraits<TermT>::family == discretization::functional_space::Family::l2
+        ));
+
+    template <typename OperandT>
+    concept ScalarValueOperand =
+        discretization::forms::IsOperand<OperandT> &&
+        ScalarValueTerm<typename OperandT::term_type> &&
+        std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Value>;
+
+    /// A conforming scalar gradient; discontinuous fields also require face terms.
+    template <typename OperandT>
+    concept ScalarGradientOperand =
+        discretization::forms::IsOperand<OperandT> &&
+        discretization::blocks::SpatialTerm<typename OperandT::term_type> &&
+        discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::rank == 0 &&
+        discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::family == discretization::functional_space::Family::h1 &&
+        std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Gradient>;
+
+    /** @brief Full vector gradients are supported for componentwise H1 fields. */
+    template <typename OperandT>
+    concept VectorGradientOperand =
+        discretization::forms::IsOperand<OperandT> &&
+        discretization::blocks::SpatialTerm<typename OperandT::term_type> &&
+        discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::rank == 1 &&
+        discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::family == discretization::functional_space::Family::h1 &&
+        std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Gradient>;
+
+    template <typename EquationT, typename OperandsT>
+    struct HasScalarVolumeOperands : std::false_type {};
+
+    template <typename EquationT, typename... Operands>
+    struct HasScalarVolumeOperands<EquationT, utils::types::TypeList<Operands...>> :
+        std::bool_constant<
+            sizeof...(Operands) != 0 &&
+            (ScalarValueOperand<Operands> && ...) &&
+            (requires(const EquationT& equation, Operands operand) {
+                {equation.additional_order(operand, 1)}->std::same_as<int>;
+            } && ...)> {};
+
+    /// Equations executable by the current scalar-value volume operator.
+    template <typename EquationT>
+    concept ScalarVolumeEquation =
+        equations::Equation<EquationT> &&
+        ScalarValueTerm<typename EquationT::Test> &&
+        HasScalarVolumeOperands<EquationT, typename EquationT::Operands>::value;
+    /// Mathematical operand combinations implemented by the current volume backend.
+    template <typename OperandT>
+    concept MappedOperand =
+        discretization::forms::IsOperand<OperandT> &&
+        (
+            ScalarValueOperand<OperandT> || ScalarGradientOperand<OperandT> || VectorGradientOperand<OperandT> ||
+            (
+                discretization::blocks::SpatialTerm<typename OperandT::term_type> &&
+                discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::rank == 1 &&
+                discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::family == discretization::functional_space::Family::h1 &&
+                (
+                    std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Value> ||
+                    std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Divergence>
+                )
+            )
+        );
+
+    template <typename C, typename ListT> struct VolumeOperands : std::false_type {};
+    template <typename C, typename... O>
+    struct VolumeOperands<C, utils::types::TypeList<O...>> : std::bool_constant<
+        sizeof...(O) != 0 && (MappedOperand<O> && ...) &&
+        (requires(const C& contribution, O operand) {
+        {contribution.additional_order(operand, 1)}->std::same_as<int>; } && ...)> {};
+
+    template <typename TupleT> struct VolumeContributions : std::false_type {};
+    template <typename... C>
+    struct VolumeContributions<std::tuple<C...>>
+        : std::bool_constant<((discretization::domain::IsDomainOrSet<equations::SupportOf<C>> &&
+                               MappedOperand<equations::TestOperandOf<C>> && VolumeOperands<C, typename C::Operands>::value) &&
+                              ...)> {};
+
+    /// One or more supported weak volume contributions assembled into a single residual row.
+    template <typename T>
+    concept VolumeEquation = equations::Equation<T> && VolumeContributions<equations::ContributionsOf<T>>::value;
+
+} // namespace serif::operators
diff --git a/src/include/serif/operators/contributions.hpp b/src/include/serif/operators/contributions.hpp
new file mode 100644
index 00000000..e35b7328
--- /dev/null
+++ b/src/include/serif/operators/contributions.hpp
@@ -0,0 +1,85 @@
+#pragma once
+
+#include <cstddef>
+#include <memory>
+#include <type_traits>
+#include <utility>
+#include <tuple>
+
+#include "serif/equations/concepts.hpp"
+
+namespace serif::operators {
+    /// Retain one leaf equation and the constant coefficient of its arithmetic expression.
+    template <typename EquationT>
+    struct ContributionBinding {
+        std::shared_ptr<const EquationT> equation;
+        double coefficient;
+    };
+
+    namespace detail {
+        /// A returned tuple may own generated equations or borrow an equation's existing children.
+        /// Retaining its parent keeps reference tuples valid without copying move-only equations.
+        template <typename EquationT>
+        struct ContributionChildren {
+            std::shared_ptr<const EquationT> parent;
+            decltype(std::declval<const EquationT&>().contributions()) values;
+
+            /// Evaluate contributions once, retaining either its returned value or its returned reference.
+            explicit ContributionChildren(std::shared_ptr<const EquationT> equation) :
+                parent(std::move(equation)),
+                values(parent->contributions()) {
+            }
+        };
+
+        /// Flatten equation arithmetic while preserving every generated tuple and borrowed parent.
+        template <typename EquationT>
+        [[nodiscard]] auto bind_contributions(std::shared_ptr<const EquationT> equation, double coefficient) {
+            if constexpr (requires { equation->contributions(); }) {
+                const auto children = std::make_shared<const ContributionChildren<EquationT>>(std::move(equation));
+                return std::apply([&](const auto&... child) {
+                    // Aliasing pointers address each child but retain the tuple and its parent together.
+                    return std::tuple_cat(bind_contributions(
+                        std::shared_ptr<const std::remove_cvref_t<decltype(child)>>(children, &child), coefficient
+                    )...);
+                }, children->values);
+            } else {
+                return std::tuple{ContributionBinding<EquationT>{
+                    .equation = std::move(equation),
+                    .coefficient = coefficient
+                }};
+            }
+        }
+
+        /// Visit arithmetic leaves without retaining them beyond the supplied callback.
+        template <typename EquationT, typename FunctionT>
+        void visit_contributions(
+            const EquationT& equation,
+            FunctionT& function,
+            std::size_t& index,
+            double coefficient
+        ) {
+            if constexpr (requires { equation.contributions(); }) {
+                std::apply([&](const auto&... child) {(visit_contributions(child, function, index, coefficient), ...);},equation.contributions());
+            } else {
+                function(equation, index++, coefficient);
+            }
+        }
+    } // namespace detail
+
+    /// Give each retained leaf a compile-time tuple position for independently cached contribution results.
+    template <equations::Equation EquationT>
+    [[nodiscard]] auto bind_contributions(std::shared_ptr<const EquationT> equation) {
+        return detail::bind_contributions(std::move(equation), 1.0);
+    }
+
+    /// Visit leaves in declaration order, multiplying constant coefficients along each path.
+    /// The callback receives (leaf equation, unique contribution index, coefficient).
+    template <equations::Equation EquationT, typename FunctionT>
+    void for_each_contribution(
+        const EquationT& equation,
+        FunctionT&& function
+    ) {
+        std::size_t index = 0;
+        detail::visit_contributions(equation, function, index, 1.0);
+    }
+} // namespace serif::operators
diff --git a/src/include/serif/operators/linearized_block.hpp b/src/include/serif/operators/linearized_block.hpp
new file mode 100644
index 00000000..60428d36
--- /dev/null
+++ b/src/include/serif/operators/linearized_block.hpp
@@ -0,0 +1,8 @@
+#pragma once
+
+#include "serif/backend/mfem/operators/linearized_block.hpp"
+
+namespace serif::operators {
+    using backend::mfem::LinearizedBlock;
+    using backend::mfem::linearized_block;
+}
diff --git a/src/include/serif/operators/linearized_system_view.hpp b/src/include/serif/operators/linearized_system_view.hpp
new file mode 100644
index 00000000..d8a4a42b
--- /dev/null
+++ b/src/include/serif/operators/linearized_system_view.hpp
@@ -0,0 +1,64 @@
+#pragma once
+
+#include <type_traits>
+
+#include "serif/discretization/blocks/system_layout.hpp"
+#include "serif/utils/dependencies/dependency_context.hpp"
+
+namespace serif::operators
+{
+    /** @brief Identify a physical Jacobian block by its equation rows and unknown coordinates. */
+    template <typename RowsT, typename ValuesT>
+    struct PhysicalBlock final
+    {
+        using Rows = RowsT;
+        using Values = ValuesT;
+    };
+
+    /** @brief Borrow an accepted collection of physical blocks and its independently ordered system layout. */
+    template <discretization::blocks::SystemLayoutType Layout, utils::dependencies::IsDependencyContext Blocks>
+    class LinearizedSystemView final
+    {
+    public:
+        using LayoutType = std::remove_cvref_t<Layout>;
+        using BlockSources = std::remove_cvref_t<Blocks>;
+
+        /** @brief Both owners must outlive the view and remain at this accepted linearization while it is used. */
+        LinearizedSystemView(const LayoutType& layout, const BlockSources& blocks) : m_layout(layout), m_blocks(blocks)
+        {
+        }
+
+        /** @brief A borrowed system cannot retain temporary coordinate metadata. */
+        LinearizedSystemView(const LayoutType&&, const BlockSources&) = delete;
+        /** @brief A borrowed system cannot retain temporary block sources. */
+        LinearizedSystemView(const LayoutType&, const BlockSources&&) = delete;
+
+        /** @brief Access a declared physical coupling; absent blocks are compile-time errors. */
+        template <typename Rows, typename Values>
+            requires utils::dependencies::BindingsWithBinding<typename BlockSources::Bindings, PhysicalBlock<Rows, Values>>
+        [[nodiscard]] const auto& block() const
+        {
+            return m_blocks.template get<PhysicalBlock<Rows, Values>>();
+        }
+
+        /** @brief Access the shared coordinate description without taking ownership. */
+        [[nodiscard]] const LayoutType& layout() const noexcept
+        {
+            return m_layout;
+        }
+
+        /** @brief Supply accepted block sources to ordinary dependency preparation. */
+        [[nodiscard]] const BlockSources& blocks() const noexcept
+        {
+            return m_blocks;
+        }
+
+    private:
+        const LayoutType& m_layout;
+        const BlockSources& m_blocks;
+    };
+
+    /** @brief Infer canonical owner types while retaining the explicit borrowing lifetime checks. */
+    template <discretization::blocks::SystemLayoutType Layout, utils::dependencies::IsDependencyContext Blocks>
+    LinearizedSystemView(const Layout&, const Blocks&) -> LinearizedSystemView<Layout, Blocks>;
+} // namespace serif::operators
diff --git a/src/include/serif/operators/plan.hpp b/src/include/serif/operators/plan.hpp
new file mode 100644
index 00000000..4ff219f0
--- /dev/null
+++ b/src/include/serif/operators/plan.hpp
@@ -0,0 +1,85 @@
+#pragma once
+
+#include <tuple>
+#include "serif/utils/dependencies/key_list_operations.hpp"
+#include "serif/discretization/forms/forms.hpp"
+#include "serif/operators/concepts.hpp"
+
+namespace serif::operators {
+    namespace detail {
+        template <typename ListT> struct AsKeyList;
+        template <typename... T> struct AsKeyList<utils::types::TypeList<T...>> {
+            using type = utils::dependencies::CombinedKeyList<utils::dependencies::KeyList<T>...>;
+        };
+        template <typename... Lists> struct Join {
+            using type = typename utils::dependencies::CombinedKeyList<typename AsKeyList<Lists>::type...>::Keys;
+        };
+        template <typename TupleT> struct ContributionOperands;
+        template <typename... C> struct ContributionOperands<std::tuple<C...>> : Join<typename C::Operands...> {};
+        template <typename ListT> struct Columns;
+        template <typename... O>
+        struct Columns<utils::types::TypeList<O...>> {
+            using type = typename utils::dependencies::CombinedKeyList<utils::dependencies::KeyList<discretization::forms::Operand<typename O::term_type>>...>::Keys;
+        };
+
+        template <typename Input, typename TestOperand> struct CouplingForm {
+            using Trial = typename Input::term_type;
+            using Test = typename TestOperand::term_type;
+            using type = std::conditional_t<
+                std::same_as<typename Input::operation_type, discretization::forms::Operation::Divergence>,
+                discretization::forms::divergence<Trial, Test>,
+                std::conditional_t<std::same_as<typename TestOperand::operation_type, discretization::forms::Operation::Divergence>,
+                                   discretization::forms::divergence<Test, Trial>,
+                                   std::conditional_t<std::same_as<Trial, Test>, discretization::forms::mass<Trial, Test>,
+                                                      discretization::forms::source<Trial, Test>>>>;
+        };
+    } // namespace detail
+
+    template <VolumeEquation EquationT,
+              typename OperandsT = typename detail::ContributionOperands<equations::ContributionsOf<EquationT>>::type>
+    struct EquationPlan;
+
+    template <VolumeEquation EquationT, typename... O> struct EquationPlan<EquationT, utils::types::TypeList<O...>> {
+        using Test = typename EquationT::Test;
+        using Support = discretization::blocks::SupportOf<Test>;
+        using Operands = utils::types::TypeList<O...>;
+        using Columns = typename detail::Columns<Operands>::type;
+        static constexpr std::size_t operand_count = sizeof...(O);
+        static constexpr std::size_t column_count = equations::OperandValues<Columns>::count;
+        static constexpr int component_count = (discretization::forms::component_count<discretization::forms::OperandValue<O>> + ... + 0);
+
+        template <typename OperandT, typename ContributionT = EquationT>
+        using Form = typename detail::CouplingForm<OperandT, equations::TestOperandOf<ContributionT>>::type;
+
+        /// Locate the unique coefficient column shared by every operation on a term.
+        template <typename OperandT> static consteval std::size_t column_index() {
+            return equations::OperandValues<Columns>::template index<typename OperandT::term_type>();
+        }
+        /// Locate the first scalar component of an operand in pointwise derivative storage.
+        template <typename OperandT> static consteval int component_offset() {
+            constexpr std::array counts{discretization::forms::component_count<discretization::forms::OperandValue<O>>...};
+            constexpr auto index = equations::OperandValues<Operands>::template index<OperandT>();
+            int offset = 0;
+            for (std::size_t i = 0; i < index; ++i) {
+                offset += counts[i];
+            }
+            return offset;
+        }
+        /// Visit all sampled operations in stable first-appearance order.
+        template <typename FunctionT> static void for_each_operand(FunctionT&& function) {
+            (function(O{}, equations::OperandValues<Operands>::template index<O>()), ...);
+        }
+        /// Visit only the operands needed by this contribution, retaining global sample indices.
+        template <typename ContributionT, typename FunctionT> static void for_each_input(FunctionT&& function) {
+            [&]<typename... Input>(utils::types::TypeList<Input...>) {
+                (function(Input{}, equations::OperandValues<Operands>::template index<Input>()), ...);
+            }(typename ContributionT::Operands{});
+        }
+        /// Visit each physical coefficient vector once for binding and packed actions.
+        template <typename FunctionT> static void for_each_column(FunctionT&& function) {
+            [&]<typename... Column>(utils::types::TypeList<Column...>) {
+                (function(Column{}, equations::OperandValues<Columns>::template index<Column>()), ...);
+            }(Columns{});
+        }
+    };
+} // namespace serif::operators