Commit 8abbe339 authored by Stefan Pfeifer's avatar Stefan Pfeifer

Add stiffness matrix test for bar element

parent 7b2b01a4
......@@ -185,7 +185,7 @@ public: // Todo: private
}
}
void get_mass_matrix(VectorXd& M)
void get_mass_matrix(VectorXd& M) const
{
M.setZero();
......@@ -207,7 +207,7 @@ public: // Todo: private
}
}
void get_internal_forces(VectorXd& q)
void get_internal_forces(VectorXd& q) const
{
q.setZero();
......@@ -229,7 +229,7 @@ public: // Todo: private
}
}
void get_tangent_stiffness(MatrixXd& K)
void get_tangent_stiffness(MatrixXd& K) const
{
K.setZero();
......@@ -244,8 +244,6 @@ public: // Todo: private
{
e->get_tangent_stiffness(view);
}
std::cout << "K = " << K << "\n";
}
// Numeric stiffness matrix via central difference quotient. Only for testing purposes.
......@@ -261,9 +259,13 @@ public: // Todo: private
double ui = u(i);
u(i) = ui + h;
update_element_states();
get_internal_forces(q_fwd);
u(i) = ui - h;
update_element_states();
get_internal_forces(q_bwd);
u(i) = ui;
K.col(i) = (q_fwd - q_bwd)/(2.0*h);
......
......@@ -36,8 +36,8 @@ public:
dx = u(nodes[1].x) - u(nodes[0].x);
dy = u(nodes[1].y) - u(nodes[0].y);
L_new = std::hypot(dx, dy); // Actual length
L_dot = 1.0/L*(dx*(v(nodes[1].x) - v(nodes[0].x)) + dy*(v(nodes[1].y) - v(nodes[0].y)));
L_new = std::hypot(dx, dy);
L_dot = 1.0/L_new*(dx*(v(nodes[1].x) - v(nodes[0].x)) + dy*(v(nodes[1].y) - v(nodes[0].y)));
}
virtual void get_masses(VectorView<Dof> M) const override
......@@ -54,16 +54,16 @@ public:
{
double N = get_normal_force();
q(nodes[0].x) -= N*dx/L;
q(nodes[0].y) -= N*dy/L;
q(nodes[1].x) += N*dx/L;
q(nodes[1].y) += N*dy/L;
q(nodes[0].x) -= N*dx/L_new;
q(nodes[0].y) -= N*dy/L_new;
q(nodes[1].x) += N*dx/L_new;
q(nodes[1].y) += N*dy/L_new;
}
virtual void get_tangent_stiffness(MatrixView<Dof> K) const override
{
double c0 = EA*(L_new - L)/(L_new*L);
double c1 = EA/(L*L*L);
double c1 = EA/std::pow(L_new, 3);
// col 0
K(nodes[0].x, nodes[0].x) += c1*dx*dx + c0;
......
......@@ -6,5 +6,6 @@ CONFIG -= app_bundle
CONFIG -= qt
SOURCES += tests/main.cpp \
tests/harmonic_oscillator.cpp
tests/harmonic_oscillator.cpp \
tests/tangent_stiffness.cpp
#include "../src/fem/System.hpp"
#include "../src/fem/elements/BarElement.hpp"
#include <catch.hpp>
#include <iostream>
void check_stiffness_matrix(System& system)
{
Eigen::MatrixXd K_sym = Eigen::MatrixXd::Zero(system.dofs(), system.dofs());
Eigen::MatrixXd K_num = Eigen::MatrixXd::Zero(system.dofs(), system.dofs());
system.get_tangent_stiffness(K_sym);
system.get_numeric_tangent_stiffness(K_num, 1e-8);
REQUIRE(K_sym.isApprox(K_num, 1e-5)); // Equal to numeric derivatives
REQUIRE(K_sym.isApprox(K_num.transpose(), 1e-5)); // Symmetric
}
TEST_CASE("Tangent stiffness matrix: BarElement")
{
auto check_at = [](double dx0, double dy0, double dx1, double dy1)
{
double L = 1.0;
double EA = 100.0;
System system;
Node node0 = system.create_node({{ dx0, dy0, 0.0}}, {{Dof::Type::Active, Dof::Type::Active, Dof::Type::Fixed}});
Node node1 = system.create_node({{L + dx1, dy1, 0.0}}, {{Dof::Type::Active, Dof::Type::Active, Dof::Type::Fixed}});
BarElement element01(node0, node1, L, EA, 0.0, 0.0);
system.add_element(element01);
check_stiffness_matrix(system);
};
check_at(0.0, 0.0, 0.0, 0.0);
check_at(0.0, 0.0, 1.5, 1.5);
check_at(0.0, 0.0, -1.5, 1.5);
check_at(0.0, 0.0, 1.5, -1.5);
check_at(0.0, 0.0, -1.5, -1.5);
check_at( 1.5, 1.5, 0.0, 0.0);
check_at(-1.5, 1.5, 0.0, 0.0);
check_at( 1.5, -1.5, 0.0, 0.0);
check_at(-1.5, -1.5, 0.0, 0.0);
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment