feat(euler-beam-3d): add local beam kernel

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"

src/fesa/elements/euler_beam_3d.cpp

@@ -0,0 +1,113 @@
+#include <fesa/elements/euler_beam_3d.hpp>
+
+#include <cmath>
+#include <stdexcept>
+
+namespace fesa::elements {
+namespace {
+
+constexpr int matrix_size = 12;
+
+std::size_t index(int row, int column)
+{
+    return static_cast<std::size_t>(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<std::size_t>(row)] +=
+                stiffness[index(row, column)] * local_displacements[static_cast<std::size_t>(column)];
+        }
+    }
+    return forces;
+}
+
+} // namespace fesa::elements

src/fesa/elements/euler_beam_3d.hpp

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <array>
+
+namespace fesa::elements {
+
+using Vector12 = std::array<double, 12>;
+using Matrix12 = std::array<double, 144>;
+
+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

tests/unit/euler_beam_3d_local_stiffness_test.cpp

@@ -0,0 +1,125 @@
+#include <fesa/elements/euler_beam_3d.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <stdexcept>
+
+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<std::size_t>(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<std::size_t>(row)] +=
+                entry(matrix, row, column) * vector[static_cast<std::size_t>(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;
+}