Commit 9fef71b9 authored by Stefan Pfeifer's avatar Stefan Pfeifer

Implement bow model setup and bracing

parent c38e0968
......@@ -75,12 +75,12 @@
"string": {
"strand_stiffness": 3100.0,
"strand_density": 0.05,
"n_strands": 10,
"n_strands": 50,
"center_mass": 0.0,
"end_mass": 0.0
},
"operation": {
"brace_height": 0.2,
"brace_height": 0.46,
"draw_length": 0.7,
"arrow_mass": 0.02
},
......
......@@ -40,4 +40,7 @@ HEADERS += \
../src/model/BowParameters.hpp \
../src/gui/Document.hpp \
../src/model/InputData.hpp \
../src/model/DiscreteLimb.hpp
../src/model/DiscreteLimb.hpp \
../src/model/BowModel.hpp \
../src/model/OutputData.hpp \
../src/numerics/SecantMethod.hpp
......@@ -124,6 +124,7 @@ public:
}
// 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)
{
if(dof.active)
......@@ -147,7 +148,7 @@ public:
return std::hypot(dx, dy);
}
// Angle between the x-axis and the line connecting node_a and node_b
// 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);
......@@ -156,13 +157,14 @@ public:
return std::atan2(dy, dx);
}
// Todo: Rename statics methods to 'solve_equilibrium' and 'solve_equilibrium_path'?
void solve_statics_lc()
{
v.setZero();
a.setZero();
unsigned max_iter = 50; // Todo: Magic number
double epsilon = 1e-8; // Todo: Magic number
double epsilon = 1e-6; // Todo: Magic number
Eigen::LDLT<Eigen::MatrixXd> stiffness_dec;
Eigen::MatrixXd K(dofs(), dofs());
......@@ -272,12 +274,20 @@ public:
};
double init_displacement = u(dof.index);
double delta_displacement = target_displacement - init_displacement;
for(unsigned i = 0; i < n_steps; ++i)
// Todo: Separate into different methods
// Todo: Think about meaning of 'n_steps'. Number of times the newton method is employed? Or number of intervals in between?
if(n_steps <= 1)
{
solve_equilibrium(init_displacement);
}
else
{
double displacement = init_displacement + delta_displacement*double(i)/double(n_steps - 1);
solve_equilibrium(displacement);
double delta_displacement = target_displacement - init_displacement;
for(unsigned i = 0; i < n_steps; ++i)
{
double displacement = init_displacement + delta_displacement*double(i)/double(n_steps - 1);
solve_equilibrium(displacement);
}
}
}
......@@ -322,7 +332,51 @@ public:
}
}
public: // Todo: private
// Todo: tangent stiffness methods only need to be public for tests. Wat do?
void get_tangent_stiffness(MatrixXd& K) const
{
K.setZero();
MatrixView<Dof> view(
[&](Dof dof_row, Dof dof_col, double val)
{
if(dof_row.active && dof_col.active)
K(dof_row.index, dof_col.index) += val;
});
for(auto e: elements)
{
e->get_tangent_stiffness(view);
}
}
// 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;
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);
}
}
private:
void update_element_states()
{
for(auto e: elements)
......@@ -374,47 +428,4 @@ public: // Todo: private
e->get_internal_forces(view);
}
}
void get_tangent_stiffness(MatrixXd& K) const
{
K.setZero();
MatrixView<Dof> view(
[&](Dof dof_row, Dof dof_col, double val)
{
if(dof_row.active && dof_col.active)
K(dof_row.index, dof_col.index) += val;
});
for(auto e: elements)
{
e->get_tangent_stiffness(view);
}
}
// 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;
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);
}
}
};
......@@ -18,7 +18,7 @@ private:
double dx;
double dy;
double L_new; // Actual length
double L_dot; // Time derivative of length
double L_dot; // Time derivative of actual length
public:
BarElement(Node nd0, Node nd1, double L, double EA, double etaA, double rhoA)
......@@ -31,6 +31,11 @@ public:
}
void set_length(double val)
{
L = val;
}
virtual void set_state(const VectorView<Dof> u, const VectorView<Dof> v)
{
dx = u(nodes[1].x) - u(nodes[0].x);
......
......@@ -6,15 +6,17 @@
#include "gui/Document.hpp"
#include "model/InputData.hpp"
#include "model/DiscreteLimb.hpp"
#include "model/BowModel.hpp"
int main()
{
//InputData data;
//data.limb.layers.push_back(Limb::Layer());
//data.save("../examples/default.bow");
InputData input;
input.load("../examples/layers.bow");
OutputData output = BowModel::simulate(input, true, false);
/*
InputData data;
data.load("../examples/layers.bow");
......@@ -24,6 +26,7 @@ int main()
{
qInfo() << limb.s[i] << ", " << limb.hc[i];
}
*/
return 0;
}
......
#pragma once
#include "InputData.hpp"
#include "OutputData.hpp"
#include "DiscreteLimb.hpp"
#include "../fem/System.hpp"
#include "../fem/elements/BeamElement.hpp"
#include "../fem/elements/BarElement.hpp"
#include "../numerics/SecantMethod.hpp"
#include <QtCore>
class BowModel
{
public:
static OutputData simulate(const InputData& input, bool statics, bool dynamics)
{
OutputData output;
BowModel model(input, output);
// Todo: Omit detailed static simulation if statics == false
model.simulate_setup();
model.simulate_statics();
if(dynamics)
model.simulate_dynamics();
return output;
}
private:
const InputData& input;
OutputData& output;
DiscreteLimb limb;
System system;
std::vector<BeamElement> elements_limb;
std::vector<BarElement> elements_string;
std::vector<Node> nodes_limb;
std::vector<Node> nodes_string;
BowModel(const InputData& input, OutputData& output)
: input(input),
output(output),
limb(input)
{
}
void simulate_setup()
{
size_t n = input.settings.n_elements_limb;
size_t k = input.settings.n_elements_string;
// Create limb nodes
for(size_t i = 0; i < n+1; ++i)
{
bool active = (i != 0);
Node node = system.create_node({{limb.x[i], limb.y[i], limb.phi[i]}}, {{active, active, active}});
nodes_limb.push_back(node);
}
// Create limb elements
for(size_t i = 0; i < n; ++i)
{
double rhoA = 0.5*(limb.rhoA[i] + limb.rhoA[i+1]);
double L = system.get_distance(nodes_limb[i], nodes_limb[i+1]);
double Cee = 0.5*(limb.Cee[i] + limb.Cee[i+1]);
double Ckk = 0.5*(limb.Ckk[i] + limb.Ckk[i+1]);
double Cek = 0.5*(limb.Cek[i] + limb.Cek[i+1]);
// 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);
BeamElement element(nodes_limb[i], nodes_limb[i+1], rhoA, L);
element.set_reference_angles(phi0, phi1);
element.set_stiffness(Cee, Ckk, Cek);
elements_limb.push_back(element);
}
// Limb tip
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;
double yc = 0.0;
// Create string nodes
for(size_t i = 0; i < k; ++i)
{
double p = double(i)/double(k);
double x = xc*(1.0 - p) + xt*p;
double y = yc*(1.0 - p) + yt*p;
Node node = system.create_node({{x, y, 0.0}}, {{true, (i != 0), false}});
nodes_string.push_back(node);
}
nodes_string.push_back(nodes_limb.back());
// Create string elements
double rhoA = double(input.string.n_strands)*input.string.strand_density;
double EA = double(input.string.n_strands)*input.string.strand_stiffness;
for(size_t i = 0; i < k; ++i)
{
BarElement element(nodes_string[i], nodes_string[i+1], 0.0, EA, 0.0, rhoA); // Element lengths are reset later when string length is determined
elements_string.push_back(element);
}
// Todo: Add elements inside loops, prevent iterator invalidation by using unique_ptr or shared_ptr
// or maybe even let system own elements. Look at boost pointer container library.
system.add_elements(elements_limb);
system.add_elements(elements_string);
// Takes a string length, iterates to equilibrium with the constraint of the brace height
// and returns the angle of the string center
auto try_string_length = [&](double string_length)
{
double L = 0.5*string_length/double(k);
for(auto& element: elements_string)
{
element.set_length(L);
}
Dof dof = nodes_string[0].x;
system.solve_statics_dc(dof, xc, 1, [](){});
return system.get_angle(nodes_string[0], nodes_string[1]) - M_PI/2;
};
// Todo: Perhaps limit the step size of the root finding algorithm ti increase robustness.
double string_length = 2.0*std::hypot(xc - xt, yc - yt);
string_length = secant_method(try_string_length, 0.95*string_length, 0.9*string_length, 1e-8, 50); // Todo: Magic numbers
// Store setup results
output.setup.limb = limb;
output.setup.string_length = string_length;
/*
for(auto& node: nodes_limb)
{
qInfo() << system.get_u()(node.x) << ", " << system.get_u()(node.y);
}
qInfo() << "=============";
for(auto& node: nodes_string)
{
qInfo() << system.get_u()(node.x) << ", " << system.get_u()(node.y);
}
*/
}
void simulate_statics()
{
qInfo() << "Simulate Statics";
}
void simulate_dynamics()
{
qInfo() << "Simulate Dynamics";
}
};
......@@ -21,7 +21,7 @@ struct DiscreteLimb
std::vector<double> y;
std::vector<double> phi;
// Properties
// Section properties
std::vector<double> Cee;
std::vector<double> Ckk;
std::vector<double> Cek;
......@@ -31,7 +31,12 @@ struct DiscreteLimb
DiscreteLimb(const InputData& input)
{
calculate_nodes(input);
calculate_properties(input);
calculate_section_properties(input);
}
DiscreteLimb() // Todo: Remove this, find better solution
{
}
private:
......@@ -67,10 +72,11 @@ private:
integrate_n_steps(stepper, system, z0, s0, ds, n, observer);
}
void calculate_properties(const InputData& input)
void calculate_section_properties(const InputData& input)
{
// Construct splines for layer geometry
SplineFunction layer_width(input.limb.width);
std::vector<SplineFunction> layer_heights;
for(size_t i = 0; i < input.limb.layers.size(); ++i)
{
......@@ -78,7 +84,7 @@ private:
}
// Calculate section properties at every node
for(size_t i = 0; i < s.size(); ++i) // Iterate over nodes
for(size_t i = 0; i < s.size(); ++i)
{
Cee.push_back(0.0);
Ckk.push_back(0.0);
......@@ -86,12 +92,12 @@ private:
rhoA.push_back(0.0);
hc.push_back(0.0);
double yu = 0.0; // Upper Position layer top
double yl = 0.0; // Position layer bottom
double p = s[i]/s.back();
double p = s[i]/s.back(); // Relative position along the limb curve
double yu = 0.0; // Upper border of the layer
double yl = 0.0; // Lower border of the layer
for(size_t j = 0; j < input.limb.layers.size(); ++j) // Iterate over layers
// Iterate over layers and sum up section properties
for(size_t j = 0; j < input.limb.layers.size(); ++j)
{
if(p >= layer_heights[j].arg_min() && p <= layer_heights[j].arg_max()) // Todo: Abstract this?
{
......
#pragma once
#include "DiscreteLimb.hpp"
struct BowOutput
{
/*
pub time: Vec<f64>,
pub draw_force: Vec<f64>,
pub pos_limb: Vec<Vec<[f64; 2]>>,
pub vel_limb: Vec<Vec<[f64; 2]>>,
pub acc_limb: Vec<Vec<[f64; 2]>>,
pub pos_string: Vec<Vec<[f64; 2]>>,
pub vel_string: Vec<Vec<[f64; 2]>>,
pub acc_string: Vec<Vec<[f64; 2]>>,
*/
};
struct BowSetup
{
DiscreteLimb limb;
double string_length;
};
struct OutputData
{
BowSetup setup;
BowOutput statics;
BowOutput dynamics;
};
#pragma once
#include <functional>
#include <stdexcept>
// Find the root of the function f by using the secant method [1] with the initial values xa, xb
// [1] https://en.wikipedia.org/wiki/Secant_method
// Todo: Can the lambda be called statically?
template<class F>
double secant_method(const F& f, double x0, double x1, double epsilon, unsigned iter)
{
double f0 = f(x0);
double f1 = f(x1);
for(unsigned i = 0; i < iter; ++i)
{
if(std::abs(f1) < epsilon)
{
return x1;
}
double x2 = x1 - (x1 - x0)/(f1 - f0)*f1;
x0 = x1; f0 = f1;
x1 = x2; f1 = f(x2);
}
throw std::runtime_error("Secant Method: Maximum number of iterations exceeded");
}
......@@ -63,10 +63,6 @@ TEST_CASE("Tangent stiffness matrix: BeamElement")
element01.set_stiffness(EA, EI, 0.0);
system.add_element(element01);
system.update_element_states();
system.update_element_states();
system.update_element_states();
check_stiffness_matrix(system);
};
......
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