From 7ddd8e690e2201a0067869d21bf1e095bf67ca1d Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 31 Jul 2026 02:08:27 +0900 Subject: [PATCH] =?UTF-8?q?feat(fem-and-beam-kernel):=20step=202=20?= =?UTF-8?q?=E2=80=94=20beam-local-frame?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 1 + include/fesa/fem/beam_frame.hpp | 32 ++++ src/fesa/fem/beam_frame.cpp | 146 +++++++++++++++ tests/CMakeLists.txt | 13 ++ tests/unit/fem/beam_frame_test.cpp | 278 +++++++++++++++++++++++++++++ 5 files changed, 470 insertions(+) create mode 100644 include/fesa/fem/beam_frame.hpp create mode 100644 src/fesa/fem/beam_frame.cpp create mode 100644 tests/unit/fem/beam_frame_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 98f165b..5942524 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/fem/beam_frame.cpp src/fesa/fem/dof_manager.cpp src/fesa/fem/gauss_rule.cpp src/fesa/fem/line2_shape.cpp diff --git a/include/fesa/fem/beam_frame.hpp b/include/fesa/fem/beam_frame.hpp new file mode 100644 index 0000000..790af2f --- /dev/null +++ b/include/fesa/fem/beam_frame.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include +#include +#include + +#include +#include + +namespace fesa { + +using Matrix12 = std::array, 12>; + +struct BeamFrame final { + Vec3 ex; + Vec3 ey; + Vec3 ez; +}; + +struct BeamFrameResult final { + std::optional frame; + std::vector diagnostics; +}; + +[[nodiscard]] BeamFrameResult make_beam_frame( + const Vec3& first, + const Vec3& second, + const Vec3& orientation); + +[[nodiscard]] Matrix12 beam_transformation(const BeamFrame& frame); + +} // namespace fesa diff --git a/src/fesa/fem/beam_frame.cpp b/src/fesa/fem/beam_frame.cpp new file mode 100644 index 0000000..6006c23 --- /dev/null +++ b/src/fesa/fem/beam_frame.cpp @@ -0,0 +1,146 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace fesa { +namespace { + +constexpr double kParallelTolerance = + 64.0 * std::numeric_limits::epsilon(); + +double length(const Vec3 value) { + return std::hypot(value.x, value.y, value.z); +} + +double max_abs_component(const Vec3 value) { + return std::max( + std::abs(value.x), + std::max(std::abs(value.y), std::abs(value.z))); +} + +double dot(const Vec3 first, const Vec3 second) { + return first.x * second.x + first.y * second.y + first.z * second.z; +} + +Vec3 cross(const Vec3 first, const 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, + }; +} + +Vec3 normalized(const Vec3 value) { + const double scale = max_abs_component(value); + const Vec3 scaled{ + value.x / scale, + value.y / scale, + value.z / scale, + }; + const double scaled_length = length(scaled); + return { + scaled.x / scaled_length, + scaled.y / scaled_length, + scaled.z / scaled_length, + }; +} + +BeamFrameResult error_result(std::string code, std::string message) { + BeamFrameResult result; + result.diagnostics.push_back({ + DiagnosticStage::model, + Severity::error, + std::move(code), + std::move(message), + std::nullopt, + }); + return result; +} + +} // namespace + +BeamFrameResult make_beam_frame( + const Vec3& first, + const Vec3& second, + const Vec3& orientation) { + if (!is_finite(first) || !is_finite(second)) { + return error_result( + "model.nonfinite_value", + "Beam frame requires finite node coordinates."); + } + + const Vec3 axis{ + second.x - first.x, + second.y - first.y, + second.z - first.z, + }; + if (!is_finite(axis)) { + return error_result( + "model.nonfinite_value", + "Beam frame axis is not representable as a finite vector."); + } + + if (max_abs_component(axis) == 0.0) { + return error_result( + "model.zero_length_element", + "Beam frame requires distinct node coordinates."); + } + + if (!is_finite(orientation)) { + return error_result( + "model.nonfinite_value", + "Beam frame requires a finite orientation vector."); + } + + if (max_abs_component(orientation) == 0.0) { + return error_result( + "model.invalid_orientation", + "Beam frame requires a nonzero orientation vector."); + } + + const Vec3 ex = normalized(axis); + const Vec3 orientation_unit = normalized(orientation); + const double axial_projection = dot(orientation_unit, ex); + const Vec3 transverse{ + orientation_unit.x - axial_projection * ex.x, + orientation_unit.y - axial_projection * ex.y, + orientation_unit.z - axial_projection * ex.z, + }; + const double transverse_length = length(transverse); + if (transverse_length <= kParallelTolerance) { + return error_result( + "model.invalid_orientation", + "Beam orientation is parallel to the element axis."); + } + + const Vec3 ey = normalized(transverse); + const Vec3 ez = cross(ex, ey); + return { + BeamFrame{ex, ey, ez}, + {}, + }; +} + +Matrix12 beam_transformation(const BeamFrame& frame) { + Matrix12 transformation{}; + const std::array basis{frame.ex, frame.ey, frame.ez}; + constexpr std::array block_offsets{0, 3, 6, 9}; + + for (const std::size_t offset : block_offsets) { + for (std::size_t row = 0; row < basis.size(); ++row) { + transformation[offset + row][offset] = basis[row].x; + transformation[offset + row][offset + 1] = basis[row].y; + transformation[offset + row][offset + 2] = basis[row].z; + } + } + + return transformation; +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index be416c7..2dbd613 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -172,6 +172,7 @@ add_test( ) add_executable(fesa_fem_primitives_tests + unit/fem/beam_frame_test.cpp unit/fem/dof_manager_test.cpp unit/fem/gauss_rule_test.cpp unit/fem/line2_shape_test.cpp @@ -221,3 +222,15 @@ add_test( COMMAND "$" --gtest_filter=EquationNumbering.* ) + +add_test( + NAME BeamFrame + COMMAND "$" + --gtest_filter=BeamFrame.* +) + +add_test( + NAME BeamTransformation + COMMAND "$" + --gtest_filter=BeamTransformation.* +) diff --git a/tests/unit/fem/beam_frame_test.cpp b/tests/unit/fem/beam_frame_test.cpp new file mode 100644 index 0000000..5bd422f --- /dev/null +++ b/tests/unit/fem/beam_frame_test.cpp @@ -0,0 +1,278 @@ +#include + +#include +#include +#include +#include +#include + +#include + +namespace { + +constexpr double kUnitTolerance = + 128.0 * std::numeric_limits::epsilon(); + +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), + }; +} + +void expect_vec_near( + const fesa::Vec3 actual, + const fesa::Vec3 expected, + const double tolerance = kUnitTolerance) { + EXPECT_NEAR(actual.x, expected.x, tolerance); + EXPECT_NEAR(actual.y, expected.y, tolerance); + EXPECT_NEAR(actual.z, expected.z, tolerance); +} + +void expect_error( + const fesa::BeamFrameResult& result, + const std::string_view code) { + EXPECT_FALSE(result.frame.has_value()); + ASSERT_EQ(result.diagnostics.size(), 1); + EXPECT_EQ(result.diagnostics[0].stage, fesa::DiagnosticStage::model); + EXPECT_EQ(result.diagnostics[0].severity, fesa::Severity::error); + EXPECT_EQ(result.diagnostics[0].code, code); +} + +std::array multiply( + const fesa::Matrix12& matrix, + const std::array& vector) { + std::array product{}; + for (std::size_t row = 0; row < product.size(); ++row) { + for (std::size_t column = 0; column < vector.size(); ++column) { + product[row] += matrix[row][column] * vector[column]; + } + } + return product; +} + +double squared_norm(const std::array& vector) { + double result = 0.0; + for (const double value : vector) { + result += value * value; + } + return result; +} + +TEST(BeamFrame, AlignsLocalAxesForGlobalXAxisBeam) { + const auto result = fesa::make_beam_frame( + {0.0, 0.0, 0.0}, + {2.0, 0.0, 0.0}, + {0.0, 3.0, 0.0}); + + ASSERT_TRUE(result.frame.has_value()); + EXPECT_TRUE(result.diagnostics.empty()); + expect_vec_near(result.frame->ex, {1.0, 0.0, 0.0}); + expect_vec_near(result.frame->ey, {0.0, 1.0, 0.0}); + expect_vec_near(result.frame->ez, {0.0, 0.0, 1.0}); +} + +TEST(BeamFrame, ProducesRightHandedOrthonormalBasisForArbitraryAxis) { + const auto result = fesa::make_beam_frame( + {1.0, 2.0, 3.0}, + {3.0, 4.0, 4.0}, + {0.0, 1.0, 0.0}); + + ASSERT_TRUE(result.frame.has_value()); + const fesa::BeamFrame& frame = *result.frame; + const double inverse_sqrt_five = 1.0 / std::sqrt(5.0); + expect_vec_near(frame.ex, {2.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0}); + expect_vec_near( + frame.ey, + { + -4.0 * inverse_sqrt_five / 3.0, + 5.0 * inverse_sqrt_five / 3.0, + -2.0 * inverse_sqrt_five / 3.0, + }); + expect_vec_near( + frame.ez, + {-inverse_sqrt_five, 0.0, 2.0 * inverse_sqrt_five}); + + EXPECT_NEAR(dot(frame.ex, frame.ex), 1.0, kUnitTolerance); + EXPECT_NEAR(dot(frame.ey, frame.ey), 1.0, kUnitTolerance); + EXPECT_NEAR(dot(frame.ez, frame.ez), 1.0, kUnitTolerance); + EXPECT_NEAR(dot(frame.ex, frame.ey), 0.0, kUnitTolerance); + EXPECT_NEAR(dot(frame.ey, frame.ez), 0.0, kUnitTolerance); + EXPECT_NEAR(dot(frame.ez, frame.ex), 0.0, kUnitTolerance); + EXPECT_NEAR(dot(frame.ex, cross(frame.ey, frame.ez)), 1.0, kUnitTolerance); +} + +TEST(BeamFrame, RotatesCovariantlyWithBeamGeometry) { + constexpr fesa::Vec3 first{1.0, -2.0, 0.5}; + constexpr fesa::Vec3 second{4.0, 1.0, 2.5}; + constexpr fesa::Vec3 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; + + const auto original = + fesa::make_beam_frame(first, second, orientation); + const auto rotated = fesa::make_beam_frame( + rotate_about_axis(first, rotation_axis, angle), + rotate_about_axis(second, rotation_axis, angle), + rotate_about_axis(orientation, rotation_axis, angle)); + + ASSERT_TRUE(original.frame.has_value()); + ASSERT_TRUE(rotated.frame.has_value()); + expect_vec_near( + rotated.frame->ex, + rotate_about_axis(original.frame->ex, rotation_axis, angle), + 256.0 * std::numeric_limits::epsilon()); + expect_vec_near( + rotated.frame->ey, + rotate_about_axis(original.frame->ey, rotation_axis, angle), + 256.0 * std::numeric_limits::epsilon()); + expect_vec_near( + rotated.frame->ez, + rotate_about_axis(original.frame->ez, rotation_axis, angle), + 256.0 * std::numeric_limits::epsilon()); +} + +TEST(BeamFrame, IsInvariantToFiniteOrientationScale) { + const double maximum = std::numeric_limits::max(); + const auto unit_scale = fesa::make_beam_frame( + {0.0, 0.0, 0.0}, + {1.0, 1.0, 1.0}, + {1.0, -1.0, 0.0}); + const auto maximum_scale = fesa::make_beam_frame( + {0.0, 0.0, 0.0}, + {1.0, 1.0, 1.0}, + {maximum, -maximum, 0.0}); + + ASSERT_TRUE(unit_scale.frame.has_value()); + ASSERT_TRUE(maximum_scale.frame.has_value()); + expect_vec_near(maximum_scale.frame->ex, unit_scale.frame->ex); + expect_vec_near(maximum_scale.frame->ey, unit_scale.frame->ey); + expect_vec_near(maximum_scale.frame->ez, unit_scale.frame->ez); +} + +TEST(BeamFrame, RejectsZeroLengthAxis) { + expect_error( + fesa::make_beam_frame( + {1.0, 2.0, 3.0}, + {1.0, 2.0, 3.0}, + {0.0, 1.0, 0.0}), + "model.zero_length_element"); +} + +TEST(BeamFrame, RejectsZeroOrientation) { + expect_error( + fesa::make_beam_frame( + {0.0, 0.0, 0.0}, + {1.0, 0.0, 0.0}, + {0.0, 0.0, 0.0}), + "model.invalid_orientation"); +} + +TEST(BeamFrame, RejectsOrientationParallelToElementAxis) { + expect_error( + fesa::make_beam_frame( + {0.0, 0.0, 0.0}, + {2.0, 0.0, 0.0}, + {5.0, 0.0, 0.0}), + "model.invalid_orientation"); +} + +TEST(BeamFrame, RejectsNumericallyParallelOrientationAtLargeScale) { + constexpr double transverse_ratio = + 32.0 * std::numeric_limits::epsilon(); + expect_error( + fesa::make_beam_frame( + {0.0, 0.0, 0.0}, + {2.0, 0.0, 0.0}, + {1.0e200, transverse_ratio * 1.0e200, 0.0}), + "model.invalid_orientation"); +} + +TEST(BeamTransformation, IsIdentityForGlobalAlignedFrame) { + const fesa::Matrix12 transformation = fesa::beam_transformation({ + {1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}, + }); + + for (std::size_t row = 0; row < transformation.size(); ++row) { + for (std::size_t column = 0; + column < transformation[row].size(); + ++column) { + EXPECT_DOUBLE_EQ( + transformation[row][column], + row == column ? 1.0 : 0.0); + } + } +} + +TEST(BeamTransformation, MapsGlobalComponentsIntoLocalDofOrder) { + const fesa::Matrix12 transformation = fesa::beam_transformation({ + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}, + {1.0, 0.0, 0.0}, + }); + const std::array global{ + 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, + 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, + }; + + EXPECT_EQ( + multiply(transformation, global), + (std::array{ + 2.0, 3.0, 1.0, 5.0, 6.0, 4.0, + 8.0, 9.0, 7.0, 11.0, 12.0, 10.0, + })); +} + +TEST(BeamTransformation, PreservesVectorEnergyForArbitraryFrame) { + const auto frame_result = fesa::make_beam_frame( + {-2.0, 1.0, 3.0}, + {1.0, 5.0, 5.0}, + {2.0, -1.0, 4.0}); + ASSERT_TRUE(frame_result.frame.has_value()); + const fesa::Matrix12 transformation = + fesa::beam_transformation(*frame_result.frame); + const std::array global{ + 3.0, -1.0, 2.0, 0.5, 4.0, -2.0, + -3.0, 6.0, 1.5, 2.5, -0.25, 5.0, + }; + + const std::array local = + multiply(transformation, global); + const double global_energy = squared_norm(global); + const double tolerance = + 256.0 * std::numeric_limits::epsilon() * + std::max(1.0, global_energy); + EXPECT_NEAR(squared_norm(local), global_energy, tolerance); +} + +} // namespace