feat(fem-and-beam-kernel): step 2 — beam-local-frame

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

include/fesa/fem/beam_frame.hpp

@@ -0,0 +1,32 @@
+#pragma once
+
+#include <array>
+#include <optional>
+#include <vector>
+
+#include <fesa/core/diagnostic.hpp>
+#include <fesa/core/vec3.hpp>
+
+namespace fesa {
+
+using Matrix12 = std::array<std::array<double, 12>, 12>;
+
+struct BeamFrame final {
+    Vec3 ex;
+    Vec3 ey;
+    Vec3 ez;
+};
+
+struct BeamFrameResult final {
+    std::optional<BeamFrame> frame;
+    std::vector<Diagnostic> diagnostics;
+};
+
+[[nodiscard]] BeamFrameResult make_beam_frame(
+    const Vec3& first,
+    const Vec3& second,
+    const Vec3& orientation);
+
+[[nodiscard]] Matrix12 beam_transformation(const BeamFrame& frame);
+
+}  // namespace fesa

src/fesa/fem/beam_frame.cpp

@@ -0,0 +1,146 @@
+#include <fesa/fem/beam_frame.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <string>
+#include <utility>
+
+namespace fesa {
+namespace {
+
+constexpr double kParallelTolerance =
+    64.0 * std::numeric_limits<double>::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<Vec3, 3> basis{frame.ex, frame.ey, frame.ez};
+    constexpr std::array<std::size_t, 4> 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

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 "$<TARGET_FILE:fesa_fem_primitives_tests>"
         --gtest_filter=EquationNumbering.*
 )
+
+add_test(
+    NAME BeamFrame
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=BeamFrame.*
+)
+
+add_test(
+    NAME BeamTransformation
+    COMMAND "$<TARGET_FILE:fesa_fem_primitives_tests>"
+        --gtest_filter=BeamTransformation.*
+)

tests/unit/fem/beam_frame_test.cpp

@@ -0,0 +1,278 @@
+#include <fesa/fem/beam_frame.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <limits>
+#include <string_view>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+constexpr double kUnitTolerance =
+    128.0 * std::numeric_limits<double>::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<double, 12> multiply(
+    const fesa::Matrix12& matrix,
+    const std::array<double, 12>& vector) {
+    std::array<double, 12> 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<double, 12>& 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<double>::epsilon());
+    expect_vec_near(
+        rotated.frame->ey,
+        rotate_about_axis(original.frame->ey, rotation_axis, angle),
+        256.0 * std::numeric_limits<double>::epsilon());
+    expect_vec_near(
+        rotated.frame->ez,
+        rotate_about_axis(original.frame->ez, rotation_axis, angle),
+        256.0 * std::numeric_limits<double>::epsilon());
+}
+
+TEST(BeamFrame, IsInvariantToFiniteOrientationScale) {
+    const double maximum = std::numeric_limits<double>::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<double>::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<double, 12> 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<double, 12>{
+            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<double, 12> 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<double, 12> local =
+        multiply(transformation, global);
+    const double global_energy = squared_norm(global);
+    const double tolerance =
+        256.0 * std::numeric_limits<double>::epsilon() *
+        std::max(1.0, global_energy);
+    EXPECT_NEAR(squared_norm(local), global_energy, tolerance);
+}
+
+}  // namespace