Commit 35fc35be authored by Stefan Pfeifer's avatar Stefan Pfeifer

Add static solution method and bar truss test

parent 8abbe339
......@@ -2,12 +2,7 @@
struct Dof
{
enum class Type{
Active,
Fixed,
};
Type type;
bool active;
size_t index;
bool operator!=(const Dof& rhs) const
......
......@@ -28,18 +28,18 @@ public:
}
Node create_node(std::array<double, 3> u_node, std::array<Dof::Type, 3> dof_types)
Node create_node(std::array<double, 3> u_node, std::array<bool, 3> dof_types)
{
return Node{create_dof(u_node[0], dof_types[0]),
create_dof(u_node[1], dof_types[1]),
create_dof(u_node[2], dof_types[2])};
}
Dof create_dof(double u_dof, Dof::Type type)
Dof create_dof(double u_dof, bool type)
{
switch(type)
{
case Dof::Type::Active:
case true:
{
size_t n = u.size() + 1;
u = (VectorXd(n) << u, u_dof).finished();
......@@ -49,7 +49,7 @@ public:
return Dof{type, n-1};
}
case Dof::Type::Fixed:
case false:
{
size_t n = uf.size() + 1;
uf = (VectorXd(n) << uf, u_dof).finished();
......@@ -80,7 +80,7 @@ public:
return VectorView<Dof>(
[&](Dof dof)
{
if(dof.type == Dof::Type::Active)
if(dof.active)
return u(dof.index);
else
return uf(dof.index);
......@@ -94,7 +94,7 @@ public:
return VectorView<Dof>(
[&](Dof dof)
{
if(dof.type == Dof::Type::Active)
if(dof.active)
return v(dof.index);
else
return 0.0;
......@@ -108,7 +108,7 @@ public:
return VectorView<Dof>(
[&](Dof dof)
{
if(dof.type == Dof::Type::Active)
if(dof.active)
return a(dof.index);
else
return 0.0;
......@@ -117,22 +117,56 @@ public:
nullptr);
}
VectorView<Dof> get_p()
// Todo: Should this be done via a View instead? Problem: Add vs set
void set_external_force(Dof dof, double val)
{
return VectorView<Dof>(
[&](Dof dof)
if(dof.active)
p(dof.index) = val;
}
void solve_statics_lc()
{
v.setZero();
a.setZero();
unsigned max_iter = 50; // Todo: Magic number
double epsilon = 1e-8; // Todo: Magic number
Eigen::LDLT<Eigen::MatrixXd> stiffness_dec;
Eigen::MatrixXd K(dofs(), dofs());
Eigen::VectorXd q(dofs());
get_internal_forces(q);
get_tangent_stiffness(K);
double norm_0 = (p - q).norm();
if(norm_0 == 0)
{
if(dof.type == Dof::Type::Active)
return p(dof.index);
else
return 0.0;
},
return; // System already in equilibrium
}
[&](Dof dof, double val)
for(unsigned i = 0; i < max_iter; ++i)
{
if(dof.type == Dof::Type::Active)
p(dof.index) += val;
});
stiffness_dec.compute(K);
if(stiffness_dec.info() != Eigen::Success)
{
throw std::runtime_error("Stiffness matrix decomposition failed");
}
u += stiffness_dec.solve(p - q);
update_element_states();
get_internal_forces(q);
get_tangent_stiffness(K);
double norm_i = (p - q).norm();
if(norm_i/norm_0 < epsilon)
{
return;
}
}
throw std::runtime_error("Maximum number of iterations exceeded");
}
void solve_dynamics(double timestep_factor, const std::function<bool()>& callback)
......@@ -197,7 +231,7 @@ public: // Todo: private
[&](Dof dof, double val)
{
if(dof.type == Dof::Type::Active)
if(dof.active)
M(dof.index) += val;
});
......@@ -219,7 +253,7 @@ public: // Todo: private
[&](Dof dof, double val)
{
if(dof.type == Dof::Type::Active)
if(dof.active)
q(dof.index) += val;
});
......@@ -236,7 +270,7 @@ public: // Todo: private
MatrixView<Dof> view(
[&](Dof dof_row, Dof dof_col, double val)
{
if(dof_row.type == Dof::Type::Active && dof_col.type == Dof::Type::Active)
if(dof_row.active && dof_col.active)
K(dof_row.index, dof_col.index) += val;
});
......
......@@ -7,5 +7,6 @@ CONFIG -= qt
SOURCES += tests/main.cpp \
tests/harmonic_oscillator.cpp \
tests/tangent_stiffness.cpp
tests/tangent_stiffness.cpp \
tests/bar_trusses.cpp
#include "../src/fem/System.hpp"
#include "../src/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({{ 0.0, 0.0, 0.0}}, {{false, false, false}});
Node node_02 = system.create_node({{ L, 0.0, 0.0}}, {{ true, true, false}});
Node node_03 = system.create_node({{2.0*L, 0.0, 0.0}}, {{ true, true, false}});
Node node_04 = system.create_node({{3.0*L, 0.0, 0.0}}, {{ true, true, false}});
Node node_05 = system.create_node({{4.0*L, 0.0, 0.0}}, {{ true, false, false}});
Node node_06 = system.create_node({{ 0.0, L, 0.0}}, {{ true, true, false}});
Node node_07 = system.create_node({{ L, L, 0.0}}, {{ true, true, false}});
Node node_08 = system.create_node({{2.0*L, L, 0.0}}, {{ true, true, false}});
Node node_09 = system.create_node({{3.0*L, L, 0.0}}, {{ true, true, false}});
Node node_10 = system.create_node({{4.0*L, L, 0.0}}, {{ true, true, false}});
BarElement element_01_02(node_01, node_02, L, EA, 0.0, 0.0);
BarElement element_02_03(node_02, node_03, L, EA, 0.0, 0.0);
BarElement element_03_04(node_03, node_04, L, EA, 0.0, 0.0);
BarElement element_04_05(node_04, node_05, L, EA, 0.0, 0.0);
BarElement element_06_07(node_06, node_07, L, EA, 0.0, 0.0);
BarElement element_07_08(node_07, node_08, L, EA, 0.0, 0.0);
BarElement element_08_09(node_08, node_09, L, EA, 0.0, 0.0);
BarElement element_09_10(node_09, node_10, L, EA, 0.0, 0.0);
BarElement element_01_06(node_01, node_06, L, EA, 0.0, 0.0);
BarElement element_02_07(node_02, node_07, L, EA, 0.0, 0.0);
BarElement element_03_08(node_03, node_08, L, EA, 0.0, 0.0);
BarElement element_04_09(node_04, node_09, L, EA, 0.0, 0.0);
BarElement element_05_10(node_05, node_10, L, EA, 0.0, 0.0);
BarElement element_01_07(node_01, node_07, std::sqrt(2)*L, EA, 0.0, 0.0);
BarElement element_07_03(node_07, node_03, std::sqrt(2)*L, EA, 0.0, 0.0);
BarElement element_03_09(node_03, node_09, std::sqrt(2)*L, EA, 0.0, 0.0);
BarElement element_09_05(node_09, node_05, std::sqrt(2)*L, EA, 0.0, 0.0);
system.add_element(element_01_02);
system.add_element(element_02_03);
system.add_element(element_03_04);
system.add_element(element_04_05);
system.add_element(element_06_07);
system.add_element(element_07_08);
system.add_element(element_08_09);
system.add_element(element_09_10);
system.add_element(element_01_06);
system.add_element(element_02_07);
system.add_element(element_03_08);
system.add_element(element_04_09);
system.add_element(element_05_10);
system.add_element(element_01_07);
system.add_element(element_07_03);
system.add_element(element_03_09);
system.add_element(element_09_05);
system.set_external_force(node_02.y, -F);
system.set_external_force(node_04.y, -F);
system.solve_statics_lc();
double s_numeric = -system.get_u()(node_03.y);
double s_analytic = (4.0+2.0*std::sqrt(2))*F*L/EA;
REQUIRE(std::abs(s_numeric - s_analytic) < 1e-6);
}
......@@ -16,8 +16,8 @@ TEST_CASE("Dynamic solution harmonic oscillator")
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}});
Node node_a = system.create_node({{ 0.0, 0.0, 0.0}}, {{false, false, false}});
Node node_b = system.create_node({{l + s0, 0.0, 0.0}}, {{true, false, false}});
BarElement element1(node_a, node_b, l, l*k, l*d, 0.0);
MassElement element2(node_b, m, 0.0);
......
......@@ -24,8 +24,8 @@ TEST_CASE("Tangent stiffness matrix: BarElement")
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}});
Node node0 = system.create_node({{ dx0, dy0, 0.0}}, {{true, true, false}});
Node node1 = system.create_node({{L + dx1, dy1, 0.0}}, {{true, true, false}});
BarElement element01(node0, node1, L, EA, 0.0, 0.0);
system.add_element(element01);
......
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