Commit 7b2b01a4 authored by Stefan Pfeifer's avatar Stefan Pfeifer

Add tests directory and change project structure

parent 833a0378
include(shared.pri)
QT += core gui
QT += widgets
TARGET = bow-design-simulator
TEMPLATE = app
SOURCES += src/main.cpp
#-------------------------------------------------
#
# Project created by QtCreator 2016-06-06T06:15:39
#
#-------------------------------------------------
QT += core gui
QT += widgets
TARGET = bow-design-simulator
TEMPLATE = app
SOURCES += src/main.cpp \
src/gui/MainWindow.cpp \
src/fem/System.cpp \
src/fem/elements/Element.cpp
HEADERS += \
src/gui/MainWindow.hpp \
src/fem/System.hpp \
src/fem/elements/Element.hpp \
src/fem/elements/BarElement.hpp \
src/fem/Node.hpp \
src/fem/elements/MassElement.hpp \
src/fem/View.hpp
INCLUDEPATH += /home/s/Libraries/Eigen-3.2.8
INCLUDEPATH += /home/s/Libraries/Eigen-3.2.8 \
/home/s/Libraries/Catch-1.5.6 \
SOURCES += src/gui/MainWindow.cpp \
src/fem/System.cpp \
src/fem/elements/Element.cpp
HEADERS += src/gui/MainWindow.hpp \
src/fem/System.hpp \
src/fem/elements/Element.hpp \
src/fem/elements/BarElement.hpp \
src/fem/Node.hpp \
src/fem/elements/MassElement.hpp \
src/fem/View.hpp
...@@ -3,6 +3,8 @@ ...@@ -3,6 +3,8 @@
#include "Node.hpp" #include "Node.hpp"
#include "elements/Element.hpp" #include "elements/Element.hpp"
#include <Eigen/Eigenvalues>
#include <vector> #include <vector>
#include <iostream> #include <iostream>
...@@ -59,6 +61,7 @@ public: ...@@ -59,6 +61,7 @@ public:
void add_element(Element& element) void add_element(Element& element)
{ {
element.set_state(get_u(), get_v());
elements.push_back(&element); elements.push_back(&element);
} }
...@@ -132,19 +135,28 @@ public: ...@@ -132,19 +135,28 @@ public:
}); });
} }
void solve_dynamics(double dt, const std::function<bool()>& callback) void solve_dynamics(double timestep_factor, const std::function<bool()>& callback)
{ {
VectorXd M(dofs()); VectorXd M(dofs());
VectorXd q(dofs()); VectorXd q(dofs());
MatrixXd K(dofs(), dofs()); MatrixXd K(dofs(), dofs());
update_element_states();
get_mass_matrix(M); get_mass_matrix(M);
get_internal_forces(q); get_internal_forces(q);
get_tangent_stiffness(K);
// Timestep estimation
Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> eigen_solver(K, M.asDiagonal(), Eigen::DecompositionOptions::EigenvaluesOnly);
if(eigen_solver.info() != Eigen::Success)
{
throw std::runtime_error("Failed to compute eigenvalues of the system");
}
double omega_max = std::sqrt(eigen_solver.eigenvalues().maxCoeff());
double dt = timestep_factor*2.0/omega_max;
// Version with a more accurate velocity based on central differences // Version with a more accurate velocity based on central differences
a = M.asDiagonal().inverse()*(p - q); a = M.asDiagonal().inverse()*(p - q);
// Previous displacements // Previous displacements
Eigen::VectorXd u_p2(dofs()); Eigen::VectorXd u_p2(dofs());
Eigen::VectorXd u_p1 = u - dt*v + dt*dt/2.0*a; Eigen::VectorXd u_p1 = u - dt*v + dt*dt/2.0*a;
...@@ -232,5 +244,29 @@ public: // Todo: private ...@@ -232,5 +244,29 @@ public: // Todo: private
{ {
e->get_tangent_stiffness(view); e->get_tangent_stiffness(view);
} }
std::cout << "K = " << K << "\n";
}
// Numeric stiffness matrix via central difference quotient. Only for testing purposes.
void get_numeric_tangent_stiffness(MatrixXd& K, double h)
{
K.setZero();
Eigen::VectorXd q_fwd = Eigen::VectorXd::Zero(dofs());
Eigen::VectorXd q_bwd = Eigen::VectorXd::Zero(dofs());
for(std::size_t i = 0; i < dofs(); ++i)
{
double ui = u(i);
u(i) = ui + h;
get_internal_forces(q_fwd);
u(i) = ui - h;
get_internal_forces(q_bwd);
u(i) = ui;
K.col(i) = (q_fwd - q_bwd)/(2.0*h);
}
} }
}; };
...@@ -4,12 +4,7 @@ double Element::get_kinetic_energy(const VectorView<Dof> v) const ...@@ -4,12 +4,7 @@ double Element::get_kinetic_energy(const VectorView<Dof> v) const
{ {
double T = 0.0; double T = 0.0;
get_masses(VectorView<Dof>([&](Dof dof) get_masses(VectorView<Dof>(nullptr, [&](Dof dof, double val)
{
return 0.0; // Todo: Unreachable, replace with assert of some sort?
},
[&](Dof dof, double val)
{ {
T += 0.5*val*std::pow(v(dof), 2); T += 0.5*val*std::pow(v(dof), 2);
})); }));
......
...@@ -2,6 +2,8 @@ ...@@ -2,6 +2,8 @@
#include "fem/elements/BarElement.hpp" #include "fem/elements/BarElement.hpp"
#include "fem/elements/MassElement.hpp" #include "fem/elements/MassElement.hpp"
#include <catch.hpp>
#include <iostream> #include <iostream>
void test_harmonic_oscillator() void test_harmonic_oscillator()
...@@ -32,7 +34,7 @@ void test_harmonic_oscillator() ...@@ -32,7 +34,7 @@ void test_harmonic_oscillator()
//assert!(delta < omega0); // Make sure the system is underdamped // Todo: Reintroduce this check //assert!(delta < omega0); // Make sure the system is underdamped // Todo: Reintroduce this check
system.solve_dynamics(1e-5, [&]() system.solve_dynamics(0.01, [&]()
{ {
// Time // Time
double t = system.get_time(); double t = system.get_time();
...@@ -51,8 +53,6 @@ void test_harmonic_oscillator() ...@@ -51,8 +53,6 @@ void test_harmonic_oscillator()
double error_s = std::abs(s_num - s_ref); double error_s = std::abs(s_num - s_ref);
double error_v = std::abs(v_num - v_ref); double error_v = std::abs(v_num - v_ref);
//std::cout << error_s << "\n";
return t < T; return t < T;
}); });
} }
......
include(shared.pri)
TEMPLATE = app
CONFIG += console c++11
CONFIG -= app_bundle
CONFIG -= qt
SOURCES += tests/main.cpp \
tests/harmonic_oscillator.cpp
#include "../src/fem/System.hpp"
#include "../src/fem/elements/BarElement.hpp"
#include "../src/fem/elements/MassElement.hpp"
#include <catch.hpp>
#include <iostream>
TEST_CASE("Dynamic solution harmonic oscillator")
{
// https://de.wikipedia.org/wiki/Schwingung#Linear_ged.C3.A4mpfte_Schwingung
double l = 1.0;
double k = 100.0;
double d = 0.0*10.0;
double m = 5.0;
double s0 = 0.1; // Initial displacement
System system;
Node node_a = system.create_node({{ 0.0, 0.0, 0.0}}, {{Dof::Type::Fixed, Dof::Type::Fixed, Dof::Type::Fixed}});
Node node_b = system.create_node({{l + s0, 0.0, 0.0}}, {{Dof::Type::Active, Dof::Type::Fixed, Dof::Type::Fixed}});
BarElement element1(node_a, node_b, l, l*k, l*d, 0.0);
MassElement element2(node_b, m, 0.0);
system.add_element(element1);
system.add_element(element2);
// Constants for the analytical solution
double delta = d/(2.0*m); // Decay constant
double omega0 = std::sqrt(k/m); // Natural frequency of the undamped system
double omega = std::sqrt(omega0*omega0 - delta*delta); // Natural frequency
double T = 2.0*M_PI/omega; // Period length
double A = delta/omega*s0;
double B = s0;
//assert!(delta < omega0); // Make sure the system is underdamped // Todo: Reintroduce this check
system.solve_dynamics(0.01, [&]()
{
// Time
double t = system.get_time();
// Numerical solution
double s_num = system.get_u()(node_b.x) - l;
double v_num = system.get_v()(node_b.x);
// Analytical solution
double s_ref = std::exp(-delta*t)*(A*std::sin(omega*t) + B*std::cos(omega*t));
double v_ref = -std::exp(-delta*t)*((omega*B + delta*A)*std::sin(omega*t) + (delta*B - omega*A)*std::cos(omega*t));
// Error
double error_s = std::abs(s_num - s_ref);
double error_v = std::abs(v_num - v_ref);
REQUIRE(error_s < 1e-4);
REQUIRE(error_v < 1e-4);
//std::cout << t << ", " << s_ref << "\n";
//std::cout << error_s << "\n";
return t < T;
});
}
#define CATCH_CONFIG_MAIN
#include <catch.hpp>
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