Implement Newmark scheme

src/enum_scheme.cc

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

src/enums.hh

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

src/one-body-sample.cc

@@ -47,6 +47,7 @@
 #include <dune/istl/bvector.hh>
 
 #include <dune/fufem/assemblers/functionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/operatorassembler.hh>
@@ -165,6 +166,40 @@ int main(int argc, char *argv[]) {
     P0Basis const p0Basis(leafView);
     P1Basis const p1Basis(leafView);
 
+    MassAssembler<GridType, P1Basis::LocalFiniteElement,
+                  P1Basis::LocalFiniteElement> const localMass;
+    MatrixType massMatrix;
+    {
+      timer.reset();
+      OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
+          .assemble(localMass, massMatrix);
+
+      // We have volume*gravity*density/area = normalstress (V*g*rho/A =
+      // sigma_n)
+      double volume = 1.0;
+      for (int i = 0; i < dim; ++i)
+        volume *= (upperRight[i] - lowerLeft[i]);
+
+      double area = 1.0;
+      for (int i = 0; i < dim; ++i)
+        if (i != 1)
+          area *= (upperRight[i] - lowerLeft[i]);
+
+      double const gravity = 9.81;
+
+      double const normalStress =
+          parset.get<double>("boundary.friction.normalstress");
+
+      // rho    = sigma     * A       / V   / g
+      // kg/m^d = N/m^(d-1) * m^(d-1) / m^d / (N/kg)
+      double const density = normalStress * area / volume / gravity;
+
+      massMatrix *= density;
+      if (parset.get<bool>("enable_timer"))
+        std::cout << "Assembled mass matrix in " << timer.elapsed() << "s"
+                  << std::endl;
+    }
+
     // Assemble elastic force on the body
     StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
                                P1Basis::LocalFiniteElement> const
@@ -200,6 +235,13 @@ int main(int argc, char *argv[]) {
 
     VectorType ud(finestSize);
     ud = 0.0;
+    VectorType ud_old(finestSize);
+    ud_old = 0.0;
+
+    VectorType udd(finestSize);
+    udd = 0.0;
+    VectorType udd_old(finestSize);
+    udd_old = 0.0;
 
     SingletonVectorType alpha_old(finestSize);
     alpha_old = parset.get<double>("boundary.friction.state.initial");
@@ -280,6 +322,12 @@ int main(int argc, char *argv[]) {
                   ell, stiffnessMatrix, u_old, u_old_old_ptr, ignoreNodes,
                   dirichletFunction, time, tau);
               break;
+            case Config::Newmark:
+              timeSteppingScheme = Dune::make_shared<Newmark<
+                  VectorType, MatrixType, decltype(dirichletFunction), dim>>(
+                  ell, stiffnessMatrix, massMatrix, u_old, ud_old, udd_old,
+                  ignoreNodes, dirichletFunction, time, tau);
+              break;
           }
         }
         timeSteppingScheme->setup(problem_rhs, problem_iterate, problem_A);
@@ -307,6 +355,11 @@ int main(int argc, char *argv[]) {
           timeSteppingScheme->extractSolution(problem_iterate, u);
           timeSteppingScheme->extractVelocity(problem_iterate, ud);
 
+          udd = 0;
+          Arithmetic::addProduct(udd, 2.0 / tau, ud);
+          Arithmetic::addProduct(udd, -2.0 / tau, ud_old);
+          Arithmetic::addProduct(udd, -1.0, udd_old);
+
           // Update the state
           for (size_t i = 0; i < frictionalNodes.size(); ++i) {
             if (frictionalNodes[i][0]) {
@@ -400,6 +453,8 @@ int main(int argc, char *argv[]) {
       }
       u_old_old = u_old;
       u_old = u;
+      ud_old = ud;
+      udd_old = udd;
       alpha_old = alpha;
 
       // Compute von Mises stress and write everything to a file

src/timestepping.cc

+#include <dune/fufem/arithmetic.hh>
+
 template <class VectorType, class MatrixType, class FunctionType, int dim>
 class TimeSteppingScheme {
 public:
@@ -149,3 +151,61 @@ class ImplicitTwoStep
 private:
   VectorType const *u_old_old_ptr;
 };
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class Newmark
+    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+  Newmark(VectorType const &_ell, MatrixType const &_A, MatrixType const &_B,
+          VectorType const &_u_old, VectorType const &_ud_old,
+          VectorType const &_udd_old,
+          Dune::BitSetVector<dim> const &_dirichletNodes,
+          FunctionType const &_dirichletFunction, double _time, double _tau)
+      : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
+            _ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
+        B(_B),
+        ud_old(_ud_old),
+        udd_old(_udd_old) {}
+
+  void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
+                     MatrixType &problem_A) const {
+    problem_rhs = this->ell;
+    /* */ B.usmv(2.0 / this->tau, ud_old, problem_rhs);
+    /* */ B.usmv(1.0, udd_old, problem_rhs);
+    this->A.usmv(-1.0, this->u_old, problem_rhs);
+    this->A.usmv(-this->tau / 2.0, ud_old, problem_rhs);
+
+    // For fixed tau, we'd only really have to do this once
+    problem_A = this->A;
+    problem_A *= this->tau / 2.0;
+    Arithmetic::addProduct(problem_A, 2.0 / this->tau, B);
+
+    // ud_old makes a good initial iterate; we could use anything, though
+    problem_iterate = ud_old;
+
+    for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
+      if (this->dirichletNodes[i].count() == dim) {
+        problem_iterate[i] = 0;
+        this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
+      } else if (this->dirichletNodes[i][1])
+        problem_iterate[i][1] = 0; // Y direction prescribed
+  }
+
+  // u1 = u0 + tau/2 (du1 + du0)
+  void virtual extractSolution(VectorType const &problem_iterate,
+                               VectorType &solution) const {
+    solution = this->u_old;
+    Arithmetic::addProduct(solution, this->tau / 2.0, problem_iterate); // ud
+    Arithmetic::addProduct(solution, this->tau / 2.0, ud_old);
+  }
+
+  void virtual extractVelocity(VectorType const &problem_iterate,
+                               VectorType &velocity) const {
+    velocity = problem_iterate;
+  }
+
+private:
+  MatrixType const &B;
+  VectorType const &ud_old;
+  VectorType const &udd_old;
+};