diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 173aa45054c36f46c97102da3aabdffc9cb0a760..264eb169b8ab712e7545e19e5d76e7522ec50bee 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -234,6 +234,8 @@ int main(int argc, char *argv[]) {
u_initial = 0.0;
solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
parset.sub("u0.solver"));
+ Vector ur_initial(fineVertexCount);
+ ur_initial = 0.0;
ScalarVector alpha_initial(fineVertexCount);
alpha_initial = parset.get<double>("boundary.friction.initialLogState");
@@ -253,6 +255,8 @@ int main(int argc, char *argv[]) {
velocityDirichletFunction.evaluate(0.0, v_initial_const);
assert(v_initial_const == 0.0);
}
+ Vector vr_initial(fineVertexCount);
+ vr_initial = 0.0;
Vector a_initial(fineVertexCount);
a_initial = 0.0;
@@ -304,7 +308,7 @@ int main(int argc, char *argv[]) {
frictionWriter.writeKinetics(_u, _v);
frictionWriter.writeOther(c, _alpha);
};
- report(u_initial, v_initial, alpha_initial);
+ report(ur_initial, vr_initial, alpha_initial);
MyVTKWriter<typename MyAssembler::VertexBasis,
typename MyAssembler::CellBasis> const
@@ -314,7 +318,7 @@ int main(int argc, char *argv[]) {
ScalarVector stress;
myAssembler.assembleVonMisesStress(
body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress);
- vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress);
+ vtkWriter.write(0, ur_initial, vr_initial, alpha_initial, stress);
}
SpecialWriter<GridView, dims> specialVelocityWriter("specialVelocities",
@@ -334,10 +338,10 @@ int main(int argc, char *argv[]) {
std::fstream iterationWriter("iterations", std::fstream::out),
relaxationWriter("relaxation", std::fstream::out);
- auto timeSteppingScheme =
- initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
- velocityDirichletFunction, dirichletNodes, M, A, C,
- u_initial, v_initial, a_initial);
+ auto timeSteppingScheme = initTimeStepper(
+ parset.get<Config::scheme>("timeSteps.scheme"),
+ velocityDirichletFunction, dirichletNodes, M, A, C, u_initial,
+ ur_initial, v_initial, vr_initial, a_initial);
auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
alpha_initial, frictionalNodes,
@@ -451,24 +455,30 @@ int main(int argc, char *argv[]) {
if (printProgress)
std::cout << std::endl;
- report(u, v, alpha);
+ Vector ur;
+ Vector vr;
+ timeSteppingScheme->postProcessRelativeQuantities();
+ timeSteppingScheme->extractRelativeDisplacement(ur);
+ timeSteppingScheme->extractRelativeVelocity(vr);
+
+ report(ur, vr, alpha);
iterationWriter << std::endl;
relaxationWriter << std::endl;
{
- BasisGridFunction<typename MyAssembler::VertexBasis, Vector> velocity(
- myAssembler.vertexBasis, v);
BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
- displacement(myAssembler.vertexBasis, u);
- specialVelocityWriter.write(velocity);
- specialDisplacementWriter.write(displacement);
+ relativeVelocity(myAssembler.vertexBasis, vr);
+ BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
+ relativeDisplacement(myAssembler.vertexBasis, ur);
+ specialVelocityWriter.write(relativeVelocity);
+ specialDisplacementWriter.write(relativeDisplacement);
}
if (parset.get<bool>("io.writeVTK")) {
ScalarVector stress;
myAssembler.assembleVonMisesStress(body.getYoungModulus(),
body.getPoissonRatio(), u, stress);
- vtkWriter.write(timeStep, u, v, alpha, stress);
+ vtkWriter.write(timeStep, ur, vr, alpha, stress);
}
}
iterationWriter.close();
diff --git a/src/sliding-block.cc b/src/sliding-block.cc
index 42a284e62512f895783338c165d7445252fbef6f..ef52fbe3c59b4e356ebcc0a4b56268c7e6214326 100644
--- a/src/sliding-block.cc
+++ b/src/sliding-block.cc
@@ -222,6 +222,8 @@ int main(int argc, char *argv[]) {
u_initial = 0.0;
solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
parset.sub("u0.solver"));
+ Vector ur_initial(fineVertexCount);
+ ur_initial = 0.0;
ScalarVector alpha_initial(fineVertexCount);
alpha_initial =
@@ -246,6 +248,9 @@ int main(int argc, char *argv[]) {
for (auto &x : v_initial)
x[0] = v_initial_const;
}
+ Vector vr_initial(fineVertexCount);
+ vr_initial = 0.0;
+
Vector a_initial(fineVertexCount);
a_initial = 0.0;
{
@@ -285,7 +290,7 @@ int main(int argc, char *argv[]) {
frictionWriter.writeKinetics(_u, _v);
frictionWriter.writeOther(c, _alpha);
};
- reportFriction(u_initial, v_initial, alpha_initial);
+ reportFriction(ur_initial, vr_initial, alpha_initial);
MyVTKWriter<typename MyAssembler::VertexBasis,
typename MyAssembler::CellBasis> const
@@ -295,7 +300,7 @@ int main(int argc, char *argv[]) {
ScalarVector stress;
myAssembler.assembleVonMisesStress(
body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress);
- vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress);
+ vtkWriter.write(0, ur_initial, vr_initial, alpha_initial, stress);
}
// Set up TNNMG solver
@@ -310,10 +315,10 @@ int main(int argc, char *argv[]) {
std::fstream iterationWriter("iterations", std::fstream::out),
relaxationWriter("relaxation", std::fstream::out);
- auto timeSteppingScheme =
- initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
- velocityDirichletFunction, dirichletNodes, M, A, C,
- u_initial, v_initial, a_initial);
+ auto timeSteppingScheme = initTimeStepper(
+ parset.get<Config::scheme>("timeSteps.scheme"),
+ velocityDirichletFunction, dirichletNodes, M, A, C, u_initial,
+ ur_initial, v_initial, vr_initial, a_initial);
auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
alpha_initial, frictionalNodes,
@@ -427,7 +432,13 @@ int main(int argc, char *argv[]) {
if (printProgress)
std::cout << std::endl;
- reportFriction(u, v, alpha);
+ Vector ur;
+ Vector vr;
+ timeSteppingScheme->postProcessRelativeQuantities();
+ timeSteppingScheme->extractRelativeDisplacement(ur);
+ timeSteppingScheme->extractRelativeVelocity(vr);
+
+ reportFriction(ur, vr, alpha);
iterationWriter << std::endl;
relaxationWriter << std::endl;
@@ -435,7 +446,7 @@ int main(int argc, char *argv[]) {
ScalarVector stress;
myAssembler.assembleVonMisesStress(body.getYoungModulus(),
body.getPoissonRatio(), u, stress);
- vtkWriter.write(timeStep, u, v, alpha, stress);
+ vtkWriter.write(timeStep, ur, vr, alpha, stress);
}
}
iterationWriter.close();
diff --git a/src/timestepping.hh b/src/timestepping.hh
index f2ee1f8396389051eee223a23c7c75e72c186473..258327637a1057e1132d5bbcc07d66ff69d298d5 100644
--- a/src/timestepping.hh
+++ b/src/timestepping.hh
@@ -12,8 +12,11 @@ class TimeSteppingScheme {
Vector &rhs, Vector &iterate, Matrix &AB) = 0;
void virtual postProcess(Vector const &iterate) = 0;
+ void virtual postProcessRelativeQuantities() = 0;
void virtual extractDisplacement(Vector &displacement) const = 0;
+ void virtual extractRelativeDisplacement(Vector &displacement) const = 0;
void virtual extractVelocity(Vector &velocity) const = 0;
+ void virtual extractRelativeVelocity(Vector &velocity) const = 0;
void virtual extractOldVelocity(Vector &velocity) const = 0;
};
@@ -27,17 +30,19 @@ initTimeStepper(Config::scheme scheme,
Dune::BitSetVector<dimension> const &velocityDirichletNodes,
Matrix const &massMatrix, Matrix const &stiffnessMatrix,
Matrix const &dampingMatrix, Vector const &u_initial,
- Vector const &v_initial, Vector const &a_initial) {
+ Vector const &ur_initial, Vector const &v_initial,
+ Vector const &vr_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);
+ stiffnessMatrix, massMatrix, dampingMatrix, u_initial, ur_initial,
+ v_initial, vr_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);
+ stiffnessMatrix, massMatrix, dampingMatrix, u_initial, ur_initial,
+ v_initial, velocityDirichletNodes, velocityDirichletFunction);
default:
assert(false);
}
diff --git a/src/timestepping/backward_euler.cc b/src/timestepping/backward_euler.cc
index 6d5f6bed7611082acca1dd2942870341a52974cb..442897bd4064c53726dc5b8491dcc4da7f9732dd 100644
--- a/src/timestepping/backward_euler.cc
+++ b/src/timestepping/backward_euler.cc
@@ -1,13 +1,14 @@
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,
+ Vector const &_u_initial, Vector const &_ur_initial,
+ Vector const &_v_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
: A(_A),
M(_M),
C(_C),
u(_u_initial),
+ ur(_ur_initial),
v(_v_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
@@ -16,6 +17,7 @@ template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::nextTimeStep() {
v_o = v;
u_o = u;
+ ur_o = ur;
}
template <class Vector, class Matrix, class Function, size_t dim>
@@ -23,6 +25,7 @@ void BackwardEuler<Vector, Matrix, Function, dim>::setup(
Vector const &ell, double _tau, double relativeTime, Vector &rhs,
Vector &iterate, Matrix &AM) {
postProcessCalled = false;
+ postProcessRelativeQuantitiesCalled = false;
dirichletFunction.evaluate(relativeTime, dirichletValue);
tau = _tau;
@@ -87,6 +90,21 @@ void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
Arithmetic::addProduct(u, tau, v);
}
+template <class Vector, class Matrix, class Function, size_t dim>
+void
+BackwardEuler<Vector, Matrix, Function, dim>::postProcessRelativeQuantities() {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+ postProcessRelativeQuantitiesCalled = true;
+
+ vr = v;
+ for (size_t i = 0; i < dirichletNodes.size(); ++i)
+ vr[i][0] -= dirichletValue;
+
+ ur = ur_o;
+ Arithmetic::addProduct(ur, tau, vr);
+}
+
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractDisplacement(
Vector &displacement) const {
@@ -96,6 +114,16 @@ void BackwardEuler<Vector, Matrix, Function, dim>::extractDisplacement(
displacement = u;
}
+template <class Vector, class Matrix, class Function, size_t dim>
+void BackwardEuler<Vector, Matrix, Function, dim>::extractRelativeDisplacement(
+ Vector &displacement) const {
+ if (!postProcessRelativeQuantitiesCalled)
+ DUNE_THROW(Dune::Exception,
+ "It seems you forgot to call postProcessRelativeQuantities!");
+
+ displacement = u;
+}
+
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractVelocity(
Vector &velocity) const {
@@ -105,6 +133,16 @@ void BackwardEuler<Vector, Matrix, Function, dim>::extractVelocity(
velocity = v;
}
+template <class Vector, class Matrix, class Function, size_t dim>
+void BackwardEuler<Vector, Matrix, Function, dim>::extractRelativeVelocity(
+ Vector &velocity) const {
+ if (!postProcessRelativeQuantitiesCalled)
+ DUNE_THROW(Dune::Exception,
+ "It seems you forgot to call postProcessRelativeQuantities!");
+
+ velocity = vr;
+}
+
template <class Vector, class MatrixType, class FunctionType, size_t dim>
void BackwardEuler<Vector, MatrixType, FunctionType, dim>::extractOldVelocity(
Vector &velocity) const {
diff --git a/src/timestepping/backward_euler.hh b/src/timestepping/backward_euler.hh
index a07b878d8f9cea2108406bf44f0455c261244f3c..657b57b5669bd2b0c599454ccd9c901bc0999be0 100644
--- a/src/timestepping/backward_euler.hh
+++ b/src/timestepping/backward_euler.hh
@@ -5,7 +5,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,
+ Vector const &_u_initial, Vector const &_ur_initial,
+ Vector const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction);
@@ -13,8 +14,11 @@ class BackwardEuler : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
void setup(Vector const &, double, double, Vector &, Vector &,
Matrix &) override;
void postProcess(Vector const &) override;
+ void postProcessRelativeQuantities() override;
void extractDisplacement(Vector &) const override;
+ void extractRelativeDisplacement(Vector &) const override;
void extractVelocity(Vector &) const override;
+ void extractRelativeVelocity(Vector &) const override;
void extractOldVelocity(Vector &) const override;
private:
@@ -22,16 +26,20 @@ class BackwardEuler : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
Matrix const &M;
Matrix const &C;
Vector u;
+ Vector ur;
Vector v;
+ Vector vr;
Dune::BitSetVector<dim> const &dirichletNodes;
Function const &dirichletFunction;
double dirichletValue;
Vector u_o;
+ Vector ur_o;
Vector v_o;
double tau;
bool postProcessCalled = false;
+ bool postProcessRelativeQuantitiesCalled = false;
};
#endif
diff --git a/src/timestepping/newmark.cc b/src/timestepping/newmark.cc
index 6b70313585fc53bc4625310f5963db837524828e..5ec3c8a09f8af774ca6e3013f2dfe1197d0185b7 100644
--- a/src/timestepping/newmark.cc
+++ b/src/timestepping/newmark.cc
@@ -1,14 +1,17 @@
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 &_u_initial, Vector const &_ur_initial,
+ Vector const &_v_initial, Vector const &_vr_initial,
Vector const &_a_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
: A(_A),
M(_M),
C(_C),
u(_u_initial),
+ ur(_ur_initial),
v(_v_initial),
+ vr(_vr_initial),
a(_a_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
@@ -17,7 +20,9 @@ template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::nextTimeStep() {
a_o = a;
v_o = v;
+ vr_o = vr;
u_o = u;
+ ur_o = ur;
}
template <class Vector, class Matrix, class Function, size_t dim>
@@ -27,6 +32,7 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
Vector &rhs, Vector &iterate,
Matrix &AM) {
postProcessCalled = false;
+ postProcessRelativeQuantitiesCalled = false;
dirichletFunction.evaluate(relativeTime, dirichletValue);
tau = _tau;
@@ -102,6 +108,22 @@ void Newmark<Vector, Matrix, Function, dim>::postProcess(
Arithmetic::subtractProduct(a, 1.0, a_o);
}
+template <class Vector, class Matrix, class Function, size_t dim>
+void Newmark<Vector, Matrix, Function, dim>::postProcessRelativeQuantities() {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+ postProcessRelativeQuantitiesCalled = true;
+
+ vr = v;
+ for (size_t i = 0; i < dirichletNodes.size(); ++i)
+ vr[i][0] -= dirichletValue;
+
+ // u1 = tau/2 ( v1 + v0 ) + u0
+ ur = ur_o;
+ Arithmetic::addProduct(ur, tau / 2.0, vr);
+ Arithmetic::addProduct(ur, tau / 2.0, vr_o);
+}
+
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::extractDisplacement(
Vector &displacement) const {
@@ -111,6 +133,16 @@ void Newmark<Vector, Matrix, Function, dim>::extractDisplacement(
displacement = u;
}
+template <class Vector, class Matrix, class Function, size_t dim>
+void Newmark<Vector, Matrix, Function, dim>::extractRelativeDisplacement(
+ Vector &displacement) const {
+ if (!postProcessRelativeQuantitiesCalled)
+ DUNE_THROW(Dune::Exception,
+ "It seems you forgot to call postProcessRelativeQuantities!");
+
+ displacement = ur;
+}
+
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::extractVelocity(Vector &velocity)
const {
@@ -120,6 +152,16 @@ void Newmark<Vector, Matrix, Function, dim>::extractVelocity(Vector &velocity)
velocity = v;
}
+template <class Vector, class Matrix, class Function, size_t dim>
+void Newmark<Vector, Matrix, Function, dim>::extractRelativeVelocity(
+ Vector &velocity) const {
+ if (!postProcessRelativeQuantitiesCalled)
+ DUNE_THROW(Dune::Exception,
+ "It seems you forgot to call postProcessRelativeQuantities!");
+
+ velocity = vr;
+}
+
template <class Vector, class MatrixType, class FunctionType, size_t dim>
void Newmark<Vector, MatrixType, FunctionType, dim>::extractOldVelocity(
Vector &velocity) const {
diff --git a/src/timestepping/newmark.hh b/src/timestepping/newmark.hh
index 78f7731a19d22c7b8eed853de827ea17ca3f0d67..a9a6e1b2e877baa03bd3cc9b134c7f37226dfc4c 100644
--- a/src/timestepping/newmark.hh
+++ b/src/timestepping/newmark.hh
@@ -5,7 +5,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 &_u_initial, Vector const &_ur_initial,
+ Vector const &_v_initial, Vector const &_vr_initial,
Vector const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction);
@@ -14,8 +15,11 @@ class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
void setup(Vector const &, double, double, Vector &, Vector &,
Matrix &) override;
void postProcess(Vector const &) override;
+ void postProcessRelativeQuantities() override;
void extractDisplacement(Vector &) const override;
+ void extractRelativeDisplacement(Vector &) const override;
void extractVelocity(Vector &) const override;
+ void extractRelativeVelocity(Vector &) const override;
void extractOldVelocity(Vector &) const override;
private:
@@ -23,18 +27,23 @@ class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
Matrix const &M;
Matrix const &C;
Vector u;
+ Vector ur;
Vector v;
+ Vector vr;
Vector a;
Dune::BitSetVector<dim> const &dirichletNodes;
Function const &dirichletFunction;
double dirichletValue;
Vector u_o;
+ Vector ur_o;
Vector v_o;
+ Vector vr_o;
Vector a_o;
double tau;
bool postProcessCalled = false;
+ bool postProcessRelativeQuantitiesCalled = false;
};
#endif