[Cleanup] Make wc into a parameter

src/one-body-sample.cc

@@ -100,18 +100,19 @@ initTimeStepper(Config::scheme scheme,
                 FunctionType const &velocityDirichletFunction,
                 Dune::BitSetVector<dimension> const &velocityDirichletNodes,
                 MatrixType const &massMatrix, MatrixType const &stiffnessMatrix,
+                MatrixType const &dampingMatrix, double wc,
                 VectorType const &u_initial, VectorType const &v_initial,
                 VectorType const &a_initial) {
   switch (scheme) {
     case Config::Newmark:
       return Dune::make_shared<
           Newmark<VectorType, MatrixType, FunctionType, dims>>(
-          stiffnessMatrix, massMatrix, u_initial, v_initial, a_initial,
-          velocityDirichletNodes, velocityDirichletFunction);
+          stiffnessMatrix, massMatrix, dampingMatrix, wc, u_initial, v_initial,
+          a_initial, velocityDirichletNodes, velocityDirichletFunction);
     case Config::BackwardEuler:
       return Dune::make_shared<
           BackwardEuler<VectorType, MatrixType, FunctionType, dims>>(
-          stiffnessMatrix, massMatrix, u_initial, v_initial,
+          stiffnessMatrix, massMatrix, dampingMatrix, wc, u_initial, v_initial,
           velocityDirichletNodes, velocityDirichletFunction);
     default:
       assert(false);
@@ -314,6 +315,9 @@ int main(int argc, char *argv[]) {
       p1Assembler.assembleOperator(localStiffness, A);
     }
 
+    auto const wc = parset.get<double>("problem.damping");
+    MatrixType const &C = A;
+
     SingletonMatrixType frictionalBoundaryMassMatrix;
     EnergyNorm<SingletonMatrixType, SingletonVectorType> const stateEnergyNorm(
         frictionalBoundaryMassMatrix);
@@ -398,8 +402,6 @@ int main(int argc, char *argv[]) {
       /* We solve Au + Cv + Ma + Psi(v) = ell, thus
                                      Ma = - (Au + Cv + Psi(v) - ell)
       */
-      MatrixType const &C = A;
-      double const wc = 0.0; // watch out -- do not choose 1.0 here!
       VectorType accelerationRHS(finestSize);
       {
         accelerationRHS = 0.0;
@@ -462,7 +464,7 @@ int main(int argc, char *argv[]) {
     auto timeSteppingScheme =
         initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
                         velocityDirichletFunction, velocityDirichletNodes, M, A,
-                        u_initial, v_initial, a_initial);
+                        C, wc, u_initial, v_initial, a_initial);
     auto stateUpdater = initStateUpdater<SingletonVectorType, VectorType>(
         parset.get<Config::stateModel>("boundary.friction.stateModel"),
         alpha_initial, frictionalNodes, frictionData);

src/one-body-sample.parset

@@ -6,6 +6,7 @@ writeVTK      = false
 
 [problem]
 finalTime = 15
+damping   = 0.0 # Needs to lie in [0,1)
 
 [body]
 E       = 5e7

src/timestepping/backward_euler.cc

 template <class VectorType, class MatrixType, class FunctionType, size_t dim>
 BackwardEuler<VectorType, MatrixType, FunctionType, dim>::BackwardEuler(
-    MatrixType const &_A, MatrixType const &_M, VectorType const &_u_initial,
-    VectorType const &_v_initial,
+    MatrixType const &_A, MatrixType const &_M, MatrixType const &_C,
+    double _wc, VectorType const &_u_initial, VectorType const &_v_initial,
     Dune::BitSetVector<dim> const &_dirichletNodes,
     FunctionType const &_dirichletFunction)
     : A(_A),
       M(_M),
+      C(_C),
+      wc(_wc),
       u(_u_initial),
       v(_v_initial),
       dirichletNodes(_dirichletNodes),
@@ -25,8 +27,6 @@ void BackwardEuler<VectorType, MatrixType, FunctionType, dim>::setup(
 
   tau = _tau;
 
-  MatrixType const &C = A;
-  double const wc = 0.0;
   /* We start out with the formulation
 
        M a + C v + A u = ell

src/timestepping/backward_euler.hh

@@ -6,7 +6,8 @@ class BackwardEuler
     : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
 public:
   BackwardEuler(MatrixType const &_A, MatrixType const &_M,
-                VectorType const &_u_initial, VectorType const &_v_initial,
+                MatrixType const &_C, double _wc, VectorType const &_u_initial,
+                VectorType const &_v_initial,
                 Dune::BitSetVector<dim> const &_dirichletNodes,
                 FunctionType const &_dirichletFunction);
 
@@ -20,6 +21,8 @@ class BackwardEuler
 private:
   MatrixType const &A;
   MatrixType const &M;
+  MatrixType const &C;
+  double wc;
   VectorType u;
   VectorType v;
   Dune::BitSetVector<dim> const &dirichletNodes;

src/timestepping/newmark.cc

 template <class VectorType, class MatrixType, class FunctionType, size_t dim>
 Newmark<VectorType, MatrixType, FunctionType, dim>::Newmark(
-    MatrixType const &_A, MatrixType const &_M, VectorType const &_u_initial,
-    VectorType const &_v_initial, VectorType const &_a_initial,
+    MatrixType const &_A, MatrixType const &_M, MatrixType const &_C,
+    double _wc, VectorType const &_u_initial, VectorType const &_v_initial,
+    VectorType const &_a_initial,
     Dune::BitSetVector<dim> const &_dirichletNodes,
     FunctionType const &_dirichletFunction)
     : A(_A),
       M(_M),
+      C(_C),
+      wc(_wc),
       u(_u_initial),
       v(_v_initial),
       a(_a_initial),
@@ -27,8 +30,6 @@ void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
 
   tau = _tau;
 
-  MatrixType const &C = A;
-  double const wc = 0.0;
   /* We start out with the formulation
 
        M a + C v + A u = ell

src/timestepping/newmark.hh

@@ -5,9 +5,9 @@ template <class VectorType, class MatrixType, class FunctionType, size_t dim>
 class Newmark
     : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
 public:
-  Newmark(MatrixType const &_A, MatrixType const &_M,
-          VectorType const &_u_initial, VectorType const &_v_initial,
-          VectorType const &_a_initial,
+  Newmark(MatrixType const &_A, MatrixType const &_M, MatrixType const &_C,
+          double _wc, VectorType const &_u_initial,
+          VectorType const &_v_initial, VectorType const &_a_initial,
           Dune::BitSetVector<dim> const &_dirichletNodes,
           FunctionType const &_dirichletFunction);
 
@@ -21,6 +21,8 @@ class Newmark
 private:
   MatrixType const &A;
   MatrixType const &M;
+  MatrixType const &C;
+  double wc;
   VectorType u;
   VectorType v;
   VectorType a;