diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 5571ac5c3230395a870e3642fdfc094915af861b..4747485f483853b2dbb6194169914b44cbf59add 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -56,6 +56,7 @@
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/viscosityassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/dunepython.hh>
@@ -106,19 +107,18 @@ 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) {
+ MatrixType const &dampingMatrix, 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, dampingMatrix, wc, u_initial, v_initial,
+ stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
a_initial, velocityDirichletNodes, velocityDirichletFunction);
case Config::BackwardEuler:
return Dune::make_shared<
BackwardEuler<VectorType, MatrixType, FunctionType, dims>>(
- stiffnessMatrix, massMatrix, dampingMatrix, wc, u_initial, v_initial,
+ stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
velocityDirichletNodes, velocityDirichletFunction);
default:
assert(false);
@@ -312,9 +312,14 @@ int main(int argc, char *argv[]) {
localStiffness(E, nu);
p1Assembler.assembleOperator(localStiffness, A);
}
-
- auto const wc = parset.get<double>("problem.damping");
- MatrixType const &C = A;
+ MatrixType C;
+ {
+ ViscosityAssembler<GridType, P1Basis::LocalFiniteElement,
+ P1Basis::LocalFiniteElement> const
+ localViscosity(parset.get<double>("body.shearViscosity"),
+ parset.get<double>("body.bulkViscosity"));
+ p1Assembler.assembleOperator(localViscosity, C);
+ }
SingletonMatrixType frictionalBoundaryMassMatrix;
EnergyNorm<SingletonMatrixType, SingletonVectorType> const stateEnergyNorm(
@@ -404,13 +409,12 @@ int main(int argc, char *argv[]) {
{
accelerationRHS = 0.0;
Arithmetic::addProduct(accelerationRHS, A, u_initial);
- Arithmetic::addProduct(accelerationRHS, wc, C, v_initial);
+ Arithmetic::addProduct(accelerationRHS, C, v_initial);
myGlobalNonlinearity->updateState(alpha_initial);
// NOTE: We assume differentiability of Psi at v0 here!
myGlobalNonlinearity->addGradient(v_initial, accelerationRHS);
accelerationRHS -= ell;
- // instead of multiplying M by (1.0 - wc), we divide the RHS
- accelerationRHS *= -1.0 / (1.0 - wc);
+ accelerationRHS *= -1.0;
}
LinearFactoryType accelerationFactory(parset.sub("solver.tnnmg"), // FIXME
refinements, *grid,
@@ -460,7 +464,7 @@ int main(int argc, char *argv[]) {
auto timeSteppingScheme =
initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunction, velocityDirichletNodes, M, A,
- C, wc, u_initial, v_initial, a_initial);
+ C, 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 7eb60fab0c0a6441526dc04134c0e832d0fd40e9..4f663a2fafd60496faa6e117456c4272e5589eb6 100644
--- a/src/one-body-sample.parset
+++ b/src/one-body-sample.parset
@@ -5,7 +5,6 @@ writeVTK = false
[problem]
finalTime = 15
-damping = 0.0 # Needs to lie in [0,1)
[body]
E = 5e7
@@ -13,6 +12,8 @@ nu = 0.3 # The closer we get to 0.5, the more wiggly everything gets
density = 5000
height = 1
width = 5
+shearViscosity = 0
+bulkViscosity = 0
[boundary.friction]
mu0 = 0.6
diff --git a/src/timestepping/backward_euler.cc b/src/timestepping/backward_euler.cc
index d9828616cd137f2e9274744d55ade3a87a80de24..288ddb6094f658c54752f90fbf6f322511c74ee3 100644
--- a/src/timestepping/backward_euler.cc
+++ b/src/timestepping/backward_euler.cc
@@ -1,13 +1,12 @@
template <class VectorType, class MatrixType, class FunctionType, size_t dim>
BackwardEuler<VectorType, MatrixType, FunctionType, dim>::BackwardEuler(
MatrixType const &_A, MatrixType const &_M, MatrixType const &_C,
- double _wc, VectorType const &_u_initial, VectorType const &_v_initial,
+ 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),
@@ -57,14 +56,14 @@ void BackwardEuler<VectorType, MatrixType, FunctionType, dim>::setup(
indices.exportIdx(AM);
}
AM = 0.0;
- Arithmetic::addProduct(AM, (1.0 - wc) / tau, M);
- Arithmetic::addProduct(AM, wc, C);
+ Arithmetic::addProduct(AM, 1.0 / tau, M);
+ Arithmetic::addProduct(AM, 1.0, C);
Arithmetic::addProduct(AM, tau, A);
// set up RHS
{
rhs = ell;
- Arithmetic::addProduct(rhs, (1.0 - wc) / tau, M, v_o);
+ Arithmetic::addProduct(rhs, 1.0 / tau, M, v_o);
Arithmetic::subtractProduct(rhs, A, u_o);
}
diff --git a/src/timestepping/backward_euler.hh b/src/timestepping/backward_euler.hh
index 02e8ddc62255b6d5a928deed22ab3efeae6eda53..50d73d0a101b419d2856c9f32baba86f8558ba9f 100644
--- a/src/timestepping/backward_euler.hh
+++ b/src/timestepping/backward_euler.hh
@@ -6,7 +6,7 @@ class BackwardEuler
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
BackwardEuler(MatrixType const &_A, MatrixType const &_M,
- MatrixType const &_C, double _wc, VectorType const &_u_initial,
+ MatrixType const &_C, VectorType const &_u_initial,
VectorType const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction);
@@ -22,7 +22,6 @@ class BackwardEuler
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 9e7fb91f9ba2baeceb4173285934ad7e05609c54..d1ce9f3962aeaacb29cfe5532aeea3bac28ab63b 100644
--- a/src/timestepping/newmark.cc
+++ b/src/timestepping/newmark.cc
@@ -1,14 +1,13 @@
template <class VectorType, class MatrixType, class FunctionType, size_t dim>
Newmark<VectorType, MatrixType, FunctionType, dim>::Newmark(
MatrixType const &_A, MatrixType const &_M, MatrixType const &_C,
- double _wc, VectorType const &_u_initial, VectorType const &_v_initial,
+ 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),
@@ -61,15 +60,15 @@ void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
indices.exportIdx(AM);
}
AM = 0.0;
- Arithmetic::addProduct(AM, (1.0 - wc) * 2.0 / tau, M);
- Arithmetic::addProduct(AM, wc, C);
+ Arithmetic::addProduct(AM, 2.0 / tau, M);
+ Arithmetic::addProduct(AM, 1.0, C);
Arithmetic::addProduct(AM, tau / 2.0, A);
// set up RHS
{
rhs = ell;
- Arithmetic::addProduct(rhs, (1.0 - wc) * 2.0 / tau, M, v_o);
- Arithmetic::addProduct(rhs, 1.0 - wc, M, a_o);
+ 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);
}
diff --git a/src/timestepping/newmark.hh b/src/timestepping/newmark.hh
index cfd3b5f1744f85d2086852eba29bba827cdd2129..ce00de305a579a4fc36256e28e5cb861f9210171 100644
--- a/src/timestepping/newmark.hh
+++ b/src/timestepping/newmark.hh
@@ -6,8 +6,8 @@ class Newmark
: public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
public:
Newmark(MatrixType const &_A, MatrixType const &_M, MatrixType const &_C,
- double _wc, VectorType const &_u_initial,
- VectorType const &_v_initial, VectorType const &_a_initial,
+ VectorType const &_u_initial, VectorType const &_v_initial,
+ VectorType const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
FunctionType const &_dirichletFunction);
@@ -22,7 +22,6 @@ class Newmark
MatrixType const &A;
MatrixType const &M;
MatrixType const &C;
- double wc;
VectorType u;
VectorType v;
VectorType a;