dune-tectonic.pc.in

@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
 Version: @VERSION@
 Description: dune-tectonic module
 URL: http://dune-project.org/
-Requires: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Requires: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
 Libs: -L${libdir}
 Cflags: -I${includedir}

dune.module

@@ -5,4 +5,4 @@
 Module: dune-tectonic
 Version: 2.5-1
 Maintainer: elias.pipping@fu-berlin.de
-Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Depends: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid

dune/tectonic/globalratestatefriction.hh

@@ -32,7 +32,7 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
       : restrictions(), localToGlobal(), zeroFriction() {
     auto const gridView = frictionalBoundary.gridView();
     Dune::MultipleCodimMultipleGeomTypeMapper<
-        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
+        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
     for (auto it = gridView.template begin<block_size>();
          it != gridView.template end<block_size>(); ++it) {
       auto const i = vertexMapper.index(*it);
@@ -44,8 +44,8 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
       restrictions.emplace_back(weights[i], weightedNormalStress[i],
                                 frictionInfo(geoToPoint(it->geometry())));
     }
-    assert(restrictions.size() == frictionalBoundary.numVertices());
-    assert(localToGlobal.size() == frictionalBoundary.numVertices());
+    assert(restrictions.size() == (size_t) frictionalBoundary.numVertices());
+    assert(localToGlobal.size() == (size_t) frictionalBoundary.numVertices());
   }
 
   void updateAlpha(ScalarVector const &alpha) override {

dune/tectonic/myblockproblem.hh

@@ -57,7 +57,7 @@ class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
   MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
       : Base(parset, problem),
         maxEigenvalues_(problem.f.size()),
-        localBisection(0.0, 1.0, 0.0, true, 0.0) {
+        localBisection(0.0, 1.0, 0.0, 0.0) {
     for (size_t i = 0; i < problem.f.size(); ++i) {
       LocalVectorType eigenvalues;
       Dune::FMatrixHelp::eigenValues(problem.A[i][i], eigenvalues);

src/CMakeLists.txt

@@ -19,6 +19,29 @@ set(SW_SOURCE_FILES
   time-stepping/state.cc
   vtk.cc
 )
+
+set(MSW_SOURCE_FILES
+  assemblers.cc
+  enumparser.cc
+  hdf5/frictionalboundary-writer.cc
+  hdf5/iteration-writer.cc
+  hdf5/patchinfo-writer.cc
+  hdf5/restart-io.cc
+  hdf5/surface-writer.cc
+  hdf5/time-writer.cc
+  multi-body-problem-data/cuboidgeometry.cc
+  multi-body-problem-data/mygrids.cc
+  multi-body-problem.cc
+  spatial-solving/fixedpointiterator.cc
+  spatial-solving/solverfactory.cc
+  time-stepping/adaptivetimestepper.cc
+  time-stepping/coupledtimestepper.cc
+  time-stepping/rate.cc
+  time-stepping/rate/rateupdater.cc
+  time-stepping/state.cc
+  vtk.cc
+)
+
 set(UGW_SOURCE_FILES
   assemblers.cc # FIXME
   one-body-problem-data/mygrid.cc
@@ -28,15 +51,20 @@ set(UGW_SOURCE_FILES
 
 foreach(_dim 2 3)
   set(_sw_target one-body-problem-${_dim}D)
+  set(_msw_target multi-body-problem-${_dim}D)
   set(_ugw_target uniform-grid-writer-${_dim}D)
 
   add_executable(${_sw_target} ${SW_SOURCE_FILES})
+  add_executable(${_msw_target} ${MSW_SOURCE_FILES})
   add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
   add_dune_ug_flags(${_sw_target})
+  add_dune_ug_flags(${_msw_target})
   add_dune_ug_flags(${_ugw_target})
 
   add_dune_hdf5_flags(${_sw_target})
+  add_dune_hdf5_flags(${_msw_target})
 
   set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
   set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
 endforeach()

src/assemblers_tmpl.cc

@@ -4,4 +4,4 @@
 
 #include "explicitgrid.hh"
 
-template class MyAssembler<GridView, MY_DIM>;
+template class MyAssembler<DeformedGrid::LeafGridView, MY_DIM>;

src/explicitgrid.hh

@@ -2,6 +2,13 @@
 #define SRC_EXPLICITGRID_HH
 
 #include "gridselector.hh"
+#include "explicitvectors.hh"
 
-using GridView = Grid::LeafGridView;
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+
+using LeafGridView = Grid::LeafGridView;
+using LevelGridView = Grid::LevelGridView;
+
+using DeformedGridComplex = typename Dune::Contact::DeformedContinuaComplex<Grid, Vector>;
+using DeformedGrid = DeformedGridComplex::DeformedGridType;
 #endif

src/hdf5/frictionalboundary-writer.cc

@@ -35,23 +35,27 @@ template <class ProgramState, class GridView>
 template <class Friction>
 void FrictionalBoundaryWriter<ProgramState, GridView>::write(
     ProgramState const &programState, Friction &friction) {
+    // TODO
   auto const frictionalBoundaryDisplacements =
-      restrictToSurface(programState.u, frictionalBoundary_);
+      restrictToSurface(programState.u[0], frictionalBoundary_);
   addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
            frictionalBoundaryDisplacements);
 
   auto const frictionalBoundaryVelocities =
-      restrictToSurface(programState.v, frictionalBoundary_);
+          // TODO
+      restrictToSurface(programState.v[0], frictionalBoundary_);
   addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
            frictionalBoundaryVelocities);
 
   auto const frictionalBoundaryStates =
-      restrictToSurface(programState.alpha, frictionalBoundary_);
+      // TODO
+      restrictToSurface(programState.alpha[0], frictionalBoundary_);
   addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
            frictionalBoundaryStates);
 
-  friction.updateAlpha(programState.alpha);
-  auto const c = friction.coefficientOfFriction(programState.v);
+  // TODO
+  friction.updateAlpha(programState.alpha[0]);
+  auto const c = friction.coefficientOfFriction(programState.v[0]);
   auto const frictionalBoundaryCoefficient =
       restrictToSurface(c, frictionalBoundary_);
   addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,

src/hdf5/frictionalboundary-writer_tmpl.cc

-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using MyFriction = GlobalFriction<Matrix, Vector>;
-
-template class FrictionalBoundaryWriter<MyProgramState, GridView>;
-template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
-    MyProgramState const &programState, MyFriction &friction);

src/hdf5/patchinfo-writer.cc

@@ -87,7 +87,8 @@ PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
 template <class ProgramState, class VertexBasis, class GridView>
 void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
     ProgramState const &programState) {
-  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
+    // TODO
+  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[0]);
   auto const gridVelocity = gridEvaluator_.evaluate(velocity);
   addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
 }

src/hdf5/patchinfo-writer_tmpl.cc

@@ -4,11 +4,12 @@
 #include "../program_state.hh"
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
-using P1Basis = P1NodalBasis<GridView, double>;
+using P1Basis = P1NodalBasis<DeformedGrid::LeafGridView, double>;
 using MyFunction = BasisGridFunction<P1Basis, Vector>;
 
-template class GridEvaluator<LocalVector, GridView>;
+template class GridEvaluator<LocalVector, DeformedGrid::LeafGridView>;
+
 template Dune::Matrix<typename MyFunction::RangeType>
-GridEvaluator<LocalVector, GridView>::evaluate(MyFunction const &f) const;
+GridEvaluator<LocalVector, DeformedGrid::LeafGridView>::evaluate(MyFunction const &f) const;
 
-template class PatchInfoWriter<MyProgramState, P1Basis, GridView>;
+template class PatchInfoWriter<MyProgramState, P1Basis, DeformedGrid::LeafGridView>;

src/hdf5/restart-io.cc

@@ -5,26 +5,66 @@
 #include "restart-io.hh"
 
 template <class ProgramState>
-RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
-    : displacementWriter_(group, "displacement", vertexCount,
-                          Vector::block_type::dimension),
-      velocityWriter_(group, "velocity", vertexCount,
-                      Vector::block_type::dimension),
-      accelerationWriter_(group, "acceleration", vertexCount,
-                          Vector::block_type::dimension),
-      stateWriter_(group, "state", vertexCount),
-      weightedNormalStressWriter_(group, "weightedNormalStress", vertexCount),
+RestartBodyIO<ProgramState>::RestartBodyIO(HDF5::Grouplike &group, const size_t vertexCount, const size_t id)
+    : id_(id),
+      displacementWriter_(group, "displacement"+std::to_string(id_), vertexCount,
+                          ProgramState::Vector::block_type::dimension),
+      velocityWriter_(group, "velocity"+std::to_string(id_), vertexCount,
+                      ProgramState::Vector::block_type::dimension),
+      accelerationWriter_(group, "acceleration"+std::to_string(id_), vertexCount,
+                          ProgramState::Vector::block_type::dimension),
+      stateWriter_(group, "state"+std::to_string(id_), vertexCount),
+      weightedNormalStressWriter_(group, "weightedNormalStress"+std::to_string(id_), vertexCount) {}
+
+template <class ProgramState>
+void RestartBodyIO<ProgramState>::write(const ProgramState& programState) {
+
+  addEntry(displacementWriter_, programState.timeStep, programState.u[id_]);
+  addEntry(velocityWriter_, programState.timeStep, programState.v[id_]);
+  addEntry(accelerationWriter_, programState.timeStep, programState.a[id_]);
+  addEntry(stateWriter_, programState.timeStep, programState.alpha[id_]);
+  addEntry(weightedNormalStressWriter_, programState.timeStep,
+           programState.weightedNormalStress[id_]);
+}
+
+template <class ProgramState>
+void RestartBodyIO<ProgramState>::read(size_t timeStep,
+                                   ProgramState& programState) {
+
+  readEntry(displacementWriter_, timeStep, programState.u[id_]);
+  readEntry(velocityWriter_, timeStep, programState.v[id_]);
+  readEntry(accelerationWriter_, timeStep, programState.a[id_]);
+  readEntry(stateWriter_, timeStep, programState.alpha[id_]);
+  readEntry(weightedNormalStressWriter_, timeStep,
+            programState.weightedNormalStress[id_]);
+}
+
+template <class ProgramState>
+RestartIO<ProgramState>::RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts)
+    : bodyIO_(vertexCounts.size()),
       relativeTimeWriter_(group, "relativeTime"),
-      relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {}
+      relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {
+
+  for (size_t i=0; i<bodyIO_.size(); i++) {
+    bodyIO_[i] = new RestartBodyIO<ProgramState>(group, vertexCounts[i], i);
+  }
+}
+
+template <class ProgramState>
+RestartIO<ProgramState>::~RestartIO() {
+  for (size_t i=0; i<bodyIO_.size(); i++) {
+    delete bodyIO_[i];
+  }
+}
 
 template <class ProgramState>
 void RestartIO<ProgramState>::write(ProgramState const &programState) {
-  addEntry(displacementWriter_, programState.timeStep, programState.u);
-  addEntry(velocityWriter_, programState.timeStep, programState.v);
-  addEntry(accelerationWriter_, programState.timeStep, programState.a);
-  addEntry(stateWriter_, programState.timeStep, programState.alpha);
-  addEntry(weightedNormalStressWriter_, programState.timeStep,
-           programState.weightedNormalStress);
+  assert(programState.size() == bodyIO_.size());
+
+  for (size_t i=0; i<bodyIO_.size(); i++) {
+    bodyIO_[i]->write(programState);
+  }
+
   addEntry(relativeTimeWriter_, programState.timeStep,
            programState.relativeTime);
   addEntry(relativeTimeIncrementWriter_, programState.timeStep,
@@ -33,14 +73,15 @@ void RestartIO<ProgramState>::write(ProgramState const &programState) {
 
 template <class ProgramState>
 void RestartIO<ProgramState>::read(size_t timeStep,
-                                   ProgramState &programState) {
+                                   ProgramState& programState) {
+  assert(programState.size() == bodyIO_.size());
+
   programState.timeStep = timeStep;
-  readEntry(displacementWriter_, timeStep, programState.u);
-  readEntry(velocityWriter_, timeStep, programState.v);
-  readEntry(accelerationWriter_, timeStep, programState.a);
-  readEntry(stateWriter_, timeStep, programState.alpha);
-  readEntry(weightedNormalStressWriter_, timeStep,
-            programState.weightedNormalStress);
+
+  for (size_t i=0; i<bodyIO_.size(); i++) {
+    bodyIO_[i]->read(timeStep, programState);
+  }
+
   readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
   readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
 }

src/hdf5/restart-io.hh

@@ -4,23 +4,37 @@
 #include <dune/fufem/hdf5/file.hh>
 #include <dune/fufem/hdf5/sequenceio.hh>
 
-template <class ProgramState> class RestartIO {
-  using ScalarVector = typename ProgramState::ScalarVector;
-  using Vector = typename ProgramState::Vector;
-
+template <class ProgramState> class RestartBodyIO {
 public:
-  RestartIO(HDF5::Grouplike &group, size_t vertexCount);
+  RestartBodyIO(HDF5::Grouplike& group, const size_t vertexCount, const size_t id);
 
-  void write(ProgramState const &programState);
+  void write(const ProgramState& programState);
 
-  void read(size_t timeStep, ProgramState &programState);
+  void read(size_t timeStep, ProgramState& programState);
 
 private:
+  const size_t id_;
+
   HDF5::SequenceIO<2> displacementWriter_;
   HDF5::SequenceIO<2> velocityWriter_;
   HDF5::SequenceIO<2> accelerationWriter_;
   HDF5::SequenceIO<1> stateWriter_;
   HDF5::SequenceIO<1> weightedNormalStressWriter_;
+};
+
+template <class ProgramState> class RestartIO {
+
+public:
+  RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts);
+  ~RestartIO();
+
+  void write(const ProgramState& programState);
+
+  void read(size_t timeStep, ProgramState& programState);
+
+private:
+  std::vector<RestartBodyIO<ProgramState>* > bodyIO_;
+
   HDF5::SequenceIO<0> relativeTimeWriter_;
   HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
 };

src/hdf5/restart-io_tmpl.cc

@@ -4,4 +4,5 @@
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
 
+template class RestartBodyIO<MyProgramState>;
 template class RestartIO<MyProgramState>;

src/hdf5/surface-writer.cc

@@ -24,11 +24,13 @@ SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
 template <class ProgramState, class GridView>
 void SurfaceWriter<ProgramState, GridView>::write(
     ProgramState const &programState) {
-  auto const surfaceDisplacements = restrictToSurface(programState.u, surface_);
+    // TODO
+  auto const surfaceDisplacements = restrictToSurface(programState.u[0], surface_);
   addEntry(surfaceDisplacementWriter_, programState.timeStep,
            surfaceDisplacements);
 
-  auto const surfaceVelocities = restrictToSurface(programState.v, surface_);
+  // TODO
+  auto const surfaceVelocities = restrictToSurface(programState.v[0], surface_);
   addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
 }
 

src/hdf5/surface-writer_tmpl.cc

@@ -5,4 +5,4 @@
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
 
-template class SurfaceWriter<MyProgramState, GridView>;
+template class SurfaceWriter<MyProgramState, LeafGridView>;

src/matrices.hh

 #ifndef SRC_MATRICES_HH
 #define SRC_MATRICES_HH
 
-template <class Matrix> struct Matrices {
-  Matrix elasticity;
-  Matrix damping;
-  Matrix mass;
+template <class Matrix> class Matrices {
+public:
+  std::vector<std::shared_ptr<Matrix>> elasticity;
+  std::vector<std::shared_ptr<Matrix>> damping;
+  std::vector<std::shared_ptr<Matrix>> mass;
+
+  Matrices(size_t n) {
+      elasticity.resize(n);
+      damping.resize(n);
+      mass.resize(n);
+  }
 };
 #endif

src/multi-body-problem-2D.cfg

+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter= 2e-3  # [m]
+
+[timeSteps]
+refinementTolerance = 1e-5
+
+[u0.solver]
+tolerance         = 1e-8
+
+[a0.solver]
+tolerance         = 1e-8
+
+[v.solver]
+tolerance         = 1e-8
+
+[v.fpi]
+tolerance         = 1e-5
+
+[solver.tnnmg.linear]
+tolerance          = 1e-10

src/multi-body-problem-3D.cfg

+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter= 2e-2  # [m]
+
+[boundary.friction.weakening]
+patchType       = Trapezoidal
+
+[timeSteps]
+refinementTolerance = 1e-5
+
+[u0.solver]
+tolerance         = 1e-6
+
+[a0.solver]
+tolerance         = 1e-6
+
+[v.solver]
+tolerance         = 1e-6
+
+[v.fpi]
+tolerance         = 1e-5
+
+[solver.tnnmg.linear]
+tolerance          = 1e-10

src/multi-body-problem-data/bc.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+
+class VelocityDirichletCondition
+    : public Dune::VirtualFunction<double, double> {
+  void evaluate(double const &relativeTime, double &y) const {
+    // Assumed to vanish at time zero
+    double const finalVelocity = -5e-5;
+    y = (relativeTime <= 0.1)
+            ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
+            : finalVelocity;
+  }
+};
+
+class NeumannCondition : public Dune::VirtualFunction<double, double> {
+  void evaluate(double const &relativeTime, double &y) const { y = 0.0; }
+};
+#endif