Commit 833a0378 authored by Stefan Pfeifer's avatar Stefan Pfeifer

Implement tangent stiffness matrix

parent 2f16b46c
......@@ -19,9 +19,9 @@ HEADERS += \
src/gui/MainWindow.hpp \
src/fem/System.hpp \
src/fem/elements/Element.hpp \
src/fem/View.hpp \
src/fem/elements/BarElement.hpp \
src/fem/Node.hpp \
src/fem/elements/MassElement.hpp
src/fem/elements/MassElement.hpp \
src/fem/View.hpp
INCLUDEPATH += /home/s/Libraries/Eigen-3.2.8
......@@ -9,6 +9,11 @@ struct Dof
Type type;
size_t index;
bool operator!=(const Dof& rhs) const
{
return rhs.index != index;
}
};
struct Node
......
......@@ -83,10 +83,7 @@ public:
return uf(dof.index);
},
[&](Dof, double)
{
// Todo: Unreachable
});
nullptr);
}
const VectorView<Dof> get_v() const
......@@ -100,10 +97,7 @@ public:
return 0.0;
},
[&](Dof, double)
{
// Todo: Unreachable
});
nullptr);
}
const VectorView<Dof> get_a() const
......@@ -117,65 +111,26 @@ public:
return 0.0;
},
[&](Dof, double)
{
// Todo: Unreachable
});
nullptr);
}
/*
VectorView<Dof> get_u()
VectorView<Dof> get_p()
{
return VectorView<Dof>(
[&](Dof dof)
{
if(dof.type == Dof::Type::Active)
return u(dof.index);
return p(dof.index);
else
return uf(dof.index);
return 0.0;
},
[&](Dof dof, double val)
{
if(dof.type == Dof::Type::Active)
u(dof.index) = val;
else
uf(dof.index) = val;
p(dof.index) += val;
});
}
*/
/*
const VectorView<const VectorXd> get_u() const
{
VectorAccess<const VectorXd> access(&u, &uf);
return VectorView<const VectorXd>(access);
}
const VectorView<const VectorXd> get_v() const
{
VectorAccess<const VectorXd> access(&v, nullptr);
return VectorView<const VectorXd>(access);
}
const VectorView<const VectorXd> get_a() const
{
VectorAccess<const VectorXd> access(&a, nullptr);
return VectorView<const VectorXd>(access);
}
VectorView<VectorXd> get_p()
{
VectorAccess<VectorXd> access(&p, nullptr);
return VectorView<VectorXd>(access);
}
const VectorView<const VectorXd> get_p() const
{
VectorAccess<const VectorXd> access(&p, nullptr);
return VectorView<const VectorXd>(access);
}
*/
void solve_dynamics(double dt, const std::function<bool()>& callback)
{
......@@ -209,7 +164,7 @@ public:
}
}
private:
public: // Todo: private
void update_element_states()
{
for(auto e: elements)
......@@ -231,7 +186,7 @@ private:
[&](Dof dof, double val)
{
if(dof.type == Dof::Type::Active)
M(dof.index) = val;
M(dof.index) += val;
});
for(auto e: elements)
......@@ -253,7 +208,7 @@ private:
[&](Dof dof, double val)
{
if(dof.type == Dof::Type::Active)
q(dof.index) = val;
q(dof.index) += val;
});
for(auto e: elements)
......@@ -261,4 +216,21 @@ private:
e->get_internal_forces(view);
}
}
void get_tangent_stiffness(MatrixXd& K)
{
K.setZero();
MatrixView<Dof> view(
[&](Dof dof_row, Dof dof_col, double val)
{
if(dof_row.type == Dof::Type::Active && dof_col.type == Dof::Type::Active)
K(dof_row.index, dof_col.index) += val;
});
for(auto e: elements)
{
e->get_tangent_stiffness(view);
}
}
};
......@@ -14,6 +14,8 @@ using Matrix = Eigen::Matrix<double, n, n>;
template<std::size_t n>
using Vector = Eigen::Matrix<double, n, 1>;
//================================================================================================
template<class Key>
class LocalVectorEntry
{
......@@ -128,98 +130,24 @@ public:
}
};
/*
template<class Key>
class VectorAccess
{
virtual double get(Key) const = 0;
virtual void add(Key, double) = 0;
};
//================================================================================================
template<class Key>
class MatrixAccess
{
virtual double get(Key, Key) const = 0;
virtual void add(Key, Key, double) = 0;
};
*/
/*
template<class Access, class Key>
class LocalVectorEntry
class LocalMatrixEntry
{
private:
Access access;
Key key;
public:
LocalVectorEntry(Access access, Key key)
: access(access), key(key)
{
}
// Read
operator double() const
{
return access.get(key);
}
std::function<void(Key, Key, double)>& add;
Key key_row;
Key key_col;
// Write
// Write symmetrically
void operator+=(double rhs)
{
access.add(key, rhs);
}
add(key_row, key_col, rhs);
void operator-=(double rhs)
{
operator+=(-rhs);
}
};
*/
/*
template<class Access, class Key, size_t N>
class LocalVectorView
{
private:
Access access;
const std::array<Key, N>& keys;
public:
LocalVectorView(Access access, const std::array<Key, N>& keys)
: access(access), keys(keys)
{
}
// Read
operator Vector<N>() const
{
Vector<N> res;
for(size_t i = 0; i < N; ++i)
if(key_row != key_col)
{
res(i) = access.get(keys[i]);
}
return res;
}
// Write
template<class T>
void operator+=(const Eigen::MatrixBase<T>& rhs)
{
for(size_t i = 0; i < N; ++i)
{
access.add(keys[i], rhs(i));
add(key_col, key_row, rhs);
}
}
......@@ -227,123 +155,15 @@ public:
{
operator+=(-rhs);
}
// Index
LocalVectorEntry<Access, Key> operator()(size_t i)
{
return LocalVectorEntry<Access, Key>(access, keys[i]);
}
const LocalVectorEntry<Access, Key> operator()(size_t i) const
{
return LocalVectorEntry<Access, Key>(access, keys[i]);
}
};
*/
/*
template<class Access, class Key>
class VectorViewBase
{
private:
Access access;
public:
VectorViewBase(Access access)
: access(access)
{
}
// Index
template<size_t N>
LocalVectorView<Access, Key, N> operator()(const std::array<Key, N>& keys)
{
return LocalVectorView<Access, Key, N>(access, keys);
}
template<size_t N>
const LocalVectorView<Access, Key, N> operator()(const std::array<Key, N>& keys) const
{
return LocalVectorView<Access, Key, N>(access, keys);
}
LocalVectorEntry<Access, Key> operator()(Key key)
{
return LocalVectorEntry<Access, Key>(access, key);
}
const LocalVectorEntry<Access, Key> operator()(Key key) const
{
return LocalVectorEntry<Access, Key>(access, key);
}
};
*/
/*
template<class Access, class Key>
class LocalMatrixEntry
{
private:
Access access;
Key key_row;
Key key_col;
public:
LocalMatrixEntry(Access access, Key key_row, Key key_col)
: access(access), key_row(key_row), key_col(key_col)
{
}
// Read
operator double() const
{
return access.get(key_row, key_col);
}
// Write
void operator+=(double rhs)
{
access.add(key_row, key_col, rhs);
}
void operator-=(double rhs)
{
operator+=(-rhs);
}
};
*/
/*
template<class Access, class Key, size_t N>
template<class Key, size_t N>
class LocalMatrixView
{
private:
Access access;
const std::array<Key, N>& keys;
public:
LocalMatrixView(Access access, const std::array<Key, N>& keys)
: access(access), keys(keys)
{
}
// Read
operator Matrix<N>() const
{
Matrix<N> res;
for(size_t j = 0; j < N; ++j)
{
for(size_t i = 0; i < N; ++i)
{
res(i, j) = access.get(keys[i], keys[j]);
}
}
return res;
}
std::function<void(Key, Key, double)>& add;
const std::array<Key, N>& keys;
// Write
template<class T>
......@@ -353,7 +173,7 @@ public:
{
for(size_t i = 0; i < N; ++i)
{
access.add(keys[i], keys[j], rhs(i, j));
add(keys[i], keys[j], rhs(i, j));
}
}
}
......@@ -364,53 +184,50 @@ public:
}
// Index
LocalMatrixEntry<Access, Key> operator()(size_t i, size_t j)
LocalMatrixEntry<Key> operator()(size_t i, size_t j)
{
return LocalMatrixEntry<Access, Key>(access, keys[i], keys[j]);
return LocalMatrixEntry<Key>{add, keys[i], keys[j]};
}
const LocalMatrixEntry<Access, Key> operator()(size_t i, size_t j) const
const LocalMatrixEntry<Key> operator()(size_t i, size_t j) const
{
return LocalMatrixEntry<Access, Key>(access, keys[i], keys[j]);
return LocalMatrixEntry<Key>{add, keys[i], keys[j]};
}
};
*/
/*
template<class Access, class Key>
class MatrixViewBase
template<class Key>
class MatrixView
{
private:
Access access;
mutable std::function<void(Key, Key, double)> add;
public:
MatrixViewBase(Access access)
: access(access)
MatrixView(std::function<void(Key, Key, double)> add)
: add(add)
{
}
// Index
template<size_t N>
LocalMatrixView<Access, Key, N> operator()(const std::array<Key, N>& keys)
LocalMatrixView<Key, N> operator()(const std::array<Key, N>& keys)
{
return LocalMatrixView<Access, Key, N>(access, keys);
return LocalMatrixView<Key, N>{add, keys};
}
template<size_t N>
const LocalMatrixView<Access, Key, N> operator()(const std::array<Key, N>& keys) const
const LocalMatrixView<Key, N> operator()(const std::array<Key, N>& keys) const
{
return LocalMatrixView<Access, Key, N>(access, keys);
return LocalMatrixView<Key, N>{add, keys};
}
LocalMatrixEntry<Access, Key> operator()(Key key_row, Key key_col)
LocalMatrixEntry<Key> operator()(Key key_row, Key key_col)
{
return LocalMatrixEntry<Access, Key>(access, key_row, key_col);
return LocalMatrixEntry<Key>{add, key_row, key_col};
}
const LocalMatrixEntry<Access, Key> operator()(Key key_row, Key key_col) const
const LocalMatrixEntry<Key> operator()(Key key_row, Key key_col) const
{
return LocalMatrixEntry<Access, Key>(access, key_row, key_col);
return LocalMatrixEntry<Key>{add, key_row, key_col};
}
};
*/
......@@ -60,6 +60,30 @@ public:
q(nodes[1].y) += N*dy/L;
}
virtual void get_tangent_stiffness(MatrixView<Dof> K) const override
{
double c0 = EA*(L_new - L)/(L_new*L);
double c1 = EA/(L*L*L);
// col 0
K(nodes[0].x, nodes[0].x) += c1*dx*dx + c0;
// col 1
K(nodes[0].x, nodes[0].y) += c1*dx*dy;
K(nodes[0].y, nodes[0].y) += c1*dy*dy + c0;
// col 2
K(nodes[0].x, nodes[1].x) += -c1*dx*dx - c0;
K(nodes[0].y, nodes[1].x) += -c1*dx*dy;
K(nodes[1].x, nodes[1].x) += c1*dx*dx + c0;
// col 3
K(nodes[0].x, nodes[1].y) += -c1*dx*dy;
K(nodes[0].y, nodes[1].y) += -c1*dy*dy - c0;
K(nodes[1].x, nodes[1].y) += c1*dx*dy;
K(nodes[1].y, nodes[1].y) += c1*dy*dy + c0;
}
virtual double get_potential_energy() const override
{
return 0.5*EA/L*std::pow(L_new - L, 2);
......@@ -69,107 +93,4 @@ public:
{
return EA/L*(L_new - L) + etaA/L*L_dot;
}
};
/*
pub struct BarElement
{
nodes: [RcMut<Node2d>; 2],
L: f64, // Length
rhoA: f64,
etaA: f64,
EA: f64, // Longitudinial stiffness
}
impl BarElement
{
pub fn new(node0: &RcMut<Node2d>, node1: &RcMut<Node2d>, L: f64, rhoA: f64, etaA: f64, EA: f64) -> BarElement
{
assert!(L != 0.0); // Todo: Use more asserts when there are numbers with restricted domains
BarElement
{
nodes: [node0.clone(), node1.clone()],
L: L,
rhoA: rhoA,
etaA: etaA,
EA: EA,
}
}
pub fn get_length(&self) -> f64
{
self.L
}
pub fn set_length(&mut self, L: f64)
{
self.L = L;
}
}
impl Element for BarElement
{
fn mass_matrix(&mut self)
{
let m = 0.5*self.rhoA*self.L;
self.nodes[0].dofs[0].m += m;
self.nodes[0].dofs[1].m += m;
self.nodes[1].dofs[0].m += m;
self.nodes[1].dofs[1].m += m;
}
fn internal_forces(&mut self, K: Option<&mut MatrixView>)
{
let dx = self.nodes[1].dofs[0].u - self.nodes[0].dofs[0].u;
let dy = self.nodes[1].dofs[1].u - self.nodes[0].dofs[1].u;
let L = dx.hypot(dy);
let dot_L = 1.0/L*(dx*(self.nodes[1].dofs[0].v - self.nodes[0].dofs[0].v) + dy*(self.nodes[1].dofs[1].v - self.nodes[0].dofs[1].v));
let N = self.EA/self.L*(L - self.L) + self.etaA/self.L*dot_L; // Axial force
self.nodes[0].dofs[0].q -= N*dx/L;
self.nodes[0].dofs[1].q -= N*dy/L;
self.nodes[1].dofs[0].q += N*dx/L;
self.nodes[1].dofs[1].q += N*dy/L;
if let Some(K) = K
{
let c0 = self.EA*(L - self.L)/(L*self.L); // Todo: Unify computation with axial force
let c1 = self.EA/L.powi(3);
// col 0
K.add_sym(&self.nodes[0].dofs[0], &self.nodes[0].dofs[0], c1*dx*dx + c0);
// col 1
K.add_sym(&self.nodes[0].dofs[0], &self.nodes[0].dofs[1], c1*dx*dy);
K.add_sym(&self.nodes[0].dofs[1], &self.nodes[0].dofs[1], c1*dy*dy + c0);
// col 2
K.add_sym(&self.nodes[0].dofs[0], &self.nodes[1].dofs[0], -c1*dx*dx - c0);
K.add_sym(&self.nodes[0].dofs[1], &self.nodes[1].dofs[0], -c1*dx*dy);
K.add_sym(&self.nodes[1].dofs[0], &self.nodes[1].dofs[0], c1*dx*dx + c0);
// col 3
K.add_sym(&self.nodes[0].dofs[0], &self.nodes[1].dofs[1], -c1*dx*dy);
K.add_sym(&self.nodes[0].dofs[1], &self.nodes[1].dofs[1], -c1*dy*dy - c0);
K.add_sym(&self.nodes[1].dofs[0], &self.nodes[1].dofs[1], c1*dx*dy);
K.add_sym(&self.nodes[1].dofs[1], &self.nodes[1].dofs[1], c1*dy*dy + c0);
}
}
fn potential_energy(&self, V: &mut f64)
{
// Todo: Perhaps store these instead of recalculating or at least make a function to prevent code repetition
let dx = self.nodes[1].dofs[0].u - self.nodes[0].dofs[0].u;
let dy = self.nodes[1].dofs[1].u - self.nodes[0].dofs[1].u;
let L = dx.hypot(dy);
*V += 0.5*self.EA/self.L*(L - self.L).powi(2);
}
}
*/
......@@ -13,7 +13,7 @@ public:
virtual void get_masses(VectorView<Dof> M) const = 0;
virtual void get_internal_forces(VectorView<Dof> q) const = 0;
//virtual void get_tangent_stiffness(MatrixView<Dof> q) const = 0;
virtual void get_tangent_stiffness(MatrixView<Dof> K) const = 0;
virtual double get_potential_energy() const = 0;
double get_kinetic_energy(const VectorView<Dof> v) const;
......
......@@ -35,6 +35,11 @@ public:
}
virtual void get_tangent_stiffness(MatrixView<Dof> K) const override
{
}
virtual double get_potential_energy() const override
{
return 0.0;
......
......@@ -14,7 +14,6 @@ void test_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}});
......@@ -33,7 +32,7 @@ void test_harmonic_oscillator()
//assert!(delta < omega0); // Make sure the system is underdamped // Todo: Reintroduce this check
system.solve_dynamics(1e-7, [&]()
system.solve_dynamics(1e-5, [&]()
{
// Time
double t = system.get_time();
......@@ -52,7 +51,7 @@ void test_harmonic_oscillator()
double error_s = std::abs(s_num - s_ref);
double error_v = std::abs(v_num - v_ref);
std::cout << error_s << "\n";
//std::cout << error_s << "\n";
return t < T;
});
......
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