src/time-stepping/state_tmpl.cc

+#include "../explicitvectors.hh"
+
+template std::shared_ptr<StateUpdater<ScalarVector, Vector>>
+initStateUpdater<ScalarVector, Vector>(
+    Config::stateModel model, ScalarVector const &alpha_initial,
+    Dune::BitSetVector<1> const &frictionalNodes, double L, double V0);

src/time-stepping/updaters.hh

+#ifndef SRC_TIME_STEPPING_UPDATERS_HH
+#define SRC_TIME_STEPPING_UPDATERS_HH
+
+template <class RateUpdaterTEMPLATE, class StateUpdaterTEMPLATE>
+struct Updaters {
+  using RateUpdater = RateUpdaterTEMPLATE;
+  using StateUpdater = StateUpdaterTEMPLATE;
+
+  Updaters() {}
+
+  Updaters(std::shared_ptr<RateUpdaterTEMPLATE> rateUpdater,
+           std::shared_ptr<StateUpdaterTEMPLATE> stateUpdater)
+      : rate_(rateUpdater), state_(stateUpdater) {}
+
+  bool operator==(Updaters const &other) const {
+    return other.rate_ == rate_ and other.state_ == state_;
+  }
+
+  Updaters<RateUpdater, StateUpdater> clone() const {
+    return Updaters<RateUpdater, StateUpdater>(rate_->clone(), state_->clone());
+  }
+
+  std::shared_ptr<RateUpdaterTEMPLATE> rate_;
+  std::shared_ptr<StateUpdaterTEMPLATE> state_;
+};
+
+#endif

src/timestepping.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/istl/matrixindexset.hh>
-
-#include "timestepping.hh"
-
-#include "timestepping/backward_euler.cc"
-#include "timestepping/newmark.cc"
-
-#include "timestepping_tmpl.cc"

src/timestepping.hh

-#ifndef SRC_TIMESTEPPING_HH
-#define SRC_TIMESTEPPING_HH
-
-#include <memory>
-
-#include "enums.hh"
-#include <dune/common/bitsetvector.hh>
-
-template <class Vector, class Matrix, class Function, size_t dim>
-class TimeSteppingScheme {
-public:
-  void virtual nextTimeStep() = 0;
-  void virtual setup(Vector const &ell, double _tau, double relativeTime,
-                     Vector &rhs, Vector &iterate, Matrix &AB) = 0;
-
-  void virtual postProcess(Vector const &iterate) = 0;
-  void virtual extractDisplacement(Vector &displacement) const = 0;
-  void virtual extractVelocity(Vector &velocity) const = 0;
-  void virtual extractOldVelocity(Vector &velocity) const = 0;
-
-  std::shared_ptr<
-      TimeSteppingScheme<Vector, Matrix, Function, dim>> virtual clone()
-      const = 0;
-};
-
-#include "timestepping/newmark.hh"
-#include "timestepping/backward_euler.hh"
-
-template <class Vector, class Matrix, class Function, int dimension>
-std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dimension>>
-initTimeStepper(Config::scheme scheme,
-                Function const &velocityDirichletFunction,
-                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) {
-  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);
-    case Config::BackwardEuler:
-      return std::make_shared<
-          BackwardEuler<Vector, Matrix, Function, dimension>>(
-          stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
-          velocityDirichletNodes, velocityDirichletFunction);
-    default:
-      assert(false);
-  }
-}
-#endif

src/timestepping/backward_euler.hh

-#ifndef SRC_TIMESTEPPING_BACKWARD_EULER_HH
-#define SRC_TIMESTEPPING_BACKWARD_EULER_HH
-
-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,
-                Dune::BitSetVector<dim> const &_dirichletNodes,
-                Function const &_dirichletFunction);
-
-  void nextTimeStep() override;
-  void setup(Vector const &, double, double, Vector &, Vector &,
-             Matrix &) override;
-  void postProcess(Vector const &) override;
-  void extractDisplacement(Vector &) const override;
-  void extractVelocity(Vector &) const override;
-  void extractOldVelocity(Vector &) const override;
-
-  std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>> clone()
-      const;
-
-private:
-  Matrix const &A;
-  Matrix const &M;
-  Matrix const &C;
-  Vector u;
-  Vector v;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  Function const &dirichletFunction;
-  double dirichletValue;
-
-  Vector u_o;
-  Vector v_o;
-
-  double tau;
-
-  bool postProcessCalled = true;
-};
-#endif

src/timestepping/newmark.hh

-#ifndef SRC_TIMESTEPPING_NEWMARK_HH
-#define SRC_TIMESTEPPING_NEWMARK_HH
-
-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 &_a_initial,
-          Dune::BitSetVector<dim> const &_dirichletNodes,
-          Function const &_dirichletFunction);
-
-  void nextTimeStep() override;
-  void setup(Vector const &, double, double, Vector &, Vector &,
-             Matrix &) override;
-  void postProcess(Vector const &) override;
-  void extractDisplacement(Vector &) const override;
-  void extractVelocity(Vector &) const override;
-  void extractOldVelocity(Vector &) const override;
-
-  std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>> clone()
-      const;
-
-private:
-  Matrix const &A;
-  Matrix const &M;
-  Matrix const &C;
-  Vector u;
-  Vector v;
-  Vector a;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  Function const &dirichletFunction;
-  double dirichletValue;
-
-  Vector u_o;
-  Vector v_o;
-  Vector a_o;
-
-  double tau;
-
-  bool postProcessCalled = true;
-};
-#endif

src/uniform-grid-writer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cmath>
+#include <exception>
+#include <iostream>
+
+#include <dune/common/exceptions.hh>
+// #include <dune/common/parametertree.hh>
+// #include <dune/common/parametertreeparser.hh>
+
+#include "assemblers.hh"
+#include "diameter.hh"
+#include "gridselector.hh"
+#include "one-body-problem-data/mygrid.hh"
+#include "vtk.hh"
+
+size_t const dims = MY_DIM;
+size_t const refinements = 5; // FIXME?
+
+int main(int argc, char *argv[]) {
+  try {
+    // Dune::ParameterTree parset;
+    // Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+
+    using GridView = Grid::LeafGridView;
+    using MyAssembler = MyAssembler<GridView, dims>;
+
+    GridConstructor<Grid> gridConstructor;
+    auto grid = gridConstructor.getGrid();
+
+    // refine uniformly!
+    for (size_t refinement = 0; refinement < refinements; ++refinement)
+      grid->globalRefine(1);
+
+    double minDiameter = std::numeric_limits<double>::infinity();
+    double maxDiameter = 0.0;
+    for (auto &&e : elements(grid->leafGridView())) {
+      auto const geometry = e.geometry();
+      auto const diam = diameter(geometry);
+      minDiameter = std::min(minDiameter, diam);
+      maxDiameter = std::max(maxDiameter, diam);
+    }
+    std::cout << "min diameter: " << minDiameter << std::endl;
+    std::cout << "max diameter: " << maxDiameter << std::endl;
+
+    auto const leafView = grid->leafGridView();
+    auto const leafVertexCount = leafView.size(dims);
+
+    std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
+
+    MyAssembler const myAssembler(leafView);
+    MyVTKWriter<typename MyAssembler::VertexBasis,
+                typename MyAssembler::CellBasis> const
+        vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+    vtkWriter.writeGrid();
+  } catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  }
+}

src/vtk.cc

@@ -43,7 +43,16 @@ void MyVTKWriter<VertexBasis, CellBasis>::write(
   writer.addCellData(stressPointer);
 
   std::string const filename = prefix + std::to_string(record);
-  writer.write(filename.c_str());
+  writer.write(filename.c_str(), Dune::VTK::appendedraw);
+}
+
+template <class VertexBasis, class CellBasis>
+void MyVTKWriter<VertexBasis, CellBasis>::writeGrid() const {
+  Dune::VTKWriter<typename VertexBasis::GridView> writer(
+      vertexBasis.getGridView());
+
+  std::string const filename = prefix + "_grid";
+  writer.write(filename.c_str(), Dune::VTK::appendedraw);
 }
 
 #include "vtk_tmpl.cc"

src/vtk.hh

@@ -15,5 +15,7 @@ template <class VertexBasis, class CellBasis> class MyVTKWriter {
   template <class Vector, class ScalarVector>
   void write(size_t record, Vector const &u, Vector const &v,
              ScalarVector const &alpha, ScalarVector const &stress) const;
+
+  void writeGrid() const;
 };
 #endif