diff --git a/src/matrices.hh b/src/matrices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fc51d472d5b2a359d37ebccef9b45ecb6411e137
--- /dev/null
+++ b/src/matrices.hh
@@ -0,0 +1,9 @@
+#ifndef SRC_MATRICES_HH
+#define SRC_MATRICES_HH
+
+template <class Matrix> struct Matrices {
+ Matrix elasticity;
+ Matrix damping;
+ Matrix mass;
+};
+#endif
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 6c75bb1711730dc35cc82eacba0428b01bffd201..da973e581b769509697e3b9e6f3584a089522ae7 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -59,6 +59,7 @@
#include "enums.hh"
#include "friction_writer.hh"
#include "gridselector.hh"
+#include "matrices.hh"
#include "sand-wedge-data/mybody.hh"
#include "sand-wedge-data/mygeometry.hh"
#include "sand-wedge-data/myglobalfrictiondata.hh"
@@ -188,12 +189,13 @@ int main(int argc, char *argv[]) {
MyBody<dims> const body(parset);
- Matrix A, C, M;
+ Matrices<Matrix> matrices;
myAssembler.assembleElasticity(body.getYoungModulus(),
- body.getPoissonRatio(), A);
+ body.getPoissonRatio(), matrices.elasticity);
myAssembler.assembleViscosity(body.getShearViscosityField(),
- body.getBulkViscosityField(), C);
- myAssembler.assembleMass(body.getDensityField(), M);
+ body.getBulkViscosityField(),
+ matrices.damping);
+ myAssembler.assembleMass(body.getDensityField(), matrices.mass);
ScalarMatrix frictionalBoundaryMass;
myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
@@ -252,7 +254,7 @@ int main(int argc, char *argv[]) {
// automatically attained in the case [v = 0 = a]. Assuming
// dPhi(v = 0) = 0, we thus only have to solve Au = ell0
Vector u_initial(leafVertexCount);
- solveLinearProblem(dirichletNodes, A, ell0, u_initial,
+ solveLinearProblem(dirichletNodes, matrices.elasticity, ell0, u_initial,
parset.sub("u0.solver"));
ScalarVector normalStress(leafVertexCount);
@@ -271,8 +273,9 @@ int main(int argc, char *argv[]) {
// constraints), again assuming dPhi(v = 0) = 0
Vector a_initial(leafVertexCount);
Vector accelerationRHS = ell0;
- Arithmetic::subtractProduct(accelerationRHS, A, u_initial);
- solveLinearProblem(noNodes, M, accelerationRHS, a_initial,
+ Arithmetic::subtractProduct(accelerationRHS, matrices.elasticity,
+ u_initial);
+ solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a_initial,
parset.sub("a0.solver"));
// }}}
@@ -356,7 +359,7 @@ int main(int argc, char *argv[]) {
parset.get<double>("boundary.friction.L"),
parset.get<double>("boundary.friction.V0")),
initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
- velocityDirichletFunction, dirichletNodes, M, A, C,
+ velocityDirichletFunction, dirichletNodes, matrices,
u_initial, v_initial, a_initial));
auto const refinementTolerance =
diff --git a/src/timestepping.hh b/src/timestepping.hh
index 524e76997c2ad92582e5c7dd370d751b5673079b..5547007512a3b1c7d14a9c810d5dd608f2cb10fc 100644
--- a/src/timestepping.hh
+++ b/src/timestepping.hh
@@ -6,6 +6,7 @@
#include <dune/common/bitsetvector.hh>
#include "enums.hh"
+#include "matrices.hh"
template <class Vector, class Matrix, class Function, size_t dim>
class TimeSteppingScheme {
@@ -31,19 +32,18 @@ std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dimension>>
initTimeStepper(Config::scheme scheme,
Function const &velocityDirichletFunction,
Dune::BitSetVector<dimension> const &velocityDirichletNodes,
- Matrix const &massMatrix, Matrix const &stiffnessMatrix,
- Matrix const &dampingMatrix, Vector const &u_initial,
+ Matrices<Matrix> const &matrices, Vector const &u_initial,
Vector const &v_initial, Vector const &a_initial) {
switch (scheme) {
case Config::Newmark:
return std::make_shared<Newmark<Vector, Matrix, Function, dimension>>(
- stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
- a_initial, velocityDirichletNodes, velocityDirichletFunction);
+ matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
+ velocityDirichletFunction);
case Config::BackwardEuler:
return std::make_shared<
BackwardEuler<Vector, Matrix, Function, dimension>>(
- stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
- velocityDirichletNodes, velocityDirichletFunction);
+ matrices, u_initial, v_initial, velocityDirichletNodes,
+ velocityDirichletFunction);
default:
assert(false);
}
diff --git a/src/timestepping/backward_euler.cc b/src/timestepping/backward_euler.cc
index 73ecfbd39faf530178a329f9c4557e05efa82e2d..f3e5158848ad2ad17e1578ad550f61c5b22b75fb 100644
--- a/src/timestepping/backward_euler.cc
+++ b/src/timestepping/backward_euler.cc
@@ -2,13 +2,10 @@
template <class Vector, class Matrix, class Function, size_t dim>
BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
- Matrix const &_A, Matrix const &_M, Matrix const &_C,
- Vector const &_u_initial, Vector const &_v_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
+ Matrices<Matrix> const &_matrices, Vector const &_u_initial,
+ Vector const &_v_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
- : A(_A),
- M(_M),
- C(_C),
+ : matrices(_matrices),
u(_u_initial),
v(_v_initial),
dirichletNodes(_dirichletNodes),
@@ -51,22 +48,23 @@ void BackwardEuler<Vector, Matrix, Function, dim>::setup(
// set up LHS (for fixed tau, we'd only really have to do this once)
{
- Dune::MatrixIndexSet indices(A.N(), A.M());
- indices.import(A);
- indices.import(M);
- indices.import(C);
+ Dune::MatrixIndexSet indices(matrices.elasticity.N(),
+ matrices.elasticity.M());
+ indices.import(matrices.elasticity);
+ indices.import(matrices.mass);
+ indices.import(matrices.damping);
indices.exportIdx(AM);
}
AM = 0.0;
- Arithmetic::addProduct(AM, 1.0 / tau, M);
- Arithmetic::addProduct(AM, 1.0, C);
- Arithmetic::addProduct(AM, tau, A);
+ Arithmetic::addProduct(AM, 1.0 / tau, matrices.mass);
+ Arithmetic::addProduct(AM, 1.0, matrices.damping);
+ Arithmetic::addProduct(AM, tau, matrices.elasticity);
// set up RHS
{
rhs = ell;
- Arithmetic::addProduct(rhs, 1.0 / tau, M, v_o);
- Arithmetic::subtractProduct(rhs, A, u_o);
+ Arithmetic::addProduct(rhs, 1.0 / tau, matrices.mass, v_o);
+ Arithmetic::subtractProduct(rhs, matrices.elasticity, u_o);
}
iterate = v_o;
diff --git a/src/timestepping/backward_euler.hh b/src/timestepping/backward_euler.hh
index 1e17547616bbda18d56e07ee4b3eda7ef1e7ba8a..db1a8402d279469663b48108cd571aadfd284eea 100644
--- a/src/timestepping/backward_euler.hh
+++ b/src/timestepping/backward_euler.hh
@@ -4,8 +4,8 @@
template <class Vector, class Matrix, class Function, size_t dim>
class BackwardEuler : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
public:
- BackwardEuler(Matrix const &_A, Matrix const &_M, Matrix const &_C,
- Vector const &_u_initial, Vector const &_v_initial,
+ BackwardEuler(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
+ Vector const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction);
@@ -21,9 +21,7 @@ class BackwardEuler : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
const;
private:
- Matrix const &A;
- Matrix const &M;
- Matrix const &C;
+ Matrices<Matrix> const &matrices;
Vector u;
Vector v;
Dune::BitSetVector<dim> const &dirichletNodes;
diff --git a/src/timestepping/newmark.cc b/src/timestepping/newmark.cc
index 867f808b90b3baa128b390d69992c15a58adde10..ee81b31005e4e66c9bfb568adfc8e7e4bc6d8f00 100644
--- a/src/timestepping/newmark.cc
+++ b/src/timestepping/newmark.cc
@@ -2,13 +2,11 @@
template <class Vector, class Matrix, class Function, size_t dim>
Newmark<Vector, Matrix, Function, dim>::Newmark(
- Matrix const &_A, Matrix const &_M, Matrix const &_C,
- Vector const &_u_initial, Vector const &_v_initial,
- Vector const &_a_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
+ Matrices<Matrix> const &_matrices, Vector const &_u_initial,
+ Vector const &_v_initial, Vector const &_a_initial,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
- : A(_A),
- M(_M),
- C(_C),
+ : matrices(_matrices),
u(_u_initial),
v(_v_initial),
a(_a_initial),
@@ -56,24 +54,25 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
// set up LHS (for fixed tau, we'd only really have to do this once)
{
- Dune::MatrixIndexSet indices(A.N(), A.M());
- indices.import(A);
- indices.import(M);
- indices.import(C);
+ Dune::MatrixIndexSet indices(matrices.elasticity.N(),
+ matrices.elasticity.M());
+ indices.import(matrices.elasticity);
+ indices.import(matrices.mass);
+ indices.import(matrices.damping);
indices.exportIdx(AM);
}
AM = 0.0;
- Arithmetic::addProduct(AM, 2.0 / tau, M);
- Arithmetic::addProduct(AM, 1.0, C);
- Arithmetic::addProduct(AM, tau / 2.0, A);
+ Arithmetic::addProduct(AM, 2.0 / tau, matrices.mass);
+ Arithmetic::addProduct(AM, 1.0, matrices.damping);
+ Arithmetic::addProduct(AM, tau / 2.0, matrices.elasticity);
// set up RHS
{
rhs = ell;
- Arithmetic::addProduct(rhs, 2.0 / tau, M, v_o);
- Arithmetic::addProduct(rhs, M, a_o);
- Arithmetic::subtractProduct(rhs, tau / 2.0, A, v_o);
- Arithmetic::subtractProduct(rhs, A, u_o);
+ Arithmetic::addProduct(rhs, 2.0 / tau, matrices.mass, v_o);
+ Arithmetic::addProduct(rhs, matrices.mass, a_o);
+ Arithmetic::subtractProduct(rhs, tau / 2.0, matrices.elasticity, v_o);
+ Arithmetic::subtractProduct(rhs, matrices.elasticity, u_o);
}
iterate = v_o;
diff --git a/src/timestepping/newmark.hh b/src/timestepping/newmark.hh
index 2ed4a1b4ffadc7416648196496f9f175fe8d244c..72ae85433a0f7bdbd4539317ece85e6e8d87fb73 100644
--- a/src/timestepping/newmark.hh
+++ b/src/timestepping/newmark.hh
@@ -4,9 +4,8 @@
template <class Vector, class Matrix, class Function, size_t dim>
class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
public:
- Newmark(Matrix const &_A, Matrix const &_M, Matrix const &_C,
- Vector const &_u_initial, Vector const &_v_initial,
- Vector const &_a_initial,
+ Newmark(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
+ Vector const &_v_initial, Vector const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction);
@@ -22,9 +21,7 @@ class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
const;
private:
- Matrix const &A;
- Matrix const &M;
- Matrix const &C;
+ Matrices<Matrix> const &matrices;
Vector u;
Vector v;
Vector a;