From 95ca95180a09bb5fe2cebd172f7132581eb1cce6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=EA=B9=80=EA=B2=BD=EC=A2=85?= Date: Fri, 12 Jun 2026 18:02:05 +0900 Subject: [PATCH] feat(euler-beam-3d): add local beam kernel --- phases/euler-beam-3d/index.json | 3 +- src/fesa/elements/euler_beam_3d.cpp | 113 ++++++++++++++++ src/fesa/elements/euler_beam_3d.hpp | 25 ++++ .../euler_beam_3d_local_stiffness_test.cpp | 125 ++++++++++++++++++ 4 files changed, 265 insertions(+), 1 deletion(-) create mode 100644 src/fesa/elements/euler_beam_3d.cpp create mode 100644 src/fesa/elements/euler_beam_3d.hpp create mode 100644 tests/unit/euler_beam_3d_local_stiffness_test.cpp diff --git a/phases/euler-beam-3d/index.json b/phases/euler-beam-3d/index.json index 3bd9e07..e804418 100644 --- a/phases/euler-beam-3d/index.json +++ b/phases/euler-beam-3d/index.json @@ -73,7 +73,8 @@ { "step": 7, "name": "local-stiffness-kernel", - "status": "pending", + "status": "completed", + "summary": "local 3D Euler beam stiffness and local end-force kernel added", "allowed_paths": [ "src/fesa/elements/", "tests/unit/euler_beam_3d_*_test.cpp" diff --git a/src/fesa/elements/euler_beam_3d.cpp b/src/fesa/elements/euler_beam_3d.cpp new file mode 100644 index 0000000..8bd4337 --- /dev/null +++ b/src/fesa/elements/euler_beam_3d.cpp @@ -0,0 +1,113 @@ +#include + +#include +#include + +namespace fesa::elements { +namespace { + +constexpr int matrix_size = 12; + +std::size_t index(int row, int column) +{ + return static_cast(row * matrix_size + column); +} + +void require_positive_finite(double value, const char* name) +{ + if (!std::isfinite(value) || value <= 0.0) { + throw std::invalid_argument(name); + } +} + +void validate_section(double length, const EulerBeam3DSection& section) +{ + require_positive_finite(length, "Euler beam length must be positive and finite"); + require_positive_finite(section.young_modulus, "Euler beam Young modulus must be positive and finite"); + require_positive_finite(section.shear_modulus, "Euler beam shear modulus must be positive and finite"); + require_positive_finite(section.area, "Euler beam area must be positive and finite"); + require_positive_finite(section.torsion_constant, "Euler beam torsion constant must be positive and finite"); + require_positive_finite(section.second_moment_y, "Euler beam Iy must be positive and finite"); + require_positive_finite(section.second_moment_z, "Euler beam Iz must be positive and finite"); +} + +void set_symmetric(Matrix12& matrix, int row, int column, double value) +{ + matrix[index(row, column)] = value; + matrix[index(column, row)] = value; +} + +} // namespace + +Matrix12 euler_beam_3d_local_stiffness(double length, const EulerBeam3DSection& section) +{ + validate_section(length, section); + + const double length2 = length * length; + const double length3 = length2 * length; + const double axial = section.young_modulus * section.area / length; + const double torsion = section.shear_modulus * section.torsion_constant / length; + + const double eiy = section.young_modulus * section.second_moment_y; + const double cy1 = 12.0 * eiy / length3; + const double cy2 = 6.0 * eiy / length2; + const double cy3 = 4.0 * eiy / length; + const double cy4 = 2.0 * eiy / length; + + const double eiz = section.young_modulus * section.second_moment_z; + const double cz1 = 12.0 * eiz / length3; + const double cz2 = 6.0 * eiz / length2; + const double cz3 = 4.0 * eiz / length; + const double cz4 = 2.0 * eiz / length; + + Matrix12 stiffness{}; + + set_symmetric(stiffness, 0, 0, axial); + set_symmetric(stiffness, 0, 6, -axial); + set_symmetric(stiffness, 6, 6, axial); + + set_symmetric(stiffness, 3, 3, torsion); + set_symmetric(stiffness, 3, 9, -torsion); + set_symmetric(stiffness, 9, 9, torsion); + + set_symmetric(stiffness, 1, 1, cz1); + set_symmetric(stiffness, 1, 5, cz2); + set_symmetric(stiffness, 1, 7, -cz1); + set_symmetric(stiffness, 1, 11, cz2); + set_symmetric(stiffness, 5, 5, cz3); + set_symmetric(stiffness, 5, 7, -cz2); + set_symmetric(stiffness, 5, 11, cz4); + set_symmetric(stiffness, 7, 7, cz1); + set_symmetric(stiffness, 7, 11, -cz2); + set_symmetric(stiffness, 11, 11, cz3); + + set_symmetric(stiffness, 2, 2, cy1); + set_symmetric(stiffness, 2, 4, -cy2); + set_symmetric(stiffness, 2, 8, -cy1); + set_symmetric(stiffness, 2, 10, -cy2); + set_symmetric(stiffness, 4, 4, cy3); + set_symmetric(stiffness, 4, 8, cy2); + set_symmetric(stiffness, 4, 10, cy4); + set_symmetric(stiffness, 8, 8, cy1); + set_symmetric(stiffness, 8, 10, cy2); + set_symmetric(stiffness, 10, 10, cy3); + + return stiffness; +} + +Vector12 euler_beam_3d_local_end_forces(double length, + const EulerBeam3DSection& section, + const Vector12& local_displacements) +{ + const auto stiffness = euler_beam_3d_local_stiffness(length, section); + Vector12 forces{}; + for (int row = 0; row < matrix_size; ++row) { + for (int column = 0; column < matrix_size; ++column) { + forces[static_cast(row)] += + stiffness[index(row, column)] * local_displacements[static_cast(column)]; + } + } + return forces; +} + +} // namespace fesa::elements diff --git a/src/fesa/elements/euler_beam_3d.hpp b/src/fesa/elements/euler_beam_3d.hpp new file mode 100644 index 0000000..94af430 --- /dev/null +++ b/src/fesa/elements/euler_beam_3d.hpp @@ -0,0 +1,25 @@ +#pragma once + +#include + +namespace fesa::elements { + +using Vector12 = std::array; +using Matrix12 = std::array; + +struct EulerBeam3DSection { + double young_modulus; + double shear_modulus; + double area; + double torsion_constant; + double second_moment_y; + double second_moment_z; +}; + +Matrix12 euler_beam_3d_local_stiffness(double length, const EulerBeam3DSection& section); + +Vector12 euler_beam_3d_local_end_forces(double length, + const EulerBeam3DSection& section, + const Vector12& local_displacements); + +} // namespace fesa::elements diff --git a/tests/unit/euler_beam_3d_local_stiffness_test.cpp b/tests/unit/euler_beam_3d_local_stiffness_test.cpp new file mode 100644 index 0000000..5c2f1dc --- /dev/null +++ b/tests/unit/euler_beam_3d_local_stiffness_test.cpp @@ -0,0 +1,125 @@ +#include + +#include +#include +#include + +namespace { + +constexpr double tolerance = 1.0e-10; + +bool close(double actual, double expected, double tol = tolerance) +{ + return std::abs(actual - expected) <= tol; +} + +double entry(const fesa::elements::Matrix12& matrix, int row, int column) +{ + return matrix[static_cast(row * 12 + column)]; +} + +fesa::elements::Vector12 multiply(const fesa::elements::Matrix12& matrix, + const fesa::elements::Vector12& vector) +{ + fesa::elements::Vector12 result{}; + for (int row = 0; row < 12; ++row) { + for (int column = 0; column < 12; ++column) { + result[static_cast(row)] += + entry(matrix, row, column) * vector[static_cast(column)]; + } + } + return result; +} + +bool throws_invalid_length(const fesa::elements::EulerBeam3DSection& section) +{ + try { + (void)fesa::elements::euler_beam_3d_local_stiffness(0.0, section); + } catch (const std::invalid_argument&) { + return true; + } + return false; +} + +bool throws_invalid_section(fesa::elements::EulerBeam3DSection section) +{ + section.second_moment_z = 0.0; + try { + (void)fesa::elements::euler_beam_3d_local_stiffness(2.0, section); + } catch (const std::invalid_argument&) { + return true; + } + return false; +} + +} // namespace + +int main() +{ + const double length = 2.0; + const fesa::elements::EulerBeam3DSection section{ + 210.0, + 80.0, + 3.0, + 4.0, + 5.0, + 6.0 + }; + + const auto stiffness = fesa::elements::euler_beam_3d_local_stiffness(length, section); + + if (!close(entry(stiffness, 0, 0), 315.0) || + !close(entry(stiffness, 0, 6), -315.0) || + !close(entry(stiffness, 6, 6), 315.0)) { + return 1; + } + + if (!close(entry(stiffness, 3, 3), 160.0) || + !close(entry(stiffness, 3, 9), -160.0) || + !close(entry(stiffness, 9, 9), 160.0)) { + return 2; + } + + if (!close(entry(stiffness, 1, 1), 1890.0) || + !close(entry(stiffness, 1, 5), 1890.0) || + !close(entry(stiffness, 5, 5), 2520.0)) { + return 3; + } + + if (!close(entry(stiffness, 2, 2), 1575.0) || + !close(entry(stiffness, 2, 4), -1575.0) || + !close(entry(stiffness, 4, 4), 2100.0)) { + return 4; + } + + for (int row = 0; row < 12; ++row) { + for (int column = 0; column < 12; ++column) { + if (!close(entry(stiffness, row, column), entry(stiffness, column, row))) { + return 5; + } + } + } + + fesa::elements::Vector12 displacement{}; + displacement[0] = 0.1; + displacement[5] = 0.2; + displacement[8] = -0.05; + + const auto expected_forces = multiply(stiffness, displacement); + const auto recovered_forces = + fesa::elements::euler_beam_3d_local_end_forces(length, section, displacement); + for (std::size_t index = 0; index < expected_forces.size(); ++index) { + if (!close(recovered_forces[index], expected_forces[index])) { + return 6; + } + } + + if (!throws_invalid_length(section)) { + return 7; + } + if (!throws_invalid_section(section)) { + return 8; + } + + return 0; +}