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<maxIndex_.size(); i++) {
+         if (globalIdx < maxIndex_[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();
+
+    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;
+    }
+  }
+
+  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];
+  }
+
+private:
+  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/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 {
+    //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];
@@ -108,6 +122,8 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
       double const entry = tensorweight * x[k] * x[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] = (*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();
+    }
+}
+
+#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;
+};
+#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"

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

+#ifndef SRC_DATA_STRUCTURES_LEVELCONTACTNETWORK_HH
+#define SRC_DATA_STRUCTURES_LEVELCONTACTNETWORK_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/istl/bvector.hh>
+
+#include <dune/grid-glue/gridglue.hh>
+#include <dune/grid-glue/extractors/codim1extractor.hh>
+#include <dune/grid-glue/merging/merger.hh>
+
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/assemblers/nbodyassembler.hh>
+//#include <dune/contact/assemblers/dualmortarcoupling.hh>
+#include <dune/contact/projections/normalprojection.hh>
+
+#include <dune/fufem/boundarypatch.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"
+
+template <class GridTypeTEMPLATE, class FrictionCouplingPair, class field_type>
+class LevelContactNetwork {
+public:
+    using GridType = GridTypeTEMPLATE;
+
+    enum {dim = GridType::dimension};
+    enum {dimworld = GridType::dimensionworld};
+
+    using ctype = typename GridType::ctype;
+
+    using GridView = typename GridType::LevelGridView;
+    using Body = Body<GridView>;
+
+    using Function = Dune::VirtualFunction<double, double>;
+
+    using BoundaryFunctions = std::vector<typename Body::BoundaryFunctions>;
+    using BoundaryNodes = std::vector<typename Body::BoundaryNodes>;
+    using BoundaryPatchNodes = std::vector<typename Body::BoundaryPatchNodes>;
+    using BoundaryPatches = std::vector<typename Body::BoundaryPatches>;
+    using BoundaryPatch = BoundaryPatch<GridView>;
+
+    using Extractor = Dune::GridGlue::Codim1Extractor<GridView>;
+    using Glue = Dune::GridGlue::GridGlue<Extractor, Extractor>;
+
+public:
+    LevelContactNetwork(int nBodies, int nCouplings, int level = 0);
+
+    void setBodies(const std::vector<std::shared_ptr<Body>> bodies);
+    void setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings);
+
+    void constructBody(
+            const std::shared_ptr<BodyData<dimworld>>& bodyData,
+            const std::shared_ptr<GridType>& grid,
+            const size_t level,
+            std::shared_ptr<Body>& body) const;
+
+    void build(field_type overlap = 1e-2);
+
+    // getter
+    void totalNodes(
+            const std::string& tag,
+            Dune::BitSetVector<dimworld>& 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;
+
+    auto level() const -> int;
+
+    auto nBodies() const -> size_t;
+    auto nCouplings() const -> size_t;
+
+    auto body(int i) const -> const std::shared_ptr<Body>&;
+
+    auto coupling(int i) const -> const std::shared_ptr<FrictionCouplingPair>&;
+    auto couplings() const -> const std::vector<std::shared_ptr<FrictionCouplingPair>>&;
+
+    auto glue(int i) const -> const std::shared_ptr<Glue>&;
+
+private:
+    void prolong(const BoundaryPatch& coarsePatch, BoundaryPatch& finePatch, const size_t fineLevel);
+
+    const int level_;
+
+    std::vector<std::shared_ptr<Body>> bodies_;
+    std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
+
+    std::vector<BoundaryPatch> nonmortarPatches_;
+    std::vector<BoundaryPatch> mortarPatches_;
+    std::shared_ptr< Dune::GridGlue::Merger<field_type, dim-1, dim-1, dim> > gridGlueBackend_;
+    std::vector<std::shared_ptr<Glue>> glues_;
+};
+#endif

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

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "../../explicitvectors.hh"
+
+#include "../friction/frictioncouplingpair.hh"
+#include "contactnetwork.hh"
+#include "levelcontactnetwork.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+template class LevelContactNetwork<typename MyContactNetwork::GridType, typename MyContactNetwork::FrictionCouplingPair, typename MyContactNetwork::field_type>;

dune/tectonic/data-structures/program_state.hh

+#ifndef SRC_PROGRAM_STATE_HH
+#define SRC_PROGRAM_STATE_HH
+
+#include <dune/common/parametertree.hh>
+
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+
+#include "body/bodydata.hh"
+
+#include "../assemblers.hh"
+#include "network/contactnetwork.hh"
+#include "matrices.hh"
+#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
+#include "../spatial-solving/solverfactory.hh"
+#include "../spatial-solving/solverfactory.cc"
+#include "../spatial-solving/tnnmg/functional.hh"
+#include "../spatial-solving/tnnmg/zerononlinearity.hh"
+#include "../utils/debugutils.hh"
+
+
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include "../spatial-solving/contact/nbodycontacttransfer.hh"
+
+template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class BodyState {
+public:
+  using Vector = VectorTEMPLATE;
+  using ScalarVector = ScalarVectorTEMPLATE;
+
+  BodyState(Vector * _u, Vector * _v, Vector * _a, ScalarVector * _alpha, ScalarVector * _weightedNormalStress)
+    : u(_u),
+      v(_v),
+      a(_a),
+      alpha(_alpha),
+      weightedNormalStress(_weightedNormalStress) {}
+
+public:
+  Vector * const u;
+  Vector * const v;
+  Vector * const a;
+  ScalarVector * const alpha;
+  ScalarVector * const weightedNormalStress;
+};
+
+
+template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
+public:
+  using Vector = VectorTEMPLATE;
+  using ScalarVector = ScalarVectorTEMPLATE;
+  using BodyState = BodyState<Vector, ScalarVector>;
+
+private:
+    using LocalVector = typename Vector::block_type;
+    //using LocalMatrix = typename Matrix::block_type;
+    const static int dims = LocalVector::dimension;
+    using BitVector = Dune::BitSetVector<dims>;
+
+public:
+  ProgramState(const std::vector<size_t>& leafVertexCounts)
+    : bodyCount_(leafVertexCounts.size()),
+      bodyStates(bodyCount_),
+      u(bodyCount_),
+      v(bodyCount_),
+      a(bodyCount_),
+      alpha(bodyCount_),
+      weightedNormalStress(bodyCount_) {
+    for (size_t i=0; i<bodyCount_; i++) {
+      size_t leafVertexCount = leafVertexCounts[i];
+
+      u[i].resize(leafVertexCount),
+      v[i].resize(leafVertexCount),
+      a[i].resize(leafVertexCount),
+      alpha[i].resize(leafVertexCount),
+      weightedNormalStress[i].resize(leafVertexCount),
+
+      bodyStates[i] = new BodyState(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]);
+    }
+  }
+
+  ~ProgramState() {
+    for (size_t i=0; i<bodyCount_; i++) {
+      delete bodyStates[i];
+    }
+  }
+
+  size_t size() const {
+      return bodyCount_;
+  }
+
+  // Set up initial conditions
+  template <class ContactNetwork>
+  void setupInitialConditions(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) {
+    using Matrix = typename ContactNetwork::Matrix;
+    const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+
+    std::cout << "-- setupInitialConditions --" << std::endl;
+    std::cout << "----------------------------" << std::endl;
+
+    // make linear solver for linear correction in TNNMGStep
+    using Norm =  EnergyNorm<Matrix, Vector>;
+    using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+    using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+ /*   const auto& preconditionerParset = parset.sub("solver.tnnmg.preconditioner");
+
+    Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+    Preconditioner preconditioner(preconditionerParset, contactNetwork, activeLevels);
+    preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
+
+    std::cout << "Building preconditioner..." << std::endl;
+    preconditioner.build();
+
+    auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
+    cgStep->setPreconditioner(preconditioner);
+
+    Norm norm(*cgStep);
+
+    auto linearSolver = std::make_shared<LinearSolver>(cgStep, parset.get<int>("solver.tnnmg.main.multi"),
+                                                       parset.get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"),
+                                                       norm, Solver::QUIET);
+
+*/
+    // set multigrid solver
+    auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
+
+    using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+    using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+    TransferOperators transfer(contactNetwork.nLevels()-1);
+    for (size_t i=0; i<transfer.size(); i++) {
+        transfer[i] = std::make_shared<TransferOperator>();
+        transfer[i]->setup(contactNetwork, i, i+1);
+    }
+
+    // Remove any recompute filed so that initially the full transferoperator is assembled
+    for (size_t i=0; i<transfer.size(); i++)
+        std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+    auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+    linearMultigridStep->setMGType(1, 3, 3);
+    linearMultigridStep->setSmoother(smoother);
+    linearMultigridStep->setTransferOperators(transfer);
+
+    Norm norm(*linearMultigridStep);
+
+    auto linearSolver = std::make_shared<LinearSolver>(linearMultigridStep, parset.get<int>("solver.tnnmg.main.multi"), parset.get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+
+    // Solving a linear problem with a multigrid solver
+    auto const solveLinearProblem = [&](
+        const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
+        const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
+        Dune::ParameterTree const &_localParset) {
+
+      std::vector<const Matrix*> matrices_ptr(_matrices.size());
+      for (size_t i=0; i<matrices_ptr.size(); i++) {
+            matrices_ptr[i] = _matrices[i].get();
+      }
+
+      /*std::vector<Matrix> matrices(velocityMatrices.size());
+          std::vector<Vector> rhs(velocityRHSs.size());
+          for (size_t i=0; i<globalFriction_.size(); i++) {
+            matrices[i] = velocityMatrices[i];
+            rhs[i] = velocityRHSs[i];
+
+            globalFriction_[i]->addHessian(v_rel[i], matrices[i]);
+            globalFriction_[i]->addGradient(v_rel[i], rhs[i]);
+
+            matrices_ptr[i] = &matrices[i];
+          }*/
+
+      // assemble full global contact problem
+      Matrix bilinearForm;
+
+      nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+      Vector totalRhs, oldTotalRhs;
+      nBodyAssembler.assembleRightHandSide(_rhs, totalRhs);
+      oldTotalRhs = totalRhs;
+
+      Vector totalX, oldTotalX;
+      nBodyAssembler.nodalToTransformed(_x, totalX);
+      oldTotalX = totalX;
+
+      // get lower and upper obstacles
+      const auto& totalObstacles = nBodyAssembler.totalObstacles_;
+      Vector lower(totalObstacles.size());
+      Vector upper(totalObstacles.size());
+
+      for (size_t j=0; j<totalObstacles.size(); ++j) {
+          const auto& totalObstaclesj = totalObstacles[j];
+          auto& lowerj = lower[j];
+          auto& upperj = upper[j];
+        for (size_t d=0; d<dims; ++d) {
+            lowerj[d] = totalObstaclesj[d][0];
+            upperj[d] = totalObstaclesj[d][1];
+        }
+      }
+
+      // print problem
+     /* print(bilinearForm, "bilinearForm");
+      print(totalRhs, "totalRhs");
+      print(_dirichletNodes, "ignore");
+      print(totalObstacles, "totalObstacles");
+      print(lower, "lower");
+      print(upper, "upper");*/
+
+      // set up functional
+      using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
+      Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper); //TODO
+
+      // set up TNMMG solver
+      using Factory = SolverFactory<Functional, BitVector>;
+      //Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
+      Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
+      factory.build(linearSolver);
+
+     /* std::vector<BitVector> bodyDirichletNodes;
+      nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
+      for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
+        print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
+      }*/
+
+     /* print(bilinearForm, "matrix: ");
+      print(totalX, "totalX: ");
+      print(totalRhs, "totalRhs: ");*/
+
+      auto tnnmgStep = factory.step();
+      factory.setProblem(totalX);
+
+      const EnergyNorm<Matrix, Vector> norm(bilinearForm);
+
+      LoopSolver<Vector> solver(
+          *tnnmgStep.get(), _localParset.get<size_t>("maximumIterations"),
+          _localParset.get<double>("tolerance"), norm,
+          _localParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+
+      solver.preprocess();
+      solver.solve();
+
+      std::cout << "ProgramState: Energy of TNNMG solution: " << J(tnnmgStep->getSol()) << std::endl;
+
+      nBodyAssembler.postprocess(tnnmgStep->getSol(), _x);
+
+      Vector res = totalRhs;
+      bilinearForm.mmv(tnnmgStep->getSol(), res);
+      std::cout << "TNNMG Res - energy norm: " << norm.operator ()(res) << std::endl;
+
+      // TODO: remove after debugging
+  /*    using DeformedGridType = typename LevelContactNetwork<GridType, dims>::DeformedGridType;
+      using OldLinearFactory = SolverFactoryOld<DeformedGridType, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+      OldLinearFactory oldFactory(parset.sub("solver.tnnmg"), nBodyAssembler, _dirichletNodes);
+
+      auto oldStep = oldFactory.getStep();
+      oldStep->setProblem(bilinearForm, oldTotalX, oldTotalRhs);
+
+          LoopSolver<Vector> oldSolver(
+              oldStep.get(), _localParset.get<size_t>("maximumIterations"),
+              _localParset.get<double>("tolerance"), &norm,
+              _localParset.get<Solver::VerbosityMode>("verbosity"),
+              false); // absolute error
+
+          oldSolver.preprocess();
+          oldSolver.solve();
+
+      Vector oldRes = totalRhs;
+      bilinearForm.mmv(oldStep->getSol(), oldRes);
+      std::cout << "Old Res - energy norm: " << norm.operator ()(oldRes) << std::endl;*/
+
+   //   print(tnnmgStep->getSol(), "TNNMG Solution: ");
+   /*   print(oldStep->getSol(), "Old Solution: ");
+      auto diff = tnnmgStep->getSol();
+      diff -= oldStep->getSol();
+      std::cout << "Energy norm: " << norm.operator ()(diff) << std::endl;*/
+    };
+
+    timeStep = parset.get<size_t>("initialTime.timeStep");
+    relativeTime = parset.get<double>("initialTime.relativeTime");
+    relativeTau = parset.get<double>("initialTime.relativeTau");
+
+    std::vector<Vector> ell0(bodyCount_);
+    for (size_t i=0; i<bodyCount_; i++) {
+      // Initial velocity
+      u[i] = 0.0;
+      v[i] = 0.0;
+
+      ell0[i].resize(u[i].size());
+      ell0[i] = 0.0;
+
+      contactNetwork.body(i)->externalForce()(relativeTime, ell0[i]);
+    }
+
+    // Initial displacement: Start from a situation of minimal stress,
+    // which is automatically attained in the case [v = 0 = a].
+    // Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
+    BitVector dirichletNodes;
+    contactNetwork.totalNodes("dirichlet", dirichletNodes);
+    /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+        bool val = false;
+        for (size_t d=0; d<dims; d++) {
+            val = val || dirichletNodes[i][d];
+        }
+
+        dirichletNodes[i] = val;
+        for (size_t d=0; d<dims; d++) {
+            dirichletNodes[i][d] = val;
+        }
+    }*/
+
+    std::cout << "solving linear problem for u..." << std::endl;
+
+    solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, u,
+                       parset.sub("u0.solver"));
+
+    //print(u, "initial u:");
+
+    // Initial acceleration: Computed in agreement with Ma = ell0 - Au
+    // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
+    std::vector<Vector> accelerationRHS = ell0;
+    for (size_t i=0; i<bodyCount_; i++) {
+      // Initial state
+      alpha[i] = parset.get<double>("boundary.friction.initialAlpha");
+
+      // Initial normal stress
+
+      const auto& body = contactNetwork.body(i);
+      /*std::vector<std::shared_ptr<typename ContactNetwork::LeafBody::BoundaryCondition>> frictionBoundaryConditions;
+      body->boundaryConditions("friction", frictionBoundaryConditions);
+      for (size_t j=0; j<frictionBoundaryConditions.size(); j++) {
+          ScalarVector frictionBoundaryStress(weightedNormalStress[i].size());
+
+          body->assembler()->assembleWeightedNormalStress(
+            *frictionBoundaryConditions[j]->boundaryPatch(), frictionBoundaryStress, body->data()->getYoungModulus(),
+            body->data()->getPoissonRatio(), u[i]);
+
+          weightedNormalStress[i] += frictionBoundaryStress;
+      }*/
+
+      Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
+    }
+
+    for (size_t i=0; i<contactNetwork.nCouplings(); i++) {
+      const auto& coupling = contactNetwork.coupling(i);
+      const auto& contactCoupling = contactNetwork.nBodyAssembler().getContactCouplings()[i];
+
+      const auto nonmortarIdx = coupling->gridIdx_[0];
+      const auto& body = contactNetwork.body(nonmortarIdx);
+
+      ScalarVector frictionBoundaryStress(weightedNormalStress[nonmortarIdx].size());
+
+      body->assembler()->assembleWeightedNormalStress(
+        contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
+        body->data()->getPoissonRatio(), u[nonmortarIdx]);
+
+      weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
+    }
+
+    std::cout << "solving linear problem for a..." << std::endl;
+
+    BitVector noNodes(dirichletNodes.size(), false);
+    solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, a,
+                       parset.sub("a0.solver"));
+
+    //print(a, "initial a:");
+  }
+
+private:
+  const size_t bodyCount_;
+
+public:
+  std::vector<BodyState* > bodyStates;
+  std::vector<Vector> u;
+  std::vector<Vector> v;
+  std::vector<Vector> a;
+  std::vector<ScalarVector> alpha;
+  std::vector<ScalarVector> weightedNormalStress;
+  double relativeTime;
+  double relativeTau;
+  size_t timeStep;
+
+
+};
+#endif

dune/tectonic/explicitgrid.hh

+#ifndef SRC_EXPLICITGRID_HH
+#define SRC_EXPLICITGRID_HH
+
+#include "gridselector.hh"
+#include "explicitvectors.hh"
+
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+
+using LeafGridView = Grid::LeafGridView;
+using LevelGridView = Grid::LevelGridView;
+
+using DeformedGridComplex = typename Dune::Contact::DeformedContinuaComplex<Grid, Vector>;
+using DeformedGrid = DeformedGridComplex::DeformedGridType;
+
+using DefLeafGridView = DeformedGrid::LeafGridView;
+using DefLevelGridView = DeformedGrid::LevelGridView;
+
+#endif

src/explicitvectors.hh → dune/tectonic/explicitvectors.hh

@@ -3,13 +3,18 @@
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
-using LocalVector = Dune::FieldVector<double, MY_DIM>;
-using LocalMatrix = Dune::FieldMatrix<double, MY_DIM, MY_DIM>;
+using ctype = double;
+
+using LocalVector = Dune::FieldVector<ctype, MY_DIM>;
+using LocalMatrix = Dune::FieldMatrix<ctype, MY_DIM, MY_DIM>;
 using Vector = Dune::BlockVector<LocalVector>;
 using Matrix = Dune::BCRSMatrix<LocalMatrix>;
-using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
-using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+using ScalarVector = Dune::BlockVector<Dune::FieldVector<ctype, 1>>;
+using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<ctype, 1, 1>>;
+using BitVector = Dune::BitSetVector<MY_DIM>;
+
 #endif

dune/tectonic/factories/CMakeLists.txt

+add_custom_target(tectonic_dune_factories SOURCES
+  cantorfactory.hh
+  cantorfactory.cc 
+  contactnetworkfactory.hh
+  levelcontactnetworkfactory.hh
+  stackedblocksfactory.hh
+  stackedblocksfactory.cc 
+  threeblocksfactory.hh
+  threeblocksfactory.cc 
+  twoblocksfactory.hh
+  twoblocksfactory.cc 
+)
+
+#install headers
+install(FILES
+  cantorfactory.hh
+  contactnetworkfactory.hh
+  levelcontactnetworkfactory.hh
+  stackedblocksfactory.hh
+  threeblocksfactory.hh
+  twoblocksfactory.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)