Implement eulerPair scheme (implicitEuler for u,v)

src/enum_scheme.cc

@@ -11,6 +11,9 @@ template <> struct StringToEnum<Config::scheme> {
     if (s == "newmark")
       return Config::Newmark;
 
+    if (s == "eulerPair")
+      return Config::EulerPair;
+
     DUNE_THROW(Dune::Exception, "failed to parse enum");
   }
 };

src/enums.hh

@@ -9,7 +9,8 @@ struct Config {
   enum scheme {
     ImplicitTwoStep,
     ImplicitEuler,
-    Newmark
+    Newmark,
+    EulerPair
   };
 };
 #endif

src/one-body-sample.cc

@@ -108,6 +108,11 @@ initTimeStepper(Config::scheme scheme, FunctionType const &dirichletFunction,
           Newmark<VectorType, MatrixType, FunctionType, dims>>(
           stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
           ignoreNodes, dirichletFunction);
+    case Config::EulerPair:
+      return Dune::make_shared<
+          EulerPair<VectorType, MatrixType, FunctionType, dims>>(
+          stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
+          ignoreNodes, dirichletFunction);
   }
 }
 template <class SingletonVectorType, class VectorType>

src/timestepping.cc

@@ -337,5 +337,111 @@ TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
 Newmark<VectorType, MatrixType, FunctionType, dim>::clone() {
   return new Newmark<VectorType, MatrixType, FunctionType, dim>(*this);
 }
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+EulerPair<VectorType, MatrixType, FunctionType, dim>::EulerPair(
+    MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
+    VectorType const &_ud_initial, VectorType const &_udd_initial,
+    Dune::BitSetVector<dim> const &_dirichletNodes,
+    FunctionType const &_dirichletFunction)
+    : A(_A),
+      B(_B),
+      u(_u_initial),
+      ud(_ud_initial),
+      udd(_udd_initial),
+      dirichletNodes(_dirichletNodes),
+      dirichletFunction(_dirichletFunction) {}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
+  udd_old = udd;
+  ud_old = ud;
+  u_old = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::setup(
+    VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
+    VectorType &problem_iterate, MatrixType &problem_A) {
+  postProcessCalled = false;
+
+  tau = _tau;
+
+  problem_rhs = ell;
+  B.usmv(1.0 / tau, ud_old, problem_rhs);
+  A.usmv(-1.0, u_old, problem_rhs);
+
+  // For fixed tau, we'd only really have to do this once
+  Dune::MatrixIndexSet indices(A.N(), A.M());
+  indices.import(A);
+  indices.import(B);
+  indices.exportIdx(problem_A);
+  problem_A = 0.0;
+  Arithmetic::addProduct(problem_A, tau, A);
+  Arithmetic::addProduct(problem_A, 1.0 / tau, B);
+
+  // ud_old makes a good initial iterate; we could use anything, though
+  problem_iterate = ud_old;
+
+  for (size_t i = 0; i < dirichletNodes.size(); ++i)
+    switch (dirichletNodes[i].count()) {
+      case 0:
+        continue;
+      case dim:
+        problem_iterate[i] = 0;
+        dirichletFunction.evaluate(time, problem_iterate[i][0]);
+        break;
+      case 1:
+        if (dirichletNodes[i][0]) {
+          dirichletFunction.evaluate(time, problem_iterate[i][0]);
+          break;
+        }
+        if (dirichletNodes[i][1]) {
+          problem_iterate[i][1] = 0;
+          break;
+        }
+        assert(false);
+      default:
+        assert(false);
+    }
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::postProcess(
+    VectorType const &problem_iterate) {
+  postProcessCalled = true;
+
+  ud = problem_iterate;
+
+  u = u_old;
+  Arithmetic::addProduct(u, tau, ud);
+
+  udd = 0;
+  Arithmetic::addProduct(udd, 1.0 / tau, ud);
+  Arithmetic::addProduct(udd, -1.0 / tau, ud_old);
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::extractDisplacement(
+    VectorType &displacement) const {
+  if (!postProcessCalled)
+    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+  displacement = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
+    VectorType &velocity) const {
+  if (!postProcessCalled)
+    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+  velocity = ud;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+EulerPair<VectorType, MatrixType, FunctionType, dim>::clone() {
+  return new EulerPair<VectorType, MatrixType, FunctionType, dim>(*this);
+}
 
 #include "timestepping_tmpl.cc"

src/timestepping.hh

@@ -112,6 +112,43 @@ class Newmark
   virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
   clone();
 
+private:
+  MatrixType const &A;
+  MatrixType const &B;
+  VectorType u;
+  VectorType ud;
+  VectorType udd;
+  Dune::BitSetVector<dim> const &dirichletNodes;
+  FunctionType const &dirichletFunction;
+
+  VectorType u_old;
+  VectorType ud_old;
+  VectorType udd_old;
+
+  double tau;
+
+  bool postProcessCalled = false;
+};
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class EulerPair
+    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+  EulerPair(MatrixType const &_A, MatrixType const &_B,
+            VectorType const &_u_initial, VectorType const &_ud_initial,
+            VectorType const &_udd_initial,
+            Dune::BitSetVector<dim> const &_dirichletNodes,
+            FunctionType const &_dirichletFunction);
+
+  void virtual nextTimeStep();
+  void virtual setup(VectorType const &, double, double, VectorType &,
+                     VectorType &, MatrixType &);
+  void virtual postProcess(VectorType const &);
+  void virtual extractDisplacement(VectorType &) const;
+  void virtual extractVelocity(VectorType &) const;
+
+  virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+  clone();
+
 private:
   MatrixType const &A;
   MatrixType const &B;

src/timestepping_tmpl.cc

@@ -17,3 +17,4 @@ typedef Dune::VirtualFunction<double, double> FunctionType;
 template class ImplicitEuler<VectorType, MatrixType, FunctionType, DIM>;
 template class ImplicitTwoStep<VectorType, MatrixType, FunctionType, DIM>;
 template class Newmark<VectorType, MatrixType, FunctionType, DIM>;
+template class EulerPair<VectorType, MatrixType, FunctionType, DIM>;