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;
+        case Config::Tresca:
+          globalFriction_ = std::make_shared<GlobalRateStateFriction<
+            Matrix, VectorType, Tresca, GridType>>(
+            nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
+        case Config::None:
+          globalFriction_ = std::make_shared<GlobalRateStateFriction<
+            Matrix, VectorType, ZeroFunction, GridType>>(
+            nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
+        default:
+          assert(false);
+          break;
+      }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::setBodies(const std::vector<std::shared_ptr<LeafBody>> bodies) {
+    assert(nBodies()==bodies.size());
+    bodies_ = bodies;
+
+    matrices_.elasticity.resize(nBodies());
+    matrices_.damping.resize(nBodies());
+    matrices_.mass.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        matrices_.elasticity[i] = bodies_[i]->matrices().elasticity;
+        matrices_.damping[i] = bodies_[i]->matrices().damping;
+        matrices_.mass[i] = bodies_[i]->matrices().mass;
+    }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings) {
+    assert(this->nCouplings()==couplings.size());
+    couplings_ = couplings;
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::constructBody(
+        const std::shared_ptr<BodyData<dim>>& bodyData,
+        const std::shared_ptr<HostGridType>& grid,
+        std::shared_ptr<LeafBody>& body) const {
+
+    body = std::make_shared<LeafBody>(bodyData, grid);
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::setDeformation(const std::vector<VectorType>& totalDeformation) {
+    assert(totalDeformation.size() == nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        body(i)->setDeformation(totalDeformation[i]);
+    }
+
+    nBodyAssembler_.assembleTransferOperator();
+    nBodyAssembler_.assembleObstacle();
+
+    for (size_t i=0; i<levelContactNetworks_.size(); i++) {
+        levelContactNetworks_[i]->build();
+    }
+}
+
+template <class HostGridType, class VectorType>
+template <class LevelBoundaryPatch>
+void ContactNetwork<HostGridType, VectorType>::constructCoupling(
+        int nonMortarBodyIdx, int mortarBodyIdx,
+        const std::shared_ptr<LevelBoundaryPatch>& nonmortarPatch,
+        const std::shared_ptr<LevelBoundaryPatch>& mortarPatch,
+        const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
+        const std::shared_ptr<GlobalFrictionData<dim>>& frictionData,
+        std::shared_ptr<FrictionCouplingPair>& coupling) const {
+
+    coupling = std::make_shared<FrictionCouplingPair>();
+
+    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<BoundaryPatch>>();
+    std::shared_ptr<typename FrictionCouplingPair::BackEndType> backend = nullptr;
+    coupling->set(nonMortarBodyIdx, mortarBodyIdx, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
+    coupling->setWeakeningPatch(weakeningPatch);
+    coupling->setFrictionData(frictionData);
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::level(size_t level) const
+-> const std::shared_ptr<LevelContactNetwork>& {
+
+    return levelContactNetworks_[level];
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nLevels() const
+-> size_t {
+
+    return levelContactNetworks_.size();
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodies() const
+-> size_t {
+
+    return bodies_.size();
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nCouplings() const
+-> size_t {
+
+    return couplings_.size();
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::body(int i) const
+-> const std::shared_ptr<LeafBody>& {
+
+    return bodies_[i];
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::coupling(int i) const
+-> const std::shared_ptr<FrictionCouplingPair>& {
+
+    return couplings_[i];
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::couplings() const
+-> const std::vector<std::shared_ptr<FrictionCouplingPair>>& {
+
+    return couplings_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler()
+-> NBodyAssembler& {
+
+    return nBodyAssembler_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler() const
+-> const NBodyAssembler& {
+
+    return nBodyAssembler_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::matrices() const
+-> const Matrices<Matrix,2>& {
+
+    return matrices_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::stateEnergyNorms() const
+-> const StateEnergyNorms& {
+
+    return stateEnergyNorms_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::globalFriction() const
+-> GlobalFriction& {
+
+    return *globalFriction_;
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::externalForces(ExternalForces& externalForces) const {
+    externalForces.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        externalForces[i] = std::make_unique<typename LeafBody::ExternalForce>(bodies_[i]->externalForce());
+    }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::totalNodes(
+        const std::string& tag,
+        Dune::BitSetVector<dim>& totalNodes) const {
+
+    int totalSize = 0;
+    for (size_t i=0; i<nBodies(); i++) {
+        totalSize += this->body(i)->nVertices();
+    }
+
+    totalNodes.resize(totalSize);
+
+    int idx=0;
+    for (size_t i=0; i<nBodies(); i++) {
+        const auto& body = this->body(i);
+        std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
+        body->boundaryConditions(tag, boundaryConditions);
+
+        if (boundaryConditions.size()>0) {
+            const int idxBackup = idx;
+            for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
+                const auto& nodes = boundaryConditions[bc]->boundaryNodes();
+                for (size_t j=0; j<nodes->size(); j++, idx++)
+                    for (int k=0; k<dim; k++)
+                        totalNodes[idx][k] = totalNodes[idx][k] or (*nodes)[j][k];
+                idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+            }
+        } else {
+            idx += body->nVertices();
+        }
+    }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::boundaryPatchNodes(
+        const std::string& tag,
+        BoundaryPatchNodes& nodes) const {
+
+    nodes.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryPatchNodes(tag, nodes[i]);
+    }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::boundaryNodes(
+        const std::string& tag,
+        BoundaryNodes& nodes) const {
+
+    nodes.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryNodes(tag, nodes[i]);
+    }
+}
+
+// collects all leaf boundary functions from the different bodies identified by "tag"
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::boundaryFunctions(
+        const std::string& tag,
+        BoundaryFunctions& functions) const {
+
+    functions.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryFunctions(tag, functions[i]);
+    }
+}
+
+// collects all leaf boundary patches from the different bodies identified by "tag"
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::boundaryPatches(
+        const std::string& tag,
+        BoundaryPatches& patches) const {
+
+    patches.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        this->body(i)->boundaryPatches(tag, patches[i]);
+    }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const {
+    nodes.resize(nBodies());
+    for (size_t i=0; i<nBodies(); i++) {
+        nodes[i] = frictionBoundaries_[i]->getVertices();
+    }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::frictionPatches(std::vector<const BoundaryPatch*>& patches) const {
+    patches.resize(nBodies());
+    for (size_t i=0; i<nBodies(); i++) {
+        patches[i] = frictionBoundaries_[i].get();
+    }
+}
+
+#include "contactnetwork_tmpl.cc"

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

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

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

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

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

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

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/makelinearsolver.hh"
+#include "../spatial-solving/nonlinearsolver.hh"
+#include "../spatial-solving/solverfactory.hh"
+#include "../spatial-solving/functionalfactory.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/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_),
+      weights(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),
+      weights[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) {
+
+    std::cout << "-- setupInitialConditions --" << std::endl;
+    std::cout << "----------------------------" << std::endl;
+
+    using Matrix = typename ContactNetwork::Matrix;
+    const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+
+    auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(parset, contactNetwork);
+    auto nonlinearity = ZeroNonlinearity();
+
+    // 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) {
+
+        Vector totalX;
+        nBodyAssembler.nodalToTransformed(_x, totalX);
+
+        FunctionalFactory<std::decay_t<decltype(nBodyAssembler)>, decltype(nonlinearity), Matrix, Vector> functionalFactory(nBodyAssembler, nonlinearity, _matrices, _rhs);
+        functionalFactory.build();
+        auto functional = functionalFactory.functional();
+
+        NonlinearSolver<std::remove_reference_t<decltype(*functional)>, BitVector> solver(parset.sub("solver.tnnmg"), linearSolver, functional, _dirichletNodes);
+        solver.solve(_localParset, totalX);
+
+        nBodyAssembler.postprocess(totalX, _x);
+    };
+
+    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);
+
+    size_t dof = 0;
+    for (size_t i=0; i<bodyCount_; i++) {
+        const auto& body = *contactNetwork.body(i);
+
+        if (body.data()->getYoungModulus() == 0.0) {
+            for (size_t j=0; j<body.nVertices(); j++) {
+                dirichletNodes[dof] = true;
+                dof++;
+            }
+        } else {
+            dof += body.nVertices();
+        }
+    }
+
+    std::cout << "solving linear problem for u..." << std::endl;
+
+    solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, u,
+                       parset.sub("u0.solver"));
+
+    // 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());
+      ScalarVector nmWeights(frictionBoundaryStress);
+
+      body->assembler()->assembleWeightedNormalStress(
+        contactCoupling->nonmortarBoundary(), frictionBoundaryStress, nmWeights, body->data()->getYoungModulus(),
+        body->data()->getPoissonRatio(), u[nonmortarIdx]);
+
+      weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
+      weights[nonmortarIdx] += nmWeights;
+    }
+
+    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> weights;
+  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 
+  strikeslipfactory.hh
+  strikeslipfactory.cc
+  threeblocksfactory.hh
+  threeblocksfactory.cc 
+  twoblocksfactory.hh
+  twoblocksfactory.cc 
+)
+
+#install headers
+install(FILES
+  cantorfactory.hh
+  contactnetworkfactory.hh
+  levelcontactnetworkfactory.hh
+  stackedblocksfactory.hh
+  strikeslipfactory.hh
+  threeblocksfactory.hh
+  twoblocksfactory.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/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"

dune/tectonic/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/assemblers_tmpl.cc → dune/tectonic/factories/cantorfactory_tmpl.cc

@@ -2,6 +2,6 @@
 #error MY_DIM unset
 #endif
 
-#include "explicitgrid.hh"
+#include "../explicitgrid.hh"
 
-template class MyAssembler<GridView, MY_DIM>;
+template class CantorFactory<Grid, MY_DIM>;

dune/tectonic/factories/contactnetworkfactory.hh

+#ifndef SRC_CONTACTNETWORKFACTORY_HH
+#define SRC_CONTACTNETWORKFACTORY_HH
+
+#include <dune/common/parametertree.hh>
+
+#include "../data-structures/network/contactnetwork.hh"
+
+template <class HostGridType, class VectorType>
+class ContactNetworkFactory {
+public:
+    using ContactNetwork = ContactNetwork<HostGridType, VectorType>;
+
+protected:
+    using LeafBody = typename ContactNetwork::LeafBody;
+    using FrictionCouplingPair = typename ContactNetwork::FrictionCouplingPair;
+
+    const Dune::ParameterTree& parset_;
+    const size_t bodyCount_;
+    const size_t couplingCount_;
+
+    std::vector<std::shared_ptr<LeafBody>> bodies_;
+    std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
+
+    ContactNetwork contactNetwork_;
+
+private:
+    virtual void setBodies() = 0;
+    virtual void setLevelBodies() = 0;
+    virtual void setCouplings() = 0;
+    virtual void setLevelCouplings() = 0;
+    virtual void setBoundaryConditions() = 0;
+
+public:
+    ContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
+        parset_(parset),
+        bodyCount_(bodyCount),
+        couplingCount_(couplingCount),
+        bodies_(bodyCount_),
+        couplings_(couplingCount_),
+        contactNetwork_(bodyCount_, couplingCount_) {}
+
+    void build() {
+        setBodies();
+        contactNetwork_.setBodies(bodies_);
+
+        setCouplings();
+        contactNetwork_.setCouplings(couplings_);
+
+        setLevelBodies();
+        setLevelCouplings();
+
+        setBoundaryConditions();      
+    }
+
+    ContactNetwork& contactNetwork() {
+        return contactNetwork_;
+    }
+};
+#endif

dune/tectonic/factories/levelcontactnetworkfactory.hh

+#ifndef SRC_CONTACTNETWORKFACTORY_HH
+#define SRC_CONTACTNETWORKFACTORY_HH
+
+#include <dune/common/parametertree.hh>
+
+#include "../data-structures/levelcontactnetwork.hh"
+
+template <class HostGridType>
+class ContactNetworkFactory {
+public:
+    using LevelContactNetwork = LevelContactNetwork<GridType, dims>;
+
+protected:
+    using Body = typename LevelContactNetwork::Body;
+    using FrictionCouplingPair = typename LevelContactNetwork::FrictionCouplingPair;
+
+    const Dune::ParameterTree& parset_;
+    const size_t bodyCount_;
+    const size_t couplingCount_;
+
+    std::vector<std::shared_ptr<Body>> bodies_;
+    std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
+
+    LevelContactNetwork levelContactNetwork_;
+
+private:
+    virtual void setBodies() = 0;
+    virtual void setCouplings() = 0;
+    virtual void setBoundaryConditions() = 0;
+
+public:
+    LevelContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
+        parset_(parset),
+        bodyCount_(bodyCount),
+        couplingCount_(couplingCount),
+        bodies_(bodyCount_),
+        couplings_(couplingCount_),
+        levelContactNetwork_(bodyCount_, couplingCount_) {}
+
+    void build() {
+        setBodies();
+        levelContactNetwork_.setBodies(bodies_);
+
+        setCouplings();
+        levelContactNetwork_.setCouplings(couplings_);
+
+        setBoundaryConditions();
+    }
+
+    LevelContactNetwork& levelContactNetwork() {
+        return levelContactNetwork_;
+    }
+};
+#endif

dune/tectonic/factories/stackedblocksfactory.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../problem-data/bc.hh"
+#include "../problem-data/grid/cuboidgeometry.hh"
+#include "../problem-data/myglobalfrictiondata.hh"
+
+#include "../utils/diameter.hh"
+
+#include "stackedblocksfactory.hh"
+
+
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {
+
+    // set up cuboid geometries
+    double const length = 1.00;
+    double const height = 0.27;
+    double const weakLength = 0.60;
+
+    std::vector<GlobalCoords> origins(this->bodyCount_);
+
+    const auto& subParset = this->parset_.sub("boundary.friction.weakening");
+
+#if MY_DIM == 3
+        double const depth = 0.60;
+
+        auto weakBound = [&] (const GlobalCoords& origin) {
+            GlobalCoords h = {origin[0] + weakLength, origin[1], origin[2]};
+            return h;
+        };
+
+        origins[0] = {0,0,0};
+        GlobalCoords lowerWeakOrigin = {0.2, 0, 0};
+        GlobalCoords upperWeakOrigin = {0.2, height, 0};
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height, depth);
+        cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));
+
+        for (size_t i=1; i<this->bodyCount_-1; i++) {
+            origins[i] = cuboidGeometries_[i-1]->upperLeft();
+            GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
+            GlobalCoords upperWeakOrigin_i = upperWeakOrigin + origins[i];
+
+            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height, depth);
+            cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
+            cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
+        }
+
+        const size_t idx = this->bodyCount_-1;
+        origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
+        lowerWeakOrigin += origins[idx];
+
+        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height, depth);
+        cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
+#elif MY_DIM == 2
+        auto weakBound = [&] (const GlobalCoords& origin) {
+            GlobalCoords h = {origin[0] + weakLength, origin[1]};
+            return h;
+        };
+
+        origins[0] = {0,0};
+        GlobalCoords lowerWeakOrigin = {0.2, 0};
+        GlobalCoords upperWeakOrigin = {0.2, height};
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height);
+        cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));
+
+        for (size_t i=1; i<this->bodyCount_-1; i++) {
+            origins[i] = cuboidGeometries_[i-1]->upperLeft();
+            auto height_i = height/3.0;
+
+            GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
+            GlobalCoords upperWeakOrigin_i = {upperWeakOrigin[0], height_i};
+            upperWeakOrigin_i += origins[i];
+
+            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height_i);
+
+            cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
+            cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));  
+        }
+
+        const size_t idx = this->bodyCount_-1;
+        origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
+        lowerWeakOrigin += origins[idx];
+
+        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height);
+        cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
+#else
+#error CuboidGeometry only supports 2D and 3D!"
+#endif
+
+    // set up reference grids
+    gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(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
+        const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
+        for (size_t j=0; j<weakeningRegions.size(); j++) {
+            refine(*grids[i], weakeningRegions[j], this->parset_.template get<double>("boundary.friction.smallestDiameter"), CuboidGeometry::lengthScale());
+        }
+
+        // determine minDiameter and maxDiameter
+        double minDiameter = std::numeric_limits<double>::infinity();
+        double maxDiameter = 0.0;
+        for (auto &&e : elements(grids[i]->leafGridView())) {
+          auto const geometry = e.geometry();
+          auto const diam = diameter(geometry);
+          minDiameter = std::min(minDiameter, diam);
+          maxDiameter = std::max(maxDiameter, diam);
+        }
+        std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
+        std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
+    }
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        this->bodies_[i] = std::make_shared<typename Base::LeafBody>(bodyData_, grids[i]);
+    }
+}
+
+
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setCouplings() {
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& cuboidGeometry = *cuboidGeometries_[i];
+
+        leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
+        levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
+    }
+
+    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+    std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;
+
+
+    for (size_t i=0; i<this->couplings_.size(); i++) {
+      auto& coupling = this->couplings_[i];
+      coupling = std::make_shared<typename Base::FrictionCouplingPair>();
+
+      nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
+      mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
+      weakPatches_[i] = std::make_shared<WeakeningRegion>(cuboidGeometries_[i+1]->weakeningRegions()[0]);
+
+      coupling->set(i+1, i, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
+      coupling->setWeakeningPatch(weakPatches_[i]);
+      coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[i], CuboidGeometry::lengthScale()));
+    }
+}
+
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setLevelBodies() {
+    const size_t nBodies = this->bodyCount_;
+
+    size_t maxLevel = 0;
+    std::vector<size_t> maxLevels(nBodies);
+    for (size_t i=0; i<nBodies; i++) {
+        const size_t bodyMaxLevel = this->bodies_[i]->grid()->maxLevel();
+        maxLevels[i] = bodyMaxLevel;
+        maxLevel = std::max(maxLevel, bodyMaxLevel);
+    }
+
+    for (int l=maxLevel; l>=0; l--) {
+        this->contactNetwork_.addLevel(maxLevels, l);
+
+        for (size_t i=0; i<nBodies; i++) {
+            maxLevels[i]--;
+        }
+    }
+}
+
+
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setBoundaryConditions() {
+    using LeafBoundaryCondition = BoundaryCondition<LeafGridView, 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>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"));
+
+    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->nVertices();
+
+        // Neumann boundary
+        std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->gridView()), neumannFunction, "neumann");
+        body->addBoundaryCondition(neumannBoundary);
+
+        // upper Dirichlet Boundary
+        // identify Dirichlet nodes on leaf level
+        if (i==this->bodyCount_-1) {
+            std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
+            for (int j=0; j<leafVertexCount; j++) {
+                if (leafFaces_[i]->upper.containsVertex(j))
+                    (*velocityDirichletNodes)[j][0] = true;
+            }
+
+            std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+
+            velocityDirichletBoundary->setBoundaryPatch(body->gridView(), velocityDirichletNodes);
+            velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+            body->addBoundaryCondition(velocityDirichletBoundary);
+        } else {
+
+        }
+    }
+
+    // lower Dirichlet Boundary
+    const auto& firstBody = this->contactNetwork_.body(0);
+    const auto& firstLeafVertexCount = firstBody->nVertices();
+    std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(firstLeafVertexCount);
+    for (int j=0; j<firstLeafVertexCount; j++) {
+        if (leafFaces_[0]->lower.containsVertex(j)) {
+            for (size_t d=0; d<dim; d++) {
+              (*zeroDirichletNodes)[j][d] = true;
+            }
+        }
+
+        #if MY_DIM == 3 //TODO: wrong, needs revision
+        if (leafFaces_[0]->front.containsVertex(j) || leafFaces_[0]->back.containsVertex(j))
+            (*zeroDirichletNodes)[j][2] = true;
+        #endif
+    }
+
+    std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+    zeroDirichletBoundary->setBoundaryPatch(firstBody->gridView(), zeroDirichletNodes);
+
+    std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
+    zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
+    firstBody->addBoundaryCondition(zeroDirichletBoundary);
+}
+
+#include "stackedblocksfactory_tmpl.cc"

dune/tectonic/factories/stackedblocksfactory.hh

+#ifndef SRC_STACKEDBLOCKSFACTORY_HH
+#define SRC_STACKEDBLOCKSFACTORY_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 "../problem-data/mybody.hh"
+#include "../problem-data/grid/mygrids.hh"
+#include "../problem-data/grid/cuboidgeometry.hh"
+
+template <class HostGridType, class VectorType> class StackedBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
+private:
+    using Base = ContactNetworkFactory<HostGridType, VectorType>;
+
+public:
+    using ContactNetwork = typename Base::ContactNetwork;
+
+private:
+    using GlobalCoords = typename ContactNetwork::LocalVector;
+
+    using LeafGridView = typename ContactNetwork::GridView;
+    using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
+
+    using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+    using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
+
+    using LeafFaces = MyFaces<LeafGridView>;
+    using LevelFaces = MyFaces<LevelGridView>;
+
+    using CuboidGeometry= CuboidGeometry<typename GlobalCoords::field_type>;
+    using WeakeningRegion = typename CuboidGeometry::WeakeningRegion;
+
+    static const int dim = ContactNetwork::dim;
+
+    const std::shared_ptr<MyBodyData<dim>> bodyData_;     // material properties of bodies
+
+    std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
+
+    std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
+
+    std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
+    std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
+
+    std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
+
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
+
+public:
+    StackedBlocksFactory(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_.sub("body"), this->parset_.template get<double>("gravity"), zenith_())),
+        cuboidGeometries_(this->bodyCount_),
+        leafFaces_(this->bodyCount_),
+        levelFaces_(this->bodyCount_),
+        weakPatches_(this->bodyCount_),
+        nonmortarPatch_(this->couplingCount_),
+        mortarPatch_(this->couplingCount_)
+    {}
+
+    void setBodies();
+    void setLevelBodies();
+    void setCouplings();
+    void setLevelCouplings() {}
+    void setBoundaryConditions();
+
+    auto& weakPatches() {
+        return weakPatches_;
+    }
+
+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
+
+