[Cleanup] Introduce Matrices class

src/matrices.hh

+#ifndef SRC_MATRICES_HH
+#define SRC_MATRICES_HH
+
+template <class Matrix> struct Matrices {
+  Matrix elasticity;
+  Matrix damping;
+  Matrix mass;
+};
+#endif

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 =

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);
   }

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;

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;

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;

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;