feat(cpp-object-oriented-modular-refactoring): step 18 - load-hierarchy

include/fesa/analysis/analysis_model.h

@@ -21,7 +21,7 @@ class AnalysisModel {
   const Domain& GetDomain() const noexcept;
 
   /// @brief Returns the sole active static step.
-  const StaticStepDefinition& Step() const noexcept;
+  const StepDefinition& Step() const noexcept;
 
   /// @brief Returns active element-definition indices in stable Domain order.
   const std::vector<EntityIndex>& ActiveElements() const noexcept;

include/fesa/assembly/load_assembler.h

@@ -3,6 +3,7 @@
 
 #include "fesa/analysis/analysis_model.h"
 #include "fesa/fem/dof_manager.h"
+#include "fesa/loads/load.h"
 #include "fesa/math/sparse_matrix.h"
 #include "fesa/math/vector.h"
 
@@ -16,6 +17,14 @@ class LoadAssembler {
   static Result<Vector> AssembleFullNodalLoad(const AnalysisModel& model,
                                               const DofManager& dofs);
 
+  /// @brief Accumulates explicitly supplied loads in their view order.
+  /// @param loads Non-owning loads whose contribution source orders must match
+  /// their view positions.
+  /// @return A candidate committed only after all contributions validate.
+  static Result<Vector> AssembleFullNodalLoad(const AnalysisModel& model,
+                                              const DofManager& dofs,
+                                              const LoadView& loads);
+
   /// @brief Forms Ff-Kfc*dc in stable free/constrained order.
   /// @note This operation neither factorizes nor invokes a solver.
   static Result<Vector> EffectiveFreeRhs(const Vector& full_load,

include/fesa/loads/concentrated_nodal_load.h

@@ -0,0 +1,52 @@
+#ifndef FESA_LOADS_CONCENTRATED_NODAL_LOAD_H_
+#define FESA_LOADS_CONCENTRATED_NODAL_LOAD_H_
+
+#include <array>
+#include <cstddef>
+
+#include "fesa/core/diagnostic.h"
+#include "fesa/loads/load.h"
+#include "fesa/model/source_target_resolver.h"
+
+namespace fesa {
+
+/// @brief Emits global concentrated nodal components for one source target.
+class ConcentratedNodalLoad final : public Load {
+ public:
+  /// @brief Creates a six-component global concentrated nodal load.
+  ConcentratedNodalLoad(SourceTargetQuery target,
+                        std::array<double, 6> global_components,
+                        std::size_t source_order);
+
+  /// @brief Creates one parsed CLOAD component while preserving diagnostics.
+  ConcentratedNodalLoad(SourceTargetQuery target, int source_dof,
+                        double magnitude, std::size_t source_order,
+                        SourceLocation location);
+
+  /// @brief Computes target-major, component-minor full-DOF contributions.
+  Result<std::vector<LoadContribution>> ComputeContributions(
+      const LoadContext& context) const override;
+
+  /// @brief Returns the immutable source target query.
+  const SourceTargetQuery& Target() const noexcept;
+
+  /// @brief Returns six global force/moment components without reordering.
+  const std::array<double, 6>& GlobalComponents() const noexcept;
+
+  /// @brief Returns the stable CLOAD declaration order.
+  std::size_t SourceOrder() const noexcept;
+
+  /// @brief Returns the source location used by structured diagnostics.
+  const SourceLocation& Location() const noexcept;
+
+ private:
+  SourceTargetQuery target_;
+  std::array<double, 6> global_components_{};
+  std::size_t source_order_;
+  SourceLocation location_{};
+  int source_dof_{0};
+};
+
+}  // namespace fesa
+
+#endif  // FESA_LOADS_CONCENTRATED_NODAL_LOAD_H_

include/fesa/loads/load.h

@@ -0,0 +1,48 @@
+#ifndef FESA_LOADS_LOAD_H_
+#define FESA_LOADS_LOAD_H_
+
+#include <cstddef>
+#include <functional>
+#include <vector>
+
+#include "fesa/core/status.h"
+
+namespace fesa {
+
+class DofManager;
+class Domain;
+class SourceTargetResolver;
+
+/// @brief Describes one ordered contribution to the full load vector.
+struct LoadContribution {
+  std::size_t source_order;
+  std::size_t full_dof_index;
+  double value;
+};
+
+/// @brief Provides immutable semantic and equation context to a Load.
+/// @note Every referenced object must outlive a contribution request.
+struct LoadContext {
+  const Domain& domain;
+  const DofManager& dof_manager;
+  const SourceTargetResolver& target_resolver;
+};
+
+/// @brief Produces local ordered load contributions without global mutation.
+class Load {
+ public:
+  virtual ~Load() = default;
+
+  /// @brief Computes finite full-DOF contributions in stable target order.
+  /// @param context Non-owning semantic and equation context for this call.
+  /// @return Ordered contributions or a structured model failure.
+  virtual Result<std::vector<LoadContribution>> ComputeContributions(
+      const LoadContext& context) const = 0;
+};
+
+/// @brief Holds non-owning loads in an explicitly supplied source order.
+using LoadView = std::vector<std::reference_wrapper<const Load>>;
+
+}  // namespace fesa
+
+#endif  // FESA_LOADS_LOAD_H_

include/fesa/model/domain.h

@@ -8,10 +8,14 @@
 #include <vector>
 
 #include "fesa/core/status.h"
+#include "fesa/loads/concentrated_nodal_load.h"
+#include "fesa/loads/load.h"
 #include "fesa/model/model_types.h"
 
 namespace fesa {
 
+class StepDefinition;
+
 /// @brief Exposes immutable references without transferring Domain ownership.
 /// @tparam T Base or concrete semantic type stored by the Domain.
 template <class T>
@@ -33,6 +37,7 @@ class DomainCollectionView {
 
  private:
   friend class Domain;
+  friend class StepDefinition;
 
   /// @brief Adds one reference while the owning Domain candidate is built.
   void Add(const T& entry) { entries_.push_back(&entry); }
@@ -40,6 +45,62 @@ class DomainCollectionView {
   std::vector<const T*> entries_;
 };
 
+/// @brief Owns one immutable static-step semantic definition.
+/// @note Loads retain source order and are owned polymorphically by unique
+/// pointers.
+class StepDefinition {
+ public:
+  StepDefinition(const StepDefinition&) = delete;
+  StepDefinition& operator=(const StepDefinition&) = delete;
+  StepDefinition(StepDefinition&&) noexcept = default;
+  StepDefinition& operator=(StepDefinition&&) noexcept = default;
+
+  /// @brief Returns the source step name.
+  const std::string& Name() const noexcept;
+
+  /// @brief Returns current boundary records in declaration order.
+  const std::vector<BoundaryCondition>& Boundaries() const noexcept;
+
+  /// @brief Returns polymorphic loads in stable source order.
+  const LoadView& Loads() const noexcept;
+
+  /// @brief Returns current concentrated loads in stable source order.
+  const DomainCollectionView<ConcentratedNodalLoad>& ConcentratedLoads()
+      const noexcept;
+
+  /// @brief Returns the static initial increment provenance value.
+  double InitialIncrement() const noexcept;
+
+  /// @brief Returns the static time-period provenance value.
+  double TimePeriod() const noexcept;
+
+  /// @brief Returns the static minimum-increment provenance value.
+  double MinimumIncrement() const noexcept;
+
+  /// @brief Returns the static maximum-increment provenance value.
+  double MaximumIncrement() const noexcept;
+
+  /// @brief Returns the source location of the step declaration.
+  const SourceLocation& Location() const noexcept;
+
+ private:
+  friend class Domain;
+
+  /// @brief Converts one parsed static-step record to owned semantic objects.
+  explicit StepDefinition(StaticStepDefinition definition);
+
+  std::string name_;
+  std::vector<BoundaryCondition> boundaries_;
+  std::vector<std::unique_ptr<Load>> loads_;
+  LoadView loads_view_;
+  DomainCollectionView<ConcentratedNodalLoad> concentrated_loads_view_;
+  double initial_increment_;
+  double time_period_;
+  double minimum_increment_;
+  double maximum_increment_;
+  SourceLocation location_;
+};
+
 /// @brief Owns the complete immutable semantic model definition.
 /// @note Collection positions remain stable internal indices after
 /// construction.
@@ -96,7 +157,7 @@ class Domain {
   const std::vector<ElementSet>& ElementSets() const noexcept;
 
   /// @brief Returns static steps in stable declaration order.
-  const std::vector<StaticStepDefinition>& Steps() const noexcept;
+  const DomainCollectionView<StepDefinition>& Steps() const noexcept;
 
   /// @brief Returns sorted nonfatal mapping diagnostics.
   const std::vector<Diagnostic>& Warnings() const noexcept;
@@ -115,6 +176,7 @@ class Domain {
   std::vector<std::unique_ptr<ElementDefinition>> element_definitions_;
   std::vector<std::unique_ptr<ElementProperty>> element_properties_;
   std::vector<std::unique_ptr<Material>> materials_;
+  std::vector<std::unique_ptr<StepDefinition>> step_definitions_;
   DomainCollectionView<ElementDefinition> elements_view_;
   DomainCollectionView<EulerBeam3DDefinition> beam_elements_view_;
   DomainCollectionView<Mitc4ShellDefinition> shell_elements_view_;
@@ -123,6 +185,7 @@ class Domain {
   DomainCollectionView<ShellSection> shell_sections_view_;
   DomainCollectionView<Material> materials_view_;
   DomainCollectionView<LinearElasticMaterial> linear_materials_view_;
+  DomainCollectionView<StepDefinition> steps_view_;
 };
 
 }  // namespace fesa

src/fesa/CMakeLists.txt

@@ -20,6 +20,7 @@ add_library(
     io/abaqus/domain_mapper.cpp
     io/abaqus/input_reader.cpp
     io/hdf5/hdf5_results_writer.cpp
+    loads/concentrated_nodal_load.cpp
     materials/isotropic_linear_elastic_material.cpp
     math/dense_blas_internal.cpp
     math/matrix.cpp

src/fesa/analysis/analysis_model.cpp

@@ -6,7 +6,7 @@
 namespace fesa {
 
 Result<AnalysisModel> AnalysisModel::Create(const Domain& domain) {
-  if (domain.Steps().empty()) {
+  if (domain.Steps().Empty()) {
     return Result<AnalysisModel>::Failure(
         Status::Failure(FailureCategory::kInput,
                         {{Severity::kError,
@@ -16,12 +16,12 @@ Result<AnalysisModel> AnalysisModel::Create(const Domain& domain) {
                           "0",
                           "AnalysisModel requires exactly one static step."}}));
   }
-  if (domain.Steps().size() > 1U) {
+  if (domain.Steps().Size() > 1U) {
     const auto& second_step = domain.Steps()[1];
     return Result<AnalysisModel>::Failure(
         Status::Failure(FailureCategory::kInput,
                         {{Severity::kError, "unsupported-multiple-step",
-                          second_step.location, "STEP", second_step.name,
+                          second_step.Location(), "STEP", second_step.Name(),
                           "AnalysisModel does not support multiple steps."}}));
   }
   return Result<AnalysisModel>::Success(AnalysisModel{domain});
@@ -29,8 +29,8 @@ Result<AnalysisModel> AnalysisModel::Create(const Domain& domain) {
 
 const Domain& AnalysisModel::GetDomain() const noexcept { return *domain_; }
 
-const StaticStepDefinition& AnalysisModel::Step() const noexcept {
-  return domain_->Steps().front();
+const StepDefinition& AnalysisModel::Step() const noexcept {
+  return domain_->Steps()[0U];
 }
 
 const std::vector<EntityIndex>& AnalysisModel::ActiveElements() const noexcept {
@@ -101,10 +101,10 @@ AnalysisModel::AnalysisModel(const Domain& domain) : domain_{&domain} {
     }
   }
 
-  for (std::size_t index = 0U; index < Step().boundaries.size(); ++index) {
+  for (std::size_t index = 0U; index < Step().Boundaries().size(); ++index) {
     active_boundary_conditions_.push_back(static_cast<EntityIndex>(index));
   }
-  for (std::size_t index = 0U; index < Step().loads.size(); ++index) {
+  for (std::size_t index = 0U; index < Step().Loads().size(); ++index) {
     active_loads_.push_back(static_cast<EntityIndex>(index));
   }
 }

src/fesa/assembly/load_assembler.cpp

@@ -9,6 +9,7 @@
 #include <vector>
 
 #include "fesa/constraints/essential_constraints.h"
+#include "fesa/loads/load.h"
 #include "fesa/model/source_target_resolver.h"
 
 namespace fesa {
@@ -25,33 +26,6 @@ Status LoadFailure(const std::string& code, const SourceLocation& location,
       {{Severity::kError, code, location, keyword, identity, message}});
 }
 
-Result<std::vector<EntityIndex>> ResolveTarget(
-    const SourceTargetResolver& resolver, const Domain& domain,
-    const NodalLoad& load) {
-  auto resolved = resolver.Resolve({SourceEntityKind::kNode, "", load.target});
-  if (!resolved.HasValue()) {
-    return Result<std::vector<EntityIndex>>::Failure(
-        LoadFailure("invalid-load-target", load.location, "CLOAD", load.target,
-                    "The load target must resolve unambiguously to one node or "
-                    "one expanded node set."));
-  }
-  std::vector<EntityIndex> nodes;
-  nodes.reserve(resolved.Value().size());
-  std::vector<unsigned char> seen(domain.Nodes().size(), 0U);
-  for (const auto& target : resolved.Value()) {
-    const EntityIndex node = target.entity_index;
-    if (node >= domain.Nodes().size() || seen[node] != 0U) {
-      return Result<std::vector<EntityIndex>>::Failure(LoadFailure(
-          "invalid-load-target", load.location, "CLOAD", load.target,
-          "The expanded node set must contain unique in-range stable node "
-          "identities."));
-    }
-    seen[node] = 1U;
-    nodes.push_back(node);
-  }
-  return Result<std::vector<EntityIndex>>::Success(std::move(nodes));
-}
-
 Status ValidateFiniteVector(const Vector& values,
                             const SourceLocation& location,
                             const std::string& identity) {
@@ -125,6 +99,35 @@ Status ValidateShellMoments(const Domain& domain, const Vector& full_load) {
 
 Result<Vector> LoadAssembler::AssembleFullNodalLoad(const AnalysisModel& model,
                                                     const DofManager& dofs) {
+  const auto& active_loads = model.ActiveLoads();
+  const auto& owned_loads = model.Step().Loads();
+  if (active_loads.size() != owned_loads.size()) {
+    return Result<Vector>::Failure(LoadFailure(
+        "invalid-load-order", model.Step().Location(), "CLOAD",
+        model.Step().Name(),
+        "The active load view must include every sole-step load once."));
+  }
+
+  LoadView loads;
+  loads.reserve(active_loads.size());
+  for (std::size_t source_order = 0U; source_order < active_loads.size();
+       ++source_order) {
+    const EntityIndex load_index = active_loads[source_order];
+    if (static_cast<std::size_t>(load_index) != source_order ||
+        load_index >= owned_loads.size()) {
+      return Result<Vector>::Failure(LoadFailure(
+          "invalid-load-order", model.Step().Location(), "CLOAD",
+          std::to_string(source_order),
+          "Active loads must retain complete stable source order."));
+    }
+    loads.push_back(owned_loads[load_index]);
+  }
+  return AssembleFullNodalLoad(model, dofs, loads);
+}
+
+Result<Vector> LoadAssembler::AssembleFullNodalLoad(const AnalysisModel& model,
+                                                    const DofManager& dofs,
+                                                    const LoadView& loads) {
   const Domain& domain = model.GetDomain();
   if (domain.Nodes().size() >
       (std::numeric_limits<std::size_t>::max)() / kDofsPerNode) {
@@ -165,52 +168,55 @@ Result<Vector> LoadAssembler::AssembleFullNodalLoad(const AnalysisModel& model,
     }
   }
 
-  const auto& active_loads = model.ActiveLoads();
-  const auto& loads = model.Step().loads;
-  if (active_loads.size() != loads.size()) {
-    return Result<Vector>::Failure(LoadFailure(
-        "invalid-load-order", model.Step().location, "CLOAD", model.Step().name,
-        "The active load view must include every sole-step load once."));
-  }
   const SourceTargetIndex target_index = SourceTargetIndex::FromDomain(domain);
   const SourceTargetResolver target_resolver{target_index};
+  const LoadContext context{domain, dofs, target_resolver};
+
+  std::vector<std::vector<LoadContribution>> contributions_by_load;
+  contributions_by_load.reserve(loads.size());
+  for (std::size_t source_order = 0U; source_order < loads.size();
+       ++source_order) {
+    auto contributions =
+        loads[source_order].get().ComputeContributions(context);
+    if (!contributions.HasValue()) {
+      return Result<Vector>::Failure(contributions.GetStatus());
+    }
+    for (const auto& contribution : contributions.Value()) {
+      if (contribution.source_order != source_order) {
+        return Result<Vector>::Failure(LoadFailure(
+            "invalid-load-order", model.Step().Location(), "LOAD_ASSEMBLER",
+            std::to_string(contribution.source_order),
+            "Every contribution must retain its supplying load source "
+            "order."));
+      }
+      if (contribution.full_dof_index >= expected_full_count) {
+        return Result<Vector>::Failure(LoadFailure(
+            "invalid-load-index", model.Step().Location(), "LOAD_ASSEMBLER",
+            std::to_string(contribution.full_dof_index),
+            "A load contribution full-DOF index is outside the active "
+            "model."));
+      }
+      if (!std::isfinite(contribution.value)) {
+        return Result<Vector>::Failure(
+            LoadFailure("nonfinite-load-value", model.Step().Location(),
+                        "LOAD_ASSEMBLER", std::to_string(source_order),
+                        "A load contribution value must be finite."));
+      }
+    }
+    contributions_by_load.push_back(std::move(contributions.Value()));
+  }
 
   Vector full_load{expected_full_count};
-  // Active load indices are required to be the original source order; this
-  // loop is therefore also the fixed floating-point accumulation order.
-  for (std::size_t source_order = 0U; source_order < active_loads.size();
-       ++source_order) {
-    const EntityIndex load_index = active_loads[source_order];
-    if (static_cast<std::size_t>(load_index) != source_order ||
-        load_index >= loads.size()) {
-      return Result<Vector>::Failure(LoadFailure(
-          "invalid-load-order", model.Step().location, "CLOAD",
-          std::to_string(source_order),
-          "Active loads must retain complete stable source order."));
-    }
-    const auto& load = loads[load_index];
-    if (load.dof < 1 || load.dof > static_cast<int>(kDofsPerNode)) {
-      return Result<Vector>::Failure(LoadFailure(
-          "invalid-load-dof", load.location, "CLOAD", load.target,
-          "A nodal load component must be in the range 1 through 6."));
-    }
-    if (!std::isfinite(load.magnitude)) {
-      return Result<Vector>::Failure(
-          LoadFailure("nonfinite-load-value", load.location, "CLOAD",
-                      load.target, "A nodal load magnitude must be finite."));
-    }
-
-    auto target = ResolveTarget(target_resolver, domain, load);
-    if (!target.HasValue()) {
-      return Result<Vector>::Failure(target.GetStatus());
-    }
-    const auto component = static_cast<DofComponent>(load.dof - 1);
-    for (const EntityIndex node : target.Value()) {
-      const std::size_t full_dof = dofs.FullDof(node, component);
-      const double accumulated = full_load[full_dof] + load.magnitude;
+  // Contribution validation completes before this fixed-order candidate
+  // accumulation begins, so failures cannot expose a partial global vector.
+  for (const auto& contributions : contributions_by_load) {
+    for (const auto& contribution : contributions) {
+      const std::size_t full_dof = contribution.full_dof_index;
+      const double accumulated = full_load[full_dof] + contribution.value;
       if (!std::isfinite(accumulated)) {
         return Result<Vector>::Failure(LoadFailure(
-            "nonfinite-load-accumulation", load.location, "CLOAD", load.target,
+            "nonfinite-load-accumulation", model.Step().Location(),
+            "LOAD_ASSEMBLER", std::to_string(full_dof),
             "Source-order load accumulation produced a nonfinite value."));
       }
       full_load[full_dof] = accumulated;

src/fesa/fem/dof_manager.cpp

@@ -180,7 +180,8 @@ Status DofManager::BuildLayouts(const AnalysisModel& analysis_model,
   std::vector<std::optional<double>> prescribed_by_full_dof(full_count);
   for (const EntityIndex boundary_index :
        analysis_model.ActiveBoundaryConditions()) {
-    const auto& boundary = analysis_model.Step().boundaries.at(boundary_index);
+    const auto& boundary =
+        analysis_model.Step().Boundaries().at(boundary_index);
     const auto target = ExpandBoundaryTarget(target_resolver, boundary);
     for (const EntityIndex node : target) {
       for (int component = boundary.first_dof; component <= boundary.last_dof;

src/fesa/loads/concentrated_nodal_load.cpp

@@ -0,0 +1,124 @@
+#include "fesa/loads/concentrated_nodal_load.h"
+
+#include <cmath>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/fem/dof_manager.h"
+#include "fesa/model/domain.h"
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t kDofsPerNode = 6U;
+
+Status LoadFailure(const std::string& code, const SourceLocation& location,
+                   const std::string& identity, const std::string& message) {
+  return Status::Failure(
+      FailureCategory::kModel,
+      {{Severity::kError, code, location, "CLOAD", identity, message}});
+}
+
+}  // namespace
+
+ConcentratedNodalLoad::ConcentratedNodalLoad(
+    SourceTargetQuery target, std::array<double, 6> global_components,
+    const std::size_t source_order)
+    : target_{std::move(target)},
+      global_components_{global_components},
+      source_order_{source_order} {}
+
+ConcentratedNodalLoad::ConcentratedNodalLoad(SourceTargetQuery target,
+                                             const int source_dof,
+                                             const double magnitude,
+                                             const std::size_t source_order,
+                                             SourceLocation location)
+    : target_{std::move(target)},
+      source_order_{source_order},
+      location_{std::move(location)},
+      source_dof_{source_dof} {
+  if (source_dof >= 1 && source_dof <= static_cast<int>(kDofsPerNode)) {
+    global_components_[static_cast<std::size_t>(source_dof - 1)] = magnitude;
+  }
+}
+
+Result<std::vector<LoadContribution>>
+ConcentratedNodalLoad::ComputeContributions(const LoadContext& context) const {
+  const std::string& identity = target_.target_name_or_label;
+  if (source_dof_ != 0 &&
+      (source_dof_ < 1 || source_dof_ > static_cast<int>(kDofsPerNode))) {
+    return Result<std::vector<LoadContribution>>::Failure(LoadFailure(
+        "invalid-load-dof", location_, identity,
+        "A nodal load component must be in the range 1 through 6."));
+  }
+  for (const double component : global_components_) {
+    if (!std::isfinite(component)) {
+      return Result<std::vector<LoadContribution>>::Failure(
+          LoadFailure("nonfinite-load-value", location_, identity,
+                      "A nodal load magnitude must be finite."));
+    }
+  }
+  if (target_.entity_kind != SourceEntityKind::kNode) {
+    return Result<std::vector<LoadContribution>>::Failure(
+        LoadFailure("invalid-load-target", location_, identity,
+                    "A concentrated nodal load requires a node target."));
+  }
+
+  auto resolved = context.target_resolver.Resolve(target_);
+  if (!resolved.HasValue()) {
+    return Result<std::vector<LoadContribution>>::Failure(LoadFailure(
+        "invalid-load-target", location_, identity,
+        "The load target must resolve unambiguously to one node or one "
+        "expanded node set."));
+  }
+
+  std::vector<unsigned char> seen(context.domain.Nodes().size(), 0U);
+  std::vector<LoadContribution> contributions;
+  contributions.reserve(resolved.Value().size() * kDofsPerNode);
+  for (const auto& target : resolved.Value()) {
+    const EntityIndex node = target.entity_index;
+    if (node >= context.domain.Nodes().size() || seen[node] != 0U) {
+      return Result<std::vector<LoadContribution>>::Failure(LoadFailure(
+          "invalid-load-target", location_, identity,
+          "The expanded node set must contain unique in-range stable node "
+          "identities."));
+    }
+    seen[node] = 1U;
+    for (std::size_t component = 0U; component < kDofsPerNode; ++component) {
+      try {
+        contributions.push_back(
+            {source_order_,
+             context.dof_manager.FullDof(node,
+                                         static_cast<DofComponent>(component)),
+             global_components_[component]});
+      } catch (const std::out_of_range&) {
+        return Result<std::vector<LoadContribution>>::Failure(LoadFailure(
+            "invalid-load-target", location_, identity,
+            "The resolved load target must have all six full DOFs."));
+      }
+    }
+  }
+  return Result<std::vector<LoadContribution>>::Success(
+      std::move(contributions));
+}
+
+const SourceTargetQuery& ConcentratedNodalLoad::Target() const noexcept {
+  return target_;
+}
+
+const std::array<double, 6>& ConcentratedNodalLoad::GlobalComponents()
+    const noexcept {
+  return global_components_;
+}
+
+std::size_t ConcentratedNodalLoad::SourceOrder() const noexcept {
+  return source_order_;
+}
+
+const SourceLocation& ConcentratedNodalLoad::Location() const noexcept {
+  return location_;
+}
+
+}  // namespace fesa

src/fesa/model/domain.cpp

@@ -1,11 +1,66 @@
 #include "fesa/model/domain.h"
 
 #include <cstddef>
+#include <functional>
 #include <memory>
 #include <utility>
 
 namespace fesa {
 
+StepDefinition::StepDefinition(StaticStepDefinition definition)
+    : name_{std::move(definition.name)},
+      boundaries_{std::move(definition.boundaries)},
+      initial_increment_{definition.initial_increment},
+      time_period_{definition.time_period},
+      minimum_increment_{definition.minimum_increment},
+      maximum_increment_{definition.maximum_increment},
+      location_{std::move(definition.location)} {
+  loads_.reserve(definition.loads.size());
+  loads_view_.reserve(definition.loads.size());
+  for (std::size_t source_order = 0U; source_order < definition.loads.size();
+       ++source_order) {
+    auto& load = definition.loads[source_order];
+    auto owned = std::make_unique<ConcentratedNodalLoad>(
+        SourceTargetQuery{SourceEntityKind::kNode, "", load.target}, load.dof,
+        load.magnitude, source_order, std::move(load.location));
+    loads_view_.push_back(std::cref(*owned));
+    concentrated_loads_view_.Add(*owned);
+    loads_.push_back(std::move(owned));
+  }
+}
+
+const std::string& StepDefinition::Name() const noexcept { return name_; }
+
+const std::vector<BoundaryCondition>& StepDefinition::Boundaries()
+    const noexcept {
+  return boundaries_;
+}
+
+const LoadView& StepDefinition::Loads() const noexcept { return loads_view_; }
+
+const DomainCollectionView<ConcentratedNodalLoad>&
+StepDefinition::ConcentratedLoads() const noexcept {
+  return concentrated_loads_view_;
+}
+
+double StepDefinition::InitialIncrement() const noexcept {
+  return initial_increment_;
+}
+
+double StepDefinition::TimePeriod() const noexcept { return time_period_; }
+
+double StepDefinition::MinimumIncrement() const noexcept {
+  return minimum_increment_;
+}
+
+double StepDefinition::MaximumIncrement() const noexcept {
+  return maximum_increment_;
+}
+
+const SourceLocation& StepDefinition::Location() const noexcept {
+  return location_;
+}
+
 Result<Domain> Domain::Create(ModelDefinition definition) {
   return Result<Domain>::Success(Domain{std::move(definition)});
 }
@@ -66,8 +121,8 @@ const std::vector<ElementSet>& Domain::ElementSets() const noexcept {
   return definition_.element_sets;
 }
 
-const std::vector<StaticStepDefinition>& Domain::Steps() const noexcept {
-  return definition_.steps;
+const DomainCollectionView<StepDefinition>& Domain::Steps() const noexcept {
+  return steps_view_;
 }
 
 const std::vector<Diagnostic>& Domain::Warnings() const noexcept {
@@ -131,6 +186,15 @@ Domain::Domain(ModelDefinition definition)
     element_definitions_.push_back(std::move(owned));
   }
   definition_.shell_elements.clear();
+
+  step_definitions_.reserve(definition_.steps.size());
+  for (auto& step : definition_.steps) {
+    auto owned =
+        std::unique_ptr<StepDefinition>(new StepDefinition(std::move(step)));
+    steps_view_.Add(*owned);
+    step_definitions_.push_back(std::move(owned));
+  }
+  definition_.steps.clear();
 }
 
 }  // namespace fesa

src/fesa/results/result_recovery.cpp

@@ -15,6 +15,7 @@
 #include <vector>
 
 #include "fesa/elements/element_factory.h"
+#include "fesa/loads/concentrated_nodal_load.h"
 #include "fesa/math/vector3.h"
 #include "fesa/model/source_target_resolver.h"
 
@@ -345,11 +346,12 @@ Status AppendShellRows(const SourceLocation& location,
 
 Result<std::vector<EntityIndex>> ResolveLoadTarget(
     const SourceTargetResolver& resolver, const Domain& domain,
-    const NodalLoad& load) {
-  auto resolved = resolver.Resolve({SourceEntityKind::kNode, "", load.target});
+    const ConcentratedNodalLoad& load) {
+  auto resolved = resolver.Resolve(load.Target());
   if (!resolved.HasValue()) {
     return RecoveryResultFailure<std::vector<EntityIndex>>(
-        "invalid-node-station-entity", load.location, load.target,
+        "invalid-node-station-entity", load.Location(),
+        load.Target().target_name_or_label,
         "A station-eligibility load target must resolve unambiguously.");
   }
   std::vector<EntityIndex> nodes;
@@ -359,7 +361,8 @@ Result<std::vector<EntityIndex>> ResolveLoadTarget(
     const EntityIndex node = target.entity_index;
     if (node >= domain.Nodes().size() || seen[node] != 0U) {
       return RecoveryResultFailure<std::vector<EntityIndex>>(
-          "invalid-node-station-entity", load.location, load.target,
+          "invalid-node-station-entity", load.Location(),
+          load.Target().target_name_or_label,
           "A station-eligibility node set must contain unique valid nodes.");
     }
     seen[node] = 1U;
@@ -789,33 +792,46 @@ ResultRecovery::NormalizeSectionResultantsToNodeStations(
   }
 
   std::vector<unsigned char> loaded_nodes(domain.Nodes().size(), 0U);
-  if (model.ActiveLoads().size() != model.Step().loads.size()) {
+  const auto& loads = model.Step().ConcentratedLoads();
+  if (model.ActiveLoads().size() != loads.Size()) {
     return RecoveryResultFailure<std::vector<NodeStationResultRow>>(
-        "invalid-node-station-entity", model.Step().location, model.Step().name,
+        "invalid-node-station-entity", model.Step().Location(),
+        model.Step().Name(),
         "The active load view must preserve every sole-step load.");
   }
   const SourceTargetIndex target_index = SourceTargetIndex::FromDomain(domain);
   const SourceTargetResolver target_resolver{target_index};
   for (std::size_t order = 0U; order < model.ActiveLoads().size(); ++order) {
     const EntityIndex load_index = model.ActiveLoads()[order];
-    if (load_index != order || load_index >= model.Step().loads.size()) {
+    if (load_index != order || load_index >= loads.Size()) {
       return RecoveryResultFailure<std::vector<NodeStationResultRow>>(
-          "invalid-node-station-entity", model.Step().location,
+          "invalid-node-station-entity", model.Step().Location(),
           std::to_string(load_index),
           "Active loads must remain in stable source order.");
     }
-    const auto& load = model.Step().loads[load_index];
-    if (!std::isfinite(load.magnitude)) {
+    const auto& load = loads[load_index];
+    if (load.SourceOrder() != order) {
       return RecoveryResultFailure<std::vector<NodeStationResultRow>>(
-          "nonfinite-node-station-value", load.location, load.target,
-          "Station eligibility requires finite concentrated loads.");
+          "invalid-node-station-entity", load.Location(),
+          std::to_string(load.SourceOrder()),
+          "Station eligibility requires stable load source order.");
+    }
+    bool has_nonzero_component = false;
+    for (const double component : load.GlobalComponents()) {
+      if (!std::isfinite(component)) {
+        return RecoveryResultFailure<std::vector<NodeStationResultRow>>(
+            "nonfinite-node-station-value", load.Location(),
+            load.Target().target_name_or_label,
+            "Station eligibility requires finite concentrated loads.");
+      }
+      has_nonzero_component = has_nonzero_component || component != 0.0;
     }
     auto targets = ResolveLoadTarget(target_resolver, domain, load);
     if (!targets.HasValue()) {
       return Result<std::vector<NodeStationResultRow>>::Failure(
           targets.GetStatus());
     }
-    if (load.magnitude != 0.0) {
+    if (has_nonzero_component) {
       for (const EntityIndex node : targets.Value()) {
         loaded_nodes[node] = 1U;
       }

tests/CMakeLists.txt

@@ -26,6 +26,7 @@ add_executable(
     unit/io/abaqus/input_reader_test.cpp
     unit/io/abaqus/input_syntax_test.cpp
     unit/io/hdf5/hdf5_results_writer_test.cpp
+    unit/loads/load_test.cpp
     unit/model/domain_test.cpp
     unit/model/model_types_test.cpp
     unit/model/shell_geometry_test.cpp

tests/unit/analysis/analysis_model_test.cpp

@@ -109,6 +109,9 @@ TEST(AnalysisModel, ClassifiesActiveEntitiesInStableOrder) {
   EXPECT_EQ(model.ActiveBoundaryConditions(),
             (std::vector<fesa::EntityIndex>{0U, 1U}));
   EXPECT_EQ(model.ActiveLoads(), (std::vector<fesa::EntityIndex>{0U, 1U, 2U}));
+  ASSERT_EQ(model.Step().Loads().size(), 3U);
+  EXPECT_EQ(&model.Step().Loads()[0U].get(),
+            &domain_result.Value().Steps()[0U].Loads()[0U].get());
 }
 
 TEST(AnalysisModel, ReferencesWithoutCopyingOrMutatingDomain) {
@@ -118,8 +121,9 @@ TEST(AnalysisModel, ReferencesWithoutCopyingOrMutatingDomain) {
   const auto* const element_address = &domain.Elements()[0U];
   const auto* const material_address = &domain.Materials()[0U];
   const auto* const section_address = &domain.Sections()[0U];
-  const std::string step_name = domain.Steps()[0].name;
-  const double first_load_magnitude = domain.Steps()[0].loads[0].magnitude;
+  const std::string step_name = domain.Steps()[0].Name();
+  const double first_load_magnitude =
+      domain.Steps()[0].ConcentratedLoads()[0U].GlobalComponents()[0U];
 
   auto model_result = fesa::AnalysisModel::Create(domain);
   ASSERT_TRUE(model_result.HasValue());
@@ -136,8 +140,10 @@ TEST(AnalysisModel, ReferencesWithoutCopyingOrMutatingDomain) {
             &domain.Materials()[2]);
   EXPECT_EQ(&model.GetDomain().Sections()[model.ActiveSections()[1]],
             &domain.Sections()[1]);
-  EXPECT_EQ(domain.Steps()[0].name, step_name);
-  EXPECT_DOUBLE_EQ(domain.Steps()[0].loads[0].magnitude, first_load_magnitude);
+  EXPECT_EQ(domain.Steps()[0].Name(), step_name);
+  EXPECT_DOUBLE_EQ(
+      domain.Steps()[0].ConcentratedLoads()[0U].GlobalComponents()[0U],
+      first_load_magnitude);
 }
 
 // C-MODEL-002

tests/unit/assembly/load_assembler_test.cpp

@@ -16,6 +16,7 @@
 
 #include "fesa/analysis/analysis_model.h"
 #include "fesa/fem/dof_manager.h"
+#include "fesa/loads/load.h"
 #include "fesa/model/domain.h"
 
 namespace {
@@ -167,8 +168,66 @@ void ExpectFailureCode(const fesa::Result<fesa::Vector>& result,
   EXPECT_EQ(result.GetStatus().Diagnostics()[0U].code, code);
 }
 
+class FakeLoad final : public fesa::Load {
+ public:
+  explicit FakeLoad(std::vector<fesa::LoadContribution> contributions)
+      : contributions_{std::move(contributions)} {}
+
+  fesa::Result<std::vector<fesa::LoadContribution>> ComputeContributions(
+      const fesa::LoadContext&) const override {
+    return fesa::Result<std::vector<fesa::LoadContribution>>::Success(
+        contributions_);
+  }
+
+ private:
+  std::vector<fesa::LoadContribution> contributions_;
+};
+
 }  // namespace
 
+// C-LOAD-001
+TEST(LoadAssembly, AccumulatesGenericContributionsInSuppliedSourceOrder) {
+  auto fixture = MakeFixture(1U, {}, {}, {});
+  const FakeLoad first(
+      {{0U, 0U, 1.0e16}, {0U, 0U, -1.0e16}, {0U, 0U, 1.0}, {0U, 2U, -4.0}});
+  const FakeLoad second({{1U, 2U, 1.5}, {1U, 0U, -2.0}});
+  const fesa::LoadView loads{std::cref(first), std::cref(second)};
+
+  const auto result = fesa::LoadAssembler::AssembleFullNodalLoad(
+      *fixture.model, *fixture.dofs, loads);
+
+  ASSERT_TRUE(result.HasValue());
+  EXPECT_DOUBLE_EQ(result.Value()[0U], -1.0);
+  EXPECT_DOUBLE_EQ(result.Value()[2U], -2.5);
+}
+
+TEST(LoadAssembly, RejectsGenericContributionErrorsBeforeCandidateCommit) {
+  auto fixture = MakeFixture(1U, {}, {}, {});
+  const FakeLoad valid({{0U, 1U, 3.0}});
+  const FakeLoad nonfinite(
+      {{1U, 2U, std::numeric_limits<double>::quiet_NaN()}});
+  const FakeLoad out_of_range({{1U, 6U, 4.0}});
+  const FakeLoad wrong_order({{3U, 2U, 4.0}});
+
+  ExpectFailureCode(fesa::LoadAssembler::AssembleFullNodalLoad(
+                        *fixture.model, *fixture.dofs,
+                        {std::cref(valid), std::cref(nonfinite)}),
+                    "nonfinite-load-value");
+  ExpectFailureCode(fesa::LoadAssembler::AssembleFullNodalLoad(
+                        *fixture.model, *fixture.dofs,
+                        {std::cref(valid), std::cref(out_of_range)}),
+                    "invalid-load-index");
+  ExpectFailureCode(fesa::LoadAssembler::AssembleFullNodalLoad(
+                        *fixture.model, *fixture.dofs,
+                        {std::cref(valid), std::cref(wrong_order)}),
+                    "invalid-load-order");
+
+  const auto repeated = fesa::LoadAssembler::AssembleFullNodalLoad(
+      *fixture.model, *fixture.dofs, {std::cref(valid)});
+  ASSERT_TRUE(repeated.HasValue());
+  EXPECT_DOUBLE_EQ(repeated.Value()[1U], 3.0);
+}
+
 TEST(LoadAssembly, AssemblesNodeSetAndSixComponentLoads) {
   const std::filesystem::path source{"models/load-assembly.inp"};
   auto fixture =

tests/unit/io/abaqus/domain_mapper_test.cpp

@@ -285,23 +285,23 @@ TEST(InpDomainMapping, MapsEverySupportedKeywordAndLegacyDeck) {
   EXPECT_EQ(domain.BeamElements()[0U].material_index, 0U);
   EXPECT_EQ(domain.BeamElements()[0U].section_index, 0U);
 
-  ASSERT_EQ(domain.Steps().size(), 1U);
+  ASSERT_EQ(domain.Steps().Size(), 1U);
   const auto& step = domain.Steps()[0];
-  EXPECT_EQ(step.name, "Step-1");
-  EXPECT_DOUBLE_EQ(step.initial_increment, 0.25);
-  EXPECT_DOUBLE_EQ(step.time_period, 1.5);
-  EXPECT_DOUBLE_EQ(step.minimum_increment, 0.01);
-  EXPECT_DOUBLE_EQ(step.maximum_increment, 1.5);
-  ASSERT_EQ(step.boundaries.size(), 2U);
-  EXPECT_EQ(step.boundaries[0].target, "rootassembly");
-  EXPECT_EQ(step.boundaries[0].first_dof, 1);
-  EXPECT_DOUBLE_EQ(step.boundaries[0].value, 0.125);
-  EXPECT_EQ(step.boundaries[1].last_dof, 2);
-  EXPECT_DOUBLE_EQ(step.boundaries[1].value, 0.0);
-  ASSERT_EQ(step.loads.size(), 1U);
-  EXPECT_EQ(step.loads[0].target, "RootAssembly");
-  EXPECT_EQ(step.loads[0].dof, 6);
-  EXPECT_DOUBLE_EQ(step.loads[0].magnitude, -12.5);
+  EXPECT_EQ(step.Name(), "Step-1");
+  EXPECT_DOUBLE_EQ(step.InitialIncrement(), 0.25);
+  EXPECT_DOUBLE_EQ(step.TimePeriod(), 1.5);
+  EXPECT_DOUBLE_EQ(step.MinimumIncrement(), 0.01);
+  EXPECT_DOUBLE_EQ(step.MaximumIncrement(), 1.5);
+  ASSERT_EQ(step.Boundaries().size(), 2U);
+  EXPECT_EQ(step.Boundaries()[0].target, "rootassembly");
+  EXPECT_EQ(step.Boundaries()[0].first_dof, 1);
+  EXPECT_DOUBLE_EQ(step.Boundaries()[0].value, 0.125);
+  EXPECT_EQ(step.Boundaries()[1].last_dof, 2);
+  EXPECT_DOUBLE_EQ(step.Boundaries()[1].value, 0.0);
+  ASSERT_EQ(step.ConcentratedLoads().Size(), 1U);
+  EXPECT_EQ(step.ConcentratedLoads()[0U].Target().target_name_or_label,
+            "RootAssembly");
+  EXPECT_DOUBLE_EQ(step.ConcentratedLoads()[0U].GlobalComponents()[5U], -12.5);
   EXPECT_EQ(domain.Warnings().size(), 8U);
 
   const auto legacy_path = RepositoryRoot() / "reference" / "cantilever beam" /
@@ -316,7 +316,7 @@ TEST(InpDomainMapping, MapsEverySupportedKeywordAndLegacyDeck) {
   EXPECT_EQ(legacy.Value().Elements().Size(), 10U);
   EXPECT_EQ(legacy.Value().Materials().Size(), 1U);
   EXPECT_EQ(legacy.Value().Sections().Size(), 1U);
-  EXPECT_EQ(legacy.Value().Steps().size(), 1U);
+  EXPECT_EQ(legacy.Value().Steps().Size(), 1U);
   EXPECT_EQ(legacy.Value().Warnings().size(), 7U);
   EXPECT_EQ(ReadExactBytes(legacy_path), bytes_before);
   EXPECT_EQ(std::filesystem::last_write_time(legacy_path), timestamp_before);
@@ -410,17 +410,24 @@ OnlySecond, 2, 5.
   EXPECT_EQ(domain.ElementSets()[2].element_indices,
             (std::vector<fesa::EntityIndex>{1U}));
 
-  ASSERT_EQ(domain.Steps().size(), 1U);
-  EXPECT_EQ(domain.Steps()[0].boundaries[0].target, "OnlySecond");
-  EXPECT_EQ(domain.Steps()[0].loads[0].target, "OnlySecond");
+  ASSERT_EQ(domain.Steps().Size(), 1U);
+  EXPECT_EQ(domain.Steps()[0].Boundaries()[0].target, "OnlySecond");
+  EXPECT_EQ(
+      domain.Steps()[0].ConcentratedLoads()[0U].Target().target_name_or_label,
+      "OnlySecond");
 
   auto direct =
       MapText("direct-node-labels",
               ReplaceOnce(ReplaceOnce(MinimalDeck(), "Root, 1, 6", "1, 1, 6"),
                           "Tip, 2, -1.", "2, 2, -1."));
   ASSERT_TRUE(direct.HasValue());
-  EXPECT_EQ(direct.Value().Steps()[0].boundaries[0].target, "1");
-  EXPECT_EQ(direct.Value().Steps()[0].loads[0].target, "2");
+  EXPECT_EQ(direct.Value().Steps()[0].Boundaries()[0].target, "1");
+  EXPECT_EQ(direct.Value()
+                .Steps()[0]
+                .ConcentratedLoads()[0U]
+                .Target()
+                .target_name_or_label,
+            "2");
 
   auto above_thresholds = MapText(
       "above-geometry-thresholds",
@@ -543,11 +550,11 @@ TEST(InpDomainMapping, NoOpAllowlistWarnsWithoutSemanticEffect) {
             plain.Value().NodeSets().size());
   EXPECT_EQ(with_no_ops.Value().ElementSets().size(),
             plain.Value().ElementSets().size());
-  EXPECT_EQ(with_no_ops.Value().Steps().size(), plain.Value().Steps().size());
-  EXPECT_EQ(with_no_ops.Value().Steps()[0].boundaries.size(),
-            plain.Value().Steps()[0].boundaries.size());
-  EXPECT_EQ(with_no_ops.Value().Steps()[0].loads.size(),
-            plain.Value().Steps()[0].loads.size());
+  EXPECT_EQ(with_no_ops.Value().Steps().Size(), plain.Value().Steps().Size());
+  EXPECT_EQ(with_no_ops.Value().Steps()[0].Boundaries().size(),
+            plain.Value().Steps()[0].Boundaries().size());
+  EXPECT_EQ(with_no_ops.Value().Steps()[0].Loads().size(),
+            plain.Value().Steps()[0].Loads().size());
 }
 
 // MITC4-MAP-001
@@ -605,11 +612,12 @@ TEST(InpDomainMapping, MapsS4AndS4rThroughOneMitc4Identity) {
   EXPECT_EQ(domain.ShellNodeInitialFrames()[0].tangent_b,
             (std::array<double, 3>{0.0, 1.0, 0.0}));
 
-  ASSERT_EQ(domain.Steps().size(), 1U);
-  ASSERT_EQ(domain.Steps()[0].boundaries.size(), 1U);
-  EXPECT_EQ(domain.Steps()[0].boundaries[0].last_dof, 6);
-  ASSERT_EQ(domain.Steps()[0].loads.size(), 1U);
-  EXPECT_EQ(domain.Steps()[0].loads[0].dof, 6);
+  ASSERT_EQ(domain.Steps().Size(), 1U);
+  ASSERT_EQ(domain.Steps()[0].Boundaries().size(), 1U);
+  EXPECT_EQ(domain.Steps()[0].Boundaries()[0].last_dof, 6);
+  ASSERT_EQ(domain.Steps()[0].ConcentratedLoads().Size(), 1U);
+  EXPECT_DOUBLE_EQ(
+      domain.Steps()[0].ConcentratedLoads()[0U].GlobalComponents()[5U], 1.0);
 }
 
 // MITC4-MAP-002

tests/unit/loads/load_test.cpp

@@ -0,0 +1,107 @@
+#include <gtest/gtest.h>
+
+#include <array>
+#include <cstddef>
+#include <filesystem>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+#include "fesa/analysis/analysis_model.h"
+#include "fesa/fem/dof_manager.h"
+#include "fesa/loads/concentrated_nodal_load.h"
+#include "fesa/model/domain.h"
+#include "fesa/model/source_target_resolver.h"
+
+namespace {
+
+fesa::Domain MakeDomain() {
+  const std::filesystem::path source{"models/load-hierarchy.inp"};
+  fesa::ModelDefinition definition{};
+  definition.source_path = source;
+  definition.source_content_identity = "fnv1a64:loadhierarchy";
+  definition.nodes = {{{"Part-1-1", 10, "10"}, {0.0, 0.0, 0.0}, {source, 2U}},
+                      {{"Part-1-1", 20, "20"}, {1.0, 0.0, 0.0}, {source, 3U}}};
+  definition.node_sets = {
+      {"Pair", std::string{"Part-1-1"}, {0U, 1U}, {source, 4U}}};
+  definition.steps = {{"Step-1", {}, {}, 0.1, 1.0, 0.01, 1.0, {source, 5U}}};
+  auto domain = fesa::Domain::Create(std::move(definition));
+  if (!domain.HasValue()) {
+    throw std::runtime_error{"Load hierarchy Domain construction failed."};
+  }
+  return std::move(domain.Value());
+}
+
+void ExpectFailureCode(
+    const fesa::Result<std::vector<fesa::LoadContribution>>& result,
+    const std::string& code) {
+  ASSERT_FALSE(result.HasValue());
+  EXPECT_EQ(result.GetStatus().Category(), fesa::FailureCategory::kModel);
+  ASSERT_EQ(result.GetStatus().Diagnostics().size(), 1U);
+  EXPECT_EQ(result.GetStatus().Diagnostics()[0U].code, code);
+}
+
+}  // namespace
+
+// C-LOAD-001
+TEST(Load, ConcentratedNodalLoadEmitsStableFullDofContributions) {
+  static_assert(std::has_virtual_destructor_v<fesa::Load>);
+  static_assert(!std::is_copy_constructible_v<fesa::StepDefinition>);
+  static_assert(std::is_move_constructible_v<fesa::StepDefinition>);
+
+  fesa::Domain domain = MakeDomain();
+  auto model = fesa::AnalysisModel::Create(domain);
+  ASSERT_TRUE(model.HasValue());
+  auto dofs = fesa::DofManager::Create(model.Value());
+  ASSERT_TRUE(dofs.HasValue());
+  const fesa::SourceTargetIndex target_index =
+      fesa::SourceTargetIndex::FromDomain(domain);
+  const fesa::SourceTargetResolver resolver{target_index};
+  const fesa::LoadContext context{domain, dofs.Value(), resolver};
+  const std::array<double, 6> components{1.0, -2.0, 0.0, 4.0, 0.0, -6.0};
+  const fesa::ConcentratedNodalLoad load(
+      {fesa::SourceEntityKind::kNode, "Part-1-1", "pair"}, components, 7U);
+
+  const auto result = load.ComputeContributions(context);
+
+  ASSERT_TRUE(result.HasValue());
+  ASSERT_EQ(result.Value().size(), 12U);
+  for (std::size_t node = 0U; node < 2U; ++node) {
+    for (std::size_t component = 0U; component < components.size();
+         ++component) {
+      const auto& contribution =
+          result.Value()[node * components.size() + component];
+      EXPECT_EQ(contribution.source_order, 7U);
+      EXPECT_EQ(contribution.full_dof_index,
+                node * components.size() + component);
+      EXPECT_DOUBLE_EQ(contribution.value, components[component]);
+    }
+  }
+}
+
+TEST(Load, ConcentratedNodalLoadRejectsInvalidInputWithoutContributions) {
+  fesa::Domain domain = MakeDomain();
+  auto model = fesa::AnalysisModel::Create(domain);
+  ASSERT_TRUE(model.HasValue());
+  auto dofs = fesa::DofManager::Create(model.Value());
+  ASSERT_TRUE(dofs.HasValue());
+  const fesa::SourceTargetIndex target_index =
+      fesa::SourceTargetIndex::FromDomain(domain);
+  const fesa::SourceTargetResolver resolver{target_index};
+  const fesa::LoadContext context{domain, dofs.Value(), resolver};
+
+  std::array<double, 6> nonfinite{};
+  nonfinite[2U] = std::numeric_limits<double>::infinity();
+  const fesa::ConcentratedNodalLoad invalid_value(
+      {fesa::SourceEntityKind::kNode, "Part-1-1", "10"}, nonfinite, 0U);
+  ExpectFailureCode(invalid_value.ComputeContributions(context),
+                    "nonfinite-load-value");
+
+  const fesa::ConcentratedNodalLoad invalid_target(
+      {fesa::SourceEntityKind::kNode, "Part-1-1", "Missing"}, {}, 0U);
+  ExpectFailureCode(invalid_target.ComputeContributions(context),
+                    "invalid-load-target");
+}

tests/unit/model/domain_test.cpp

@@ -191,15 +191,18 @@ TEST(DomainModel, ImmutableOwnershipPreservesStableOrder) {
   ASSERT_EQ(domain.ElementSets().size(), 1U);
   EXPECT_EQ(domain.ElementSets()[0].element_indices[0], 0U);
 
-  ASSERT_EQ(domain.Steps().size(), 1U);
-  EXPECT_EQ(domain.Steps()[0].name, "Step-1");
-  EXPECT_DOUBLE_EQ(domain.Steps()[0].initial_increment, 0.1);
-  EXPECT_DOUBLE_EQ(domain.Steps()[0].time_period, 1.0);
-  EXPECT_DOUBLE_EQ(domain.Steps()[0].minimum_increment, 1.0e-5);
-  EXPECT_DOUBLE_EQ(domain.Steps()[0].maximum_increment, 1.0);
-  ASSERT_EQ(domain.Steps()[0].boundaries.size(), 1U);
-  ASSERT_EQ(domain.Steps()[0].loads.size(), 1U);
-  EXPECT_DOUBLE_EQ(domain.Steps()[0].loads[0].magnitude, -1000.0);
+  ASSERT_EQ(domain.Steps().Size(), 1U);
+  EXPECT_EQ(domain.Steps()[0].Name(), "Step-1");
+  EXPECT_DOUBLE_EQ(domain.Steps()[0].InitialIncrement(), 0.1);
+  EXPECT_DOUBLE_EQ(domain.Steps()[0].TimePeriod(), 1.0);
+  EXPECT_DOUBLE_EQ(domain.Steps()[0].MinimumIncrement(), 1.0e-5);
+  EXPECT_DOUBLE_EQ(domain.Steps()[0].MaximumIncrement(), 1.0);
+  ASSERT_EQ(domain.Steps()[0].Boundaries().size(), 1U);
+  ASSERT_EQ(domain.Steps()[0].Loads().size(), 1U);
+  ASSERT_EQ(domain.Steps()[0].ConcentratedLoads().Size(), 1U);
+  EXPECT_DOUBLE_EQ(
+      domain.Steps()[0].ConcentratedLoads()[0U].GlobalComponents()[2U],
+      -1000.0);
 
   ASSERT_EQ(domain.Warnings().size(), 1U);
   EXPECT_EQ(domain.Warnings()[0].code, "ignored-output-request");