feat(linear-static-3d-euler-beam): step 16 - euler-beam-element

docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md

@@ -513,3 +513,70 @@
 - handoff: the exact backend-free result records and `AnalysisState` ledger API
   provide full-space mutable solution vectors, step/frame identity, and stable
   recovery-row collections to Step 16 and later result recovery/output tasks.
+
+## Step 16 — euler-beam-element
+
+- task_id: `TASK-16`
+- status: `completed`
+- changed_files: `include/fesa/elements/euler_beam_3d.hpp`,
+  `src/fesa/elements/euler_beam_3d.cpp`,
+  `tests/unit/elements/euler_beam_3d_test.cpp`, `src/fesa/CMakeLists.txt`,
+  `tests/CMakeLists.txt`,
+  `docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md`,
+  `phases/linear-static-3d-euler-beam/index.json`,
+  `.superpowers/sdd/linear-static-3d-euler-beam/task-16-report.md`
+- requirement_ids: `FESA-REQ-LS3DEB-002`, `FESA-REQ-LS3DEB-004`,
+  `FESA-REQ-LS3DEB-008`, `FESA-REQ-LS3DEB-012`,
+  `FESA-REQ-LS3DEB-013`, `FESA-REQ-LS3DEB-014`,
+  `FESA-REQ-LS3DEB-031`, `FESA-REQ-LS3DEB-032`,
+  `FESA-REQ-LS3DEB-035`, `FESA-REQ-LS3DEB-044`
+- test_ids: `T16-BEAM-001`, `T16-BEAM-002`, `T16-BEAM-003`,
+  `T16-BEAM-004`, `T16-BEAM-005`, `T16-BEAM-006`, `T16-BEAM-007`,
+  `T16-BEAM-008`, `T16-BEAM-009`, `T16-BEAM-010`
+
+| stage | exact command | exit_code | expected_or_observed_result | evidence_tail |
+| --- | --- | ---: | --- | --- |
+| RED-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 1 | Exactly ten planned tests were registered before production and the EulerBeam3D public API was absent | MSVC C1083 reported missing `fesa/elements/euler_beam_3d.hpp` from `euler_beam_3d_test.cpp` after successful CMake regeneration |
+| GREEN-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 0 | The minimum beam kernel, ten tests, solver library, and unit executable compile and link | `euler_beam_3d.cpp`, `euler_beam_3d_test.cpp`, `fesa_solver.lib`, and `fesa_unit_tests.exe` built without a FESA warning under `/W4 /WX` |
+| GREEN-test | `ctest --test-dir .harness/build -C Debug -R EulerBeam3D --output-on-failure` | 0 | Reviewed signs, Gauss integration, rank, transform, loads, analytical modes, validation, recovery, patches, and negative control pass | 10/10 exact `EulerBeam3D` tests passed |
+| VERIFY-configure | `cmake -S . -B .harness/build -A x64 -DFESA_GTEST_SOURCE_DIR=C:/git/googletest "-DMKL_DIR=C:/Program Files (x86)/Intel/oneAPI/mkl/2026.1/lib/cmake/mkl" "-DTBB_DIR=C:/Program Files (x86)/Intel/oneAPI/tbb/2023.1/lib/cmake/tbb" "-DHDF5_DIR=C:/Program Files/HDF_Group/HDF5/2.1.1/cmake"` | 0 | Approved explicit-dependency MSVC x64 build tree generates | Windows SDK, oneMKL 2026.1, oneTBB, and HDF5 resolved; configure and generate completed |
+| VERIFY-build | `cmake --build .harness/build --config Debug` | 0 | Full Debug build passes without a new FESA warning | `fesa_solver.lib` and `fesa_unit_tests.exe` built under `/W4 /WX` |
+| VERIFY-targeted | `ctest --test-dir .harness/build -C Debug -R EulerBeam3D --output-on-failure` | 0 | Focused Step 16 suite remains green | 10/10 exact `EulerBeam3D` tests passed |
+| VERIFY-discovery | `ctest --test-dir .harness/build -C Debug --show-only=json-v1` | 0 | CTest discovers the accumulated suite and all ten exact EulerBeam3D names | 39 tests discovered, including 10 `EulerBeam3D` tests, with feature/unit labels |
+| VERIFY-full | `ctest --test-dir .harness/build -C Debug --output-on-failure` | 0 | Full accumulated C++ suite has zero failures | 39/39 tests passed |
+| VERIFY-contract-scans | Backend public-header, upward-dependency, out-of-scope formulation, production one-point path, two-point Gauss/invariant, exact-test-count, and CMake-registration scans using fail-on-match `rg` wrappers | 0 | The kernel implements only the approved exact ledger and numerical integration contract | backend leaks 0; upward dependencies 0; out-of-scope couplings 0; production one-point paths 0; two-point Gauss consumers 3; normalized closed-form invariant `1e-12`; tests 10; registrations 1/1 |
+| VERIFY-diff | `git diff --check` plus trailing-whitespace scan over the three new files | 0 | Tracked and untracked Step 16 files have no whitespace errors | Diff check exit 0; new-file trailing whitespace matches 0 |
+| VERIFY-reference | `git diff --exit-code -- reference/`; `git status --short -- reference/` | 0 | Approved legacy reference artifacts remain read-only and unchanged | Reference diff exit 0 and reference status empty |
+
+- green_triage: the first implementation build passed and the first focused
+  run was 9/10. The single failure was isolated to a large-coordinate test
+  fixture whose requested equality rounded to a represented length above the
+  threshold. The test fixture was replaced by an exactly representable unit
+  threshold equality plus below/above large-coordinate cases; the production
+  strict `>` comparison and formulation were unchanged before the 10/10 run.
+- contract_checks: local DOF order is exactly
+  `[u1,v1,w1,rx1,ry1,rz1,u2,v2,w2,rx2,ry2,rz2]`. Rotation rows are
+  `(ex,ey,ez)`, `dl=T*dg`, and `Kg=T^T*Kl*T`. The B-matrix follows the reviewed
+  `theta_y=-w'`, `theta_z=v'` signs; two-point Gauss `B^T D B` is the sole
+  production stiffness path and is checked against the independent closed
+  matrix with normalized tolerance `1e-12`. The test-only one-point rule has
+  rank four while the production rule has rank six and six rigid modes.
+- contract_checks: length and guide-vector tests use the exact scale-aware
+  strict thresholds with no fallback axis. `E`, derived `G`, `A`, `Iy`, `Iz`,
+  and `J` must be finite and positive and `I12` must be exactly zero. Constant
+  local `px,py,pz,mx` receives the consistent signed load only; no parser,
+  `*DLOAD`, B31, Timoshenko, shear correction, or transverse/torsional stress
+  support was added.
+- contract_checks: recovery keeps local outward equilibrium end actions,
+  endpoint section resultants, and two-Gauss generalized strain/resultants
+  distinct. Axial stress is
+  `S11=E(epsilon0+x2*kappa_y-x1*kappa_z)` in input section-point order; absent
+  input points produce section point 0 at the centroid with source
+  `fesa-default`.
+- generated_evidence: `.harness/build/src/fesa/Debug/fesa_solver.lib`,
+  `.harness/build/tests/Debug/fesa_unit_tests.exe`
+- reference_diff: unchanged; `git diff --exit-code -- reference/` exit 0
+- handoff: the exact backend-free `EulerBeam3D` ledger API supplies global and
+  local stiffness, constant local equivalent load, and distinct recovery data
+  to Step 18 sparse assembly and Step 22 result recovery without owning element
+  identity, equation numbering, parser coupling, or result persistence.

include/fesa/elements/euler_beam_3d.hpp

@@ -0,0 +1,74 @@
+#pragma once
+
+#include "fesa/core/status.hpp"
+#include "fesa/math/matrix.hpp"
+#include "fesa/math/vector.hpp"
+#include "fesa/model/model_types.hpp"
+
+#include <array>
+#include <cstddef>
+#include <string>
+#include <vector>
+
+namespace fesa {
+
+struct ConstantLocalLineLoad {
+    double px;
+    double py;
+    double pz;
+    double mx;
+};
+
+struct BeamStressPoint {
+    int gaussPoint;
+    std::size_t sectionPoint;
+    double x1;
+    double x2;
+    double s11;
+    std::string source;
+};
+
+struct BeamRecovery {
+    std::array<std::array<double, 6>, 2> equilibriumEndActions;
+    std::array<std::array<double, 4>, 2> endpointSectionResultants;
+    std::array<std::array<double, 4>, 2> gaussGeneralizedStrains;
+    std::array<std::array<double, 4>, 2> gaussGeneralizedResultants;
+    std::vector<BeamStressPoint> stressPoints;
+};
+
+// Implements the approved two-node straight prismatic B33 Euler-Bernoulli
+// kernel. Equation numbering and element identity remain outside this type.
+class EulerBeam3D {
+public:
+    static Result<EulerBeam3D> create(const Node& firstNode,
+                                      const Node& secondNode,
+                                      const GeneralBeamSection& section,
+                                      const LinearElasticMaterial& material);
+    Matrix localStiffness() const;
+    Matrix globalStiffness() const;
+    Vector localEquivalentLoad(const ConstantLocalLineLoad& load) const;
+    BeamRecovery recover(const Vector& globalElementDisplacement) const;
+
+private:
+    EulerBeam3D(double length,
+                double youngsModulus,
+                double shearModulus,
+                double area,
+                double iy,
+                double iz,
+                double torsionalConstant,
+                std::array<double, 9> rotation,
+                std::vector<std::array<double, 2>> sectionPoints);
+
+    double length_;
+    double youngsModulus_;
+    double shearModulus_;
+    double area_;
+    double iy_;
+    double iz_;
+    double torsionalConstant_;
+    std::array<double, 9> rotation_;
+    std::vector<std::array<double, 2>> sectionPoints_;
+};
+
+} // namespace fesa

src/fesa/CMakeLists.txt

@@ -6,6 +6,7 @@ add_library(
     build_info.cpp
     core/diagnostic.cpp
     core/status.cpp
+    elements/euler_beam_3d.cpp
     fem/dof_manager.cpp
     io/abaqus/domain_mapper.cpp
     io/abaqus/input_reader.cpp

src/fesa/elements/euler_beam_3d.cpp

@@ -0,0 +1,458 @@
+#include "fesa/elements/euler_beam_3d.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t kElementDofCount = 12U;
+constexpr std::size_t kGeneralizedComponentCount = 4U;
+constexpr double kGeometryTolerance = 1.0e-12;
+constexpr double kStiffnessInvariantTolerance = 1.0e-12;
+
+using Vector3 = std::array<double, 3>;
+
+double norm(const Vector3& value) {
+    return std::hypot(value[0], value[1], value[2]);
+}
+
+double dot(const Vector3& lhs, const Vector3& rhs) {
+    return lhs[0] * rhs[0] + lhs[1] * rhs[1] + lhs[2] * rhs[2];
+}
+
+Vector3 cross(const Vector3& lhs, const Vector3& rhs) {
+    return {
+        lhs[1] * rhs[2] - lhs[2] * rhs[1],
+        lhs[2] * rhs[0] - lhs[0] * rhs[2],
+        lhs[0] * rhs[1] - lhs[1] * rhs[0]};
+}
+
+bool isFinite(const Vector3& value) {
+    return std::isfinite(value[0]) && std::isfinite(value[1]) &&
+           std::isfinite(value[2]);
+}
+
+std::string elementIdentity(const Node& firstNode, const Node& secondNode) {
+    return firstNode.sourceId.instanceName + ":" +
+           firstNode.sourceId.sourceLabelText + "-" +
+           secondNode.sourceId.sourceLabelText;
+}
+
+Result<EulerBeam3D> modelFailure(const std::string& code,
+                                 const SourceLocation& location,
+                                 const std::string& identity,
+                                 const std::string& message) {
+    return Result<EulerBeam3D>::failure(Status::failure(
+        FailureCategory::model,
+        {{Severity::error, code, location, "*ELEMENT", identity, message}}));
+}
+
+Matrix transformation(const std::array<double, 9>& rotation) {
+    Matrix result{kElementDofCount, kElementDofCount};
+    // Blocks preserve [translation, rotation] at node 1 then node 2.
+    for (std::size_t block = 0; block < 4U; ++block) {
+        for (std::size_t row = 0; row < 3U; ++row) {
+            for (std::size_t column = 0; column < 3U; ++column) {
+                result(block * 3U + row, block * 3U + column) =
+                    rotation[row * 3U + column];
+            }
+        }
+    }
+    return result;
+}
+
+Matrix strainDisplacement(double xi, double length) {
+    Matrix b{kGeneralizedComponentCount, kElementDofCount};
+    const double r = 0.5 * (1.0 + xi);
+    const double inverseLength = 1.0 / length;
+    const double inverseLengthSquared = inverseLength * inverseLength;
+
+    b(0U, 0U) = -inverseLength;
+    b(0U, 6U) = inverseLength;
+    b(1U, 3U) = -inverseLength;
+    b(1U, 9U) = inverseLength;
+
+    // theta_y=-w' makes kappa_y=-w''; theta_z=v' makes kappa_z=v''.
+    b(2U, 2U) = (6.0 - 12.0 * r) * inverseLengthSquared;
+    b(2U, 4U) = (-4.0 + 6.0 * r) * inverseLength;
+    b(2U, 8U) = (-6.0 + 12.0 * r) * inverseLengthSquared;
+    b(2U, 10U) = (-2.0 + 6.0 * r) * inverseLength;
+
+    b(3U, 1U) = (-6.0 + 12.0 * r) * inverseLengthSquared;
+    b(3U, 5U) = (-4.0 + 6.0 * r) * inverseLength;
+    b(3U, 7U) = (6.0 - 12.0 * r) * inverseLengthSquared;
+    b(3U, 11U) = (-2.0 + 6.0 * r) * inverseLength;
+    return b;
+}
+
+std::array<double, kGeneralizedComponentCount> constitutiveDiagonal(
+    double youngsModulus,
+    double shearModulus,
+    double area,
+    double iy,
+    double iz,
+    double torsionalConstant) {
+    return {
+        youngsModulus * area,
+        shearModulus * torsionalConstant,
+        youngsModulus * iy,
+        youngsModulus * iz};
+}
+
+Matrix closedStiffness(double length,
+                       const std::array<double, kGeneralizedComponentCount>& diagonal) {
+    Matrix closed{kElementDofCount, kElementDofCount};
+    const auto addBlock = [&closed](const std::array<std::size_t, 2>& indices,
+                                    double coefficient) {
+        closed(indices[0], indices[0]) = coefficient;
+        closed(indices[0], indices[1]) = -coefficient;
+        closed(indices[1], indices[0]) = -coefficient;
+        closed(indices[1], indices[1]) = coefficient;
+    };
+    addBlock({0U, 6U}, diagonal[0U] / length);
+    addBlock({3U, 9U}, diagonal[1U] / length);
+
+    const auto addBendingBlock = [&closed, length](
+                                     const std::array<std::size_t, 4>& indices,
+                                     double flexuralRigidity,
+                                     double rotationSign) {
+        const double value = 12.0 * flexuralRigidity /
+                             (length * length * length);
+        const double coupling = rotationSign * 6.0 * flexuralRigidity /
+                                (length * length);
+        const double diagonalRotation = 4.0 * flexuralRigidity / length;
+        const double offDiagonalRotation = 2.0 * flexuralRigidity / length;
+        const std::array<double, 16> block = {
+            value, coupling, -value, coupling,
+            coupling, diagonalRotation, -coupling, offDiagonalRotation,
+            -value, -coupling, value, -coupling,
+            coupling, offDiagonalRotation, -coupling, diagonalRotation};
+        for (std::size_t row = 0; row < indices.size(); ++row) {
+            for (std::size_t column = 0; column < indices.size(); ++column) {
+                closed(indices[row], indices[column]) = block[row * indices.size() + column];
+            }
+        }
+    };
+    addBendingBlock({1U, 5U, 7U, 11U}, diagonal[3U], 1.0);
+    addBendingBlock({2U, 4U, 8U, 10U}, diagonal[2U], -1.0);
+    return closed;
+}
+
+double normalizedMatrixError(const Matrix& lhs, const Matrix& rhs) {
+    double maximumDifference = 0.0;
+    double scale = 1.0;
+    for (std::size_t row = 0; row < lhs.rows(); ++row) {
+        for (std::size_t column = 0; column < lhs.columns(); ++column) {
+            maximumDifference = (std::max)(
+                maximumDifference,
+                std::abs(lhs(row, column) - rhs(row, column)));
+            scale = (std::max)(scale, std::abs(lhs(row, column)));
+            scale = (std::max)(scale, std::abs(rhs(row, column)));
+        }
+    }
+    return maximumDifference / scale;
+}
+
+std::array<double, kGeneralizedComponentCount> generalizedStrain(
+    const Matrix& b,
+    const Vector& localDisplacement) {
+    std::array<double, kGeneralizedComponentCount> strain{};
+    for (std::size_t component = 0; component < strain.size(); ++component) {
+        for (std::size_t dof = 0; dof < localDisplacement.size(); ++dof) {
+            strain[component] += b(component, dof) * localDisplacement[dof];
+        }
+    }
+    return strain;
+}
+
+std::array<double, kGeneralizedComponentCount> generalizedResultant(
+    const std::array<double, kGeneralizedComponentCount>& strain,
+    const std::array<double, kGeneralizedComponentCount>& diagonal) {
+    std::array<double, kGeneralizedComponentCount> resultant{};
+    for (std::size_t component = 0; component < resultant.size(); ++component) {
+        resultant[component] = diagonal[component] * strain[component];
+    }
+    return resultant;
+}
+
+Matrix kinematicInterpolation(double xi, double length) {
+    Matrix interpolation{4U, kElementDofCount};
+    const double r = 0.5 * (1.0 + xi);
+    const double rSquared = r * r;
+    const double rCubed = rSquared * r;
+    const double n1 = 1.0 - r;
+    const double n2 = r;
+    const double h1 = 1.0 - 3.0 * rSquared + 2.0 * rCubed;
+    const double h2 = length * (r - 2.0 * rSquared + rCubed);
+    const double h3 = 3.0 * rSquared - 2.0 * rCubed;
+    const double h4 = length * (-rSquared + rCubed);
+
+    interpolation(0U, 0U) = n1;
+    interpolation(0U, 6U) = n2;
+    interpolation(1U, 1U) = h1;
+    interpolation(1U, 5U) = h2;
+    interpolation(1U, 7U) = h3;
+    interpolation(1U, 11U) = h4;
+    interpolation(2U, 2U) = h1;
+    interpolation(2U, 4U) = -h2;
+    interpolation(2U, 8U) = h3;
+    interpolation(2U, 10U) = -h4;
+    interpolation(3U, 3U) = n1;
+    interpolation(3U, 9U) = n2;
+    return interpolation;
+}
+
+} // namespace
+
+Result<EulerBeam3D> EulerBeam3D::create(
+    const Node& firstNode,
+    const Node& secondNode,
+    const GeneralBeamSection& section,
+    const LinearElasticMaterial& material) {
+    const std::string identity = elementIdentity(firstNode, secondNode);
+    const Vector3& first = firstNode.coordinates;
+    const Vector3& second = secondNode.coordinates;
+    const Vector3 delta = {
+        second[0] - first[0], second[1] - first[1], second[2] - first[2]};
+    const double length = norm(delta);
+    const double coordinateScale =
+        (std::max)({1.0, norm(first), norm(second)});
+    if (!isFinite(first) || !isFinite(second) || !isFinite(delta) ||
+        !std::isfinite(length) || !std::isfinite(coordinateScale) ||
+        !(length > kGeometryTolerance * coordinateScale)) {
+        return modelFailure(
+            "invalid-beam-length",
+            firstNode.location,
+            identity,
+            "Beam length must exceed the scale-aware geometry threshold.");
+    }
+
+    const Vector3 ex = {delta[0] / length, delta[1] / length, delta[2] / length};
+    const Vector3& guide = section.firstAxis;
+    const double guideNorm = norm(guide);
+    const double guideProjection = dot(guide, ex);
+    const Vector3 eyTrial = {
+        guide[0] - guideProjection * ex[0],
+        guide[1] - guideProjection * ex[1],
+        guide[2] - guideProjection * ex[2]};
+    const double eyTrialNorm = norm(eyTrial);
+    if (!isFinite(guide) || !std::isfinite(guideNorm) || !isFinite(eyTrial) ||
+        !std::isfinite(eyTrialNorm) ||
+        !(eyTrialNorm > kGeometryTolerance * (std::max)(1.0, guideNorm))) {
+        return modelFailure(
+            "invalid-beam-guide-vector",
+            section.location,
+            identity,
+            "Beam guide vector must define a scale-aware transverse direction.");
+    }
+
+    if (!std::isfinite(section.i12)) {
+        return modelFailure(
+            "invalid-beam-property",
+            section.location,
+            identity,
+            "Beam section properties must be finite and positive.");
+    }
+    if (section.i12 != 0.0) {
+        return modelFailure(
+            "unsupported-coupled-section",
+            section.location,
+            identity,
+            "The Euler beam kernel requires exact I12=0.");
+    }
+
+    const double shearModulus =
+        material.youngsModulus / (2.0 * (1.0 + material.poissonRatio));
+    const std::array<double, 6> positiveProperties = {
+        material.youngsModulus,
+        shearModulus,
+        section.area,
+        section.i11,
+        section.i22,
+        section.torsionalConstant};
+    if (!std::isfinite(material.poissonRatio) ||
+        std::any_of(
+            positiveProperties.begin(),
+            positiveProperties.end(),
+            [](double property) { return !std::isfinite(property) || !(property > 0.0); })) {
+        return modelFailure(
+            "invalid-beam-property",
+            section.location,
+            identity,
+            "E, G, A, Iy, Iz, and J must be finite and positive.");
+    }
+
+    const Vector3 ey = {
+        eyTrial[0] / eyTrialNorm,
+        eyTrial[1] / eyTrialNorm,
+        eyTrial[2] / eyTrialNorm};
+    const Vector3 ez = cross(ex, ey);
+    // Rows map global vectors to the approved right-handed local (ex,ey,ez) basis.
+    const std::array<double, 9> rotation = {
+        ex[0], ex[1], ex[2],
+        ey[0], ey[1], ey[2],
+        ez[0], ez[1], ez[2]};
+
+    return Result<EulerBeam3D>::success(EulerBeam3D{
+        length,
+        material.youngsModulus,
+        shearModulus,
+        section.area,
+        section.i11,
+        section.i22,
+        section.torsionalConstant,
+        rotation,
+        section.sectionPoints});
+}
+
+Matrix EulerBeam3D::localStiffness() const {
+    const auto diagonal = constitutiveDiagonal(
+        youngsModulus_, shearModulus_, area_, iy_, iz_, torsionalConstant_);
+    Matrix stiffness{kElementDofCount, kElementDofCount};
+    const double inverseSqrtThree = 1.0 / std::sqrt(3.0);
+    const std::array<double, 2> gaussPoints = {-inverseSqrtThree, inverseSqrtThree};
+    const double jacobian = 0.5 * length_;
+
+    // Both Gauss points are required: a one-point bending rule loses two ranks.
+    for (const double xi : gaussPoints) {
+        const Matrix b = strainDisplacement(xi, length_);
+        for (std::size_t row = 0; row < kElementDofCount; ++row) {
+            for (std::size_t column = 0; column < kElementDofCount; ++column) {
+                for (std::size_t component = 0; component < diagonal.size(); ++component) {
+                    stiffness(row, column) +=
+                        b(component, row) * diagonal[component] *
+                        b(component, column) * jacobian;
+                }
+            }
+        }
+    }
+
+    const Matrix closed = closedStiffness(length_, diagonal);
+    if (normalizedMatrixError(stiffness, closed) > kStiffnessInvariantTolerance) {
+        throw std::logic_error{
+            "Two-point Euler beam stiffness violated the closed-form invariant."};
+    }
+    return stiffness;
+}
+
+Matrix EulerBeam3D::globalStiffness() const {
+    const Matrix local = localStiffness();
+    const Matrix transform = transformation(rotation_);
+    const Matrix localTimesTransform = local.multiply(transform);
+    Matrix global{kElementDofCount, kElementDofCount};
+    // Kg=T^T*Kl*T while dl=T*dg.
+    for (std::size_t row = 0; row < kElementDofCount; ++row) {
+        for (std::size_t column = 0; column < kElementDofCount; ++column) {
+            for (std::size_t inner = 0; inner < kElementDofCount; ++inner) {
+                global(row, column) +=
+                    transform(inner, row) * localTimesTransform(inner, column);
+            }
+        }
+    }
+    return global;
+}
+
+Vector EulerBeam3D::localEquivalentLoad(const ConstantLocalLineLoad& load) const {
+    const std::array<double, 4> components = {load.px, load.py, load.pz, load.mx};
+    Vector equivalent{kElementDofCount};
+    const double inverseSqrtThree = 1.0 / std::sqrt(3.0);
+    const std::array<double, 2> gaussPoints = {-inverseSqrtThree, inverseSqrtThree};
+    const double jacobian = 0.5 * length_;
+    for (const double xi : gaussPoints) {
+        const Matrix interpolation = kinematicInterpolation(xi, length_);
+        for (std::size_t dof = 0; dof < equivalent.size(); ++dof) {
+            for (std::size_t component = 0; component < components.size(); ++component) {
+                equivalent[dof] +=
+                    interpolation(component, dof) * components[component] * jacobian;
+            }
+        }
+    }
+    return equivalent;
+}
+
+BeamRecovery EulerBeam3D::recover(const Vector& globalElementDisplacement) const {
+    const Matrix transform = transformation(rotation_);
+    const Vector localDisplacement = transform.multiply(globalElementDisplacement);
+    const auto diagonal = constitutiveDiagonal(
+        youngsModulus_, shearModulus_, area_, iy_, iz_, torsionalConstant_);
+    BeamRecovery recovery{};
+
+    // With parser/CLI distributed loading excluded, Kl*dl is the local outward end action.
+    const Vector endAction = localStiffness().multiply(localDisplacement);
+    for (std::size_t endpoint = 0; endpoint < 2U; ++endpoint) {
+        for (std::size_t component = 0; component < 6U; ++component) {
+            recovery.equilibriumEndActions[endpoint][component] =
+                endAction[endpoint * 6U + component];
+        }
+
+        const double xi = endpoint == 0U ? -1.0 : 1.0;
+        recovery.endpointSectionResultants[endpoint] = generalizedResultant(
+            generalizedStrain(strainDisplacement(xi, length_), localDisplacement),
+            diagonal);
+    }
+
+    const double inverseSqrtThree = 1.0 / std::sqrt(3.0);
+    const std::array<double, 2> gaussPoints = {-inverseSqrtThree, inverseSqrtThree};
+    for (std::size_t point = 0; point < gaussPoints.size(); ++point) {
+        recovery.gaussGeneralizedStrains[point] = generalizedStrain(
+            strainDisplacement(gaussPoints[point], length_), localDisplacement);
+        recovery.gaussGeneralizedResultants[point] = generalizedResultant(
+            recovery.gaussGeneralizedStrains[point], diagonal);
+
+        if (sectionPoints_.empty()) {
+            recovery.stressPoints.push_back({
+                static_cast<int>(point + 1U),
+                0U,
+                0.0,
+                0.0,
+                youngsModulus_ * recovery.gaussGeneralizedStrains[point][0U],
+                "fesa-default"});
+            continue;
+        }
+
+        for (std::size_t sectionPoint = 0; sectionPoint < sectionPoints_.size();
+             ++sectionPoint) {
+            const double x1 = sectionPoints_[sectionPoint][0U];
+            const double x2 = sectionPoints_[sectionPoint][1U];
+            const auto& strain = recovery.gaussGeneralizedStrains[point];
+            // x1=y and x2=z: S11=E(epsilon0+x2*kappa_y-x1*kappa_z).
+            recovery.stressPoints.push_back({
+                static_cast<int>(point + 1U),
+                sectionPoint + 1U,
+                x1,
+                x2,
+                youngsModulus_ *
+                    (strain[0U] + x2 * strain[2U] - x1 * strain[3U]),
+                "input"});
+        }
+    }
+    return recovery;
+}
+
+EulerBeam3D::EulerBeam3D(
+    double length,
+    double youngsModulus,
+    double shearModulus,
+    double area,
+    double iy,
+    double iz,
+    double torsionalConstant,
+    std::array<double, 9> rotation,
+    std::vector<std::array<double, 2>> sectionPoints)
+    : length_{length},
+      youngsModulus_{youngsModulus},
+      shearModulus_{shearModulus},
+      area_{area},
+      iy_{iy},
+      iz_{iz},
+      torsionalConstant_{torsionalConstant},
+      rotation_{rotation},
+      sectionPoints_{std::move(sectionPoints)} {}
+
+} // namespace fesa

tests/CMakeLists.txt

@@ -8,6 +8,7 @@ add_executable(
     unit/core/diagnostic_test.cpp
     unit/core/source_identity_test.cpp
     unit/core/status_test.cpp
+    unit/elements/euler_beam_3d_test.cpp
     unit/fem/dof_manager_test.cpp
     unit/math/matrix_test.cpp
     unit/math/vector_test.cpp

tests/unit/elements/euler_beam_3d_test.cpp

@@ -0,0 +1,885 @@
+#include "fesa/elements/euler_beam_3d.hpp"
+
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t kElementDofCount = 12U;
+constexpr double kMatrixTolerance = 1.0e-12;
+constexpr double kRigidTolerance = 1.0e-10;
+constexpr double kAnalyticalTolerance = 1.0e-9;
+
+Node makeNode(std::array<double, 3> coordinates, std::size_t line) {
+    return {{"Beam-1", static_cast<std::int64_t>(line), std::to_string(line)},
+            coordinates,
+            {"beam-test.inp", line}};
+}
+
+LinearElasticMaterial makeMaterial(double youngsModulus = 210.0e9,
+                                   double poissonRatio = 0.3) {
+    return {"Steel", youngsModulus, poissonRatio, {"beam-test.inp", 20U}};
+}
+
+GeneralBeamSection makeSection(
+    std::array<double, 3> firstAxis = {0.0, 1.0, 0.0},
+    std::vector<std::array<double, 2>> sectionPoints = {}) {
+    return {"Section-1",
+            0.012,
+            2.5e-5,
+            0.0,
+            4.0e-5,
+            1.5e-5,
+            firstAxis,
+            std::move(sectionPoints),
+            {"beam-test.inp", 30U}};
+}
+
+EulerBeam3D requireBeam(const Node& firstNode,
+                        const Node& secondNode,
+                        const GeneralBeamSection& section,
+                        const LinearElasticMaterial& material) {
+    auto result = EulerBeam3D::create(firstNode, secondNode, section, material);
+    if (!result.hasValue()) {
+        throw std::runtime_error{"Expected a valid EulerBeam3D fixture."};
+    }
+    return std::move(result.value());
+}
+
+EulerBeam3D alignedBeam(double length,
+                        const GeneralBeamSection& section = makeSection(),
+                        const LinearElasticMaterial& material = makeMaterial()) {
+    return requireBeam(
+        makeNode({0.0, 0.0, 0.0}, 1U),
+        makeNode({length, 0.0, 0.0}, 2U),
+        section,
+        material);
+}
+
+double maximumAbsoluteEntry(const Matrix& matrix) {
+    double maximum = 0.0;
+    for (std::size_t row = 0; row < matrix.rows(); ++row) {
+        for (std::size_t column = 0; column < matrix.columns(); ++column) {
+            maximum = (std::max)(maximum, std::abs(matrix(row, column)));
+        }
+    }
+    return maximum;
+}
+
+double normalizedMatrixError(const Matrix& actual, const Matrix& expected) {
+    if (actual.rows() != expected.rows() || actual.columns() != expected.columns()) {
+        throw std::invalid_argument{"Matrix comparison requires equal shapes."};
+    }
+
+    double maximumDifference = 0.0;
+    for (std::size_t row = 0; row < actual.rows(); ++row) {
+        for (std::size_t column = 0; column < actual.columns(); ++column) {
+            maximumDifference = (std::max)(
+                maximumDifference,
+                std::abs(actual(row, column) - expected(row, column)));
+        }
+    }
+
+    const double scale = (std::max)(
+        1.0,
+        (std::max)(maximumAbsoluteEntry(actual), maximumAbsoluteEntry(expected)));
+    return maximumDifference / scale;
+}
+
+double vectorNorm(const Vector& vector) {
+    double sum = 0.0;
+    for (std::size_t index = 0; index < vector.size(); ++index) {
+        sum += vector[index] * vector[index];
+    }
+    return std::sqrt(sum);
+}
+
+double quadraticEnergy(const Matrix& matrix, const Vector& vector) {
+    const Vector product = matrix.multiply(vector);
+    double value = 0.0;
+    for (std::size_t index = 0; index < vector.size(); ++index) {
+        value += vector[index] * product[index];
+    }
+    return value;
+}
+
+void expectScaledNear(double actual, double expected, double relativeTolerance) {
+    const double scale = (std::max)(1.0, std::abs(expected));
+    EXPECT_LE(std::abs(actual - expected), relativeTolerance * scale);
+}
+
+Matrix expectedClosedStiffness(double length,
+                               const GeneralBeamSection& section,
+                               const LinearElasticMaterial& material) {
+    Matrix expected{kElementDofCount, kElementDofCount};
+    const double shearModulus =
+        material.youngsModulus / (2.0 * (1.0 + material.poissonRatio));
+
+    const auto addBlock = [&expected](const std::vector<std::size_t>& indices,
+                                      const std::vector<double>& values) {
+        const std::size_t width = indices.size();
+        for (std::size_t row = 0; row < width; ++row) {
+            for (std::size_t column = 0; column < width; ++column) {
+                expected(indices[row], indices[column]) = values[row * width + column];
+            }
+        }
+    };
+
+    const double axial = material.youngsModulus * section.area / length;
+    addBlock({0U, 6U}, {axial, -axial, -axial, axial});
+
+    const double torsion = shearModulus * section.torsionalConstant / length;
+    addBlock({3U, 9U}, {torsion, -torsion, -torsion, torsion});
+
+    const auto bendingBlock = [length](double flexuralRigidity, double rotationSign) {
+        const double v = 12.0 * flexuralRigidity / (length * length * length);
+        const double c = rotationSign * 6.0 * flexuralRigidity / (length * length);
+        const double d = 4.0 * flexuralRigidity / length;
+        const double e = 2.0 * flexuralRigidity / length;
+        return std::vector<double>{
+            v, c, -v, c,
+            c, d, -c, e,
+            -v, -c, v, -c,
+            c, e, -c, d};
+    };
+
+    addBlock(
+        {1U, 5U, 7U, 11U},
+        bendingBlock(material.youngsModulus * section.i22, 1.0));
+    addBlock(
+        {2U, 4U, 8U, 10U},
+        bendingBlock(material.youngsModulus * section.i11, -1.0));
+    return expected;
+}
+
+std::array<double, kElementDofCount> symmetricEigenvalues(Matrix matrix) {
+    for (std::size_t iteration = 0; iteration < 100U * kElementDofCount; ++iteration) {
+        std::size_t p = 0U;
+        std::size_t q = 1U;
+        double maximumOffDiagonal = 0.0;
+        for (std::size_t row = 0; row < kElementDofCount; ++row) {
+            for (std::size_t column = row + 1U; column < kElementDofCount; ++column) {
+                const double candidate = std::abs(matrix(row, column));
+                if (candidate > maximumOffDiagonal) {
+                    maximumOffDiagonal = candidate;
+                    p = row;
+                    q = column;
+                }
+            }
+        }
+        if (maximumOffDiagonal <=
+            1.0e-14 * (std::max)(1.0, maximumAbsoluteEntry(matrix))) {
+            break;
+        }
+
+        const double app = matrix(p, p);
+        const double aqq = matrix(q, q);
+        const double apq = matrix(p, q);
+        const double angle = 0.5 * std::atan2(2.0 * apq, aqq - app);
+        const double cosine = std::cos(angle);
+        const double sine = std::sin(angle);
+
+        for (std::size_t index = 0; index < kElementDofCount; ++index) {
+            if (index == p || index == q) {
+                continue;
+            }
+            const double aip = matrix(index, p);
+            const double aiq = matrix(index, q);
+            matrix(index, p) = cosine * aip - sine * aiq;
+            matrix(p, index) = matrix(index, p);
+            matrix(index, q) = sine * aip + cosine * aiq;
+            matrix(q, index) = matrix(index, q);
+        }
+
+        matrix(p, p) = cosine * cosine * app - 2.0 * sine * cosine * apq +
+                       sine * sine * aqq;
+        matrix(q, q) = sine * sine * app + 2.0 * sine * cosine * apq +
+                       cosine * cosine * aqq;
+        matrix(p, q) = 0.0;
+        matrix(q, p) = 0.0;
+    }
+
+    std::array<double, kElementDofCount> eigenvalues{};
+    for (std::size_t index = 0; index < kElementDofCount; ++index) {
+        eigenvalues[index] = matrix(index, index);
+    }
+    return eigenvalues;
+}
+
+std::size_t symmetricRank(const Matrix& matrix, double relativeTolerance) {
+    const auto eigenvalues = symmetricEigenvalues(matrix);
+    double maximum = 0.0;
+    for (const double value : eigenvalues) {
+        maximum = (std::max)(maximum, std::abs(value));
+    }
+    return static_cast<std::size_t>(std::count_if(
+        eigenvalues.begin(),
+        eigenvalues.end(),
+        [maximum, relativeTolerance](double value) {
+            return std::abs(value) > relativeTolerance * maximum;
+        }));
+}
+
+Matrix testOnlyOnePointStiffness(double length,
+                                 const GeneralBeamSection& section,
+                                 const LinearElasticMaterial& material) {
+    // At xi=0 the two bending curvature rows retain only the nodal rotations.
+    // This deliberately under-integrated negative control is independent of production.
+    Matrix b{4U, kElementDofCount};
+    b(0U, 0U) = -1.0 / length;
+    b(0U, 6U) = 1.0 / length;
+    b(1U, 3U) = -1.0 / length;
+    b(1U, 9U) = 1.0 / length;
+    b(2U, 4U) = -1.0 / length;
+    b(2U, 10U) = 1.0 / length;
+    b(3U, 5U) = -1.0 / length;
+    b(3U, 11U) = 1.0 / length;
+
+    const double shearModulus =
+        material.youngsModulus / (2.0 * (1.0 + material.poissonRatio));
+    const std::array<double, 4> diagonal = {
+        material.youngsModulus * section.area,
+        shearModulus * section.torsionalConstant,
+        material.youngsModulus * section.i11,
+        material.youngsModulus * section.i22};
+
+    Matrix stiffness{kElementDofCount, kElementDofCount};
+    for (std::size_t row = 0; row < kElementDofCount; ++row) {
+        for (std::size_t column = 0; column < kElementDofCount; ++column) {
+            for (std::size_t component = 0; component < diagonal.size(); ++component) {
+                stiffness(row, column) +=
+                    b(component, row) * diagonal[component] * b(component, column) * length;
+            }
+        }
+    }
+    return stiffness;
+}
+
+Vector solveFixedFirstNode(const Matrix& stiffness,
+                           const std::array<double, 6>& freeEndLoad) {
+    std::array<std::array<double, 7>, 6> augmented{};
+    for (std::size_t row = 0; row < 6U; ++row) {
+        for (std::size_t column = 0; column < 6U; ++column) {
+            augmented[row][column] = stiffness(row + 6U, column + 6U);
+        }
+        augmented[row][6U] = freeEndLoad[row];
+    }
+
+    for (std::size_t pivot = 0; pivot < 6U; ++pivot) {
+        std::size_t pivotRow = pivot;
+        for (std::size_t row = pivot + 1U; row < 6U; ++row) {
+            if (std::abs(augmented[row][pivot]) >
+                std::abs(augmented[pivotRow][pivot])) {
+                pivotRow = row;
+            }
+        }
+        if (std::abs(augmented[pivotRow][pivot]) <=
+            std::numeric_limits<double>::min()) {
+            throw std::runtime_error{"Cantilever fixture is singular."};
+        }
+        std::swap(augmented[pivot], augmented[pivotRow]);
+
+        const double pivotValue = augmented[pivot][pivot];
+        for (std::size_t column = pivot; column < 7U; ++column) {
+            augmented[pivot][column] /= pivotValue;
+        }
+        for (std::size_t row = 0; row < 6U; ++row) {
+            if (row == pivot) {
+                continue;
+            }
+            const double factor = augmented[row][pivot];
+            for (std::size_t column = pivot; column < 7U; ++column) {
+                augmented[row][column] -= factor * augmented[pivot][column];
+            }
+        }
+    }
+
+    Vector displacement{kElementDofCount};
+    for (std::size_t component = 0; component < 6U; ++component) {
+        displacement[component + 6U] = augmented[component][6U];
+    }
+    return displacement;
+}
+
+Matrix transformationFromKnownRows(
+    const std::array<std::array<double, 3>, 3>& rotation) {
+    Matrix transformation{kElementDofCount, kElementDofCount};
+    for (std::size_t block = 0; block < 4U; ++block) {
+        for (std::size_t row = 0; row < 3U; ++row) {
+            for (std::size_t column = 0; column < 3U; ++column) {
+                transformation(block * 3U + row, block * 3U + column) =
+                    rotation[row][column];
+            }
+        }
+    }
+    return transformation;
+}
+
+Vector transposeMultiply(const Matrix& matrix, const Vector& vector) {
+    if (matrix.rows() != vector.size()) {
+        throw std::invalid_argument{"Transpose multiply dimension mismatch."};
+    }
+    Vector result{matrix.columns()};
+    for (std::size_t column = 0; column < matrix.columns(); ++column) {
+        for (std::size_t row = 0; row < matrix.rows(); ++row) {
+            result[column] += matrix(row, column) * vector[row];
+        }
+    }
+    return result;
+}
+
+double determinant(const std::array<std::array<double, 3>, 3>& matrix) {
+    return matrix[0][0] *
+               (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]) -
+           matrix[0][1] *
+               (matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0]) +
+           matrix[0][2] *
+               (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]);
+}
+
+TEST(EulerBeam3D, HermiteAndBMatrixMatchReviewedSigns) {
+    const double length = 2.5;
+    const auto section = makeSection();
+    const auto material = makeMaterial();
+    const auto beam = alignedBeam(length, section, material);
+
+    const double axialStrain = 0.012;
+    const double twist = -0.021;
+    const std::array<double, 4> v = {1.2, -0.4, 0.3, -0.07};
+    const std::array<double, 4> w = {-0.8, 0.6, -0.2, 0.05};
+    const auto value = [](const std::array<double, 4>& coefficients, double x) {
+        return coefficients[0] + coefficients[1] * x + coefficients[2] * x * x +
+               coefficients[3] * x * x * x;
+    };
+    const auto slope = [](const std::array<double, 4>& coefficients, double x) {
+        return coefficients[1] + 2.0 * coefficients[2] * x +
+               3.0 * coefficients[3] * x * x;
+    };
+    const auto curvature = [](const std::array<double, 4>& coefficients, double x) {
+        return 2.0 * coefficients[2] + 6.0 * coefficients[3] * x;
+    };
+
+    Vector displacement{kElementDofCount};
+    displacement[0U] = 0.2;
+    displacement[1U] = value(v, 0.0);
+    displacement[2U] = value(w, 0.0);
+    displacement[3U] = -0.1;
+    displacement[4U] = -slope(w, 0.0);
+    displacement[5U] = slope(v, 0.0);
+    displacement[6U] = displacement[0U] + axialStrain * length;
+    displacement[7U] = value(v, length);
+    displacement[8U] = value(w, length);
+    displacement[9U] = displacement[3U] + twist * length;
+    displacement[10U] = -slope(w, length);
+    displacement[11U] = slope(v, length);
+
+    const BeamRecovery recovery = beam.recover(displacement);
+    const double inverseSqrtThree = 1.0 / std::sqrt(3.0);
+    const std::array<double, 2> gaussXi = {-inverseSqrtThree, inverseSqrtThree};
+    for (std::size_t point = 0; point < gaussXi.size(); ++point) {
+        const double x = 0.5 * length * (1.0 + gaussXi[point]);
+        EXPECT_NEAR(recovery.gaussGeneralizedStrains[point][0U], axialStrain, 1.0e-14);
+        EXPECT_NEAR(recovery.gaussGeneralizedStrains[point][1U], twist, 1.0e-14);
+        EXPECT_NEAR(
+            recovery.gaussGeneralizedStrains[point][2U],
+            -curvature(w, x),
+            1.0e-13);
+        EXPECT_NEAR(
+            recovery.gaussGeneralizedStrains[point][3U],
+            curvature(v, x),
+            1.0e-13);
+    }
+
+    EXPECT_NEAR(
+        recovery.endpointSectionResultants[0U][2U],
+        material.youngsModulus * section.i11 * -curvature(w, 0.0),
+        1.0e-5);
+    EXPECT_NEAR(
+        recovery.endpointSectionResultants[1U][2U],
+        material.youngsModulus * section.i11 * -curvature(w, length),
+        1.0e-5);
+    EXPECT_NEAR(
+        recovery.endpointSectionResultants[0U][3U],
+        material.youngsModulus * section.i22 * curvature(v, 0.0),
+        1.0e-5);
+    EXPECT_NEAR(
+        recovery.endpointSectionResultants[1U][3U],
+        material.youngsModulus * section.i22 * curvature(v, length),
+        1.0e-5);
+}
+
+TEST(EulerBeam3D, TwoPointGaussMatchesClosedStiffness) {
+    const double length = 3.7;
+    auto section = makeSection();
+    section.area = 0.019;
+    section.i11 = 3.1e-5;
+    section.i22 = 7.4e-5;
+    section.torsionalConstant = 2.2e-5;
+    const auto material = makeMaterial(73.0e9, 0.27);
+    const auto beam = alignedBeam(length, section, material);
+
+    const Matrix actual = beam.localStiffness();
+    const Matrix closed = expectedClosedStiffness(length, section, material);
+    EXPECT_LE(normalizedMatrixError(actual, closed), kMatrixTolerance);
+
+    Matrix transpose{actual.rows(), actual.columns()};
+    for (std::size_t row = 0; row < actual.rows(); ++row) {
+        for (std::size_t column = 0; column < actual.columns(); ++column) {
+            transpose(row, column) = actual(column, row);
+        }
+    }
+    EXPECT_LE(normalizedMatrixError(actual, transpose), kMatrixTolerance);
+}
+
+TEST(EulerBeam3D, HasSixRigidModesRankSixAndPositiveDeformationEnergy) {
+    const double length = 2.0;
+    auto section = makeSection();
+    section.area = 1.4;
+    section.i11 = 0.8;
+    section.i22 = 1.1;
+    section.torsionalConstant = 0.6;
+    const auto material = makeMaterial(5.0, 0.25);
+    const Matrix stiffness = alignedBeam(length, section, material).localStiffness();
+
+    std::array<Vector, 6> rigidModes = {
+        Vector{kElementDofCount}, Vector{kElementDofCount}, Vector{kElementDofCount},
+        Vector{kElementDofCount}, Vector{kElementDofCount}, Vector{kElementDofCount}};
+    rigidModes[0U][0U] = rigidModes[0U][6U] = 1.0;
+    rigidModes[1U][1U] = rigidModes[1U][7U] = 1.0;
+    rigidModes[2U][2U] = rigidModes[2U][8U] = 1.0;
+    rigidModes[3U][3U] = rigidModes[3U][9U] = 1.0;
+    rigidModes[4U][4U] = rigidModes[4U][10U] = 1.0;
+    rigidModes[4U][8U] = -length;
+    rigidModes[5U][5U] = rigidModes[5U][11U] = 1.0;
+    rigidModes[5U][7U] = length;
+
+    const double stiffnessScale = (std::max)(1.0, maximumAbsoluteEntry(stiffness));
+    for (const Vector& mode : rigidModes) {
+        const double normalizedResidual =
+            vectorNorm(stiffness.multiply(mode)) /
+            (stiffnessScale * (std::max)(1.0, vectorNorm(mode)));
+        EXPECT_LE(normalizedResidual, kRigidTolerance);
+    }
+
+    const std::array<double, kElementDofCount> q = {
+        1.0, 1.0, 1.0, length, length, length,
+        1.0, 1.0, 1.0, length, length, length};
+    Matrix scaled{kElementDofCount, kElementDofCount};
+    for (std::size_t row = 0; row < kElementDofCount; ++row) {
+        for (std::size_t column = 0; column < kElementDofCount; ++column) {
+            scaled(row, column) = stiffness(row, column) / (q[row] * q[column]);
+        }
+    }
+    const auto eigenvalues = symmetricEigenvalues(scaled);
+    double maximumSingularValue = 0.0;
+    for (const double value : eigenvalues) {
+        maximumSingularValue = (std::max)(maximumSingularValue, std::abs(value));
+    }
+    const auto positiveCount = std::count_if(
+        eigenvalues.begin(), eigenvalues.end(), [maximumSingularValue](double value) {
+            return std::abs(value) > kRigidTolerance * maximumSingularValue;
+        });
+    EXPECT_EQ(positiveCount, 6);
+    for (const double value : eigenvalues) {
+        EXPECT_GE(value, -kRigidTolerance * maximumSingularValue);
+    }
+
+    for (std::size_t component = 0; component < 6U; ++component) {
+        Vector deformation{kElementDofCount};
+        deformation[6U + component] = 1.0;
+        EXPECT_GT(quadraticEnergy(stiffness, deformation), 0.0);
+    }
+}
+
+TEST(EulerBeam3D, RotatedTransformPreservesWorkAndEnergy) {
+    const double inverseSqrtTwo = 1.0 / std::sqrt(2.0);
+    const std::array<std::array<double, 3>, 3> rotation = {{
+        {{2.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0}},
+        {{-inverseSqrtTwo, inverseSqrtTwo, 0.0}},
+        {{-inverseSqrtTwo / 3.0, -inverseSqrtTwo / 3.0,
+          4.0 * inverseSqrtTwo / 3.0}}}};
+    const Matrix transformation = transformationFromKnownRows(rotation);
+
+    for (std::size_t row = 0; row < 3U; ++row) {
+        for (std::size_t column = 0; column < 3U; ++column) {
+            double dot = 0.0;
+            for (std::size_t component = 0; component < 3U; ++component) {
+                dot += rotation[row][component] * rotation[column][component];
+            }
+            EXPECT_NEAR(dot, row == column ? 1.0 : 0.0, 1.0e-14);
+        }
+    }
+    EXPECT_NEAR(determinant(rotation), 1.0, 1.0e-14);
+
+    const auto section = makeSection({-2.0, 2.0, 0.0});
+    const auto material = makeMaterial();
+    const auto beam = requireBeam(
+        makeNode({1.0, -2.0, 0.5}, 1U),
+        makeNode({3.0, 0.0, 1.5}, 2U),
+        section,
+        material);
+    const Matrix local = beam.localStiffness();
+    const Matrix global = beam.globalStiffness();
+
+    Matrix expectedGlobal{kElementDofCount, kElementDofCount};
+    const Matrix localTimesTransform = local.multiply(transformation);
+    for (std::size_t row = 0; row < kElementDofCount; ++row) {
+        for (std::size_t column = 0; column < kElementDofCount; ++column) {
+            for (std::size_t inner = 0; inner < kElementDofCount; ++inner) {
+                expectedGlobal(row, column) +=
+                    transformation(inner, row) * localTimesTransform(inner, column);
+            }
+        }
+    }
+    EXPECT_LE(normalizedMatrixError(global, expectedGlobal), kMatrixTolerance);
+
+    Vector localDisplacement{kElementDofCount};
+    for (std::size_t index = 0; index < localDisplacement.size(); ++index) {
+        localDisplacement[index] = 0.01 * static_cast<double>(index + 1U) - 0.04;
+    }
+    const Vector globalDisplacement = transposeMultiply(transformation, localDisplacement);
+    const Vector localForce = local.multiply(localDisplacement);
+    const Vector globalForce = global.multiply(globalDisplacement);
+    const Vector expectedGlobalForce = transposeMultiply(transformation, localForce);
+    for (std::size_t index = 0; index < kElementDofCount; ++index) {
+        expectScaledNear(globalForce[index], expectedGlobalForce[index], kMatrixTolerance);
+    }
+    expectScaledNear(
+        quadraticEnergy(global, globalDisplacement),
+        quadraticEnergy(local, localDisplacement),
+        kMatrixTolerance);
+
+    Vector globalVariation{kElementDofCount};
+    for (std::size_t index = 0; index < globalVariation.size(); ++index) {
+        globalVariation[index] = 0.03 - 0.002 * static_cast<double>(index);
+    }
+    const Vector localVariation = transformation.multiply(globalVariation);
+    expectScaledNear(
+        globalVariation.dot(globalForce),
+        localVariation.dot(localForce),
+        kMatrixTolerance);
+
+    const BeamRecovery recovery = beam.recover(globalDisplacement);
+    EXPECT_NEAR(
+        recovery.gaussGeneralizedStrains[0U][0U],
+        (localDisplacement[6U] - localDisplacement[0U]) / 3.0,
+        1.0e-14);
+}
+
+TEST(EulerBeam3D, ConstantLineLoadMatchesAllSignedComponents) {
+    const double length = 4.0;
+    const ConstantLocalLineLoad load{2.5, -3.0, 5.5, -7.0};
+    const Vector equivalent = alignedBeam(length).localEquivalentLoad(load);
+    const std::array<double, kElementDofCount> expected = {
+        5.0, -6.0, 11.0, -14.0, -22.0 / 3.0, -4.0,
+        5.0, -6.0, 11.0, -14.0, 22.0 / 3.0, 4.0};
+    ASSERT_EQ(equivalent.size(), expected.size());
+    for (std::size_t index = 0; index < expected.size(); ++index) {
+        expectScaledNear(equivalent[index], expected[index], kMatrixTolerance);
+    }
+}
+
+TEST(EulerBeam3D, AnalyticalAxialTorsionAndTwoPlaneBendingRecover) {
+    const double length = 3.0;
+    const auto section = makeSection();
+    const auto material = makeMaterial();
+    const auto beam = alignedBeam(length, section, material);
+    const Matrix stiffness = beam.localStiffness();
+    const double shearModulus =
+        material.youngsModulus / (2.0 * (1.0 + material.poissonRatio));
+
+    const double axialForce = 1250.0;
+    const Vector axial = solveFixedFirstNode(stiffness, {axialForce, 0.0, 0.0, 0.0, 0.0, 0.0});
+    expectScaledNear(
+        axial[6U],
+        axialForce * length / (material.youngsModulus * section.area),
+        kAnalyticalTolerance);
+    const BeamRecovery axialRecovery = beam.recover(axial);
+    expectScaledNear(axialRecovery.equilibriumEndActions[0U][0U], -axialForce, kMatrixTolerance);
+    expectScaledNear(axialRecovery.equilibriumEndActions[1U][0U], axialForce, kMatrixTolerance);
+    expectScaledNear(axialRecovery.endpointSectionResultants[0U][0U], axialForce, kMatrixTolerance);
+    expectScaledNear(axialRecovery.endpointSectionResultants[1U][0U], axialForce, kMatrixTolerance);
+
+    const double torque = -870.0;
+    const Vector torsion = solveFixedFirstNode(stiffness, {0.0, 0.0, 0.0, torque, 0.0, 0.0});
+    expectScaledNear(
+        torsion[9U],
+        torque * length / (shearModulus * section.torsionalConstant),
+        kAnalyticalTolerance);
+    const BeamRecovery torsionRecovery = beam.recover(torsion);
+    expectScaledNear(torsionRecovery.equilibriumEndActions[0U][3U], -torque, kMatrixTolerance);
+    expectScaledNear(torsionRecovery.equilibriumEndActions[1U][3U], torque, kMatrixTolerance);
+    expectScaledNear(torsionRecovery.endpointSectionResultants[0U][1U], torque, kMatrixTolerance);
+
+    const double localYForce = 640.0;
+    const Vector localY = solveFixedFirstNode(stiffness, {0.0, localYForce, 0.0, 0.0, 0.0, 0.0});
+    expectScaledNear(
+        localY[7U],
+        localYForce * length * length * length /
+            (3.0 * material.youngsModulus * section.i22),
+        kAnalyticalTolerance);
+    expectScaledNear(
+        localY[11U],
+        localYForce * length * length /
+            (2.0 * material.youngsModulus * section.i22),
+        kAnalyticalTolerance);
+    const BeamRecovery localYRecovery = beam.recover(localY);
+    expectScaledNear(localYRecovery.equilibriumEndActions[0U][1U], -localYForce, kMatrixTolerance);
+    expectScaledNear(localYRecovery.equilibriumEndActions[1U][1U], localYForce, kMatrixTolerance);
+    expectScaledNear(
+        localYRecovery.equilibriumEndActions[0U][5U],
+        -localYForce * length,
+        kMatrixTolerance);
+    expectScaledNear(
+        localYRecovery.endpointSectionResultants[0U][3U],
+        localYForce * length,
+        kMatrixTolerance);
+    EXPECT_NEAR(localYRecovery.endpointSectionResultants[1U][3U], 0.0, 1.0e-8);
+
+    const double localZForce = -510.0;
+    const Vector localZ = solveFixedFirstNode(stiffness, {0.0, 0.0, localZForce, 0.0, 0.0, 0.0});
+    expectScaledNear(
+        localZ[8U],
+        localZForce * length * length * length /
+            (3.0 * material.youngsModulus * section.i11),
+        kAnalyticalTolerance);
+    expectScaledNear(
+        localZ[10U],
+        -localZForce * length * length /
+            (2.0 * material.youngsModulus * section.i11),
+        kAnalyticalTolerance);
+    const BeamRecovery localZRecovery = beam.recover(localZ);
+    expectScaledNear(localZRecovery.equilibriumEndActions[0U][2U], -localZForce, kMatrixTolerance);
+    expectScaledNear(localZRecovery.equilibriumEndActions[1U][2U], localZForce, kMatrixTolerance);
+    expectScaledNear(
+        localZRecovery.equilibriumEndActions[0U][4U],
+        localZForce * length,
+        kMatrixTolerance);
+    expectScaledNear(
+        localZRecovery.endpointSectionResultants[0U][2U],
+        -localZForce * length,
+        kMatrixTolerance);
+    EXPECT_NEAR(localZRecovery.endpointSectionResultants[1U][2U], 0.0, 1.0e-8);
+}
+
+TEST(EulerBeam3D, RejectsInvalidGeometryAndProperties) {
+    const Node origin = makeNode({0.0, 0.0, 0.0}, 1U);
+    const Node unitX = makeNode({1.0, 0.0, 0.0}, 2U);
+    const auto validSection = makeSection();
+    const auto validMaterial = makeMaterial();
+
+    const auto expectFailure = [](const Result<EulerBeam3D>& result,
+                                  const std::string& code) {
+        ASSERT_FALSE(result.hasValue());
+        EXPECT_EQ(result.status().failureCategory(), FailureCategory::model);
+        ASSERT_EQ(result.status().diagnostics().size(), 1U);
+        EXPECT_EQ(result.status().diagnostics()[0U].code, code);
+    };
+
+    expectFailure(
+        EulerBeam3D::create(origin, origin, validSection, validMaterial),
+        "invalid-beam-length");
+    expectFailure(
+        EulerBeam3D::create(
+            origin,
+            makeNode({1.0e-12, 0.0, 0.0}, 2U),
+            validSection,
+            validMaterial),
+        "invalid-beam-length");
+
+    const double coordinate = 1048576.0;
+    const Node scaledFirst = makeNode({coordinate, 0.0, 0.0}, 1U);
+    const Node belowThreshold = makeNode({coordinate + coordinate * 0.5e-12, 0.0, 0.0}, 2U);
+    const Node aboveThreshold = makeNode({coordinate + coordinate * 2.0e-12, 0.0, 0.0}, 2U);
+    expectFailure(
+        EulerBeam3D::create(scaledFirst, belowThreshold, validSection, validMaterial),
+        "invalid-beam-length");
+    EXPECT_TRUE(EulerBeam3D::create(
+                    scaledFirst, aboveThreshold, validSection, validMaterial)
+                    .hasValue());
+
+    auto parallelGuide = validSection;
+    parallelGuide.firstAxis = {1.0, 0.0, 0.0};
+    expectFailure(
+        EulerBeam3D::create(origin, unitX, parallelGuide, validMaterial),
+        "invalid-beam-guide-vector");
+    auto guideAtThreshold = validSection;
+    guideAtThreshold.firstAxis = {1.0, 1.0e-12, 0.0};
+    expectFailure(
+        EulerBeam3D::create(origin, unitX, guideAtThreshold, validMaterial),
+        "invalid-beam-guide-vector");
+    auto guideAboveThreshold = validSection;
+    guideAboveThreshold.firstAxis = {1.0, 2.0e-12, 0.0};
+    EXPECT_TRUE(EulerBeam3D::create(origin, unitX, guideAboveThreshold, validMaterial).hasValue());
+
+    auto invalidMaterial = validMaterial;
+    invalidMaterial.youngsModulus = 0.0;
+    expectFailure(
+        EulerBeam3D::create(origin, unitX, validSection, invalidMaterial),
+        "invalid-beam-property");
+    invalidMaterial = validMaterial;
+    invalidMaterial.poissonRatio = -2.0;
+    expectFailure(
+        EulerBeam3D::create(origin, unitX, validSection, invalidMaterial),
+        "invalid-beam-property");
+
+    for (std::size_t property = 0; property < 4U; ++property) {
+        auto invalidSection = validSection;
+        double* properties[] = {
+            &invalidSection.area,
+            &invalidSection.i11,
+            &invalidSection.i22,
+            &invalidSection.torsionalConstant};
+        *properties[property] = 0.0;
+        expectFailure(
+            EulerBeam3D::create(origin, unitX, invalidSection, validMaterial),
+            "invalid-beam-property");
+    }
+
+    auto coupledSection = validSection;
+    coupledSection.i12 = 1.0e-9;
+    expectFailure(
+        EulerBeam3D::create(origin, unitX, coupledSection, validMaterial),
+        "unsupported-coupled-section");
+}
+
+TEST(EulerBeam3D, RecoversSectionPointAndDefaultCentroidS11) {
+    const double length = 2.0;
+    const double epsilon = 0.01;
+    const double kappaY = 0.02;
+    const double kappaZ = -0.03;
+    const auto material = makeMaterial();
+    auto section = makeSection({0.0, 1.0, 0.0}, {{0.25, -0.5}, {-0.4, 0.3}});
+    const auto beam = alignedBeam(length, section, material);
+
+    Vector displacement{kElementDofCount};
+    displacement[6U] = epsilon * length;
+    displacement[7U] = 0.5 * kappaZ * length * length;
+    displacement[8U] = -0.5 * kappaY * length * length;
+    displacement[10U] = kappaY * length;
+    displacement[11U] = kappaZ * length;
+
+    const BeamRecovery recovery = beam.recover(displacement);
+    ASSERT_EQ(recovery.stressPoints.size(), 4U);
+    for (std::size_t gaussPoint = 0; gaussPoint < 2U; ++gaussPoint) {
+        for (std::size_t point = 0; point < section.sectionPoints.size(); ++point) {
+            const BeamStressPoint& stress =
+                recovery.stressPoints[gaussPoint * section.sectionPoints.size() + point];
+            const double x1 = section.sectionPoints[point][0U];
+            const double x2 = section.sectionPoints[point][1U];
+            EXPECT_EQ(stress.gaussPoint, static_cast<int>(gaussPoint + 1U));
+            EXPECT_EQ(stress.sectionPoint, point + 1U);
+            EXPECT_DOUBLE_EQ(stress.x1, x1);
+            EXPECT_DOUBLE_EQ(stress.x2, x2);
+            expectScaledNear(
+                stress.s11,
+                material.youngsModulus * (epsilon + x2 * kappaY - x1 * kappaZ),
+                kMatrixTolerance);
+        }
+    }
+
+    const auto defaultBeam = alignedBeam(length, makeSection(), material);
+    const BeamRecovery defaultRecovery = defaultBeam.recover(displacement);
+    ASSERT_EQ(defaultRecovery.stressPoints.size(), 2U);
+    for (std::size_t gaussPoint = 0; gaussPoint < 2U; ++gaussPoint) {
+        const BeamStressPoint& stress = defaultRecovery.stressPoints[gaussPoint];
+        EXPECT_EQ(stress.gaussPoint, static_cast<int>(gaussPoint + 1U));
+        EXPECT_EQ(stress.sectionPoint, 0U);
+        EXPECT_DOUBLE_EQ(stress.x1, 0.0);
+        EXPECT_DOUBLE_EQ(stress.x2, 0.0);
+        EXPECT_EQ(stress.source, "fesa-default");
+        expectScaledNear(
+            stress.s11,
+            material.youngsModulus * epsilon,
+            kMatrixTolerance);
+    }
+}
+
+TEST(EulerBeam3D, ReproducesConstantStrainTwistAndCurvaturePatches) {
+    const double length = 2.8;
+    const double epsilon = -0.014;
+    const double twist = 0.023;
+    const double kappaY = -0.031;
+    const double kappaZ = 0.047;
+    const auto section = makeSection();
+    const auto material = makeMaterial();
+    const double shearModulus =
+        material.youngsModulus / (2.0 * (1.0 + material.poissonRatio));
+    const auto beam = alignedBeam(length, section, material);
+
+    Vector displacement{kElementDofCount};
+    displacement[6U] = epsilon * length;
+    displacement[7U] = 0.5 * kappaZ * length * length;
+    displacement[8U] = -0.5 * kappaY * length * length;
+    displacement[9U] = twist * length;
+    displacement[10U] = kappaY * length;
+    displacement[11U] = kappaZ * length;
+
+    const std::array<double, 4> expectedStrain = {epsilon, twist, kappaY, kappaZ};
+    const std::array<double, 4> expectedResultant = {
+        material.youngsModulus * section.area * epsilon,
+        shearModulus * section.torsionalConstant * twist,
+        material.youngsModulus * section.i11 * kappaY,
+        material.youngsModulus * section.i22 * kappaZ};
+    const BeamRecovery recovery = beam.recover(displacement);
+    for (std::size_t point = 0; point < 2U; ++point) {
+        for (std::size_t component = 0; component < 4U; ++component) {
+            expectScaledNear(
+                recovery.gaussGeneralizedStrains[point][component],
+                expectedStrain[component],
+                kMatrixTolerance);
+            expectScaledNear(
+                recovery.gaussGeneralizedResultants[point][component],
+                expectedResultant[component],
+                kMatrixTolerance);
+            expectScaledNear(
+                recovery.endpointSectionResultants[point][component],
+                expectedResultant[component],
+                kMatrixTolerance);
+        }
+    }
+
+    const std::array<std::size_t, 4> endActionComponents = {0U, 3U, 4U, 5U};
+    for (std::size_t component = 0; component < expectedResultant.size(); ++component) {
+        expectScaledNear(
+            recovery.equilibriumEndActions[0U][endActionComponents[component]],
+            -expectedResultant[component],
+            kMatrixTolerance);
+        expectScaledNear(
+            recovery.equilibriumEndActions[1U][endActionComponents[component]],
+            expectedResultant[component],
+            kMatrixTolerance);
+    }
+    EXPECT_NEAR(recovery.equilibriumEndActions[0U][1U], 0.0, 1.0e-8);
+    EXPECT_NEAR(recovery.equilibriumEndActions[0U][2U], 0.0, 1.0e-8);
+    EXPECT_NEAR(recovery.equilibriumEndActions[1U][1U], 0.0, 1.0e-8);
+    EXPECT_NEAR(recovery.equilibriumEndActions[1U][2U], 0.0, 1.0e-8);
+}
+
+TEST(EulerBeam3D, OnePointNegativeControlHasRankFour) {
+    const double length = 3.7;
+    auto section = makeSection();
+    section.area = 1.0;
+    section.i11 = 0.7;
+    section.i22 = 1.2;
+    section.torsionalConstant = 0.9;
+    const auto material = makeMaterial(4.0, 0.25);
+
+    const Matrix onePoint = testOnlyOnePointStiffness(length, section, material);
+    const Matrix production = alignedBeam(length, section, material).localStiffness();
+    EXPECT_EQ(symmetricRank(onePoint, kRigidTolerance), 4U);
+    EXPECT_EQ(symmetricRank(production, kRigidTolerance), 6U);
+}
+
+} // namespace
+} // namespace fesa