feat(cpp-object-oriented-modular-refactoring): step 14 - runtime-element-factory

include/fesa/elements/element.h

@@ -0,0 +1,85 @@
+#ifndef FESA_ELEMENTS_ELEMENT_H_
+#define FESA_ELEMENTS_ELEMENT_H_
+
+#include <cstdint>
+#include <functional>
+#include <variant>
+#include <vector>
+
+#include "fesa/core/source_identity.h"
+#include "fesa/core/status.h"
+#include "fesa/math/matrix.h"
+#include "fesa/math/vector.h"
+#include "fesa/results/result_records.h"
+
+namespace fesa {
+
+/// @brief Identifies one component in the stable six-DOF node layout.
+enum class DofComponent : std::uint8_t {
+  kUx,
+  kUy,
+  kUz,
+  kUrx,
+  kUry,
+  kUrz,
+};
+
+/// @brief Describes one runtime element's stable node and component order.
+struct ElementDofLayout {
+  SourceEntityId source_id;
+  std::vector<EntityIndex> node_indices;
+  std::vector<DofComponent> components_per_node;
+};
+
+/// @brief Carries one element-local stiffness in its declared DOF order.
+struct ElementStiffnessContribution {
+  ElementDofLayout layout;
+  Matrix values;
+};
+
+/// @brief Keeps beam recovery locations in their distinct row collections.
+struct BeamElementResultRows {
+  std::vector<EndpointResultRow> endpoint_rows;
+  std::vector<GaussResultRow> gauss_rows;
+  std::vector<StressS11Row> stress_rows;
+};
+
+/// @brief Keeps physical shell rows separate from numerical drilling data.
+struct ShellElementResultRows {
+  std::vector<ShellResultRow> rows;
+  double physical_strain_energy{0.0};
+};
+
+/// @brief Selects the physical recovery shape of one runtime element.
+using ElementResultPayload =
+    std::variant<BeamElementResultRows, ShellElementResultRows>;
+
+/// @brief Carries stable source identity with one typed recovery payload.
+struct ElementResultBundle {
+  SourceEntityId source_id;
+  ElementResultPayload payload;
+};
+
+/// @brief Defines the numerical contract consumed by solver pipeline owners.
+class Element {
+ public:
+  virtual ~Element() = default;
+
+  /// @brief Returns the stable element-local DOF ordering.
+  virtual const ElementDofLayout& DofLayout() const noexcept = 0;
+
+  /// @brief Computes the finite stiffness in the declared local ordering.
+  virtual Result<ElementStiffnessContribution> ComputeStiffness() const = 0;
+
+  /// @brief Recovers typed physical rows from element-local global DOFs.
+  virtual Result<ElementResultBundle> Recover(
+      const Vector& element_displacement) const = 0;
+};
+
+/// @brief Provides non-owning runtime elements in stable owner order.
+/// @note Every referenced element must outlive this view.
+using ElementView = std::vector<std::reference_wrapper<const Element>>;
+
+}  // namespace fesa
+
+#endif  // FESA_ELEMENTS_ELEMENT_H_

include/fesa/elements/element_factory.h

@@ -0,0 +1,36 @@
+#ifndef FESA_ELEMENTS_ELEMENT_FACTORY_H_
+#define FESA_ELEMENTS_ELEMENT_FACTORY_H_
+
+#include <memory>
+
+#include "fesa/core/status.h"
+#include "fesa/elements/element.h"
+
+namespace fesa {
+
+class Domain;
+class ElementDefinition;
+
+/// @brief Creates checked numerical elements from Domain-owned definitions.
+class ElementFactory {
+ public:
+  /// @brief Creates the supported runtime kernel for one semantic definition.
+  /// @param definition Definition owned by domain for the returned operation.
+  /// @param domain Immutable owner of referenced nodes, property, and material.
+  /// @return A non-null element or a structured model failure.
+  Result<std::unique_ptr<Element>> Create(const ElementDefinition& definition,
+                                          const Domain& domain) const;
+
+ private:
+  /// @brief Creates one checked B33 runtime candidate.
+  Result<std::unique_ptr<Element>> CreateBeam(
+      const ElementDefinition& definition, const Domain& domain) const;
+
+  /// @brief Creates one checked MITC4 runtime candidate.
+  Result<std::unique_ptr<Element>> CreateShell(
+      const ElementDefinition& definition, const Domain& domain) const;
+};
+
+}  // namespace fesa
+
+#endif  // FESA_ELEMENTS_ELEMENT_FACTORY_H_

include/fesa/elements/euler_beam_3d.h

@@ -8,6 +8,7 @@
 #include <vector>
 
 #include "fesa/core/status.h"
+#include "fesa/elements/element.h"
 #include "fesa/elements/element_definition.h"
 #include "fesa/materials/isotropic_linear_elastic_material.h"
 #include "fesa/math/matrix.h"
@@ -79,7 +80,7 @@ struct BeamRecovery {
 
 /// @brief Implements the approved two-node prismatic B33 Euler beam kernel.
 /// @note Equation numbering and semantic element identity remain external.
-class EulerBeam3D {
+class EulerBeam3D final : public Element {
  public:
   /// @brief Creates a validated beam kernel and right-handed local frame.
   /// @param first_node First source node in the element connectivity.
@@ -92,6 +93,16 @@ class EulerBeam3D {
                                     const GeneralBeamSection& section,
                                     const LinearElasticMaterial& material);
 
+  /// @brief Returns the factory-bound B33 DOF order.
+  const ElementDofLayout& DofLayout() const noexcept override;
+
+  /// @brief Computes global B33 stiffness through the runtime contract.
+  Result<ElementStiffnessContribution> ComputeStiffness() const override;
+
+  /// @brief Recovers typed B33 rows through the runtime contract.
+  Result<ElementResultBundle> Recover(
+      const Vector& element_displacement) const override;
+
   /// @brief Computes the 12-by-12 stiffness in local DOF order.
   /// @note Uses the approved two-point Gauss operation order.
   Matrix LocalStiffness() const;
@@ -107,9 +118,16 @@ class EulerBeam3D {
   /// @brief Recovers signed physical quantities at their distinct locations.
   /// @param global_element_displacement Twelve global element DOF values.
   /// @return Beam recovery rows in deterministic location order.
-  BeamRecovery Recover(const Vector& global_element_displacement) const;
+  BeamRecovery RecoverBeam(const Vector& global_element_displacement) const;
 
  private:
+  friend class ElementFactory;
+
+  /// @brief Binds Domain identity after the numerical candidate is valid.
+  void BindRuntime(ElementDofLayout layout, EntityIndex element_index,
+                   std::vector<SourceEntityId> node_source_ids,
+                   SourceLocation location);
+
   /// @brief Stores already validated geometry, material, and section state.
   EulerBeam3D(double length, double youngs_modulus, double shear_modulus,
               double area, double iy, double iz, double torsional_constant,
@@ -125,6 +143,10 @@ class EulerBeam3D {
   double torsional_constant_;
   std::array<double, 9> rotation_;
   std::vector<std::array<double, 2>> section_points_;
+  ElementDofLayout dof_layout_;
+  EntityIndex element_index_{0U};
+  std::vector<SourceEntityId> node_source_ids_;
+  SourceLocation runtime_location_;
 };
 
 }  // namespace fesa

include/fesa/elements/mitc4_shell.h

@@ -7,6 +7,7 @@
 #include <vector>
 
 #include "fesa/core/status.h"
+#include "fesa/elements/element.h"
 #include "fesa/elements/element_definition.h"
 #include "fesa/materials/isotropic_linear_elastic_material.h"
 #include "fesa/math/matrix.h"
@@ -114,7 +115,7 @@ struct Mitc4PhysicalRecovery {
 
 /// @brief Implements the approved small-rotation FESA-MITC4 shell kernel.
 /// @note Physical and numerical drilling contributions remain separate.
-class Mitc4Shell {
+class Mitc4Shell final : public Element {
  public:
   /// @brief Creates a validated shell kernel from four non-owning node
   /// pointers.
@@ -128,6 +129,16 @@ class Mitc4Shell {
       std::array<std::array<double, 3>, 4> initial_directors,
       const ShellSection& section, const LinearElasticMaterial& material);
 
+  /// @brief Returns the factory-bound MITC4 DOF order.
+  const ElementDofLayout& DofLayout() const noexcept override;
+
+  /// @brief Computes stabilized MITC4 stiffness through the runtime contract.
+  Result<ElementStiffnessContribution> ComputeStiffness() const override;
+
+  /// @brief Recovers physical MITC4 rows through the runtime contract.
+  Result<ElementResultBundle> Recover(
+      const Vector& element_displacement) const override;
+
   /// @brief Evaluates bilinear shape functions and derivatives.
   static Mitc4ShapeFunctions ShapeFunctions(double xi, double eta) noexcept;
 
@@ -183,6 +194,12 @@ class Mitc4Shell {
       const Vector& global_element_displacement24) const;
 
  private:
+  friend class ElementFactory;
+
+  /// @brief Binds Domain identity after the numerical candidate is valid.
+  void BindRuntime(ElementDofLayout layout, EntityIndex element_index,
+                   SourceLocation location);
+
   /// @brief Stores covariant, reciprocal, frame, and Jacobian data at one
   /// point.
   struct GeometryData {
@@ -221,6 +238,9 @@ class Mitc4Shell {
   double poisson_ratio_;
   SourceLocation source_location_;
   std::string identity_;
+  ElementDofLayout dof_layout_;
+  EntityIndex element_index_{0U};
+  SourceLocation runtime_location_;
 };
 
 }  // namespace fesa

include/fesa/fem/dof_manager.h

@@ -3,25 +3,15 @@
 
 #include <array>
 #include <cstddef>
-#include <cstdint>
 #include <optional>
 #include <vector>
 
 #include "fesa/analysis/analysis_model.h"
+#include "fesa/elements/element.h"
 #include "fesa/math/vector.h"
 
 namespace fesa {
 
-/// @brief Identifies one component in the stable six-DOF node layout.
-enum class DofComponent : std::uint8_t {
-  kUx,
-  kUy,
-  kUz,
-  kUrx,
-  kUry,
-  kUrz,
-};
-
 /// @brief Stores the stable structural CSR pattern.
 struct SparsePattern {
   std::vector<std::size_t> row_offsets;

include/fesa/results/result_records.h

@@ -6,7 +6,7 @@
 #include <string>
 #include <vector>
 
-#include "fesa/model/model_types.h"
+#include "fesa/core/source_identity.h"
 
 namespace fesa {
 

src/fesa/CMakeLists.txt

@@ -13,6 +13,7 @@ add_library(
     core/ascii.cpp
     core/diagnostic.cpp
     core/status.cpp
+    elements/element_factory.cpp
     elements/euler_beam_3d.cpp
     elements/mitc4_shell.cpp
     fem/dof_manager.cpp

src/fesa/elements/element_factory.cpp

@@ -0,0 +1,240 @@
+#include "fesa/elements/element_factory.h"
+
+#include <array>
+#include <cstddef>
+#include <memory>
+#include <optional>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/elements/euler_beam_3d.h"
+#include "fesa/elements/mitc4_shell.h"
+#include "fesa/materials/material.h"
+#include "fesa/model/domain.h"
+#include "fesa/properties/element_property.h"
+
+namespace fesa {
+namespace {
+
+/// @brief Creates one deterministic factory failure for a semantic element.
+Result<std::unique_ptr<Element>> FactoryFailure(
+    const ElementDefinition& definition, const Domain& domain, std::string code,
+    std::string message) {
+  return Result<std::unique_ptr<Element>>::Failure(Status::Failure(
+      FailureCategory::kModel, {{Severity::kError,
+                                 std::move(code),
+                                 {domain.SourcePath(), 0U},
+                                 "*ELEMENT",
+                                 definition.SourceId().source_label_text,
+                                 std::move(message)}}));
+}
+
+/// @brief Returns the shared stable six-component node ordering.
+std::vector<DofComponent> FullNodeComponents() {
+  return {DofComponent::kUx,  DofComponent::kUy,  DofComponent::kUz,
+          DofComponent::kUrx, DofComponent::kUry, DofComponent::kUrz};
+}
+
+/// @brief Finds one element's unified stable Domain collection position.
+std::optional<EntityIndex> FindElementIndex(const ElementDefinition& definition,
+                                            const Domain& domain) {
+  for (std::size_t index = 0U; index < domain.Elements().Size(); ++index) {
+    if (&domain.Elements()[index] == &definition) {
+      return static_cast<EntityIndex>(index);
+    }
+  }
+  return std::nullopt;
+}
+
+/// @brief Finds a Domain-owned B33 definition without performing a cast.
+std::optional<EntityIndex> FindBeamIndex(const ElementDefinition& definition,
+                                         const Domain& domain) {
+  for (std::size_t index = 0U; index < domain.BeamElements().Size(); ++index) {
+    const ElementDefinition* candidate = &domain.BeamElements()[index];
+    if (candidate == &definition) {
+      return static_cast<EntityIndex>(index);
+    }
+  }
+  return std::nullopt;
+}
+
+/// @brief Finds a Domain-owned MITC4 definition without performing a cast.
+std::optional<EntityIndex> FindShellIndex(const ElementDefinition& definition,
+                                          const Domain& domain) {
+  for (std::size_t index = 0U; index < domain.ShellElements().Size(); ++index) {
+    const ElementDefinition* candidate = &domain.ShellElements()[index];
+    if (candidate == &definition) {
+      return static_cast<EntityIndex>(index);
+    }
+  }
+  return std::nullopt;
+}
+
+/// @brief Verifies every topology entry before indexing Domain nodes.
+bool HasValidNodes(const ElementDefinition& definition, const Domain& domain,
+                   const std::size_t expected_count) {
+  if (definition.NodeIndices().size() != expected_count) {
+    return false;
+  }
+  for (const EntityIndex node : definition.NodeIndices()) {
+    if (node >= domain.Nodes().size()) {
+      return false;
+    }
+  }
+  return true;
+}
+
+}  // namespace
+
+Result<std::unique_ptr<Element>> ElementFactory::CreateBeam(
+    const ElementDefinition& definition, const Domain& domain) const {
+  const auto element_index = FindElementIndex(definition, domain);
+  const auto beam_index = FindBeamIndex(definition, domain);
+  if (!element_index.has_value() || !beam_index.has_value()) {
+    return FactoryFailure(
+        definition, domain, "invalid-element-definition",
+        "The B33 definition must be owned by the supplied Domain.");
+  }
+  // The concrete access is valid only after the kind and Domain-owned
+  // collection identity checks performed by Create and FindBeamIndex.
+  const auto& beam_definition =
+      static_cast<const EulerBeam3DDefinition&>(definition);
+  if (!HasValidNodes(definition, domain, 2U)) {
+    return FactoryFailure(definition, domain, "invalid-element-topology",
+                          "B33 requires two valid Domain node indices.");
+  }
+  if (definition.PropertyIndex() >= domain.Properties().Size() ||
+      definition.MaterialIndex() >= domain.Materials().Size()) {
+    return FactoryFailure(
+        definition, domain, "invalid-element-assignment",
+        "B33 property and material indices must resolve in the Domain.");
+  }
+  const ElementProperty& property =
+      domain.Properties()[definition.PropertyIndex()];
+  const Material& material = domain.Materials()[definition.MaterialIndex()];
+  if (property.Kind() != ElementPropertyKind::kGeneralBeamSection ||
+      material.Kind() != MaterialKind::kIsotropicLinearElastic) {
+    return FactoryFailure(
+        definition, domain, "incompatible-element-assignment",
+        "B33 requires a general beam section and isotropic elastic material.");
+  }
+  // These casts are centralized after explicit discriminator validation.
+  const auto& section = static_cast<const GeneralBeamSection&>(property);
+  const auto& elastic =
+      static_cast<const IsotropicLinearElasticMaterial&>(material);
+  auto candidate = EulerBeam3D::Create(
+      domain.Nodes()[definition.NodeIndices()[0U]],
+      domain.Nodes()[definition.NodeIndices()[1U]], section, elastic);
+  if (!candidate.HasValue()) {
+    return Result<std::unique_ptr<Element>>::Failure(candidate.GetStatus());
+  }
+
+  ElementDofLayout layout{definition.SourceId(), definition.NodeIndices(),
+                          FullNodeComponents()};
+  std::vector<SourceEntityId> node_source_ids;
+  node_source_ids.reserve(2U);
+  for (const EntityIndex node : definition.NodeIndices()) {
+    node_source_ids.push_back(domain.Nodes()[node].source_id);
+  }
+  auto runtime = std::make_unique<EulerBeam3D>(std::move(candidate.Value()));
+  runtime->BindRuntime(std::move(layout), *element_index,
+                       std::move(node_source_ids), beam_definition.location);
+  std::unique_ptr<Element> result = std::move(runtime);
+  return Result<std::unique_ptr<Element>>::Success(std::move(result));
+}
+
+Result<std::unique_ptr<Element>> ElementFactory::CreateShell(
+    const ElementDefinition& definition, const Domain& domain) const {
+  const auto element_index = FindElementIndex(definition, domain);
+  const auto shell_index = FindShellIndex(definition, domain);
+  if (!element_index.has_value() || !shell_index.has_value()) {
+    return FactoryFailure(
+        definition, domain, "invalid-element-definition",
+        "The MITC4 definition must be owned by the supplied Domain.");
+  }
+  // The concrete access is valid only after the kind and Domain-owned
+  // collection identity checks performed by Create and FindShellIndex.
+  const auto& shell_definition =
+      static_cast<const Mitc4ShellDefinition&>(definition);
+  if (!HasValidNodes(definition, domain, 4U)) {
+    return FactoryFailure(definition, domain, "invalid-element-topology",
+                          "MITC4 requires four valid Domain node indices.");
+  }
+  if (definition.PropertyIndex() >= domain.Properties().Size() ||
+      definition.MaterialIndex() >= domain.Materials().Size()) {
+    return FactoryFailure(
+        definition, domain, "invalid-element-assignment",
+        "MITC4 property and material indices must resolve in the Domain.");
+  }
+  const ElementProperty& property =
+      domain.Properties()[definition.PropertyIndex()];
+  const Material& material = domain.Materials()[definition.MaterialIndex()];
+  if (property.Kind() != ElementPropertyKind::kShellSection ||
+      material.Kind() != MaterialKind::kIsotropicLinearElastic) {
+    return FactoryFailure(
+        definition, domain, "incompatible-element-assignment",
+        "MITC4 requires a shell section and isotropic elastic material.");
+  }
+  const auto& section = static_cast<const ShellSection&>(property);
+  if (section.MaterialIndex() != definition.MaterialIndex()) {
+    return FactoryFailure(
+        definition, domain, "incompatible-element-assignment",
+        "MITC4 definition and shell section must reference one material.");
+  }
+  // These casts are centralized after explicit discriminator validation.
+  const auto& elastic =
+      static_cast<const IsotropicLinearElasticMaterial&>(material);
+
+  std::vector<const ShellNodeInitialFrame*> frame_by_node(domain.Nodes().size(),
+                                                          nullptr);
+  for (const auto& frame : domain.ShellNodeInitialFrames()) {
+    if (frame.node_index >= frame_by_node.size() ||
+        frame_by_node[frame.node_index] != nullptr) {
+      return FactoryFailure(
+          definition, domain, "invalid-shell-frame",
+          "Shell initial frames must map uniquely to valid Domain nodes.");
+    }
+    frame_by_node[frame.node_index] = &frame;
+  }
+
+  std::array<const Node*, 4> nodes{};
+  std::array<std::array<double, 3>, 4> directors{};
+  for (std::size_t position = 0U; position < nodes.size(); ++position) {
+    const EntityIndex node = definition.NodeIndices()[position];
+    const ShellNodeInitialFrame* frame = frame_by_node[node];
+    if (frame == nullptr) {
+      return FactoryFailure(
+          definition, domain, "missing-shell-frame",
+          "MITC4 requires one initial frame for every element node.");
+    }
+    nodes[position] = &domain.Nodes()[node];
+    directors[position] = frame->director;
+  }
+
+  auto candidate = Mitc4Shell::Create(nodes, directors, section, elastic);
+  if (!candidate.HasValue()) {
+    return Result<std::unique_ptr<Element>>::Failure(candidate.GetStatus());
+  }
+  ElementDofLayout layout{definition.SourceId(), definition.NodeIndices(),
+                          FullNodeComponents()};
+  auto runtime = std::make_unique<Mitc4Shell>(std::move(candidate.Value()));
+  runtime->BindRuntime(std::move(layout), *element_index,
+                       shell_definition.location);
+  std::unique_ptr<Element> result = std::move(runtime);
+  return Result<std::unique_ptr<Element>>::Success(std::move(result));
+}
+
+Result<std::unique_ptr<Element>> ElementFactory::Create(
+    const ElementDefinition& definition, const Domain& domain) const {
+  switch (definition.Kind()) {
+    case ElementDefinitionKind::kEulerBeam3D:
+      return CreateBeam(definition, domain);
+    case ElementDefinitionKind::kMitc4Shell:
+      return CreateShell(definition, domain);
+  }
+  return FactoryFailure(definition, domain, "unsupported-element-definition",
+                        "The element definition kind is not supported.");
+}
+
+}  // namespace fesa

src/fesa/elements/euler_beam_3d.cpp

@@ -6,6 +6,7 @@
 #include <limits>
 #include <stdexcept>
 #include <string>
+#include <type_traits>
 #include <utility>
 #include <vector>
 
@@ -77,6 +78,29 @@ Result<EulerBeam3D> ModelFailure(const std::string& code,
       {{Severity::kError, code, location, "*ELEMENT", identity, message}}));
 }
 
+/// @brief Creates a structured failure at the runtime element boundary.
+template <class T>
+Result<T> RuntimeFailure(const SourceLocation& location,
+                         const SourceEntityId& source_id,
+                         const std::string& message) {
+  return Result<T>::Failure(
+      Status::Failure(FailureCategory::kModel,
+                      {{Severity::kError, "invalid-runtime-element", location,
+                        "*ELEMENT", source_id.source_label_text, message}}));
+}
+
+/// @brief Reports whether every fixed-size recovery component is finite.
+template <class T, std::size_t Size>
+bool IsFinite(const std::array<T, Size>& values) {
+  return std::all_of(values.begin(), values.end(), [](const auto& value) {
+    if constexpr (std::is_arithmetic_v<T>) {
+      return std::isfinite(value);
+    } else {
+      return IsFinite(value);
+    }
+  });
+}
+
 /// @brief Repeats the local-axis rotation in stable nodal translation and
 /// rotation blocks.
 Matrix Transformation(const std::array<double, 9>& rotation) {
@@ -428,7 +452,7 @@ Vector EulerBeam3D::LocalEquivalentLoad(
   return equivalent;
 }
 
-BeamRecovery EulerBeam3D::Recover(
+BeamRecovery EulerBeam3D::RecoverBeam(
     const Vector& global_element_displacement) const {
   const Matrix transform = Transformation(rotation_);
   const Vector local_displacement =
@@ -484,6 +508,105 @@ BeamRecovery EulerBeam3D::Recover(
   return recovery;
 }
 
+const ElementDofLayout& EulerBeam3D::DofLayout() const noexcept {
+  return dof_layout_;
+}
+
+Result<ElementStiffnessContribution> EulerBeam3D::ComputeStiffness() const {
+  if (dof_layout_.node_indices.size() != 2U ||
+      dof_layout_.components_per_node.size() != 6U) {
+    return RuntimeFailure<ElementStiffnessContribution>(
+        runtime_location_, dof_layout_.source_id,
+        "B33 runtime stiffness requires a bound two-node six-DOF layout.");
+  }
+  try {
+    return Result<ElementStiffnessContribution>::Success(
+        {dof_layout_, GlobalStiffness()});
+  } catch (const std::exception& error) {
+    return RuntimeFailure<ElementStiffnessContribution>(
+        runtime_location_, dof_layout_.source_id, error.what());
+  }
+}
+
+Result<ElementResultBundle> EulerBeam3D::Recover(
+    const Vector& element_displacement) const {
+  if (dof_layout_.node_indices.size() != 2U ||
+      dof_layout_.components_per_node.size() != 6U ||
+      node_source_ids_.size() != 2U ||
+      element_displacement.Size() != kElementDofCount) {
+    return RuntimeFailure<ElementResultBundle>(
+        runtime_location_, dof_layout_.source_id,
+        "B33 runtime recovery requires a bound two-node layout and 12 DOFs.");
+  }
+  for (std::size_t dof = 0U; dof < element_displacement.Size(); ++dof) {
+    if (!std::isfinite(element_displacement[dof])) {
+      return RuntimeFailure<ElementResultBundle>(
+          runtime_location_, dof_layout_.source_id,
+          "B33 runtime recovery displacement must be finite.");
+    }
+  }
+
+  BeamRecovery recovered{};
+  try {
+    recovered = RecoverBeam(element_displacement);
+  } catch (const std::exception& error) {
+    return RuntimeFailure<ElementResultBundle>(
+        runtime_location_, dof_layout_.source_id, error.what());
+  }
+
+  BeamElementResultRows rows{};
+  rows.endpoint_rows.reserve(2U);
+  rows.gauss_rows.reserve(2U);
+  rows.stress_rows.reserve(recovered.stress_points.size());
+  for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) {
+    if (!IsFinite(recovered.equilibrium_end_actions[endpoint]) ||
+        !IsFinite(recovered.endpoint_section_resultants[endpoint])) {
+      return RuntimeFailure<ElementResultBundle>(
+          runtime_location_, dof_layout_.source_id,
+          "B33 endpoint recovery values must be finite.");
+    }
+    rows.endpoint_rows.push_back(
+        {element_index_, static_cast<int>(endpoint), node_source_ids_[endpoint],
+         recovered.equilibrium_end_actions[endpoint],
+         recovered.endpoint_section_resultants[endpoint]});
+  }
+  for (std::size_t gauss = 0U; gauss < 2U; ++gauss) {
+    if (!IsFinite(recovered.gauss_generalized_strains[gauss]) ||
+        !IsFinite(recovered.gauss_generalized_resultants[gauss])) {
+      return RuntimeFailure<ElementResultBundle>(
+          runtime_location_, dof_layout_.source_id,
+          "B33 Gauss recovery values must be finite.");
+    }
+    rows.gauss_rows.push_back({element_index_, static_cast<int>(gauss + 1U),
+                               recovered.gauss_generalized_strains[gauss],
+                               recovered.gauss_generalized_resultants[gauss]});
+  }
+  for (const auto& point : recovered.stress_points) {
+    if ((point.gauss_point != 1 && point.gauss_point != 2) ||
+        !std::isfinite(point.x1) || !std::isfinite(point.x2) ||
+        !std::isfinite(point.s11)) {
+      return RuntimeFailure<ElementResultBundle>(
+          runtime_location_, dof_layout_.source_id,
+          "B33 stress recovery identity and values must be finite.");
+    }
+    rows.stress_rows.push_back({element_index_, point.gauss_point,
+                                point.section_point, point.x1, point.x2,
+                                point.s11, point.source});
+  }
+  return Result<ElementResultBundle>::Success(
+      {dof_layout_.source_id, std::move(rows)});
+}
+
+void EulerBeam3D::BindRuntime(ElementDofLayout layout,
+                              const EntityIndex element_index,
+                              std::vector<SourceEntityId> node_source_ids,
+                              SourceLocation location) {
+  dof_layout_ = std::move(layout);
+  element_index_ = element_index;
+  node_source_ids_ = std::move(node_source_ids);
+  runtime_location_ = std::move(location);
+}
+
 EulerBeam3D::EulerBeam3D(double length, double youngs_modulus,
                          double shear_modulus, double area, double iy,
                          double iz, double torsional_constant,

src/fesa/elements/mitc4_shell.cpp

@@ -288,6 +288,16 @@ Result<Mitc4PhysicalRecovery> RecoveryFailure(const SourceLocation& location,
                         "*ELEMENT", identity, std::move(message)}}));
 }
 
+/// @brief Creates a structured failure at the runtime element boundary.
+template <class T>
+Result<T> RuntimeFailure(const SourceLocation& location,
+                         const SourceEntityId& source_id, std::string message) {
+  return Result<T>::Failure(Status::Failure(
+      FailureCategory::kModel,
+      {{Severity::kError, "invalid-runtime-element", location, "*ELEMENT",
+        source_id.source_label_text, std::move(message)}}));
+}
+
 }  // namespace
 
 Result<Mitc4Shell> Mitc4Shell::Create(
@@ -794,6 +804,95 @@ Result<Mitc4PhysicalRecovery> Mitc4Shell::RecoverPhysical(
   return Result<Mitc4PhysicalRecovery>::Success(std::move(recovery));
 }
 
+const ElementDofLayout& Mitc4Shell::DofLayout() const noexcept {
+  return dof_layout_;
+}
+
+Result<ElementStiffnessContribution> Mitc4Shell::ComputeStiffness() const {
+  if (dof_layout_.node_indices.size() != kNodeCount ||
+      dof_layout_.components_per_node.size() != kGlobalDofsPerNode) {
+    return RuntimeFailure<ElementStiffnessContribution>(
+        runtime_location_, dof_layout_.source_id,
+        "MITC4 runtime stiffness requires a bound four-node six-DOF layout.");
+  }
+  try {
+    auto stiffness = Stiffness();
+    if (!stiffness.HasValue()) {
+      return Result<ElementStiffnessContribution>::Failure(
+          stiffness.GetStatus());
+    }
+    return Result<ElementStiffnessContribution>::Success(
+        {dof_layout_, std::move(stiffness.Value().stabilized_global24)});
+  } catch (const std::exception& error) {
+    return RuntimeFailure<ElementStiffnessContribution>(
+        runtime_location_, dof_layout_.source_id, error.what());
+  }
+}
+
+Result<ElementResultBundle> Mitc4Shell::Recover(
+    const Vector& element_displacement) const {
+  if (dof_layout_.node_indices.size() != kNodeCount ||
+      dof_layout_.components_per_node.size() != kGlobalDofsPerNode ||
+      element_displacement.Size() != kGlobalDofCount) {
+    return RuntimeFailure<ElementResultBundle>(
+        runtime_location_, dof_layout_.source_id,
+        "MITC4 runtime recovery requires a bound four-node layout and 24 "
+        "DOFs.");
+  }
+
+  Result<Mitc4PhysicalRecovery> recovered = [&]() {
+    try {
+      return RecoverPhysical(element_displacement);
+    } catch (const std::exception& error) {
+      return RuntimeFailure<Mitc4PhysicalRecovery>(
+          runtime_location_, dof_layout_.source_id, error.what());
+    }
+  }();
+  if (!recovered.HasValue()) {
+    return Result<ElementResultBundle>::Failure(recovered.GetStatus());
+  }
+
+  constexpr std::array<ShellMidsurfaceLocation, kNodeCount> locations{
+      ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
+      ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
+  constexpr std::array<ShellSectionPosition, 3> positions{
+      ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle,
+      ShellSectionPosition::kTop};
+  constexpr std::array<double, 3> zeta{-1.0, 0.0, 1.0};
+
+  ShellElementResultRows rows{};
+  rows.rows.reserve(kNodeCount);
+  rows.physical_strain_energy = recovered.Value().strain_energy;
+  for (std::size_t point = 0U; point < recovered.Value().points.size();
+       ++point) {
+    const auto& physical_point = recovered.Value().points[point];
+    ShellResultRow row{};
+    row.element = element_index_;
+    row.location = locations[point];
+    row.natural_coordinates = physical_point.natural_coordinates;
+    row.local_frame = {physical_point.local_frame.e1,
+                       physical_point.local_frame.e2,
+                       physical_point.local_frame.e3};
+    row.generalized_strain = physical_point.generalized_strain;
+    row.section_resultant = physical_point.section_resultant;
+    for (std::size_t position = 0U; position < positions.size(); ++position) {
+      row.stress[position] = {positions[position], zeta[position],
+                              physical_point.in_plane_stress[position]};
+    }
+    rows.rows.push_back(std::move(row));
+  }
+  return Result<ElementResultBundle>::Success(
+      {dof_layout_.source_id, std::move(rows)});
+}
+
+void Mitc4Shell::BindRuntime(ElementDofLayout layout,
+                             const EntityIndex element_index,
+                             SourceLocation location) {
+  dof_layout_ = std::move(layout);
+  element_index_ = element_index;
+  runtime_location_ = std::move(location);
+}
+
 Mitc4Shell::Mitc4Shell(std::array<Vector3, 4> coordinates,
                        std::array<Vector3, 4> directors,
                        std::array<Vector3, 4> tangent_a,

src/fesa/results/result_recovery.cpp

@@ -607,7 +607,8 @@ Status ResultRecovery::Recover(const AnalysisModel& model,
       element_displacement[local_dof] =
           state.Displacement()[scatter[local_dof]];
     }
-    const BeamRecovery recovered = beam.Value().Recover(element_displacement);
+    const BeamRecovery recovered =
+        beam.Value().RecoverBeam(element_displacement);
     for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) {
       if (!IsFinite(recovered.equilibrium_end_actions[endpoint]) ||
           !IsFinite(recovered.endpoint_section_resultants[endpoint])) {

tests/CMakeLists.txt

@@ -13,6 +13,7 @@ add_executable(
     unit/core/diagnostic_test.cpp
     unit/core/source_identity_test.cpp
     unit/core/status_test.cpp
+    unit/elements/element_factory_test.cpp
     unit/elements/euler_beam_3d_test.cpp
     unit/elements/mitc4_shell_test.cpp
     unit/fem/dof_manager_test.cpp

tests/unit/elements/element_factory_test.cpp

@@ -0,0 +1,346 @@
+#include "fesa/elements/element_factory.h"
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <functional>
+#include <memory>
+#include <stdexcept>
+#include <string_view>
+#include <type_traits>
+#include <utility>
+#include <variant>
+#include <vector>
+
+#include "fesa/elements/element.h"
+#include "fesa/elements/euler_beam_3d.h"
+#include "fesa/elements/mitc4_shell.h"
+#include "fesa/fem/dof_manager.h"
+#include "fesa/model/domain.h"
+
+namespace {
+
+fesa::ModelDefinition BeamDefinition(
+    std::array<fesa::EntityIndex, 2> node_indices = {0U, 1U},
+    fesa::EntityIndex material_index = 0U,
+    fesa::EntityIndex property_index = 0U) {
+  fesa::ModelDefinition definition{};
+  definition.source_path = "element-factory-beam.inp";
+  definition.source_content_identity = "fnv1a64:beam";
+  definition.nodes = {
+      {{"Beam-1", 1, "1"}, {0.0, 0.0, 0.0}, {definition.source_path, 1U}},
+      {{"Beam-1", 2, "2"}, {2.0, 0.0, 0.0}, {definition.source_path, 2U}}};
+  definition.materials = {
+      {"Steel", 100.0, 0.25, {definition.source_path, 10U}}};
+  definition.sections = {{"BeamSection",
+                          2.0,
+                          1.0,
+                          0.0,
+                          1.0,
+                          1.0,
+                          {0.0, 1.0, 0.0},
+                          {},
+                          {definition.source_path, 20U}}};
+  definition.elements = {{{"Beam-1", 7, "7"},
+                          node_indices,
+                          material_index,
+                          property_index,
+                          {definition.source_path, 30U}}};
+  return definition;
+}
+
+fesa::ModelDefinition ShellDefinition(
+    fesa::EntityIndex definition_material_index = 0U,
+    fesa::EntityIndex section_material_index = 0U) {
+  fesa::ModelDefinition definition{};
+  definition.source_path = "element-factory-shell.inp";
+  definition.source_content_identity = "fnv1a64:shell";
+  definition.nodes = {
+      {{"Shell-1", 1, "1"}, {-1.0, -1.0, 0.0}, {definition.source_path, 1U}},
+      {{"Shell-1", 2, "2"}, {1.0, -1.0, 0.0}, {definition.source_path, 2U}},
+      {{"Shell-1", 3, "3"}, {1.0, 1.0, 0.0}, {definition.source_path, 3U}},
+      {{"Shell-1", 4, "4"}, {-1.0, 1.0, 0.0}, {definition.source_path, 4U}}};
+  definition.materials = {
+      {"Steel-A", 120.0, 0.25, {definition.source_path, 10U}},
+      {"Steel-B", 200.0, 0.30, {definition.source_path, 11U}}};
+  definition.shell_sections = {{"ShellSection",
+                                0.2,
+                                section_material_index,
+                                {definition.source_path, 20U}}};
+  definition.shell_elements = {{{"Shell-1", 9, "9"},
+                                fesa::ShellSourceElementType::kS4,
+                                {0U, 1U, 2U, 3U},
+                                definition_material_index,
+                                0U,
+                                {definition.source_path, 30U}}};
+  for (fesa::EntityIndex node = 0U; node < 4U; ++node) {
+    definition.shell_node_initial_frames.push_back(
+        {node, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}});
+  }
+  return definition;
+}
+
+fesa::Domain MakeDomain(fesa::ModelDefinition definition) {
+  auto candidate = fesa::Domain::Create(std::move(definition));
+  if (!candidate.HasValue()) {
+    throw std::runtime_error{"Expected a valid Domain fixture."};
+  }
+  return std::move(candidate.Value());
+}
+
+void ExpectSameMatrix(const fesa::Matrix& actual,
+                      const fesa::Matrix& expected) {
+  ASSERT_EQ(actual.Rows(), expected.Rows());
+  ASSERT_EQ(actual.Columns(), expected.Columns());
+  for (std::size_t row = 0U; row < actual.Rows(); ++row) {
+    for (std::size_t column = 0U; column < actual.Columns(); ++column) {
+      EXPECT_DOUBLE_EQ(actual(row, column), expected(row, column));
+    }
+  }
+}
+
+void ExpectModelFailure(
+    const fesa::Result<std::unique_ptr<fesa::Element>>& candidate) {
+  ASSERT_FALSE(candidate.HasValue());
+  EXPECT_EQ(candidate.GetStatus().Category(), fesa::FailureCategory::kModel);
+  ASSERT_EQ(candidate.GetStatus().Diagnostics().size(), 1U);
+  EXPECT_EQ(candidate.GetStatus().Diagnostics()[0].severity,
+            fesa::Severity::kError);
+}
+
+class UnknownElementDefinition final : public fesa::ElementDefinition {
+ public:
+  fesa::ElementDefinitionKind Kind() const noexcept override {
+    return static_cast<fesa::ElementDefinitionKind>(99);
+  }
+
+  const fesa::SourceEntityId& SourceId() const noexcept override {
+    return source_id_;
+  }
+
+  std::string_view SourceElementType() const noexcept override {
+    return "UNKNOWN";
+  }
+
+  const std::vector<fesa::EntityIndex>& NodeIndices() const noexcept override {
+    return node_indices_;
+  }
+
+  fesa::EntityIndex MaterialIndex() const noexcept override { return 0U; }
+
+  fesa::EntityIndex PropertyIndex() const noexcept override { return 0U; }
+
+ private:
+  fesa::SourceEntityId source_id_{"Unknown-1", 99, "99"};
+  std::vector<fesa::EntityIndex> node_indices_{0U, 1U};
+};
+
+class DestructionProbeElement final : public fesa::Element {
+ public:
+  explicit DestructionProbeElement(bool& destroyed) : destroyed_{&destroyed} {}
+
+  ~DestructionProbeElement() override { *destroyed_ = true; }
+
+  const fesa::ElementDofLayout& DofLayout() const noexcept override {
+    return layout_;
+  }
+
+  fesa::Result<fesa::ElementStiffnessContribution> ComputeStiffness()
+      const override {
+    return fesa::Result<fesa::ElementStiffnessContribution>::Success(
+        {layout_, fesa::Matrix{0U, 0U}});
+  }
+
+  fesa::Result<fesa::ElementResultBundle> Recover(
+      const fesa::Vector&) const override {
+    return fesa::Result<fesa::ElementResultBundle>::Success(
+        {layout_.source_id, fesa::BeamElementResultRows{}});
+  }
+
+ private:
+  bool* destroyed_;
+  fesa::ElementDofLayout layout_{};
+};
+
+}  // namespace
+
+// C-ELEMENT-001
+TEST(ElementFactory, CreatesBeamWithStableLayoutStiffnessAndRecovery) {
+  static_assert(std::has_virtual_destructor_v<fesa::Element>);
+
+  bool destroyed = false;
+  {
+    std::unique_ptr<fesa::Element> probe =
+        std::make_unique<DestructionProbeElement>(destroyed);
+    EXPECT_NE(probe, nullptr);
+  }
+  EXPECT_TRUE(destroyed);
+
+  const fesa::Domain domain = MakeDomain(BeamDefinition());
+  const fesa::ElementFactory factory{};
+  auto candidate = factory.Create(domain.Elements()[0U], domain);
+  ASSERT_TRUE(candidate.HasValue());
+  ASSERT_NE(candidate.Value(), nullptr);
+  const std::unique_ptr<fesa::Element>& element = candidate.Value();
+
+  const auto& layout = element->DofLayout();
+  EXPECT_EQ(layout.source_id.source_label_text, "7");
+  EXPECT_EQ(layout.node_indices, (std::vector<fesa::EntityIndex>{0U, 1U}));
+  EXPECT_EQ(layout.components_per_node,
+            (std::vector<fesa::DofComponent>{
+                fesa::DofComponent::kUx, fesa::DofComponent::kUy,
+                fesa::DofComponent::kUz, fesa::DofComponent::kUrx,
+                fesa::DofComponent::kUry, fesa::DofComponent::kUrz}));
+
+  auto stiffness = element->ComputeStiffness();
+  ASSERT_TRUE(stiffness.HasValue());
+  EXPECT_EQ(stiffness.Value().layout.node_indices, layout.node_indices);
+  auto direct = fesa::EulerBeam3D::Create(
+      domain.Nodes()[0U], domain.Nodes()[1U], domain.Sections()[0U],
+      domain.LinearElasticMaterials()[0U]);
+  ASSERT_TRUE(direct.HasValue());
+  ExpectSameMatrix(stiffness.Value().values, direct.Value().GlobalStiffness());
+
+  fesa::Vector displacement{12U};
+  displacement[6U] = 0.1;
+  auto result = element->Recover(displacement);
+  ASSERT_TRUE(result.HasValue());
+  EXPECT_EQ(result.Value().source_id.source_label_text, "7");
+  ASSERT_TRUE(std::holds_alternative<fesa::BeamElementResultRows>(
+      result.Value().payload));
+  const auto& rows =
+      std::get<fesa::BeamElementResultRows>(result.Value().payload);
+  ASSERT_EQ(rows.endpoint_rows.size(), 2U);
+  EXPECT_DOUBLE_EQ(rows.endpoint_rows[0U].end_action[0U], -10.0);
+  EXPECT_DOUBLE_EQ(rows.endpoint_rows[1U].end_action[0U], 10.0);
+  EXPECT_DOUBLE_EQ(rows.endpoint_rows[0U].section_resultant[0U], 10.0);
+  EXPECT_EQ(rows.endpoint_rows[0U].node.source_label_text, "1");
+  EXPECT_EQ(rows.endpoint_rows[1U].node.source_label_text, "2");
+  EXPECT_EQ(rows.gauss_rows.size(), 2U);
+  EXPECT_EQ(rows.stress_rows.size(), 2U);
+}
+
+TEST(ElementFactory, CreatesMitc4WithStableLayoutStiffnessAndRecovery) {
+  auto mixed_definition = ShellDefinition();
+  mixed_definition.sections = {{"BeamSection",
+                                2.0,
+                                1.0,
+                                0.0,
+                                1.0,
+                                1.0,
+                                {0.0, 1.0, 0.0},
+                                {},
+                                {mixed_definition.source_path, 21U}}};
+  mixed_definition.elements = {{{"Shell-1", 8, "8"},
+                                {0U, 1U},
+                                0U,
+                                0U,
+                                {mixed_definition.source_path, 29U}}};
+  const fesa::Domain domain = MakeDomain(std::move(mixed_definition));
+  const fesa::ElementFactory factory{};
+  ASSERT_EQ(domain.Elements().Size(), 2U);
+  auto candidate = factory.Create(domain.Elements()[1U], domain);
+  ASSERT_TRUE(candidate.HasValue());
+  ASSERT_NE(candidate.Value(), nullptr);
+  const std::unique_ptr<fesa::Element>& element = candidate.Value();
+
+  EXPECT_EQ(element->DofLayout().node_indices,
+            (std::vector<fesa::EntityIndex>{0U, 1U, 2U, 3U}));
+  EXPECT_EQ(element->DofLayout().components_per_node.size(), 6U);
+  auto stiffness = element->ComputeStiffness();
+  ASSERT_TRUE(stiffness.HasValue());
+
+  std::array<const fesa::Node*, 4> nodes{
+      &domain.Nodes()[0U], &domain.Nodes()[1U], &domain.Nodes()[2U],
+      &domain.Nodes()[3U]};
+  std::array<std::array<double, 3>, 4> directors{
+      std::array<double, 3>{0.0, 0.0, 1.0},
+      {0.0, 0.0, 1.0},
+      {0.0, 0.0, 1.0},
+      {0.0, 0.0, 1.0}};
+  auto direct =
+      fesa::Mitc4Shell::Create(nodes, directors, domain.ShellSections()[0U],
+                               domain.LinearElasticMaterials()[0U]);
+  ASSERT_TRUE(direct.HasValue());
+  auto direct_stiffness = direct.Value().Stiffness();
+  ASSERT_TRUE(direct_stiffness.HasValue());
+  ExpectSameMatrix(stiffness.Value().values,
+                   direct_stiffness.Value().stabilized_global24);
+
+  auto result = element->Recover(fesa::Vector{24U});
+  ASSERT_TRUE(result.HasValue());
+  ASSERT_TRUE(std::holds_alternative<fesa::ShellElementResultRows>(
+      result.Value().payload));
+  const auto& rows =
+      std::get<fesa::ShellElementResultRows>(result.Value().payload);
+  ASSERT_EQ(rows.rows.size(), 4U);
+  EXPECT_DOUBLE_EQ(rows.physical_strain_energy, 0.0);
+  EXPECT_EQ(rows.rows[0U].location, fesa::ShellMidsurfaceLocation::kGp1);
+  EXPECT_EQ(rows.rows[3U].location, fesa::ShellMidsurfaceLocation::kGp4);
+  for (const auto& row : rows.rows) {
+    EXPECT_EQ(row.element, fesa::EntityIndex{1U});
+  }
+}
+
+TEST(ElementFactory, OwnerBuildsNonOwningViewInStableElementOrder) {
+  const fesa::Domain beam_domain = MakeDomain(BeamDefinition());
+  const fesa::Domain shell_domain = MakeDomain(ShellDefinition());
+  const fesa::ElementFactory factory{};
+  auto beam = factory.Create(beam_domain.Elements()[0U], beam_domain);
+  auto shell = factory.Create(shell_domain.Elements()[0U], shell_domain);
+  ASSERT_TRUE(beam.HasValue());
+  ASSERT_TRUE(shell.HasValue());
+
+  std::vector<std::unique_ptr<fesa::Element>> owner;
+  owner.push_back(std::move(beam.Value()));
+  owner.push_back(std::move(shell.Value()));
+  fesa::ElementView view;
+  view.reserve(owner.size());
+  for (const auto& element : owner) {
+    view.push_back(std::cref(*element));
+  }
+
+  ASSERT_EQ(view.size(), 2U);
+  EXPECT_EQ(&view[0U].get(), owner[0U].get());
+  EXPECT_EQ(&view[1U].get(), owner[1U].get());
+  EXPECT_EQ(view[0U].get().DofLayout().source_id.source_label_text, "7");
+  EXPECT_EQ(view[1U].get().DofLayout().source_id.source_label_text, "9");
+}
+
+TEST(ElementFactory, RejectsUnknownAndIncompatibleDefinitions) {
+  const fesa::ElementFactory factory{};
+  const fesa::Domain beam_domain = MakeDomain(BeamDefinition());
+  const UnknownElementDefinition unknown{};
+  ExpectModelFailure(factory.Create(unknown, beam_domain));
+
+  auto property_mismatch_definition = BeamDefinition();
+  property_mismatch_definition.shell_sections = {
+      {"ShellSection",
+       0.2,
+       0U,
+       {property_mismatch_definition.source_path, 21U}}};
+  property_mismatch_definition.elements[0U].section_index = 1U;
+  const fesa::Domain property_mismatch =
+      MakeDomain(std::move(property_mismatch_definition));
+  ExpectModelFailure(
+      factory.Create(property_mismatch.Elements()[0U], property_mismatch));
+
+  const fesa::Domain material_mismatch = MakeDomain(ShellDefinition(0U, 1U));
+  ExpectModelFailure(
+      factory.Create(material_mismatch.Elements()[0U], material_mismatch));
+
+  const fesa::Domain invalid_node = MakeDomain(BeamDefinition({0U, 2U}));
+  ExpectModelFailure(factory.Create(invalid_node.Elements()[0U], invalid_node));
+
+  auto missing_frame_definition = ShellDefinition();
+  missing_frame_definition.shell_node_initial_frames.pop_back();
+  const fesa::Domain missing_frame =
+      MakeDomain(std::move(missing_frame_definition));
+  auto missing_frame_result =
+      factory.Create(missing_frame.Elements()[0U], missing_frame);
+  ExpectModelFailure(missing_frame_result);
+  EXPECT_EQ(missing_frame_result.GetStatus().Diagnostics()[0U].code,
+            "missing-shell-frame");
+}

tests/unit/elements/euler_beam_3d_test.cpp

@@ -535,7 +535,7 @@ TEST(EulerBeam3D, HermiteAndBMatrixMatchReviewedSigns) {
   displacement[10U] = -slope(w, length);
   displacement[11U] = slope(v, length);
 
-  const BeamRecovery recovery = beam.Recover(displacement);
+  const BeamRecovery recovery = beam.RecoverBeam(displacement);
   const double inverse_sqrt_three = 1.0 / std::sqrt(3.0);
   const std::array<double, 2> gauss_xi = {-inverse_sqrt_three,
                                           inverse_sqrt_three};
@@ -718,7 +718,7 @@ TEST(EulerBeam3D, RotatedTransformPreservesWorkAndEnergy) {
   ExpectScaledNear(global_variation.Dot(global_force),
                    local_variation.Dot(local_force), kMatrixTolerance);
 
-  const BeamRecovery recovery = beam.Recover(global_displacement);
+  const BeamRecovery recovery = beam.RecoverBeam(global_displacement);
   EXPECT_NEAR(recovery.gauss_generalized_strains[0U][0U],
               (local_displacement[6U] - local_displacement[0U]) / 3.0, 1.0e-14);
 }
@@ -738,7 +738,7 @@ TEST(EulerBeam3D, PreservesExactRotatedResultsAcrossVector3Migration) {
   for (std::size_t index = 0U; index < displacement.Size(); ++index) {
     displacement[index] = 0.01 * static_cast<double>(index + 1U) - 0.04;
   }
-  const BeamRecovery recovery = beam.Recover(displacement);
+  const BeamRecovery recovery = beam.RecoverBeam(displacement);
 
   EXPECT_DOUBLE_EQ(global(0U, 0U), 0x1.65f135da12f68p+28);
   EXPECT_DOUBLE_EQ(global(0U, 1U), 0x1.6261c084bda12p+28);
@@ -829,7 +829,7 @@ TEST(EulerBeam3D, AnalyticalAxialTorsionAndTwoPlaneBendingRecover) {
       axial[6U],
       axial_force * length / (material.youngs_modulus * section.area),
       kAnalyticalTolerance);
-  const BeamRecovery axial_recovery = beam.Recover(axial);
+  const BeamRecovery axial_recovery = beam.RecoverBeam(axial);
   ExpectRelativeNear(axial_recovery.equilibrium_end_actions[0U][0U],
                      -axial_force, kAnalyticalTolerance);
   ExpectRelativeNear(axial_recovery.equilibrium_end_actions[1U][0U],
@@ -846,7 +846,7 @@ TEST(EulerBeam3D, AnalyticalAxialTorsionAndTwoPlaneBendingRecover) {
       torsion[9U],
       torque * length / (shear_modulus * section.torsional_constant),
       kAnalyticalTolerance);
-  const BeamRecovery torsion_recovery = beam.Recover(torsion);
+  const BeamRecovery torsion_recovery = beam.RecoverBeam(torsion);
   ExpectRelativeNear(torsion_recovery.equilibrium_end_actions[0U][3U], -torque,
                      kAnalyticalTolerance);
   ExpectRelativeNear(torsion_recovery.equilibrium_end_actions[1U][3U], torque,
@@ -865,7 +865,7 @@ TEST(EulerBeam3D, AnalyticalAxialTorsionAndTwoPlaneBendingRecover) {
                      local_yforce * length * length /
                          (2.0 * material.youngs_modulus * section.i22),
                      kAnalyticalTolerance);
-  const BeamRecovery local_yrecovery = beam.Recover(local_y);
+  const BeamRecovery local_yrecovery = beam.RecoverBeam(local_y);
   ExpectRelativeNear(local_yrecovery.equilibrium_end_actions[0U][1U],
                      -local_yforce, kAnalyticalTolerance);
   ExpectRelativeNear(local_yrecovery.equilibrium_end_actions[1U][1U],
@@ -887,7 +887,7 @@ TEST(EulerBeam3D, AnalyticalAxialTorsionAndTwoPlaneBendingRecover) {
                      -local_zforce * length * length /
                          (2.0 * material.youngs_modulus * section.i11),
                      kAnalyticalTolerance);
-  const BeamRecovery local_zrecovery = beam.Recover(local_z);
+  const BeamRecovery local_zrecovery = beam.RecoverBeam(local_z);
   ExpectRelativeNear(local_zrecovery.equilibrium_end_actions[0U][2U],
                      -local_zforce, kAnalyticalTolerance);
   ExpectRelativeNear(local_zrecovery.equilibrium_end_actions[1U][2U],
@@ -1028,7 +1028,7 @@ TEST(EulerBeam3D, RecoversSectionPointAndDefaultCentroidS11) {
   displacement[10U] = kappa_y * length;
   displacement[11U] = kappa_z * length;
 
-  const BeamRecovery recovery = beam.Recover(displacement);
+  const BeamRecovery recovery = beam.RecoverBeam(displacement);
   ASSERT_EQ(recovery.stress_points.size(), 4U);
   for (std::size_t gauss_point = 0; gauss_point < 2U; ++gauss_point) {
     for (std::size_t point = 0; point < section.section_points.size();
@@ -1050,7 +1050,7 @@ TEST(EulerBeam3D, RecoversSectionPointAndDefaultCentroidS11) {
   }
 
   const auto default_beam = AlignedBeam(length, MakeSection(), material);
-  const BeamRecovery default_recovery = default_beam.Recover(displacement);
+  const BeamRecovery default_recovery = default_beam.RecoverBeam(displacement);
   ASSERT_EQ(default_recovery.stress_points.size(), 2U);
   for (std::size_t gauss_point = 0; gauss_point < 2U; ++gauss_point) {
     const BeamStressPoint& stress = default_recovery.stress_points[gauss_point];
@@ -1091,7 +1091,7 @@ TEST(EulerBeam3D, ReproducesConstantStrainTwistAndCurvaturePatches) {
       shear_modulus * section.torsional_constant * twist,
       material.youngs_modulus * section.i11 * kappa_y,
       material.youngs_modulus * section.i22 * kappa_z};
-  const BeamRecovery recovery = beam.Recover(displacement);
+  const BeamRecovery recovery = beam.RecoverBeam(displacement);
   for (std::size_t point = 0; point < 2U; ++point) {
     for (std::size_t component = 0; component < 4U; ++component) {
       ExpectScaledNear(recovery.gauss_generalized_strains[point][component],