CMakeLists.txt

@@ -17,7 +17,7 @@ include(DuneMacros)
 
 # start a dune project with information from dune.module
 dune_project()
-
+dune_enable_all_packages()
 find_package(HDF5 COMPONENTS C REQUIRED)
 
 add_subdirectory("src")

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/CMakeLists.txt

 install(FILES
-  body.hh
+  bodydata.hh
   frictiondata.hh
   frictionpotential.hh
   globalfrictiondata.hh
@@ -12,4 +12,5 @@ install(FILES
   mydirectionalconvexfunction.hh
   quadraticenergy.hh
   tectonic.hh
+  transformedglobalratestatefriction.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/body.hh → dune/tectonic/bodydata.hh

-#ifndef DUNE_TECTONIC_BODY_HH
-#define DUNE_TECTONIC_BODY_HH
+#ifndef DUNE_TECTONIC_BODY_DATA_HH
+#define DUNE_TECTONIC_BODY_DATA_HH
 
-template <int dimension> struct Body {
+template <int dimension> struct BodyData {
   using ScalarFunction =
       Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
                             Dune::FieldVector<double, 1>>;

dune/tectonic/frictionpotential.hh

@@ -106,6 +106,9 @@ class RegularisedRateState : public FrictionPotential {
 
 class ZeroFunction : public FrictionPotential {
 public:
+  template <typename... Args>
+  ZeroFunction(Args... args) {}
+
   double evaluate(double) const override { return 0; }
 
   double coefficientOfFriction(double s) const override { return 0; }

dune/tectonic/globalfriction.hh

@@ -9,9 +9,12 @@
 
 #include <dune/solvers/common/interval.hh>
 
+#include <dune/tnnmg/functionals/nonsmoothconvexfunctional.hh>
+
 #include <dune/tectonic/localfriction.hh>
 
-template <class Matrix, class Vector> class GlobalFriction {
+template <class Matrix, class Vector>
+class GlobalFriction { //: public Dune::TNNMG::NonsmoothConvexFunctional<> {
 protected:
   using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
 
@@ -81,6 +84,7 @@ template <class Matrix, class Vector> class GlobalFriction {
     return ret;
   }
 
-  void virtual updateAlpha(ScalarVector const &alpha) = 0;
+  void virtual updateAlpha(const std::vector<ScalarVector>& alpha) = 0;
+
 };
 #endif

dune/tectonic/globalfrictiondata.hh

@@ -27,6 +27,7 @@ template <int dimension> class GlobalFrictionData {
       Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
                             Dune::FieldVector<double, 1>>;
 
+public:
   double virtual const &C() const = 0;
   double virtual const &L() const = 0;
   double virtual const &V0() const = 0;

dune/tectonic/globalratestatefriction.hh

@@ -10,62 +10,123 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+
 #include <dune/tectonic/geocoordinate.hh>
 #include <dune/tectonic/globalfrictiondata.hh>
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/index-in-sorted-range.hh>
 
-template <class Matrix, class Vector, class ScalarFriction, class GridView>
+#include "../../src/frictioncouplingpair.hh"
+
+template <class Matrix, class Vector, class ScalarFriction, class GridType>
 class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
 public:
   using GlobalFriction<Matrix, Vector>::block_size;
   using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
 
 private:
-  using typename GlobalFriction<Matrix, Vector>::ScalarVector;
+  using Base = GlobalFriction<Matrix, Vector>;
+
+  using field_type = typename Vector::field_type;
+  using typename Base::ScalarVector;
+  using typename Base::LocalVectorType;
+
+  using FrictionCouplingPair = FrictionCouplingPair<GridType, LocalVectorType, field_type>;
+  using ContactCoupling =  Dune::Contact::DualMortarCoupling<field_type, GridType>;
+
+  size_t bodyIndex(const size_t globalIdx) {
+     size_t i=0;
+
+     for (; i<maxIndex_.size(); i++) {
+         if (globalIdx < maxIndex_[i])
+             break;
+     }
+     return i;
+  }
 
 public:
-  GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
-                          GlobalFrictionData<block_size> const &frictionInfo,
-                          ScalarVector const &weights,
-                          ScalarVector const &weightedNormalStress)
-      : restrictions(), localToGlobal(), zeroFriction() {
-    auto const gridView = frictionalBoundary.gridView();
-    Dune::MultipleCodimMultipleGeomTypeMapper<
-        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
-    for (auto it = gridView.template begin<block_size>();
-         it != gridView.template end<block_size>(); ++it) {
-      auto const i = vertexMapper.index(*it);
-
-      if (not frictionalBoundary.containsVertex(i))
-        continue;
-
-      localToGlobal.emplace_back(i);
-      restrictions.emplace_back(weights[i], weightedNormalStress[i],
-                                frictionInfo(geoToPoint(it->geometry())));
+  GlobalRateStateFriction(const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+                          const std::vector<std::shared_ptr<FrictionCouplingPair>>& couplings, // contains frictionInfo
+                          const std::vector<ScalarVector>& weights,
+                          const std::vector<ScalarVector>& weightedNormalStress)
+      : restrictions_(), localToGlobal_(), zeroFriction_() {
+
+    assert(contactCouplings.size() == couplings.size());
+    assert(weights.size() == weightedNormalStress.size());
+
+    const auto nBodies = weights.size();
+
+    std::vector<std::vector<int>> nonmortarBodies(nBodies); // first index body, second index coupling
+    for (size_t i=0; i<contactCouplings.size(); i++) {
+        const auto nonmortarIdx = couplings[i]->gridIdx_[0];
+        nonmortarBodies[nonmortarIdx].emplace_back(i);
+    }
+
+    maxIndex_.resize(nBodies);
+    size_t globalIdx = 0;
+
+    for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+        const auto& couplingIndices = nonmortarBodies[bodyIdx];
+
+        if (couplingIndices.size()==0)
+            continue;
+
+        const auto gridView = contactCouplings[couplingIndices[0]]->nonmortarBoundary().gridView();
+
+        Dune::MultipleCodimMultipleGeomTypeMapper<
+            decltype(gridView), Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
+
+        for (auto it = gridView.template begin<block_size>(); it != gridView.template end<block_size>(); ++it) {
+            const auto i = vertexMapper.index(*it);
+
+            for (size_t j=0; j<couplingIndices.size(); j++) {
+                const auto couplingIdx = couplingIndices[j];
+
+                if (not contactCouplings[couplingIdx]->nonmortarBoundary().containsVertex(i))
+                  continue;
+
+                localToGlobal_.emplace_back(i);
+                restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
+                                          couplings[j]->frictionData()(geoToPoint(it->geometry())));
+                break;
+            }
+
+            globalIdx++;
+        }
+        maxIndex_[bodyIdx] = globalIdx;
     }
-    assert(restrictions.size() == frictionalBoundary.numVertices());
-    assert(localToGlobal.size() == frictionalBoundary.numVertices());
   }
 
-  void updateAlpha(ScalarVector const &alpha) override {
-    for (size_t j = 0; j < restrictions.size(); ++j)
-      restrictions[j].updateAlpha(alpha[localToGlobal[j]]);
+  void updateAlpha(const std::vector<ScalarVector>& alpha) override {
+    for (size_t j = 0; j < restrictions_.size(); ++j) {
+      const auto bodyIdx = bodyIndex(j);
+      const auto globalDof = localToGlobal_[j];
+
+      size_t bodyDof;
+      if (bodyIdx>0) {
+          bodyDof = globalDof - maxIndex_[bodyIdx-1];
+      } else {
+          bodyDof = globalDof;
+      }
+      restrictions_[j].updateAlpha(alpha[bodyIdx][bodyDof]);
+    }
   }
 
   /*
     Return a restriction of the outer function to the i'th node.
   */
   LocalNonlinearity const &restriction(size_t i) const override {
-    auto const index = indexInSortedRange(localToGlobal, i);
-    if (index == localToGlobal.size())
-      return zeroFriction;
-    return restrictions[index];
+    auto const index = indexInSortedRange(localToGlobal_, i);
+    if (index == localToGlobal_.size())
+      return zeroFriction_;
+    return restrictions_[index];
   }
 
 private:
-  std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions;
-  std::vector<size_t> localToGlobal;
-  WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction;
+  std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions_;
+  std::vector<size_t> maxIndex_; // max global index per body
+  std::vector<size_t> localToGlobal_;
+  WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction_;
 };
 #endif

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);

dune/tectonic/mydirectionalconvexfunction.hh

-#ifndef DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
-#define DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
-
-// Copied from dune/tnnmg/problem-classes/directionalconvexfunction.hh
-// Allows phi to be const
-
-#include <dune/fufem/arithmetic.hh>
-#include <dune/solvers/common/interval.hh>
-
-template <class Nonlinearity> class MyDirectionalConvexFunction {
-public:
-  using Vector = typename Nonlinearity::VectorType;
-  using Matrix = typename Nonlinearity::MatrixType;
-
-  MyDirectionalConvexFunction(double A, double b, Nonlinearity const &phi,
-                              Vector const &u, Vector const &v)
-      : A(A), b(b), phi(phi), u(u), v(v) {
-    phi.directionalDomain(u, v, dom);
-  }
-
-  double quadraticPart() const { return A; }
-
-  double linearPart() const { return b; }
-
-  void subDiff(double x, Dune::Solvers::Interval<double> &D) const {
-    Vector uxv = u;
-    Arithmetic::addProduct(uxv, x, v);
-    phi.directionalSubDiff(uxv, v, D);
-    auto const Axmb = A * x - b;
-    D[0] += Axmb;
-    D[1] += Axmb;
-  }
-
-  void domain(Dune::Solvers::Interval<double> &domain) const {
-    domain[0] = this->dom[0];
-    domain[1] = this->dom[1];
-  }
-
-  double A;
-  double b;
-
-private:
-  Nonlinearity const &phi;
-  Vector const &u;
-  Vector const &v;
-
-  Dune::Solvers::Interval<double> dom;
-};
-
-/*
-  1/2 <A(u + hv),u + hv> - <b, u + hv>
-  = 1/2 <Av,v> h^2 - <b - Au, v> h + const.
-
-  localA = <Av,v>
-  localb = <b - Au, v>
-*/
-
-template <class Matrix, class Vector, class Nonlinearity>
-MyDirectionalConvexFunction<Nonlinearity> restrict(Matrix const &A,
-                                                   Vector const &b,
-                                                   Vector const &u,
-                                                   Vector const &v,
-                                                   Nonlinearity const &phi) {
-  return MyDirectionalConvexFunction<Nonlinearity>(
-      Arithmetic::Axy(A, v, v), Arithmetic::bmAxy(A, b, u, v), phi, u, v);
-}
-
-#endif

dune/tectonic/tectonic.hh

-#ifndef DUNE_TECTONIC_TECTONIC_HH
-#define DUNE_TECTONIC_TECTONIC_HH
-
-// add your classes here
-#endif

program_structure.txt

+--------------------------
+--  Program structure:  --
+--------------------------
+
+
+1. build n-body system 
+    - contains grids, couplings, gridViews, assembler 
+    
+    Data structure: LevelContactNetwork
+    Factories:      StackedBlocksFactory
+    
+2. initialize/set up program state
+    - holds bodyStates, u, v, a, alpha for each body
+    - defines time, timeStep
+    - computes initial conditions
+    
+    Data structure: ProgramState
+
+-- tested until here
+
+3. assemble RSD friction
+
+4. set up TNNMG solver
+    - rate updater
+    - state updater
+    
+5. adaptive time stepper
+
+ 

src/CMakeLists.txt

+add_subdirectory("tests")
+add_subdirectory("spatial-solving")
+
 set(SW_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
-  one-body-problem-data/mygeometry.cc
-  one-body-problem-data/mygrid.cc
-  one-body-problem.cc
+  data-structures/body.cc
+  data-structures/enumparser.cc
+  io/vtk.cc
+  io/hdf5/frictionalboundary-writer.cc
+  io/hdf5/iteration-writer.cc
+  io/hdf5/patchinfo-writer.cc
+  io/hdf5/restart-io.cc
+  io/hdf5/surface-writer.cc
+  io/hdf5/time-writer.cc
+#  one-body-problem-data/mygeometry.cc
+#  one-body-problem-data/mygrid.cc
+#  one-body-problem.cc
   spatial-solving/fixedpointiterator.cc
   spatial-solving/solverfactory.cc
   time-stepping/adaptivetimestepper.cc
@@ -17,26 +22,64 @@ set(SW_SOURCE_FILES
   time-stepping/rate.cc
   time-stepping/rate/rateupdater.cc
   time-stepping/state.cc
-  vtk.cc
 )
+
+set(MSW_SOURCE_FILES
+  assemblers.cc
+  #nodalweights.cc
+  data-structures/body.cc
+  data-structures/levelcontactnetwork.cc
+  data-structures/contactnetwork.cc
+  data-structures/enumparser.cc
+  #factories/cantorfactory.cc
+  factories/threeblocksfactory.cc
+  factories/stackedblocksfactory.cc
+  io/vtk.cc
+  io/hdf5/frictionalboundary-writer.cc
+  io/hdf5/iteration-writer.cc
+  io/hdf5/patchinfo-writer.cc
+  io/hdf5/restart-io.cc
+  io/hdf5/surface-writer.cc
+  io/hdf5/time-writer.cc
+  #multi-body-problem-data/grid/cube.cc
+  #multi-body-problem-data/grid/cubefaces.cc
+  multi-body-problem-data/grid/cuboidgeometry.cc
+  multi-body-problem-data/grid/mygrids.cc
+  multi-body-problem-data/grid/simplexmanager.cc
+  multi-body-problem.cc
+  #spatial-solving/solverfactory.cc
+  spatial-solving/fixedpointiterator.cc
+  #spatial-solving/solverfactory_old.cc
+  time-stepping/adaptivetimestepper.cc
+  time-stepping/coupledtimestepper.cc
+  time-stepping/rate.cc
+  time-stepping/rate/rateupdater.cc
+  time-stepping/state.cc
+)
+
 set(UGW_SOURCE_FILES
   assemblers.cc # FIXME
+  io/uniform-grid-writer.cc
+  io/vtk.cc
   one-body-problem-data/mygrid.cc
-  uniform-grid-writer.cc
-  vtk.cc
 )
 
 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.cc

@@ -12,7 +12,6 @@
 #include <dune/fufem/assemblers/localassemblers/variablecoefficientviscosityassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
-#include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
 #include <dune/fufem/functiontools/p0p1interpolation.hh>
@@ -24,17 +23,28 @@
 #include "assemblers.hh"
 
 template <class GridView, int dimension>
-MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
-    : cellBasis(_gridView),
+MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView) :
+      cellBasis(_gridView),
       vertexBasis(_gridView),
       gridView(_gridView),
       cellAssembler(cellBasis, cellBasis),
       vertexAssembler(vertexBasis, vertexBasis) {}
 
+template <class GridView, int dimension>
+template <class LocalBoundaryFunctionalAssemblerType, class GlobalVectorType>
+void MyAssembler<GridView, dimension>::assembleBoundaryFunctional(LocalBoundaryFunctionalAssemblerType& localAssembler,
+                                GlobalVectorType& b,
+                                const BoundaryPatch<GridView>& boundaryPatch,
+                                bool initializeVector) const {
+
+    vertexAssembler.assembleBoundaryFunctional(localAssembler, b, boundaryPatch, initializeVector);
+}
+
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
-    BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarMatrix &frictionalBoundaryMass) const {
+        BoundaryPatch<GridView> const &frictionalBoundary,
+        ScalarMatrix &frictionalBoundaryMass) const {
+
   BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
                         LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
       frictionalBoundaryMassAssembler(frictionalBoundary);
@@ -44,9 +54,9 @@ void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleMass(
-    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-        densityFunction,
-    Matrix &M) const {
+        Dune::VirtualFunction<LocalVector, LocalScalarVector> const & densityFunction,
+        Matrix &M) const {
+
   // NOTE: We treat the weight as a constant function
   QuadratureRuleKey quadKey(dimension, 0);
 
@@ -58,8 +68,11 @@ void MyAssembler<GridView, dimension>::assembleMass(
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
-                                                          Matrix &A) const {
+void MyAssembler<GridView, dimension>::assembleElasticity(
+        double E,
+        double nu,
+        Matrix &A) const {
+
   StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
       localStiffness(E, nu);
   vertexAssembler.assembleOperator(localStiffness, A);
@@ -67,9 +80,10 @@ void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleViscosity(
-    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
-    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
-    Matrix &C) const {
+        Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
+        Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
+        Matrix &C) const {
+
   // NOTE: We treat the weights as constant functions
   QuadratureRuleKey shearViscosityKey(dimension, 0);
   QuadratureRuleKey bulkViscosityKey(dimension, 0);
@@ -83,8 +97,9 @@ void MyAssembler<GridView, dimension>::assembleViscosity(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleBodyForce(
-    Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-    Vector &f) const {
+        Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+        Vector &f) const {
+
   L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
       gravityFunctionalAssembler(gravityField);
   vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
@@ -92,9 +107,11 @@ void MyAssembler<GridView, dimension>::assembleBodyForce(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleNeumann(
-    BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
-    Dune::VirtualFunction<double, double> const &neumann,
-    double relativeTime) const {
+        BoundaryPatch<GridView> const &neumannBoundary,
+        Vector &f,
+        Dune::VirtualFunction<double, double> const &neumann,
+        double relativeTime) const {
+
   LocalVector localNeumann(0);
   neumann.evaluate(relativeTime, localNeumann[0]);
   NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
@@ -106,9 +123,11 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
-    BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarVector &weightedNormalStress, double youngModulus,
-    double poissonRatio, Vector const &displacement) const {
+        BoundaryPatch<GridView> const &frictionalBoundary,
+        ScalarVector &weightedNormalStress,
+        double youngModulus,
+        double poissonRatio,
+        Vector const &displacement) const {
 
   BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
                                                               displacement);
@@ -136,41 +155,13 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
         std::fmin(normals[i] * nodalTractionAverage[i], 0) * weights[i];
 }
 
-template <class GridView, int dimension>
-auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
-    Config::FrictionModel frictionModel,
-    BoundaryPatch<GridView> const &frictionalBoundary,
-    GlobalFrictionData<dimension> const &frictionInfo,
-    ScalarVector const &weightedNormalStress) const
-    -> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
-  // Lumping of the nonlinearity
-  ScalarVector weights;
-  {
-    NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
-        frictionalBoundaryAssembler(std::make_shared<
-            ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
-            1));
-    vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
-                                               weights, frictionalBoundary);
-  }
-  switch (frictionModel) {
-    case Config::Truncated:
-      return std::make_shared<GlobalRateStateFriction<
-          Matrix, Vector, TruncatedRateState, GridView>>(
-          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
-    case Config::Regularised:
-      return std::make_shared<GlobalRateStateFriction<
-          Matrix, Vector, RegularisedRateState, GridView>>(
-          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
-    default:
-      assert(false);
-  }
-}
-
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleVonMisesStress(
-    double youngModulus, double poissonRatio, Vector const &u,
-    ScalarVector &stress) const {
+        double youngModulus,
+        double poissonRatio,
+        Vector const &u,
+        ScalarVector &stress) const {
+
   auto const gridDisplacement =
       std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
           vertexBasis, u);

src/assemblers.hh

@@ -12,79 +12,96 @@
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #pragma clang diagnostic pop
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/boundarypatch.hh>
 
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/globalfrictiondata.hh>
 
-#include "enums.hh"
+#include "data-structures/enums.hh"
 
-template <class GridView, int dimension> class MyAssembler {
+template <class GridView, int dimension>
+class MyAssembler {
 public:
-  using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
-  using Matrix =
+    using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using Matrix =
       Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
-  using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
-  using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
+    using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
 
-  using CellBasis = P0Basis<GridView, double>;
-  using VertexBasis = P1NodalBasis<GridView, double>;
+    using CellBasis = P0Basis<GridView, double>;
+    using VertexBasis = P1NodalBasis<GridView, double>;
 
-  CellBasis const cellBasis;
-  VertexBasis const vertexBasis;
-  GridView const &gridView;
+    CellBasis const cellBasis;
+    VertexBasis const vertexBasis;
+    GridView const &gridView;
 
 private:
-  using Grid = typename GridView::Grid;
-  using LocalVector = typename Vector::block_type;
-  using LocalScalarVector = typename ScalarVector::block_type;
+    using Grid = typename GridView::Grid;
+    using LocalVector = typename Vector::block_type;
+    using LocalScalarVector = typename ScalarVector::block_type;
 
-  using LocalCellBasis = typename CellBasis::LocalFiniteElement;
-  using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
+    using LocalCellBasis = typename CellBasis::LocalFiniteElement;
+    using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
 
-  Assembler<CellBasis, CellBasis> cellAssembler;
-  Assembler<VertexBasis, VertexBasis> vertexAssembler;
+    Assembler<CellBasis, CellBasis> cellAssembler;
+    Assembler<VertexBasis, VertexBasis> vertexAssembler;
 
 public:
-  MyAssembler(GridView const &gridView);
-
-  void assembleFrictionalBoundaryMass(
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarMatrix &frictionalBoundaryMass) const;
-
-  void assembleMass(Dune::VirtualFunction<
-                        LocalVector, LocalScalarVector> const &densityFunction,
-                    Matrix &M) const;
-
-  void assembleElasticity(double E, double nu, Matrix &A) const;
-
-  void assembleViscosity(
-      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-          shearViscosity,
-      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-          bulkViscosity,
-      Matrix &C) const;
-
-  void assembleBodyForce(
-      Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-      Vector &f) const;
-
-  void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
-                       Vector &f,
-                       Dune::VirtualFunction<double, double> const &neumann,
-                       double relativeTime) const;
-
-  void assembleWeightedNormalStress(
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarVector &weightedNormalStress, double youngModulus,
-      double poissonRatio, Vector const &displacement) const;
-
-  std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
-      Config::FrictionModel frictionModel,
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      GlobalFrictionData<dimension> const &frictionInfo,
-      ScalarVector const &weightedNormalStress) const;
-
-  void assembleVonMisesStress(double youngModulus, double poissonRatio,
-                              Vector const &u, ScalarVector &stress) const;
+    MyAssembler(GridView const &gridView);
+
+    template <class LocalBoundaryFunctionalAssemblerType, class GlobalVectorType>
+    void assembleBoundaryFunctional(LocalBoundaryFunctionalAssemblerType& localAssembler,
+                                    GlobalVectorType& b,
+                                    const BoundaryPatch<GridView>& boundaryPatch,
+                                    bool initializeVector=true) const;
+
+    void assembleFrictionalBoundaryMass(
+            BoundaryPatch<GridView> const &frictionalBoundary,
+            ScalarMatrix &frictionalBoundaryMass) const;
+
+    void assembleMass(
+            Dune::VirtualFunction<
+            LocalVector, LocalScalarVector> const &densityFunction,
+            Matrix &M) const;
+
+    void assembleElasticity(
+            double E,
+            double nu,
+            Matrix &A) const;
+
+    void assembleViscosity(
+            Dune::VirtualFunction<LocalVector, LocalScalarVector> const & shearViscosity,
+            Dune::VirtualFunction<LocalVector, LocalScalarVector> const & bulkViscosity,
+            Matrix &C) const;
+
+    void assembleBodyForce(
+            Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+            Vector &f) const;
+
+    void assembleNeumann(
+            BoundaryPatch<GridView> const &neumannBoundary,
+            Vector &f,
+            Dune::VirtualFunction<double, double> const &neumann,
+            double relativeTime) const;
+
+    void assembleWeightedNormalStress(
+            BoundaryPatch<GridView> const &frictionalBoundary,
+            ScalarVector &weightedNormalStress,
+            double youngModulus,
+            double poissonRatio,
+            Vector const &displacement) const;
+
+    auto assembleFrictionNonlinearity(
+            Config::FrictionModel frictionModel,
+            BoundaryPatch<GridView> const &frictionalBoundary,
+            GlobalFrictionData<dimension> const &frictionInfo,
+            ScalarVector const &weightedNormalStress) const
+    -> std::shared_ptr<GlobalFriction<Matrix, Vector>>;
+
+    void assembleVonMisesStress(
+            double youngModulus,
+            double poissonRatio,
+            Vector const &u,
+            ScalarVector &stress) const;
 };
 #endif

src/assemblers_tmpl.cc

@@ -3,5 +3,20 @@
 #endif
 
 #include "explicitgrid.hh"
+#include "explicitvectors.hh"
 
-template class MyAssembler<GridView, MY_DIM>;
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/boundarypatch.hh>
+
+#include "assemblers.hh"
+
+template class MyAssembler<DefLeafGridView, MY_DIM>;
+
+
+using MyNeumannBoundaryAssembler = NeumannBoundaryAssembler<DeformedGrid, typename ScalarVector::block_type>;
+
+template void MyAssembler<DefLeafGridView, MY_DIM>::assembleBoundaryFunctional<MyNeumannBoundaryAssembler, ScalarVector>(
+        MyNeumannBoundaryAssembler& localAssembler,
+        ScalarVector& b,
+        const BoundaryPatch<DefLeafGridView>& boundaryPatch,
+        bool initializeVector) const;

src/boundarycondition.hh

+#ifndef SRC_BOUNDARYCONDITION_HH
+#define SRC_BOUNDARYCONDITION_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+
+template <class GridView, int dims>
+class BoundaryCondition {
+public:
+  using BoundaryPatch = BoundaryPatch<GridView>;
+
+private:
+  using Function = Dune::VirtualFunction<double, double>;
+
+  const std::string tag_; // store type of boundary condition, e.g. dirichlet, neumann, friction, etc
+
+  std::shared_ptr<BoundaryPatch> boundaryPatch_;
+  std::shared_ptr<Dune::BitSetVector<dims>> boundaryNodes_;
+  std::shared_ptr<Function> boundaryFunction_;
+
+public:
+  BoundaryCondition(const std::string& tag = "") :
+      tag_(tag)
+  {}
+
+  BoundaryCondition(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function, const std::string& tag = "") :
+      tag_(tag),
+      boundaryPatch_(patch),
+      boundaryFunction_(function)
+  {}
+
+  void setBoundaryPatch(const GridView& gridView, std::shared_ptr<Dune::BitSetVector<dims>> nodes) {
+      boundaryNodes_ = nodes;
+      boundaryPatch_ = std::make_shared<BoundaryPatch>(gridView, *nodes);
+  }
+
+  void setBoundaryPatch(std::shared_ptr<BoundaryPatch> patch) {
+      boundaryPatch_ = patch;
+  }
+
+  void setBoundaryPatch(const BoundaryPatch& patch) {
+      boundaryPatch_ = std::make_shared<BoundaryPatch>(patch);
+  }
+
+  void setBoundaryFunction(std::shared_ptr<Function> function) {
+      boundaryFunction_ = function;
+  }
+
+  void set(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function) {
+      boundaryPatch_ = patch;
+      boundaryFunction_ = function;
+  }
+
+  const std::string& tag() const {
+      return tag_;
+  }
+
+  const std::shared_ptr<BoundaryPatch>& boundaryPatch() const {
+      return boundaryPatch_;
+  }
+
+  const std::shared_ptr<Dune::BitSetVector<dims>>& boundaryNodes() const {
+      return boundaryNodes_;
+  }
+
+  const std::shared_ptr<Function>& boundaryFunction() const {
+      return boundaryFunction_;
+  }
+};
+
+#endif

src/data-structures/.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+
+template class ContactNetwork<Grid, MY_DIM>;

src/data-structures/body.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/hybridutilities.hh>
+#include "body.hh"
+
+// -------------------------------
+// --- LeafBody Implementation ---
+// -------------------------------
+template <class GridTEMPLATE, class VectorType>
+LeafBody<GridTEMPLATE, VectorType>::LeafBody(
+        const std::shared_ptr<BodyData<dim>>& bodyData,
+        const std::shared_ptr<GridTEMPLATE>& hostGrid) :
+    bodyData_(bodyData),
+    hostGrid_(hostGrid),
+    deformation_(std::make_shared<DeformationFunction>(hostGrid_->leafGridView())),
+    grid_(std::make_shared<GridType>(*hostGrid_, *deformation_)),
+    gridView_(grid_->leafGridView()),
+    assembler_(std::make_shared<Assembler>(gridView_)),
+    matrices_()
+{
+
+    boundaryConditions_.resize(0);
+
+    externalForce_ = [&](double relativeTime, VectorType &ell) {
+        // Problem formulation: right-hand side
+        std::vector<std::shared_ptr<BoundaryCondition>> leafNeumannConditions;
+        boundaryConditions("neumann", leafNeumannConditions);
+        ell.resize(gravityFunctional_.size());
+        ell = gravityFunctional_;
+
+        for (size_t i=0; i<leafNeumannConditions.size(); i++) {
+            const auto& leafNeumannCondition = leafNeumannConditions[i];
+
+            VectorType b;
+            assembler_->assembleNeumann(*leafNeumannCondition->boundaryPatch(), b, *leafNeumannCondition->boundaryFunction(),
+                                    relativeTime);
+            ell += b;
+        }
+    };
+}
+
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::assemble() {
+
+    // assemble matrices
+    assembler_->assembleElasticity(bodyData_->getYoungModulus(), bodyData_->getPoissonRatio(), *matrices_.elasticity);
+    assembler_->assembleViscosity(bodyData_->getShearViscosityField(), bodyData_->getBulkViscosityField(), *matrices_.damping);
+    assembler_->assembleMass(bodyData_->getDensityField(), *matrices_.mass);
+
+    // assemble forces
+    assembler_->assembleBodyForce(bodyData_->getGravityField(), gravityFunctional_);
+}
+
+// setter and getter
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::data() const
+-> const std::shared_ptr<BodyData<dim>>& {
+
+    return bodyData_;
+}
+
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::setDeformation(const VectorType& def) {
+    deformation_->setDeformation(def);
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::deformation() const
+-> DeformationFunction& {
+
+    return *deformation_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::grid() const
+-> std::shared_ptr<GridType> {
+
+    return grid_;
+}
+
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::gridView() const
+-> const GridView& {
+
+    return gridView_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::nVertices() const
+-> size_t {
+    return gridView_.size(dim);
+}
+
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::assembler() const
+-> const std::shared_ptr<Assembler>& {
+
+    return assembler_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::matrices() const
+-> const Matrices<Matrix,1>& {
+
+    return matrices_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::externalForce() const
+-> const ExternalForce& {
+
+    return externalForce_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {
+    boundaryConditions_.emplace_back(boundaryCondition);
+}
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryConditions(
+        const std::string& tag,
+        std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {
+
+    selectedConditions.resize(0);
+    for (size_t i=0; i<boundaryConditions_.size(); i++) {
+        if (boundaryConditions_[i]->tag() == tag)
+            selectedConditions.emplace_back(boundaryConditions_[i]);
+    }
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::boundaryConditions() const
+-> const std::vector<std::shared_ptr<BoundaryCondition>>& {
+
+    return boundaryConditions_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryPatchNodes(
+        const std::string& tag,
+        BoundaryPatchNodes& nodes) const {
+
+    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    nodes.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<nodes.size(); i++) {
+        nodes[i] = _boundaryConditions[i]->boundaryPatch()->getVertices();
+    }
+}
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryNodes(
+        const std::string& tag,
+        BoundaryNodes& nodes) const {
+
+    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    nodes.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<nodes.size(); i++) {
+        nodes[i] = _boundaryConditions[i]->boundaryNodes().get();
+    }
+}
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryFunctions(
+        const std::string& tag,
+        BoundaryFunctions& functions) const {
+
+    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    functions.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<functions.size(); i++) {
+        functions[i] = _boundaryConditions[i]->boundaryFunction().get();
+    }
+}
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryPatches(
+        const std::string& tag,
+        BoundaryPatches& patches) const {
+
+    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    patches.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<patches.size(); i++) {
+        patches[i] = _boundaryConditions[i]->boundaryPatch().get();
+    }
+}
+
+
+// ---------------------------
+// --- Body Implementation ---
+// ---------------------------
+
+// setter and getter
+
+template <class GridView>
+auto Body<GridView>::data() const
+-> const std::shared_ptr<BodyData<dim>>& {
+
+    return bodyData_;
+}
+
+
+template <class GridView>
+auto Body<GridView>::grid() const
+-> std::shared_ptr<GridType> {
+
+    return grid_;
+}
+
+template <class GridView>
+auto Body<GridView>::gridView() const
+-> const GridView& {
+
+    return gridView_;
+}
+
+template <class GridView>
+auto Body<GridView>::nVertices() const
+-> size_t {
+    return gridView_.size(dim);
+}
+
+
+template <class GridView>
+void Body<GridView>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {
+
+    boundaryConditions_.push_back(boundaryCondition);
+}
+
+
+template <class GridView>
+void Body<GridView>::boundaryConditions(
+        const std::string& tag,
+        std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {
+
+    selectedConditions.resize(0);
+    for (size_t i=0; i<boundaryConditions_.size(); i++) {
+        if (boundaryConditions_[i]->tag() == tag)
+            selectedConditions.push_back(boundaryConditions_[i]);
+    }
+}
+
+
+template <class GridView>
+auto Body<GridView>::boundaryConditions() const
+-> const std::vector<std::shared_ptr<BoundaryCondition>>& {
+
+    return boundaryConditions_;
+}
+
+
+template <class GridView>
+void Body<GridView>::boundaryPatchNodes(
+        const std::string& tag,
+        BoundaryPatchNodes& nodes) const {
+
+    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    nodes.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<nodes.size(); i++) {
+        nodes[i] = _boundaryConditions[i]->boundaryPatch()->getVertices();
+    }
+}
+
+template <class GridView>
+void Body<GridView>::boundaryNodes(
+        const std::string& tag,
+        BoundaryNodes& nodes) const {
+
+    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    nodes.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<nodes.size(); i++) {
+        nodes[i] = _boundaryConditions[i]->boundaryNodes().get();
+    }
+}
+
+template <class GridView>
+void Body<GridView>::boundaryFunctions(
+        const std::string& tag,
+        BoundaryFunctions& functions) const {
+
+    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    functions.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<functions.size(); i++) {
+        functions[i] = _boundaryConditions[i]->boundaryFunction().get();
+    }
+}
+
+template <class GridView>
+void Body<GridView>::boundaryPatches(
+        const std::string& tag,
+        BoundaryPatches& patches) const {
+
+    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    patches.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<patches.size(); i++) {
+        patches[i] = _boundaryConditions[i]->boundaryPatch().get();
+    }
+}
+
+#include "body_tmpl.cc"