src/data-structures/body.hh

+#ifndef SRC_DATA_STRUCTURES_BODY_HH
+#define SRC_DATA_STRUCTURES_BODY_HH
+
+#include <functional>
+#include <type_traits>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+//#include <dune/fufem/assemblers/assembler.hh>
+//#pragma clang diagnostic push
+//#pragma clang diagnostic ignored "-Wsign-compare"
+//#include <dune/fufem/functionspacebases/p0basis.hh>
+//#pragma clang diagnostic pop
+//#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+
+#include <dune/grid/geometrygrid/grid.hh>
+
+#include <dune/fufem/functions/deformationfunction.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+
+#include <dune/tectonic/bodydata.hh>
+
+#include "../assemblers.hh"
+#include "../boundarycondition.hh"
+#include "enums.hh"
+#include "matrices.hh"
+
+template <class HostGridType, class VectorType>
+class LeafBody {
+public:
+    static const int dim = HostGridType::dimensionworld;
+
+    using DeformationFunction = DeformationFunction<typename HostGridType::LeafGridView, VectorType>;
+    using GridType = Dune::GeometryGrid<HostGridType, DeformationFunction>;
+    using GridView = typename GridType::LeafGridView;
+
+    using Function = Dune::VirtualFunction<double, double>;
+    using ExternalForce = std::function<void(double, VectorType &)>;
+
+    using Assembler = MyAssembler<GridView, dim>;
+    using Matrix = typename Assembler::Matrix;
+    using LocalVector = typename VectorType::block_type;
+
+    using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+
+    using BoundaryCondition = BoundaryCondition<GridView, dim>;
+
+    using BoundaryFunctions = std::vector<const Function* >;
+    using BoundaryNodes = std::vector<const Dune::BitSetVector<dim>* >;
+    using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
+    using BoundaryPatches = std::vector<const typename BoundaryCondition::BoundaryPatch* >;
+
+    using StateEnergyNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+
+private:
+    std::shared_ptr<BodyData<dim>> bodyData_;
+
+    std::shared_ptr<HostGridType> hostGrid_;
+    std::shared_ptr<DeformationFunction> deformation_;
+    std::shared_ptr<GridType> grid_;
+    GridView gridView_;
+
+    std::shared_ptr<Assembler> assembler_;
+    Matrices<Matrix,1> matrices_;
+
+    ExternalForce externalForce_;
+
+    VectorType gravityFunctional_;
+
+    // boundary conditions
+    std::vector<std::shared_ptr<BoundaryCondition>> boundaryConditions_;
+
+public:
+    LeafBody(
+            const std::shared_ptr<BodyData<dim>>& bodyData,
+            const std::shared_ptr<HostGridType>& hostGrid);
+
+    void assemble();
+
+    // setter and getter
+    auto data() const -> const std::shared_ptr<BodyData<dim>>&;
+
+    void setDeformation(const VectorType& def);
+    auto deformation() const -> DeformationFunction&;
+    auto grid() const -> std::shared_ptr<GridType>;
+
+    auto gridView() const -> const GridView&;
+    auto nVertices() const -> size_t;
+
+    auto assembler() const -> const std::shared_ptr<Assembler>&;
+    auto matrices() const -> const Matrices<Matrix,1>&;
+
+    auto externalForce() const -> const ExternalForce&;
+
+    void addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition);
+    void boundaryConditions(
+            const std::string& tag,
+            std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const;
+    auto boundaryConditions() const -> const std::vector<std::shared_ptr<BoundaryCondition>>&;
+
+    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 boundaryPatches(
+            const std::string& tag,
+            BoundaryPatches& patches) const;
+};
+
+template <class GridView>
+class Body {
+public:
+    enum {dim = GridView::dimensionworld};
+
+    using GridType = typename GridView::Grid;
+
+    using Function = Dune::VirtualFunction<double, double>;
+
+    using BoundaryCondition = BoundaryCondition<GridView, dim>;
+
+    using BoundaryFunctions = std::vector<const Function* >;
+    using BoundaryNodes = std::vector<const Dune::BitSetVector<dim>* >;
+    using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
+    using BoundaryPatches = std::vector<const typename BoundaryCondition::BoundaryPatch* >;
+
+private:
+    std::shared_ptr<BodyData<dim>> bodyData_;
+
+    std::shared_ptr<GridType> grid_;
+    const size_t level_;
+    GridView gridView_;
+
+    // boundary conditions
+    std::vector<std::shared_ptr<BoundaryCondition>> boundaryConditions_;
+
+public:
+    template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LeafGridView>::value, int> = 0>
+    Body(
+            const std::shared_ptr<BodyData<dim>>& bodyData,
+            const std::shared_ptr<GridType>& grid) :
+        bodyData_(bodyData),
+        grid_(grid),
+        level_(grid_->maxLevel()),
+        gridView_(grid_->leafGridView()) {
+
+        boundaryConditions_.resize(0);
+    }
+
+    template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LevelGridView>::value, int> = 0>
+    Body(
+            const std::shared_ptr<BodyData<dim>>& bodyData,
+            const std::shared_ptr<GridType>& grid,
+            const size_t level) :
+        bodyData_(bodyData),
+        grid_(grid),
+        level_(level),
+        gridView_(grid_->levelGridView(level_)) {
+
+        boundaryConditions_.resize(0);
+    }
+
+
+    // setter and getter
+    auto data() const -> const std::shared_ptr<BodyData<dim>>&;
+    auto grid() const -> std::shared_ptr<GridType>;
+
+    //template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LevelGridView>::value, int> = 0>
+    auto level() const
+    -> size_t {
+
+        return level_;
+    }
+
+    auto gridView() const -> const GridView&;
+    auto nVertices() const -> size_t;
+
+    void addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition);
+    void boundaryConditions(
+            const std::string& tag,
+            std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const;
+    auto boundaryConditions() const -> const std::vector<std::shared_ptr<BoundaryCondition>>&;
+
+    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 boundaryPatches(
+            const std::string& tag,
+            BoundaryPatches& patches) const;
+};
+#endif

src/data-structures/body_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+template class LeafBody<Grid, Vector>;
+
+template class Body<DefLeafGridView>;
+template class Body<DefLevelGridView>;

src/data-structures/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 <dune/tectonic/globalratestatefriction.hh>
+#include <dune/tectonic/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)
+{}
+
+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);
+    }
+
+    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()
+-> Dune::Contact::NBodyAssembler<GridType, VectorType>& {
+
+    return nBodyAssembler_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler() const
+-> const Dune::Contact::NBodyAssembler<GridType, VectorType>& {
+
+    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);
+
+        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);
+        }
+    }
+}
+
+// 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"

src/data-structures/contactnetwork.hh

+#ifndef SRC_CONTACTNETWORK_HH
+#define SRC_CONTACTNETWORK_HH
+
+#include <memory>
+
+#include <dune/common/promotiontraits.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
+#include <dune/tectonic/bodydata.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "../assemblers.hh"
+#include "../frictioncouplingpair.hh"
+#include "body.hh"
+#include "enums.hh"
+#include "matrices.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 = Dune::Contact::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

src/data-structures/contactnetwork_tmpl.cc

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

src/enumparser.hh → src/data-structures/enumparser.hh

-#ifndef SRC_ENUMPARSER_HH
-#define SRC_ENUMPARSER_HH
+#ifndef SRC_DATA_STRUCTURES_ENUMPARSER_HH
+#define SRC_DATA_STRUCTURES_ENUMPARSER_HH
 
 // Copyright Carsten Graeser 2012
 

src/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"
+*/

src/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

src/data-structures/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 <dune/tectonic/globalratestatefriction.hh>
+#include <dune/tectonic/frictionpotential.hh>
+
+#include "levelcontactnetwork.hh"
+
+#include "../assemblers.hh"
+#include "../utils/tobool.hh"
+#include "../utils/debugutils.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"

src/data-structures/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 <dune/tectonic/bodydata.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "../assemblers.hh"
+#include "../frictioncouplingpair.hh"
+#include "body.hh"
+#include "enums.hh"
+#include "matrices.hh"
+//#include "multi-body-problem-data/myglobalfrictiondata.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

src/data-structures/levelcontactnetwork_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include "../frictioncouplingpair.hh"
+#include "contactnetwork_tmpl.cc"
+#include "levelcontactnetwork.hh"
+
+template class LevelContactNetwork<typename MyContactNetwork::GridType, typename MyContactNetwork::FrictionCouplingPair, typename MyContactNetwork::field_type>;

src/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

src/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/tnnmg/problem-classes/blocknonlineartnnmgproblem.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 <dune/tectonic/bodydata.hh>
+
+#include "../assemblers.hh"
+#include "contactnetwork.hh"
+#include "matrices.hh"
+#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
+#include "../spatial-solving/solverfactory.hh"
+#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/preconditioners/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();
+
+    const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
+
+    std::cout << "-- setupInitialConditions --" << std::endl;
+    std::cout << "----------------------------" << std::endl;
+
+    /*
+    std::cout << "Building preconditioner..." << std::endl;
+    using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+    Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+    /*Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), false);
+    activeLevels.back() = true;
+    activeLevels[activeLevels.size()-2] = true;*/
+    /*print(activeLevels, "activeLevels:");
+    Preconditioner preconditioner(preconditionerParset, contactNetwork, activeLevels);
+    preconditioner.build();
+
+    using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
+
+    LinearSolverStep linearSolverStep;
+    linearSolverStep.setPreconditioner(preconditioner);
+
+    EnergyNorm<Matrix, Vector> energyNorm(linearSolverStep);
+    LinearSolver linearSolver(linearSolverStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), energyNorm, Solver::FULL);
+    */
+    using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    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 smoother = TruncatedBlockGSStep<Matrix, Vector>{};
+    auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+    linearMultigridStep->setMGType(1, 3, 3);
+    linearMultigridStep->setSmoother(&smoother);
+    linearMultigridStep->setTransferOperators(transfer);
+
+    EnergyNorm<Matrix, Vector> mgNorm(*linearMultigridStep);
+    LinearSolver mgSolver(linearMultigridStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), mgNorm, Solver::FULL);
+
+
+    // 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, mgSolver, _dirichletNodes);
+
+     /* 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"),
+          false); // absolute error
+
+      solver.preprocess();
+      solver.solve();
+
+      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]);
+    }
+
+    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

src/explicitgrid.hh

@@ -2,6 +2,17 @@
 #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;
 
-using GridView = Grid::LeafGridView;
 #endif

src/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

src/factories/cantorfactory.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../utils/almostequal.hh"
+
+#include "cantorfactory.hh"
+
+#include "../multi-body-problem-data/bc.hh"
+#include "../multi-body-problem-data/cuboidgeometry.hh"
+#include "../multi-body-problem-data/myglobalfrictiondata.hh"
+#include "../multi-body-problem-data/weakpatch.hh"
+#include "../multi-body-problem-data/cubegridconstructor.hh"
+
+#include "../utils/diameter.hh"
+
+template <class GridType, int dim>
+void CantorFactory<GridType, dim>::setBodies() {
+
+
+
+    for (size_t level=0; level<=maxLevel_; level++) {
+        const LevelCubes& levelCubes = cubes_[level];
+
+        for (size_t i=0; i<levelCubes.size(); i++) {
+            const std::shared_ptr<Cube>& cube = levelCubes[i];
+
+            CubeGridConstructor
+        }
+    }
+
+    std::vector<LocalVector> origins(this->bodyCount_);
+    std::vector<LocalVector> lowerWeakOrigins(this->bodyCount_);
+    std::vector<LocalVector> upperWeakOrigins(this->bodyCount_);
+
+#if MY_DIM == 3
+        double const depth = 0.60;
+
+        origins[0] = {0,0,0};
+        lowerWeakOrigins[0] = {0.2,0,0};
+        upperWeakOrigins[0] = {0.2,width,0};
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width, depth);
+        for (size_t i=1; i<this->bodyCount_-1; i++) {
+            origins[i] = cuboidGeometries_[i-1]->D;
+            lowerWeakOrigins[i] = {0.6,i*width,0};
+            upperWeakOrigins[i] = {0.6,(i+1)*width,0};
+
+            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width, depth);
+        }
+        const size_t idx = this->bodyCount_-1;
+        origins[idx] = cuboidGeometries_[idx-1]->D;
+        lowerWeakOrigins[idx] = {0.6,idx*width,0};
+        upperWeakOrigins[idx] = {0.6,(idx+1)*width,0};
+
+        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, length, width, depth);
+#elif MY_DIM == 2
+        origins[0] = {0,0};
+        lowerWeakOrigins[0] = {0.2,0};
+        upperWeakOrigins[0] = {0.2,width};
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width);
+        for (size_t i=1; i<this->bodyCount_-1; i++) {
+            origins[i] = cuboidGeometries_[i-1]->D;
+            lowerWeakOrigins[i] = {0.6,i*width};
+            upperWeakOrigins[i] = {0.6,(i+1)*width};
+
+            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width);
+        }
+        const size_t idx = this->bodyCount_-1;
+        origins[idx] = cuboidGeometries_[idx-1]->D;
+        lowerWeakOrigins[idx] = {0.6,idx*width};
+        upperWeakOrigins[idx] = {0.6,(idx+1)*width};
+
+        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, length, width);
+#else
+#error CuboidGeometry only supports 2D and 3D!"
+#endif
+
+    // set up reference grids
+    gridConstructor_ = new GridsConstructor<GridType>(cuboidGeometries_);
+    auto& grids = gridConstructor_->getGrids();
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& cuboidGeometry = *cuboidGeometries_[i];
+
+        // define weak patch and refine grid
+        auto lowerWeakPatch = lowerWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
+        auto upperWeakPatch = upperWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
+        getWeakPatch<LocalVector>(this->parset_.sub("boundary.friction.weakening"), cuboidGeometry, *lowerWeakPatch, *upperWeakPatch);
+        refine(*grids[i], *lowerWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
+        refine(*grids[i], *upperWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
+    }
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        this->bodies_[i] = std::make_shared<typename Base::Body>(bodyData_, grids[i]);
+    }
+}
+
+template <class GridType, int dim>
+void CantorFactory<GridType, dim>::setCouplings() {
+    const auto deformedGrids = this->contactNetwork_.deformedGrids();
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& cuboidGeometry = *cuboidGeometries_[i];
+        leafFaces_[i] = std::make_shared<LeafFaces>(this->contactNetwork_.leafView(i), cuboidGeometry);
+        levelFaces_[i] = std::make_shared<LevelFaces>(this->contactNetwork_.levelView(i), cuboidGeometry);
+    }
+
+    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+    std::shared_ptr<typename Base::CouplingPair::BackEndType> backend = nullptr;
+    for (size_t i=0; i<this->couplings_.size(); i++) {
+      auto& coupling = this->couplings_[i];
+      coupling = std::make_shared<typename Base::CouplingPair>();
+
+      nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
+      mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
+
+      coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::CouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
+    }
+}
+
+template <class GridType, int dim>
+void CantorFactory<GridType, dim>::setBoundaryConditions() {
+    using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dim>;
+    using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dim>;
+
+    using Function = Dune::VirtualFunction<double, double>;
+    std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+    std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>();
+
+    const double lengthScale = CuboidGeometry::lengthScale;
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+      const auto& body = this->contactNetwork_.body(i);
+      const auto& leafVertexCount = body->leafVertexCount();
+
+      std::cout << "Grid" << i << " Number of DOFs: " << leafVertexCount << std::endl;
+
+      // friction boundary
+      if (i<this->bodyCount_-1 && upperWeakPatches_[i]->vertices.size()>0) {
+        std::shared_ptr<FrictionBoundaryCondition> frictionBoundary = std::make_shared<FrictionBoundaryCondition>();
+        frictionBoundary->setBoundaryPatch(leafFaces_[i]->upper);
+        frictionBoundary->setWeakeningPatch(upperWeakPatches_[i]);
+        frictionBoundary->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *upperWeakPatches_[i], lengthScale));
+        body->addBoundaryCondition(frictionBoundary);
+      }
+      if (i>0 && lowerWeakPatches_[i]->vertices.size()>0) {
+        std::shared_ptr<FrictionBoundaryCondition> frictionBoundary = std::make_shared<FrictionBoundaryCondition>();
+        frictionBoundary->setBoundaryPatch(leafFaces_[i]->lower);
+        frictionBoundary->setWeakeningPatch(lowerWeakPatches_[i]);
+        frictionBoundary->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *lowerWeakPatches_[i], lengthScale));
+        body->addBoundaryCondition(frictionBoundary);
+      }
+
+      // Neumann boundary
+      std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->leafView()), neumannFunction, "neumann");
+      body->addBoundaryCondition(neumannBoundary);
+
+      // Dirichlet Boundary
+      // identify Dirichlet nodes on leaf level
+      std::shared_ptr<Dune::BitSetVector<dim>> dirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
+      for (int j=0; j<leafVertexCount; j++) {
+        if (leafFaces_[i]->right.containsVertex(j))
+          (*dirichletNodes)[j][0] = true;
+
+        if (leafFaces_[i]->lower.containsVertex(j))
+          (*dirichletNodes)[j][0] = true;
+
+        #if MY_DIM == 3
+        if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
+            dirichletNodes->at(j)[2] = true;
+        #endif
+      }
+
+      std::shared_ptr<LeafBoundaryCondition> dirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+      dirichletBoundary->setBoundaryPatch(body->leafView(), dirichletNodes);
+      dirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+      body->addBoundaryCondition(dirichletBoundary);
+    }
+}
+
+
+class CantorIndexHierarchy {
+    public:
+        typedef std::map<std::pair<int,int>, bool> LevelCantorIndices;
+
+    private:
+        const size_t maxLevel_;
+        std::vector<size_t> levelN_;
+        std::vector<LevelCantorIndices> cantorIndexHierarchy_;
+
+        void prolongIndices(size_t currentLevel, size_t newLevel) {
+            const LevelCantorIndices& indices = cantorIndexHierarchy_[currentLevel];
+            LevelCantorIndices& newIndices = cantorIndexHierarchy_[newLevel];
+
+            size_t offset = levelN_[currentLevel];
+
+            std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin();
+            std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end();
+            for (; it!=endIt; ++it) {
+                int xID = it->first.first;
+                int yID = it->first.second;
+                newIndices[std::make_pair(xID, yID)] = true;
+                newIndices[std::make_pair(xID+offset, yID)] = true;
+                newIndices[std::make_pair(xID, yID+offset)] = true;
+            }
+
+            size_t doubleOffset = 2*offset;
+            for (size_t i=offset+1; i<=doubleOffset; i=i+2) {
+                newIndices[std::make_pair(doubleOffset, i)] = true;
+                newIndices[std::make_pair(i, doubleOffset)] = true;
+            }
+        }
+
+        void print(const LevelCantorIndices& indices) const {
+            std::cout << "LevelCantorIndices: " << std::endl;
+
+            std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin();
+            std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end();
+            for (; it!=endIt; ++it) {
+                std::cout << "(" << it->first.first << ", " << it->first.second << ")"<< std::endl;
+            }
+        }
+
+    public:
+        CantorIndexHierarchy(size_t maxLevel) :
+            maxLevel_(maxLevel)
+        {
+            levelN_.resize(maxLevel_+1);
+            cantorIndexHierarchy_.resize(maxLevel_+1);
+
+            // init levelCantorIndices on level 0
+            LevelCantorIndices& initCantorIndices = cantorIndexHierarchy_[0];
+            initCantorIndices[std::make_pair(0, 1)] = true;
+            initCantorIndices[std::make_pair(1, 0)] = true;
+            initCantorIndices[std::make_pair(1, 2)] = true;
+            initCantorIndices[std::make_pair(2, 1)] = true;
+
+            for (size_t i=0; i<maxLevel_; i++) {
+                levelN_[i] = std::pow(2, i+1);
+                prolongIndices(i, i+1);
+            }
+
+            levelN_[maxLevel_] = std::pow(2, maxLevel_+1);
+        }
+
+        const LevelCantorIndices& levelCantorIndices(size_t i) const {
+                return cantorIndexHierarchy_[i];
+        }
+
+        size_t levelN(const size_t i) const {
+            return levelN_[i];
+        }
+};
+
+
+template <class GridType>
+class CantorFaultFactory
+{
+    //! Parameter for mapper class
+    template<int dim>
+    struct FaceMapperLayout
+    {
+        bool contains (Dune::GeometryType gt)
+        {
+            return gt.dim() == dim-1;
+        }
+    };
+
+    protected:
+        typedef OscUnitCube<GridType, 2> GridOb;
+        static const int dimworld = GridType::dimensionworld;
+        static const int dim = GridType::dimension;
+        typedef typename GridType::ctype ctype;
+        typedef typename GridType::LevelGridView GV;
+
+        typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;
+
+    private:
+        const std::map<size_t, size_t>& levelToCantorLevel_;
+        const int coarseResolution_;
+        const size_t maxLevel_;
+        const size_t maxCantorLevel_;
+        const CantorIndexHierarchy cantorIndexHierarchy_;
+        const int coarseGridN_;
+
+        std::shared_ptr<GridType> grid_;
+        std::vector<double> levelResolutions_;
+
+        std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_;
+
+private:
+
+
+    bool computeID(Dune::FieldVector<ctype, dimworld> vertex, const size_t gridN, std::pair<size_t, size_t>& IDs) const {
+        ctype xID = vertex[0]*gridN;
+        ctype yID = vertex[1]*gridN;
+
+        ctype x = 0;
+        ctype y = 0;
+
+        bool xIsInt = almost_equal(std::modf(xID, &x), 0.0, 2);
+        bool yIsInt = almost_equal(std::modf(yID, &y), 0.0, 2);
+
+        if (xIsInt && yIsInt) {
+            IDs = std::make_pair((size_t) x, (size_t) y);
+            return true;
+        }
+
+        if (xIsInt) {
+            y = std::ceil(yID);
+            if ((size_t) y % 2==0)
+                y = y-1;
+            IDs = std::make_pair((size_t) x, (size_t) y);
+            return true;
+        }
+
+        if (yIsInt) {
+            x = std::ceil(xID);
+            if ((size_t) x % 2==0)
+                x = x-1;
+            IDs = std::make_pair((size_t) x, (size_t) y);
+            return true;
+        }
+
+        return false;
+    }
+
+public:
+        //setup
+    CantorFaultFactory(const std::map<size_t, size_t>& levelToCantorLevel, const int coarseResolution, const size_t maxLevel, const size_t maxCantorLevel) :
+        levelToCantorLevel_(levelToCantorLevel),
+        coarseResolution_(coarseResolution),
+        maxLevel_(maxLevel),
+        maxCantorLevel_(maxCantorLevel),
+        cantorIndexHierarchy_(maxCantorLevel_),
+        coarseGridN_(std::pow(2, coarseResolution_)),
+        interfaceNetwork_(nullptr)
+        {
+            Dune::UGGrid<dim>::setDefaultHeapSize(4000);
+            GridOb unitCube(coarseGridN_);
+            grid_ = unitCube.grid();
+
+            grid_->globalRefine(maxLevel_);
+
+            levelResolutions_.resize(maxLevel_+1);
+
+            // init interface network
+            interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);
+
+            for (size_t i=0; i<=maxLevel_; i++) {
+                if (i>0) {
+                    interfaceNetwork_->prolongLevel(i-1, i);
+                }
+
+                if (levelToCantorLevel_.count(i)) {
+
+                    levelResolutions_[i] = std::pow(2, coarseResolution_+i);
+                    const size_t levelResolution = levelResolutions_[i];
+                    const size_t cantorResolution = cantorIndexHierarchy_.levelN(levelToCantorLevel_.at(i));
+
+                    if (2*levelResolution<cantorResolution)
+                        DUNE_THROW(Dune::Exception, "Level " + std::to_string(i) + " does not support cantor set with resolution " + std::to_string(cantorResolution));
+
+                    const typename CantorIndexHierarchy::LevelCantorIndices& levelCantorIndices = cantorIndexHierarchy_.levelCantorIndices(levelToCantorLevel_.at(i));
+                    std::shared_ptr<LevelInterfaceNetwork<GV>> levelInterfaceNetwork = interfaceNetwork_->levelInterfaceNetworkPtr(i);
+
+                    const GV& gridView = levelInterfaceNetwork->levelGridView();
+                    FaceMapper intersectionMapper(gridView);
+                    std::vector<bool> intersectionHandled(intersectionMapper.size(),false);
+
+                    for (const auto& elem:elements(gridView)) {
+                        for (const auto& isect:intersections(gridView, elem)) {
+
+                            if (intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)])
+                                continue;
+
+                            intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)]=true;
+
+                            if (isect.boundary())
+                                continue;
+
+                            std::pair<size_t, size_t> intersectionID;
+                            bool isAdmissibleID = computeID(isect.geometry().center(), cantorResolution, intersectionID);
+
+                            if (isAdmissibleID && levelCantorIndices.count(intersectionID)) {
+                                levelInterfaceNetwork->addIntersection(isect, i);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+
+
+#include "cantorfactory_tmpl.cc"

src/factories/cantorfactory.hh

+#ifndef SRC_CANTORFACTORY_HH
+#define SRC_CANTORFACTORY_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "contactnetworkfactory.hh"
+
+#include "../multi-body-problem-data/mybody.hh"
+#include "../multi-body-problem-data/grid/mygrids.hh"
+
+
+
+
+template <class GridType, int dim> class CantorFactory : public ContactNetworkFactory<GridType, dim>{
+private:
+    using Base = ContactNetworkFactory<GridType, dim>;
+
+    using LevelCubes = std::vector<std::shared_ptr<Cube>>;
+
+    using LocalVector = typename Base::ContactNetwork::LocalVector;
+
+    using DeformedLeafGridView = typename Base::ContactNetwork::DeformedLeafGridView;
+    using DeformedLevelGridView = typename Base::ContactNetwork::DeformedLevelGridView;
+
+    using LevelBoundaryPatch = BoundaryPatch<DeformedLevelGridView>;
+    using LeafBoundaryPatch = BoundaryPatch<DeformedLeafGridView>;
+
+    using LeafFaces = MyFaces<DeformedLeafGridView>;
+    using LevelFaces = MyFaces<DeformedLevelGridView>;
+
+private:
+    const std::shared_ptr<MyBodyData<dim>> bodyData_;     // material properties of bodies
+
+    GridsConstructor<GridType>* gridConstructor_;
+
+    std::vector<LevelCubes> cubes_;
+
+    std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
+    std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
+
+    std::vector<std::shared_ptr<ConvexPolyhedron<LocalVector>>> lowerWeakPatches_;
+    std::vector<std::shared_ptr<ConvexPolyhedron<LocalVector>>> upperWeakPatches_;
+
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
+
+public:
+    CantorFactory(const Dune::ParameterTree& parset) :
+        Base(parset, parset.get<size_t>("problem.bodyCount"), parset.get<size_t>("problem.bodyCount")-1),
+        bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_, zenith_())),
+        cuboidGeometries_(this->bodyCount_),
+        leafFaces_(this->bodyCount_),
+        levelFaces_(this->bodyCount_),
+        lowerWeakPatches_(this->bodyCount_),
+        upperWeakPatches_(this->bodyCount_),
+        nonmortarPatch_(this->couplingCount_),
+        mortarPatch_(this->couplingCount_)
+    {}
+
+    ~CantorFactory() {
+        delete gridConstructor_;
+    }
+
+    void setBodies();
+    void setCouplings();
+    void setBoundaryConditions();
+
+private:
+    static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
+        #if MY_DIM == 2
+        return {0, 1};
+        #elif MY_DIM == 3
+        return {0, 1, 0};
+        #endif
+    }
+};
+#endif
+
+

src/factories/cantorfactory_tmpl.cc

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