src/one-body-problem-data/mybody.hh → dune/tectonic/problem-data/mybody.hh

-#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYBODY_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_MYBODY_HH
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYBODYDATA_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_MYBODYDATA_HH
 
 #include <dune/common/fvector.hh>
 
 #include <dune/fufem/functions/constantfunction.hh>
 
-#include <dune/tectonic/body.hh>
-#include <dune/tectonic/gravity.hh>
-
-#include "mygeometry.hh"
+#include "../data-structures/body/bodydata.hh"
+#include "gravity.hh"
+#include "grid/cuboidgeometry.hh"
 #include "segmented-function.hh"
 
-template <int dimension> class MyBody : public Body<dimension> {
-  using typename Body<dimension>::ScalarFunction;
-  using typename Body<dimension>::VectorField;
+template <int dimension> class MyBodyData : public BodyData<dimension> {
+  using typename BodyData<dimension>::ScalarFunction;
+  using typename BodyData<dimension>::VectorField;
 
 public:
-  MyBody(Dune::ParameterTree const &parset)
-      : poissonRatio_(parset.get<double>("body.poissonRatio")),
-        youngModulus_(3.0 * parset.get<double>("body.bulkModulus") *
+  MyBodyData(Dune::ParameterTree const &parset, const double gravity, const Dune::FieldVector<double, dimension>& zenith)
+      : poissonRatio_(parset.get<double>("poissonRatio")),
+        youngModulus_(3.0 * parset.get<double>("bulkModulus") *
                       (1.0 - 2.0 * poissonRatio_)),
+        zenith_(zenith),
         shearViscosityField_(
-            parset.get<double>("body.elastic.shearViscosity"),
-            parset.get<double>("body.viscoelastic.shearViscosity")),
+            parset.get<double>("elastic.shearViscosity"),
+            parset.get<double>("viscoelastic.shearViscosity")),
         bulkViscosityField_(
-            parset.get<double>("body.elastic.bulkViscosity"),
-            parset.get<double>("body.viscoelastic.bulkViscosity")),
-        densityField_(parset.get<double>("body.elastic.density"),
-                      parset.get<double>("body.viscoelastic.density")),
-        gravityField_(densityField_, MyGeometry::zenith,
-                      parset.get<double>("gravity")) {}
+            parset.get<double>("elastic.bulkViscosity"),
+            parset.get<double>("viscoelastic.bulkViscosity")),
+        densityField_(parset.get<double>("elastic.density"),
+                      parset.get<double>("viscoelastic.density")),
+        gravityField_(densityField_, zenith_,
+                      gravity) {}
 
   double getPoissonRatio() const override { return poissonRatio_; }
   double getYoungModulus() const override { return youngModulus_; }
@@ -47,6 +47,8 @@ template <int dimension> class MyBody : public Body<dimension> {
 private:
   double const poissonRatio_;
   double const youngModulus_;
+  Dune::FieldVector<double, dimension> zenith_;
+
   SegmentedFunction const shearViscosityField_;
   SegmentedFunction const bulkViscosityField_;
   SegmentedFunction const densityField_;

src/one-body-problem-data/myglobalfrictiondata.hh → dune/tectonic/problem-data/myglobalfrictiondata.hh

-#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
 
 #include <dune/common/function.hh>
 
-#include <dune/tectonic/globalfrictiondata.hh>
+#include "../data-structures/friction/globalfrictiondata.hh"
 
 #include "patchfunction.hh"
 
@@ -14,14 +14,15 @@ class MyGlobalFrictionData : public GlobalFrictionData<LocalVector::dimension> {
 
 public:
   MyGlobalFrictionData(Dune::ParameterTree const &parset,
-                       ConvexPolyhedron<LocalVector> const &segment)
+                       ConvexPolyhedron<LocalVector> const &segment,
+                       double lengthScale = 1.0)
       : C_(parset.get<double>("C")),
         L_(parset.get<double>("L")),
         V0_(parset.get<double>("V0")),
         a_(parset.get<double>("strengthening.a"),
-           parset.get<double>("weakening.a"), segment),
+           parset.get<double>("weakening.a"), segment, lengthScale),
         b_(parset.get<double>("strengthening.b"),
-           parset.get<double>("weakening.b"), segment),
+           parset.get<double>("weakening.b"), segment, lengthScale),
         mu0_(parset.get<double>("mu0")) {}
 
   double const &C() const override { return C_; }

src/one-body-problem-data/patchfunction.hh → dune/tectonic/problem-data/patchfunction.hh

-#ifndef SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
 
 #include <dune/common/function.hh>
 #include <dune/common/fvector.hh>
@@ -17,13 +17,15 @@ class PatchFunction
   double const v2_;
   Polyhedron const &segment_;
 
+  double const lengthScale_;
+
 public:
-  PatchFunction(double v1, double v2, Polyhedron const &segment)
-      : v1_(v1), v2_(v2), segment_(segment) {}
+  PatchFunction(double v1, double v2, Polyhedron const &segment, double lengthScale = 1.0)
+      : v1_(v1), v2_(v2), segment_(segment), lengthScale_(lengthScale) {}
 
   void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
                 Dune::FieldVector<double, 1> &y) const {
-    y = distance(x, segment_, 1e-6 * MyGeometry::lengthScale) <= 1e-5 ? v2_
+    y = distance(x, segment_, 1e-6 * lengthScale_) <= 1e-5 ? v2_
                                                                       : v1_;
   }
 };

src/one-body-problem-data/segmented-function.hh → dune/tectonic/problem-data/segmented-function.hh

@@ -5,7 +5,7 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/parametertree.hh>
 
-#include "mygeometry.hh"
+#include "grid/cuboidgeometry.hh"
 
 class SegmentedFunction
     : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
@@ -15,10 +15,10 @@ class SegmentedFunction
                  Dune::FieldVector<double, MY_DIM> const &y,
                  Dune::FieldVector<double, MY_DIM> const &z) const {
     return x[1] + (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1];
-  };
+  }
   bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {
-    return liesBelow(MyGeometry::K, MyGeometry::M, z);
-  };
+    return false; //TODO liesBelow(MyGeometry::K, MyGeometry::M, z);
+  }
 
   double const _v1;
   double const _v2;

dune/tectonic/problem-data/strikeslipbody.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYBODYDATA_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_MYBODYDATA_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include "../data-structures/body/bodydata.hh"
+#include "gravity.hh"
+#include "grid/cuboidgeometry.hh"
+
+
+class SegmentedFunction
+    : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
+                                   Dune::FieldVector<double, 1>> {
+private:
+  bool distFromDiagonal(const Dune::FieldVector<double, MY_DIM>& z) const {
+    return std::abs(z[0] - z[1])/std::sqrt(2);
+  }
+
+  bool insideViscousRegion(Dune::FieldVector<double, MY_DIM> const &z) const {
+    return distFromDiagonal(z) > distFromDiag_;
+  }
+
+  const double distFromDiag_;
+
+  const double v1_;
+  const double v2_;
+
+public:
+  SegmentedFunction(double v1, double v2, double distFromDiag = 0.1) :
+      distFromDiag_(distFromDiag) , v1_(v1), v2_(v2) {}
+
+  void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
+                Dune::FieldVector<double, 1> &y) const {
+    y = insideViscousRegion(x) ? v2_ : v1_;
+  }
+};
+
+
+template <int dimension> class MyBodyData : public BodyData<dimension> {
+  using typename BodyData<dimension>::ScalarFunction;
+  using typename BodyData<dimension>::VectorField;
+
+public:
+  MyBodyData(Dune::ParameterTree const &parset, const double gravity, const Dune::FieldVector<double, dimension>& zenith)
+      : poissonRatio_(parset.get<double>("poissonRatio")),
+        youngModulus_(3.0 * parset.get<double>("bulkModulus") *
+                      (1.0 - 2.0 * poissonRatio_)),
+        zenith_(zenith),
+        shearViscosityField_(
+            parset.get<double>("elastic.shearViscosity"),
+            parset.get<double>("viscoelastic.shearViscosity"),
+            parset.get<double>("elastic.distFromDiag")),
+        bulkViscosityField_(
+            parset.get<double>("elastic.bulkViscosity"),
+            parset.get<double>("viscoelastic.bulkViscosity"),
+            parset.get<double>("elastic.distFromDiag")),
+        densityField_(parset.get<double>("elastic.density"),
+                      parset.get<double>("viscoelastic.density"),
+                      parset.get<double>("elastic.distFromDiag")),
+        gravityField_(densityField_, zenith_,
+                      gravity) {}
+
+  double getPoissonRatio() const override { return poissonRatio_; }
+  double getYoungModulus() const override { return youngModulus_; }
+  ScalarFunction const &getShearViscosityField() const override {
+    return shearViscosityField_;
+  }
+  ScalarFunction const &getBulkViscosityField() const override {
+    return bulkViscosityField_;
+  }
+  ScalarFunction const &getDensityField() const override {
+    return densityField_;
+  }
+  VectorField const &getGravityField() const override { return gravityField_; }
+
+private:
+  double const poissonRatio_;
+  double const youngModulus_;
+  Dune::FieldVector<double, dimension> zenith_;
+
+  SegmentedFunction const shearViscosityField_;
+  SegmentedFunction const bulkViscosityField_;
+  SegmentedFunction const densityField_;
+  Gravity<dimension> const gravityField_;
+};
+#endif

dune/tectonic/quadraticenergy.hh

-#ifndef DUNE_TECTONIC_QUADRATICENERGY_HH
-#define DUNE_TECTONIC_QUADRATICENERGY_HH
-
-#include <memory>
-
-template <class NonlinearityTEMPLATE> class QuadraticEnergy {
-public:
-  using Nonlinearity = NonlinearityTEMPLATE;
-  using LocalVector = typename Nonlinearity::VectorType;
-
-  QuadraticEnergy(double alpha, LocalVector const &b, Nonlinearity const &phi)
-      : alpha(alpha), b(b), phi(phi) {}
-
-  double const alpha;
-  LocalVector const &b;
-  Nonlinearity const &phi;
-};
-#endif

dune/tectonic/spatial-solving/CMakeLists.txt

+add_subdirectory("contact")
+add_subdirectory("tnnmg")
+add_subdirectory("preconditioners")
+
+add_custom_target(tectonic_dune_spatial-solving SOURCES
+  fixedpointiterator.hh
+  fixedpointiterator.cc
+  functionalfactory.hh
+  makelinearsolver.hh
+  nonlinearsolver.hh
+  solverfactory.hh
+  solverfactory.cc
+)
+
+#install headers
+install(FILES
+  fixedpointiterator.hh
+  functionalfactory.hh.hh
+  makelinearsolver.hh
+  nonlinearsolver.hh
+  solverfactory.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/spatial-solving/contact/CMakeLists.txt

+add_custom_target(tectonic_dune_spatial-solving_contact SOURCES
+  dualmortarcoupling.hh
+  dualmortarcoupling.cc
+  nbodyassembler.hh
+  nbodyassembler.cc
+  nbodycontacttransfer.hh
+  nbodycontacttransfer.cc
+)
+
+#install headers
+install(FILES
+  dualmortarcoupling.hh
+  nbodyassembler.hh
+  nbodycontacttransfer.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#include <dune/common/exceptions.hh>
+
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+
+#include <dune/grid-glue/merging/contactmerge.hh>
+//#include <dune/grid-glue/adapter/gridgluevtkwriter.hh>
+//#include <dune/grid-glue/adapter/gridglueamirawriter.hh>
+
+#include <dune/contact/common/dualbasisadapter.hh>
+#include <dune/matrix-vector/addtodiagonal.hh>
+
+template <class field_type, class GridType0, class GridType1>
+void DualMortarCoupling<field_type, GridType0, GridType1>::reducePatches() {
+    typedef typename GridType0::LeafGridView GridView0;
+    typedef typename GridType1::LeafGridView GridView1;
+
+    using Element0 = typename GridView0::template Codim<0>::Entity;
+    using Element1 = typename GridView1::template Codim<0>::Entity;
+
+    auto desc0 = [&] (const Element0& e, unsigned int face) {
+        return nonmortarBoundary_.contains(e,face);
+    };
+
+    auto desc1 = [&] (const Element1& e, unsigned int face) {
+        return mortarBoundary_.contains(e,face);
+    };
+
+    GridView0 gridView0 = grid0_->leafGridView();
+    GridView1 gridView1 = grid1_->leafGridView();
+
+    auto extract0 = std::make_shared<Extractor0>(gridView0,desc0);
+    auto extract1 = std::make_shared<Extractor1>(gridView1,desc1);
+
+    if (!gridGlueBackend_)
+        gridGlueBackend_ = std::make_shared< Dune::GridGlue::ContactMerge<dimworld, ctype> >(overlap_);
+    glue_ = std::make_shared<Glue>(extract0, extract1, gridGlueBackend_);
+    auto& glue = *glue_;
+    glue.build();
+
+    //std::cout << glue.size() << " remote intersections found." << std::endl;
+    //GridGlueAmiraWriter::write<GlueType>(glue,debugPath_);
+
+    // Restrict the hasObstacle fields to the part of the nonmortar boundary
+    // where we could actually create a contact mapping
+    LeafBoundaryPatch0 boundaryWithMapping(gridView0);
+
+    const auto& indexSet0 = gridView0.indexSet();
+    const auto& indexSet1 = gridView1.indexSet();
+
+    ///////////////////////////////////
+    //  reducing nonmortar boundary
+    /////////////////////////////////
+
+    // Get all fine grid boundary segments that are totally covered by the grid-glue segments
+    typedef std::pair<int,int> Pair;
+    std::map<Pair,ctype> coveredArea, fullArea;
+
+    // initialize with area of boundary faces
+    for (const auto& bIt : nonmortarBoundary_) {
+        const Pair p(indexSet0.index(bIt.inside()),bIt.indexInInside());
+        fullArea[p] = bIt.geometry().volume();
+        coveredArea[p] = 0;
+    }
+
+    // sum up the remote intersection areas to find out which are totally covered
+    for (const auto& rIs : intersections(glue))
+        coveredArea[Pair(indexSet0.index(rIs.inside()),rIs.indexInInside())] += rIs.geometry().volume();
+
+    // add all fine grid faces that are totally covered by the contact mapping
+    for (const auto& bIt : nonmortarBoundary_) {
+        const auto& inside = bIt.inside();
+        if(coveredArea[Pair(indexSet0.index(inside),bIt.indexInInside())]/
+            fullArea[Pair(indexSet0.index(inside),bIt.indexInInside())] >= coveredArea_)
+            boundaryWithMapping.addFace(inside, bIt.indexInInside());
+    }
+
+    //writeBoundary(boundaryWithMapping,debugPath_ + "relevantNonmortar");
+
+
+    /** \todo replace by all fine grid segments which are totally covered by the RemoteIntersections. */
+    //for (const auto& rIs : intersections(glue))
+    //    boundaryWithMapping.addFace(rIs.inside(),rIs.indexInInside());
+
+    /*printf("contact mapping could be built for %d of %d boundary segments.\n",
+           boundaryWithMapping.numFaces(), nonmortarBoundary_.numFaces()); */
+
+    nonmortarBoundary_ = boundaryWithMapping;
+    mortarBoundary_.setup(gridView1);
+    for (const auto& rIs : intersections(glue))
+        if (nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+            mortarBoundary_.addFace(rIs.outside(),rIs.indexInOutside());
+}
+
+template <class field_type, class GridType0, class GridType1>
+void DualMortarCoupling<field_type, GridType0, GridType1>::setCrosspoints(const std::set<size_t>& allCrosspoints) {
+    crosspoints_.clear();
+
+    if (!allCrosspoints.empty()) {
+        Dune::BitSetVector<1> nonmortarPatchBoundary;
+        nonmortarBoundary_.getPatchBoundaryVertices(nonmortarPatchBoundary);
+
+        for (auto p : allCrosspoints) {
+            if (nonmortarPatchBoundary[p][0]) {
+                crosspoints_.emplace(p);
+            }
+        }
+    }
+}
+
+template <class field_type, class GridType0, class GridType1>
+void DualMortarCoupling<field_type, GridType0, GridType1>::assembleNonmortarLagrangeMatrix()
+{  
+    // clear matrix
+    const auto nNonmortarVertices = nonmortarBoundary_.numVertices(); // - crosspoints_.size();
+    nonmortarLagrangeMatrix_ = Dune::BDMatrix<MatrixBlock>(nNonmortarVertices);
+    nonmortarLagrangeMatrix_ = 0;
+
+    const auto& indexSet = grid0_->leafIndexSet();
+
+    // loop over all faces of the boundary patch
+    for (const auto& nIt : nonmortarBoundary_) {
+
+        const auto& geometry = nIt.geometry();
+        const field_type numCorners = geometry.corners();
+
+        ctype intElem = geometry.integrationElement(Dune::FieldVector<ctype,dim-1>(0));
+
+        field_type sfI = (numCorners==3) ? intElem/6.0 : intElem/numCorners;
+
+        // turn scalar element mass matrix into vector-valued one
+        // and add element mass matrix to the global matrix
+
+        // Get global node ids
+        const auto& inside = nIt.inside();
+        const auto& refElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+
+        // Add localMatrix to nonmortarLagrangeMat
+        for (int i=0; i<refElement.size(nIt.indexInInside(),1,dim); i++) {
+            // we can use subEntity here because we add all indices anyway
+            int v = globalToLocal_[indexSet.subIndex(inside,refElement.subEntity(nIt.indexInInside(),1,
+                                                                          i,dim),dim)];
+            if (v == -1)
+               continue;
+
+            Dune::MatrixVector::addToDiagonal(nonmortarLagrangeMatrix_[v][v],sfI);
+        }
+
+    }
+
+}
+
+template <class field_type, class GridType0, class GridType1>
+void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
+{
+    typedef typename GridType0::LeafGridView GridView0;
+    typedef typename GridType1::LeafGridView GridView1;
+
+    using FiniteElementCache1 = Dune::PQkLocalFiniteElementCache<typename GridType1::ctype, field_type, GridType1::dimension,1>;
+
+    // cache for the dual functions on the boundary
+    using DualCache = Dune::Contact::DualBasisAdapter<GridView0, field_type>;
+    using ConstantDualCache = Dune::PQkLocalFiniteElementCache<typename GridType0::ctype, field_type, GridType0::dimension, 0>;
+
+    GridView0 gridView0 = grid0_->leafGridView();
+    GridView1 gridView1 = grid1_->leafGridView();
+
+    const auto& indexSet0 = gridView0.indexSet();
+    const auto& indexSet1 = gridView1.indexSet();
+
+    // Create mapping from the global set of block dofs to the ones on the contact boundary
+    const auto& vertices = *nonmortarBoundary_.getVertices();
+
+    globalToLocal_.resize(vertices.size());
+    hasObstacle_.resize(vertices.size());
+    int idx = 0;
+    for (size_t i=0; i<globalToLocal_.size(); i++) {
+        bool nonmortarVertex = vertices[i][0]; // and crosspoints_.count(i)==0;
+        globalToLocal_[i] = nonmortarVertex ? idx++ : -1;
+        hasObstacle_[i][0] = nonmortarVertex ? true : false;
+    }
+
+    // Assemble the diagonal matrix coupling the nonmortar side with the lagrange multiplyers there
+    assembleNonmortarLagrangeMatrix();
+
+    // The weak obstacle vector
+    const auto nNonmortarVertices = nonmortarBoundary_.numVertices();// - crosspoints_.size();
+    weakObstacle_.resize(nNonmortarVertices);
+    weakObstacle_ = 0;
+
+    // ///////////////////////////////////////////////////////////
+    //   Get the occupation structure for the mortar matrix
+    // ///////////////////////////////////////////////////////////
+
+    /** \todo Also restrict mortar indices and don't use the whole grid level. */
+    Dune::MatrixIndexSet mortarIndices(nNonmortarVertices, grid1_->size(dim));
+
+    // loop over all intersections
+    for (const auto& rIs : intersections(*glue_)) {
+
+        if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+            continue;
+
+        const auto& inside = rIs.inside();
+        const auto& outside = rIs.outside();
+
+        const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+        const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+
+        int nDomainVertices = domainRefElement.size(dim);
+        int nTargetVertices = targetRefElement.size(dim);
+
+        for (int j=0; j<nDomainVertices; j++) {
+
+            int localDomainIdx = globalToLocal_[indexSet0.subIndex(inside,j,dim)];
+
+            // if the vertex is not contained in the restricted contact boundary or a crosspoint, then dismiss it
+            if (localDomainIdx == -1)
+                continue;
+
+            for (int k=0; k<nTargetVertices; k++) {
+                int globalTargetIdx = indexSet1.subIndex(outside,k,dim);
+                if (!mortarBoundary_.containsVertex(globalTargetIdx))
+                    continue;
+
+                mortarIndices.add(localDomainIdx, globalTargetIdx);
+            }
+        }
+    }
+
+    mortarIndices.exportIdx(mortarLagrangeMatrix_);
+
+    // Clear it
+    mortarLagrangeMatrix_ = 0;
+
+    //cache of local bases
+    FiniteElementCache1 cache1;
+
+    std::unique_ptr<DualCache> dualCache;
+    dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView0, field_type> >();
+
+    // adapt dual basis close to crosspoints using ConstantDualCache
+    ConstantDualCache constantDualCache;
+
+    std::vector<Dune::FieldVector<ctype,dim> > avNormals;
+    avNormals = nonmortarBoundary_.getNormals();
+
+    // loop over all intersections and compute the matrix entriescrosspoints_
+    for (const auto& rIs : intersections(*glue_)) {
+        const auto& inside = rIs.inside();
+
+        if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+            continue;
+
+        const auto& outside = rIs.outside();
+
+        // types of the elements supporting the boundary segments in question
+        Dune::GeometryType nonmortarEType = inside.type();
+        Dune::GeometryType mortarEType    = outside.type();
+
+        const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(nonmortarEType);
+
+        const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(mortarEType);
+
+        int noOfMortarVec = targetRefElement.size(dim);
+
+        Dune::GeometryType nmFaceType = domainRefElement.type(rIs.indexInInside(),1);
+        Dune::GeometryType mFaceType  = targetRefElement.type(rIs.indexInOutside(),1);
+
+        // Select a quadrature rule
+        // 2 in 2d and for integration over triangles in 3d.  If one (or both) of the two faces involved
+        // are quadrilaterals, then the quad order has to be risen to 3 (4).
+        int quadOrder = 2 + (!nmFaceType.isSimplex()) + (!mFaceType.isSimplex());
+        const auto& quadRule = Dune::QuadratureRules<ctype, dim-1>::rule(rIs.type(), quadOrder);
+
+        //const typename FiniteElementCache0::FiniteElementType& nonmortarFiniteElement = cache0.get(nonmortarEType);
+        const auto& mortarFiniteElement = cache1.get(mortarEType);
+
+        const auto& constantDualFE = constantDualCache.get(nonmortarEType);
+        dualCache->bind(inside, rIs.indexInInside());
+
+        std::vector<Dune::FieldVector<field_type,1> > mortarQuadValues, dualQuadValues;
+
+        const auto& rGeom = rIs.geometry();
+        const auto& rGeomOutside = rIs.geometryOutside();
+        const auto& rGeomInInside = rIs.geometryInInside();
+        const auto& rGeomInOutside = rIs.geometryInOutside();
+
+        bool isCrosspointFace = false;
+        std::vector<int> nonmortarFaceNodes;
+        for (int i=0; i<domainRefElement.size(rIs.indexInInside(),1,dim); i++) {
+          int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
+          int globalIdx = indexSet0.subIndex(inside,faceIdxi,dim);
+
+          /*if (crosspoints_.count(globalIdx)>0) {
+            isCrosspointFace = true;
+          } else {*/
+            nonmortarFaceNodes.push_back(faceIdxi);
+          //}
+        }
+
+
+
+        for (const auto& quadPt : quadRule) {
+
+            // compute integration element of overlap
+            ctype integrationElement = rGeom.integrationElement(quadPt.position());
+
+            // quadrature point positions on the reference element
+            Dune::FieldVector<ctype,dim> nonmortarQuadPos = rGeomInInside.global(quadPt.position());
+            Dune::FieldVector<ctype,dim> mortarQuadPos    = rGeomInOutside.global(quadPt.position());
+
+            // The current quadrature point in world coordinates
+            Dune::FieldVector<field_type,dim> nonmortarQpWorld = rGeom.global(quadPt.position());
+            Dune::FieldVector<field_type,dim> mortarQpWorld    = rGeomOutside.global(quadPt.position());;
+
+            // the gap direction (normal * gapValue)
+            Dune::FieldVector<field_type,dim> gapVector = mortarQpWorld  - nonmortarQpWorld;
+
+            //evaluate all shapefunctions at the quadrature point
+            //nonmortarFiniteElement.localBasis().evaluateFunction(nonmortarQuadPos,nonmortarQuadValues);
+            mortarFiniteElement.localBasis().evaluateFunction(mortarQuadPos,mortarQuadValues);
+
+            /*if (isCrosspointFace) {
+              constantDualFE.localBasis().evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+              // loop over all Lagrange multiplier shape functions
+              for (size_t j=0; j<nonmortarFaceNodes.size(); j++) {
+
+                  int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+                  int rowIdx = globalToLocal_[globalDomainIdx];
+
+                  weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
+                      * dualQuadValues[0] * (gapVector*avNormals[globalDomainIdx]);
+
+                  // loop over all mortar shape functions
+                  for (int k=0; k<noOfMortarVec; k++) {
+
+                      int colIdx  = indexSet1.subIndex(outside, k, dim);
+                      if (!mortarBoundary_.containsVertex(colIdx))
+                          continue;
+
+                      // Integrate over the product of two shape functions
+                      field_type mortarEntry =  integrationElement* quadPt.weight()* dualQuadValues[0]* mortarQuadValues[k];
+
+                      Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
+
+                  }
+
+               }
+            } else {*/
+              dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+              // loop over all Lagrange multiplier shape functions
+              for (size_t j=0; j<nonmortarFaceNodes.size(); j++) {
+
+                  int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+                  int rowIdx = globalToLocal_[globalDomainIdx];
+
+                  weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
+                      * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);
+
+                  // loop over all mortar shape functions
+                  for (int k=0; k<noOfMortarVec; k++) {
+
+                      int colIdx  = indexSet1.subIndex(outside, k, dim);
+                      if (!mortarBoundary_.containsVertex(colIdx))
+                          continue;
+
+                      // Integrate over the product of two shape functions
+                      field_type mortarEntry =  integrationElement* quadPt.weight()* dualQuadValues[nonmortarFaceNodes[j]]* mortarQuadValues[k];
+
+                      Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
+
+                  }
+
+               }
+            //}
+        }
+
+    }
+
+    // ///////////////////////////////////////
+    //    Compute M = D^{-1} \hat{M}
+    // ///////////////////////////////////////
+
+    Dune::BCRSMatrix<MatrixBlock>& M  = mortarLagrangeMatrix_;
+    Dune::BDMatrix<MatrixBlock>& D    = nonmortarLagrangeMatrix_;
+
+    // First compute D^{-1}
+    D.invert();
+
+    // Then the matrix product D^{-1} \hat{M}
+    for (auto rowIt = M.begin(); rowIt != M.end(); ++rowIt) {
+        const auto rowIndex = rowIt.index();
+        for (auto& entry : *rowIt)
+            entry.leftmultiply(D[rowIndex][rowIndex]);
+    }
+
+    // weakObstacles in transformed basis = D^{-1}*weakObstacle_
+    for(size_t rowIdx=0; rowIdx<weakObstacle_.size(); rowIdx++)
+        weakObstacle_[rowIdx] *= D[rowIdx][rowIdx][0][0];
+
+    gridGlueBackend_->clear();
+}

dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_DUAL_MORTAR_COUPLING_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_DUAL_MORTAR_COUPLING_HH
+
+#include <memory>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bdmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/boundarypatchprolongator.hh>
+
+#include <dune/grid-glue/gridglue.hh>
+#include <dune/grid-glue/extractors/codim1extractor.hh>
+#include <dune/grid-glue/merging/merger.hh>
+
+/** \brief Assembles the transfer operator for two-body contact.
+ * Extension of Dune::Contact::DualMortarCoupling by crosspoints.
+ * Note: so war only works in 2D!
+ */
+template<class field_type, class GridType0, class GridType1=GridType0>
+class DualMortarCoupling {
+
+    // /////////////////////
+    // private types
+    // /////////////////////
+
+    enum {dim = GridType0::dimension};
+
+    enum {dimworld = GridType0::dimensionworld};
+
+  protected:
+    typedef Dune::FieldMatrix<field_type, dim, dim> MatrixBlock;
+    typedef Dune::BCRSMatrix<MatrixBlock> MatrixType;
+    using ObstacleVector = Dune::BlockVector<Dune::FieldVector<field_type, 1> >;
+
+    typedef typename GridType0::ctype ctype;
+
+    typedef BoundaryPatch<typename GridType0::LeafGridView> LeafBoundaryPatch0;
+    typedef BoundaryPatch<typename GridType1::LeafGridView> LeafBoundaryPatch1;
+
+    typedef BoundaryPatch<typename GridType0::LevelGridView> LevelBoundaryPatch0;
+    typedef BoundaryPatch<typename GridType1::LevelGridView> LevelBoundaryPatch1;
+
+    using Extractor0 = Dune::GridGlue::Codim1Extractor<typename GridType0::LeafGridView>;
+    using Extractor1 = Dune::GridGlue::Codim1Extractor<typename GridType1::LeafGridView>;
+
+public:
+
+    using Glue = Dune::GridGlue::GridGlue<Extractor0, Extractor1>;
+
+    DualMortarCoupling(field_type overlap = 1e-2, field_type coveredArea = 1.0 - 1e-2)
+        : overlap_(overlap), gridGlueBackend_(nullptr), coveredArea_(coveredArea)
+    {}
+
+    DualMortarCoupling(const GridType0& grid0, const GridType1& grid1,
+                       field_type overlap = 1e-2, field_type coveredArea = 1.0 - 1e-2)
+        : grid0_(&grid0), grid1_(&grid1), overlap_(overlap),
+          gridGlueBackend_(nullptr), coveredArea_(coveredArea)
+    {}
+
+    virtual ~DualMortarCoupling() {
+        /* Nothing. */
+    }
+
+    void reducePatches();
+
+    void setCrosspoints(const std::set<size_t>& allCrosspoints);
+
+    /** \brief Sets up the contact coupling operator on the grid leaf level */
+    virtual void setup();
+
+    /** \brief Assemble nonmortar matrix on the leaf level. */
+    void assembleNonmortarLagrangeMatrix();
+
+    /** \brief Get the dual mortar coupling matrix of the leaf level. */
+    virtual const MatrixType& mortarLagrangeMatrix() const {
+        return mortarLagrangeMatrix_;
+    }
+
+    /** \brief Get the dual mortar coupling matrix of the leaf level. */
+    virtual MatrixType& mortarLagrangeMatrix() {
+        return mortarLagrangeMatrix_;
+    }
+
+    /** \brief Get the non-mortar matrix. */
+    const Dune::BDMatrix<MatrixBlock>& nonmortarLagrangeMatrix() const {
+        return nonmortarLagrangeMatrix_;
+    }
+
+    /** \brief Setup leaf nonmortar and mortar patches. */
+    virtual void setupContactPatch(const LevelBoundaryPatch0& coarseNonmortar,
+                           const LevelBoundaryPatch1& coarseMortar) {
+
+        BoundaryPatchProlongator<GridType0>::prolong(coarseNonmortar,nonmortarBoundary_);
+        BoundaryPatchProlongator<GridType1>::prolong(coarseMortar,mortarBoundary_);
+    }
+
+    /** \brief Return the leafnonmortar boundary. */
+    const LeafBoundaryPatch0& nonmortarBoundary() const {
+        return nonmortarBoundary_;
+    }
+
+    /** \brief Return the leafmortar boundary. */
+    const LeafBoundaryPatch1& mortarBoundary() const {
+        return mortarBoundary_;
+    }
+
+    /** \brief Set the non-mortar and mortar grid. */
+    void setGrids(const GridType0& grid0, const GridType1& grid1) {
+        grid0_ = &grid0;
+        grid1_ = &grid1;
+    }
+
+    /** \brief Get non-mortar grid. */
+    const GridType0* getGrid0() const {
+        return grid0_;
+    }
+
+    /** \brief Get mortar grid. */
+    const GridType1* getGrid1() const {
+        return grid1_;
+    }
+
+    /** \brief Set the percentage of covered area each face must at least have
+     *         by the grid-glue projection.
+     */
+    void setCoveredArea(field_type coveredArea) {
+        coveredArea_ = coveredArea;
+    }
+
+    void setOverlap(field_type overlap) {
+        overlap_ = overlap;
+    }
+
+    /** \brief Get the GridGlue. */
+    std::shared_ptr<Glue> getGlue() {
+        return glue_;
+    }
+
+    const ObstacleVector& getWeakObstacle() const {
+        return weakObstacle_;
+    }
+
+    const Dune::BitSetVector<1>* hasObstacle() const {
+        return &hasObstacle_;
+    }
+
+    const std::vector<int>& globalToLocal() const {
+        return globalToLocal_;
+    }
+
+   /* const auto& crosspoints() const {
+        return crosspoints_;
+    }*/
+
+    // /////////////////////////////////////////////
+    //   Data members
+    // /////////////////////////////////////////////
+
+    /** \brief The two grids involved in the two-body contact problem. */
+    const GridType0* grid0_;
+    const GridType1* grid1_;
+
+    /** \brief For each dof a bit specifying whether the dof carries an obstacle or not. */
+    LeafBoundaryPatch0 nonmortarBoundary_;
+
+    /** \brief The mortar boundary. */
+    LeafBoundaryPatch1 mortarBoundary_;
+
+    /** \brief The crosspoints of the nonmortar body. */
+    std::set<size_t> crosspoints_;
+
+    /** \brief Nonmortar boundary patch global to local map. */
+    std::vector<int> globalToLocal_;
+
+    /** \brief The matrix coupling mortar side and Lagrange multipliers. */
+    MatrixType mortarLagrangeMatrix_;
+
+    /** \brief The diagonal matrix coupling nonmortar side and Lagrange multipliers */
+    Dune::BDMatrix<MatrixBlock> nonmortarLagrangeMatrix_;
+
+    /** \brief The weak obstacles. */
+    ObstacleVector weakObstacle_;
+
+    Dune::BitSetVector<1> hasObstacle_;
+
+    /** \brief Allow some small penetration of the bodies. */
+    field_type overlap_;
+
+    /** \brief A grid-glue merger. If this is not set ContactMerge is used. */
+    std::shared_ptr< Dune::GridGlue::Merger<field_type,dim-1,dim-1,dim> > gridGlueBackend_;
+
+    // Path where to write the merged grids and the reduced contact boundaries
+    //std::string debugPath_;
+
+    std::shared_ptr<Glue> glue_;
+
+    /** \brief Dismiss all faces that are not at least covered by the grid-glue projection for this
+     *         much percentage ranging between one - for total coverage and zero for taking all faces.
+     */
+    field_type coveredArea_;
+};
+
+#include "dualmortarcoupling.cc"
+
+#endif

dune/tectonic/spatial-solving/contact/nbodyassembler.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#include <memory>
+
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/matrix-vector/transpose.hh>
+#include <dune/matrix-vector/blockmatrixview.hh>
+
+#include <dune/fufem/boundarypatchprolongator.hh>
+#include <dune/contact/projections/normalprojection.hh>
+
+
+#include "../../utils/debugutils.hh"
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::setCrosspoints() {
+    /*std::vector<Dune::BitSetVector<1>> bodywiseBoundaries(nGrids());
+    for (size_t i=0; i<nGrids(); i++) {
+        bodywiseBoundaries[i] = *dirichletVertices_[i];
+    }
+
+    for (size_t i=0; i<nCouplings(); i++) {
+        auto& coupling = coupling_[i];
+        const auto& contactCoupling = contactCoupling_[i]; // dual mortar object holding boundary patches
+        const auto nonmortarIdx = coupling.gridIdx_[0];
+        //const auto mortarIdx = coupling->gridIdx_[1];
+
+        auto& bodywiseBoundary = bodywiseBoundaries[nonmortarIdx];
+        auto& crosspoints = crosspoints_[nonmortarIdx];
+        const auto& nmBoundaryVertices = *contactCoupling->nonmortarBoundary().getVertices();
+        for (size_t j=0; j<nmBoundaryVertices.size(); j++) {
+            if (bodywiseBoundary[j][0] and nmBoundaryVertices[j][0]) {
+                crosspoints.emplace(j);
+            }
+            bodywiseBoundary[j][0] = bodywiseBoundary[j][0] or nmBoundaryVertices[j][0];
+        }
+
+    }*/
+
+    for (int i=0; i<nCouplings(); i++) {
+        auto& coupling = coupling_[i];
+        const auto& contactCoupling = contactCoupling_[i]; // dual mortar object holding boundary patches
+        const auto nonmortarIdx = coupling.gridIdx_[0];
+
+        auto& crosspoints = crosspoints_[nonmortarIdx];
+
+        Dune::BitSetVector<1> nonmortarPatchBoundary;
+        contactCoupling->nonmortarBoundary().getPatchBoundaryVertices(nonmortarPatchBoundary);
+
+        for (size_t j=0; j<nonmortarPatchBoundary.size(); j++) {
+            if (nonmortarPatchBoundary[j][0]) {
+                crosspoints.emplace(j);
+            }
+        }
+
+    }
+
+
+    /*for (size_t i=0; i<crosspoints_.size(); i++) {
+        print(crosspoints_[i], "crosspoints " + std::to_string(i));
+    }*/
+}
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::assembleTransferOperator()
+{
+
+    // /////////////////////////////////////////////////////////////////
+    //   Check if contact data is present
+    // /////////////////////////////////////////////////////////////////
+
+    for (int i=0; i<nCouplings(); i++) {
+
+
+         if (!coupling_[i].obsPatch_)
+             DUNE_THROW(Dune::Exception, "You have to supply a nonmortar patch for the " << i << "-th coupling!");
+
+         if (!coupling_[i].mortarPatch_)
+             DUNE_THROW(Dune::Exception, "You have to supply a mortar patch for the " << i << "-th coupling!");
+    }
+
+    // ////////////////////////////////////////////////////
+    //   Set up Mortar element transfer operators
+    // ///////////////////////////////////////////////////
+
+    //std::cout<<"Setup mortar transfer operators\n";
+
+    for (size_t i=0; i<contactCoupling_.size(); i++) {
+        contactCoupling_[i]->setGrids(*grids_[coupling_[i].gridIdx_[0]],*grids_[coupling_[i].gridIdx_[1]]);
+        contactCoupling_[i]->setupContactPatch(*coupling_[i].patch0(),*coupling_[i].patch1());
+        contactCoupling_[i]->gridGlueBackend_ = coupling_[i].backend();
+        contactCoupling_[i]->setCoveredArea(coveredArea_);
+        contactCoupling_[i]->reducePatches();
+    }
+
+    //setCrosspoints();
+
+    /*for (int i=0; i<nGrids(); i++) {
+        print(crosspoints_[i], "crosspoints " + std::to_string(i));
+    }*/
+
+    for (size_t i=0; i<contactCoupling_.size(); i++) {
+        contactCoupling_[i]->setCrosspoints(crosspoints_[coupling_[i].gridIdx_[0]]);
+        contactCoupling_[i]->setup();
+
+        //print(contactCoupling_[i]->nonmortarLagrangeMatrix() , "nonmortarLagrangeMatrix " + std::to_string(i) + ": ");
+        //print(contactCoupling_[i]->mortarLagrangeMatrix() , "mortarLagrangeMatrix " + std::to_string(i) + ": ");
+    }
+
+    // setup Householder reflections
+    assembleHouseholderReflections();
+
+    // compute the mortar transformation matrix
+    computeTransformationMatrix();
+}
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::assembleHouseholderReflections()
+{
+    // //////////////////////////////////////////////////////////////////
+    //   Compute local coordinate systems for the dofs with obstacles
+    // //////////////////////////////////////////////////////////////////
+
+    // first canonical basis vector
+    using CoordinateType = Dune::FieldVector<field_type,dim>;
+    CoordinateType e0(0);
+    e0[0] = 1;
+
+    // local identity
+    Dune::ScaledIdentityMatrix<field_type,dim> id(1);
+
+    // bodywise local coordinate systems: initialise with identity
+    std::vector<std::vector<MatrixBlock> > coordinateSystems(nGrids());
+    for (size_t i=0; i<coordinateSystems.size(); i++) {
+        coordinateSystems[i].resize(grids_[i]->size(dim));
+        for (size_t j=0; j<coordinateSystems[i].size(); j++) {
+            coordinateSystems[i][j] = id;
+        }
+    }
+
+    std::vector<std::map<size_t, CoordinateType>> crosspointDirections(nGrids());
+
+    for (int i=0; i<nCouplings(); i++) {
+        auto& coupling = coupling_[i];
+
+        double dist = coupling.obsDistance_;
+        auto projection = coupling.projection();
+
+        if (!projection)
+            DUNE_THROW(Dune::Exception, "You have to supply a contact projection for the " << i << "-th coupling!");
+
+        std::vector<CoordinateType> directions;
+
+        const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
+        const auto& mortarBoundary    = contactCoupling_[i]->mortarBoundary();
+
+        projection->project(nonmortarBoundary, mortarBoundary,dist);
+        projection->getObstacleDirections(directions);
+
+        const auto nonmortarIdx = coupling_[i].gridIdx_[0];
+        auto globalToLocal = nonmortarBoundary.makeGlobalToLocal();
+
+        for (size_t j=0; j<globalToLocal.size(); j++) {
+            if (globalToLocal[j] > -1) {
+                // There is an obstacle at this vertex --> the coordinate system
+                // will consist of the surface normal and tangential vectors
+
+                if (crosspoints_[nonmortarIdx].count(j)>0) {
+                    crosspointDirections[nonmortarIdx][j] += directions[globalToLocal[j]];
+                } else {
+                    this->householderReflection(e0, directions[globalToLocal[j]], coordinateSystems[nonmortarIdx][j]);
+                }
+            }
+        }
+    }
+
+    // local coordinate systems at crosspoints
+    for (size_t i=0; i<crosspoints_.size(); i++) {
+        for (auto c : crosspoints_[i]) {
+            auto& direction = crosspointDirections[i][c];
+            direction /= direction.two_norm();
+
+            this->householderReflection(e0, direction, coordinateSystems[i][c]);
+        }
+    }
+
+    // ///////////////////////////////////////////////////////////////
+    // Combine the coordinate systems for each grid to one long array
+    // ///////////////////////////////////////////////////////////////
+    this->localCoordSystems_.resize(numDofs());
+
+    size_t offSet = 0;
+    for (size_t i=0; i<nGrids(); i++) {
+        for (std::size_t k=0; k<coordinateSystems[i].size(); k++) {
+            this->localCoordSystems_[offSet+k] = coordinateSystems[i][k];
+        }
+
+        offSet += grids_[i]->size(dim);
+    }
+
+    //print(this->localCoordSystems_, "localCoordSystems:");
+    //std::cout << "done" << std::endl;
+}
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::assembleObstacle()
+{
+
+    std::vector<std::vector<int> > globalToLocals(nCouplings());
+    for (std::size_t i = 0; i < globalToLocals.size(); ++i)
+        globalToLocals[i] = contactCoupling_[i]->globalToLocal();
+
+    ///////////////////////////////////////
+    //  Set the obstacle values
+    ///////////////////////////////////////
+
+    totalHasObstacle_.resize(numDofs(),false);
+
+    for (int j=0; j<nCouplings(); j++) {
+
+        int grid0Idx = coupling_[j].gridIdx_[0];
+        const auto& globalToLocal = globalToLocals[j];
+
+        // compute offset
+        int idx = 0;
+        for (int k=0;k<grid0Idx;k++)
+            idx += grids_[k]->size(dim);
+
+        for (int k=0; k<grids_[grid0Idx]->size(dim); k++)
+            if (globalToLocal[k] > -1)
+                totalHasObstacle_[idx+k] = true;
+    }
+
+    // //////////////////////////////////////////////////////////////////
+    //   Set the obstacle distances
+    // //////////////////////////////////////////////////////////////////
+
+    // Combine the obstacle values for each grid to one long array
+    totalObstacles_.resize(numDofs());
+
+    for (size_t j=0; j<totalObstacles_.size(); j++)
+        totalObstacles_[j].clear();
+
+    // Init the obstacle values
+    for (int i=0; i<nCouplings(); i++) {
+
+        int grid0Idx = coupling_[i].gridIdx_[0];
+
+        // compute offset
+        int idx = 0;
+        for (int k=0;k<grid0Idx;k++)
+            idx += grids_[k]->size(dim);
+
+        // The grids involved in this coupling
+        const auto& indexSet = grids_[grid0Idx]->leafIndexSet();
+
+        // the obstacles are stored using local indices
+        const std::vector<int>& globalToLocal = globalToLocals[i];
+
+        // the weak obstacles
+        const auto& obs = contactCoupling_[i]->weakObstacle_;
+
+        // get strong obstacles, the projection method of the correct level has already been called
+        //std::vector<field_type> obs;
+        //coupling_[i].projection()->getObstacles(obs);
+
+        const auto& leafView = grids_[grid0Idx]->leafGridView();
+
+        for (const auto& v : vertices(leafView)) {
+
+            int vIdx = indexSet.index(v);
+            if (globalToLocal[vIdx]<0) // and crosspoints_[grid0Idx].count(vIdx)<1)
+                continue;
+
+            // Set the obstacle
+            switch (coupling_[i].type_) {
+
+            case Dune::Contact::CouplingPairBase::CONTACT:
+                totalObstacles_[idx+vIdx].upper(0) = obs[globalToLocal[vIdx]];
+                break;
+
+            case Dune::Contact::CouplingPairBase::STICK_SLIP:
+                totalObstacles_[idx+vIdx].lower(0) = 0;
+                totalObstacles_[idx+vIdx].upper(0) = 0;
+
+                /*if (crosspoints_[grid0Idx].count(vIdx) > 0) {
+                    totalObstacles_[idx+vIdx].lower(1) = 0;
+                    totalObstacles_[idx+vIdx].upper(1) = 0;
+                }*/
+
+                break;
+
+            case Dune::Contact::CouplingPairBase::GLUE:
+                for (int j=0; j<dim; j++) {
+                    totalObstacles_[idx+vIdx].lower(j) = 0;
+                    totalObstacles_[idx+vIdx].upper(j) = 0;
+                }
+                break;
+
+            case Dune::Contact::CouplingPairBase::NONE:
+                break;
+            default:
+                DUNE_THROW(Dune::NotImplemented, "Coupling type not handled yet!");
+            }
+        }
+    }
+
+    //totalHasObstacle_[35] = true;
+    //totalObstacles_[35].lower(0) = 0;
+   // totalObstacles_[35].lower(1) = 0;
+    //totalObstacles_[35].upper(0) = 0;
+   // totalObstacles_[35].upper(1) = 0;
+}
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::
+assembleJacobian(const std::vector<const MatrixType*>& submat,
+                 MatrixType& totalMatrix) const
+{
+    // Create a block view of the grand matrix
+    Dune::MatrixVector::BlockMatrixView<MatrixType> blockView(submat);
+
+    int nRows = blockView.nRows();
+    int nCols = blockView.nCols();
+
+    // Create the index set of B \hat{A} B^T
+    Dune::MatrixIndexSet indicesBABT(nRows, nCols);
+
+    for (int i=0; i<nGrids(); i++) {
+
+        for (size_t iRow = 0; iRow<submat[i]->N(); iRow++) {
+
+            const RowType& row = (*submat[i])[iRow];
+
+            // Loop over all columns of the stiffness matrix
+            ConstColumnIterator j    = row.begin();
+            ConstColumnIterator jEnd = row.end();
+
+            for (; j!=jEnd; ++j) {
+
+                ConstColumnIterator k    = BT_[blockView.row(i,iRow)].begin();
+                ConstColumnIterator kEnd = BT_[blockView.row(i,iRow)].end();
+                for (; k!=kEnd; ++k) {
+
+                    ConstColumnIterator l    = BT_[blockView.col(i,j.index())].begin();
+                    ConstColumnIterator lEnd = BT_[blockView.col(i,j.index())].end();
+
+                    for (; l!=lEnd; ++l)
+                        indicesBABT.add(k.index(), l.index());
+
+                }
+
+            }
+
+        }
+
+    }
+
+    // ////////////////////////////////////////////////////////////////////
+    //   Multiply transformation matrix to the global stiffness matrix
+    // ////////////////////////////////////////////////////////////////////
+
+    indicesBABT.exportIdx(totalMatrix);
+    totalMatrix = 0;
+
+    for (int i=0; i<nGrids(); i++) {
+
+        for (size_t iRow = 0; iRow<submat[i]->N(); iRow++) {
+
+            const RowType& row = (*submat[i])[iRow];
+
+            // Loop over all columns of the stiffness matrix
+            ConstColumnIterator j    = row.begin();
+            ConstColumnIterator jEnd = row.end();
+
+            for (; j!=jEnd; ++j) {
+
+                ConstColumnIterator k    = BT_[blockView.row(i,iRow)].begin();
+                ConstColumnIterator kEnd = BT_[blockView.row(i,iRow)].end();
+
+                for (; k!=kEnd; ++k) {
+
+                    ConstColumnIterator l    = BT_[blockView.col(i,j.index())].begin();
+                    ConstColumnIterator lEnd = BT_[blockView.col(i,j.index())].end();
+
+                    for (; l!=lEnd; ++l) {
+
+                        //BABT[k][l] += BT[i][k] * mat_[i][j] * BT[j][l];
+                        MatrixBlock m_ij = *j;
+
+                        MatrixBlock BT_ik_trans;
+                        Dune::MatrixVector::transpose(*k, BT_ik_trans);
+
+                        m_ij.leftmultiply(BT_ik_trans);
+                        m_ij.rightmultiply(*l);
+
+                        totalMatrix[k.index()][l.index()] += m_ij;
+
+                    }
+
+                }
+
+            }
+
+        }
+
+    }
+}
+
+template <class GridType, class VectorType>
+template <class VectorTypeContainer>
+void NBodyAssembler<GridType, VectorType>::
+assembleRightHandSide(const VectorTypeContainer& rhs,
+               VectorType& totalRhs) const
+{
+    // Concatenate the two rhs vectors to a large one
+    VectorType untransformedTotalRhs(BT_.M());
+
+    int idx = 0;
+
+    for (size_t i=0; i<rhs.size(); i++) {
+
+        for (size_t j=0; j<rhs[i].size(); j++)
+            untransformedTotalRhs[idx++] = rhs[i][j];
+
+    }
+
+    if ((int) BT_.M() != idx)
+        DUNE_THROW(Dune::Exception, "assembleRightHandSide: vector size and matrix size don't match!");
+
+    // Transform the basis of the ansatz space
+    totalRhs.resize(untransformedTotalRhs.size());
+    totalRhs = 0;
+    BT_.umtv(untransformedTotalRhs, totalRhs);
+}
+
+
+template <class GridType, class VectorType>
+template <class VectorTypeContainer>
+void NBodyAssembler<GridType, VectorType>::
+postprocess(const VectorType& totalX, VectorTypeContainer& x) const
+{
+    // ///////////////////////////////////////////////////////
+    //   Transform the solution vector to the nodal basis
+    // ///////////////////////////////////////////////////////
+
+    VectorType nodalBasisTotalX(totalX.size());
+    BT_.mv(totalX, nodalBasisTotalX);
+
+
+    // ///////////////////////////////////////////////////////
+    //   Split the total solution vector into the parts
+    //   corresponding to the grids.
+    // ///////////////////////////////////////////////////////
+
+    int idx = 0;
+    for (int i=0; i<nGrids(); i++) {
+
+        x[i].resize(grids_[i]->size(dim));
+
+        for (size_t j=0; j<x[i].size(); j++, idx++)
+            x[i][j] = nodalBasisTotalX[idx];
+    }
+}
+
+template <class GridType, class VectorType>
+template <class VectorTypeContainer>
+void NBodyAssembler<GridType, VectorType>::
+concatenateVectors(const VectorTypeContainer& parts, VectorType& whole)
+{
+    int totalSize = 0;
+    for (size_t i=0; i<parts.size(); i++)
+        totalSize += parts[i].size();
+
+    whole.resize(totalSize);
+
+    int idx = 0;
+    for (size_t i=0; i<parts.size(); i++)
+        for (size_t j=0; j<parts[i].size(); j++)
+            whole[idx++] = parts[i][j];
+
+}
+
+// ////////////////////////////////////////////////////////////////////////////
+//   Turn the initial solutions from the nodal basis to the transformed basis
+// ////////////////////////////////////////////////////////////////////////////
+template <class GridType, class VectorType>
+template <class VectorTypeContainer>
+void NBodyAssembler<GridType, VectorType>::
+nodalToTransformed(const VectorTypeContainer& x,
+                   VectorType& transformedX) const
+{
+    VectorType canonicalTotalX;
+
+    concatenateVectors(x, canonicalTotalX);
+
+    // Make small cg solver
+    Dune::MatrixAdapter<MatrixType,VectorType,VectorType> op(getTransformationMatrix());
+
+    // A preconditioner
+    Dune::SeqILU<MatrixType,VectorType,VectorType> ilu0(getTransformationMatrix(),1.0);
+
+    // A cg solver for nonsymmetric matrices
+    Dune::BiCGSTABSolver<VectorType> bicgstab(op,ilu0,1E-4,100,0);
+
+    // Object storing some statistics about the solving process
+    Dune::InverseOperatorResult statistics;
+
+    // Solve!
+    transformedX = canonicalTotalX;  // seems to be a good initial value
+    bicgstab.apply(transformedX, canonicalTotalX, statistics);
+}
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::
+computeTransformationMatrix()
+{
+    //std::cout<<"Setup transformation matrix...";
+    // compute offsets for the different grids
+    std::vector<int> offsets(grids_.size());
+    offsets[0] = 0;
+
+    // P1 elements are hard-wired here
+    size_t k;
+    for (k=0; k<grids_.size()-1; k++)
+        offsets[k+1] = offsets[k] + grids_[k]->size(dim);
+
+    int nRows = offsets[k] + grids_[k]->size(dim);
+    int nCols = nRows;
+
+    // /////////////////////////////////////////////////////////////////
+    //   First create the index structure
+    // /////////////////////////////////////////////////////////////////
+
+    Dune::MatrixIndexSet indicesBT(nRows, nCols);
+
+    // BT_ is the identity plus some off-diagonal elements
+    for (size_t i=0; i<indicesBT.rows(); i++)
+        indicesBT.add(i,i);
+
+    std::vector<std::vector<int> > nonmortarToGlobal(nCouplings());
+    // Enter all the off-diagonal entries
+    for (int i=0; i<nCouplings(); i++) {
+
+        const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
+        const auto& globalToLocal = contactCoupling_[i]->globalToLocal();
+        // If the contact mapping could not be built at all then skip this coupling
+        if (nonmortarBoundary.numVertices() == 0)
+            continue;
+
+        // The grids involved in this coupling
+        int grid0Idx = coupling_[i].gridIdx_[0];
+        int grid1Idx = coupling_[i].gridIdx_[1];
+
+        // The mapping from nonmortar vertex indices to grid indices
+        nonmortarToGlobal[i].resize(nonmortarBoundary.numVertices());
+        int idx = 0;
+        for (int j=0; j<grids_[grid0Idx]->size(dim); j++)
+            if (globalToLocal[j] > -1)
+                nonmortarToGlobal[i][idx++] = j;
+        nonmortarToGlobal[i].resize(idx);
+
+        // the off-diagonal part
+        const MatrixType& M = contactCoupling_[i]->mortarLagrangeMatrix();
+        for (size_t j=0; j<M.N(); j++) {
+
+            const RowType& row = M[j];
+
+            ConstColumnIterator cIt    = row.begin();
+            ConstColumnIterator cEndIt = row.end();
+
+            for (; cIt!=cEndIt; ++cIt)
+                indicesBT.add(offsets[grid0Idx] + nonmortarToGlobal[i][j],
+                              offsets[grid1Idx] + cIt.index());
+
+        }
+
+    }
+
+    // ////////////////////////////////////////////////////////////
+    //   Enter the values of the different couplings
+    // ////////////////////////////////////////////////////////////
+
+    indicesBT.exportIdx(BT_);
+    BT_ = 0;
+
+    // Enter identity part
+    for (int i=0; i<nRows; i++)
+        for (int j=0; j<dim; j++)
+            BT_[i][i][j][j] = 1;
+
+    for (int i=0; i<nCouplings(); i++) {
+
+        const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
+        const auto& globalToLocal = contactCoupling_[i]->globalToLocal();
+        if (nonmortarBoundary.numVertices() == 0)
+            continue;
+
+        // The grids involved in this coupling
+        int grid0Idx = coupling_[i].gridIdx_[0];
+        int grid1Idx = coupling_[i].gridIdx_[1];
+
+        // the off-diagonal part
+        const MatrixType& M = contactCoupling_[i]->mortarLagrangeMatrix();
+        for (size_t j=0; j<M.N(); j++) {
+
+            const RowType& row = M[j];
+
+            ConstColumnIterator cIt    = row.begin();
+            ConstColumnIterator cEndIt = row.end();
+
+            for (; cIt!=cEndIt; ++cIt)
+                BT_[offsets[grid0Idx] + nonmortarToGlobal[i][j]][offsets[grid1Idx] +cIt.index()] = *cIt;
+
+        }
+
+        int offset = offsets[grid0Idx];
+
+        // modify non-mortar dofs and rotate them in normal and tangential coordinates
+        for (int k=0; k<grids_[grid0Idx]->size(dim); k++)
+            if(globalToLocal[k] > -1)
+                BT_[offset + k][offset+k] =this->localCoordSystems_[offset + k];
+    }
+}

dune/tectonic/spatial-solving/contact/nbodyassembler.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_N_BODY_MORTAR_ASSEMBLER_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_N_BODY_MORTAR_ASSEMBLER_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/promotiontraits.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/contact/assemblers/contactassembler.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/projections/contactprojection.hh>
+
+#include "dualmortarcoupling.hh"
+
+#ifdef HAVE_DUNE_GRID_GLUE
+#include <dune/grid-glue/merging/merger.hh>
+#endif
+
+
+/**  \brief Assembler for mortar discretized contact problems with arbitrary number of bodies.
+ *
+ *  \tparam GridType Type of the corresponding grids.
+ *  \tparam VectorType The vector type for the displacements.
+ */
+template <class GridType, class VectorType>
+class NBodyAssembler : public Dune::Contact::ContactAssembler<GridType::dimension, typename VectorType::field_type>
+{
+protected:
+    enum {dim = GridType::dimension};
+
+    using field_type = typename Dune::PromotionTraits<typename VectorType::field_type,
+                                                typename GridType::ctype>::PromotedType;
+
+    using CouplingType = DualMortarCoupling<field_type, GridType>;
+
+    using MatrixBlock = Dune::FieldMatrix<field_type, dim, dim>;
+    using MatrixType = Dune::BCRSMatrix<MatrixBlock>;
+
+    using RowType =  typename MatrixType::row_type;
+    using ConstColumnIterator = typename RowType::ConstIterator;
+
+    using LeafBoundaryPatch = BoundaryPatch<typename GridType::LeafGridView>;
+
+public:
+    /** \brief Construct assembler from number of couplings and grids.
+     *
+     *  \param nGrids The number of involved bodies.
+     *  \param nCouplings The number of couplings.
+     */
+    NBodyAssembler(int nGrids, int nCouplings, field_type coveredArea = 0.99) :
+      coveredArea_(coveredArea)
+    {
+        grids_.resize(nGrids);
+
+        crosspoints_.resize(nGrids);
+
+        coupling_.resize(nCouplings);
+        contactCoupling_.resize(nCouplings);
+        // initialise the couplings with DualMortarCoupling
+        for (int i=0; i<nCouplings; i++)
+            contactCoupling_[i] = std::make_shared<CouplingType>();
+    }
+
+    /** \brief Get the number of couplings.*/
+    int nCouplings() const {return coupling_.size();}
+
+    /** \brief Get the number of involved bodies.*/
+    int nGrids() const {return grids_.size();}
+
+    /** \brief Get the total number of degrees of freedom of the leaf grids. */
+    int numDofs() const {
+        int n=0;
+        for (int i=0; i<nGrids(); i++)
+            n += grids_[i]->size(dim);
+        return n;
+    }
+
+    void setCrosspoints();
+
+    /** \brief Setup the patches for the contact boundary and compute the obstacles. */
+    void assembleObstacle();
+
+    /** \brief Assemble the mortar matrices and compute the basis transformation.*/
+    void assembleTransferOperator();
+
+    /** \brief Turn initial solutions from nodal basis to the transformed basis.
+     *         i.e. transformedX = O*B^{-T}nodalX
+     *  \param x Initial solution vector containing nodal coefficients.
+     *  \param transformedX Coefficient vector for all bodies in mortar transformed coordinates.
+     */
+    template <class VectorTypeContainer>
+    void nodalToTransformed(const VectorTypeContainer& x,
+                            VectorType& transformedX) const;
+
+    /** \brief Compute stiffness matrices in mortar transformed coordinates.
+     *         i.e. transformedA = O*B*nodalA*B^T*O^T
+     *  \param submat Stiffness matrices w.r.t the nodal finite element basis.
+     *  \param totalMatrix Reference to mortar transformed stiffness matrix for all bodies.
+     */
+    void assembleJacobian(const std::vector<const MatrixType*>& submat,
+                          MatrixType& totalMatrix) const;
+
+    /** \brief Compute the right hand side in mortar transformed coordinates.
+     *         i.e. totalRhs = O*B*nodalRhs
+     *  \param  rhs Right hand side coefficients w.r.t. the nodal finite element basis.
+     *  \param  totalRhs Reference to coefficient vector w.r.t. the transformed basis.
+     */
+    template <class VectorTypeContainer>
+    void assembleRightHandSide(const VectorTypeContainer& rhs,
+                               VectorType& totalRhs) const;
+
+    /** \brief Transform a vector from local coordinates to canonical coordinates.
+     *
+     *  \param totalX Coefficient vector of the mortar transformed basis.
+     *  \param x    Reference to target vector for the standard nodal coefficients of each body.
+     */
+    template <class VectorTypeContainer>
+    void postprocess(const VectorType& totalX, VectorTypeContainer& x) const;
+
+    /** \brief Concatenate a family of vectors .
+     *
+     * \param parts A vector of vectors.
+     * \param whole The vector to contain the concatenated family
+     */
+    template <class VectorTypeContainer>
+    static void concatenateVectors(const VectorTypeContainer& parts, VectorType& whole);
+
+    /** \brief Get the contact couplings. */
+    const auto& getContactCouplings() const {
+        return contactCoupling_;
+    }
+
+    /** \brief Get the contact couplings. */
+    auto& getContactCouplings() {
+        return contactCoupling_;
+    }
+
+    /** \brief Set the contact couplings. */
+    void setContactCouplings(std::vector<std::shared_ptr<CouplingType> >& contactCouplings) {
+
+        contactCoupling_.resize(contactCouplings.size());
+        for (size_t i=0; i<contactCouplings.size(); i++)
+            contactCoupling_[i] = contactCouplings[i];
+    }
+
+    /** \brief Get the transposed of the mortar transformation matrix B^T.*/
+    const MatrixType& getTransformationMatrix() const {return BT_;}
+
+    /** \brief Get the grids. */
+    const std::vector<const GridType*> getGrids() const { return grids_; }
+
+    /** \brief Set the grids. */
+    void setGrids(const std::vector<const GridType*>& grids) {grids_ = grids;}
+
+    /** \brief Set grids. */
+    void setCoupling(const Dune::Contact::CouplingPair<GridType, GridType, field_type>& coupling, size_t i) {
+        coupling_[i] = coupling;
+    }
+
+    /** \brief Get reference to i'th coupling. */
+    const auto& getCoupling(size_t i) const {return coupling_[i];}
+
+    void setDirichletVertices(std::vector<std::shared_ptr<Dune::BitSetVector<1>>> dirichletVertices) {
+        dirichletVertices_ = dirichletVertices;
+    }
+protected:
+    /** \brief Compute the transposed mortar transformation matrix. */
+    void computeTransformationMatrix();
+
+    /** \brief Setup the Householder reflections. */
+    void assembleHouseholderReflections();
+
+public:
+    std::vector<std::shared_ptr<Dune::BitSetVector<1>>> dirichletVertices_;
+
+    /** \brief Bitvector for all bodies whose flags are set if a dof has an obstacle.*/
+    Dune::BitSetVector<dim> totalHasObstacle_;
+
+    /** \brief Obstacles for all bodies on the leafview.*/
+    std::vector<BoxConstraint<field_type,dim> > totalObstacles_;
+
+    std::vector<std::set<size_t>> crosspoints_;
+
+    /** \brief The mortar couplings between grid pairs */
+    std::vector<std::shared_ptr<CouplingType> > contactCoupling_;
+
+    /** \brief The coupling pairs. */
+    std::vector<Dune::Contact::CouplingPair<GridType,GridType,field_type> > coupling_;
+
+    /** \brief Vector containing the involved grids. */
+    std::vector<const GridType*> grids_;
+
+    /** \brief The transposed of the mortar transformation matrix. */
+    MatrixType BT_;
+
+    /** \brief Dismiss all faces that are not at least covered by the grid-glue projection for this
+     *         much percentage ranging between one - for total coverage and zero for taking all faces.
+     */
+    field_type coveredArea_;
+};
+
+#include "nbodyassembler.cc"
+
+#endif

dune/tectonic/spatial-solving/contact/nbodycontacttransfer.cc

+#include <dune/common/fvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/istl/bdmatrix.hh>
+#include <dune/istl/matrixindexset.hh>
+
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
+#include <dune/matrix-vector/blockmatrixview.hh>
+
+#include <dune/matrix-vector/addtodiagonal.hh>
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/solvers/common/numproc.hh>
+#include <dune/solvers/transferoperators/densemultigridtransfer.hh>
+
+template<class ContactNetwork, class VectorType>
+void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwork& contactNetwork, const size_t coarseLevel, const size_t fineLevel) {
+    const size_t nBodies     = contactNetwork.nBodies();
+    const size_t nCouplings = contactNetwork.nCouplings();
+
+    const auto& coarseContactNetwork = *contactNetwork.level(coarseLevel);
+    const auto& fineContactNetwork = *contactNetwork.level(fineLevel);
+
+    std::vector<size_t> maxL(nBodies), cL(nBodies), fL(nBodies);
+
+    for (size_t i=0; i<nBodies; i++) {
+        maxL[i] = contactNetwork.body(i)->grid()->maxLevel();
+        cL[i] = coarseContactNetwork.body(i)->level();
+        fL[i] = fineContactNetwork.body(i)->level();
+    }
+
+    // ////////////////////////////////////////////////////////////
+    //   Create the standard prolongation matrices for each grid
+    // ////////////////////////////////////////////////////////////
+
+    std::vector<TruncatedDenseMGTransfer<VectorType>* > gridTransfer(nBodies);
+    std::vector<const MatrixType*> submat(nBodies);
+
+    for (size_t i=0; i<nBodies; i++) {
+        if (fL[i] > cL[i]) {
+            gridTransfer[i] = new TruncatedDenseMGTransfer<VectorType>;
+            gridTransfer[i]->setup(*contactNetwork.body(i)->grid(), cL[i], fL[i]);
+            submat[i] = &gridTransfer[i]->getMatrix();
+
+        } else {
+            // gridTransfer is identity if coarse and fine level coincide for body
+            Dune::BDMatrix<MatrixBlock>* newMatrix = new Dune::BDMatrix<MatrixBlock>(coarseContactNetwork.body(i)->nVertices());
+
+            for (size_t j=0; j<newMatrix->N(); j++)
+                for (int k=0; k<blocksize; k++)
+                    for (int l=0; l<blocksize; l++)
+                        (*newMatrix)[j][j][k][l] = (k==l);
+
+            submat[i] = newMatrix;
+        }
+    }
+
+    // ///////////////////////////////////
+    //   Combine the submatrices
+    // ///////////////////////////////////
+
+    /*
+            (const std::vector<const GridType*>& grids, int colevel,
+                                                             const std::vector<const MatrixType*>& mortarTransferOperator,
+                                                             const CoordSystemVector& fineLocalCoordSystems,
+                                                             const std::vector<const BitSetVector<1>*>& fineHasObstacle,
+                                                             const std::vector<std::array<int,2> >& gridIdx)
+
+    transfer->setup(*grids[0], *grids[1], i,
+                            contactAssembler.contactCoupling_[0]->mortarLagrangeMatrix(),
+                            contactAssembler.localCoordSystems_,
+                            *);
+
+    */
+
+    if (fineLevel == contactNetwork.nLevels()-1) {
+        const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+        const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+
+        std::vector<const MatrixType*> mortarTransferOperators(nCouplings);
+        std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nCouplings);
+        std::vector<std::array<int,2> > gridIdx(nCouplings);
+
+        for (size_t i=0; i<nCouplings; i++) {
+            mortarTransferOperators[i] = &contactCouplings[i]->mortarLagrangeMatrix();
+            fineHasObstacle[i] = contactCouplings[i]->hasObstacle();
+            gridIdx[i] = nBodyAssembler.getCoupling(i).gridIdx_;
+        }
+
+        combineSubMatrices(submat, mortarTransferOperators, nBodyAssembler.getLocalCoordSystems(), fineHasObstacle, gridIdx);
+    } else {
+        Dune::MatrixVector::BlockMatrixView<MatrixType>::setupBlockMatrix(submat, this->matrix_);
+    }
+
+    for (size_t i=0; i<nBodies; i++) {
+        if (fL[i] <= cL[i]) {
+            delete(submat[i]);
+        }
+        delete(gridTransfer[i]);
+    }
+}
+
+/*
+template<class VectorType>
+template<class GridType>
+void NBodyContactTransfer<VectorType>::setupHierarchy(std::vector<std::shared_ptr<NonSmoothNewtonContactTransfer<VectorType> > >& mgTransfers,
+                                    const std::vector<const GridType*> grids,
+                                    const std::vector<const MatrixType*> mortarTransferOperator,
+                                    const CoordSystemVector& fineLocalCoordSystems,
+                                    const std::vector<const BitSetVector<1>*>& fineHasObstacle,
+                                    const std::vector<std::array<int,2> >& gridIdx)
+{
+    const size_t nGrids     = grids.size();
+
+    std::vector<std::vector<TruncatedDenseMGTransfer<VectorType>* > > gridTransfer(nGrids);
+    std::vector<const MatrixType*> submat(nGrids);
+
+    // ////////////////////////////////////////////////////////////
+    //   Create the standard prolongation matrices for each grid
+    // ////////////////////////////////////////////////////////////
+
+    for (size_t i=0; i<nGrids; i++) {
+
+        gridTransfer[i].resize(grids[i]->maxLevel());
+        for (size_t j=0; j<gridTransfer[i].size(); j++)
+            gridTransfer[i][j] = new TruncatedDenseMGTransfer<VectorType>;
+
+        // Assemble standard transfer operator
+        TransferOperatorAssembler<GridType> transferOperatorAssembler(*grids[i]);
+        transferOperatorAssembler.assembleOperatorPointerHierarchy(gridTransfer[i]);
+    }
+
+    // ////////////////////////////////////////////////////////////////////////////
+    //      Combine matrices in one matrix and add mortar entries on the fine level
+    // ///////////////////////////////////////////////////////////////////////////
+    int toplevel = mgTransfers.size();
+
+    for (size_t colevel=0; colevel<mgTransfers.size(); colevel++) {
+
+        std::vector<int> fL(nGrids);
+
+        for (size_t i=0; i<nGrids; i++) {
+            fL[i] = std::max(size_t(0), grids[i]->maxLevel() - colevel);
+
+            // If the prolongation matrix exists, take it
+            if (fL[i] > 0) {
+
+                submat[i] =&gridTransfer[i][fL[i]-1]->getMatrix();
+            } else {
+                // when the maxLevels of the grids differ then we add copys of the coarsest level to the "smaller" grid
+                BDMatrix<MatrixBlock>* newMatrix = new BDMatrix<MatrixBlock>(grids[i]->size(0,GridType::dimension));
+                *newMatrix = 0;
+
+                for (size_t j=0; j<newMatrix->N(); j++)
+                    for (int k=0; k<blocksize; k++)
+                            (*newMatrix)[j][j][k][k] = 1.0;
+
+                submat[i] = newMatrix;
+            }
+        }
+
+        if (colevel == 0)
+            mgTransfers[toplevel-colevel-1]->combineSubMatrices(submat, mortarTransferOperator, fineLocalCoordSystems, fineHasObstacle, gridIdx);
+        else
+          Dune::MatrixVector::BlockMatrixView<MatrixType>::setupBlockMatrix(submat, mgTransfers[toplevel-colevel-1]->matrix_);
+
+        for (size_t i=0; i<nGrids; i++)
+            if (fL[i]==0)
+                delete(submat[i]);
+    }
+
+    // Delete separate transfer objects
+    for (size_t i=0; i<nGrids; i++)
+        for (size_t j=0; j<gridTransfer[i].size(); j++)
+            delete(gridTransfer[i][j]);
+}
+*/
+
+template<class ContactNetwork, class VectorType>
+void NBodyContactTransfer<ContactNetwork, VectorType>::combineSubMatrices(const std::vector<const MatrixType*>& submat,
+                                                                      const std::vector<const MatrixType*>& mortarTransferOperator,
+                                                                      const CoordSystemVector& fineLocalCoordSystems,
+                                                                      const std::vector<const Dune::BitSetVector<1>*>& fineHasObstacle,
+                                                                      const std::vector<std::array<int,2> >& gridIdx)
+{
+    // ///////////////////////////////////
+    //   Combine the submatrices
+    // ///////////////////////////////////
+
+    const size_t nGrids     = submat.size();
+    const size_t nCouplings = mortarTransferOperator.size();
+
+    Dune::MatrixVector::BlockMatrixView<MatrixType> view(submat);
+
+    Dune::MatrixIndexSet totalIndexSet(view.nRows(), view.nCols());
+
+    // import indices of canonical transfer operator
+    for (size_t i=0; i<nGrids; i++)
+        totalIndexSet.import(*submat[i], view.row(i, 0), view.col(i, 0));
+
+    // ///////////////////////////////////////////////////////////////
+    //   Add additional matrix entries  $ -D^{-1}M I_{mm} $
+    // ///////////////////////////////////////////////////////////////
+
+    typedef typename OperatorType::row_type RowType;
+    typedef typename RowType::ConstIterator ConstIterator;
+
+    std::vector<std::vector<int> > localToGlobal(nCouplings);
+
+    for (size_t i=0; i<nCouplings; i++) {
+
+      if (fineHasObstacle[i]->size() != submat[gridIdx[i][0]]->N())
+        DUNE_THROW(Dune::Exception,
+                   "fineHasObstacle[" << i << "] doesn't have the proper length!");
+
+      localToGlobal[i].resize(mortarTransferOperator[i]->N());
+      size_t idx = 0;
+      for (size_t j=0; j<fineHasObstacle[i]->size(); j++)
+        if ((*fineHasObstacle[i])[j][0])
+          localToGlobal[i][idx++] = j;
+
+      assert(idx==localToGlobal[i].size());
+
+      for (size_t j=0; j<mortarTransferOperator[i]->N(); j++) {
+
+        ConstIterator cTIt = (*mortarTransferOperator[i])[j].begin();
+        ConstIterator cEndIt = (*mortarTransferOperator[i])[j].end();
+        for (; cTIt != cEndIt; ++cTIt) {
+
+          int k = cTIt.index();
+
+          ConstIterator cMIt = (*submat[gridIdx[i][1]])[k].begin();
+          ConstIterator cMEndIt = (*submat[gridIdx[i][1]])[k].end();
+          for (; cMIt != cMEndIt; ++cMIt)
+            totalIndexSet.add(view.row(gridIdx[i][0], localToGlobal[i][j]),
+                view.col(gridIdx[i][1], cMIt.index()));
+
+        }
+
+      }
+
+    }
+
+    totalIndexSet.exportIdx(this->matrix_);
+    this->matrix_ = 0;
+
+    // Copy matrices
+    for (size_t i=0; i<nGrids; i++) {
+
+        for(size_t rowIdx=0; rowIdx<submat[i]->N(); rowIdx++) {
+
+            const RowType& row = (*submat[i])[rowIdx];
+
+            ConstIterator cIt    = row.begin();
+            ConstIterator cEndIt = row.end();
+
+            for(; cIt!=cEndIt; ++cIt)
+                this->matrix_[view.row(i, rowIdx)][view.col(i, cIt.index())] = *cIt;
+
+        }
+
+    }
+
+    // ///////////////////////////////////////////////////////////////
+    //   Add additional matrix entries  $ -D^{-1}M I_{mm} $
+    // ///////////////////////////////////////////////////////////////
+
+    for (size_t i=0; i<nCouplings; i++) {
+
+        for (size_t j=0; j<mortarTransferOperator[i]->N(); j++) {
+
+            ConstIterator cTIt = (*mortarTransferOperator[i])[j].begin();
+            ConstIterator cTEndIt = (*mortarTransferOperator[i])[j].end();
+
+            for (; cTIt != cTEndIt; ++cTIt) {
+
+                int k = cTIt.index();
+
+                ConstIterator cMIt = (*submat[gridIdx[i][1]])[k].begin();
+                ConstIterator cMEndIt = (*submat[gridIdx[i][1]])[k].end();
+
+                for (; cMIt != cMEndIt; ++cMIt) {
+
+                    auto& currentMatrixBlock = this->matrix_[view.row(gridIdx[i][0], localToGlobal[i][j])][view.col(gridIdx[i][1], cMIt.index())];
+
+                    // the entry in the prolongation matrix of the mortar grid
+                    const auto& subMatBlock = this->matrix_[view.row(gridIdx[i][1], k)][view.col(gridIdx[i][1], cMIt.index())];
+
+                    // the - is due to the negation formula for BT
+                    Dune::MatrixVector::subtractProduct(currentMatrixBlock,(*cTIt), subMatBlock);
+                }
+            }
+        }
+    }
+
+    // ///////////////////////////////////////////////////////////////
+    // Transform matrix to account for the local coordinate systems
+    // ///////////////////////////////////////////////////////////////
+
+    /** \todo Hack.  Those should really be equal */
+    assert(fineLocalCoordSystems.size() <= this->matrix_.N());
+
+    for(size_t rowIdx=0; rowIdx<fineLocalCoordSystems.size(); rowIdx++) {
+
+        auto& row = this->matrix_[rowIdx];
+
+        for(auto& col : row)
+            col.leftmultiply(fineLocalCoordSystems[rowIdx]);
+
+    }
+
+}

dune/tectonic/spatial-solving/contact/nbodycontacttransfer.hh

+#ifndef N_BODY_CONTACT_TRANSFER_HH
+#define N_BODY_CONTACT_TRANSFER_HH
+
+#include <vector>
+#include <dune/common/bitsetvector.hh>
+#include <dune/solvers/transferoperators/truncateddensemgtransfer.hh>
+#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
+
+
+/** \brief Multigrid restriction and prolongation operator for contact problems using the TNNMG method.
+ *
+ * Currently only works for first-order Lagrangian elements!
+ *
+ *   This class sets up the multigrid transfer operator for n-body contact
+ *   problems.  It first constructs the standard prolongation matrices for the
+ *   grids and then combines them into a single matrix.  Then, this matrix
+ *   is modified using the mortar transformation basis of Wohlmuth, Krause: "Monotone
+ *   Methods on Nonmatching Grids for Nonlinear Contact Problems". In the TNNMG method
+ *   the nonlinear obstacle problem is only solved on the finest level using a smoother.
+ *   As coarse grid correction a linear multigrid step is applied which is later projected
+ *   back onto the admissible set. Consequently the mortar basis is only needed on the finest
+ *   level. This transfer operator thus additionally to the restriction and truncation,
+ *   transformes the mortar basis back to the nodal basis on the finest level and otherwise
+ *   just does standard restriction/prolongation. On the finest level the prolongation matrix
+ *   is of the form
+
+ *   \f[ \begin{pmatrix}
+ *        (I_l^{l+1})_{ii} & (I_l^{l+1})_{im}       & (I_l^{l+1})_{is} \\
+ *        0                & (I_l^{l+1})_{mm}       & 0 \\
+ *        0                & -M_{l+1}(I_l^{l+1})_{mm} & (I_l^{l+1})_{ss}
+ *       \end{pmatrix} \f].
+ *
+ *   where \f[ I_l^{l+1} \f] denote the standard transfer operators
+ *     and \f[ M_{l+1} \f] denotes the weighted mortar matrix
+ */
+template<class ContactNetwork, class VectorType>
+class NBodyContactTransfer : public TruncatedDenseMGTransfer<VectorType> {
+
+    enum {blocksize = VectorType::block_type::dimension};
+
+    using field_type = typename VectorType::field_type;
+
+    using MatrixBlock = Dune::FieldMatrix<field_type, blocksize, blocksize>;
+    using MatrixType = Dune::BCRSMatrix<MatrixBlock>;
+
+    using CoordSystemVector = std::vector<MatrixBlock>;
+
+public:
+
+    typedef typename MultigridTransfer<VectorType>::OperatorType OperatorType;
+
+    /** \brief Setup transfer operator for a N-body contact problem.
+     *
+     *  \param grids The involved grids
+     *  \param colevel The co-level of the fine grid level to which the transfer operator belongs
+     *  \param mortarTransferOperator The weighted mortar matrices of the couplings
+     *  \param fineLocalCoordSystems A vector containing the local coordinate systems (householder transformations) of all grids
+     *  \param fineHasObstacle  Bitfields determining for each coupling which fine grid nodes belong to the nonmortar boundary.
+     *  \param gridIdx  For each coupling store the indices of the nonmortar and mortar grid.
+    */
+    void setup(const ContactNetwork& contactNetwork, const size_t coarseLevel, const size_t fineLevel);
+
+
+protected:
+    /** \brief Combine all the standard prolongation and mortar operators to one big matrix.
+     *
+     *  \param submat The standard transfer operators
+     *  \param mortarTransferOperator The weighted mortar matrices of the couplings
+     *  \param fineLocalCoordSystems A vector containing the local coordinate systems (householder transformations) of all grids
+     *  \param fineHasObstacle  Bitfields determining for each coupling which fine grid nodes belong to the nonmortar boundary.
+     *  \param gridIdx  For each coupling store the indices of the nonmortar and mortar grid.
+
+    */
+    void combineSubMatrices(const std::vector<const MatrixType*>& submat,
+                            const std::vector<const MatrixType*>& mortarTransferOperator,
+                            const CoordSystemVector& fineLocalCoordSystems,
+                            const std::vector<const Dune::BitSetVector<1>*>& fineHasObstacle,
+                            const std::vector<std::array<int,2> >& gridIdx);
+
+};
+
+#include "nbodycontacttransfer.cc"
+
+#endif

dune/tectonic/spatial-solving/fixedpointiterator.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/solvers/norms/sumnorm.hh>
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+
+#include "functionalfactory.hh"
+#include "nonlinearsolver.hh"
+
+#include "../utils/debugutils.hh"
+
+#include "fixedpointiterator.hh"
+
+void FixedPointIterationCounter::operator+=(
+    FixedPointIterationCounter const &other) {
+  iterations += other.iterations;
+  multigridIterations += other.multigridIterations;
+}
+
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIterator(
+    Dune::ParameterTree const &parset,
+    const NBodyAssembler& nBodyAssembler,
+    const IgnoreVector& ignoreNodes,
+    GlobalFriction& globalFriction,
+    const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+    const ErrorNorms& errorNorms)
+    : parset_(parset),
+      nBodyAssembler_(nBodyAssembler),
+      ignoreNodes_(ignoreNodes),
+      globalFriction_(globalFriction),
+      bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
+      fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
+      fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
+      lambda_(parset.get<double>("v.fpi.lambda")),
+      solverParset_(parset.sub("v.solver")),
+      errorNorms_(errorNorms) {}
+
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+bool FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::displacementCriterion(const Updaters& updaters, std::vector<Vector>& last_u) const {
+    bool criterion = true;
+
+    std::vector<Vector> u;
+    updaters.rate_->extractDisplacement(u);
+
+    const auto& mats = updaters.rate_->getMatrices();
+
+    for (int i=0; i<u.size(); i++) {
+        const EnergyNorm<Matrix, Vector> ANorm(*mats.elasticity[i]), MNorm(*mats.mass[i]);
+        const SumNorm<Vector> AMNorm(1.0, ANorm, 1.0, MNorm);
+
+        if (u[i].size()==0 || last_u[i].size()==0)
+            continue;
+
+        if (AMNorm.diff(u[i], last_u[i]) >= fixedPointTolerance_) {
+            criterion = false;
+            break;
+        }
+    }
+
+    last_u = u;
+
+    return criterion;
+}
+
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+bool FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::stateCriterion(const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& newAlpha) const {
+    bool criterion = true;
+    for (int i=0; i<alpha.size(); i++) {
+        if (alpha[i].size()==0 || newAlpha[i].size()==0)
+            continue;
+
+        if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
+            criterion = false;
+            break;
+        }
+    }
+
+    return criterion;
+}
+
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class LinearSolver>
+FixedPointIterationCounter
+FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
+    Updaters updaters, std::shared_ptr<LinearSolver>& linearSolver,
+    const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+    std::vector<Vector>& velocityIterates) {
+
+    // debugging
+    /*const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
+    for (size_t i=0; i<contactCouplings.size(); i++) {
+      print(*contactCouplings[i]->nonmortarBoundary().getVertices(), "nonmortarBoundaries:");
+    }*/
+
+  const auto nBodies = nBodyAssembler_.nGrids();
+
+  FunctionalFactory<std::decay_t<decltype(nBodyAssembler_)>, decltype(globalFriction_), Matrix, Vector> functionalFactory(nBodyAssembler_, globalFriction_, velocityMatrices, velocityRHSs);
+  functionalFactory.build();
+  auto functional = functionalFactory.functional();
+
+  NonlinearSolver<std::remove_reference_t<decltype(*functional)>, IgnoreVector> solver(parset_.sub("solver.tnnmg"), linearSolver, functional, ignoreNodes_);
+
+  size_t fixedPointIteration;
+  size_t multigridIterations = 0;
+
+  std::vector<Vector> last_u;
+  updaters.rate_->extractOldDisplacement(last_u);
+
+  std::vector<ScalarVector> alpha(nBodies);
+  updaters.state_->extractAlpha(alpha);
+
+  std::vector<Vector> v_o;
+  updaters.rate_->extractOldVelocity(v_o);
+
+  Vector total_v_o;
+  nBodyAssembler_.nodalToTransformed(v_o, total_v_o);
+
+  Vector total_v;
+  nBodyAssembler_.nodalToTransformed(velocityIterates, total_v);
+
+  print(velocityIterates, "velocityIterates start:");
+
+  for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
+       ++fixedPointIteration) {
+
+    // update friction
+    globalFriction_.updateAlpha(alpha);
+
+    // solve rate problem
+    auto res = solver.solve(solverParset_, total_v);
+    multigridIterations += res.iterations;
+
+    Vector v_m = total_v_o;
+    v_m *= 1.0 - lambda_;
+    Dune::MatrixVector::addProduct(v_m, lambda_, total_v);
+
+    // extract relative velocities in mortar basis
+    std::vector<Vector> v_rel;
+    split(v_m, v_rel);
+
+    // solve state problem
+    updaters.state_->solve(v_rel);
+    std::vector<ScalarVector> newAlpha(nBodies);
+    updaters.state_->extractAlpha(newAlpha);
+
+    nBodyAssembler_.postprocess(total_v, velocityIterates);
+    //Rprint(velocityIterates, "velocityIterates loop:");
+    updaters.rate_->postProcess(velocityIterates);
+    bool breakCriterion = true; //displacementCriterion(updaters, last_u); //stateCriterion(alpha, newAlpha);
+
+    printRegularityTruncation(globalFriction_, total_v);
+
+    if (lambda_ < 1e-12 or breakCriterion) {
+      fixedPointIteration++;
+      break;
+    }
+
+    if (res.iterations>200) {
+        std::cout << "FixedPointIterator: TNNMG did not converge, iterations: " << res.iterations << std::endl;
+        //DUNE_THROW(Dune::Exception, "FixedPointIterator: TNNMG did not converge");
+        break;
+    }
+
+    alpha = newAlpha;
+  }
+
+  if (fixedPointIteration == fixedPointMaxIterations_)
+    DUNE_THROW(Dune::Exception, "FPI failed to converge");
+
+  nBodyAssembler_.postprocess(total_v, velocityIterates);
+  updaters.rate_->postProcess(velocityIterates);
+
+  print(velocityIterates, "velocityIterates end:");
+
+  // Cannot use return { fixedPointIteration, multigridIterations };
+  // with gcc 4.9.2, see also http://stackoverflow.com/a/37777814/179927
+  FixedPointIterationCounter ret;
+  ret.iterations = fixedPointIteration;
+  ret.multigridIterations = multigridIterations;
+  return ret;
+}
+
+/*std::ostream &operator<<(std::ostream &stream,
+                         FixedPointIterationCounter const &fpic) {
+  return stream << "(" << fpic.iterations << "," << fpic.multigridIterations
+                << ")";
+}*/
+
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::split(
+    const Vector& v,
+    std::vector<Vector>& splitV) const {
+
+    const size_t nBodies = nBodyAssembler_.nGrids();
+
+    size_t globalIdx = 0;
+    splitV.resize(nBodies);
+    for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+        auto& splitV_ = splitV[bodyIdx];
+
+        splitV_.resize(nBodyAssembler_.grids_[bodyIdx]->size(dims));
+
+        for (size_t i=0; i<splitV_.size(); i++) {
+            splitV_[i] = v[globalIdx];
+            globalIdx++;
+        }
+    }
+
+
+   /*
+        boundaryNodes =
+
+        const auto gridView = contactCouplings[couplingIndices[0]]->nonmortarBoundary().gridView();
+
+        Dune::MultipleCodimMultipleGeomTypeMapper<
+            decltype(gridView), Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
+
+        for (auto it = gridView.template begin<block_size>(); it != gridView.template end<block_size>(); ++it) {
+            const auto i = vertexMapper.index(*it);
+
+            for (size_t j=0; j<couplingIndices.size(); j++) {
+                const auto couplingIdx = couplingIndices[j];
+
+                if (not contactCouplings[couplingIdx]->nonmortarBoundary().containsVertex(i))
+                  continue;
+
+                localToGlobal_.emplace_back(i);
+                restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
+                                          couplings[i]->frictionData()(geoToPoint(it->geometry())));
+                break;
+            }
+
+            globalIdx++;
+        }
+        maxIndex_[bodyIdx] = globalIdx;
+    }*/
+}
+
+#include "fixedpointiterator_tmpl.cc"

src/spatial-solving/fixedpointiterator.hh → dune/tectonic/spatial-solving/fixedpointiterator.hh

@@ -4,10 +4,17 @@
 #include <memory>
 
 #include <dune/common/parametertree.hh>
+#include <dune/common/bitsetvector.hh>
 
 #include <dune/solvers/norms/norm.hh>
 #include <dune/solvers/solvers/solver.hh>
 
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
+#include "solverfactory.hh"
+
 struct FixedPointIterationCounter {
   size_t iterations = 0;
   size_t multigridIterations = 0;
@@ -18,36 +25,59 @@ struct FixedPointIterationCounter {
 std::ostream &operator<<(std::ostream &stream,
                          FixedPointIterationCounter const &fpic);
 
-template <class Factory, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
 class FixedPointIterator {
   using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
   using Vector = typename Factory::Vector;
   using Matrix = typename Factory::Matrix;
-  using ConvexProblem = typename Factory::ConvexProblem;
-  using BlockProblem = typename Factory::BlockProblem;
-  using Nonlinearity = typename ConvexProblem::NonlinearityType;
+
+  using Nonlinearity = typename Factory::Functional::Nonlinearity;
+
+  const static int dims = Vector::block_type::dimension;
+
+  using IgnoreVector = typename Factory::BitVector;
+ // using Nonlinearity = typename ConvexProblem::NonlinearityType;
+ //  using DeformedGrid = typename Factory::DeformedGrid;
+
+public:
+  using GlobalFriction = Nonlinearity;
+  using BitVector = Dune::BitSetVector<1>;
+
+private:
+  void split(const Vector& v, std::vector<Vector>& v_rel) const;
+  bool displacementCriterion(const Updaters& updaters, std::vector<Vector>& last_u) const;
+  bool stateCriterion(const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& newAlpha) const;
 
 public:
-  FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
-                     std::shared_ptr<Nonlinearity> globalFriction,
-                     ErrorNorm const &errorNorm_);
+  FixedPointIterator(const Dune::ParameterTree& parset,
+                     const NBodyAssembler& nBodyAssembler,
+                     const IgnoreVector& ignoreNodes,
+                     GlobalFriction& globalFriction,
+                     const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+                     const ErrorNorms& errorNorms);
 
+  template <class LinearSolver>
   FixedPointIterationCounter run(Updaters updaters,
-                                 Matrix const &velocityMatrix,
-                                 Vector const &velocityRHS,
-                                 Vector &velocityIterate);
+                                 std::shared_ptr<LinearSolver>& linearSolver,
+                                 const std::vector<Matrix>& velocityMatrices,
+                                 const std::vector<Vector>& velocityRHSs,
+                                 std::vector<Vector>& velocityIterates);
 
 private:
-  std::shared_ptr<typename Factory::Step> step_;
-  Dune::ParameterTree const &parset_;
-  std::shared_ptr<Nonlinearity> globalFriction_;
+  const Dune::ParameterTree& parset_;
+  const NBodyAssembler& nBodyAssembler_;
+  const IgnoreVector& ignoreNodes_;
+
+  GlobalFriction& globalFriction_;
+  const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
+
 
   size_t fixedPointMaxIterations_;
   double fixedPointTolerance_;
   double lambda_;
-  size_t velocityMaxIterations_;
-  double velocityTolerance_;
-  Solver::VerbosityMode verbosity_;
-  ErrorNorm const &errorNorm_;
+  const Dune::ParameterTree& solverParset_;
+
+  const ErrorNorms& errorNorms_;
 };
+
 #endif

dune/tectonic/spatial-solving/fixedpointiterator_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "../spatial-solving/solverfactory.hh"
+#include "../spatial-solving/tnnmg/zerononlinearity.hh"
+#include "../data-structures/network/contactnetwork.hh"
+#include "../data-structures/friction/globalfriction.hh"
+
+#include "tnnmg/functional.hh"
+
+#include "../time-stepping/rate/rateupdater.hh"
+#include "../time-stepping/state/stateupdater.hh"
+#include "../time-stepping/updaters.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+
+using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
+
+using LinearSolver = Dune::Solvers::LoopSolver<Vector>;
+using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
+using MyNBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+
+using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
+using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
+using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
+
+template class FixedPointIterator<MySolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>;
+template FixedPointIterationCounter FixedPointIterator<MySolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>::run<LinearSolver>(
+    MyUpdaters, std::shared_ptr<LinearSolver>&, const std::vector<Matrix>&, const std::vector<Vector>&, std::vector<Vector>&);
+
+using NoFriction = ZeroNonlinearity;
+using NoFrictionFunctional = Functional<Matrix&, Vector&, NoFriction&, Vector&, Vector&, double>;
+using NoFrictionSolverFactory = SolverFactory<NoFrictionFunctional, BitVector>;
+
+template class FixedPointIterator<NoFrictionSolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>;
+template FixedPointIterationCounter FixedPointIterator<NoFrictionSolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>::run<LinearSolver>(
+    MyUpdaters, std::shared_ptr<LinearSolver>&, const std::vector<Matrix>&, const std::vector<Vector>&, std::vector<Vector>&);

dune/tectonic/spatial-solving/functionalfactory.hh

+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_GLOBALPROBLEM_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_GLOBALPROBLEM_HH
+
+#include <vector>
+#include <memory>
+
+#include "tnnmg/functional.hh"
+
+template <class NBodyAssembler, class Nonlinearity, class MatrixType, class VectorType>
+class FunctionalFactory {
+    using Vector = VectorType;
+    using Matrix = MatrixType;
+
+    const static int dims = Vector::block_type::dimension;
+
+public:
+    using Functional = Functional<Matrix&, Vector&, Nonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
+
+    FunctionalFactory(const NBodyAssembler& nBodyAssembler, const Nonlinearity& nonlinearity, const std::vector<Matrix>& matrices, const std::vector<Vector>& rhs) :
+        nBodyAssembler_(nBodyAssembler),
+        nonlinearity_(nonlinearity),
+        rhs_(rhs) {
+
+        const auto nBodies = nBodyAssembler_.nGrids();
+
+        matrices_.resize(nBodies);
+        for (int i=0; i<nBodies; i++) {
+            matrices_[i] = &matrices[i];
+        }
+    }
+
+    FunctionalFactory(const NBodyAssembler& nBodyAssembler, const Nonlinearity& nonlinearity, const std::vector<std::shared_ptr<Matrix>>& matrices, const std::vector<Vector>& rhs) :
+        nBodyAssembler_(nBodyAssembler),
+        nonlinearity_(nonlinearity),
+        rhs_(rhs) {
+
+        const auto nBodies = nBodyAssembler_.nGrids();
+
+        matrices_.resize(nBodies);
+        for (size_t i=0; i<nBodies; i++) {
+            matrices_[i] = matrices[i].get();
+      }
+    }
+
+    void build() {
+        // assemble full global contact problem
+        nBodyAssembler_.assembleJacobian(matrices_, bilinearForm_);
+        nBodyAssembler_.assembleRightHandSide(rhs_, totalRhs_);
+
+        // get lower and upper obstacles
+        const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
+        lower_.resize(totalObstacles.size());
+        upper_.resize(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];
+          }
+        }
+
+        // set up functional
+        functional_ = std::make_shared<Functional>(bilinearForm_, totalRhs_, nonlinearity_, lower_, upper_);
+    }
+
+    std::shared_ptr<Functional> functional() const {
+        return functional_;
+    }
+
+    const Matrix& bilinearForm() const {
+        return bilinearForm_;
+    }
+
+    const Vector& rhs() const {
+        return totalRhs_;
+    }
+
+    const Vector& lower() const {
+        return lower_;
+    }
+
+    const Vector& upper() const {
+        return upper_;
+    }
+
+private:
+    const NBodyAssembler nBodyAssembler_;
+    const Nonlinearity& nonlinearity_;
+    std::vector<const Matrix*> matrices_;
+    const std::vector<Vector>& rhs_;
+
+    Matrix bilinearForm_;
+    Vector totalRhs_;
+    Vector lower_;
+    Vector upper_;
+
+    std::shared_ptr<Functional> functional_;
+};
+#endif

dune/tectonic/spatial-solving/makelinearsolver.hh

+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH
+
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "preconditioners/multilevelpatchpreconditioner.hh"
+
+template <class ContactNetwork, class VectorType>
+auto makeLinearSolver(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) {
+    // make linear solver for linear correction in TNNMGStep
+    using Vector = VectorType;
+    using Matrix = typename ContactNetwork::Matrix;
+
+        /* old, pre multi threading, was unused
+        using Norm =  EnergyNorm<Matrix, Vector>;
+        using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+        const auto& preconditionerParset = parset_.sub("solver.tnnmg.linear.preconditioner");
+
+        Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
+        Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
+        preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
+        preconditioner.build();
+
+        auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
+        cgStep->setPreconditioner(preconditioner);
+
+        Norm norm(*cgStep);
+
+        return std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.linear.tolerance"), norm, Solver::QUIET);
+        */
+
+        // patch preconditioner only needs to be computed once per advance()
+        // make linear solver for linear correction in TNNMGStep
+        using Norm = EnergyNorm<Matrix, Vector>;
+        using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+        /*const auto& preconditionerParset = parset_.sub("solver.tnnmg.preconditioner");
+
+        Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
+        Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
+        preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
+        preconditioner.build();
+
+        auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
+        cgStep->setPreconditioner(preconditioner);
+
+        Norm norm(*cgStep);
+
+        auto linearSolver = std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+        */
+
+
+
+
+        // transfer operators need to be recomputed on change due to setDeformation()
+        using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+        using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+        TransferOperators transfer(contactNetwork.nLevels()-1);
+        for (size_t i=0; i<transfer.size(); i++) {
+            transfer[i] = std::make_shared<TransferOperator>();
+            transfer[i]->setup(contactNetwork, i, i+1);
+        }
+
+        // Remove any recompute filed so that initially the full transferoperator is assembled
+        for (size_t i=0; i<transfer.size(); i++)
+            std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+        auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+        linearMultigridStep->setMGType(1, 3, 3);
+        linearMultigridStep->setSmoother(TruncatedBlockGSStep<Matrix, Vector>());
+        linearMultigridStep->setTransferOperators(transfer);
+
+        return std::make_shared<LinearSolver>(linearMultigridStep, parset.template get<int>("solver.tnnmg.main.multi"), parset.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), Norm(*linearMultigridStep), Solver::QUIET);
+    }
+#endif

dune/tectonic/spatial-solving/nonlinearsolver.hh

+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_SOLVENONLINEARPROBLEM_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_SOLVENONLINEARPROBLEM_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/contact/common/dualbasisadapter.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/functions/gridfunctions/gridfunction.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
+#include "../data-structures/enums.hh"
+#include "../data-structures/enumparser.hh"
+
+
+
+#include "../utils/tobool.hh"
+#include "../utils/debugutils.hh"
+
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include "../spatial-solving/contact/nbodycontacttransfer.hh"
+
+#include "solverfactory.hh"
+#include "solverfactory.cc"
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+
+#include "fixedpointiterator.hh"
+
+
+template <class Functional, class BitVector>
+class NonlinearSolver {
+protected:
+    using Matrix = typename Functional::Matrix;
+    using Vector = typename Functional::Vector;
+    using Factory = SolverFactory<Functional, BitVector>;
+
+    using Norm = EnergyNorm<Matrix, Vector>;
+    using SolverType = Dune::Solvers::LoopSolver<Vector>;
+
+    const static int dims = Vector::block_type::dimension;
+
+public:
+    template <class LinearSolver>
+    NonlinearSolver(const Dune::ParameterTree& tnnmgParset, std::shared_ptr<LinearSolver> linearSolver, std::shared_ptr<Functional> functional, const BitVector& dirichletNodes) :
+        functional_(functional),
+        norm_(functional_->quadraticPart()) {
+
+        // set up TNMMG step
+        solverFactory_ = std::make_shared<Factory>(tnnmgParset, functional_, dirichletNodes);
+        solverFactory_->build(linearSolver);
+    }
+
+    auto solve(const Dune::ParameterTree& solverParset, Vector& x) {
+        auto tnnmgStep = solverFactory_->step();
+
+        SolverType solver(*tnnmgStep.get(), solverParset.get<size_t>("maximumIterations"),
+            solverParset.get<double>("tolerance"), norm_,
+            solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+
+        const auto& lower = functional_->lowerObstacle();
+        const auto& upper = functional_->upperObstacle();
+
+        // project in onto admissible set
+        for (size_t i=0; i<x.size(); i++) {
+            for (size_t j=0; j<dims; j++) {
+                if (x[i][j] < lower[i][j]) {
+                    x[i][j] = lower[i][j];
+                }
+
+                if (x[i][j] > upper[i][j]) {
+                    x[i][j] = upper[i][j];
+                }
+            }
+        }
+
+        solverFactory_->setProblem(x);
+
+        solver.preprocess();
+        solver.solve();
+
+        return solver.getResult();
+    }
+
+private:
+    std::shared_ptr<Functional> functional_;
+    std::shared_ptr<Factory> solverFactory_;
+
+    Norm norm_;
+};
+
+#endif