Commit 045ee37b authored by Stefan Pfeifer's avatar Stefan Pfeifer

Merge remote-tracking branch 'origin/fix-tests' into develop

parents 9aa1882b e15ca78e
......@@ -6,10 +6,10 @@ set(CMAKE_CXX_STANDARD 14)
# Parameters used for configuring files
set(APPLICATION_NAME "bow-simulator")
set(APPLICATION_TITLE "Bow Simulator")
set(APPLICATION_VERSION "0.5")
set(APPLICATION_VERSION "0.6")
set(APPLICATION_WEBSITE "bow-simulator.org")
set(APPLICATION_MAINTAINER "Stefan Pfeifer")
set(APPLICATION_COPYRIGHT "Copyright (C) 2016-2018 Stefan Pfeifer")
set(APPLICATION_COPYRIGHT "Copyright (C) 2016-2019 Stefan Pfeifer")
set(APPLICATION_LICENSE "GNU General Public License v3.0")
set(APPLICATION_DESCRIPTION_SHORT "Bow and arrow physics simulation")
set(APPLICATION_DESCRIPTION_LONG "Software tool for simulating the static and dynamic performance of bow designs")
......@@ -81,18 +81,13 @@ include_directories(${SOURCE_DIR})
# System thread library (https://stackoverflow.com/a/39547577)
find_package(Threads REQUIRED)
# Main executable
# Build application library
configure_file(source/config.hpp.in ${CMAKE_BINARY_DIR}/config.hpp)
include_directories(${PROJECT_BINARY_DIR})
include_directories(source)
add_executable(
bow-simulator
source/main.cpp
${CMAKE_BINARY_DIR}/config.hpp
configure_file(source/config.hpp ${CMAKE_BINARY_DIR}/config.hpp)
include_directories(source ${PROJECT_BINARY_DIR})
add_library(
main-library
source/fem/elements/BarElement.cpp
source/fem/elements/BeamElement.cpp
source/fem/elements/MassElement.cpp
......@@ -104,16 +99,13 @@ add_executable(
source/fem/StaticSolver.cpp
source/fem/DynamicSolver.cpp
source/fem/System.cpp
source/bow/LimbProperties.cpp
source/bow/LayerProperties.cpp
source/bow/input/InputData.cpp
source/bow/BowModel.cpp
source/numerics/ArcCurve.cpp
source/numerics/CubicSpline.cpp
source/numerics/Series.cpp
source/gui/Application.cpp
source/gui/Application.hpp
source/gui/MainWindow.cpp
......@@ -128,7 +120,6 @@ add_executable(
source/gui/ProgressDialog.hpp
source/gui/EditableTabBar.cpp
source/gui/EditableTabBar.hpp
source/gui/input/dialogs/CommentDialog.cpp
source/gui/input/dialogs/CommentDialog.hpp
source/gui/input/dialogs/GroupDialog.cpp
......@@ -147,7 +138,6 @@ add_executable(
source/gui/input/dialogs/StringDialog.hpp
source/gui/input/dialogs/WidthDialog.cpp
source/gui/input/dialogs/WidthDialog.hpp
source/gui/input/editors/BowEditor.cpp
source/gui/input/editors/BowEditor.hpp
source/gui/input/editors/DoubleEditor.cpp
......@@ -161,7 +151,6 @@ add_executable(
source/gui/input/editors/TreeEditor.cpp
source/gui/input/editors/TreeEditor.hpp
source/gui/input/editors/TreeItem.hpp
source/gui/input/views/LimbView.cpp
source/gui/input/views/LimbView.hpp
source/gui/input/views/LimbMesh.cpp
......@@ -172,7 +161,6 @@ add_executable(
source/gui/input/views/ProfileView.hpp
source/gui/input/views/SplineView.cpp
source/gui/input/views/SplineView.hpp
source/gui/output/ComboPlot.cpp
source/gui/output/ComboPlot.hpp
source/gui/output/EnergyPlot.cpp
......@@ -187,40 +175,31 @@ add_executable(
source/gui/output/Slider.hpp
source/gui/output/StressPlot.cpp
source/gui/output/StressPlot.hpp
source/external/qcustomplot/qcustomplot.cpp
resources/resources.qrc
)
add_dependencies(main-library eigen catch json boost)
target_link_libraries(main-library Qt5::Widgets Qt5::PrintSupport Threads::Threads)
# Declare dependencies, link libraries
# Build application executable
add_dependencies(bow-simulator eigen catch json boost)
target_link_libraries(bow-simulator Qt5::Widgets Qt5::PrintSupport Threads::Threads)
add_executable(bow-simulator source/main.cpp resources/resources.qrc)
target_link_libraries(bow-simulator main-library)
# Test executable
# Build test executable
#add_executable(
# bow-simulator-tests
# source/tests/main.cpp
# source/tests/tangent_stiffness.cpp
# source/tests/large_deformation_beams.cpp
# source/tests/harmonic_oscillator.cpp
# source/tests/bar_trusses.cpp
add_executable(
bow-simulator-test
tests/main.cpp
tests/fem/BarTrusses.cpp
tests/fem/HarmonicOscillator.cpp
tests/fem/LargeDeformationBeams.cpp
tests/fem/TangentStiffness.cpp
)
# source/fem/elements/BarElement.cpp
# source/fem/elements/BeamElement.cpp
# source/fem/elements/MassElement.cpp
# source/fem/elements/ContactElement.cpp
# source/fem/elements/ContactSurface.cpp
# source/fem/elements/ConstraintElement.cpp
# source/fem/Node.cpp
# source/fem/DynamicSolver.cpp
# source/fem/System.cpp
#)
target_link_libraries(bow-simulator-test main-library)
# Release packages for different platforms
# Build release packages for different platforms
if(UNIX)
add_subdirectory(platforms/linux)
......
......@@ -245,7 +245,7 @@ void BowModel::init_string(const Callback& callback)
// Find a element length at which the brace height difference is zero
// Todo: Perhaps limit the step size of the root finding algorithm to increase robustness.
double l = (points[1] - points[0]).norm();
if(!try_element_length(l) > 0)
if(try_element_length(l) <= 0.0)
throw std::runtime_error("Invalid input: Brace height is too low");
double dl = 1e-3*l; // Initial step length, later adjusted by the algorithm // Magic number
......
#include "DynamicSolver.hpp"
DynamicSolver::DynamicSolver(System& system, double dt, double f, const StopFn& stop)
DynamicSolver::DynamicSolver(System& system, double dt, double f_sample, const StopFn& stop)
: system(system),
stop(stop),
dt(dt),
n(std::max(std::ceil(1.0/(f*dt)), 1.0))
n(std::max(std::ceil(1.0/(f_sample*dt)), 1.0))
{
// Initialise previous displacement
u_p2 = system.get_u() - dt*system.get_v() + dt*dt/2.0*system.get_a();
......
......@@ -8,7 +8,7 @@ class DynamicSolver
public:
using StopFn = std::function<bool()>;
DynamicSolver(System& system, double dt, double f, const StopFn& stop);
DynamicSolver(System& system, double dt, double f_sample, const StopFn& stop);
static double estimate_timestep(const System& system, double factor);
bool step();
......
......@@ -4,6 +4,7 @@
#include "numerics/Optimization.hpp"
#include <Eigen/Core>
// Todo: Make constraint function a member with templated type instead of using inheritance.
class StaticSolver
{
public:
......
#include "System.hpp"
#include <QtCore>
System::System()
: t(0.0), n_a(0), n_f(0)
{
......
......@@ -5,8 +5,6 @@
#include "fem/System.hpp"
#include <map>
#include <QtCore>
// Holds a collection of segments (two nodes, two distances) and points (one node) and creates/removes
// contact elements for them as needed. It uses the Sweep and Prune broadphase algorithm [1],[2] along
// both axes for determining which points and segments can actually come into contact.
......
#include "fem/System.hpp"
#include "fem/StaticSolver.hpp"
#include "fem/elements/BarElement.hpp"
#include <catch.hpp>
#include <iostream>
TEST_CASE("small-deformation-bar-truss")
{
double L = 0.5;
double EA = 21000.0; // Steel rod, 1cm x 1cm
double F = 10.0;
System system;
Node node_01 = system.create_node({DofType::Fixed, DofType::Fixed, DofType::Fixed}, { 0.0, 0.0, 0.0});
Node node_02 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, { L, 0.0, 0.0});
Node node_03 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, {2.0*L, 0.0, 0.0});
Node node_04 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, {3.0*L, 0.0, 0.0});
Node node_05 = system.create_node({DofType::Active, DofType::Fixed, DofType::Fixed}, {4.0*L, 0.0, 0.0});
Node node_06 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, { 0.0, L, 0.0});
Node node_07 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, { L, L, 0.0});
Node node_08 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, {2.0*L, L, 0.0});
Node node_09 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, {3.0*L, L, 0.0});
Node node_10 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, {4.0*L, L, 0.0});
system.add_element(BarElement(node_01, node_02, L, EA, 0.0));
system.add_element(BarElement(node_02, node_03, L, EA, 0.0));
system.add_element(BarElement(node_03, node_04, L, EA, 0.0));
system.add_element(BarElement(node_04, node_05, L, EA, 0.0));
system.add_element(BarElement(node_06, node_07, L, EA, 0.0));
system.add_element(BarElement(node_07, node_08, L, EA, 0.0));
system.add_element(BarElement(node_08, node_09, L, EA, 0.0));
system.add_element(BarElement(node_09, node_10, L, EA, 0.0));
system.add_element(BarElement(node_01, node_06, L, EA, 0.0));
system.add_element(BarElement(node_02, node_07, L, EA, 0.0));
system.add_element(BarElement(node_03, node_08, L, EA, 0.0));
system.add_element(BarElement(node_04, node_09, L, EA, 0.0));
system.add_element(BarElement(node_05, node_10, L, EA, 0.0));
system.add_element(BarElement(node_01, node_07, std::sqrt(2)*L, EA, 0.0));
system.add_element(BarElement(node_07, node_03, std::sqrt(2)*L, EA, 0.0));
system.add_element(BarElement(node_03, node_09, std::sqrt(2)*L, EA, 0.0));
system.add_element(BarElement(node_09, node_05, std::sqrt(2)*L, EA, 0.0));
node_02[1].p_mut() = -F;
node_04[1].p_mut() = -F;
StaticSolverLC solver(system);
solver.solve();
double s_numeric = -node_03[1].u();
double s_analytic = (4.0+2.0*std::sqrt(2))*F*L/EA;
REQUIRE(std::abs(s_numeric - s_analytic) < 1e-6);
}
// Todo: Why does the displacement control not allow passing the point 0.5*H?
// Read section on displacement control in 'Nonlinear Finite Element Analysis of Solids and Structures (René De Borst,Mike A. Crisfield,Joris J. C.)
TEST_CASE("large-deformation-bar-truss")
{
double H = 1.0;
double EA = 10000.0;
double F = 1500.0;
System system;
Node node01 = system.create_node({DofType::Fixed, DofType::Fixed, DofType::Fixed}, { 0.0, 0.0, 0.0});
Node node02 = system.create_node({DofType::Fixed, DofType::Active, DofType::Fixed}, { H, H, 0.0});
Node node03 = system.create_node({DofType::Fixed, DofType::Fixed, DofType::Fixed}, {2.0*H, 0.0, 0.0});
system.add_element(BarElement(node01, node02, M_SQRT2*H, EA, 0.0));
system.add_element(BarElement(node02, node03, M_SQRT2*H, EA, 0.0));
StaticSolverDC solver(system, node02[1]);//, 0.6*H, 100);
/*
while(solver.step())
{
double s = node02[1].u();
double f_num = node02[1].p();
double L0 = M_SQRT2*H;
double L = hypot(H, s);
double f_ref = 2.0*EA*s*(L - L0)/(L*L0);
double error = std::abs(f_num - f_ref); // Absolute error because f_ref = 0 initially
REQUIRE(error < 1e-9);
}
*/
}
#include "fem/System.hpp"
#include "fem/StaticSolver.hpp"
#include "fem/elements/BarElement.hpp"
#include <catch.hpp>
#include <iostream>
TEST_CASE("small-deformation-bar-truss")
{
double L = 0.5;
double EA = 21000.0; // Steel rod, 1cm x 1cm
double F = 10.0;
System system;
Node node_01 = system.create_node({false, false, false}, { 0.0, 0.0, 0.0});
Node node_02 = system.create_node({ true, true, false}, { L, 0.0, 0.0});
Node node_03 = system.create_node({ true, true, false}, {2.0*L, 0.0, 0.0});
Node node_04 = system.create_node({ true, true, false}, {3.0*L, 0.0, 0.0});
Node node_05 = system.create_node({ true, false, false}, {4.0*L, 0.0, 0.0});
Node node_06 = system.create_node({ true, true, false}, { 0.0, L, 0.0});
Node node_07 = system.create_node({ true, true, false}, { L, L, 0.0});
Node node_08 = system.create_node({ true, true, false}, {2.0*L, L, 0.0});
Node node_09 = system.create_node({ true, true, false}, {3.0*L, L, 0.0});
Node node_10 = system.create_node({ true, true, false}, {4.0*L, L, 0.0});
system.mut_elements().add(BarElement(system, node_01, node_02, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_02, node_03, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_03, node_04, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_04, node_05, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_06, node_07, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_07, node_08, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_08, node_09, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_09, node_10, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_01, node_06, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_02, node_07, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_03, node_08, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_04, node_09, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_05, node_10, L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_01, node_07, std::sqrt(2)*L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_07, node_03, std::sqrt(2)*L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_03, node_09, std::sqrt(2)*L, EA, 0.0));
system.mut_elements().add(BarElement(system, node_09, node_05, std::sqrt(2)*L, EA, 0.0));
system.set_p(node_02.y, -F);
system.set_p(node_04.y, -F);
StaticSolverLC solver(system);
solver.solve();
double s_num = -system.get_u(node_03.y);
double s_ref = (4.0 + 2.0*std::sqrt(2))*F*L/EA;
REQUIRE(std::abs(s_num - s_ref) < 1e-6);
}
// Todo: Why does the displacement control not allow passing the point 0.5*H?
// Read section on displacement control in 'Nonlinear Finite Element Analysis of Solids and Structures (René De Borst,Mike A. Crisfield,Joris J. C.)
TEST_CASE("large-deformation-bar-truss")
{
double H = 1.0;
double EA = 10000.0;
double F = 1500.0;
System system;
Node node01 = system.create_node({false, false, false}, { 0.0, 0.0, 0.0});
Node node02 = system.create_node({false, true, false}, { H, H, 0.0});
Node node03 = system.create_node({false, false, false}, {2.0*H, 0.0, 0.0});
system.mut_elements().add(BarElement(system, node01, node02, M_SQRT2*H, EA, 0.0));
system.mut_elements().add(BarElement(system, node02, node03, M_SQRT2*H, EA, 0.0));
StaticSolverDC solver(system, node02.y);
for(double d = H; d > 0.6*H; d -= 0.05*H)
{
// Numerical solution
solver.solve(d);
double s = system.get_u(node02.y);
double f_num = system.get_p(node02.y);
// Analytical solution
double L0 = M_SQRT2*H;
double L = hypot(H, s);
double f_ref = 2.0*EA*s*(L - L0)/(L*L0);
// Error
REQUIRE(std::abs(f_num - f_ref) < 1e-9);
}
}
......@@ -8,7 +8,7 @@
TEST_CASE("harmonic-oscillator")
{
// https://de.wikipedia.org/wiki/Schwingung#Linear_ged.C3.A4mpfte_Schwingung
// https://en.wikipedia.org/wiki/Harmonic_oscillator
double l = 1.0;
double k = 100.0;
double m = 5.0;
......@@ -16,38 +16,31 @@ TEST_CASE("harmonic-oscillator")
double s0 = 0.1; // Initial displacement
System system;
Node node_a = system.create_node({DofType::Fixed, DofType::Fixed, DofType::Fixed}, { 0.0, 0.0, 0.0});
Node node_b = system.create_node({DofType::Active, DofType::Fixed, DofType::Fixed}, {l + s0, 0.0, 0.0});
Node node_a = system.create_node({false, false, false}, { 0.0, 0.0, 0.0});
Node node_b = system.create_node({true, false, false}, {l + s0, 0.0, 0.0});
system.add_element(BarElement(node_a, node_b, l, l*k, 0.0));
system.add_element(MassElement(node_b, m, 0.0));
system.mut_elements().add(BarElement(system, node_a, node_b, l, l*k, 0.0));
system.mut_elements().add(MassElement(system, node_b, m, 0.0));
// Constants for the analytical solution
double omega = std::sqrt(k/m); // Natural frequency
double T = 2.0*M_PI/omega; // Period length
DynamicSolver solver(system, 0.01, 1e-15, [&]{ return system.get_t() >= T; });
DynamicSolver solver(system, 0.001, 1000.0, [&]{ return system.get_t() >= T; });
while(solver.step())
{
// Time
double t = system.get_t();
//std::cout << "t = " << t << "\n";
// Numerical solution
double s_num = node_b[0].u() - l;
double v_num = node_b[0].v();
double s_num = system.get_u(node_b.x) - l;
double v_num = system.get_v(node_b.x);
// Analytical solution
double t = system.get_t();
double s_ref = s0*cos(omega*t);
double v_ref = -s0*omega*sin(omega*t);
// Error
double error_s = std::abs(s_num - s_ref);
double error_v = std::abs(v_num - v_ref);
REQUIRE(error_s < 0.62e-5);
REQUIRE(error_v < 0.50e-4);
REQUIRE(std::abs(s_num - s_ref) < 0.62e-5);
REQUIRE(std::abs(v_num - v_ref) < 0.50e-4);
}
}
......@@ -30,40 +30,39 @@ TEST_CASE("large-deformation-cantilever")
// Create nodes
for(unsigned i = 0; i < N+1; ++i)
{
DofType type = (i == 0) ? DofType::Fixed : DofType::Active;
nodes.push_back(system.create_node({type, type, type}, {double(i)/double(N)*L, 0.0, 0.0}));
bool active = (i != 0);
nodes.push_back(system.create_node({active, active, active}, {double(i)/double(N)*L, 0.0, 0.0}));
}
// Create elements
for(unsigned i = 0; i < N; ++i)
{
BeamElement element(nodes[i], nodes[i+1], 0.0, L/double(N));
BeamElement element(system, nodes[i], nodes[i+1], 0.0, L/double(N));
element.set_stiffness(E*A, E*I, 0.0);
system.add_element(element);
system.mut_elements().add(element);
}
nodes[N][1].p_mut() = F0;
system.set_p(nodes[N].y, F0);
StaticSolverLC solver(system);
solver.solve();
// Tip displacements
double ux_num = L - nodes[N][0].u();
double uy_num = nodes[N][1].u();
// Numerical tip displacements
double ux_num = L - system.get_u(nodes[N].x);
double uy_num = system.get_u(nodes[N].y);
// Reference displacements according to [1]
// Reference tip displacements according to [1]
double ux_ref = 0.5089228704;
double uy_ref = 1.2083981311;
double error_x = std::abs((ux_num - ux_ref)/ux_ref);
double error_y = std::abs((uy_num - uy_ref)/uy_ref);
// Error
// Todo: Why is the precision not better?
REQUIRE(error_x < 1.08e-3);
REQUIRE(error_y < 5.79e-4);
REQUIRE(std::abs((ux_num - ux_ref)/ux_ref) < 1.08e-3);
REQUIRE(std::abs((uy_num - uy_ref)/uy_ref) < 5.79e-4);
}
TEST_CASE("large-deformation-circular-beam")
{
/*
// Static test "Bending of a pre-curved beam into a full circle" from [1]
// [1] On the correct representation of bending and axial deformation in the absolute nodal coordinate formulation with an elastic line approach
// Johannes Gerstmayr, Hans Irschik. Journal of Sound and Vibration 318 (2008) 461-487
......@@ -81,38 +80,43 @@ TEST_CASE("large-deformation-circular-beam")
// Create nodes
for(unsigned i = 0; i < N+1; ++i)
{
DofType type = (i == 0) ? DofType::Fixed : DofType::Active;
bool active = (i != 0);
double phi = double(i)/double(N)*M_PI;
nodes.push_back(system.create_node({type, type, type}, {R*sin(phi), R*(cos(phi) - 1.0), 0.0}));
nodes.push_back(system.create_node({active, active, active}, {R*sin(phi), R*(cos(phi) - 1.0), 0.0}));
}
// Create elements
for(unsigned i = 0; i < N; ++i)
{
double dist = Node::distance(nodes[i], nodes[i+1]);
double angle = Node::angle(nodes[i], nodes[i+1]);
double dist = system.get_distance(nodes[i], nodes[i+1]);
double angle = system.get_angle(nodes[i], nodes[i+1]);
BeamElement element(nodes[i], nodes[i+1], 0.0, dist);
BeamElement element(system, nodes[i], nodes[i+1], 0.0, dist);
element.set_stiffness(EA, EI, 0.0);
element.set_reference_angles(angle - nodes[i][2].u(), angle - nodes[i+1][2].u());
element.set_reference_angles(angle - system.get_u(nodes[i].phi), angle - system.get_u(nodes[i+1].phi));
system.add_element(element);
system.mut_elements().add(element);
}
StaticSolverDC solver(system, nodes[N][2]);
StaticSolverDC solver(system, nodes[N].phi);
for(unsigned i = 0; i < 15; ++i)
solver.solve(-double(i)/15*M_PI);
double M_num = -nodes[N][2].p();
double M_ref = EI*R;
double error_M = std::abs((M_num - M_ref)/M_ref);
double error_x = std::abs(nodes[N][0].u());
double error_y = std::abs(nodes[N][1].u());
// Numerical solution
double M_num = -system.get_p(nodes[N].phi);
double x_num = system.get_u(nodes[N].x);
double y_num = system.get_u(nodes[N].y);
// Todo: Why is the error in the torque so much higher than the displacement error?
REQUIRE(error_M < 1.03e-03);
REQUIRE(error_x < 7.82e-16);
REQUIRE(error_y < 3.58e-14);
// Reference solution
double M_ref = EI*R;
double x_ref = 0.0;
double y_ref = 0.0;
// Error
// Todo: Why is the torque so much more off than the displacements?
REQUIRE(std::abs((M_num - M_ref)/M_ref) < 1.03e-03);
REQUIRE(std::abs(x_num - x_ref) < 7.82e-16);
REQUIRE(std::abs(y_num - y_ref) < 3.58e-14);
*/
}
......@@ -16,15 +16,19 @@ MatrixXd numeric_tangent_stiffness(System& system, double h = 1e-8)
for(std::size_t i = 0; i < system.dofs(); ++i)
{
double ui = system.u()(i);
VectorXd u = system.get_u();
double ui = u(i);
system.u_mut()(i) = ui + h;
q_fwd = system.q();
u(i) = ui + h;
system.set_u(u);
q_fwd = system.get_q();
system.u_mut()(i) = ui - h;
q_bwd = system.q();
u(i) = ui - h;
system.set_u(u);
q_bwd = system.get_q();
system.u_mut()(i) = ui;
u(i) = ui;
system.set_u(u);
K.col(i) = (q_fwd - q_bwd)/(2.0*h);
}
......@@ -33,11 +37,11 @@ MatrixXd numeric_tangent_stiffness(System& system, double h = 1e-8)
void check_stiffness_matrix(System& system)
{
MatrixXd K_ana = system.K();
MatrixXd K_ana = system.get_K();
MatrixXd K_num = numeric_tangent_stiffness(system);
REQUIRE(K_ana.isApprox(K_num, 1e-5)); // Equal to numeric derivatives
REQUIRE(K_ana.isApprox(K_num.transpose(), 1e-5)); // Symmetric
REQUIRE(K_ana.isApprox(K_num, 1e-5)); // Compare to numeric derivatives
REQUIRE(K_ana.isApprox(K_ana.transpose(), 1e-5)); // Check symmetry
}
TEST_CASE("tangent-stiffness-bar-element")
......@@ -48,9 +52,9 @@ TEST_CASE("tangent-stiffness-bar-element")
double EA = 100.0;
System system;
Node node0 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, { dx0, dy0, 0.0});
Node node1 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed}, {L + dx1, dy1, 0.0});
system.add_element(BarElement(node0, node1, L, EA, 0.0));
Node node0 = system.create_node({true, true, false}, { dx0, dy0, 0.0});
Node node1 = system.create_node({true, true, false}, {L + dx1, dy1, 0.0});
system.mut_elements().add(BarElement(system, node0, node1, L, EA, 0.0));
check_stiffness_matrix(system);
};
......@@ -77,12 +81,12 @@ TEST_CASE("tangent-stiffness-beam-element")
double EI = 10.0;
System system;
Node node0 = system.create_node({DofType::Active, DofType::Active, DofType::Active}, {x0, y0, phi0});
Node node1 = system.create_node({DofType::Active, DofType::Active, DofType::Active}, {x1, y1, phi1});
Node node0 = system.create_node({true, true, true}, {x0, y0, phi0});
Node node1 = system.create_node({true, true, true}, {x1, y1, phi1});
BeamElement element01(node0, node1, 0.0, L);
BeamElement element01(system, node0, node1, 0.0, L);
element01.set_stiffness(EA, EI, 0.0);
system.add_element(element01);
system.mut_elements().add(element01);
check_stiffness_matrix(system);
};
......@@ -99,23 +103,24 @@ TEST_CASE("tangent-stiffness-beam-element")
TEST_CASE("tangent-stiffness-contact-element")
{
/*
auto check_at = [](double x0, double y0, double phi0, double x1, double y1, double phi1, double x2, double y2)
{
double L = 1.0;
double EA = 100.0;
double EI = 10.0;
System system;
Node node0 = system.create_node({DofType::Active, DofType::Active, DofType::Active}, {x0, y0, phi0});
Node node1 = system.create_node({DofType::Active, DofType::Active, DofType::Active}, {x1, y1, phi1});
Node node2 = system.create_node({DofType::Active, DofType::Active, DofType::Fixed }, {x2, y2, 0.0});
Node node0 = system.create_node({true, true, true}, {x0, y0, phi0});
Node node1 = system.create_node({true, true, true}, {x1, y1, phi1});
Node node2 = system.create_node({true, true, false}, {x2, y2, 0.0});
ContactElement element(node0, node1, node2, 0.5, 0.25, 5000.0);
system.add_element(element);
ContactElement element(system, node0, node1, node2, 0.5, 0.25, {5000.0, 1e-4});
system.mut_elements().add(element);
check_stiffness_matrix(system);
};
check_at(0.0, 0.0, M_PI_2, 0.0, 1.0, M_PI_2, 1.25, 0.5); // No contact
check_at(0.0, 0.0, M_PI_2, 0.0, 1.0, M_PI_2, 0.25, 0.5); // Contact
*/
}
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