From fa0910dd4e031aeb3dcd6d29c7a45dcc912f0bc9 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 14 Aug 2014 21:48:32 +0200
Subject: [PATCH] [Output ] Kill relative output

---
 src/fixedpointiterator.cc          |  1 -
 src/sand-wedge.cc                  | 30 ++++++++------------
 src/timestepping.hh                | 15 ++++------
 src/timestepping/backward_euler.cc | 42 ++--------------------------
 src/timestepping/backward_euler.hh | 10 +------
 src/timestepping/newmark.cc        | 44 +-----------------------------
 src/timestepping/newmark.hh        | 11 +-------
 7 files changed, 22 insertions(+), 131 deletions(-)

diff --git a/src/fixedpointiterator.cc b/src/fixedpointiterator.cc
index c013d991..5b7a394f 100644
--- a/src/fixedpointiterator.cc
+++ b/src/fixedpointiterator.cc
@@ -72,7 +72,6 @@ int FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::run(
     DUNE_THROW(Dune::Exception, "FPI failed to converge");
 
   velocityUpdater->postProcess(velocityIterate);
-  velocityUpdater->postProcessRelativeQuantities();
 
   return fixedPointIteration;
 }
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index f53ecbc8..e7273b23 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -311,8 +311,6 @@ int main(int argc, char *argv[]) {
     u_initial = 0.0;
     solveLinearProblem(dirichletNodes, A, ell0, u_initial, ANorm,
                        parset.sub("u0.solver"));
-    Vector ur_initial(leafVertexCount);
-    ur_initial = 0.0;
 
     ScalarVector alpha_initial(leafVertexCount);
     alpha_initial = parset.get<double>("boundary.friction.initialAlpha");
@@ -333,8 +331,6 @@ int main(int argc, char *argv[]) {
       velocityDirichletFunction.evaluate(0.0, v_initial_const);
       assert(std::abs(v_initial_const) < 1e-14);
     }
-    Vector vr_initial(leafVertexCount);
-    vr_initial = 0.0;
 
     Vector a_initial(leafVertexCount);
     a_initial = 0.0;
@@ -386,7 +382,7 @@ int main(int argc, char *argv[]) {
       frictionWriter.writeKinetics(_u, _v);
       frictionWriter.writeOther(c, _alpha);
     };
-    report(ur_initial, vr_initial, alpha_initial);
+    report(u_initial, v_initial, alpha_initial);
 
     MyVTKWriter<typename MyAssembler::VertexBasis,
                 typename MyAssembler::CellBasis> const
@@ -396,7 +392,7 @@ int main(int argc, char *argv[]) {
       ScalarVector stress;
       myAssembler.assembleVonMisesStress(
           body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress);
-      vtkWriter.write(0, ur_initial, vr_initial, alpha_initial, stress);
+      vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress);
     }
 
     SpecialWriter<GridView, dims> specialVelocityWriter("specialVelocities",
@@ -423,8 +419,7 @@ int main(int argc, char *argv[]) {
             parset.get<double>("boundary.friction.V0")),
         initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
                         velocityDirichletFunction, dirichletNodes, M, A, C,
-                        u_initial, ur_initial, v_initial, vr_initial,
-                        a_initial));
+                        u_initial, v_initial, a_initial));
 
     auto const refinementTolerance =
         parset.get<double>("timeSteps.refinementTolerance");
@@ -449,28 +444,27 @@ int main(int argc, char *argv[]) {
       reportTimeStep(adaptiveTimeStepper.getRelativeTime(),
                      adaptiveTimeStepper.getRelativeTau());
 
-      Vector u, ur, vr;
+      Vector u, v;
       ScalarVector alpha;
       current.second->extractDisplacement(u);
-      current.second->extractRelativeDisplacement(ur);
-      current.second->extractRelativeVelocity(vr);
+      current.second->extractVelocity(v);
       current.first->extractAlpha(alpha);
 
-      report(ur, vr, alpha);
+      report(u, v, alpha);
       {
+        BasisGridFunction<typename MyAssembler::VertexBasis, Vector> velocity(
+            myAssembler.vertexBasis, v);
         BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
-        relativeVelocity(myAssembler.vertexBasis, vr);
-        BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
-        relativeDisplacement(myAssembler.vertexBasis, ur);
-        specialVelocityWriter.write(relativeVelocity);
-        specialDisplacementWriter.write(relativeDisplacement);
+        displacement(myAssembler.vertexBasis, u);
+        specialVelocityWriter.write(velocity);
+        specialDisplacementWriter.write(displacement);
       }
 
       if (parset.get<bool>("io.writeVTK")) {
         ScalarVector stress;
         myAssembler.assembleVonMisesStress(body.getYoungModulus(),
                                            body.getPoissonRatio(), u, stress);
-        vtkWriter.write(timeStep, ur, vr, alpha, stress);
+        vtkWriter.write(timeStep, u, v, alpha, stress);
       }
       timeStep++;
     }
diff --git a/src/timestepping.hh b/src/timestepping.hh
index 1221c585..3d54198b 100644
--- a/src/timestepping.hh
+++ b/src/timestepping.hh
@@ -14,11 +14,8 @@ 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;
 
   std::shared_ptr<
@@ -36,19 +33,17 @@ initTimeStepper(Config::scheme scheme,
                 Dune::BitSetVector<dimension> const &velocityDirichletNodes,
                 Matrix const &massMatrix, Matrix const &stiffnessMatrix,
                 Matrix const &dampingMatrix, Vector const &u_initial,
-                Vector const &ur_initial, Vector const &v_initial,
-                Vector const &vr_initial, Vector const &a_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, ur_initial,
-          v_initial, vr_initial, a_initial, velocityDirichletNodes,
-          velocityDirichletFunction);
+          stiffnessMatrix, massMatrix, dampingMatrix, 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, ur_initial,
-          v_initial, velocityDirichletNodes, velocityDirichletFunction);
+          stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
+          velocityDirichletNodes, velocityDirichletFunction);
     default:
       assert(false);
   }
diff --git a/src/timestepping/backward_euler.cc b/src/timestepping/backward_euler.cc
index 6f810d59..73ecfbd3 100644
--- a/src/timestepping/backward_euler.cc
+++ b/src/timestepping/backward_euler.cc
@@ -3,14 +3,13 @@
 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 &_ur_initial,
-    Vector const &_v_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
+    Vector const &_u_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) {}
@@ -19,7 +18,6 @@ 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>
@@ -27,7 +25,6 @@ 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;
 
@@ -91,21 +88,6 @@ 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 {
@@ -115,16 +97,6 @@ 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 = ur;
-}
-
 template <class Vector, class Matrix, class Function, size_t dim>
 void BackwardEuler<Vector, Matrix, Function, dim>::extractVelocity(
     Vector &velocity) const {
@@ -134,16 +106,6 @@ 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 Matrix, class Function, size_t dim>
 void BackwardEuler<Vector, Matrix, Function, dim>::extractOldVelocity(
     Vector &velocity) const {
diff --git a/src/timestepping/backward_euler.hh b/src/timestepping/backward_euler.hh
index 510a149f..e6c0101d 100644
--- a/src/timestepping/backward_euler.hh
+++ b/src/timestepping/backward_euler.hh
@@ -5,8 +5,7 @@ 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 &_ur_initial,
-                Vector const &_v_initial,
+                Vector const &_u_initial, Vector const &_v_initial,
                 Dune::BitSetVector<dim> const &_dirichletNodes,
                 Function const &_dirichletFunction);
 
@@ -14,11 +13,8 @@ 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;
 
   std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>> clone()
@@ -29,20 +25,16 @@ 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 12e39abc..bb546a1f 100644
--- a/src/timestepping/newmark.cc
+++ b/src/timestepping/newmark.cc
@@ -3,17 +3,14 @@
 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 &_ur_initial,
-    Vector const &_v_initial, Vector const &_vr_initial,
+    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),
       u(_u_initial),
-      ur(_ur_initial),
       v(_v_initial),
-      vr(_vr_initial),
       a(_a_initial),
       dirichletNodes(_dirichletNodes),
       dirichletFunction(_dirichletFunction) {}
@@ -22,9 +19,7 @@ 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>
@@ -34,7 +29,6 @@ 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;
 
@@ -109,22 +103,6 @@ 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 {
@@ -134,16 +112,6 @@ 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 {
@@ -153,16 +121,6 @@ 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 Matrix, class Function, size_t dim>
 void Newmark<Vector, Matrix, Function, dim>::extractOldVelocity(
     Vector &velocity) const {
diff --git a/src/timestepping/newmark.hh b/src/timestepping/newmark.hh
index 0d711447..89b257c0 100644
--- a/src/timestepping/newmark.hh
+++ b/src/timestepping/newmark.hh
@@ -5,8 +5,7 @@ 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 &_ur_initial,
-          Vector const &_v_initial, Vector const &_vr_initial,
+          Vector const &_u_initial, Vector const &_v_initial,
           Vector const &_a_initial,
           Dune::BitSetVector<dim> const &_dirichletNodes,
           Function const &_dirichletFunction);
@@ -15,11 +14,8 @@ 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;
 
   std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>> clone()
@@ -30,23 +26,18 @@ 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
-- 
GitLab