feat(fem-and-beam-kernel): step 3 — timoshenko-stiffness-kernel

CMakeLists.txt

@@ -27,6 +27,7 @@ include(cmake/FesaDependencies.cmake)
 
 add_library(fesa_core STATIC
     src/fesa/core/version.cpp
+    src/fesa/elements/beam/beam3d2.cpp
     src/fesa/fem/beam_frame.cpp
     src/fesa/fem/dof_manager.cpp
     src/fesa/fem/gauss_rule.cpp

docs/formulation/timoshenko-beam-3d.md

@@ -0,0 +1,237 @@
+# 2절점 3D Timoshenko Beam 정식화
+
+## 1. 범위와 가정
+
+FESA Phase 1 Beam은 소변형, 선형 탄성의 2절점 직선
+isoparametric Timoshenko 요소다. 단면은 도심과 전단중심이 일치하는 주축 단면이며
+\(I_{yz}=0\)이다. 단면 offset, warping, 기하·재료 비선형은 포함하지 않는다.
+FESA는 단위 변환을 하지 않으므로 모든 입력은 하나의 일관 단위계를 사용해야 한다.
+
+정식화의 유한요소 이산화와 수치적분은 Bathe,
+*Finite Element Procedures*, 2판과 Hughes,
+*The Finite Element Method: Linear Static and Dynamic Finite Element
+Analysis*를 따른다. 3차원 isoparametric Beam의 국부 좌표계는 Bathe와
+Bolourchi(1979)의 선형화된 부분을 사용한다. 선택적 감차적분은 Hughes, Taylor와
+Kanoknukulchai(1977)의 원칙을 Beam 전단항에 적용한다. Plate 운동학은 사용하지
+않는다.
+
+## 2. 국부 좌표계와 자유도
+
+절점 1에서 절점 2로 향하는 단위벡터를 \(\mathbf e_x\)로 둔다. 입력 orientation
+벡터에서 \(\mathbf e_x\) 성분을 Gram-Schmidt로 제거하고 정규화한 벡터를
+\(\mathbf e_y\)로 두며,
+
+\[
+\mathbf e_z=\mathbf e_x\times\mathbf e_y
+\]
+
+로 정의한다. 따라서 \((\mathbf e_x,\mathbf e_y,\mathbf e_z)\)는 오른손 직교
+기저다. 길이가 0인 요소, 길이가 0인 orientation, 요소축과 평행한 orientation은
+오류다.
+
+국부 요소 자유도 벡터의 순서는 다음과 같다.
+
+\[
+\mathbf q_l =
+\begin{bmatrix}
+u_1&v_1&w_1&\theta_{x1}&\theta_{y1}&\theta_{z1}&
+u_2&v_2&w_2&\theta_{x2}&\theta_{y2}&\theta_{z2}
+\end{bmatrix}^{T}.
+\]
+
+\(u,v,w\)는 각각 국부 \(x,y,z\) 병진이고
+\(\theta_x,\theta_y,\theta_z\)는 오른손 규칙의 회전 성분이다. 전역 자유도도 절점별
+\([u_x,u_y,u_z,r_x,r_y,r_z]\) 순서를 사용한다.
+
+기준축에서 떨어진 점 \((0,y,z)\)의 선형화된 변위장은
+
+\[
+U_x=u+z\theta_y-y\theta_z,\qquad
+U_y=v-z\theta_x,\qquad
+U_z=w+y\theta_x
+\]
+
+이다.
+
+## 3. 변형률과 부호
+
+일반화 변형률과 대응 단면력의 component 순서는
+
+\[
+\boldsymbol\varepsilon_s =
+\begin{bmatrix}
+\epsilon&\gamma_y&\gamma_z&\kappa_x&\kappa_y&\kappa_z
+\end{bmatrix}^{T},
+\qquad
+\mathbf s =
+\begin{bmatrix}
+N&V_y&V_z&T&M_y&M_z
+\end{bmatrix}^{T}
+\]
+
+로 고정한다. 위 변위장과 오른손 회전 부호로부터
+
+\[
+\begin{aligned}
+\epsilon &= \frac{du}{dx},&
+\gamma_y &= \frac{dv}{dx}-\theta_z,&
+\gamma_z &= \frac{dw}{dx}+\theta_y,\\
+\kappa_x &= \frac{d\theta_x}{dx},&
+\kappa_y &= \frac{d\theta_y}{dx},&
+\kappa_z &= \frac{d\theta_z}{dx}
+\end{aligned}
+\]
+
+를 사용한다. 양의 \(N,V_y,V_z,T,M_y,M_z\)는 위 일반화 변형률과 양의 내부
+일률 \(\delta\boldsymbol\varepsilon_s^T\mathbf s\)을 이루는 방향이다. 이 부호에서
+단면 축응력은
+
+\[
+\sigma_{xx}=E\left(\epsilon+z\kappa_y-y\kappa_z\right)
+\]
+
+이고 \(M_y=\int_A z\sigma_{xx}\,dA\),
+\(M_z=-\int_A y\sigma_{xx}\,dA\)이다.
+
+## 4. 재료와 단면 constitutive matrix
+
+등방성 선형 탄성의 전단계수는 입력 \(E,\nu\)로부터
+
+\[
+G=\frac{E}{2(1+\nu)}
+\]
+
+로 계산한다. 일반화 constitutive 관계는
+
+\[
+\mathbf s=\mathbf D\boldsymbol\varepsilon_s,\qquad
+\mathbf D=
+\operatorname{diag}
+\left(EA,\;GA_{sy},\;GA_{sz},\;GJ,\;EI_y,\;EI_z\right).
+\]
+
+\(A,I_y,I_z,J,A_{sy},A_{sz}\)는 양의 값이어야 한다. 명시적인 전단강성이 없는
+입력은 semantic mapper가 커널 호출 전에
+\(A_{sy}=A_{sz}=5A/6\), `SCF=0`으로 정규화한다. 커널은 Abaqus의 slenderness
+compensation을 추가하지 않는다.
+
+## 5. 자연좌표, 형상함수와 Jacobian
+
+자연좌표는 \(\xi\in[-1,1]\)이고 두 절점의 선형 형상함수는
+
+\[
+N_1(\xi)=\frac{1-\xi}{2},\qquad
+N_2(\xi)=\frac{1+\xi}{2}.
+\]
+
+모든 여섯 component를 같은 형상함수로 보간한다. 직선 요소 길이를 \(L\)이라 하면
+
+\[
+J=\frac{dx}{d\xi}=\frac{L}{2},\qquad
+\frac{dN_a}{dx}=\frac{1}{J}\frac{dN_a}{d\xi}
+\]
+
+이다.
+
+절점 \(a\)의 자유도 순서
+\([u_a,v_a,w_a,\theta_{xa},\theta_{ya},\theta_{za}]\)에 대한
+strain-displacement block은
+
+\[
+\mathbf B_a(\xi)=
+\begin{bmatrix}
+N_{a,x}&0&0&0&0&0\\
+0&N_{a,x}&0&0&0&-N_a\\
+0&0&N_{a,x}&0&N_a&0\\
+0&0&0&N_{a,x}&0&0\\
+0&0&0&0&N_{a,x}&0\\
+0&0&0&0&0&N_{a,x}
+\end{bmatrix},
+\qquad
+\mathbf B=\begin{bmatrix}\mathbf B_1&\mathbf B_2\end{bmatrix}.
+\]
+
+## 6. 선택적 감차적분과 요소 강성
+
+Constitutive matrix를 축·비틀림·굽힘 부분과 전단 부분으로 나눈다.
+
+\[
+\begin{aligned}
+\mathbf D_{ab}&=
+\operatorname{diag}(EA,0,0,GJ,EI_y,EI_z),\\
+\mathbf D_s&=
+\operatorname{diag}(0,GA_{sy},GA_{sz},0,0,0).
+\end{aligned}
+\]
+
+국부 강성은
+
+\[
+\mathbf K_l =
+\sum_{g=1}^{2}
+\mathbf B(\xi_g)^T\mathbf D_{ab}\mathbf B(\xi_g)Jw_g
++
+\sum_{g=1}^{1}
+\mathbf B(\xi_g)^T\mathbf D_s\mathbf B(\xi_g)Jw_g
+\]
+
+로 직접 적분한다. 축·비틀림·굽힘은
+\(\xi=\pm1/\sqrt{3},w=1\)인 2점 Gauss rule을 사용하고, 전단은
+\(\xi=0,w=2\)인 1점 Gauss rule만 사용한다. 전단항을 2점 적분하지 않는다.
+선형 회전장에서 1점 전단 적분은 요소 중앙의 평균 회전을 사용하므로 일정곡률
+굽힘 mode의 불필요한 전단에너지를 제거한다.
+
+`Matrix12`는 `matrix[row][column]`의 12행 12열 고정 크기 저장을 사용한다. 적분은
+각 Gauss point에서 \(\mathbf B^T\mathbf D\mathbf B\)를 대칭 양쪽 entry에 같은
+부동소수점 연산 순서로 누적한다.
+
+## 7. 전역 변환
+
+기존 `beam_transformation`이 만드는 \(\mathbf T\)는 각 병진·회전 3-vector를
+전역 성분에서 국부 성분으로 투영한다.
+
+\[
+\mathbf q_l=\mathbf T\mathbf q_g.
+\]
+
+따라서 전역 강성과 에너지 관계는
+
+\[
+\mathbf K_g=\mathbf T^T\mathbf K_l\mathbf T,\qquad
+\frac12\mathbf q_g^T\mathbf K_g\mathbf q_g
+=\frac12\mathbf q_l^T\mathbf K_l\mathbf q_l
+\]
+
+이다.
+
+## 8. 구현과 검증 invariant
+
+커널은 위 적분식으로 local stiffness를 계산하고 변환식으로 global stiffness를
+계산한다. 다음 invariant가 정식화 검증 기준이다.
+
+- local/global stiffness는 roundoff 이내에서 대칭이다.
+- 세 병진과 세 회전 강체 mode의 변형에너지는 0이다.
+- 축 및 비틀림 submatrix는 각각 \(EA/L\), \(GJ/L\)을 재현한다.
+- 일정곡률 단축·이축 굽힘 에너지는 각각
+  \(EI_y\kappa_y^2L/2\), \(EI_z\kappa_z^2L/2\)이다.
+- 일정전단 mode의 에너지는 각각
+  \(GA_{sy}\gamma_y^2L/2\), \(GA_{sz}\gamma_z^2L/2\)이다.
+- 일정곡률 mode의 전단에너지는 slenderness 변화와 무관하게 0이며 전체 에너지는
+  해석값과 일치한다.
+- 요소와 자유도를 같은 강체 회전으로 회전하면 변형에너지가 보존된다.
+
+## 9. 참고문헌과 FESA 차이
+
+- K. J. Bathe, *Finite Element Procedures*, 2nd ed.
+- T. J. R. Hughes, *The Finite Element Method: Linear Static and Dynamic
+  Finite Element Analysis*.
+- K. J. Bathe and S. Bolourchi, “Large Displacement Analysis of
+  Three-Dimensional Beam Structures,” 1979. FESA는 이 문헌의 선형화된 국부
+  좌표와 변환만 사용한다.
+- T. J. R. Hughes, R. L. Taylor, and W. Kanoknukulchai, “A Simple and
+  Efficient Finite Element for Plate Bending,” 1977. FESA는 selective
+  integration 원칙만 사용한다.
+- Abaqus 2024, “Choosing a Beam Element” 및 “BEAM GENERAL SECTION.”
+  FESA의 \(5A/6\) 생략값과 `SCF=0`은 Phase 1 명시적 계약이며 임의 일반 단면의
+  보편적인 Abaqus 기본값이 아니다. Abaqus B31 slenderness compensation도
+  구현하지 않는다.

include/fesa/elements/beam/beam3d2.hpp

@@ -0,0 +1,35 @@
+#pragma once
+
+#include <array>
+#include <optional>
+#include <vector>
+
+#include <fesa/core/diagnostic.hpp>
+#include <fesa/core/vec3.hpp>
+#include <fesa/fem/beam_frame.hpp>
+#include <fesa/model/beam_section.hpp>
+#include <fesa/model/material.hpp>
+
+namespace fesa {
+
+struct Beam3D2Input final {
+    std::array<Vec3, 2> coordinates;
+    IsotropicElastic material;
+    BeamSection section;
+};
+
+struct Beam3D2Contribution final {
+    Matrix12 local_stiffness;
+    Matrix12 global_stiffness;
+    BeamFrame frame;
+};
+
+struct BeamKernelResult final {
+    std::optional<Beam3D2Contribution> contribution;
+    std::vector<Diagnostic> diagnostics;
+};
+
+[[nodiscard]] BeamKernelResult compute_beam3d2(
+    const Beam3D2Input& input);
+
+}  // namespace fesa

src/fesa/elements/beam/beam3d2.cpp

@@ -0,0 +1,244 @@
+#include <fesa/elements/beam/beam3d2.hpp>
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <optional>
+#include <string>
+#include <utility>
+
+#include <fesa/fem/gauss_rule.hpp>
+#include <fesa/fem/line2_shape.hpp>
+
+namespace fesa {
+namespace {
+
+using StrainMatrix =
+    std::array<std::array<double, 12>, 6>;
+
+BeamKernelResult error_result(std::string code, std::string message) {
+    BeamKernelResult result;
+    result.diagnostics.push_back({
+        DiagnosticStage::model,
+        Severity::error,
+        std::move(code),
+        std::move(message),
+        std::nullopt,
+    });
+    return result;
+}
+
+bool is_positive_finite(const double value) {
+    return std::isfinite(value) && value > 0.0;
+}
+
+std::optional<BeamKernelResult> validate_properties(
+    const Beam3D2Input& input) {
+    if (!std::isfinite(input.material.young) ||
+        !std::isfinite(input.material.poisson)) {
+        return error_result(
+            "model.nonfinite_value",
+            "Beam kernel requires finite elastic constants.");
+    }
+    if (input.material.young <= 0.0 ||
+        input.material.poisson <= -1.0 ||
+        input.material.poisson >= 0.5) {
+        return error_result(
+            "model.invalid_material",
+            "Beam kernel requires E > 0 and -1 < nu < 0.5.");
+    }
+
+    if (!is_positive_finite(input.section.area) ||
+        !is_positive_finite(input.section.iy) ||
+        !is_positive_finite(input.section.iz) ||
+        !is_positive_finite(input.section.torsion_j) ||
+        !is_positive_finite(input.section.shear_area_y) ||
+        !is_positive_finite(input.section.shear_area_z)) {
+        return error_result(
+            "model.invalid_section",
+            "Beam kernel requires positive finite section properties.");
+    }
+    return std::nullopt;
+}
+
+StrainMatrix strain_matrix(
+    const double xi,
+    const double jacobian) {
+    StrainMatrix strain{};
+    const auto shape = line2_shape(xi);
+    const auto natural_derivative = line2_shape_derivative(xi);
+
+    for (std::size_t node = 0; node < shape.size(); ++node) {
+        const std::size_t offset = node * 6;
+        const double derivative = natural_derivative[node] / jacobian;
+        strain[0][offset] = derivative;
+        strain[1][offset + 1] = derivative;
+        strain[1][offset + 5] = -shape[node];
+        strain[2][offset + 2] = derivative;
+        strain[2][offset + 4] = shape[node];
+        strain[3][offset + 3] = derivative;
+        strain[4][offset + 4] = derivative;
+        strain[5][offset + 5] = derivative;
+    }
+    return strain;
+}
+
+template <std::size_t ComponentCount>
+void integrate_components(
+    Matrix12& stiffness,
+    const std::span<const GaussPoint1D> rule,
+    const std::array<std::size_t, ComponentCount>& components,
+    const std::array<double, 6>& constitutive,
+    const double jacobian) {
+    for (const GaussPoint1D& point : rule) {
+        const StrainMatrix strain = strain_matrix(point.xi, jacobian);
+        const double integration_weight = jacobian * point.weight;
+        for (std::size_t row = 0; row < stiffness.size(); ++row) {
+            for (std::size_t column = row;
+                 column < stiffness[row].size();
+                 ++column) {
+                double entry = 0.0;
+                for (const std::size_t component : components) {
+                    entry +=
+                        strain[component][row] *
+                        constitutive[component] *
+                        strain[component][column];
+                }
+                stiffness[row][column] +=
+                    entry * integration_weight;
+            }
+        }
+    }
+}
+
+void mirror_upper_triangle(Matrix12& matrix) {
+    for (std::size_t row = 0; row < matrix.size(); ++row) {
+        for (std::size_t column = row + 1;
+             column < matrix[row].size();
+             ++column) {
+            matrix[column][row] = matrix[row][column];
+        }
+    }
+}
+
+Matrix12 transform_stiffness(
+    const Matrix12& local,
+    const Matrix12& transformation) {
+    Matrix12 global{};
+    for (std::size_t row = 0; row < global.size(); ++row) {
+        for (std::size_t column = row;
+             column < global[row].size();
+             ++column) {
+            double entry = 0.0;
+            for (std::size_t local_row = 0;
+                 local_row < local.size();
+                 ++local_row) {
+                for (std::size_t local_column = 0;
+                     local_column < local[local_row].size();
+                     ++local_column) {
+                    entry +=
+                        transformation[local_row][row] *
+                        local[local_row][local_column] *
+                        transformation[local_column][column];
+                }
+            }
+            global[row][column] = entry;
+        }
+    }
+    mirror_upper_triangle(global);
+    return global;
+}
+
+bool is_finite(const Matrix12& matrix) {
+    for (const auto& row : matrix) {
+        for (const double value : row) {
+            if (!std::isfinite(value)) {
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+}  // namespace
+
+BeamKernelResult compute_beam3d2(const Beam3D2Input& input) {
+    if (const auto validation = validate_properties(input);
+        validation.has_value()) {
+        return *validation;
+    }
+
+    BeamFrameResult frame_result = make_beam_frame(
+        input.coordinates[0],
+        input.coordinates[1],
+        input.section.orientation);
+    if (!frame_result.frame.has_value()) {
+        return {std::nullopt, std::move(frame_result.diagnostics)};
+    }
+
+    const Vec3 axis{
+        input.coordinates[1].x - input.coordinates[0].x,
+        input.coordinates[1].y - input.coordinates[0].y,
+        input.coordinates[1].z - input.coordinates[0].z,
+    };
+    const double length = std::hypot(axis.x, axis.y, axis.z);
+    const double jacobian = length / 2.0;
+    if (!std::isfinite(jacobian) || jacobian <= 0.0) {
+        return error_result(
+            "model.zero_length_element",
+            "Beam kernel requires a representable positive Jacobian.");
+    }
+
+    const double shear_modulus =
+        input.material.young / (2.0 * (1.0 + input.material.poisson));
+    const std::array<double, 6> constitutive{
+        input.material.young * input.section.area,
+        shear_modulus * input.section.shear_area_y,
+        shear_modulus * input.section.shear_area_z,
+        shear_modulus * input.section.torsion_j,
+        input.material.young * input.section.iy,
+        input.material.young * input.section.iz,
+    };
+    for (const double value : constitutive) {
+        if (!is_positive_finite(value)) {
+            return error_result(
+                "model.nonfinite_value",
+                "Beam constitutive stiffness is not finite and positive.");
+        }
+    }
+
+    Matrix12 local{};
+    integrate_components(
+        local,
+        gauss_rule_1d(2),
+        std::array<std::size_t, 4>{0, 3, 4, 5},
+        constitutive,
+        jacobian);
+    integrate_components(
+        local,
+        gauss_rule_1d(1),
+        std::array<std::size_t, 2>{1, 2},
+        constitutive,
+        jacobian);
+    mirror_upper_triangle(local);
+
+    const Matrix12 transformation =
+        beam_transformation(*frame_result.frame);
+    Matrix12 global = transform_stiffness(local, transformation);
+    if (!is_finite(local) || !is_finite(global)) {
+        return error_result(
+            "model.nonfinite_value",
+            "Beam stiffness contains a nonfinite value.");
+    }
+
+    return {
+        Beam3D2Contribution{
+            std::move(local),
+            std::move(global),
+            *frame_result.frame,
+        },
+        {},
+    };
+}
+
+}  // namespace fesa

tests/CMakeLists.txt

@@ -234,3 +234,40 @@ add_test(
     COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
         --gtest_filter=BeamTransformation.*
 )
+
+add_executable(fesa_beam3d2_tests
+    unit/elements/beam3d2_test.cpp
+)
+
+target_compile_features(fesa_beam3d2_tests PRIVATE cxx_std_20)
+target_compile_options(
+    fesa_beam3d2_tests
+    PRIVATE
+        /W4
+        /permissive-
+        /EHsc
+)
+
+target_link_libraries(fesa_beam3d2_tests
+    PRIVATE
+        fesa_core
+        GTest::gtest_main
+)
+
+add_test(
+    NAME Beam3D2
+    COMMAND "$<TARGET_FILE:fesa_beam3d2_tests>"
+        --gtest_filter=Beam3D2.*
+)
+
+add_test(
+    NAME Timoshenko
+    COMMAND "$<TARGET_FILE:fesa_beam3d2_tests>"
+        --gtest_filter=Timoshenko.*
+)
+
+add_test(
+    NAME RigidBody
+    COMMAND "$<TARGET_FILE:fesa_beam3d2_tests>"
+        --gtest_filter=RigidBody.*
+)

tests/unit/elements/beam3d2_test.cpp

@@ -0,0 +1,372 @@
+#include <fesa/elements/beam/beam3d2.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+constexpr double kYoung = 210.0;
+constexpr double kPoisson = 0.3;
+constexpr double kArea = 0.4;
+constexpr double kIy = 0.03;
+constexpr double kIz = 0.05;
+constexpr double kTorsionJ = 0.02;
+constexpr double kShearAreaY = 0.25;
+constexpr double kShearAreaZ = 0.2;
+
+fesa::Beam3D2Input make_x_axis_input(const double length) {
+    return {
+        {{{0.0, 0.0, 0.0}, {length, 0.0, 0.0}}},
+        {fesa::MaterialId{1}, "elastic", kYoung, kPoisson},
+        {
+            fesa::SectionId{1},
+            "section",
+            kArea,
+            kIy,
+            kIz,
+            kTorsionJ,
+            kShearAreaY,
+            kShearAreaZ,
+            fesa::ShearPropertySource::input,
+            {0.0, 1.0, 0.0},
+            {},
+        },
+    };
+}
+
+double quadratic_form(
+    const fesa::Matrix12& matrix,
+    const std::array<double, 12>& vector) {
+    double result = 0.0;
+    for (std::size_t row = 0; row < matrix.size(); ++row) {
+        for (std::size_t column = 0; column < matrix[row].size(); ++column) {
+            result += vector[row] * matrix[row][column] * vector[column];
+        }
+    }
+    return result;
+}
+
+double strain_energy(
+    const fesa::Matrix12& matrix,
+    const std::array<double, 12>& vector) {
+    return 0.5 * quadratic_form(matrix, vector);
+}
+
+double maximum_abs_entry(const fesa::Matrix12& matrix) {
+    double maximum = 0.0;
+    for (const auto& row : matrix) {
+        for (const double value : row) {
+            maximum = std::max(maximum, std::abs(value));
+        }
+    }
+    return maximum;
+}
+
+double squared_norm(const std::array<double, 12>& vector) {
+    double result = 0.0;
+    for (const double value : vector) {
+        result += value * value;
+    }
+    return result;
+}
+
+double dot(const fesa::Vec3 first, const fesa::Vec3 second) {
+    return first.x * second.x + first.y * second.y + first.z * second.z;
+}
+
+fesa::Vec3 cross(const fesa::Vec3 first, const fesa::Vec3 second) {
+    return {
+        first.y * second.z - first.z * second.y,
+        first.z * second.x - first.x * second.z,
+        first.x * second.y - first.y * second.x,
+    };
+}
+
+fesa::Vec3 rotate_about_axis(
+    const fesa::Vec3 value,
+    const fesa::Vec3 unit_axis,
+    const double angle) {
+    const double cosine = std::cos(angle);
+    const double sine = std::sin(angle);
+    const fesa::Vec3 axis_cross_value = cross(unit_axis, value);
+    const double axis_projection = dot(unit_axis, value);
+    return {
+        value.x * cosine + axis_cross_value.x * sine +
+            unit_axis.x * axis_projection * (1.0 - cosine),
+        value.y * cosine + axis_cross_value.y * sine +
+            unit_axis.y * axis_projection * (1.0 - cosine),
+        value.z * cosine + axis_cross_value.z * sine +
+            unit_axis.z * axis_projection * (1.0 - cosine),
+    };
+}
+
+std::array<double, 12> rotate_dofs(
+    const std::array<double, 12>& values,
+    const fesa::Vec3 unit_axis,
+    const double angle) {
+    std::array<double, 12> rotated{};
+    for (const std::size_t offset : std::array<std::size_t, 4>{0, 3, 6, 9}) {
+        const fesa::Vec3 value{
+            values[offset],
+            values[offset + 1],
+            values[offset + 2],
+        };
+        const fesa::Vec3 rotated_value =
+            rotate_about_axis(value, unit_axis, angle);
+        rotated[offset] = rotated_value.x;
+        rotated[offset + 1] = rotated_value.y;
+        rotated[offset + 2] = rotated_value.z;
+    }
+    return rotated;
+}
+
+void expect_relative_near(
+    const double actual,
+    const double expected,
+    const double relative_tolerance =
+        512.0 * std::numeric_limits<double>::epsilon()) {
+    const double scale = std::max(1.0, std::abs(expected));
+    EXPECT_NEAR(actual, expected, relative_tolerance * scale);
+}
+
+TEST(Beam3D2, ProducesFiniteSymmetricLocalAndGlobalStiffness) {
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(2.5));
+
+    ASSERT_TRUE(result.contribution.has_value());
+    EXPECT_TRUE(result.diagnostics.empty());
+    const auto& contribution = *result.contribution;
+    for (std::size_t row = 0; row < 12; ++row) {
+        for (std::size_t column = 0; column < 12; ++column) {
+            EXPECT_TRUE(std::isfinite(
+                contribution.local_stiffness[row][column]));
+            EXPECT_TRUE(std::isfinite(
+                contribution.global_stiffness[row][column]));
+            EXPECT_DOUBLE_EQ(
+                contribution.local_stiffness[row][column],
+                contribution.local_stiffness[column][row]);
+            EXPECT_DOUBLE_EQ(
+                contribution.global_stiffness[row][column],
+                contribution.global_stiffness[column][row]);
+        }
+    }
+}
+
+TEST(RigidBody, SixIndependentModesHaveZeroStrainEnergy) {
+    constexpr double length = 3.0;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    const auto& stiffness = result.contribution->local_stiffness;
+
+    const std::array<std::array<double, 12>, 6> modes{{
+        {1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+        {0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
+         0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
+        {0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
+         0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
+        {0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+         0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
+        {0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
+         0.0, 0.0, -length, 0.0, 1.0, 0.0},
+        {0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
+         0.0, length, 0.0, 0.0, 0.0, 1.0},
+    }};
+
+    for (std::size_t mode = 0; mode < modes.size(); ++mode) {
+        const double roundoff_bound =
+            4096.0 * std::numeric_limits<double>::epsilon() *
+            maximum_abs_entry(stiffness) * squared_norm(modes[mode]);
+        EXPECT_NEAR(
+            quadratic_form(stiffness, modes[mode]),
+            0.0,
+            roundoff_bound)
+            << "rigid mode " << mode;
+    }
+}
+
+TEST(Timoshenko, ReproducesAnalyticalAxialSubmatrix) {
+    constexpr double length = 2.5;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    const auto& stiffness = result.contribution->local_stiffness;
+    const double expected = kYoung * kArea / length;
+
+    expect_relative_near(stiffness[0][0], expected);
+    expect_relative_near(stiffness[0][6], -expected);
+    expect_relative_near(stiffness[6][0], -expected);
+    expect_relative_near(stiffness[6][6], expected);
+}
+
+TEST(Timoshenko, ReproducesAnalyticalTorsionalSubmatrix) {
+    constexpr double length = 2.5;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    const auto& stiffness = result.contribution->local_stiffness;
+    const double shear_modulus = kYoung / (2.0 * (1.0 + kPoisson));
+    const double expected = shear_modulus * kTorsionJ / length;
+
+    expect_relative_near(stiffness[3][3], expected);
+    expect_relative_near(stiffness[3][9], -expected);
+    expect_relative_near(stiffness[9][3], -expected);
+    expect_relative_near(stiffness[9][9], expected);
+}
+
+TEST(Timoshenko, ReproducesConstantCurvatureEnergyAboutLocalY) {
+    constexpr double length = 2.5;
+    constexpr double curvature = 0.4;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    std::array<double, 12> displacement{};
+    displacement[4] = -0.5 * curvature * length;
+    displacement[10] = 0.5 * curvature * length;
+
+    const double expected =
+        0.5 * kYoung * kIy * curvature * curvature * length;
+    expect_relative_near(
+        strain_energy(result.contribution->local_stiffness, displacement),
+        expected);
+}
+
+TEST(Timoshenko, ReproducesConstantCurvatureEnergyAboutLocalZ) {
+    constexpr double length = 2.5;
+    constexpr double curvature = 0.4;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    std::array<double, 12> displacement{};
+    displacement[5] = -0.5 * curvature * length;
+    displacement[11] = 0.5 * curvature * length;
+
+    const double expected =
+        0.5 * kYoung * kIz * curvature * curvature * length;
+    expect_relative_near(
+        strain_energy(result.contribution->local_stiffness, displacement),
+        expected);
+}
+
+TEST(Timoshenko, ReproducesConstantShearEnergyInLocalY) {
+    constexpr double length = 2.5;
+    constexpr double shear_strain = 0.04;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    std::array<double, 12> displacement{};
+    displacement[1] = -0.5 * shear_strain * length;
+    displacement[7] = 0.5 * shear_strain * length;
+
+    const double shear_modulus = kYoung / (2.0 * (1.0 + kPoisson));
+    const double expected =
+        0.5 * shear_modulus * kShearAreaY *
+        shear_strain * shear_strain * length;
+    expect_relative_near(
+        strain_energy(result.contribution->local_stiffness, displacement),
+        expected);
+}
+
+TEST(Timoshenko, ReproducesConstantShearEnergyInLocalZ) {
+    constexpr double length = 2.5;
+    constexpr double shear_strain = 0.04;
+    const auto result = fesa::compute_beam3d2(make_x_axis_input(length));
+    ASSERT_TRUE(result.contribution.has_value());
+    std::array<double, 12> displacement{};
+    displacement[2] = -0.5 * shear_strain * length;
+    displacement[8] = 0.5 * shear_strain * length;
+
+    const double shear_modulus = kYoung / (2.0 * (1.0 + kPoisson));
+    const double expected =
+        0.5 * shear_modulus * kShearAreaZ *
+        shear_strain * shear_strain * length;
+    expect_relative_near(
+        strain_energy(result.contribution->local_stiffness, displacement),
+        expected);
+}
+
+TEST(Timoshenko, AvoidsShearLockingAcrossSlendernessSweep) {
+    constexpr double length = 2.0;
+    constexpr double curvature = 0.2;
+    for (const double slenderness : std::array{2.0, 10.0, 100.0, 1000.0}) {
+        fesa::Beam3D2Input input = make_x_axis_input(length);
+        const double side = length / slenderness;
+        input.section.area = side * side;
+        input.section.iy = side * side * side * side / 12.0;
+        input.section.iz = input.section.iy;
+        input.section.shear_area_y = 5.0 * input.section.area / 6.0;
+        input.section.shear_area_z = input.section.shear_area_y;
+
+        const auto result = fesa::compute_beam3d2(input);
+        ASSERT_TRUE(result.contribution.has_value());
+        std::array<double, 12> displacement{};
+        displacement[5] = -0.5 * curvature * length;
+        displacement[11] = 0.5 * curvature * length;
+
+        const double expected =
+            0.5 * kYoung * input.section.iz *
+            curvature * curvature * length;
+        const double relative_tolerance =
+            4096.0 * std::numeric_limits<double>::epsilon() *
+            slenderness * slenderness;
+        EXPECT_NEAR(
+            strain_energy(
+                result.contribution->local_stiffness,
+                displacement),
+            expected,
+            relative_tolerance * expected)
+            << "L/h=" << slenderness;
+    }
+}
+
+TEST(Beam3D2, PreservesGlobalEnergyUnderRigidCoordinateRotation) {
+    fesa::Beam3D2Input original = make_x_axis_input(2.5);
+    original.coordinates = {{
+        {1.0, -2.0, 0.5},
+        {3.0, -0.5, 1.25},
+    }};
+    original.section.orientation = {-1.0, 2.0, 3.0};
+
+    const double inverse_sqrt_fourteen = 1.0 / std::sqrt(14.0);
+    const fesa::Vec3 rotation_axis{
+        inverse_sqrt_fourteen,
+        2.0 * inverse_sqrt_fourteen,
+        3.0 * inverse_sqrt_fourteen,
+    };
+    constexpr double angle = 0.73;
+    fesa::Beam3D2Input rotated = original;
+    for (std::size_t node = 0; node < rotated.coordinates.size(); ++node) {
+        rotated.coordinates[node] = rotate_about_axis(
+            original.coordinates[node],
+            rotation_axis,
+            angle);
+    }
+    rotated.section.orientation = rotate_about_axis(
+        original.section.orientation,
+        rotation_axis,
+        angle);
+
+    const auto original_result = fesa::compute_beam3d2(original);
+    const auto rotated_result = fesa::compute_beam3d2(rotated);
+    ASSERT_TRUE(original_result.contribution.has_value());
+    ASSERT_TRUE(rotated_result.contribution.has_value());
+
+    const std::array<double, 12> original_displacement{
+        0.1, -0.3, 0.2, 0.04, -0.02, 0.03,
+        -0.2, 0.5, -0.1, -0.01, 0.06, -0.05,
+    };
+    const std::array<double, 12> rotated_displacement =
+        rotate_dofs(original_displacement, rotation_axis, angle);
+    const double original_energy = strain_energy(
+        original_result.contribution->global_stiffness,
+        original_displacement);
+    const double rotated_energy = strain_energy(
+        rotated_result.contribution->global_stiffness,
+        rotated_displacement);
+
+    expect_relative_near(
+        rotated_energy,
+        original_energy,
+        4096.0 * std::numeric_limits<double>::epsilon());
+}
+
+}  // namespace