src/multi-body-problem-data/segmented-function.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include "grid/cuboidgeometry.hh"
+
+class SegmentedFunction
+    : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
+                                   Dune::FieldVector<double, 1>> {
+private:
+  bool liesBelow(Dune::FieldVector<double, MY_DIM> const &x,
+                 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 false; //TODO liesBelow(MyGeometry::K, MyGeometry::M, z);
+  };
+
+  double const _v1;
+  double const _v2;
+
+public:
+  SegmentedFunction(double v1, double v2) : _v1(v1), _v2(v2) {}
+
+  void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
+                Dune::FieldVector<double, 1> &y) const {
+    y = insideRegion2(x) ? _v2 : _v1;
+  }
+};
+#endif

src/multi-body-problem.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/formatstring.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/projections/normalprojection.hh>
+
+#include <dune/tectonic/geocoordinate.hh>
+#include <dune/tectonic/globalfriction.hh>
+
+#include "assemblers.hh"
+#include "gridselector.hh"
+#include "explicitgrid.hh"
+#include "explicitvectors.hh"
+
+#include "data-structures/enumparser.hh"
+#include "data-structures/enums.hh"
+#include "data-structures/contactnetwork.hh"
+#include "data-structures/matrices.hh"
+#include "data-structures/program_state.hh"
+
+#include "factories/stackedblocksfactory.hh"
+#include "factories/threeblocksfactory.hh"
+
+//#include "io/hdf5-levelwriter.hh"
+#include "io/hdf5/restart-io.hh"
+#include "io/vtk.hh"
+
+#include "multi-body-problem-data/bc.hh"
+#include "multi-body-problem-data/mybody.hh"
+#include "multi-body-problem-data/grid/mygrids.hh"
+
+#include "spatial-solving/tnnmg/functional.hh"
+//#include "spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
+#include "spatial-solving/tnnmg/localbisectionsolver.hh"
+#include "spatial-solving/contact/nbodyassembler.hh"
+#include "spatial-solving/solverfactory.hh"
+
+#include "time-stepping/adaptivetimestepper.hh"
+#include "time-stepping/rate.hh"
+#include "time-stepping/state.hh"
+#include "time-stepping/updaters.hh"
+
+#include "utils/debugutils.hh"
+#include "utils/diameter.hh"
+
+// for getcwd
+#include <unistd.h>
+
+//#include <tbb/tbb.h> //TODO multi threading preconditioner?
+//#include <pthread.h>
+
+size_t const dims = MY_DIM;
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+  Dune::ParameterTree parset;
+  Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem.cfg", parset);
+  Dune::ParameterTreeParser::readINITree(
+      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem-%dD.cfg", dims), parset);
+  Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+  return parset;
+}
+
+static std::atomic<bool> terminationRequested(false);
+void handleSignal(int signum) { terminationRequested = true; }
+
+int main(int argc, char *argv[]) {
+  try {
+    Dune::MPIHelper::instance(argc, argv);
+
+    char buffer[256];
+    char *val = getcwd(buffer, sizeof(buffer));
+    if (val) {
+        std::cout << buffer << std::endl;
+        std::cout << argv[0] << std::endl;
+    }
+
+    std::ofstream out("../log.txt");
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
+
+    auto const parset = getParameters(argc, argv);
+
+    using Assembler = MyAssembler<DefLeafGridView, dims>;
+    using field_type = Matrix::field_type;
+
+    // ----------------------
+    // set up contact network
+    // ----------------------
+    StackedBlocksFactory<Grid, Vector> stackedBlocksFactory(parset);
+    using ContactNetwork = typename StackedBlocksFactory<Grid, Vector>::ContactNetwork;
+    stackedBlocksFactory.build();
+
+    ContactNetwork& contactNetwork = stackedBlocksFactory.contactNetwork();
+
+    /*ThreeBlocksFactory<Grid, Vector> threeBlocksFactory(parset);
+    using ContactNetwork = typename ThreeBlocksFactory<Grid, Vector>::ContactNetwork;
+    threeBlocksFactory.build();
+
+    ContactNetwork& contactNetwork = threeBlocksFactory.contactNetwork(); */
+
+    const size_t bodyCount = contactNetwork.nBodies();
+
+    /*
+    for (size_t i=0; i<bodyCount; i++) {
+        // printDofLocation(contactNetwork.body(i)->gridView());
+
+        //Vector def(contactNetwork.deformedGrids()[i]->size(dims));
+        //def = 1;
+        //deformedGridComplex.setDeformation(def, i);
+
+        auto& levelViews = contactNetwork.levelViews(i);
+
+        for (size_t j=0; j<levelViews.size(); j++) {
+            writeToVTK(*levelViews[j], "", "body_" + std::to_string(i) + "_level_" + std::to_string(j));
+        }
+
+        writeToVTK(contactNetwork.leafView(i), "", "body_" + std::to_string(i) + "_leaf");
+    } */
+
+    // ----------------------------
+    // assemble contactNetwork
+    // ----------------------------
+    contactNetwork.assemble();
+
+    //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+    // -----------------
+    // init input/output
+    // -----------------
+    std::vector<size_t> nVertices(bodyCount);
+    for (size_t i=0; i<bodyCount; i++) {
+        nVertices[i] = contactNetwork.body(i)->nVertices();
+    }
+
+    using MyProgramState = ProgramState<Vector, ScalarVector>;
+    MyProgramState programState(nVertices);
+    auto const firstRestart = parset.get<size_t>("io.restarts.first");
+    auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
+    auto const writeRestarts = parset.get<bool>("io.restarts.write");
+    auto const writeData = parset.get<bool>("io.data.write");
+    bool const handleRestarts = writeRestarts or firstRestart > 0;
+
+
+    auto dataFile =
+        writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+
+    auto restartFile = handleRestarts
+                           ? std::make_unique<HDF5::File>(
+                                 "restarts.h5",
+                                 writeRestarts ? HDF5::Access::READWRITE
+                                               : HDF5::Access::READONLY)
+                           : nullptr;
+
+
+    auto restartIO = handleRestarts
+                         ? std::make_unique<RestartIO<MyProgramState>>(
+                               *restartFile, nVertices)
+                         : nullptr;
+
+    if (firstRestart > 0) // automatically adjusts the time and timestep
+      restartIO->read(firstRestart, programState);
+    else
+     programState.setupInitialConditions(parset, contactNetwork);
+
+
+    auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+    for (size_t i=0; i<bodyCount; i++) {
+      contactNetwork.body(i)->setDeformation(programState.u[i]);
+    }
+    nBodyAssembler.assembleTransferOperator();
+    nBodyAssembler.assembleObstacle();
+
+    // ------------------------
+    // assemble global friction
+    // ------------------------
+    contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
+
+    auto& globalFriction = contactNetwork.globalFriction();
+    globalFriction.updateAlpha(programState.alpha);
+
+    using MyVertexBasis = typename Assembler::VertexBasis;
+    using MyCellBasis = typename Assembler::CellBasis;
+    std::vector<Vector> vertexCoordinates(bodyCount);
+    std::vector<const MyVertexBasis* > vertexBases(bodyCount);
+    std::vector<const MyCellBasis* > cellBases(bodyCount);
+
+    for (size_t i=0; i<bodyCount; i++) {
+      const auto& body = contactNetwork.body(i);
+      vertexBases[i] = &(body->assembler()->vertexBasis);
+      cellBases[i] = &(body->assembler()->cellBasis);
+
+      auto& vertexCoords = vertexCoordinates[i];
+      vertexCoords.resize(nVertices[i]);
+
+      Dune::MultipleCodimMultipleGeomTypeMapper<
+          DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(body->gridView(), Dune::mcmgVertexLayout());
+      for (auto &&v : vertices(body->gridView()))
+        vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
+    }
+
+    //typename contactNetwork::BoundaryPatches frictionBoundaries;
+    //contactNetwork.boundaryPatches("friction", frictionBoundaries);
+
+    /*
+    auto dataWriter =
+        writeData ? std::make_unique<
+                        HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
+                        *dataFile, vertexCoordinates, vertexBases,
+                        frictionBoundaries) //, weakPatches)
+                  : nullptr;*/
+
+    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
+
+    IterationRegister iterationCount;
+
+    auto const report = [&](bool initial = false) {
+      /*if (writeData) {
+        dataWriter->reportSolution(programState, contactNetwork.globalFriction());
+        if (!initial)
+          dataWriter->reportIterations(programState, iterationCount);
+        dataFile->flush();
+      }
+
+      if (writeRestarts and !initial and
+          programState.timeStep % restartSpacing == 0) {
+        restartIO->write(programState);
+        restartFile->flush();
+      }*/
+
+      if (parset.get<bool>("io.printProgress"))
+        std::cout << "timeStep = " << std::setw(6) << programState.timeStep
+                  << ", time = " << std::setw(12) << programState.relativeTime
+                  << ", tau = " << std::setw(12) << programState.relativeTau
+                  << std::endl;
+
+      if (parset.get<bool>("io.vtk.write")) {
+        std::vector<ScalarVector> stress(bodyCount);
+
+        for (size_t i=0; i<bodyCount; i++) {
+          const auto& body = contactNetwork.body(i);
+          body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
+                                           body->data()->getPoissonRatio(),
+                                           programState.u[i], stress[i]);
+
+        }
+
+        vtkWriter.write(programState.timeStep, programState.u, programState.v,
+                        programState.alpha, stress);
+      }
+    };
+    report(true);
+
+    // -------------------
+    // Set up TNNMG solver
+    // -------------------
+
+    BitVector totalDirichletNodes;
+    contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+
+    /*for (size_t i=0; i<totalDirichletNodes.size(); i++) {
+        bool val = false;
+        for (size_t d=0; d<dims; d++) {
+            val = val || totalDirichletNodes[i][d];
+        }
+
+        totalDirichletNodes[i] = val;
+        for (size_t d=0; d<dims; d++) {
+            totalDirichletNodes[i][d] = val;
+        }
+    }*/
+
+    //print(totalDirichletNodes, "totalDirichletNodes:");
+
+    using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
+    using NonlinearFactory = SolverFactory<Functional, BitVector>;
+
+    using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
+    using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+    using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
+                               StateUpdater<ScalarVector, Vector>>;
+
+    BoundaryFunctions velocityDirichletFunctions;
+    contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+    BoundaryNodes dirichletNodes;
+    contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+    /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+        for (size_t j=0; j<dirichletNodes[i].size(); j++) {
+        print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
+        }
+    }*/
+
+    std::vector<const Dune::BitSetVector<1>*> frictionNodes;
+    contactNetwork.frictionNodes(frictionNodes);
+
+    /*for (size_t i=0; i<frictionNodes.size(); i++) {
+        print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
+    }*/
+
+    Updaters current(
+        initRateUpdater(
+            parset.get<Config::scheme>("timeSteps.scheme"),
+            velocityDirichletFunctions,
+            dirichletNodes,
+            contactNetwork.matrices(),
+            programState.u,
+            programState.v,
+            programState.a),
+        initStateUpdater<ScalarVector, Vector>(
+            parset.get<Config::stateModel>("boundary.friction.stateModel"),
+            programState.alpha,
+            nBodyAssembler.getContactCouplings(),
+            contactNetwork.couplings())
+            );
+
+
+    auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
+
+    const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
+
+    auto const mustRefine = [&](Updaters &coarseUpdater,
+                                Updaters &fineUpdater) {
+
+        //return false;
+      std::cout << "Step 1" << std::endl;
+
+      std::vector<ScalarVector> coarseAlpha;
+      coarseAlpha.resize(bodyCount);
+      coarseUpdater.state_->extractAlpha(coarseAlpha);
+
+      print(coarseAlpha, "coarseAlpha:");
+
+      std::vector<ScalarVector> fineAlpha;
+      fineAlpha.resize(bodyCount);
+      fineUpdater.state_->extractAlpha(fineAlpha);
+
+      print(fineAlpha, "fineAlpha:");
+
+      std::cout << "Step 3" << std::endl;
+
+      ScalarVector::field_type energyNorm = 0;
+      for (size_t i=0; i<stateEnergyNorms.size(); i++) {
+          std::cout << "for " << i << std::endl;
+
+          std::cout << not stateEnergyNorms[i] << std::endl;
+
+          if (coarseAlpha[i].size()==0 || fineAlpha[i].size()==0)
+              continue;
+
+          energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
+      }
+      std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance <<  std::endl;
+      std::cout << "must refine: " << (energyNorm > refinementTolerance) <<  std::endl;
+      return energyNorm > refinementTolerance;
+    };
+
+
+    std::signal(SIGXCPU, handleSignal);
+    std::signal(SIGINT, handleSignal);
+    std::signal(SIGTERM, handleSignal);
+
+/*
+    // set patch preconditioner for linear correction in TNNMG method
+    using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, PatchSolver, Matrix, Vector>;
+
+    const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
+
+    auto gsStep = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+    PatchSolver patchSolver(gsStep,
+                               preconditionerParset.get<size_t>("maximumIterations"),
+                               preconditionerParset.get<double>("tolerance"),
+                               nullptr,
+                               preconditionerParset.get<Solver::VerbosityMode>("verbosity"),
+                               false); // absolute error
+
+    Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+    Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<Preconditioner::Mode>("mode"));
+    preconditioner.setPatchSolver(patchSolver);
+    preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
+*/
+    // set adaptive time stepper
+    typename ContactNetwork::ExternalForces externalForces;
+    contactNetwork.externalForces(externalForces);
+
+    AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        adaptiveTimeStepper(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes, current,
+                            programState.relativeTime, programState.relativeTau,
+                            externalForces, stateEnergyNorms, mustRefine);
+
+    while (!adaptiveTimeStepper.reachedEnd()) {
+      programState.timeStep++;
+
+      //preconditioner.build();
+      iterationCount = adaptiveTimeStepper.advance();
+
+      programState.relativeTime = adaptiveTimeStepper.relativeTime_;
+      programState.relativeTau = adaptiveTimeStepper.relativeTau_;
+      current.rate_->extractDisplacement(programState.u);
+      current.rate_->extractVelocity(programState.v);
+      current.rate_->extractAcceleration(programState.a);
+      current.state_->extractAlpha(programState.alpha);
+
+      print(programState.u, "current u:");
+      print(programState.v, "current v:");
+      print(programState.a, "current a:");
+      print(programState.alpha, "current alpha:");
+
+      contactNetwork.setDeformation(programState.u);
+
+      report();
+
+      break;
+
+      if (programState.timeStep==50) {
+        std::cout << "limit of timeSteps reached!" << std::endl;
+        break; // TODO remove after debugging
+      }
+
+      if (terminationRequested) {
+        std::cerr << "Terminating prematurely" << std::endl;
+        break;
+      }
+    }
+
+
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+
+  } catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  }
+}

src/multi-body-problem.cfg

+# -*- mode:conf -*-
+gravity         = 9.81  # [m/s^2]
+
+[io]
+data.write      = false #true
+printProgress   = true
+restarts.first  = 0
+restarts.spacing= 20
+restarts.write  = false #true
+vtk.write       = true
+
+[problem]
+finalTime       = 100  # [s] #1000
+bodyCount       = 2
+
+[body]
+bulkModulus     = 0.5e5 # [Pa]
+poissonRatio    = 0.3   # [1]
+[body.elastic]
+density         = 900   # [kg/m^3]
+shearViscosity  = 1e3   # [Pas]
+bulkViscosity   = 1e3   # [Pas]
+[body.viscoelastic]
+density         = 1000  # [kg/m^3]
+shearViscosity  = 1e4   # [Pas]
+bulkViscosity   = 1e4   # [Pas]
+
+[boundary.friction]
+C               = 10    # [Pa]
+mu0             = 0.7   # [ ]
+V0              = 5e-5  # [m/s]
+L               = 2.25e-5 # [m]
+initialAlpha    = 0     # [ ]
+stateModel      = AgeingLaw
+frictionModel   = Regularised
+[boundary.friction.weakening]
+a               = 0.002 # [ ]
+b               = 0.017 # [ ]
+[boundary.friction.strengthening]
+a               = 0.020 # [ ]
+b               = 0.005 # [ ]
+
+[initialTime]
+timeStep = 0
+relativeTime = 0.0
+relativeTau = 1e-4 # 1e-6
+
+[timeSteps]
+scheme = newmark
+
+[u0.solver]
+maximumIterations = 100
+verbosity         = full
+
+[a0.solver]
+maximumIterations = 100
+verbosity         = full
+
+[v.solver]
+maximumIterations = 100
+verbosity         = full
+
+[v.fpi]
+maximumIterations = 10000
+lambda            = 0.5
+
+[solver.tnnmg.linear]
+maximumIterations = 100
+pre                = 3
+cycle              = 1  # 1 = V, 2 = W, etc.
+post               = 3
+[solver.tnnmg.linear.preconditioner]
+mode         = additive
+patchDepth   = 0
+maximumIterations = 10000
+verbosity         = quiet
+
+[solver.tnnmg.main]
+pre   = 1
+multi = 5 # number of multigrid steps
+post  = 0
+
+
+

src/nodalweights.cc

+#include "utils/almostequal.hh"
+
+#include "nodalweights.hh"
+
+template <class Basis0, class Basis1>
+NodalWeights<Basis0, Basis1>::NodalWeights(const Basis0& basis0, const Basis1& basis1)
+    : basis0_(basis0),
+      basis1_(basis1)
+    {}
+
+template <class Basis0, class Basis1>
+template <class Basis, class Element, class GlobalCoords>
+auto NodalWeights<Basis0, Basis1>::basisDof(const Basis& basis, const Element& elem, const GlobalCoords& vertex) const {
+    int dof = -1;
+
+    const typename Basis::LocalFiniteElement& lFE = basis.getLocalFiniteElement(elem);
+    const auto& geometry = elem.geometry();
+    const auto& localVertex = geometry.local(vertex);
+
+    const size_t localBasisSize = lFE.localBasis().size();
+    std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+    lFE.localBasis().evaluateFunction(localVertex, localBasisRep);
+
+    for(size_t i=0; i<localBasisSize; i++) {
+        if (almost_equal(localBasisRep[i][0], 1.0, 2)) {
+            dof = basis.index(elem, i);
+            break;
+        }
+    }
+
+    return dof;
+}
+
+template <class Basis0, class Basis1>
+template <class GridGlue, class ScalarVector>
+void NodalWeights<Basis0, Basis1>::assemble(const GridGlue& glue, ScalarVector& weights0, ScalarVector& weights1, bool initializeVector) const {
+    using ctype = typename ScalarVector::field_type;
+
+    if (initializeVector==true) {
+        weights0.resize(basis0_.size());
+        weights1.resize(basis1_.size());
+    }
+
+    // loop over all intersections
+    for (const auto& rIs : intersections(glue)) {
+        const auto& inside = rIs.inside();
+        const auto& outside = rIs.outside();
+
+        /*if (!nmBoundary.contains(inside, rIs.indexInInside())) {
+            std::cout << "it happened" << std::endl;
+            continue;
+        }*/
+
+        const auto& geometry = rIs.geometry();
+        const ctype val = 1.0/(geometry.mydimension + 1)*geometry.volume();
+
+        std::cout << geometry.mydimension  << " " << geometry.volume() << " " << val << std::endl;
+
+        for (int i=0; i<geometry.corners(); i++) {
+            const auto& vertex = geometry.corner(i);
+
+            const int inIdx = basisDof(basis0_, inside, vertex);
+            if (inIdx >= 0) {
+                weights0[inIdx] += val;
+            }
+
+            const int outIdx = basisDof(basis1_, outside, vertex);
+            if (outIdx >= 0) {
+                weights1[outIdx] += val;
+            }
+
+            std::cout << inIdx << " " << outIdx << std::endl;
+        }
+    }
+}

src/nodalweights.hh

+#ifndef SRC_NODALWEIGHTS_HH
+#define SRC_NODALWEIGHTS_HH
+
+
+/**
+ * let merged contact boundary Gamma be given by GridGlue object;
+ * computes
+ *                  int_Gamma lambda_i dx,
+ * where lambda_i is nodal basis of merged contact boundary whose
+ * dofs are given by nonmortar and mortar dofs;
+ *
+ *
+ * NOTE: works only for P1 nodal bases
+ **/
+
+template <class Basis0, class Basis1>
+class NodalWeights {
+private:
+    enum {dim = Basis0::GridView::Grid::dimension};
+
+    template <class Basis, class Element, class GlobalCoords>
+    auto basisDof(const Basis& basis, const Element& elem, const GlobalCoords& vertex) const;
+
+public:
+    NodalWeights(const Basis0& basis0, const Basis1& basis1);
+
+    template <class GridGlue, class ScalarVector>
+    void assemble(const GridGlue& glue, ScalarVector& weights0, ScalarVector& weights1, bool initializeVector = true) const;
+
+private:
+    const Basis0& basis0_;
+    const Basis1& basis1_;
+};
+
+#include "nodalweights.cc"
+
+#endif

src/one-body-problem-data/mygrid.cc

@@ -6,7 +6,7 @@
 
 #include "mygrid.hh"
 #include "midpoint.hh"
-#include "../diameter.hh"
+#include "../utils/diameter.hh"
 
 #if MY_DIM == 3
 SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}

src/one-body-problem.cc

@@ -31,9 +31,14 @@
 #include <dune/fufem/formatstring.hh>
 
 #include <dune/solvers/norms/energynorm.hh>
+
+
+/*
 #include <dune/solvers/solvers/loopsolver.hh>
 #include <dune/solvers/solvers/solver.hh>
 #include <dune/tnnmg/problem-classes/convexproblem.hh>
+*/
+
 
 #include <dune/tectonic/geocoordinate.hh>
 #include <dune/tectonic/myblockproblem.hh>
@@ -41,34 +46,46 @@
 #include <dune/fufem/hdf5/file.hh>
 
 #include "assemblers.hh"
-#include "diameter.hh"
-#include "enumparser.hh"
-#include "enums.hh"
+#include "boundarycondition.hh"
 #include "gridselector.hh"
-#include "hdf5-writer.hh"
-#include "hdf5/restart-io.hh"
-#include "matrices.hh"
-#include "program_state.hh"
+
+#include "data-structures/enumparser.hh"
+#include "data-structures/enums.hh"
+#include "data-structures/matrices.hh"
+#include "data-structures/program_state.hh"
+
+#include "io/hdf5-writer.hh"
+#include "io/hdf5/restart-io.hh"
+#include "io/vtk.hh"
+
 #include "one-body-problem-data/bc.hh"
 #include "one-body-problem-data/mybody.hh"
 #include "one-body-problem-data/mygeometry.hh"
 #include "one-body-problem-data/myglobalfrictiondata.hh"
 #include "one-body-problem-data/mygrid.hh"
 #include "one-body-problem-data/weakpatch.hh"
+
 #include "spatial-solving/solverfactory.hh"
+
 #include "time-stepping/adaptivetimestepper.hh"
 #include "time-stepping/rate.hh"
 #include "time-stepping/state.hh"
 #include "time-stepping/updaters.hh"
-#include "vtk.hh"
+
+#include "utils/diameter.hh"
+
+// for getcwd
+#include <unistd.h>
+
+#define USE_OLD_TNNMG
 
 size_t const dims = MY_DIM;
 
 Dune::ParameterTree getParameters(int argc, char *argv[]) {
   Dune::ParameterTree parset;
-  Dune::ParameterTreeParser::readINITree("one-body-problem.cfg", parset);
+  Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/one-body-problem.cfg", parset);
   Dune::ParameterTreeParser::readINITree(
-      Dune::Fufem::formatString("one-body-problem-%dD.cfg", dims), parset);
+      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/one-body-problem-%dD.cfg", dims), parset);
   Dune::ParameterTreeParser::readOptions(argc, argv, parset);
   return parset;
 }
@@ -79,6 +96,14 @@ void handleSignal(int signum) { terminationRequested = true; }
 int main(int argc, char *argv[]) {
   try {
     Dune::MPIHelper::instance(argc, argv);
+
+    char buffer[256];
+    char *val = getcwd(buffer, sizeof(buffer));
+    if (val) {
+        std::cout << buffer << std::endl;
+        std::cout << argv[0] << std::endl;
+    }
+
     auto const parset = getParameters(argc, argv);
 
     MyGeometry::render();
@@ -140,6 +165,7 @@ int main(int argc, char *argv[]) {
 #endif
     }
 
+
     // Set up functions for time-dependent boundary conditions
     using Function = Dune::VirtualFunction<double, double>;
     Function const &velocityDirichletFunction = VelocityDirichletCondition();
@@ -184,8 +210,10 @@ int main(int argc, char *argv[]) {
     auto const writeData = parset.get<bool>("io.data.write");
     bool const handleRestarts = writeRestarts or firstRestart > 0;
 
+    /*
     auto dataFile =
         writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+
     auto restartFile = handleRestarts
                            ? std::make_unique<HDF5::File>(
                                  "restarts.h5",
@@ -214,7 +242,7 @@ int main(int argc, char *argv[]) {
     Vector vertexCoordinates(leafVertexCount);
     {
       Dune::MultipleCodimMultipleGeomTypeMapper<
-          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
+          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView, Dune::mcmgVertexLayout());
       for (auto &&v : vertices(leafView))
         vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
     }
@@ -229,6 +257,7 @@ int main(int argc, char *argv[]) {
     MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
         myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
 
+
     IterationRegister iterationCount;
     auto const report = [&](bool initial = false) {
       if (writeData) {
@@ -261,6 +290,7 @@ int main(int argc, char *argv[]) {
     };
     report(true);
 
+
     // Set up TNNMG solver
     using NonlinearFactory = SolverFactory<
         dims,
@@ -321,6 +351,8 @@ int main(int argc, char *argv[]) {
         break;
       }
     }
+
+    */
   } catch (Dune::Exception &e) {
     Dune::derr << "Dune reported error: " << e << std::endl;
   } catch (std::exception &e) {

src/one-body-problem.cfg

@@ -59,7 +59,7 @@ maximumIterations = 10000
 lambda            = 0.5
 
 [solver.tnnmg.linear]
-maxiumumIterations = 100000
+maximumIterations = 100000
 pre                = 3
 cycle              = 1  # 1 = V, 2 = W, etc.
 post               = 3

src/program_state.hh

-#ifndef SRC_PROGRAM_STATE_HH
-#define SRC_PROGRAM_STATE_HH
-
-#include <dune/common/parametertree.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
-
-#include <dune/tectonic/body.hh>
-
-#include "assemblers.hh"
-#include "matrices.hh"
-#include "spatial-solving/solverfactory.hh"
-
-template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
-public:
-  using Vector = VectorTEMPLATE;
-  using ScalarVector = ScalarVectorTEMPLATE;
-
-  ProgramState(int leafVertexCount)
-      : u(leafVertexCount),
-        v(leafVertexCount),
-        a(leafVertexCount),
-        alpha(leafVertexCount),
-        weightedNormalStress(leafVertexCount) {}
-
-  // Set up initial conditions
-  template <class Matrix, class GridView>
-  void setupInitialConditions(
-      Dune::ParameterTree const &parset,
-      std::function<void(double, Vector &)> externalForces,
-      Matrices<Matrix> const matrices,
-      MyAssembler<GridView, Vector::block_type::dimension> const &myAssembler,
-      Dune::BitSetVector<Vector::block_type::dimension> const &dirichletNodes,
-      Dune::BitSetVector<Vector::block_type::dimension> const &noNodes,
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      Body<Vector::block_type::dimension> const &body) {
-
-    using LocalVector = typename Vector::block_type;
-    using LocalMatrix = typename Matrix::block_type;
-    auto constexpr dims = LocalVector::dimension;
-
-    // Solving a linear problem with a multigrid solver
-    auto const solveLinearProblem = [&](
-        Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
-        Vector const &_rhs, Vector &_x,
-        Dune::ParameterTree const &_localParset) {
-
-      using LinearFactory = SolverFactory<
-          dims, BlockNonlinearTNNMGProblem<ConvexProblem<
-                    ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-          typename GridView::Grid>;
-      ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
-
-      LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
-                            myAssembler.gridView.grid(), _dirichletNodes);
-
-      typename LinearFactory::ConvexProblem convexProblem(
-          1.0, _matrix, zeroNonlinearity, _rhs, _x);
-      typename LinearFactory::BlockProblem problem(parset, convexProblem);
-
-      auto multigridStep = factory.getStep();
-      multigridStep->setProblem(_x, problem);
-      EnergyNorm<Matrix, Vector> const norm(_matrix);
-      LoopSolver<Vector> solver(
-          multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
-          _localParset.get<double>("tolerance"), &norm,
-          _localParset.get<Solver::VerbosityMode>("verbosity"),
-          false); // absolute error
-
-      solver.preprocess();
-      solver.solve();
-    };
-
-    timeStep = 0;
-    relativeTime = 0.0;
-    relativeTau = 1e-6;
-
-    Vector ell0(u.size());
-    externalForces(relativeTime, ell0);
-
-    // Initial velocity
-    v = 0.0;
-
-    // Initial displacement: Start from a situation of minimal stress,
-    // which is automatically attained in the case [v = 0 = a].
-    // Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
-    solveLinearProblem(dirichletNodes, matrices.elasticity, ell0, u,
-                       parset.sub("u0.solver"));
-
-    // Initial acceleration: Computed in agreement with Ma = ell0 - Au
-    // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
-    Vector accelerationRHS = ell0;
-    Arithmetic::subtractProduct(accelerationRHS, matrices.elasticity, u);
-    solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a,
-                       parset.sub("a0.solver"));
-
-    // Initial state
-    alpha = parset.get<double>("boundary.friction.initialAlpha");
-
-    // Initial normal stress
-    myAssembler.assembleWeightedNormalStress(
-        frictionalBoundary, weightedNormalStress, body.getYoungModulus(),
-        body.getPoissonRatio(), u);
-  }
-
-public:
-  Vector u;
-  Vector v;
-  Vector a;
-  ScalarVector alpha;
-  ScalarVector weightedNormalStress;
-  double relativeTime;
-  double relativeTau;
-  size_t timeStep;
-};
-
-#endif

src/spatial-solving/CMakeLists.txt

+add_subdirectory("tnnmg")
+add_subdirectory("contact")
+add_subdirectory("preconditioners")

src/spatial-solving/contact/CMakeLists.txt

+add_custom_target(dune-tectonic_spatial-solving_contact_src SOURCES
+  dualmortarcoupling.cc
+  dualmortarcoupling.hh
+  nbodyassembler.cc
+  nbodyassembler.hh
+)

src/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 GridView0, class GridView1>
+void DualMortarCoupling<field_type, GridView0, GridView1>::assembleNonmortarLagrangeMatrix()
+{
+    // Create mapping from the global set of block dofs to the ones on the contact boundary
+    std::vector<int> globalToLocal;
+    nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
+
+    // clear matrix
+    nonmortarLagrangeMatrix_ = Dune::BDMatrix<MatrixBlock>(nonmortarBoundary_.numVertices());
+    nonmortarLagrangeMatrix_ = 0;
+
+    const auto& indexSet = gridView0_.indexSet();
+
+    // 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)];
+                Dune::MatrixVector::addToDiagonal(nonmortarLagrangeMatrix_[v][v],sfI);
+        }
+
+    }
+
+}
+
+template <class field_type, class GridView0, class GridView1>
+void DualMortarCoupling<field_type, GridView0, GridView1>::setup()
+{
+    typedef Dune::PQkLocalFiniteElementCache<typename GridType1::ctype, field_type, GridType1::dimension,1> FiniteElementCache1;
+
+    // cache for the dual functions on the boundary
+    using DualCache = Dune::Contact::DualBasisAdapter<GridView0, field_type>;
+
+    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);
+    };
+
+    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
+    BoundaryPatch0 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());
+
+
+    // Assemble the diagonal matrix coupling the nonmortar side with the lagrange multiplyers there
+    assembleNonmortarLagrangeMatrix();
+
+    // The weak obstacle vector
+    weakObstacle_.resize(nonmortarBoundary_.numVertices());
+    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(nonmortarBoundary_.numVertices(), grid1_->size(dim));
+
+    // Create mapping from the global set of block dofs to the ones on the contact boundary
+    std::vector<int> globalToLocal;
+    nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
+
+    // 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 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> >();
+
+    std::vector<Dune::FieldVector<ctype,dim> > avNormals;
+    avNormals = nonmortarBoundary_.getNormals();
+
+    // loop over all intersections and compute the matrix entries
+    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);
+        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();
+
+        int nNonmortarFaceNodes = domainRefElement.size(rIs.indexInInside(),1,dim);
+        std::vector<int> nonmortarFaceNodes;
+        for (int i=0; i<nNonmortarFaceNodes; i++) {
+          int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
+          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);
+            dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+            // loop over all Lagrange multiplier shape functions
+            for (int j=0; j<nNonmortarFaceNodes; 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();
+}

src/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_CONTACT_ASSEMBLERS_DUAL_MORTAR_COUPLING_HH
+#define DUNE_CONTACT_ASSEMBLERS_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 Extension of Dune::Contact::DualMortarCoupling that works on
+ *         any GridView (not just LeafGridView) and allows crosspoints.
+ *
+ *         Assembles the transfer operator for two-body contact
+ */
+template<class field_type, class GridView0, class GridView1=GridView0>
+class DualMortarCoupling {
+
+    // /////////////////////
+    // private types
+    // /////////////////////
+
+    enum {dim = GridView0::dimension};
+
+    enum {dimworld = GridView0::dimensionworld};
+
+  protected:
+    using MatrixBlock = Dune::FieldMatrix<field_type, dim, dim>;
+    using MatrixType = Dune::BCRSMatrix<MatrixBlock>;
+    using ObstacleVector = Dune::BlockVector<Dune::FieldVector<field_type, 1> >;
+
+    using ctype = typename GridView0::ctype;
+
+    using GridType0 = typename GridView0::Grid;
+    using GridType1 = typename GridView1::Grid;
+
+    using BoundaryPatch0 = BoundaryPatch<GridView0>;
+    using BoundaryPatch1 = BoundaryPatch<GridView1>;
+
+    using LevelBoundaryPatch0 = BoundaryPatch<typename GridType0::LevelGridView>;
+    using LevelBoundaryPatch1 = BoundaryPatch<typename GridType1::LevelGridView>;
+
+    using Extractor0 = Dune::GridGlue::Codim1Extractor<GridView0>;
+    using Extractor1 = Dune::GridGlue::Codim1Extractor<GridView1>;
+
+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(GridView0&& gridView0, GridView0&& gridView1,
+                       field_type overlap = 1e-2, field_type coveredArea = 1.0 - 1e-2)
+        : overlap_(overlap), gridGlueBackend_(nullptr), coveredArea_(coveredArea)
+    {
+        setGrids(gridView0, gridView1);
+    }
+
+    virtual ~DualMortarCoupling() {
+        /* Nothing. */
+    }
+
+    /** \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 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 nonmortar boundary. */
+    const BoundaryPatch0& nonmortarBoundary() const {
+        return nonmortarBoundary_;
+    }
+
+    /** \brief Return the mortar boundary. */
+    const BoundaryPatch1& mortarBoundary() const {
+        return mortarBoundary_;
+    }
+
+    /** \brief Set the non-mortar and mortar grid. */
+    void setGrids(GridView0&& gridView0, GridView1&& gridView1) {
+        gridView0_ = std::forward<GridView0>(gridView0);
+        gridView1_ = std::forward<GridView1>(gridView1);
+
+        grid0_ = &gridView0_.grid();
+        grid1_ = &gridView1_.grid();
+    }
+
+    /** \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_;
+    }
+
+
+    // /////////////////////////////////////////////
+    //   Data members
+    // /////////////////////////////////////////////
+private:
+    /** \brief The two grids involved in the two-body contact problem. */
+    const GridView0& gridView0_;
+    const GridView1& gridView1_;
+
+    const GridType0* grid0_;
+    const GridType1* grid1_;
+
+    /** \brief For each dof a bit specifying whether the dof carries an obstacle or not. */
+    BoundaryPatch0 nonmortarBoundary_;
+
+    /** \brief The mortar boundary. */
+    BoundaryPatch1 mortarBoundary_;
+
+    /** \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_;
+
+    /** \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

src/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>
+
+namespace Dune {
+namespace Contact {
+
+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]->setup();
+    }
+
+    // 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
+    // //////////////////////////////////////////////////////////////////
+
+    std::vector<std::vector<MatrixBlock> > coordinateSystems(nCouplings());
+
+    for (int i=0; i<nCouplings(); i++) {
+
+        double dist = coupling_[i].obsDistance_;
+        auto projection = coupling_[i].projection();
+
+        if (!projection)
+            DUNE_THROW(Dune::Exception, "You have to supply a contact projection for the " << i << "-th coupling!");
+
+        std::vector<Dune::FieldVector<field_type,dim> > directions;
+
+        const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
+        const auto& mortarBoundary    = contactCoupling_[i]->mortarBoundary();
+
+        projection->project(nonmortarBoundary, mortarBoundary,dist);
+        projection->getObstacleDirections(directions);
+
+        this->computeLocalCoordinateSystems(nonmortarBoundary,coordinateSystems[i], directions);
+    }
+
+    // ///////////////////////////////////////////////////////////////
+    // Combine the coordinate systems for each grid to one long array
+    // ///////////////////////////////////////////////////////////////
+    this->localCoordSystems_.resize(numDofs());
+
+    // initialise with identity
+    Dune::ScaledIdentityMatrix<field_type,dim> id(1);
+    for (size_t i=0; i<this->localCoordSystems_.size(); i++)
+        this->localCoordSystems_[i] = id;
+
+    for (int j=0; j<nCouplings(); j++) {
+
+        int grid0Idx = coupling_[j].gridIdx_[0];
+
+        // compute offset
+        int idx = 0;
+        for (int k=0;k<grid0Idx;k++)
+            idx += grids_[k]->size(dim);
+
+        const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
+
+        // if a body is non-mortar more then once, one has to be careful with not overwriting entries
+        for (std::size_t k=0; k<coordinateSystems[j].size(); k++)
+          if(nonmortarBoundary.containsVertex(k))
+            this->localCoordSystems_[idx+k] = coordinateSystems[j][k];
+    }
+}
+
+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)
+        contactCoupling_[i]->nonmortarBoundary().makeGlobalToLocal(globalToLocals[i]);
+
+    ///////////////////////////////////////
+    //  Set the obstacle values
+    ///////////////////////////////////////
+
+    totalHasObstacle_.resize(numDofs(),false);
+
+    for (int j=0; j<nCouplings(); j++) {
+
+        int grid0Idx = coupling_[j].gridIdx_[0];
+        const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
+
+        // 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 (nonmortarBoundary.containsVertex(k))
+                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)
+                continue;
+
+            // Set the obstacle
+            switch (coupling_[i].type_) {
+
+            case CouplingPairBase::CONTACT:
+                totalObstacles_[idx+vIdx].upper(0) = obs[globalToLocal[vIdx]];
+                break;
+
+            case CouplingPairBase::STICK_SLIP:
+                totalObstacles_[idx+vIdx].lower(0) = 0;
+                totalObstacles_[idx+vIdx].upper(0) = 0;
+                break;
+
+            case CouplingPairBase::GLUE:
+                for (int j=0; j<dim; j++) {
+                    totalObstacles_[idx+vIdx].lower(j) = 0;
+                    totalObstacles_[idx+vIdx].upper(j) = 0;
+                }
+                break;
+
+            case CouplingPairBase::NONE:
+                break;
+            default:
+                DUNE_THROW(Dune::NotImplemented, "Coupling type not handled yet!");
+            }
+        }
+    }
+
+}
+
+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();
+        // 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 (nonmortarBoundary.containsVertex(j))
+                nonmortarToGlobal[i][idx++] = j;
+
+        // 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();
+        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(nonmortarBoundary.containsVertex(k))
+                BT_[offset + k][offset+k] =this->localCoordSystems_[offset + k];
+    }
+}
+
+} /* namespace Contact */
+} /* namespace Dune */

src/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_CONTACT_ASSEMBLERS_N_BODY_MORTAR_ASSEMBLER_HH
+#define DUNE_CONTACT_ASSEMBLERS_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 <dune/contact/assemblers/dualmortarcoupling.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.
+     *   \param hierarchyCoupling This boolean determines if the whole hierarchy of mortar matrices should be setup
+     */
+    NBodyAssembler(int nGrids, int nCouplings, bool hierarchyCoupling=false, field_type coveredArea = 0.99) :
+      coveredArea_(coveredArea)
+    {
+        grids_.resize(nGrids);
+
+        coupling_.resize(nCouplings);
+        contactCoupling_.resize(nCouplings);
+        // initialise the couplings with DualMortarCoupling
+        for (int i=0; i<nCouplings; i++)
+            if (hierarchyCoupling)
+                contactCoupling_[i] = std::make_shared<HierarchyCouplingType>();
+            else
+                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;
+    }
+
+    /** \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 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];}
+
+protected:
+    /** \brief Compute the transposed mortar transformation matrix. */
+    void computeTransformationMatrix();
+
+    /** \brief Setup the Householder reflections. */
+    void assembleHouseholderReflections();
+
+public:
+    /** \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_;
+
+    /** \brief The mortar couplings between grid pairs */
+    std::vector<std::shared_ptr<CouplingType> > contactCoupling_;
+
+    /** \brief The coupling pairs. */
+    std::vector<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

src/spatial-solving/fixedpointiterator.cc

@@ -4,84 +4,256 @@
 
 #include <dune/common/exceptions.hh>
 
-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
 
-#include "../enums.hh"
-#include "../enumparser.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 "fixedpointiterator.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/preconditioners/nbodycontacttransfer.hh"
+
 void FixedPointIterationCounter::operator+=(
     FixedPointIterationCounter const &other) {
   iterations += other.iterations;
   multigridIterations += other.multigridIterations;
 }
 
-template <class Factory, class Updaters, class ErrorNorm>
-FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
-    Factory &factory, Dune::ParameterTree const &parset,
-    std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
-    : step_(factory.getStep()),
-      parset_(parset),
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::FixedPointIterator(
+    Dune::ParameterTree const &parset,
+    const ContactNetwork& contactNetwork,
+    const IgnoreVector& ignoreNodes,
+    GlobalFriction& globalFriction,
+    const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+    const ErrorNorms& errorNorms)
+    : parset_(parset),
+      contactNetwork_(contactNetwork),
+      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")),
       velocityMaxIterations_(parset.get<size_t>("v.solver.maximumIterations")),
       velocityTolerance_(parset.get<double>("v.solver.tolerance")),
       verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
-      errorNorm_(errorNorm) {}
+      errorNorms_(errorNorms) {}
 
-template <class Factory, class Updaters, class ErrorNorm>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
 FixedPointIterationCounter
-FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
-    Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
-    Vector &velocityIterate) {
-  EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
-  LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
-                                           velocityTolerance_, &energyNorm,
+FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
+    Updaters updaters,
+    const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+    std::vector<Vector>& velocityIterates) {
+
+  std::cout << "FixedPointIterator::run()" << std::endl;
+
+  const auto& nBodyAssembler = contactNetwork_.nBodyAssembler();
+
+  const auto nBodies = nBodyAssembler.nGrids();
+
+  std::vector<const Matrix*> matrices_ptr(nBodies);
+  for (int i=0; i<nBodies; i++) {
+      matrices_ptr[i] = &velocityMatrices[i];
+  }
+
+  // assemble full global contact problem
+  Matrix bilinearForm;
+  nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+  //print(bilinearForm, "bilinearForm:");
+
+  Vector totalRhs;
+  nBodyAssembler.assembleRightHandSide(velocityRHSs, totalRhs);
+
+  //print(totalRhs, "totalRhs:");
+
+  // get lower and upper obstacles
+  const auto& totalObstacles = nBodyAssembler.totalObstacles_;
+  Vector lower(totalObstacles.size());
+  Vector upper(totalObstacles.size());
+
+  for (size_t j=0; j<totalObstacles.size(); ++j) {
+      const auto& totalObstaclesj = totalObstacles[j];
+      auto& lowerj = lower[j];
+      auto& upperj = upper[j];
+    for (size_t d=0; d<dims; ++d) {
+        lowerj[d] = totalObstaclesj[d][0];
+        upperj[d] = totalObstaclesj[d][1];
+    }
+  }
+
+  //print(totalObstacles, "totalObstacles:");
+
+  //print(lower, "lower obstacles:");
+  //print(upper, "upper obstacles:");
+
+  // comput velocity obstacles
+  Vector vLower, vUpper;
+  std::vector<Vector> u0, v0;
+  updaters.rate_->extractOldVelocity(v0);
+  updaters.rate_->extractOldDisplacement(u0);
+
+  Vector totalu0, totalv0;
+  nBodyAssembler.concatenateVectors(u0, totalu0);
+  nBodyAssembler.concatenateVectors(v0, totalv0);
+
+  updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
+  updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);
+
+  //print(vLower, "vLower obstacles:");
+  //print(vUpper, "vUpper obstacles:");
+
+  std::cout << "- Problem assembled: success" << std::endl;
+
+  using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, IgnoreVector>;
+  using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+  using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+  TransferOperators transfer(contactNetwork_.nLevels()-1);
+  for (size_t i=0; i<transfer.size(); i++) {
+      transfer[i] = std::make_shared<TransferOperator>();
+      transfer[i]->setup(contactNetwork_, i, i+1);
+  }
+
+  // Remove any recompute filed so that initially the full transferoperator is assembled
+  for (size_t i=0; i<transfer.size(); i++)
+      std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+  auto smoother = TruncatedBlockGSStep<Matrix, Vector>{};
+  auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+  linearMultigridStep->setMGType(1, 3, 3);
+  linearMultigridStep->setSmoother(&smoother);
+  linearMultigridStep->setTransferOperators(transfer);
+
+  EnergyNorm<Matrix, Vector> mgNorm(*linearMultigridStep);
+  LinearSolver mgSolver(linearMultigridStep, parset_.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset_.get<double>("solver.tnnmg.linear.tolerance"), mgNorm, Solver::QUIET);
+
+
+  // set up functional and TNMMG solver
+  Functional J(bilinearForm, totalRhs, globalFriction_, vLower, vUpper);
+  Factory solverFactory(parset_.sub("solver.tnnmg"), J, mgSolver, ignoreNodes_);
+  auto step = solverFactory.step();
+
+  std::cout << "- Functional and TNNMG step setup: success" << std::endl;
+
+  EnergyNorm<Matrix, Vector> energyNorm(bilinearForm);
+  LoopSolver<Vector> velocityProblemSolver(*step.get(), velocityMaxIterations_,
+                                           velocityTolerance_, energyNorm,
                                            verbosity_, false); // absolute error
 
   size_t fixedPointIteration;
   size_t multigridIterations = 0;
-  ScalarVector alpha;
+  std::vector<ScalarVector> alpha(nBodies);
   updaters.state_->extractAlpha(alpha);
   for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
        ++fixedPointIteration) {
+
+    // contribution from nonlinearity
+    globalFriction_.updateAlpha(alpha);
+
+    Vector totalVelocityIterate;
+    nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+
+    //print(velocityIterates, "velocityIterates:");
+    //print(totalVelocityIterate, "totalVelocityIterate:");
+    std::cout << "- FixedPointIteration iterate" << std::endl;
+
     // solve a velocity problem
-    globalFriction_->updateAlpha(alpha);
-    ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
-                                velocityRHS, velocityIterate);
-    BlockProblem velocityProblem(parset_, convexProblem);
-    step_->setProblem(velocityIterate, velocityProblem);
+    solverFactory.setProblem(totalVelocityIterate);
+
+    std::cout << "- Velocity problem set" << std::endl;
+
     velocityProblemSolver.preprocess();
+    std::cout << "-- Preprocessed" << std::endl;
     velocityProblemSolver.solve();
+    std::cout << "-- Solved" << std::endl;
+
+    const auto& tnnmgSol = step->getSol();
+    nBodyAssembler.postprocess(tnnmgSol, velocityIterates);
+    //nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
+
+    //print(totalVelocityIterate, "totalVelocityIterate:");
+    //print(velocityIterates, "velocityIterates:");
 
     multigridIterations += velocityProblemSolver.getResult().iterations;
 
-    Vector v_m;
+    std::vector<Vector> v_m;
     updaters.rate_->extractOldVelocity(v_m);
-    v_m *= 1.0 - lambda_;
-    Arithmetic::addProduct(v_m, lambda_, velocityIterate);
+
+    for (size_t i=0; i<v_m.size(); i++) {
+      v_m[i] *= 1.0 - lambda_;
+      Dune::MatrixVector::addProduct(v_m[i], lambda_, velocityIterates[i]);
+    }
+
+    // extract relative velocities in mortar basis
+    std::vector<Vector> v_rel;
+    relativeVelocities(tnnmgSol, v_rel);
+
+    //print(v_rel, "v_rel");
+
+    std::cout << "- State problem set" << std::endl;
 
     // solve a state problem
-    updaters.state_->solve(v_m);
-    ScalarVector newAlpha;
+    updaters.state_->solve(v_rel);
+
+    std::cout << "-- Solved" << std::endl;
+
+    std::vector<ScalarVector> newAlpha(nBodies);
     updaters.state_->extractAlpha(newAlpha);
 
-    if (lambda_ < 1e-12 or
-        errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
+    bool breakCriterion = true;
+    for (size_t i=0; i<nBodies; i++) {
+        if (alpha[i].size()==0 || newAlpha[i].size()==0)
+            continue;
+
+        //print(alpha[i], "alpha i:");
+        //print(newAlpha[i], "new alpha i:");
+        if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
+            breakCriterion = false;
+            std::cout << "fixedPoint error: " << errorNorms_[i]->diff(alpha[i], newAlpha[i]) << std::endl;
+            break;
+        }
+    }
+
+    if (lambda_ < 1e-12 or breakCriterion) {
+      std::cout << "-FixedPointIteration finished! " << (lambda_ < 1e-12 ? "lambda" : "breakCriterion") << std::endl;
       fixedPointIteration++;
       break;
     }
     alpha = newAlpha;
   }
+
+  std::cout << "-FixedPointIteration finished! " << std::endl;
+
   if (fixedPointIteration == fixedPointMaxIterations_)
     DUNE_THROW(Dune::Exception, "FPI failed to converge");
 
-  updaters.rate_->postProcess(velocityIterate);
+  updaters.rate_->postProcess(velocityIterates);
 
   // Cannot use return { fixedPointIteration, multigridIterations };
   // with gcc 4.9.2, see also http://stackoverflow.com/a/37777814/179927
@@ -97,4 +269,60 @@ std::ostream &operator<<(std::ostream &stream,
                 << ")";
 }
 
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+void FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::relativeVelocities(
+    const Vector& v,
+    std::vector<Vector>& v_rel) const {
+
+    const auto& nBodyAssembler = contactNetwork_.nBodyAssembler();
+    const size_t nBodies = nBodyAssembler.nGrids();
+   // const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+
+    size_t globalIdx = 0;
+    v_rel.resize(nBodies);
+    for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+        const auto& nonmortarBoundary = *bodywiseNonmortarBoundaries_[bodyIdx];
+        auto& v_rel_ = v_rel[bodyIdx];
+
+        v_rel_.resize(nonmortarBoundary.size());
+
+        for (size_t i=0; i<v_rel_.size(); i++) {
+            if (toBool(nonmortarBoundary[i])) {
+                v_rel_[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

@@ -4,10 +4,15 @@
 #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>
+
 struct FixedPointIterationCounter {
   size_t iterations = 0;
   size_t multigridIterations = 0;
@@ -18,29 +23,48 @@ struct FixedPointIterationCounter {
 std::ostream &operator<<(std::ostream &stream,
                          FixedPointIterationCounter const &fpic);
 
-template <class Factory, class Updaters, class ErrorNorm>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
 class FixedPointIterator {
+  using Functional = typename Factory::Functional;
   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 relativeVelocities(const Vector& v, std::vector<Vector>& v_rel) const;
 
 public:
-  FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
-                     std::shared_ptr<Nonlinearity> globalFriction,
-                     ErrorNorm const &errorNorm_);
+  FixedPointIterator(const Dune::ParameterTree& parset,
+                     const ContactNetwork& contactNetwork,
+                     const IgnoreVector& ignoreNodes,
+                     GlobalFriction& globalFriction,
+                     const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+                     const ErrorNorms& errorNorms);
 
   FixedPointIterationCounter run(Updaters updaters,
-                                 Matrix const &velocityMatrix,
-                                 Vector const &velocityRHS,
-                                 Vector &velocityIterate);
+                                 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 ContactNetwork& contactNetwork_;
+  const IgnoreVector& ignoreNodes_;
+
+  GlobalFriction& globalFriction_;
+  const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
+
 
   size_t fixedPointMaxIterations_;
   double fixedPointTolerance_;
@@ -48,6 +72,6 @@ class FixedPointIterator {
   size_t velocityMaxIterations_;
   double velocityTolerance_;
   Solver::VerbosityMode verbosity_;
-  ErrorNorm const &errorNorm_;
+  const ErrorNorms& errorNorms_;
 };
 #endif

src/spatial-solving/fixedpointiterator_tmpl.cc

@@ -5,28 +5,23 @@
 #include "../explicitgrid.hh"
 #include "../explicitvectors.hh"
 
-#include <dune/common/function.hh>
-
 #include <dune/solvers/norms/energynorm.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
 
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
+#include "../spatial-solving/solverfactory_tmpl.cc"
+#include "../data-structures/contactnetwork_tmpl.cc"
 
 #include "../time-stepping/rate/rateupdater.hh"
 #include "../time-stepping/state/stateupdater.hh"
 #include "../time-stepping/updaters.hh"
-#include "solverfactory.hh"
 
-using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+
 using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
 using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
 
-using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
+
 
-template class FixedPointIterator<Factory, MyUpdaters, ErrorNorm>;
+template class FixedPointIterator<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;

src/spatial-solving/preconditioners/CMakeLists.txt

+add_custom_target(dune-tectonic_spatial-solving_preconditioners_src SOURCES
+  hierarchicleveliterator.hh
+  levelglobalpreconditioner.hh
+  levelpatchpreconditioner.hh
+  multilevelpatchpreconditioner.hh
+  nbodycontacttransfer.hh
+  supportpatchfactory.hh
+)

src/spatial-solving/preconditioners/hierarchicleveliterator.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef HIERARCHIC_LEVEL_ITERATOR_HH
+#define HIERARCHIC_LEVEL_ITERATOR_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/iteratorfacades.hh>
+
+#include <dune/geometry/multilineargeometry.hh>
+#include <dune/geometry/referenceelements.hh>
+
+/** \brief Hierarchic leaf iterator.
+ *
+ *  This iterator loops over all children of a given coarse element and only returns the ones that are leaf.
+ *  If the starting element is leaf itself, then the returned iterator returns the element itself.
+ *  This class also provides a geometry, mapping local coordinates of the children to local coordinates
+ *  in the coarse element. 
+*/
+template <class GridImp>
+class HierarchicLevelIterator :
+    public Dune::ForwardIteratorFacade<HierarchicLevelIterator<GridImp>, const typename GridImp::template Codim<0>::Entity>
+{
+    typedef typename GridImp::template Codim<0>::Entity Element;
+    typedef typename GridImp::HierarchicIterator HierarchicIterator;
+    enum {dim = GridImp::dimension};
+    enum {dimworld = GridImp::dimensionworld};
+public:
+    typedef Dune::CachedMultiLinearGeometry<typename GridImp::ctype, dim, dimworld> LocalGeometry;
+
+    enum PositionFlag {begin, end};
+
+    HierarchicLevelIterator(const GridImp& grid, const Element& element, 
+                            PositionFlag flag, int maxlevel, bool nested = true)
+        : element_(element), maxlevel_(maxlevel), hIt_(element.hend(maxlevel_)),
+            flag_(flag), nested_(nested)
+    {
+
+        // if the element itself is leaf, then we don't have to iterate over the children
+        if (flag==begin && !(element_.level()==maxlevel_)) {
+            hIt_ = element_.hbegin(maxlevel_);
+
+            //NOTE This class by now assumes that possible non-nestedness of the grid levels only arises
+            // due to boundary parametrisation
+            if (!nested_)  {
+                // Check if the element is a boundary element, and set the nested flag correspondingly
+                // If the element is not at the boundary, then we can neglect the nested flag
+                typename GridImp::LevelIntersectionIterator lIt = grid.levelGridView(element.level()).ibegin(element);
+                typename GridImp::LevelIntersectionIterator lEndIt = grid.levelGridView(element.level()).ibegin(element);
+                nested_ = true;
+                for (; lIt != lEndIt; ++lIt)
+                    if (lIt->boundary()) {
+                        nested_ = false;
+                        break;
+                    }
+            }
+
+            if (!(hIt_->level()==maxlevel_))
+                increment();
+        }
+    }
+
+    //! Copy constructor
+    HierarchicLevelIterator(const HierarchicLevelIterator<GridImp>& other) :
+        element_(other.element_),maxlevel_(other.maxlevel_), hIt_(other.hIt_),
+        flag_(other.flag_), nested_(other.nested_)
+    {}
+
+    //! Equality
+    bool equals(const HierarchicLevelIterator<GridImp>& other) const {
+        return  (element_ == other.element_)  && 
+                ((flag_==other.flag_ && flag_==end) || (hIt_ == other.hIt_ && flag_==begin && other.flag_==flag_));
+    }
+
+    //! Prefix increment
+    void increment() {
+
+        HierarchicIterator hEndIt = element_.hend(maxlevel_);
+
+        if (hIt_ == hEndIt) {
+            flag_ = end;
+            return;
+        }
+
+        // Increment until we hit a leaf child element
+        do {
+            ++hIt_;
+        } while(hIt_ != hEndIt && (!(hIt_->level()==maxlevel_)));
+    
+        if (hIt_ == hEndIt)
+            flag_ = end;
+    }
+
+    //! Dereferencing
+    const Element& dereference() const {
+        if (flag_ == end)
+            DUNE_THROW(Dune::Exception,"HierarchicLevelIterator: Cannot dereference end iterator!");
+        return (element_.level()==maxlevel_) ? element_ : *hIt_;
+    }
+
+    //! Compute the local geometry mapping from the leaf child to the starting element
+    LocalGeometry geometryInAncestor() {
+
+        //NOTE: We assume here that the type of the child and the ancestors are the same!
+        const Dune::ReferenceElement<typename GridImp::ctype, dim>& ref = Dune::ReferenceElements<typename GridImp::ctype,dim>::general(element_.type());
+
+        // store local coordinates of the leaf child element within the coarse starting element
+        std::vector<Dune::FieldVector<typename GridImp::ctype,dim> > corners(ref.size(dim));
+        for (int i=0; i<corners.size(); i++)
+            corners[i] = ref.position(i,dim);
+
+        // If the element is leaf, then return an Identity mapping
+        if (element_.level()==maxlevel_)
+            return LocalGeometry(ref,corners);
+
+        if (nested_) {
+            const typename Element::Geometry leafGeom = hIt_->geometry();
+            const typename Element::Geometry coarseGeom = element_.geometry();
+
+            for (int i=0; i<corners.size();i++) 
+                corners[i] = coarseGeom.local(leafGeom.corner(i));
+            return LocalGeometry(ref,corners);
+        }
+ 
+        Element father(*hIt_);
+        while (father != element_) {
+
+            const typename Element::LocalGeometry fatherGeom = hIt_->geometryInFather();
+            father = father->father();
+
+            for (int i=0; i<corners.size(); i++)
+                corners[i] = fatherGeom.global(corners[i]);
+        }
+
+        return LocalGeometry(ref,corners);
+    }
+
+private:
+
+    //! The starting element
+    const Element element_;
+
+    //! The grids maxlevel
+    int maxlevel_;
+
+    //! The actual hierarchic iterator
+    HierarchicIterator hIt_;
+
+    //! Position flag
+    PositionFlag flag_;
+
+    //! Flag that is set true if the grid levels are conforming (i.e. no parametrised boundaries)
+    bool nested_;
+
+};
+#endif