src/data-structures/body.hh

+#ifndef SRC_DATA_STRUCTURES_BODY_HH
+#define SRC_DATA_STRUCTURES_BODY_HH
+
+#include <functional>
+
+#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 GridTEMPLATE, class VectorTEMPLATE, int dims>
+class Body {
+public:
+    using GridType = GridTEMPLATE;
+    using VectorType = VectorTEMPLATE;
+
+    using DeformationFunction = DeformationFunction<typename GridType::LeafGridView, VectorType>;
+    using DeformedGridType = Dune::GeometryGrid<GridType, DeformationFunction>;
+    using DeformedLeafGridView = typename DeformedGridType::LeafGridView;
+    using DeformedLevelGridView = typename DeformedGridType::LevelGridView;
+
+    using Assembler = MyAssembler<DeformedLeafGridView, dims>;
+    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 Function = Dune::VirtualFunction<double, double>;
+    using ExternalForce = std::function<void(double, VectorType &)>;
+
+    using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;   
+    using LevelBoundaryCondition = BoundaryCondition<DeformedLevelGridView, dims>;
+
+    using BoundaryFunctions = std::vector<const Function* >;
+    using BoundaryNodes = std::vector<const Dune::BitSetVector<dims>* >;
+    using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
+    using BoundaryPatches = std::vector<const typename LeafBoundaryCondition::BoundaryPatch* >;
+
+    using StateEnergyNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+
+private:
+    std::shared_ptr<Body<GridType, VectorType, dims>> parent_;
+    std::vector<std::shared_ptr<Body<GridType, VectorType, dims>>> children_;
+
+    std::shared_ptr<BodyData<dims>> bodyData_;
+
+    std::shared_ptr<GridType> grid_;
+    std::shared_ptr<DeformationFunction> deformation_;
+    std::shared_ptr<DeformedGridType> deformedGrid_;
+
+    DeformedLeafGridView leafView_;
+    int leafVertexCount_;
+
+    std::shared_ptr<Assembler> assembler_;
+    Matrices<Matrix,1> matrices_;
+
+    ExternalForce externalForce_;
+
+    VectorType gravityFunctional_;
+
+    // boundary conditions
+    std::vector<std::shared_ptr<LeafBoundaryCondition>> leafBoundaryConditions_;
+    std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>> levelBoundaryConditions_; // first index: level, second index: bc on lvl
+
+public:
+    Body(
+            const std::shared_ptr<BodyData<dims>>& bodyData,
+            const std::shared_ptr<GridType>& grid);
+
+    void assemble();
+
+    // setter and getter
+    auto data() const -> const std::shared_ptr<BodyData<dims>>&;
+
+    void setDeformation(const VectorType& def);
+    auto deformation() const -> DeformationFunction&;
+    auto deformedGrid() const -> const DeformedGridType*;
+
+    auto leafView() const -> const DeformedLeafGridView&;
+    auto leafVertexCount() const -> int;
+
+    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<LeafBoundaryCondition> leafBoundaryCondition);
+    void leafBoundaryConditions(
+            const std::string& tag,
+            std::vector<std::shared_ptr<LeafBoundaryCondition>>& selectedConditions) const;
+    auto leafBoundaryConditions() const -> const std::vector<std::shared_ptr<LeafBoundaryCondition>>&;
+
+    void addBoundaryCondition(std::shared_ptr<LevelBoundaryCondition> levelBoundaryCondition, int level);
+    void levelBoundaryConditions(
+            const std::string& tag,
+            std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& selectedConditions) const;
+    void levelBoundaryConditions(
+            const std::string& tag,
+            int level,
+            std::vector<std::shared_ptr<LevelBoundaryCondition>>& selectedConditions) const;
+    auto levelBoundaryConditions(int level) const -> const std::vector<std::shared_ptr<LevelBoundaryCondition>>&;
+    auto levelBoundaryConditions() const -> const std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>&;
+
+    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/hdf5/surface-writer_tmpl.cc → src/data-structures/body_tmpl.cc

-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
 
-using MyProgramState = ProgramState<Vector, ScalarVector>;
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
 
-template class SurfaceWriter<MyProgramState, GridView>;
+template class Body<Grid, Vector, MY_DIM>;

src/data-structures/contactnetwork.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+
+#include "contactnetwork.hh"
+
+template <class GridType, int dim>
+void ContactNetwork<GridType, dim>::assemble() {
+    for (size_t i=0; i<levelContactNetworks_.size(); i++) {
+        levelContactNetworks_[i]->assemble();
+    }
+}
+
+template <class GridType, int dim>
+void ContactNetwork<GridType, dim>::addLevelContactNetwork(std::shared_ptr<LevelContactNetwork> levelContactNetwork) {
+    size_t level = levelContactNetwork->level();
+    assert(level < this->size());
+
+    levelContactNetworks_[level] = levelContactNetwork;
+}
+
+template <class GridType, int dim>
+const std::shared_ptr<typename ContactNetwork<GridType, dim>::LevelContactNetwork>& ContactNetwork<GridType, dim>::addLevelContactNetwork(int nBodies, int nCouplings, int level) {
+    assert(level < this->size());
+
+    levelContactNetworks_[level] = std::make_shared<LevelContactNetwork>(nBodies, nCouplings, level);
+    return levelContactNetworks_[level];
+}
+
+#include "contactnetwork_tmpl.cc"

src/data-structures/contactnetwork.hh

+#ifndef SRC_CONTACTNETWORK_HH
+#define SRC_CONTACTNETWORK_HH
+
+#include "levelcontactnetwork.hh"
+
+template <class GridType, int dim> class ContactNetwork {
+  public:
+    using LevelContactNetwork = LevelContactNetwork<GridType, dim>;
+
+  private:
+   std::vector<std::shared_ptr<LevelContactNetwork>> levelContactNetworks_;
+
+  public:
+        ContactNetwork() {}
+
+        void addLevelContactNetwork(std::shared_ptr<LevelContactNetwork> levelContactNetwork);
+        const std::shared_ptr<LevelContactNetwork>& addLevelContactNetwork(int nBodies, int nCouplings, int level);
+
+        void assemble();
+
+        const std::shared_ptr<LevelContactNetwork>& levelContactNetwork(size_t level) const {
+            return levelContactNetworks_[level];
+        }
+
+        void resize(size_t size) {
+            levelContactNetworks_.resize(size);
+        }
+
+        size_t size() const {
+            return levelContactNetworks_.size();
+        }
+
+        size_t maxLevel() const {
+            return size()-1;
+		} 
+};
+#endif

src/data-structures/contactnetwork_tmpl.cc

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

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/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, int dimension>
+LevelContactNetwork<GridType, dimension>::LevelContactNetwork(
+        int nBodies,
+        int nCouplings,
+        int level) :
+    level_(level),
+    bodies_(nBodies),
+    couplings_(nCouplings),
+    deformedGrids_(nBodies),
+    leafViews_(nBodies),
+    levelViews_(nBodies),
+    leafVertexCounts_(nBodies),
+    frictionBoundaries_(nCouplings),
+    stateEnergyNorms_(nBodies),
+    nBodyAssembler_(nBodies, nCouplings)
+{}
+
+template <class GridType, int dimension>
+LevelContactNetwork<GridType, dimension>::~LevelContactNetwork() {
+    for (size_t i=0; i<nBodies(); i++) {
+        delete frictionBoundaries_[i];
+        delete frictionBoundaryMass_[i];
+        delete stateEnergyNorms_[i];
+    }
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::assemble() {
+    using BoundaryPatch = DeformedLeafBoundaryPatch;
+
+    for (size_t i=0; i<nBodies(); i++) {
+        bodies_[i]->assemble();
+        //printBasisDofLocation(bodies_[i]->assembler()->vertexBasis); //TODO remove after debugging
+
+        frictionBoundaries_[i] = new BoundaryPatch(bodies_[i]->leafView(), false);
+    }
+
+    // set up dune-contact nBodyAssembler
+    nBodyAssembler_.setGrids(deformedGrids_);
+
+    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] = new ScalarMatrix();
+        bodies_[i]->assembler()->assembleFrictionalBoundaryMass(*frictionBoundaries_[i], *frictionBoundaryMass_[i]);
+        *frictionBoundaryMass_[i] /= frictionBoundaries_[i]->area(); // TODO: weight by individual friction patches?
+
+        stateEnergyNorms_[i] = new EnergyNorm<ScalarMatrix, ScalarVector>(*frictionBoundaryMass_[i]);
+    }
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::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<DeformedGridType, 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, Vector, TruncatedRateState, DeformedGridType>>(
+              nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
+        case Config::Regularised:
+          globalFriction_ = std::make_shared<GlobalRateStateFriction<
+              Matrix, Vector, RegularisedRateState, DeformedGridType>>(
+              nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
+        default:
+          assert(false);
+          break;
+      }
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::setBodies(const std::vector<std::shared_ptr<Body>> 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++) {
+        deformedGrids_[i] = bodies_[i]->deformedGrid();
+
+        leafViews_[i] = std::make_unique<DeformedLeafGridView>(deformedGrids_[i]->leafGridView());
+
+        levelViews_[i].resize(deformedGrids_[i]->maxLevel()+1);
+        for (size_t j=0; j<levelViews_[i].size(); j++) {
+            levelViews_[i][j] = std::make_unique<DeformedLevelGridView>(deformedGrids_[i]->levelGridView(j));
+        }
+
+        leafVertexCounts_[i] = leafViews_[i]->size(dimension);
+
+        matrices_.elasticity[i] = bodies_[i]->matrices().elasticity;
+        matrices_.damping[i] = bodies_[i]->matrices().damping;
+        matrices_.mass[i] = bodies_[i]->matrices().mass;
+    }
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings) {
+    assert(this->nCouplings()==couplings.size());
+    couplings_ = couplings;
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::constructBody(
+        const std::shared_ptr<BodyData<dimension>>& bodyData,
+        const std::shared_ptr<GridType>& grid,
+        std::shared_ptr<Body>& body) const {
+
+    body = std::make_shared<Body>(bodyData, grid);
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::constructCoupling(
+        int nonMortarBodyIdx, int mortarBodyIdx,
+        const std::shared_ptr<DeformedLevelBoundaryPatch>& nonmortarPatch,
+        const std::shared_ptr<DeformedLevelBoundaryPatch>& mortarPatch,
+        const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
+        const std::shared_ptr<GlobalFrictionData<dimension>>& frictionData,
+        std::shared_ptr<FrictionCouplingPair>& coupling) const {
+
+    coupling = std::make_shared<FrictionCouplingPair>();
+
+    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<DeformedLeafBoundaryPatch>>();
+    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);
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::totalNodes(
+        const std::string& tag,
+        Dune::BitSetVector<dimension>& totalNodes) const {
+
+    int totalSize = 0;
+    for (size_t i=0; i<nBodies(); i++) {
+        totalSize += this->body(i)->leafVertexCount();
+    }
+
+    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::LeafBoundaryCondition>> boundaryConditions;
+        body->leafBoundaryConditions(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<dimension; 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, int dimension>
+void LevelContactNetwork<GridType, dimension>::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, int dimension>
+void LevelContactNetwork<GridType, dimension>::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, int dimension>
+void LevelContactNetwork<GridType, dimension>::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, int dimension>
+void LevelContactNetwork<GridType, dimension>::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, int dimension>
+void LevelContactNetwork<GridType, dimension>::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 GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::stateEnergyNorms() const
+-> const StateEnergyNorms& {
+
+    return stateEnergyNorms_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::level() const
+-> int {
+
+    return level_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::nBodies() const
+-> size_t {
+
+    return bodies_.size();
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::nCouplings() const
+-> size_t {
+
+    return couplings_.size();
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::leafVertexCounts() const
+-> const std::vector<size_t>& {
+
+    return leafVertexCounts_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::body(int i) const
+-> const std::shared_ptr<Body>& {
+
+    return bodies_[i];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::coupling(int i) const
+-> const std::shared_ptr<FrictionCouplingPair>& {
+
+    return couplings_[i];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::deformedGrids() const
+-> const std::vector<const DeformedGridType*>& {
+
+    return deformedGrids_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::couplings() const
+-> const std::vector<std::shared_ptr<FrictionCouplingPair>>& {
+
+    return couplings_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::leafView(int bodyID) const
+-> const DeformedLeafGridView& {
+
+    return *leafViews_[bodyID];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::levelViews(int bodyID) const
+-> const std::vector<std::unique_ptr<const DeformedLevelGridView>>& {
+
+    return levelViews_[bodyID];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::levelView(int bodyID, int level) const
+-> const DeformedLevelGridView& {
+
+    return *levelViews_[bodyID][level];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::leafVertexCount(int bodyID) const
+-> int {
+
+    return leafVertexCounts_[bodyID];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::nBodyAssembler()
+-> Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& {
+
+    return nBodyAssembler_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::nBodyAssembler() const
+-> const Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& {
+
+    return nBodyAssembler_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::matrices() const
+-> const Matrices<Matrix,2>& {
+
+    return matrices_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::globalFriction() const
+-> GlobalFriction& {
+
+    return *globalFriction_;
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::externalForces(ExternalForces& externalForces) const {
+    externalForces.resize(nBodies());
+
+    for (size_t i=0; i<nBodies(); i++) {
+        externalForces[i] = &bodies_[i]->externalForce();
+    }
+}
+
+#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/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/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 GridType, int dimension> class LevelContactNetwork {
+    private:
+        using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
+
+    public:
+        using Body = Body<GridType, Vector, dimension>;
+
+        using DeformedGridType = typename Body::DeformedGridType;
+        using DeformedLeafGridView = typename Body::DeformedLeafGridView;
+        using DeformedLevelGridView = typename Body::DeformedLevelGridView;
+
+        using DeformedLeafBoundaryPatch = BoundaryPatch<DeformedLeafGridView>;
+        using DeformedLevelBoundaryPatch = BoundaryPatch<DeformedLevelGridView>;
+
+        using Assembler = typename Body::Assembler;
+        using Matrix = typename Assembler::Matrix;
+        using LocalVector = typename Vector::block_type;
+        using ScalarMatrix = typename Assembler::ScalarMatrix;
+        using ScalarVector = typename Assembler::ScalarVector;
+        using field_type = typename Matrix::field_type;
+    
+        using FrictionCouplingPair = FrictionCouplingPair<DeformedGridType, LocalVector, field_type>;
+    
+        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 GlobalFriction = GlobalFriction<Matrix, Vector>;
+
+        using StateEnergyNorms = std::vector<const typename Body::StateEnergyNorm *>;
+        using ExternalForces = std::vector<const typename Body::ExternalForce *>;
+
+        using NBodyAssembler = Dune::Contact::NBodyAssembler<DeformedGridType, Vector>;
+
+    public:
+        LevelContactNetwork(int nBodies, int nCouplings, int level = 0);
+        ~LevelContactNetwork();
+
+        void assemble();
+        void assembleFriction(
+                const Config::FrictionModel& frictionModel,
+                const std::vector<ScalarVector>& weightedNormalStress);
+		
+        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<dimension>>& bodyData,
+                const std::shared_ptr<GridType>& grid,
+                std::shared_ptr<Body>& body) const;
+        void constructCoupling(
+                int nonMortarBodyIdx,
+                int mortarBodyIdx,
+                const std::shared_ptr<DeformedLevelBoundaryPatch>& nonmortarPatch,
+                const std::shared_ptr<DeformedLevelBoundaryPatch>& mortarPatch,
+                const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
+                const std::shared_ptr<GlobalFrictionData<dimension>>& frictionData,
+                std::shared_ptr<FrictionCouplingPair>& coupling) const;
+
+        // getter
+        void totalNodes(
+                const std::string& tag,
+                Dune::BitSetVector<dimension>& 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;
+
+        auto stateEnergyNorms() const -> const StateEnergyNorms&;
+        auto globalFriction() const -> GlobalFriction&;
+        void externalForces(ExternalForces& externalForces) 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 deformedGrids() const -> const std::vector<const DeformedGridType*>& ;
+
+        auto leafVertexCounts() const -> const std::vector<size_t>&;
+        auto leafVertexCount(int bodyID) const -> int;
+
+        auto leafView(int bodyID) const -> const DeformedLeafGridView&;
+        auto levelView(int bodyID, int level = 0) const -> const DeformedLevelGridView&;
+        auto levelViews(int bodyID) const -> const std::vector<std::unique_ptr<const DeformedLevelGridView>>&;
+
+        auto nBodyAssembler() -> Dune::Contact::NBodyAssembler<DeformedGridType, Vector>&;
+        auto nBodyAssembler() const -> const Dune::Contact::NBodyAssembler<DeformedGridType, Vector>&;
+
+        auto matrices() const -> const Matrices<Matrix,2>&;
+
+    private:
+        const int level_;
+
+        std::vector<std::shared_ptr<Body>> bodies_;
+        std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
+
+        std::vector<const DeformedGridType*> deformedGrids_;
+        std::vector<std::unique_ptr<const DeformedLeafGridView>> leafViews_;
+        std::vector<std::vector<std::unique_ptr<const DeformedLevelGridView>>> levelViews_;
+        std::vector<size_t> leafVertexCounts_;
+
+        std::vector<DeformedLeafBoundaryPatch*> frictionBoundaries_; // union of all boundary patches per body
+        std::vector<ScalarMatrix*> frictionBoundaryMass_;
+        StateEnergyNorms stateEnergyNorms_;
+
+        NBodyAssembler nBodyAssembler_;
+
+        Matrices<Matrix,2> matrices_;
+
+        std::shared_ptr<GlobalFriction> globalFriction_;
+};
+#endif

src/data-structures/levelcontactnetwork_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "levelcontactnetwork.hh"
+
+using MyLevelContactNetwork = LevelContactNetwork<Grid, MY_DIM>;
+
+template class LevelContactNetwork<Grid, MY_DIM>;

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/tectonic/bodydata.hh>
+
+#include "../assemblers.hh"
+#include "levelcontactnetwork.hh"
+#include "matrices.hh"
+#include "../spatial-solving/solverfactory.hh"
+#include "../spatial-solving/solverfactory_old.hh"
+#include "../spatial-solving/tnnmg/functional.hh"
+#include "../spatial-solving/tnnmg/zerononlinearity.hh"
+#include "../utils/debugutils.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;
+
+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 GridType>
+  void setupInitialConditions(const Dune::ParameterTree& parset, const LevelContactNetwork<GridType, dims>& levelContactNetwork) {
+    using Matrix = typename LevelContactNetwork<GridType, dims>::Matrix;
+    const auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
+
+    // Solving a linear problem with a multigrid solver
+    auto const solveLinearProblem = [&](
+        const Dune::BitSetVector<dims>& _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 and TNMMG solver
+      using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
+      Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper);
+
+      using Factory = SolverFactory<Functional, Dune::BitSetVector<dims>>;
+      Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
+
+   /*   std::vector<Dune::BitSetVector<dims>> 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;
+
+      levelContactNetwork.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
+    Dune::BitSetVector<dims> dirichletNodes;
+    levelContactNetwork.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;
+        }
+    }*/
+
+    solveLinearProblem(dirichletNodes, levelContactNetwork.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 = levelContactNetwork.body(i);
+      std::vector<std::shared_ptr<typename LevelContactNetwork<GridType, dims>::Body::LeafBoundaryCondition>> frictionBoundaryConditions;
+      body->leafBoundaryConditions("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]);
+    }
+
+    Dune::BitSetVector<dims> noNodes(dirichletNodes.size(), false);
+    solveLinearProblem(noNodes, levelContactNetwork.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,6 +3,7 @@
 
 #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>
 
@@ -12,4 +13,6 @@ 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 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>;