Merge pull request #321 from Exawind/linear-springs

Implement R3-based geometrically nonlinear, constitutively linear spring element

src/elements/springs/create_springs.hpp

@@ -0,0 +1,38 @@
+#pragma once
+
+#include "springs.hpp"
+#include "springs_input.hpp"
+
+namespace openturbine {
+
+inline Springs CreateSprings(const SpringsInput& springs_input) {
+    Springs springs(springs_input.NumElements());
+
+    auto host_node_state_indices = Kokkos::create_mirror(springs.node_state_indices);
+    auto host_x0 = Kokkos::create_mirror(springs.x0);
+    auto host_l_ref = Kokkos::create_mirror(springs.l_ref);
+    auto host_k = Kokkos::create_mirror(springs.k);
+
+    for (size_t i_elem = 0; i_elem < springs_input.NumElements(); i_elem++) {
+        const auto& element = springs_input.elements[i_elem];
+
+        host_node_state_indices(i_elem, 0U) = static_cast<size_t>(element.nodes[0].ID);
+        host_node_state_indices(i_elem, 1U) = static_cast<size_t>(element.nodes[1].ID);
+
+        for (size_t i = 0; i < 3; i++) {
+            host_x0(i_elem, i) = element.nodes[1].x[i] - element.nodes[0].x[i];
+        }
+
+        host_l_ref(i_elem) = element.undeformed_length;
+        host_k(i_elem) = element.stiffness;
+    }
+
+    Kokkos::deep_copy(springs.node_state_indices, host_node_state_indices);
+    Kokkos::deep_copy(springs.x0, host_x0);
+    Kokkos::deep_copy(springs.l_ref, host_l_ref);
+    Kokkos::deep_copy(springs.k, host_k);
+
+    return springs;
+}
+
+}  // namespace openturbine

src/elements/springs/spring_element.hpp

@@ -0,0 +1,32 @@
+#pragma once
+
+#include <array>
+#include <cmath>
+
+#include "src/model/node.hpp"
+
+namespace openturbine {
+
+/**
+ * @brief Spring element represents a constitutively linear spring connecting two nodes and defined
+ * by its scalar stiffness and undeformed length.
+ */
+struct SpringElement {
+    std::array<Node, 2> nodes;  // 2 nodes (start and end points of spring)
+    double stiffness;           // Spring stiffness coefficient
+    double undeformed_length;   // Reference/undeformed length of spring
+
+    SpringElement(std::array<Node, 2> n, double k, double l0 = 0.)
+        : nodes(n), stiffness(k), undeformed_length(l0) {
+        // If undeformed length is not provided, compute it from the nodes
+        if (l0 == 0.) {
+            const auto& p1 = nodes[0].x;
+            const auto& p2 = nodes[1].x;
+            undeformed_length = std::sqrt(
+                std::pow(p2[0] - p1[0], 2) + std::pow(p2[1] - p1[1], 2) + std::pow(p2[2] - p1[2], 2)
+            );
+        }
+    }
+};
+
+}  // namespace openturbine

src/elements/springs/springs.hpp

@@ -0,0 +1,65 @@
+#pragma once
+
+#include <Kokkos_Core.hpp>
+
+#include "src/dof_management/freedom_signature.hpp"
+#include "src/types.hpp"
+
+namespace openturbine {
+
+/**
+ * @brief Contains field variables for spring elements to compute per-element
+ * contributions to the residual vector and system/iteration matrix
+ */
+struct Springs {
+    size_t num_elems;  //< Total number of elements
+
+    Kokkos::View<size_t*> num_nodes_per_element;  //< This is always 2 for springs
+    Kokkos::View<size_t* [2]> node_state_indices;
+    Kokkos::View<FreedomSignature* [2]> element_freedom_signature;
+    Kokkos::View<size_t* [2][3]> element_freedom_table;  //< Only translational DOFs for springs
+
+    Kokkos::View<double* [3]> x0;          //< Initial distance vector between nodes
+    Kokkos::View<double* [3]> u1;          //< Displacement of node 1
+    Kokkos::View<double* [3]> u2;          //< Displacement of node 2
+    Kokkos::View<double* [3]> r;           //< Current distance vector between nodes
+    Kokkos::View<double*> l;               //< Current length of springs
+    Kokkos::View<double*> l_ref;           //< Initial length of springs
+    Kokkos::View<double*> k;               //< Spring stiffness coefficients
+    Kokkos::View<double*> c1;              //< First coefficient for force calculation
+    Kokkos::View<double*> c2;              //< Second coefficient for force calculation
+    Kokkos::View<double* [3]> f;           //< Force components
+    Kokkos::View<double* [3][3]> a;        //< Stiffness matrices
+    Kokkos::View<double* [3][3]> r_tilde;  //< Skew-symmetric matrix of r
+
+    Kokkos::View<double* [3]> residual_vector_terms;
+    Kokkos::View<double* [2][2][3][3]> stiffness_matrix_terms;
+
+    Springs(const size_t n_spring_elems)
+        : num_elems(n_spring_elems),
+          num_nodes_per_element("num_nodes_per_element", num_elems),
+          node_state_indices("node_state_indices", num_elems),
+          element_freedom_signature("element_freedom_signature", num_elems),
+          element_freedom_table("element_freedom_table", num_elems),
+          x0("x0", num_elems),
+          u1("u1", num_elems),
+          u2("u2", num_elems),
+          r("r", num_elems),
+          l("l", num_elems),
+          l_ref("l_ref", num_elems),
+          k("k", num_elems),
+          c1("c1", num_elems),
+          c2("c2", num_elems),
+          f("f", num_elems),
+          a("a", num_elems),
+          r_tilde("r_tilde", num_elems),
+          residual_vector_terms("residual_vector_terms", num_elems),
+          stiffness_matrix_terms("stiffness_matrix_terms", num_elems) {
+        Kokkos::deep_copy(num_nodes_per_element, 2);  // Always 2 nodes per element
+        Kokkos::deep_copy(
+            element_freedom_signature, FreedomSignature::JustPosition
+        );  // Springs only have translational DOFs
+    }
+};
+
+}  // namespace openturbine

src/elements/springs/springs_input.hpp

@@ -0,0 +1,19 @@
+#pragma once
+
+#include "spring_element.hpp"
+
+namespace openturbine {
+
+/**
+ * @brief Represents the input data for creating spring elements
+ */
+struct SpringsInput {
+    std::vector<SpringElement> elements;  //< All spring elements in the system
+
+    SpringsInput(std::vector<SpringElement> elems) : elements(std::move(elems)) {}
+
+    /// Returns the total number of spring elements in the system
+    [[nodiscard]] size_t NumElements() const { return elements.size(); }
+};
+
+}  // namespace openturbine

src/system/CMakeLists.txt

@@ -1,2 +1,5 @@
+target_sources(openturbine_library PRIVATE)
+
+add_subdirectory(springs)
 add_subdirectory(masses)
 add_subdirectory(beams)

src/system/springs/CMakeLists.txt

@@ -0,0 +1 @@
+target_sources(openturbine_library PRIVATE)

src/system/springs/calculate_distance_components.hpp

@@ -0,0 +1,29 @@
+#pragma once
+
+#include <KokkosBlas.hpp>
+#include <Kokkos_Core.hpp>
+
+#include "src/math/vector_operations.hpp"
+
+namespace openturbine {
+
+struct CalculateDistanceComponents {
+    Kokkos::View<double* [3]>::const_type x0_;
+    Kokkos::View<double* [3]>::const_type u1_;
+    Kokkos::View<double* [3]>::const_type u2_;
+    Kokkos::View<double* [3]> r_;
+
+    KOKKOS_FUNCTION
+    void operator()(int i_elem) const {
+        auto x0 = Kokkos::subview(x0_, i_elem, Kokkos::ALL);
+        auto u1 = Kokkos::subview(u1_, i_elem, Kokkos::ALL);
+        auto u2 = Kokkos::subview(u2_, i_elem, Kokkos::ALL);
+        auto r = Kokkos::subview(r_, i_elem, Kokkos::ALL);
+
+        for (int i = 0; i < 3; ++i) {
+            r(i) = x0(i) + u2(i) - u1(i);
+        }
+    }
+};
+
+}  // namespace openturbine

src/system/springs/calculate_force_coefficients.hpp

@@ -0,0 +1,26 @@
+#pragma once
+
+#include <KokkosBlas.hpp>
+#include <Kokkos_Core.hpp>
+
+#include "src/math/vector_operations.hpp"
+
+namespace openturbine {
+
+struct CalculateForceCoefficients {
+    Kokkos::View<double*>::const_type k_;
+    Kokkos::View<double*>::const_type l_ref_;
+    Kokkos::View<double*>::const_type l_;
+    Kokkos::View<double*> c1_;
+    Kokkos::View<double*> c2_;
+
+    KOKKOS_FUNCTION
+    void operator()(int i_elem) const {
+        // c1 = k * (l_ref/l - 1)
+        c1_(i_elem) = k_(i_elem) * (l_ref_(i_elem) / l_(i_elem) - 1.0);
+        // c2 = k * l_ref/(l^3)
+        c2_(i_elem) = k_(i_elem) * l_ref_(i_elem) / (l_(i_elem) * l_(i_elem) * l_(i_elem));
+    }
+};
+
+}  // namespace openturbine

src/system/springs/calculate_force_vectors.hpp

@@ -0,0 +1,27 @@
+#pragma once
+
+#include <KokkosBlas.hpp>
+#include <Kokkos_Core.hpp>
+
+#include "src/math/vector_operations.hpp"
+
+namespace openturbine {
+
+struct CalculateForceVectors {
+    Kokkos::View<double* [3]>::const_type r_;
+    Kokkos::View<double*>::const_type c1_;
+    Kokkos::View<double* [3]> f_;
+
+    KOKKOS_FUNCTION
+    void operator()(int i_elem) const {
+        auto r = Kokkos::subview(r_, i_elem, Kokkos::ALL);
+        auto f = Kokkos::subview(f_, i_elem, Kokkos::ALL);
+
+        // Calculate force vector components: f = c1 * r
+        for (int i = 0; i < 3; ++i) {
+            f(i) = c1_(i_elem) * r(i);
+        }
+    }
+};
+
+}  // namespace openturbine

src/system/springs/calculate_length.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <KokkosBlas.hpp>
+#include <Kokkos_Core.hpp>
+
+#include "src/math/vector_operations.hpp"
+
+namespace openturbine {
+
+struct CalculateLength {
+    Kokkos::View<double* [3]>::const_type r_;
+    Kokkos::View<double*> l_;
+
+    KOKKOS_FUNCTION
+    void operator()(int i_elem) const {
+        auto r = Kokkos::subview(r_, i_elem, Kokkos::ALL);
+        l_(i_elem) = sqrt(r(0) * r(0) + r(1) * r(1) + r(2) * r(2));
+    }
+};
+
+}  // namespace openturbine

src/system/springs/calculate_stiffness_matrix.hpp

@@ -0,0 +1,45 @@
+#pragma once
+
+#include <iostream>
+
+#include <KokkosBatched_Gemm_Decl.hpp>
+#include <KokkosBlas.hpp>
+#include <Kokkos_Core.hpp>
+
+#include "src/math/vector_operations.hpp"
+
+namespace openturbine {
+
+struct CalculateStiffnessMatrix {
+    using NoTranspose = KokkosBatched::Trans::NoTranspose;
+    using GemmDefault = KokkosBatched::Algo::Gemm::Default;
+    using Gemm = KokkosBatched::SerialGemm<NoTranspose, NoTranspose, GemmDefault>;
+    Kokkos::View<double*>::const_type c1_;
+    Kokkos::View<double*>::const_type c2_;
+    Kokkos::View<double* [3]>::const_type r_;
+    Kokkos::View<double*>::const_type l_;
+    Kokkos::View<double* [3][3]> r_tilde_;
+    Kokkos::View<double* [3][3]> a_;
+
+    KOKKOS_FUNCTION
+    void operator()(int i_elem) const {
+        auto r = Kokkos::subview(r_, i_elem, Kokkos::ALL);
+        auto r_tilde = Kokkos::subview(r_tilde_, i_elem, Kokkos::ALL, Kokkos::ALL);
+        auto a = Kokkos::subview(a_, i_elem, Kokkos::ALL, Kokkos::ALL);
+
+        // diagonal terms: c1 - c2 * l^2
+        const double diag_term = c1_(i_elem) - c2_(i_elem) * l_(i_elem) * l_(i_elem);
+        a(0, 0) = diag_term;
+        a(1, 1) = diag_term;
+        a(2, 2) = diag_term;
+
+        // non-diagonal terms: -c2 * r_tilde * r_tilde
+        auto temp = Kokkos::Array<double, 9>{0.};
+        auto temp_view = View_3x3(temp.data());
+        VecTilde(r, r_tilde);
+        Gemm::invoke(1., r_tilde, r_tilde, 0., temp_view);
+        KokkosBlas::serial_axpy(-c2_(i_elem), temp_view, a);
+    };
+};
+
+}  // namespace openturbine

tests/regression_tests/regression/CMakeLists.txt

@@ -10,5 +10,6 @@ target_sources(
     test_rotating_beam.cpp
     test_rotor.cpp
     test_solver.cpp
+    test_springs.cpp
     test_utilities.cpp
 )

tests/regression_tests/regression/test_springs.cpp

@@ -0,0 +1,60 @@
+#include <gtest/gtest.h>
+
+#include "test_utilities.hpp"
+
+#include "src/elements/springs/create_springs.hpp"
+#include "src/elements/springs/springs.hpp"
+#include "src/model/model.hpp"
+
+namespace openturbine::tests {
+
+inline auto SetUpSprings() {
+    auto model = Model();
+    // Add two nodes for the spring element
+    model.AddNode(
+        {0, 0, 0, 1, 0, 0, 0},  // First node at origin
+        {0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}
+    );
+    model.AddNode(
+        {1, 0, 0, 1, 0, 0, 0},  // Second node at (1,0,0)
+        {0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}
+    );
+
+    const auto springs_input = SpringsInput({SpringElement(
+        std::array{model.GetNode(0), model.GetNode(1)},
+        1000.,  // Spring stiffness
+        1.      // Undeformed length
+    )});
+
+    return CreateSprings(springs_input);
+}
+
+TEST(SpringsTest, NodeStateIndices) {
+    const auto springs = SetUpSprings();
+    auto node_state_indices_host = Kokkos::create_mirror(springs.node_state_indices);
+    Kokkos::deep_copy(node_state_indices_host, springs.node_state_indices);
+    EXPECT_EQ(node_state_indices_host(0, 0), 0);
+    EXPECT_EQ(node_state_indices_host(0, 1), 1);
+}
+
+TEST(SpringsTest, InitialPositionVector) {
+    const auto springs = SetUpSprings();
+    expect_kokkos_view_2D_equal(
+        springs.x0,
+        {
+            {1., 0., 0.},  // Vector from node 0 to node 1
+        }
+    );
+}
+
+TEST(SpringsTest, ReferenceLength) {
+    const auto springs = SetUpSprings();
+    expect_kokkos_view_1D_equal(springs.l_ref, {1.});  // Undeformed length
+}
+
+TEST(SpringsTest, SpringStiffness) {
+    const auto springs = SetUpSprings();
+    expect_kokkos_view_1D_equal(springs.k, {1000.});  // Spring stiffness
+}
+
+}  // namespace openturbine::tests

tests/unit_tests/elements/CMakeLists.txt

@@ -1,5 +1,6 @@
 add_subdirectory(beams)
 add_subdirectory(masses)
+add_subdirectory(springs)
 
 # Specify the source files for the unit test executable
 target_sources(