Commit 25693a0c authored by Stefan Pfeifer's avatar Stefan Pfeifer

Change interface for getting u, v, a from System

parent 22264d88
......@@ -13,8 +13,7 @@ class System
private:
std::vector<Element*> elements;
double t;
double t; // Time
VectorXd u; // Displacements
VectorXd v; // Velocities
VectorXd a; // Accelerations
......@@ -23,7 +22,15 @@ private:
VectorXd uf; // Displacements of fixed DOFs
public:
System(): t(0.0)
const VectorView<Dof> get_u; // Todo: Rename to view_u, view_v, view_a or something better?
const VectorView<Dof> get_v;
const VectorView<Dof> get_a;
System()
: t(0.0),
get_u([&](Dof dof){ return dof.active ? u(dof.index) : uf(dof.index); }, nullptr),
get_v([&](Dof dof){ return dof.active ? v(dof.index) : 0.0; }, nullptr),
get_a([&](Dof dof){ return dof.active ? a(dof.index) : 0.0; }, nullptr)
{
}
......@@ -58,7 +65,7 @@ public:
void add_element(Element& element)
{
element.set_state(get_u(), get_v());
element.set_state(get_u, get_v);
elements.push_back(&element);
}
......@@ -81,48 +88,6 @@ public:
return t;
}
const VectorView<Dof> get_u() const
{
return VectorView<Dof>(
[&](Dof dof)
{
if(dof.active)
return u(dof.index);
else
return uf(dof.index);
},
nullptr);
}
const VectorView<Dof> get_v() const
{
return VectorView<Dof>(
[&](Dof dof)
{
if(dof.active)
return v(dof.index);
else
return 0.0;
},
nullptr);
}
const VectorView<Dof> get_a() const
{
return VectorView<Dof>(
[&](Dof dof)
{
if(dof.active)
return a(dof.index);
else
return 0.0;
},
nullptr);
}
// Todo: Should this be done via a View instead? Problem: Add vs set
// One solution: Replace get, add with get, set. Implement add as set(get+add)
double get_external_force(Dof dof)
......@@ -142,8 +107,8 @@ public:
// Euclidean distance between node_a and node_b
double get_distance(Node node_a, Node node_b)
{
double dx = get_u()(node_b.x) - get_u()(node_a.x);
double dy = get_u()(node_b.y) - get_u()(node_a.y);
double dx = get_u(node_b.x) - get_u(node_a.x);
double dy = get_u(node_b.y) - get_u(node_a.y);
return std::hypot(dx, dy);
}
......@@ -151,8 +116,8 @@ public:
// Angle between the x-axis and the line connecting node_a to node_b
double get_angle(Node node_a, Node node_b)
{
double dx = get_u()(node_b.x) - get_u()(node_a.x);
double dy = get_u()(node_b.y) - get_u()(node_a.y);
double dx = get_u(node_b.x) - get_u(node_a.x);
double dy = get_u(node_b.y) - get_u(node_a.y);
return std::atan2(dy, dx);
}
......@@ -381,7 +346,7 @@ private:
{
for(auto e: elements)
{
e->set_state(get_u(), get_v());
e->set_state(get_u, get_v);
}
}
......
......@@ -81,8 +81,8 @@ private:
// Todo: Document this
double phi = system.get_angle(nodes_limb[i], nodes_limb[i+1]);
double phi0 = phi - system.get_u()(nodes_limb[i].phi);
double phi1 = phi - system.get_u()(nodes_limb[i+1].phi);
double phi0 = phi - system.get_u(nodes_limb[i].phi);
double phi1 = phi - system.get_u(nodes_limb[i+1].phi);
BeamElement element(nodes_limb[i], nodes_limb[i+1], rhoA, L);
element.set_reference_angles(phi0, phi1);
......@@ -91,8 +91,8 @@ private:
}
// Limb tip
double xt = system.get_u()(nodes_limb.back().x);
double yt = system.get_u()(nodes_limb.back().y);
double xt = system.get_u(nodes_limb.back().x);
double yt = system.get_u(nodes_limb.back().y);
// String center at brace height
double xc = input.operation.brace_height;
......@@ -226,8 +226,8 @@ private:
{
states.time.push_back(system.get_time());
states.draw_force.push_back(system.get_external_force(nodes_string[0].x));
states.pos_string.push_back(system.get_u()(nodes_string[0].x));
states.pos_arrow.push_back(system.get_u()(node_arrow.x));
states.pos_string.push_back(system.get_u(nodes_string[0].x));
states.pos_arrow.push_back(system.get_u(node_arrow.x));
//qInfo() << states.pos_string.back() << ", " << states.draw_force.back();
qInfo() << states.time.back() << ", " << states.pos_arrow.back() - states.pos_string.back();
......
......@@ -66,7 +66,7 @@ TEST_CASE("Small deformation bar truss")
system.solve_statics_lc();
double s_numeric = -system.get_u()(node_03.y);
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);
......@@ -94,7 +94,7 @@ TEST_CASE("Large deformation bar truss")
// Solve for static equilibrium and get vertical deflection of node 02
system.solve_statics_dc(node02.y, 0.6*H, 100, [&]()
{
double s = system.get_u()(node02.y);
double s = system.get_u(node02.y);
double f_num = system.get_external_force(node02.y);
double L0 = M_SQRT2*H;
......
......@@ -40,8 +40,8 @@ TEST_CASE("Dynamic solution harmonic oscillator")
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);
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));
......
......@@ -47,8 +47,8 @@ TEST_CASE("Large deformation cantilever")
system.solve_statics_lc();
// Tip displacements
double ux_num = L - system.get_u()(nodes[N].x);
double uy_num = system.get_u()(nodes[N].y);
double ux_num = L - system.get_u(nodes[N].x);
double uy_num = system.get_u(nodes[N].y);
// Reference displacements according to [1]
double ux_ref = 0.5089228704;
......@@ -94,8 +94,8 @@ TEST_CASE("Large deformation circular beam")
BeamElement element(nodes[i], nodes[i+1], 0.0, dist);
element.set_stiffness(EA, EI, 0.0);
element.set_reference_angles(angle - system.get_u()(nodes[i].phi),
angle - system.get_u()(nodes[i+1].phi));
element.set_reference_angles(angle - system.get_u(nodes[i].phi),
angle - system.get_u(nodes[i+1].phi));
elements.push_back(element);
}
......@@ -111,8 +111,8 @@ TEST_CASE("Large deformation circular beam")
double error_M = std::abs((M_num - M_ref)/M_ref);
double error_x = std::abs(system.get_u()(nodes[N].x));
double error_y = std::abs(system.get_u()(nodes[N].y));
double error_x = std::abs(system.get_u(nodes[N].x));
double error_y = std::abs(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.1e-3);
......
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