dune/tectonic/data-structures/friction/frictioncouplingpair.hh

+#ifndef SRC_FRICTIONCOUPLINGPAIR_HH
+#define SRC_FRICTIONCOUPLINGPAIR_HH
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/contact/common/couplingpair.hh>
+
+#include "globalfrictiondata.hh"
+
+template <class GridType, class LocalVectorType, class field_type = double>
+class FrictionCouplingPair : public Dune::Contact::CouplingPair<GridType, GridType, field_type>{
+private:
+    static const int dim = GridType::dimensionworld;
+
+    using Base = Dune::Contact::CouplingPair<GridType,GridType,field_type>;
+    using LocalVector = LocalVectorType;
+
+    // friction data
+    std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch_;
+    std::shared_ptr<GlobalFrictionData<dim>> frictionData_;
+
+public:
+  void setWeakeningPatch(std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch) {
+      weakeningPatch_ = weakeningPatch;
+  }
+
+  void setFrictionData(std::shared_ptr<GlobalFrictionData<dim>> frictionData) {
+      frictionData_ = frictionData;
+  }
+
+  const auto& weakeningPatch() const {
+      return *weakeningPatch_;
+  }
+
+  const auto& frictionData() const {
+      return *frictionData_;
+  }
+};
+#endif

dune/tectonic/frictionpotential.hh → dune/tectonic/data-structures/friction/frictionpotential.hh

@@ -9,7 +9,7 @@
 #include <dune/common/exceptions.hh>
 #include <dune/common/function.hh>
 
-#include <dune/tectonic/frictiondata.hh>
+#include "frictiondata.hh"
 
 class FrictionPotential {
 public:
@@ -41,26 +41,42 @@ class TruncatedRateState : public FrictionPotential {
   }
 
   double differential(double V) const override {
+    //std::cout << "differential: " <<  weight * fd.C - weightedNormalStress * coefficientOfFriction(V) << std::endl;
     return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
   }
 
   double second_deriv(double V) const override {
+    //std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : -weightedNormalStress * (fd.a / V)) << std::endl;
+
     if (V <= Vmin)
-      return 0;
+      return 0.0;
 
     return -weightedNormalStress * (fd.a / V);
   }
 
   double regularity(double V) const override {
+    //std::cout << "regularity: " << ((std::abs(V - Vmin) < 1e-14) ? std::numeric_limits<double>::infinity() : std::abs(second_deriv(V))) << std::endl;
+
     if (std::abs(V - Vmin) < 1e-14) // TODO
       return std::numeric_limits<double>::infinity();
 
     return std::abs(second_deriv(V));
   }
 
+  double evaluate(double V) const override {
+
+    //std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1)) << std::endl;
+
+    if (V <= Vmin)
+        return 0.0;
+
+      return weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1);
+  }
+
   void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
+    //std::cout << "Vmin: " << Vmin << std::endl;
   }
 
 private:
@@ -92,6 +108,10 @@ class RegularisedRateState : public FrictionPotential {
     return std::abs(second_deriv(V));
   }
 
+  double evaluate(double V) const override {
+      return weight * fd.C * V - weightedNormalStress * 4 * fd.a * Vmin * std::pow(std::asinh(0.25 * V / Vmin), 2.0);
+  }
+
   void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
@@ -106,6 +126,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 → dune/tectonic/data-structures/friction/globalfriction.hh

@@ -9,9 +9,10 @@
 
 #include <dune/solvers/common/interval.hh>
 
-#include <dune/tectonic/localfriction.hh>
+#include "localfriction.hh"
 
-template <class Matrix, class Vector> class GlobalFriction {
+template <class Matrix, class Vector>
+class GlobalFriction {
 protected:
   using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
 
@@ -38,8 +39,9 @@ template <class Matrix, class Vector> class GlobalFriction {
   LocalNonlinearity const virtual &restriction(size_t i) const = 0;
 
   void addHessian(Vector const &v, Matrix &hessian) const {
-    for (size_t i = 0; i < v.size(); ++i)
+    for (size_t i = 0; i < v.size(); ++i) {
       restriction(i).addHessian(v[i], hessian[i][i]);
+    }
   }
 
   void directionalDomain(Vector const &, Vector const &,
@@ -81,6 +83,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 → dune/tectonic/data-structures/friction/globalfrictiondata.hh

@@ -4,7 +4,7 @@
 #include <dune/common/function.hh>
 #include <dune/common/fvector.hh>
 
-#include <dune/tectonic/frictiondata.hh>
+#include "frictiondata.hh"
 
 template <class DomainType>
 double evaluateScalarFunction(
@@ -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/data-structures/friction/globalratestatefriction.hh

+#ifndef DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
+#define DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
+
+#include <vector>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include "../../spatial-solving/contact/dualmortarcoupling.hh"
+
+#include "globalfrictiondata.hh"
+#include "globalfriction.hh"
+#include "frictioncouplingpair.hh"
+
+#include "../../utils/geocoordinate.hh"
+#include "../../utils/index-in-sorted-range.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 Base = GlobalFriction<Matrix, Vector>;
+
+  using field_type = typename Vector::field_type;
+  using typename Base::ScalarVector;
+  using typename Base::LocalVectorType;
+
+  using FrictionCoupling = FrictionCouplingPair<GridType, LocalVectorType, field_type>;
+  using ContactCoupling =  DualMortarCoupling<field_type, GridType>;
+
+  size_t bodyIndex(const size_t globalIdx) {
+     size_t i=0;
+
+     for (; i<offSet_.size()-1; i++) {
+         if (globalIdx < offSet_[i])
+             break;
+     }
+     return i;
+  }
+
+public:
+  GlobalRateStateFriction(const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+                          const std::vector<std::shared_ptr<FrictionCoupling>>& 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();
+    offSet_.resize(nBodies, 0);
+    for (size_t i=1; i<nBodies; i++) {
+        offSet_[i] = weights[i-1].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);
+    }
+
+    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(offSet_[bodyIdx] + i);
+                restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
+                                          couplings[j]->frictionData()(geoToPoint(it->geometry())));
+                break;
+            }
+        }
+    }
+  }
+
+  void updateAlpha(const std::vector<ScalarVector>& alpha) override {
+    for (size_t j = 0; j < restrictions_.size(); ++j) {
+      const auto globalDof = localToGlobal_[j];
+      const auto bodyIdx = bodyIndex(globalDof);
+
+      size_t bodyDof;
+      if (bodyIdx>0) {
+          bodyDof = globalDof - offSet_[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];
+  }
+
+private:
+  std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions_;
+  std::vector<size_t> offSet_; // index off-set by body id
+  std::vector<size_t> localToGlobal_;
+  WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction_;
+};
+#endif

dune/tectonic/localfriction.hh → dune/tectonic/data-structures/friction/localfriction.hh

@@ -10,7 +10,7 @@
 #include <dune/fufem/arithmetic.hh>
 #include <dune/solvers/common/interval.hh>
 
-#include <dune/tectonic/frictionpotential.hh>
+#include "frictionpotential.hh"
 
 template <size_t dimension> class LocalFriction {
 public:
@@ -19,6 +19,7 @@ template <size_t dimension> class LocalFriction {
   using VectorType = Dune::FieldVector<double, dimension>;
   using MatrixType = Dune::FieldMatrix<double, dimension, dimension>;
 
+  double virtual operator()(VectorType const &x) const = 0;
   void virtual updateAlpha(double alpha) = 0;
   double virtual regularity(VectorType const &x) const = 0;
   double virtual coefficientOfFriction(VectorType const &x) const = 0;
@@ -44,6 +45,10 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
   WrappedScalarFriction(Args... args)
       : func_(args...) {}
 
+  double operator()(VectorType const &x) const override {
+    return func_.evaluate(x.two_norm());
+  }
+
   void updateAlpha(double alpha) override { func_.updateAlpha(alpha); }
 
   double regularity(VectorType const &x) const override {
@@ -86,17 +91,26 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
       \f}
   */
   void addHessian(VectorType const &x, MatrixType &A) const override {
-    double const xnorm2 = x.two_norm2();
+    //std::cout << A << std::endl;
+    //std::cout << x << std::endl;
+
+ /*   double const xnorm2 = x.two_norm2();
     double const xnorm = std::sqrt(xnorm2);
     if (xnorm2 <= 0.0)
       return;
 
+    //std::cout << xnorm << " " << xnorm2 << std::endl;
+
     double const H1 = func_.differential(xnorm);
     double const H2 = func_.second_deriv(xnorm);
 
+    //std::cout << H1 << " " << H2 << std::endl;
+
     double const tensorweight = (H2 - H1 / xnorm) / xnorm2;
     double const idweight = H1 / xnorm;
 
+    //std::cout << tensorweight << " " << idweight << std::endl;
+
     for (size_t i = 0; i < dimension; ++i)
       for (size_t j = 0; j < i; ++j) {
         double const entry = tensorweight * x[i] * x[j];
@@ -107,7 +121,36 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
     for (size_t k = 0; k < dimension; ++k) {
       double const entry = tensorweight * x[k] * x[k];
       A[k][k] += entry + idweight;
+    }*/
+
+    const double xnorm = x.two_norm();
+    if (xnorm <= 0.0)
+      return;
+
+    VectorType y = x;
+    y /= xnorm;
+
+    double const H1 = func_.differential(xnorm);
+    double const H2 = func_.second_deriv(xnorm);
+
+    double const tensorweight = (H2 - H1 / xnorm);
+    double const idweight = H1 / xnorm;
+
+    //std::cout << tensorweight << " " << idweight << std::endl;
+
+    for (size_t i = 0; i < dimension; ++i)
+      for (size_t j = 0; j < i; ++j) {
+        double const entry = tensorweight * y[i] * y[j];
+        A[i][j] += entry;
+        A[j][i] += entry;
+      }
+
+    for (size_t k = 0; k < dimension; ++k) {
+      double const entry = tensorweight * y[k] * y[k];
+      A[k][k] += entry + idweight;
     }
+
+    //std::cout << A << std::endl;
   }
 
   void addGradient(VectorType const &x, VectorType &y) const override {

dune/tectonic/data-structures/globalfrictioncontainer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+/*
+#include "globalfrictioncontainer.hh"
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::operator[](std::size_t i)
+-> IndexObject& {
+    return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::operator[](std::size_t i) const
+-> const IndexObject& {
+    return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::size() const
+-> size_t {
+    return globalFriction_.size();
+}
+
+template <class BaseGlobalFriction, size_t depth>
+void GlobalFrictionContainer<BaseGlobalFriction, depth>::resize(std::list<size_t> list) {
+    assert(list.size() <= depth);
+
+    if (list.size() == 0) {
+        globalFriction_.resize(0);
+    } else {
+        globalFriction_.resize(list.front());
+        list.pop_front();
+
+        for (size_t i=0; i<size(); i++) {
+            globalFriction_[i].resize(list);
+        }
+    }
+}
+
+template <class BaseGlobalFriction, size_t depth>
+template <class VectorContainer>
+void GlobalFrictionContainer<BaseGlobalFriction, depth>::updateAlpha(const VectorContainer& newAlpha) {
+    assert(newAlpha.size() == size());
+
+    for (size_t i=0; i<size(); i++) {
+        globalFriction_[i].updateAlpha(newAlpha[i]);
+    }
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::globalFriction()
+-> GlobalFriction& {
+    return globalFriction_;
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::globalFriction() const
+-> const GlobalFriction& {
+    return globalFriction_;
+}
+
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::operator[](std::size_t i)
+-> IndexObject& {
+    return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::operator[](std::size_t i) const
+-> const IndexObject& {
+    return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::size() const
+-> size_t {
+    return globalFriction_.size();
+}
+
+template <class BaseGlobalFriction>
+void GlobalFrictionContainer<BaseGlobalFriction, 1>::resize(std::list<size_t> newSize) {
+    if (newSize.size() > 0) {
+        globalFriction_.resize(newSize.front(), nullptr);
+    } else {
+        globalFriction_.resize(0);
+    }
+}
+
+template <class BaseGlobalFriction>
+template <class Vector>
+void GlobalFrictionContainer<BaseGlobalFriction, 1>::updateAlpha(const Vector& newAlpha) {
+    for (size_t i=0; i<size(); i++) {
+        globalFriction_[i]->updateAlpha(newAlpha);
+    }
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::globalFriction()
+-> GlobalFriction& {
+    return globalFriction_;
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::globalFriction() const
+-> const GlobalFriction& {
+    return globalFriction_;
+}
+
+
+#include "globalfrictioncontainer_tmpl.cc"
+*/

dune/tectonic/data-structures/globalfrictioncontainer.hh

+#ifndef SRC_DATA_STRUCTURES_GLOBALFRICTIONCONTAINER_HH
+#define SRC_DATA_STRUCTURES_GLOBALFRICTIONCONTAINER_HH
+
+#include <vector>
+#include <cstddef>
+#include <list>
+#include <cassert>
+#include <memory>
+
+template <class BaseGlobalFriction, size_t depth>
+class GlobalFrictionContainer {
+public:
+    using IndexObject = GlobalFrictionContainer<BaseGlobalFriction, depth-1>;
+    using GlobalFriction = std::vector<IndexObject>;
+
+    GlobalFrictionContainer() {}
+
+    auto operator[](std::size_t i) -> IndexObject& {
+        return globalFriction_[i];
+    }
+
+    auto operator[](std::size_t i) const -> const IndexObject& {
+        return globalFriction_[i];
+    }
+
+    auto size() const -> size_t {
+        return globalFriction_.size();
+    }
+
+    void resize(std::list<size_t> list) {
+        assert(list.size() <= depth);
+
+        if (list.size() == 0) {
+            globalFriction_.resize(0);
+        } else {
+            globalFriction_.resize(list.front());
+            list.pop_front();
+
+            for (size_t i=0; i<size(); i++) {
+                globalFriction_[i].resize(list);
+            }
+        }
+    }
+
+    template <class VectorContainer>
+    void updateAlpha(const VectorContainer& newAlpha) {
+        assert(newAlpha.size() == size());
+
+        for (size_t i=0; i<size(); i++) {
+                globalFriction_[i].updateAlpha(newAlpha[i]);
+        }
+    }
+
+    auto globalFriction() -> GlobalFriction& {
+        return globalFriction_;
+    }
+
+    auto globalFriction() const -> const GlobalFriction& {
+        return globalFriction_;
+    }
+
+private:
+    GlobalFriction globalFriction_;
+};
+
+
+template <class BaseGlobalFriction>
+class GlobalFrictionContainer<BaseGlobalFriction, 1> {
+public:
+    using IndexObject = std::shared_ptr<BaseGlobalFriction>;
+    using GlobalFriction = std::vector<IndexObject>;
+
+    GlobalFrictionContainer() {}
+
+    auto operator[](std::size_t i) -> IndexObject& {
+        return globalFriction_[i];
+    }
+
+    auto operator[](std::size_t i) const -> const IndexObject& {
+        return globalFriction_[i];
+    }
+
+    auto size() const -> size_t {
+        return globalFriction_.size();
+    }
+
+    void resize(std::list<size_t> newSize) {
+        if (newSize.size() > 0) {
+            globalFriction_.resize(newSize.front(), nullptr);
+        } else {
+            globalFriction_.resize(0);
+        }
+    }
+
+    template <class Vector>
+    void updateAlpha(const Vector& newAlpha) {
+        for (size_t i=0; i<size(); i++) {
+                globalFriction_[i]->updateAlpha(newAlpha);
+        }
+    }
+
+    auto globalFriction() -> GlobalFriction& {
+        return globalFriction_;
+    }
+
+    auto globalFriction() const -> const GlobalFriction& {
+        return globalFriction_;
+    }
+
+private:
+    GlobalFriction globalFriction_;
+};
+
+#endif

dune/tectonic/data-structures/matrices.hh

+#ifndef SRC_MATRICES_HH
+#define SRC_MATRICES_HH
+
+template <class Matrix, size_t n>
+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() {
+    elasticity.resize(n);
+    damping.resize(n);
+    mass.resize(n);
+
+      for (size_t i=0; i<n; i++) {
+        elasticity[i] = std::make_shared<Matrix>();
+        damping[i] = std::make_shared<Matrix>();
+        mass[i] = std::make_shared<Matrix>();
+      }
+  }
+};
+
+template <class Matrix> 
+class Matrices<Matrix, 1> {
+public:
+  std::shared_ptr<Matrix> elasticity;
+  std::shared_ptr<Matrix> damping;
+  std::shared_ptr<Matrix> mass;
+
+  Matrices() :
+    elasticity(std::make_shared<Matrix>()),
+    damping(std::make_shared<Matrix>()),
+    mass(std::make_shared<Matrix>())
+  {}
+};
+#endif

dune/tectonic/data-structures/network/CMakeLists.txt

+add_custom_target(tectonic_dune_data-structures_network SOURCES
+  contactnetwork.hh
+  contactnetwork.cc
+  levelcontactnetwork.hh
+  levelcontactnetwork.cc
+)
+
+#install headers
+install(FILES
+  contactnetwork.hh
+  levelcontactnetwork.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/data-structures/network/contactnetwork.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include "../friction/globalratestatefriction.hh"
+#include "../friction/frictionpotential.hh"
+
+#include "contactnetwork.hh"
+
+#include "../../assemblers.hh"
+#include "../../utils/tobool.hh"
+
+
+template <class HostGridType, class VectorType>
+ContactNetwork<HostGridType, VectorType>::ContactNetwork(
+        size_t nBodies,
+        size_t nCouplings,
+        size_t nLevels) :
+    levelContactNetworks_(nLevels),
+    bodies_(nBodies),
+    couplings_(nCouplings),
+    frictionBoundaries_(nBodies),
+    stateEnergyNorms_(nBodies),
+    nBodyAssembler_(nBodies, nCouplings, 0.0)
+{}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::addLevel(std::shared_ptr<LevelContactNetwork> levelContactNetwork) {
+    size_t level = levelContactNetwork->level();
+    if (level >= levelContactNetworks_.size()) {
+        levelContactNetworks_.resize(level);
+    }
+
+    levelContactNetworks_[level] = levelContactNetwork;
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::addLevel(const std::vector<size_t>& bodyLevels, const size_t level) {
+    assert(bodyLevels.size() == nBodies());
+
+    if (level >= levelContactNetworks_.size()) {
+        levelContactNetworks_.resize(level+1);
+    }
+
+    auto& levelContactNetwork = levelContactNetworks_[level];
+    levelContactNetwork = std::make_shared<LevelContactNetwork>(nBodies(), nCouplings(), level);
+
+    std::vector<std::shared_ptr<Body>> bodies(nBodies());
+
+    for (size_t id=0; id<nBodies(); id++) {
+        const auto& body = this->body(id);
+        bodies[id] = std::make_shared<Body>(body->data(), body->grid(), std::min((size_t) body->grid()->maxLevel(), bodyLevels[id]));
+    }
+    levelContactNetwork->setBodies(bodies);
+    levelContactNetwork->setCouplings(this->couplings_);
+    levelContactNetwork->build();
+}
+
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::assemble() {
+    std::vector<const GridType*> grids(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        grids[i] = bodies_[i]->grid().get();
+        bodies_[i]->assemble();
+        frictionBoundaries_[i] = std::make_unique<BoundaryPatch>(bodies_[i]->gridView(), false);
+    }
+
+    // set up dune-contact nBodyAssembler
+    nBodyAssembler_.setGrids(grids);
+
+    for (size_t i=0; i<nCouplings(); i++) {
+      nBodyAssembler_.setCoupling(*couplings_[i], i);
+    }
+
+    std::vector<std::shared_ptr<Dune::BitSetVector<1>>> dirichletVertices(nBodies());
+    for (size_t i=0; i<nBodies(); i++) {
+        const auto& body = this->body(i);
+        std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
+        body->boundaryConditions("dirichlet", boundaryConditions);
+
+        dirichletVertices[i] = std::make_shared<Dune::BitSetVector<1>>(body->nVertices());
+        auto& vertices = *dirichletVertices[i];
+        vertices.unsetAll();
+
+        if (boundaryConditions.size()<1)
+            continue;
+
+        for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
+           const auto& boundaryVertices = *boundaryConditions[bc]->boundaryPatch()->getVertices();
+           for (size_t j=0; j<boundaryVertices.size(); j++) {
+               vertices[j][0] = vertices[j][0] or boundaryVertices[j][0];
+           }
+        }
+    }
+
+    nBodyAssembler_.setDirichletVertices(dirichletVertices);
+    nBodyAssembler_.assembleTransferOperator();
+    nBodyAssembler_.assembleObstacle();
+
+    for (size_t i=0; i<nCouplings(); i++) {
+      auto& coupling = couplings_[i];
+      const auto& contactCoupling = nBodyAssembler_.getContactCouplings()[i]; // dual mortar object holding boundary patches
+      const auto nonmortarIdx = coupling->gridIdx_[0];
+      //const auto mortarIdx = coupling->gridIdx_[1];
+
+      frictionBoundaries_[nonmortarIdx]->addPatch(contactCoupling->nonmortarBoundary());
+      //frictionBoundaries_[mortarIdx]->addPatch(contactCoupling->mortarBoundary());
+    }
+
+    // assemble state energy norm
+    frictionBoundaryMass_.resize(nBodies());
+    for (size_t i=0; i<nBodies(); i++) {
+        frictionBoundaryMass_[i] = std::make_unique<ScalarMatrix>();
+        bodies_[i]->assembler()->assembleFrictionalBoundaryMass(*frictionBoundaries_[i], *frictionBoundaryMass_[i]);
+        *frictionBoundaryMass_[i] /= frictionBoundaries_[i]->area(); // TODO: weight by individual friction patches?
+
+        stateEnergyNorms_[i] = std::make_unique<typename LeafBody::StateEnergyNorm>(*frictionBoundaryMass_[i]);
+    }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::assembleFriction(
+        const Config::FrictionModel& frictionModel,
+        const std::vector<ScalarVector>& weightedNormalStress) {
+
+      assert(weightedNormalStress.size() == bodies_.size());
+
+      const size_t nBodies = bodies_.size();
+
+      // Lumping of the nonlinearity
+      std::vector<ScalarVector> weights(nBodies);
+      {
+        NeumannBoundaryAssembler<GridType, typename ScalarVector::block_type>
+            frictionalBoundaryAssembler(std::make_shared<
+                ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
+                1));
+
+        for (size_t i=0; i<nBodies; i++) {
+            bodies_[i]->assembler()->assembleBoundaryFunctional(frictionalBoundaryAssembler, weights[i], *frictionBoundaries_[i]);
+        }
+      }
+
+    /*  globalFriction_ = std::make_shared<GlobalRateStateFriction<
+          Matrix, Vector, ZeroFunction, DeformedGridType>>(
+          nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress); */
+
+      switch (frictionModel) {
+        case Config::Truncated:
+          globalFriction_ = std::make_shared<GlobalRateStateFriction<
+              Matrix, VectorType, TruncatedRateState, GridType>>(
+              nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
+        case Config::Regularised:
+          globalFriction_ = std::make_shared<GlobalRateStateFriction<
+              Matrix, VectorType, RegularisedRateState, GridType>>(
+              nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
+        default:
+          assert(false);
+          break;
+      }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::setBodies(const std::vector<std::shared_ptr<LeafBody>> bodies) {
+    assert(nBodies()==bodies.size());
+    bodies_ = bodies;
+
+    matrices_.elasticity.resize(nBodies());
+    matrices_.damping.resize(nBodies());
+    matrices_.mass.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        matrices_.elasticity[i] = bodies_[i]->matrices().elasticity;
+        matrices_.damping[i] = bodies_[i]->matrices().damping;
+        matrices_.mass[i] = bodies_[i]->matrices().mass;
+    }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings) {
+    assert(this->nCouplings()==couplings.size());
+    couplings_ = couplings;
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::constructBody(
+        const std::shared_ptr<BodyData<dim>>& bodyData,
+        const std::shared_ptr<HostGridType>& grid,
+        std::shared_ptr<LeafBody>& body) const {
+
+    body = std::make_shared<LeafBody>(bodyData, grid);
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::setDeformation(const std::vector<VectorType>& totalDeformation) {
+    assert(totalDeformation.size() == nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        body(i)->setDeformation(totalDeformation[i]);
+    }
+
+    nBodyAssembler_.assembleTransferOperator();
+    nBodyAssembler_.assembleObstacle();
+
+    for (size_t i=0; i<levelContactNetworks_.size(); i++) {
+        levelContactNetworks_[i]->build();
+    }
+}
+
+template <class HostGridType, class VectorType>
+template <class LevelBoundaryPatch>
+void ContactNetwork<HostGridType, VectorType>::constructCoupling(
+        int nonMortarBodyIdx, int mortarBodyIdx,
+        const std::shared_ptr<LevelBoundaryPatch>& nonmortarPatch,
+        const std::shared_ptr<LevelBoundaryPatch>& mortarPatch,
+        const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
+        const std::shared_ptr<GlobalFrictionData<dim>>& frictionData,
+        std::shared_ptr<FrictionCouplingPair>& coupling) const {
+
+    coupling = std::make_shared<FrictionCouplingPair>();
+
+    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<BoundaryPatch>>();
+    std::shared_ptr<typename FrictionCouplingPair::BackEndType> backend = nullptr;
+    coupling->set(nonMortarBodyIdx, mortarBodyIdx, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
+    coupling->setWeakeningPatch(weakeningPatch);
+    coupling->setFrictionData(frictionData);
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::level(size_t level) const
+-> const std::shared_ptr<LevelContactNetwork>& {
+
+    return levelContactNetworks_[level];
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nLevels() const
+-> size_t {
+
+    return levelContactNetworks_.size();
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodies() const
+-> size_t {
+
+    return bodies_.size();
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nCouplings() const
+-> size_t {
+
+    return couplings_.size();
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::body(int i) const
+-> const std::shared_ptr<LeafBody>& {
+
+    return bodies_[i];
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::coupling(int i) const
+-> const std::shared_ptr<FrictionCouplingPair>& {
+
+    return couplings_[i];
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::couplings() const
+-> const std::vector<std::shared_ptr<FrictionCouplingPair>>& {
+
+    return couplings_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler()
+-> NBodyAssembler& {
+
+    return nBodyAssembler_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler() const
+-> const NBodyAssembler& {
+
+    return nBodyAssembler_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::matrices() const
+-> const Matrices<Matrix,2>& {
+
+    return matrices_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::stateEnergyNorms() const
+-> const StateEnergyNorms& {
+
+    return stateEnergyNorms_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::globalFriction() const
+-> GlobalFriction& {
+
+    return *globalFriction_;
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::externalForces(ExternalForces& externalForces) const {
+    externalForces.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        externalForces[i] = std::make_unique<typename LeafBody::ExternalForce>(bodies_[i]->externalForce());
+    }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::totalNodes(
+        const std::string& tag,
+        Dune::BitSetVector<dim>& totalNodes) const {
+
+    int totalSize = 0;
+    for (size_t i=0; i<nBodies(); i++) {
+        totalSize += this->body(i)->nVertices();
+    }
+
+    totalNodes.resize(totalSize);
+
+    int idx=0;
+    for (size_t i=0; i<nBodies(); i++) {
+        const auto& body = this->body(i);
+        std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
+        body->boundaryConditions(tag, boundaryConditions);
+
+        if (boundaryConditions.size()>0) {
+            const int idxBackup = idx;
+            for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
+                const auto& nodes = boundaryConditions[bc]->boundaryNodes();
+                for (size_t j=0; j<nodes->size(); j++, idx++)
+                    for (int k=0; k<dim; k++)
+                        totalNodes[idx][k] = totalNodes[idx][k] or (*nodes)[j][k];
+                idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+            }
+        } else {
+            idx += body->nVertices();
+        }
+    }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::boundaryPatchNodes(
+        const std::string& tag,
+        BoundaryPatchNodes& nodes) const {
+
+    nodes.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryPatchNodes(tag, nodes[i]);
+    }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::boundaryNodes(
+        const std::string& tag,
+        BoundaryNodes& nodes) const {
+
+    nodes.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryNodes(tag, nodes[i]);
+    }
+}
+
+// collects all leaf boundary functions from the different bodies identified by "tag"
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::boundaryFunctions(
+        const std::string& tag,
+        BoundaryFunctions& functions) const {
+
+    functions.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryFunctions(tag, functions[i]);
+    }
+}
+
+// collects all leaf boundary patches from the different bodies identified by "tag"
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::boundaryPatches(
+        const std::string& tag,
+        BoundaryPatches& patches) const {
+
+    patches.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryPatches(tag, patches[i]);
+    }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const {
+    nodes.resize(nBodies());
+    for (size_t i=0; i<nBodies(); i++) {
+        nodes[i] = frictionBoundaries_[i]->getVertices();
+    }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::frictionPatches(std::vector<const BoundaryPatch*>& patches) const {
+    patches.resize(nBodies());
+    for (size_t i=0; i<nBodies(); i++) {
+        patches[i] = frictionBoundaries_[i].get();
+    }
+}
+
+#include "contactnetwork_tmpl.cc"

dune/tectonic/data-structures/network/contactnetwork.hh

+#ifndef SRC_CONTACTNETWORK_HH
+#define SRC_CONTACTNETWORK_HH
+
+#include <memory>
+
+#include <dune/common/promotiontraits.hh>
+
+#include "../../spatial-solving/contact/nbodyassembler.hh"
+#include "../../assemblers.hh"
+#include "../enums.hh"
+#include "../matrices.hh"
+
+#include "../body/body.hh"
+#include "../body/bodydata.hh"
+#include "../friction/frictioncouplingpair.hh"
+#include "../friction/globalfriction.hh"
+#include "../friction/globalfrictiondata.hh"
+
+#include "levelcontactnetwork.hh"
+
+template <class HostGridType, class VectorType>
+class ContactNetwork {
+    //using Vector = Dune::BlockVector<Dune::FieldVector<double, dim>>;
+public:
+    enum {dim = HostGridType::dimensionworld};
+
+    using field_type = typename Dune::PromotionTraits<typename VectorType::field_type,
+                                                typename HostGridType::ctype>::PromotedType;
+
+    using LeafBody = LeafBody<HostGridType, VectorType>;
+    using GridType = typename LeafBody::GridType;
+    using GridView = typename LeafBody::GridView;
+    using BoundaryPatch = BoundaryPatch<GridView>;
+
+    using Assembler = typename LeafBody::Assembler;
+    using Matrix = typename Assembler::Matrix;
+
+    using ScalarMatrix = typename Assembler::ScalarMatrix;
+    using ScalarVector = typename Assembler::ScalarVector;
+
+    using LocalVector = typename VectorType::block_type;
+    using FrictionCouplingPair = FrictionCouplingPair<GridType, LocalVector, field_type>;
+
+    using Function = Dune::VirtualFunction<double, double>;
+
+    using BoundaryFunctions = std::vector<typename LeafBody::BoundaryFunctions>;
+    using BoundaryNodes = std::vector<typename LeafBody::BoundaryNodes>;
+    using BoundaryPatchNodes = std::vector<typename LeafBody::BoundaryPatchNodes>;
+    using BoundaryPatches = std::vector<typename LeafBody::BoundaryPatches>;
+
+    using GlobalFriction = GlobalFriction<Matrix, VectorType>;
+
+    using StateEnergyNorms = std::vector<std::unique_ptr<typename LeafBody::StateEnergyNorm> >;
+    using ExternalForces = std::vector<std::unique_ptr<const typename LeafBody::ExternalForce> >;
+
+    using NBodyAssembler = NBodyAssembler<GridType, VectorType>;
+
+    using LevelContactNetwork = LevelContactNetwork<GridType, FrictionCouplingPair, field_type>;
+    using Body = typename LevelContactNetwork::Body;
+
+private:
+    std::vector<std::shared_ptr<LevelContactNetwork>> levelContactNetworks_;
+
+    std::vector<std::shared_ptr<LeafBody>> bodies_;
+    std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
+
+    std::vector<std::unique_ptr<BoundaryPatch>> frictionBoundaries_; // union of all boundary patches per body
+    std::vector<std::unique_ptr<ScalarMatrix>> frictionBoundaryMass_;
+    StateEnergyNorms stateEnergyNorms_;
+
+    NBodyAssembler nBodyAssembler_;
+
+    Matrices<Matrix,2> matrices_;
+
+    std::shared_ptr<GlobalFriction> globalFriction_;
+
+public:
+    ContactNetwork(size_t nBodies, size_t nCouplings, size_t nLevels = 0);
+
+    void addLevel(std::shared_ptr<LevelContactNetwork> levelContactNetwork);
+    void addLevel(const std::vector<size_t>& bodyLevels, const size_t level);
+
+    void assemble();
+    void assembleFriction(
+            const Config::FrictionModel& frictionModel,
+            const std::vector<ScalarVector>& weightedNormalStress);
+
+    void setBodies(const std::vector<std::shared_ptr<LeafBody>> bodies);
+    void setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings);
+
+    void setDeformation(const std::vector<VectorType>& totalDeformation);
+
+    void constructBody(
+            const std::shared_ptr<BodyData<dim>>& bodyData,
+            const std::shared_ptr<HostGridType>& grid,
+            std::shared_ptr<LeafBody>& body) const;
+
+    template <class LevelBoundaryPatch>
+    void constructCoupling(
+            int nonMortarBodyIdx,
+            int mortarBodyIdx,
+            const std::shared_ptr<LevelBoundaryPatch>& nonmortarPatch,
+            const std::shared_ptr<LevelBoundaryPatch>& mortarPatch,
+            const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
+            const std::shared_ptr<GlobalFrictionData<dim>>& frictionData,
+            std::shared_ptr<FrictionCouplingPair>& coupling) const;
+
+    // getter
+    auto level(size_t level) const -> const std::shared_ptr<LevelContactNetwork>&;
+
+    auto nLevels() const -> size_t;
+    auto nBodies() const -> size_t;
+    auto nCouplings() const -> size_t;
+
+    auto body(int i) const -> const std::shared_ptr<LeafBody>&;
+
+    auto coupling(int i) const -> const std::shared_ptr<FrictionCouplingPair>&;
+    auto couplings() const -> const std::vector<std::shared_ptr<FrictionCouplingPair>>&;
+
+    auto nBodyAssembler() -> NBodyAssembler&;
+    auto nBodyAssembler() const -> const NBodyAssembler&;
+
+    auto matrices() const -> const Matrices<Matrix,2>&;
+
+    auto stateEnergyNorms() const -> const StateEnergyNorms&;
+    auto globalFriction() const -> GlobalFriction&;
+    void externalForces(ExternalForces& externalForces) const;
+
+    void totalNodes(
+            const std::string& tag,
+            Dune::BitSetVector<dim>& totalNodes) const;
+    void boundaryPatches(
+            const std::string& tag,
+            BoundaryPatches& patches) const;
+    void boundaryPatchNodes(
+            const std::string& tag,
+            BoundaryPatchNodes& nodes) const;
+    void boundaryNodes(
+            const std::string& tag,
+            BoundaryNodes& nodes) const;
+    void boundaryFunctions(
+            const std::string& tag,
+            BoundaryFunctions& functions) const;
+    void frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const;
+    void frictionPatches(std::vector<const BoundaryPatch*>& patches) const;
+};
+#endif

dune/tectonic/data-structures/network/contactnetwork_tmpl.cc

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

dune/tectonic/data-structures/network/levelcontactnetwork.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/facehierarchy.hh>
+
+#include "../../assemblers.hh"
+#include "../../utils/tobool.hh"
+#include "../../utils/debugutils.hh"
+
+#include "../friction/globalratestatefriction.hh"
+#include "../friction/frictionpotential.hh"
+
+#include "levelcontactnetwork.hh"
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::LevelContactNetwork(
+        int nBodies,
+        int nCouplings,
+        int level) :
+    level_(level),
+    bodies_(nBodies),
+    couplings_(nCouplings),
+    nonmortarPatches_(nCouplings),
+    mortarPatches_(nCouplings),
+    glues_(nCouplings)
+{}
+
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::setBodies(const std::vector<std::shared_ptr<Body>> bodies) {
+    assert(nBodies() == bodies.size());
+    bodies_ = bodies;
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings) {
+    assert(nCouplings() == couplings.size());
+    couplings_ = couplings;
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::constructBody(
+        const std::shared_ptr<BodyData<dimworld>>& bodyData,
+        const std::shared_ptr<GridType>& grid,
+        const size_t level,
+        std::shared_ptr<Body>& body) const {
+
+    body = std::make_shared<Body>(bodyData, grid, level);
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::totalNodes(
+        const std::string& tag,
+        Dune::BitSetVector<dimworld>& totalNodes) const {
+
+    int totalSize = 0;
+    for (size_t i=0; i<nBodies(); i++) {
+        totalSize += this->body(i)->nVertices();
+    }
+
+    totalNodes.resize(totalSize);
+
+    int idx=0;
+    for (size_t i=0; i<nBodies(); i++) {
+        const auto& body = this->body(i);
+        std::vector<std::shared_ptr<typename Body::BoundaryCondition>> boundaryConditions;
+        body->boundaryConditions(tag, boundaryConditions);
+
+        const int idxBackup = idx;
+        for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
+            const auto& nodes = boundaryConditions[bc]->boundaryNodes();
+            for (size_t j=0; j<nodes->size(); j++, idx++)
+                for (int k=0; k<dimworld; k++)
+                    totalNodes[idx][k] = (*nodes)[j][k];
+            idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+        }
+    }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::boundaryPatchNodes(
+        const std::string& tag,
+        BoundaryPatchNodes& nodes) const {
+
+    nodes.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryPatchNodes(tag, nodes[i]);
+    }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::boundaryNodes(
+        const std::string& tag,
+        BoundaryNodes& nodes) const {
+
+    nodes.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryNodes(tag, nodes[i]);
+    }
+}
+
+// collects all leaf boundary functions from the different bodies identified by "tag"
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::boundaryFunctions(
+        const std::string& tag,
+        BoundaryFunctions& functions) const {
+
+    functions.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryFunctions(tag, functions[i]);
+    }
+}
+
+// collects all leaf boundary patches from the different bodies identified by "tag"
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::boundaryPatches(
+        const std::string& tag,
+        BoundaryPatches& patches) const {
+
+    patches.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryPatches(tag, patches[i]);
+    }
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::level() const
+-> int {
+
+    return level_;
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::nBodies() const
+-> size_t {
+
+    return bodies_.size();
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::nCouplings() const
+-> size_t {
+
+    return couplings_.size();
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::body(int i) const
+-> const std::shared_ptr<Body>& {
+
+    return bodies_[i];
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::coupling(int i) const
+-> const std::shared_ptr<FrictionCouplingPair>& {
+
+    return couplings_[i];
+}
+
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::couplings() const
+-> const std::vector<std::shared_ptr<FrictionCouplingPair>>& {
+
+    return couplings_;
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::prolong(const BoundaryPatch& coarsePatch, BoundaryPatch& finePatch, const size_t fineLevel) {
+
+    if (not(coarsePatch.isInitialized()))
+        DUNE_THROW(Dune::Exception, "Coarse boundary patch has not been set up correctly!");
+
+    const GridType& grid = coarsePatch.gridView().grid();
+
+    // Set array sizes correctly
+    finePatch.setup(grid.levelGridView(fineLevel));
+
+    for (const auto& pIt : coarsePatch) {
+        const auto& inside = pIt.inside();
+
+        if (inside.level() == fineLevel)
+            finePatch.addFace(inside, pIt.indexInInside());
+        else {
+            Face<GridType> face(grid, inside, pIt.indexInInside());
+
+            typename Face<GridType>::HierarchicIterator it = face.hbegin(fineLevel);
+            typename Face<GridType>::HierarchicIterator end = face.hend(fineLevel);
+
+            for(; it!=end; ++it) {
+                if (it->element().level() == fineLevel)
+                    finePatch.addFace(it->element(), it->index());
+            }
+        }
+    }
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::build(field_type overlap) {
+
+    for (size_t i=0; i<couplings_.size(); i++) {
+        auto& coupling = *couplings_[i];
+
+        const size_t nmBody = coupling.gridIdx_[0];
+        const size_t mBody = coupling.gridIdx_[1];
+
+        this->prolong(*coupling.patch0(), nonmortarPatches_[i], body(nmBody)->level());
+        this->prolong(*coupling.patch1(), mortarPatches_[i], body(mBody)->level());
+
+        using Element = typename GridView::template Codim<0>::Entity;
+        auto desc0 = [&] (const Element& e, unsigned int face) {
+            return nonmortarPatches_[i].contains(e,face);
+        };
+        auto desc1 = [&] (const Element& e, unsigned int face) {
+            return mortarPatches_[i].contains(e,face);
+        };
+
+        auto extract0 = std::make_shared<Extractor>(body(nmBody)->gridView(), desc0);
+        auto extract1 = std::make_shared<Extractor>(body(mBody)->gridView(), desc1);
+
+        gridGlueBackend_ = std::make_shared< Dune::GridGlue::ContactMerge<dim, ctype> >(overlap);
+        glues_[i] = std::make_shared<Glue>(extract0, extract1, gridGlueBackend_);
+        glues_[i]->build();
+    }
+}
+
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::glue(int i) const
+-> const std::shared_ptr<Glue>& {
+
+    return glues_[i];
+}
+
+#include "levelcontactnetwork_tmpl.cc"