diff --git a/CMakeLists.txt b/CMakeLists.txt index 5942524..cb1e4a4 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/docs/formulation/timoshenko-beam-3d.md b/docs/formulation/timoshenko-beam-3d.md new file mode 100644 index 0000000..0710bfc --- /dev/null +++ b/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도 + 구현하지 않는다. diff --git a/include/fesa/elements/beam/beam3d2.hpp b/include/fesa/elements/beam/beam3d2.hpp new file mode 100644 index 0000000..265f8ce --- /dev/null +++ b/include/fesa/elements/beam/beam3d2.hpp @@ -0,0 +1,35 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace fesa { + +struct Beam3D2Input final { + std::array coordinates; + IsotropicElastic material; + BeamSection section; +}; + +struct Beam3D2Contribution final { + Matrix12 local_stiffness; + Matrix12 global_stiffness; + BeamFrame frame; +}; + +struct BeamKernelResult final { + std::optional contribution; + std::vector diagnostics; +}; + +[[nodiscard]] BeamKernelResult compute_beam3d2( + const Beam3D2Input& input); + +} // namespace fesa diff --git a/src/fesa/elements/beam/beam3d2.cpp b/src/fesa/elements/beam/beam3d2.cpp new file mode 100644 index 0000000..5b98082 --- /dev/null +++ b/src/fesa/elements/beam/beam3d2.cpp @@ -0,0 +1,244 @@ +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace fesa { +namespace { + +using StrainMatrix = + std::array, 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 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 +void integrate_components( + Matrix12& stiffness, + const std::span rule, + const std::array& components, + const std::array& 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 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{0, 3, 4, 5}, + constitutive, + jacobian); + integrate_components( + local, + gauss_rule_1d(1), + std::array{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 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2dbd613..e5cebd8 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -234,3 +234,40 @@ add_test( COMMAND "$" --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 "$" + --gtest_filter=Beam3D2.* +) + +add_test( + NAME Timoshenko + COMMAND "$" + --gtest_filter=Timoshenko.* +) + +add_test( + NAME RigidBody + COMMAND "$" + --gtest_filter=RigidBody.* +) diff --git a/tests/unit/elements/beam3d2_test.cpp b/tests/unit/elements/beam3d2_test.cpp new file mode 100644 index 0000000..98b2dd5 --- /dev/null +++ b/tests/unit/elements/beam3d2_test.cpp @@ -0,0 +1,372 @@ +#include + +#include +#include +#include +#include +#include + +#include + +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& 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& 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& 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 rotate_dofs( + const std::array& values, + const fesa::Vec3 unit_axis, + const double angle) { + std::array rotated{}; + for (const std::size_t offset : std::array{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::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, 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::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 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 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 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 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 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::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 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 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::epsilon()); +} + +} // namespace