diff --git a/include/fesa/assembly/sparse_assembler.h b/include/fesa/assembly/sparse_assembler.h
index 1bf8a76..fe1103f 100644
--- a/include/fesa/assembly/sparse_assembler.h
+++ b/include/fesa/assembly/sparse_assembler.h
@@ -2,6 +2,7 @@
 #define FESA_ASSEMBLY_SPARSE_ASSEMBLER_H_
 
 #include "fesa/core/status.h"
+#include "fesa/elements/element.h"
 #include "fesa/math/sparse_matrix.h"
 
 namespace fesa {
@@ -13,10 +14,22 @@ class ParallelFor;
 /// @brief Owns deterministic element-contribution reduction into global CSR.
 class SparseAssembler {
  public:
+  /// @brief Assembles runtime element stiffness into validated full-DOF CSR.
+  /// @param elements Non-owning elements in stable active source order.
+  /// @param dofs Owner of the matching scatter and structural pattern.
+  /// @param parallel_for Backend for index-owned element-local computation.
+  /// @return A validated matrix or a structured model failure.
+  /// @note Workers write only their element-owned COO buffers; flattening and
+  /// duplicate reduction retain fixed element and local-entry order.
+  static Result<SparseMatrix> Assemble(const ElementView& elements,
+                                       const DofManager& dofs,
+                                       const ParallelFor& parallel_for);
+
   /// @brief Assembles stiffness in stable source-element and local-entry order.
   /// @return A validated full-DOF CSR matrix or structured model failure.
   /// @note Parallel workers produce index-owned local buffers; serial reduction
-  /// remains the sole global CSR writer.
+  /// remains the sole global CSR writer. This compatibility facade creates
+  /// runtime candidates until LinearStaticAnalysis owns them directly.
   static Result<SparseMatrix> AssembleStiffness(
       const AnalysisModel& model, const DofManager& dofs,
       const ParallelFor& parallel_for);
diff --git a/src/fesa/assembly/sparse_assembler.cpp b/src/fesa/assembly/sparse_assembler.cpp
index 2a0e176..61face6 100644
--- a/src/fesa/assembly/sparse_assembler.cpp
+++ b/src/fesa/assembly/sparse_assembler.cpp
@@ -1,286 +1,97 @@
 #include "fesa/assembly/sparse_assembler.h"
 
-#include <array>
+#include <cstddef>
+#include <functional>
 #include <limits>
+#include <memory>
 #include <optional>
-#include <stdexcept>
 #include <string>
 #include <utility>
 #include <vector>
 
 #include "fesa/analysis/analysis_model.h"
 #include "fesa/assembly/parallel_for.h"
-#include "fesa/elements/euler_beam_3d.h"
-#include "fesa/elements/mitc4_shell.h"
+#include "fesa/elements/element_factory.h"
 #include "fesa/fem/dof_manager.h"
+#include "fesa/model/domain.h"
 
 namespace fesa {
 namespace {
 
-constexpr std::size_t kDofsPerNode = 6U;
-constexpr std::size_t kBeamElementDofCount = 12U;
-constexpr std::size_t kBeamContributionCount =
-    kBeamElementDofCount * kBeamElementDofCount;
-constexpr std::size_t kShellElementDofCount = 24U;
-constexpr std::size_t kShellContributionCount =
-    kShellElementDofCount * kShellElementDofCount;
-
-using BeamElementBuffer = std::array<CooContribution, kBeamContributionCount>;
-using ShellElementBuffer = std::array<CooContribution, kShellContributionCount>;
-
-Result<SparseMatrix> AssemblyFailure(const std::string& code,
-                                     const SourceLocation& location,
-                                     const std::string& identity,
-                                     const std::string& message) {
-  return Result<SparseMatrix>::Failure(Status::Failure(
+Status AssemblyFailure(const std::string& code, const SourceLocation& location,
+                       const std::string& identity,
+                       const std::string& message) {
+  return Status::Failure(
       FailureCategory::kModel,
-      {{Severity::kError, code, location, "*ELEMENT", identity, message}}));
+      {{Severity::kError, code, location, "*ELEMENT", identity, message}});
+}
+
+bool SameSourceIdentity(const SourceEntityId& left,
+                        const SourceEntityId& right) {
+  return left.instance_name == right.instance_name &&
+         left.source_label == right.source_label &&
+         left.source_label_text == right.source_label_text;
+}
+
+bool SameLayout(const ElementDofLayout& left, const ElementDofLayout& right) {
+  return SameSourceIdentity(left.source_id, right.source_id) &&
+         left.node_indices == right.node_indices &&
+         left.components_per_node == right.components_per_node;
 }
 
 }  // namespace
 
-Result<SparseMatrix> SparseAssembler::AssembleStiffness(
-    const AnalysisModel& model, const DofManager& dofs,
+Result<SparseMatrix> SparseAssembler::Assemble(
+    const ElementView& elements, const DofManager& dofs,
     const ParallelFor& parallel_for) {
-  const Domain& domain = model.GetDomain();
-  if (domain.Nodes().size() >
-          (std::numeric_limits<std::size_t>::max)() / kDofsPerNode ||
-      dofs.FullDofCount() != domain.Nodes().size() * kDofsPerNode) {
-    return AssemblyFailure(
-        "invalid-assembly-dimensions", {domain.SourcePath(), 0U},
-        std::to_string(dofs.FullDofCount()),
-        "DofManager dimensions do not match the active model nodes.");
-  }
-  if (!model.ActiveBeamElements().empty() && !domain.ShellElements().Empty()) {
-    return AssemblyFailure(
-        "unsupported-mixed-element-model", {domain.SourcePath(), 0U},
-        "B33:FESA-MITC4",
-        "Sparse assembly does not support mixed beam and shell models.");
-  }
-
-  if (!domain.ShellElements().Empty()) {
-    if (domain.ShellElements().Size() >
-        (std::numeric_limits<std::size_t>::max)() / kShellContributionCount) {
-      return AssemblyFailure(
-          "invalid-assembly-dimensions", {domain.SourcePath(), 0U},
-          std::to_string(domain.ShellElements().Size()),
-          "Shell contribution storage exceeds the addressable range.");
+  std::vector<std::vector<CooContribution>> local_buffers(elements.size());
+  std::vector<std::optional<Status>> local_failures(elements.size());
+  parallel_for.Execute(elements.size(), [&](const std::size_t element_order) {
+    const Element& element = elements[element_order].get();
+    auto stiffness = element.ComputeStiffness();
+    if (!stiffness.HasValue()) {
+      local_failures[element_order] = stiffness.GetStatus();
+      return;
     }
 
-    std::vector<std::optional<std::array<double, 3>>> directors_by_node(
-        domain.Nodes().size());
-    for (const auto& frame : domain.ShellNodeInitialFrames()) {
-      if (frame.node_index >= directors_by_node.size() ||
-          directors_by_node[frame.node_index]) {
-        return AssemblyFailure(
-            "invalid-assembly-element", {domain.SourcePath(), 0U},
-            std::to_string(frame.node_index),
-            "Shell initial frames must map uniquely to model nodes.");
-      }
-      directors_by_node[frame.node_index] = frame.director;
+    const auto& contribution = stiffness.Value();
+    const auto& layout = element.DofLayout();
+    if (!SameLayout(contribution.layout, layout)) {
+      local_failures[element_order] = AssemblyFailure(
+          "invalid-assembly-layout", {}, layout.source_id.source_label_text,
+          "Element stiffness layout must match its declared runtime layout.");
+      return;
+    }
+    auto scatter = dofs.ElementScatter(layout);
+    if (!scatter.HasValue()) {
+      local_failures[element_order] = scatter.GetStatus();
+      return;
+    }
+    const std::size_t local_dof_count = scatter.Value().size();
+    if (contribution.values.Rows() != local_dof_count ||
+        contribution.values.Columns() != local_dof_count ||
+        local_dof_count >
+            (std::numeric_limits<std::size_t>::max)() / local_dof_count) {
+      local_failures[element_order] = AssemblyFailure(
+          "invalid-assembly-dimensions", {}, layout.source_id.source_label_text,
+          "Element stiffness dimensions must match its declared DOF layout.");
+      return;
     }
 
-    struct ShellInput {
-      std::array<const Node*, 4> nodes;
-      std::array<std::array<double, 3>, 4> directors;
-      const ShellSection* section;
-      const LinearElasticMaterial* material;
-      std::array<std::size_t, kShellElementDofCount> scatter;
-    };
-    std::vector<ShellInput> inputs;
-    inputs.reserve(domain.ShellElements().Size());
-    for (std::size_t element_order = 0U;
-         element_order < domain.ShellElements().Size(); ++element_order) {
-      const auto& element = domain.ShellElements()[element_order];
-      if (element.material_index >= domain.LinearElasticMaterials().Size() ||
-          element.section_index >= domain.ShellSections().Size()) {
-        return AssemblyFailure(
-            "invalid-assembly-element", element.location,
-            element.source_id.source_label_text,
-            "Shell element references an entity outside the Domain.");
-      }
-
-      ShellInput input{};
-      input.section = &domain.ShellSections()[element.section_index];
-      input.material = &domain.LinearElasticMaterials()[element.material_index];
-      try {
-        input.scatter =
-            dofs.ShellElementScatter(static_cast<EntityIndex>(element_order));
-      } catch (const std::out_of_range&) {
-        return AssemblyFailure(
-            "invalid-assembly-scatter", element.location,
-            element.source_id.source_label_text,
-            "DofManager does not contain the active shell scatter.");
-      }
-      for (std::size_t node_position = 0U;
-           node_position < element.node_indices.size(); ++node_position) {
-        const EntityIndex node_index = element.node_indices[node_position];
-        if (node_index >= domain.Nodes().size() ||
-            !directors_by_node[node_index]) {
-          return AssemblyFailure(
-              "invalid-assembly-element", element.location,
-              element.source_id.source_label_text,
-              "Shell element requires a valid node and initial director.");
-        }
-        input.nodes[node_position] = &domain.Nodes()[node_index];
-        input.directors[node_position] = *directors_by_node[node_index];
-        for (std::size_t component = 0U; component < kDofsPerNode;
-             ++component) {
-          const std::size_t local = node_position * kDofsPerNode + component;
-          const std::size_t expected =
-              static_cast<std::size_t>(node_index) * kDofsPerNode + component;
-          if (input.scatter[local] != expected ||
-              input.scatter[local] >= dofs.FullDofCount()) {
-            return AssemblyFailure(
-                "invalid-assembly-scatter", element.location,
-                element.source_id.source_label_text,
-                "Shell scatter does not match the active model topology.");
-          }
-        }
-      }
-      inputs.push_back(input);
-    }
-
-    std::vector<ShellElementBuffer> local_buffers(inputs.size());
-    std::vector<std::optional<Status>> local_failures(inputs.size());
-    parallel_for.Execute(inputs.size(), [&](const std::size_t element_order) {
-      const auto& input = inputs[element_order];
-      const auto shell = Mitc4Shell::Create(input.nodes, input.directors,
-                                            *input.section, *input.material);
-      if (!shell.HasValue()) {
-        local_failures[element_order] = shell.GetStatus();
-        return;
-      }
-      const auto stiffness = shell.Value().Stiffness();
-      if (!stiffness.HasValue()) {
-        local_failures[element_order] = stiffness.GetStatus();
-        return;
-      }
-
-      auto& buffer = local_buffers[element_order];
-      for (std::size_t local_row = 0U; local_row < kShellElementDofCount;
-           ++local_row) {
-        for (std::size_t local_column = 0U;
-             local_column < kShellElementDofCount; ++local_column) {
-          const std::size_t local_order =
-              local_row * kShellElementDofCount + local_column;
-          buffer[local_order] = {
-              input.scatter[local_row], input.scatter[local_column],
-              stiffness.Value().stabilized_global24(local_row, local_column),
-              element_order, local_order};
-        }
-      }
-    });
-
-    for (std::size_t element_order = 0U; element_order < local_failures.size();
-         ++element_order) {
-      if (local_failures[element_order]) {
-        return Result<SparseMatrix>::Failure(*local_failures[element_order]);
+    auto& buffer = local_buffers[element_order];
+    buffer.reserve(local_dof_count * local_dof_count);
+    for (std::size_t local_row = 0U; local_row < local_dof_count; ++local_row) {
+      for (std::size_t local_column = 0U; local_column < local_dof_count;
+           ++local_column) {
+        const std::size_t local_order =
+            local_row * local_dof_count + local_column;
+        buffer.push_back({scatter.Value()[local_row],
+                          scatter.Value()[local_column],
+                          contribution.values(local_row, local_column),
+                          element_order, local_order});
       }
     }
-
-    std::vector<CooContribution> contributions;
-    contributions.reserve(local_buffers.size() * kShellContributionCount);
-    // Flatten in source-element order after workers complete. The canonical
-    // COO reduction remains the sole writer of global CSR values.
-    for (const auto& buffer : local_buffers) {
-      contributions.insert(contributions.end(), buffer.begin(), buffer.end());
-    }
-    return SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(),
-                                 std::move(contributions),
-                                 dofs.GetSparsePattern());
-  }
-
-  if (model.ActiveBeamElements().size() >
-      (std::numeric_limits<std::size_t>::max)() / kBeamContributionCount) {
-    return AssemblyFailure(
-        "invalid-assembly-dimensions", {domain.SourcePath(), 0U},
-        std::to_string(model.ActiveBeamElements().size()),
-        "Element contribution storage exceeds the addressable range.");
-  }
-
-  std::vector<std::array<std::size_t, kBeamElementDofCount>> scatters;
-  scatters.reserve(model.ActiveBeamElements().size());
-  for (const EntityIndex element_index : model.ActiveBeamElements()) {
-    if (element_index >= domain.BeamElements().Size()) {
-      return AssemblyFailure("invalid-assembly-element",
-                             {domain.SourcePath(), 0U},
-                             std::to_string(element_index),
-                             "Active element index is outside the Domain.");
-    }
-    const auto& element = domain.BeamElements()[element_index];
-    if (element.node_indices[0U] >= domain.Nodes().size() ||
-        element.node_indices[1U] >= domain.Nodes().size() ||
-        element.material_index >= domain.LinearElasticMaterials().Size() ||
-        element.section_index >= domain.Sections().Size()) {
-      return AssemblyFailure(
-          "invalid-assembly-element", element.location,
-          element.source_id.source_label_text,
-          "Element references an entity outside the Domain.");
-    }
-
-    std::array<std::size_t, kBeamElementDofCount> scatter{};
-    try {
-      scatter = dofs.ElementScatter(element_index);
-    } catch (const std::out_of_range&) {
-      return AssemblyFailure(
-          "invalid-assembly-scatter", element.location,
-          element.source_id.source_label_text,
-          "DofManager does not contain the active element scatter.");
-    }
-    for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) {
-      for (std::size_t component = 0U; component < kDofsPerNode; ++component) {
-        const std::size_t local = endpoint * kDofsPerNode + component;
-        const std::size_t expected =
-            static_cast<std::size_t>(element.node_indices[endpoint]) *
-                kDofsPerNode +
-            component;
-        if (scatter[local] != expected ||
-            scatter[local] >= dofs.FullDofCount()) {
-          return AssemblyFailure(
-              "invalid-assembly-scatter", element.location,
-              element.source_id.source_label_text,
-              "Element scatter does not match the active model topology.");
-        }
-      }
-    }
-    scatters.push_back(scatter);
-  }
-
-  std::vector<BeamElementBuffer> local_buffers(
-      model.ActiveBeamElements().size());
-  std::vector<std::optional<Status>> local_failures(
-      model.ActiveBeamElements().size());
-  parallel_for.Execute(
-      model.ActiveBeamElements().size(), [&](const std::size_t element_order) {
-        const EntityIndex element_index =
-            model.ActiveBeamElements()[element_order];
-        const auto& definition = domain.BeamElements()[element_index];
-        const auto beam = EulerBeam3D::Create(
-            domain.Nodes()[definition.node_indices[0U]],
-            domain.Nodes()[definition.node_indices[1U]],
-            domain.Sections()[definition.section_index],
-            domain.LinearElasticMaterials()[definition.material_index]);
-        if (!beam.HasValue()) {
-          local_failures[element_order] = beam.GetStatus();
-          return;
-        }
-
-        const Matrix stiffness = beam.Value().GlobalStiffness();
-        auto& buffer = local_buffers[element_order];
-        const auto& scatter = scatters[element_order];
-        for (std::size_t local_row = 0U; local_row < kBeamElementDofCount;
-             ++local_row) {
-          for (std::size_t local_column = 0U;
-               local_column < kBeamElementDofCount; ++local_column) {
-            const std::size_t local_order =
-                local_row * kBeamElementDofCount + local_column;
-            buffer[local_order] = {scatter[local_row], scatter[local_column],
-                                   stiffness(local_row, local_column),
-                                   element_order, local_order};
-          }
-        }
-      });
+  });
 
   for (std::size_t element_order = 0U; element_order < local_failures.size();
        ++element_order) {
@@ -289,9 +100,20 @@ Result<SparseMatrix> SparseAssembler::AssembleStiffness(
     }
   }
 
+  std::size_t contribution_count = 0U;
+  for (const auto& buffer : local_buffers) {
+    if (buffer.size() >
+        (std::numeric_limits<std::size_t>::max)() - contribution_count) {
+      return Result<SparseMatrix>::Failure(AssemblyFailure(
+          "invalid-assembly-dimensions", {}, std::to_string(elements.size()),
+          "Element contribution storage exceeds the addressable range."));
+    }
+    contribution_count += buffer.size();
+  }
   std::vector<CooContribution> contributions;
-  contributions.reserve(local_buffers.size() * kBeamContributionCount);
-  // Flatten only after all workers complete; workers never share CSR state.
+  contributions.reserve(contribution_count);
+  // Flatten only after all workers finish. This fixed element order remains
+  // independent of serial, reverse, or TBB task completion order.
   for (const auto& buffer : local_buffers) {
     contributions.insert(contributions.end(), buffer.begin(), buffer.end());
   }
@@ -300,4 +122,31 @@ Result<SparseMatrix> SparseAssembler::AssembleStiffness(
                                dofs.GetSparsePattern());
 }
 
+Result<SparseMatrix> SparseAssembler::AssembleStiffness(
+    const AnalysisModel& model, const DofManager& dofs,
+    const ParallelFor& parallel_for) {
+  const Domain& domain = model.GetDomain();
+  std::vector<std::unique_ptr<Element>> owned_elements;
+  ElementView elements;
+  owned_elements.reserve(model.ActiveElements().size());
+  elements.reserve(model.ActiveElements().size());
+
+  const ElementFactory factory;
+  for (const EntityIndex element_index : model.ActiveElements()) {
+    if (element_index >= domain.Elements().Size()) {
+      return Result<SparseMatrix>::Failure(
+          AssemblyFailure("invalid-assembly-element", {domain.SourcePath(), 0U},
+                          std::to_string(element_index),
+                          "Active element index is outside the Domain."));
+    }
+    auto candidate = factory.Create(domain.Elements()[element_index], domain);
+    if (!candidate.HasValue()) {
+      return Result<SparseMatrix>::Failure(candidate.GetStatus());
+    }
+    owned_elements.push_back(std::move(candidate.Value()));
+    elements.push_back(std::cref(*owned_elements.back()));
+  }
+  return Assemble(elements, dofs, parallel_for);
+}
+
 }  // namespace fesa
diff --git a/tests/unit/assembly/sparse_assembler_test.cpp b/tests/unit/assembly/sparse_assembler_test.cpp
index 5e00050..564ddbd 100644
--- a/tests/unit/assembly/sparse_assembler_test.cpp
+++ b/tests/unit/assembly/sparse_assembler_test.cpp
@@ -6,11 +6,13 @@
 #include <array>
 #include <cstring>
 #include <filesystem>
+#include <functional>
 #include <utility>
 #include <vector>
 
 #include "fesa/analysis/analysis_model.h"
 #include "fesa/assembly/parallel_for.h"
+#include "fesa/elements/element.h"
 #include "fesa/elements/mitc4_shell.h"
 #include "fesa/fem/dof_manager.h"
 #include "fesa/model/domain.h"
@@ -180,6 +182,42 @@ class ReverseParallelFor final : public fesa::ParallelFor {
   mutable std::size_t observed_count_{0U};
 };
 
+class FakeStiffnessElement final : public fesa::Element {
+ public:
+  FakeStiffnessElement(fesa::ElementDofLayout layout, fesa::Matrix stiffness)
+      : layout_{std::move(layout)}, stiffness_{std::move(stiffness)} {}
+
+  const fesa::ElementDofLayout& DofLayout() const noexcept override {
+    return layout_;
+  }
+
+  fesa::Result<fesa::ElementStiffnessContribution> ComputeStiffness()
+      const override {
+    return fesa::Result<fesa::ElementStiffnessContribution>::Success(
+        {layout_, stiffness_});
+  }
+
+  fesa::Result<fesa::ElementResultBundle> Recover(
+      const fesa::Vector&) const override {
+    return fesa::Result<fesa::ElementResultBundle>::Success(
+        {layout_.source_id, fesa::BeamElementResultRows{}});
+  }
+
+ private:
+  fesa::ElementDofLayout layout_;
+  fesa::Matrix stiffness_;
+};
+
+fesa::Matrix MatrixFromRows(const std::vector<std::vector<double>>& rows) {
+  fesa::Matrix matrix{rows.size(), rows.empty() ? 0U : rows.front().size()};
+  for (std::size_t row = 0U; row < rows.size(); ++row) {
+    for (std::size_t column = 0U; column < rows[row].size(); ++column) {
+      matrix(row, column) = rows[row][column];
+    }
+  }
+  return matrix;
+}
+
 void ExpectByteIdentical(const fesa::SparseMatrix& actual,
                          const fesa::SparseMatrix& expected) {
   EXPECT_TRUE(ByteIdentical(actual.RowOffsets(), expected.RowOffsets()));
@@ -187,6 +225,66 @@ void ExpectByteIdentical(const fesa::SparseMatrix& actual,
   EXPECT_TRUE(ByteIdentical(actual.Values(), expected.Values()));
 }
 
+// C-ASSEMBLY-001
+TEST(SparseAssembly,
+     AssemblesFakeRuntimeContributionsIntoExactDeterministicCsr) {
+  auto domain = fesa::Domain::Create(MakeDefinition());
+  ASSERT_TRUE(domain.HasValue());
+  auto model = fesa::AnalysisModel::Create(domain.Value());
+  ASSERT_TRUE(model.HasValue());
+
+  fesa::ElementDofLayout first_layout{
+      {"Beam-1", 10, "10"}, {0U, 2U}, {fesa::DofComponent::kUx}};
+  FakeStiffnessElement first{std::move(first_layout),
+                             MatrixFromRows({{1.0, 2.0}, {3.0, 4.0}})};
+  fesa::ElementDofLayout second_layout{
+      {"Beam-1", 20, "20"},
+      {2U},
+      {fesa::DofComponent::kUx, fesa::DofComponent::kUz,
+       fesa::DofComponent::kUrz}};
+  FakeStiffnessElement second{
+      std::move(second_layout),
+      MatrixFromRows({{5.0, 6.0, 7.0}, {8.0, 9.0, 10.0}, {11.0, 12.0, 13.0}})};
+  const fesa::ElementView elements{std::cref(first), std::cref(second)};
+  fesa::DofManager dofs;
+  ASSERT_TRUE(dofs.Build(model.Value(), elements).IsOk());
+
+  fesa::SerialParallelFor serial_executor;
+  fesa::TbbParallelFor tbb_executor;
+  ReverseParallelFor reverse_executor;
+  auto serial =
+      fesa::SparseAssembler::Assemble(elements, dofs, serial_executor);
+  auto tbb = fesa::SparseAssembler::Assemble(elements, dofs, tbb_executor);
+  auto reversed =
+      fesa::SparseAssembler::Assemble(elements, dofs, reverse_executor);
+  ASSERT_TRUE(serial.HasValue());
+  ASSERT_TRUE(tbb.HasValue());
+  ASSERT_TRUE(reversed.HasValue());
+
+  EXPECT_EQ(reverse_executor.Calls(), 1U);
+  EXPECT_EQ(reverse_executor.ObservedCount(), elements.size());
+  EXPECT_EQ(serial.Value().Rows(), 18U);
+  EXPECT_EQ(serial.Value().Columns(), 18U);
+  EXPECT_EQ(serial.Value().RowOffsets(),
+            (std::vector<std::size_t>{0U, 2U, 2U, 2U, 2U, 2U, 2U, 2U, 2U, 2U,
+                                      2U, 2U, 2U, 6U, 6U, 9U, 9U, 9U, 12U}));
+  EXPECT_EQ(serial.Value().ColumnIndices(),
+            (std::vector<std::size_t>{0U, 12U, 0U, 12U, 14U, 17U, 12U, 14U, 17U,
+                                      12U, 14U, 17U}));
+  EXPECT_EQ(serial.Value().Values(),
+            (std::vector<double>{1.0, 2.0, 3.0, 9.0, 6.0, 7.0, 8.0, 9.0, 10.0,
+                                 11.0, 12.0, 13.0}));
+  ExpectByteIdentical(tbb.Value(), serial.Value());
+  ExpectByteIdentical(reversed.Value(), serial.Value());
+
+  for (std::size_t repetition = 0U; repetition < 8U; ++repetition) {
+    auto repeated =
+        fesa::SparseAssembler::Assemble(elements, dofs, tbb_executor);
+    ASSERT_TRUE(repeated.HasValue());
+    ExpectByteIdentical(repeated.Value(), serial.Value());
+  }
+}
+
 TEST(SparseAssembly, SerialTbbAndRepeatedRunsAreByteIdentical) {
   auto domain_result = fesa::Domain::Create(MakeDefinition());
   ASSERT_TRUE(domain_result.HasValue());