diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 0ea095a87d9397fe33ea00f5a8701147d30c555d..24a17b2703361a99b503009a18972eff7a9c3bfa 100644
--- a/src/one-body-sample.cc
+++ b/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);
diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset
index 5bd61810f6359d84f54723fbd99a2acbf8986c43..75610beee33a5e183b50c86e6826acf22a90681f 100644
--- a/src/one-body-sample.parset
+++ b/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
diff --git a/src/timestepping/backward_euler.cc b/src/timestepping/backward_euler.cc
index eef042e35e52d896009533e2e80ea74518ad7d06..684af45838ee08b2571c3829c8335702aa63c9e4 100644
--- a/src/timestepping/backward_euler.cc
+++ b/src/timestepping/backward_euler.cc
@@ -1,11 +1,13 @@
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
diff --git a/src/timestepping/backward_euler.hh b/src/timestepping/backward_euler.hh
index 98eb362b92d994f697db52f54353c8e8a9fa9052..02e8ddc62255b6d5a928deed22ab3efeae6eda53 100644
--- a/src/timestepping/backward_euler.hh
+++ b/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;
diff --git a/src/timestepping/newmark.cc b/src/timestepping/newmark.cc
index f2190ac452f6c888a38dd7c7fe32c985c53a8da9..bd9dcbc0ae5a446108c64f1cf2596cedade188b9 100644
--- a/src/timestepping/newmark.cc
+++ b/src/timestepping/newmark.cc
@@ -1,11 +1,14 @@
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
diff --git a/src/timestepping/newmark.hh b/src/timestepping/newmark.hh
index f4dc7e600b4314b58ef401d3c28b05abea236a0a..cfd3b5f1744f85d2086852eba29bba827cdd2129 100644
--- a/src/timestepping/newmark.hh
+++ b/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;