Added a regression test with a chain of masses/springs which are in s…

…tationary equilibrium.

This test is a more complex configuration with a very simple analytic solution.  More importantly,
it is parametrically scalable, so it can be easily scaled up/down for performance testing.

tests/regression_tests/regression/test_spring_mass_system.cpp

@@ -124,4 +124,116 @@ TEST(SpringMassSystemTest, FinalDisplacement) {
     SetUpSpringMassSystem();
 }
 
+inline auto SetUpSpringMassChainSystem() {
+    auto model = Model();
+
+    // Add two nodes for the spring element
+    constexpr auto number_of_masses = 10U;
+    constexpr auto displacement = .5;
+    auto position = 0.;
+    model.AddNode(
+        {position, 0., 0., 1., 0., 0., 0.},  // First node at origin -- initial position
+        {0., 0., 0., 1., 0., 0., 0.},  // initial displacement
+        {0., 0., 0., 0., 0., 0.},      // initial velocity
+        {0., 0., 0., 0., 0., 0.}       // initial acceleration
+    );
+    for(auto mass_number = 0U; mass_number < number_of_masses; ++mass_number) {
+        position += displacement;
+        model.AddNode(
+            {position, 0., 0., 1., 0., 0., 0.},  // Second node at (2,0,0) -- initial position
+            {0., 0., 0., 1., 0., 0., 0.},  // initial displacement
+            {0., 0., 0., 0., 0., 0.},      // initial velocity
+            {0., 0., 0., 0., 0., 0.}       // initial acceleration
+        );
+    }
+    position += displacement;
+    model.AddNode(
+        {position, 0., 0., 1., 0., 0., 0.},  // Second node at (2,0,0) -- initial position
+        {0., 0., 0., 1., 0., 0., 0.},  // initial displacement
+        {0., 0., 0., 0., 0., 0.},      // initial velocity
+        {0., 0., 0., 0., 0., 0.}       // initial acceleration
+    );
+
+    // No beams
+    const auto beams_input = BeamsInput({}, {0., 0., 0.});
+    auto beams = CreateBeams(beams_input);
+
+    // We need to add a mass element with identity for mass matrix to create a spring-mass system
+    constexpr auto m = 1.;
+    constexpr auto j = 1.;
+    constexpr auto mass_matrix = std::array{
+        std::array{m, 0., 0., 0., 0., 0.},  // mass in x-direction
+        std::array{0., m, 0., 0., 0., 0.},  // mass in y-direction
+        std::array{0., 0., m, 0., 0., 0.},  // mass in z-direction
+        std::array{0., 0., 0., j, 0., 0.},  // inertia around x-axis
+        std::array{0., 0., 0., 0., j, 0.},  // inertia around y-axis
+        std::array{0., 0., 0., 0., 0., j},  // inertia around z-axis
+    };
+
+    auto mass_elements = std::vector<MassElement>{};
+    for(auto mass_number = 0U; mass_number < number_of_masses; ++mass_number) {
+        mass_elements.emplace_back(MassElement(model.GetNode(mass_number+1), mass_matrix));
+    }
+    const auto masses_input = MassesInput(mass_elements, {0., 0., 0.});
+    auto masses = CreateMasses(masses_input);
+
+    // Create springs
+    const auto k = 10.;
+    auto spring_elements = std::vector<SpringElement>{};
+    for(auto mass_number = 0U; mass_number <= number_of_masses; ++mass_number) {
+        spring_elements.emplace_back(SpringElement{std::array{model.GetNode(mass_number), model.GetNode(mass_number+1)}, k, 0.});
+    }
+    const auto springs_input = SpringsInput(spring_elements);
+    auto springs = CreateSprings(springs_input);
+
+    // Create elements
+    auto elements = Elements{beams, masses, springs};
+
+    // Add fixed BC to the first node
+    model.AddFixedBC(model.GetNode(0));
+    model.AddFixedBC(model.GetNode(number_of_masses+1));
+
+    // Set up solver components
+    auto state = model.CreateState();
+    auto constraints = Constraints(model.GetConstraints());
+    assemble_node_freedom_allocation_table(state, elements, constraints);
+    compute_node_freedom_map_table(state);
+    create_element_freedom_table(elements, state);
+    create_constraint_freedom_table(constraints, state);
+    auto solver = Solver(
+        state.ID, state.node_freedom_allocation_table, state.node_freedom_map_table,
+        elements.NumberOfNodesPerElement(), elements.NodeStateIndices(), constraints.num_dofs,
+        constraints.type, constraints.base_node_freedom_table, constraints.target_node_freedom_table,
+        constraints.row_range
+    );
+
+    // The spring-mass system should move periodically between -2 and 2 (position of the second node)
+    // with simple harmonic motion where the time period is 2 * pi * sqrt(m / k)
+    const double T = 2. * M_PI * sqrt(m / k);
+    constexpr auto num_steps = 1000;
+
+    // Set up step parameters
+    constexpr bool is_dynamic_solve(true);
+    constexpr size_t max_iter(6);
+    constexpr double rho_inf(0.);                                // No damping
+    const double step_size(T / static_cast<double>(num_steps));  // Calculate step size
+    auto parameters = StepParameters(is_dynamic_solve, max_iter, step_size, rho_inf);
+
+    auto q = Kokkos::create_mirror(state.q);
+
+    // Run simulation for T seconds
+    for (auto time_step = 1U; time_step <= num_steps; ++time_step) {
+        auto converged = Step(parameters, solver, elements, state, constraints);
+        EXPECT_TRUE(converged);
+        Kokkos::deep_copy(q, state.q);
+        for(auto node = 0U; node < q.extent(0); ++node) {
+            EXPECT_EQ(q(node, 0), 0.);
+        }
+    }
+}
+
+TEST(SpringMassChainSystemTest, FinalDisplacement) {
+    SetUpSpringMassChainSystem();
+}
+
 }  // namespace openturbine::tests