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/fixedpointiterator.cc

@@ -4,84 +4,224 @@
 
 #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"
+
+
 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 NBodyAssembler, class Updaters, class ErrorNorm>
+FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::FixedPointIterator(
+    Dune::ParameterTree const &parset,
+    NBodyAssembler& nBodyAssembler,
+    const IgnoreVector& ignoreNodes,
+    GlobalFriction& globalFriction,
+    const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+    const std::vector<const ErrorNorm* >& errorNorms)
+    : parset_(parset),
+      nBodyAssembler_(nBodyAssembler),
+      ignoreNodes_(ignoreNodes),
       globalFriction_(globalFriction),
+      bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
       fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
       fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
       lambda_(parset.get<double>("v.fpi.lambda")),
       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 NBodyAssembler, class Updaters, class ErrorNorm>
 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, NBodyAssembler, Updaters, ErrorNorm>::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 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;
+
+  // set up functional and TNMMG solver
+  Functional J(bilinearForm, totalRhs, globalFriction_, vLower, vUpper);
+  Factory solverFactory(parset_.sub("solver.tnnmg"), J, 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 +237,59 @@ std::ostream &operator<<(std::ostream &stream,
                 << ")";
 }
 
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::relativeVelocities(
+    const Vector& v,
+    std::vector<Vector>& v_rel) const {
+
+    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 NBodyAssembler, class Updaters, class ErrorNorm>
 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,
+                     NBodyAssembler& nBodyAssembler,
+                     const IgnoreVector& ignoreNodes,
+                     GlobalFriction& globalFriction,
+                     const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+                     const std::vector<const ErrorNorm* >& 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_;
+  NBodyAssembler& nBodyAssembler_;
+  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 std::vector<const ErrorNorm* >& errorNorms_;
 };
 #endif

src/spatial-solving/fixedpointiterator_tmpl.cc

@@ -5,28 +5,24 @@
 #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/levelcontactnetwork_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 NBodyAssembler = typename MyLevelContactNetwork::NBodyAssembler;
+using BoundaryNodes = typename MyLevelContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyLevelContactNetwork::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>;
 
-template class FixedPointIterator<Factory, MyUpdaters, ErrorNorm>;
+
+template class FixedPointIterator<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorm>;

src/spatial-solving/solverfactory.cc

@@ -2,58 +2,48 @@
 #include "config.h"
 #endif
 
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#include <dune/fufem/assemblers/transferoperatorassembler.hh>
-#include <dune/solvers/solvers/solver.hh>
+#include <dune/solvers/common/wrapownshare.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+#include <dune/solvers/solvers/umfpacksolver.hh>
 
 #include "solverfactory.hh"
 
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
-    Dune::ParameterTree const &parset, Grid const &grid,
-    Dune::BitSetVector<dim> const &ignoreNodes)
-    : baseEnergyNorm(linearBaseSolverStep),
-      linearBaseSolver(&linearBaseSolverStep,
-                       parset.get<size_t>("linear.maxiumumIterations"),
-                       parset.get<double>("linear.tolerance"), &baseEnergyNorm,
-                       Solver::QUIET),
-      transferOperators(grid.maxLevel()),
-      multigridStep(
-          std::make_shared<Step>(linearIterationStep, nonlinearSmoother)) {
-  // linear iteration step
-  linearIterationStep.setMGType(parset.get<int>("linear.cycle"),
-                                parset.get<int>("linear.pre"),
-                                parset.get<int>("linear.post"));
-  linearIterationStep.basesolver_ = &linearBaseSolver;
-  linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);
-
-  // transfer operators
-  for (auto &&x : transferOperators)
-    x = new CompressedMultigridTransfer<Vector>;
-  TransferOperatorAssembler<Grid>(grid)
-      .assembleOperatorPointerHierarchy(transferOperators);
-  linearIterationStep.setTransferOperators(transferOperators);
-
-  // tnnmg iteration step
-  multigridStep->setSmoothingSteps(parset.get<int>("main.pre"),
-                                   parset.get<int>("main.multi"),
-                                   parset.get<int>("main.post"));
-  multigridStep->ignoreNodes_ = &ignoreNodes;
+#include "../utils/debugutils.hh"
+
+template <class Functional, class BitVector>
+SolverFactory<Functional, BitVector>::SolverFactory(
+    const Dune::ParameterTree& parset,
+    Functional& J,
+    const BitVector& ignoreNodes) :
+      J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J)))
+      //linear solver
+      //preconditioner_(),
+      //linearSolverStep_(),
+    {
+    nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
+
+    // linearSolverStep_.setPreconditioner(preconditioner_);
+    linearSolverStep_ = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::ldlt());
+    energyNorm_ = std::make_shared<EnergyNorm<Matrix, Vector>>(*linearSolverStep_);
+    //linearSolver_ = std::make_shared<LinearSolver>(*linearSolverStep_, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), *energyNorm_, Solver::QUIET);
+    linearSolver_ = std::make_shared<Dune::Solvers::UMFPackSolver<Matrix, Vector>>();
+
+    tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_, DefectProjection(), LineSearchSolver());
+    tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
+    tnnmgStep_->setIgnore(ignoreNodes);
 }
 
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() {
-  for (auto &&x : transferOperators)
-    delete x;
+template <class Functional, class BitVector>
+void SolverFactory<Functional, BitVector>::setProblem(Vector& x) {
+    nonlinearSmoother_->setProblem(x);
+    tnnmgStep_->setProblem(x);
 }
 
-template <size_t dim, class BlockProblem, class Grid>
-auto SolverFactory<dim, BlockProblem, Grid>::getStep()
-    -> std::shared_ptr<Step> {
-  return multigridStep;
+
+template <class Functional, class BitVector>
+auto SolverFactory<Functional, BitVector>::step()
+-> std::shared_ptr<Step> {
+    return tnnmgStep_;
 }
 
 #include "solverfactory_tmpl.cc"

src/spatial-solving/solverfactory.hh

 #ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
 #define SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
 
+#define NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
+
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/parametertree.hh>
 
-#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
-#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+//#include <dune/solvers/iterationsteps/cgstep.hh>
+//#include <dune/solvers/iterationsteps/amgstep.hh>
+
+//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
+#include <dune/tnnmg/projections/obstacledefectprojection.hh>
+#include <dune/tnnmg/linearsolvers/fastamgsolver.hh>
+
+//#include "tnnmg/functional.hh"
+#include "tnnmg/tnnmgstep.hh"
+#include "tnnmg/linearization.hh"
+#include "tnnmg/linesearchsolver.hh"
+#include "tnnmg/localbisectionsolver.hh"
 
-template <size_t dim, class BlockProblemTEMPLATE, class Grid>
+template <class FunctionalTEMPLATE, class BitVectorType>
 class SolverFactory {
 public:
-  using BlockProblem = BlockProblemTEMPLATE;
-  using ConvexProblem = typename BlockProblem::ConvexProblemType;
-  using Vector = typename BlockProblem::VectorType;
-  using Matrix = typename BlockProblem::MatrixType;
+    using Functional = FunctionalTEMPLATE;
+    using Matrix = typename Functional::Matrix;
+    using Vector = typename Functional::Vector;
+    using BitVector = BitVectorType;
 
-private:
-  using NonlinearSmoother = GenericNonlinearGS<BlockProblem>;
+    using LocalSolver = LocalBisectionSolver;
+    using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, LocalSolver>;
 
-public:
-  using Step =
-      TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
+    using Linearization = Linearization<Functional, BitVector>;
+    //using Preconditioner = Dune::Solvers::Preconditioner; //TODO Platzhalter für PatchPreconditioner
+    //using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
+    //using LinearSolverStep = typename Dune::Solvers::AMGStep<Matrix, Vector>;
+    using LinearSolverStep = LinearIterationStep<Matrix, Vector>;
+    //using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using LinearSolver = typename Dune::Solvers::LinearSolver<Matrix, Vector>;
+
+    //using LinearSolver = typename Dune::Solvers::UMFPackSolver<Matrix, Vector>;
+    //using LinearSolver = typename Dune::TNNMG::FastAMGSolver<Matrix, Vector>;
 
-  SolverFactory(Dune::ParameterTree const &parset, Grid const &grid,
-                Dune::BitSetVector<dim> const &ignoreNodes);
+    using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
 
-  ~SolverFactory();
+    using Step = typename Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
 
-  std::shared_ptr<Step> getStep();
+
+  SolverFactory(Dune::ParameterTree const &,
+                Functional&,
+                const BitVector&);
+
+  void setProblem(Vector& x);
+
+  auto step() -> std::shared_ptr<Step>;
 
 private:
-  TruncatedBlockGSStep<Matrix, Vector> linearBaseSolverStep;
-  EnergyNorm<Matrix, Vector> baseEnergyNorm;
-  LoopSolver<Vector> linearBaseSolver;
-  TruncatedBlockGSStep<Matrix, Vector> linearPresmoother;
-  TruncatedBlockGSStep<Matrix, Vector> linearPostsmoother;
-  MultigridStep<Matrix, Vector> linearIterationStep;
-  std::vector<CompressedMultigridTransfer<Vector> *> transferOperators;
-  NonlinearSmoother nonlinearSmoother;
-  std::shared_ptr<Step> multigridStep;
+  Vector dummyIterate_;
+  std::shared_ptr<const Functional> J_;
+
+  // nonlinear smoother
+  std::shared_ptr<NonlinearSmoother> nonlinearSmoother_;
+
+  // linear solver
+  //Preconditioner preconditioner_;
+  std::shared_ptr<LinearSolverStep> linearSolverStep_;
+  std::shared_ptr<EnergyNorm<Matrix, Vector>> energyNorm_;
+  std::shared_ptr<LinearSolver> linearSolver_;
+
+  // TNNMGStep
+  std::shared_ptr<Step> tnnmgStep_;
 };
 #endif

src/spatial-solving/solverfactory_old.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
+#include <dune/solvers/solvers/solver.hh>
+
+#include "solverfactory_old.hh"
+
+
+template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
+SolverFactoryOld<DeformedGrid, Nonlinearity, MatrixType, VectorType>::SolverFactoryOld(
+    const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
+    const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes)
+    : nBodyAssembler_(nBodyAssembler),
+      baseEnergyNorm_(baseSolverStep_),
+      baseSolver_(&baseSolverStep_,
+                       parset.get<size_t>("linear.maximumIterations"),
+                       parset.get<double>("linear.tolerance"), &baseEnergyNorm_,
+                       Solver::QUIET),
+      transferOperators_(nBodyAssembler.getGrids().at(0)->maxLevel()) {
+
+ /* baseSolverStep_.obstacles_ = using HasObstacle = Dune::BitSetVector<BlockSize>;
+    using Obstacle = BoxConstraint<Field, BlockSize>;
+    const HasObstacle* hasObstacle_;
+    const std::vector<Obstacle>* obstacles_;*/
+
+  multigridStep_ = std::make_shared<Step>();
+  // tnnmg iteration step
+  multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
+  multigridStep_->setIgnore(ignoreNodes);
+  multigridStep_->setBaseSolver(baseSolver_);
+  multigridStep_->setSmoother(&presmoother_, &postsmoother_);
+  multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
+  multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
+  multigridStep_->setVerbosity(NumProc::QUIET);
+
+  // create the transfer operators
+  const int nCouplings = nBodyAssembler_.nCouplings();
+  const auto grids = nBodyAssembler_.getGrids();
+
+  for (size_t i=0; i<transferOperators_.size(); i++) {
+    transferOperators_[i] = std::make_shared<TransferOperator>();
+  }
+
+  std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nCouplings);
+  std::vector<const MatrixType*> mortarTransfer(nCouplings);
+  std::vector<std::array<int,2> > gridIdx(nCouplings);
+
+  for (int j=0; j<nCouplings; j++) {
+    fineHasObstacle[j]    = nBodyAssembler_.contactCoupling_[j]->nonmortarBoundary().getVertices();
+    mortarTransfer[j] = &nBodyAssembler_.contactCoupling_[j]->mortarLagrangeMatrix();
+    gridIdx[j]        = nBodyAssembler_.coupling_[j].gridIdx_;
+  }
+
+  TransferOperator::setupHierarchy(transferOperators_,
+                                   grids,
+                                   mortarTransfer,
+                                   nBodyAssembler_.localCoordSystems_,
+                                   fineHasObstacle,
+                                   gridIdx);
+
+  multigridStep_->setTransferOperators(transferOperators_);
+}
+
+template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
+auto SolverFactoryOld<DeformedGrid, Nonlinearity, MatrixType, VectorType>::getStep()
+    -> std::shared_ptr<Step> {
+  return multigridStep_;
+}
+
+#include "solverfactory_tmpl_old.cc"

src/spatial-solving/solverfactory_old.hh

+#ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORYOLD_HH
+#define SRC_SPATIAL_SOLVING_SOLVERFACTORYOLD_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/blockgsstep.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+//#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
+//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
+//#include <dune/contact/estimators/hierarchiccontactestimator.hh>
+#include <dune/contact/solvers/nsnewtonmgstep.hh>
+#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
+#include <dune/contact/solvers/contactobsrestrict.hh>
+#include <dune/contact/solvers/contacttransferoperatorassembler.hh>
+#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
+
+#define USE_OLD_TNNMG //needed for old bisection.hh in tnnmg
+
+template <class DeformedGridTEMPLATE, class NonlinearityTEMPLATE, class MatrixType, class VectorType>
+class SolverFactoryOld {
+//protected:
+//  const size_t dim = DeformedGrid::dimension;
+
+public:
+  using Matrix = MatrixType;
+  using Vector = VectorType;
+  using DeformedGrid = DeformedGridTEMPLATE;
+  using Nonlinearity = NonlinearityTEMPLATE;
+
+public:
+  using Step = Dune::Contact::NonSmoothNewtonMGStep<MatrixType, VectorType>;
+  using TransferOperator = Dune::Contact::NonSmoothNewtonContactTransfer<VectorType>;
+
+  SolverFactoryOld(Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
+                const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes);
+
+  std::shared_ptr<Step> getStep();
+
+  const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& getNBodyAssembler() const {
+    return nBodyAssembler_;
+  }
+
+private:
+  const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler_;
+
+  //ProjectedBlockGSStep<MatrixType, VectorType> baseSolverStep_;
+  BlockGSStep<MatrixType, VectorType> baseSolverStep_;
+  EnergyNorm<MatrixType, VectorType> baseEnergyNorm_;
+  LoopSolver<VectorType> baseSolver_;
+
+  ProjectedBlockGSStep<MatrixType, VectorType> presmoother_;
+  ProjectedBlockGSStep<MatrixType, VectorType> postsmoother_;
+  std::shared_ptr<Step> multigridStep_;
+
+  std::vector<std::shared_ptr<TransferOperator>> transferOperators_;
+};
+#endif

src/spatial-solving/solverfactory_tmpl.cc

@@ -5,18 +5,15 @@
 #include "../explicitgrid.hh"
 #include "../explicitvectors.hh"
 
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
+#include "../../dune/tectonic/globalfriction.hh"
+#include "tnnmg/functional.hh"
+#include "tnnmg/zerononlinearity.hh"
+#include "solverfactory.hh"
 
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
+using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
+using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
+using MyZeroFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, double>;
+using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
 
-template class SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
-template class SolverFactory<
-    MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
-                ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-    Grid>;
+template class SolverFactory<MyFunctional, BitVector>;
+template class SolverFactory<MyZeroFunctional, BitVector>;

src/spatial-solving/solverfactory_tmpl_old.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+
+#include <dune/tectonic/globalfriction.hh>
+
+template class SolverFactoryOld<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+
+/*template class SolverFactory<
+    MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
+                ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
+    Grid>;*/

src/spatial-solving/tnnmg/functional.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=8 sw=2 et sts=2:
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_FUNCTIONAL_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_FUNCTIONAL_HH
+
+#include <cstddef>
+#include <type_traits>
+
+#include <dune/solvers/common/wrapownshare.hh>
+#include <dune/solvers/common/copyorreference.hh>
+#include <dune/solvers/common/interval.hh>
+#include <dune/solvers/common/resize.hh>
+
+#include <dune/matrix-vector/algorithm.hh>
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+
+/** \brief Coordinate restriction of box constrained quadratic functional with nonlinearity;
+ *         mainly used for presmoothing in TNNMG algorithm
+ *
+ *  \tparam M global matrix type
+ *  \tparam V global vector type
+ *  \tparam V nonlinearity type
+ *  \tparam R Range type
+ */
+template<class M, class V, class N, class R>
+class ShiftedFunctional : public Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunctional<M,V,R> {
+  using Base = Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunctional<M,V,R>;
+
+public:
+  using Matrix = M;
+  using Vector = V;
+  using Nonlinearity = N;
+  using LowerObstacle = Vector;
+  using UpperObstacle = Vector;
+  using Range = R;
+
+  ShiftedFunctional(const Matrix& quadraticPart, const Vector& linearPart, const Nonlinearity& phi, const LowerObstacle& lowerObstacle, const UpperObstacle& upperObstacle, const Vector& origin) :
+    Base(quadraticPart, linearPart, lowerObstacle, upperObstacle, origin),
+    phi_(phi)
+  {}
+
+  Range operator()(const Vector& v) const
+  {
+    auto temp = Base::Base::origin();
+    temp += v;
+    if (checkInsideBox(temp, this->originalLowerObstacle(), this->originalUpperObstacle()))
+      return Base::Base::operator()(v) + phi_.operator()(temp);
+    return std::numeric_limits<Range>::max();
+  }
+
+  const auto& phi() const {
+      return phi_;
+  }
+
+protected:
+  const Nonlinearity& phi_;
+};
+
+
+
+
+
+/** \brief Coordinate restriction of box constrained quadratic functional with nonlinearity;
+ *         mainly used for line search step in TNNMG algorithm
+ *
+ *  \tparam M Global matrix type
+ *  \tparam V Global vector type
+ *  \tparam N nonlinearity type
+ *  \tparam L Global lower obstacle type
+ *  \tparam U Global upper obstacle type
+ *  \tparam R Range type
+ */
+template<class M, class V, class N, class L, class U, class R=double>
+class DirectionalRestriction :
+  public Dune::TNNMG::BoxConstrainedQuadraticDirectionalRestriction<M,V,L,U,R>
+{
+  using Base = Dune::TNNMG::BoxConstrainedQuadraticDirectionalRestriction<M,V,L,U,R>;
+  using Nonlinearity = N;
+  using Interval = typename Dune::Solvers::Interval<R>;
+
+public:
+  using GlobalMatrix = typename Base::GlobalMatrix;
+  using GlobalVector = typename Base::GlobalVector;
+  using GlobalLowerObstacle = typename Base::GlobalLowerObstacle;
+  using GlobalUpperObstacle = typename Base::GlobalUpperObstacle;
+
+  using Matrix = typename Base::Matrix;
+  using Vector = typename Base::Vector;
+
+  DirectionalRestriction(const GlobalMatrix& matrix, const GlobalVector& linearTerm, const Nonlinearity& phi, const GlobalLowerObstacle& lower, const GlobalLowerObstacle& upper, const GlobalVector& origin, const GlobalVector& direction) :
+    Base(matrix, linearTerm, lower, upper, origin, direction),
+    origin_(origin),
+    direction_(direction),
+    phi_(phi)
+  {
+    phi_.directionalDomain(origin_, direction_, domain_);
+
+    std::cout << domain_[0] << " " << domain_[1] << "Phi domain:" << std::endl;
+    std::cout << this->defectLower_ << " " << this->defectUpper_ << "defect obstacles:" << std::endl;
+    if (domain_[0] < this->defectLower_) {
+       domain_[0] = this->defectLower_;
+    }
+
+    if (domain_[1] > this->defectUpper_) {
+       domain_[1] = this->defectUpper_;
+    }
+  }
+
+    auto subDifferential(double x) const {
+      Interval D;
+      GlobalVector uxv = origin_;
+
+      Dune::MatrixVector::addProduct(uxv, x, direction_);
+      phi_.directionalSubDiff(uxv, direction_, D);
+      auto const Axmb = this->quadraticPart_ * x - this->linearPart_;
+      D[0] += Axmb;
+      D[1] += Axmb;
+
+      return D;
+    }
+
+    auto domain() const {
+      return domain_;
+    }
+
+    const GlobalVector& origin() const {
+        return origin_;
+    }
+
+    const GlobalVector& direction() const {
+        return direction_;
+    }
+
+protected:
+  const GlobalVector& origin_;
+  const GlobalVector& direction_;
+
+  const Nonlinearity& phi_;
+  Interval domain_;
+};
+
+/** \brief Box constrained quadratic functional with nonlinearity
+ *         Note: call setIgnore() to set up functional in affine subspace with ignore information
+ *
+ *  \tparam V Vector type
+ *  \tparam N Nonlinearity type
+ *  \tparam L Lower obstacle type
+ *  \tparam U Upper obstacle type
+ *  \tparam R Range type
+ */
+template<class V, class N, class L, class U, class O, class D, class R>
+class FirstOrderModelFunctional {
+    using Interval = typename Dune::Solvers::Interval<R>;
+
+public:
+    using Vector = std::decay_t<V>;
+    using Nonlinearity = std::decay_t<N>;
+    using LowerObstacle = std::decay_t<L>;
+    using UpperObstacle = std::decay_t<U>;
+    using Range = R;
+
+    using field_type = typename Vector::field_type;
+
+private:
+    void init() {
+        auto origin = origin_.get();
+        auto direction = direction_.get();
+
+        // set defect obstacles
+        defectLower_ = -std::numeric_limits<field_type>::max();
+        defectUpper_ = std::numeric_limits<field_type>::max();
+        Dune::TNNMG::directionalObstacles(origin, direction, lower_.get(), upper_.get(), defectLower_, defectUpper_);
+
+        // set domain
+        phi_.get().directionalDomain(origin, direction, domain_);
+
+        if (domain_[0] < defectLower_) {
+           domain_[0] = defectLower_;
+        }
+
+        if (domain_[1] > defectUpper_) {
+           domain_[1] = defectUpper_;
+        }
+
+        // set quadratic and linear parts of functional
+        quadraticPart_ = direction*direction;
+        quadraticPart_ *= maxEigenvalue_;
+
+        linearPart_ = linearTerm_.get()*direction - maxEigenvalue_*(direction*origin);
+    }
+
+public:
+    template <class MM, class VV1, class NN, class LL, class UU, class VV2, class VV3>
+    FirstOrderModelFunctional(
+            const MM& matrix,
+            VV1&& linearTerm,
+            NN&& phi,
+            LL&& lower,
+            UU&& upper,
+            VV2&& origin,
+            VV3&& direction) :
+        linearTerm_(std::forward<VV1>(linearTerm)),
+        lower_(std::forward<LL>(lower)),
+        upper_(std::forward<UU>(upper)),
+        origin_(std::forward<VV2>(origin)),
+        direction_(std::forward<VV3>(direction)),
+        phi_(std::forward<NN>(phi)) {
+
+        // set maxEigenvalue_ from matrix
+        Vector eigenvalues;
+        Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
+        maxEigenvalue_ = *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+
+        init();
+    }
+
+    template <class VV1, class NN, class LL, class UU, class VV2, class VV3>
+    FirstOrderModelFunctional(
+            const Range& maxEigenvalue,
+            VV1&& linearTerm,
+            NN&& phi,
+            LL&& lower,
+            UU&& upper,
+            VV2&& origin,
+            VV3&& direction) :
+        linearTerm_(std::forward<VV1>(linearTerm)),
+        lower_(std::forward<LL>(lower)),
+        upper_(std::forward<UU>(upper)),
+        origin_(std::forward<VV2>(origin)),
+        direction_(std::forward<VV3>(direction)),
+        phi_(std::forward<NN>(phi)),
+        maxEigenvalue_(maxEigenvalue) {
+
+        init();
+    }
+
+  template <class BitVector>
+  auto getIgnoreFunctional(const BitVector& ignore) const {
+      Vector direction = direction_.get();
+      Vector origin = origin_.get();
+     // Dune::Solvers::resizeInitializeZero(direction, linearPart());
+      for (size_t j = 0; j < ignore.size(); ++j)
+        if (ignore[j])
+          direction[j] = 0; // makes sure result remains in subspace after correction
+        else
+          origin[j] = 0; // shift affine offset
+
+      return FirstOrderModelFunctional<Vector, Nonlinearity&, LowerObstacle, UpperObstacle, Vector, Vector, Range>(maxEigenvalue_, linearTerm_, phi_, lower_, upper_, std::move(origin), std::move(direction));
+  }
+
+  Range operator()(const Vector& v) const
+  {
+    DUNE_THROW(Dune::NotImplemented, "Evaluation of FirstOrderModelFunctional not implemented");
+  }
+
+  auto subDifferential(double x) const {
+    Interval Di;
+    Vector uxv = origin_.get();
+
+    Dune::MatrixVector::addProduct(uxv, x, direction_.get());
+    phi_.get().directionalSubDiff(uxv, direction_.get(), Di);
+
+    const auto Axmb = quadraticPart_ * x - linearPart_;
+    Di[0] += Axmb;
+    Di[1] += Axmb;
+
+    return Di;
+  }
+
+  const Interval& domain() const {
+    return domain_;
+  }
+
+  const Vector& origin() const {
+      return origin_.get();
+  }
+
+  const Vector& direction() const {
+      return direction_.get();
+  }
+
+  const field_type& defectLower() const {
+      return defectLower_;
+  }
+
+  const field_type& defectUpper() const {
+      return defectUpper_;
+  }
+
+private:
+  Dune::Solvers::ConstCopyOrReference<V> linearTerm_;
+  Dune::Solvers::ConstCopyOrReference<L> lower_;
+  Dune::Solvers::ConstCopyOrReference<U> upper_;
+
+  Dune::Solvers::ConstCopyOrReference<O> origin_;
+  Dune::Solvers::ConstCopyOrReference<D> direction_;
+
+  Dune::Solvers::ConstCopyOrReference<N> phi_;
+
+  Interval domain_;
+
+  Range maxEigenvalue_;
+
+  field_type quadraticPart_;
+  field_type linearPart_;
+
+  field_type defectLower_;
+  field_type defectUpper_;
+};
+
+// \ToDo This should be an inline friend of ShiftedBoxConstrainedQuadraticFunctional
+// but gcc-4.9.2 shipped with debian jessie does not accept
+// inline friends with auto return type due to bug-59766.
+// Notice, that this is fixed in gcc-4.9.3.
+template<class Index, class M, class V, class Nonlinearity, class R>
+auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, const Index& i)
+{
+  using Range = R;
+  using LocalMatrix = std::decay_t<decltype(f.quadraticPart()[i][i])>;
+  using LocalVector = std::decay_t<decltype(f.originalLinearPart()[i])>;
+  using LocalLowerObstacle = std::decay_t<decltype(f.originalLowerObstacle()[i])>;
+  using LocalUpperObstacle = std::decay_t<decltype(f.originalUpperObstacle()[i])>;
+
+  using namespace Dune::MatrixVector;
+  namespace H = Dune::Hybrid;
+
+  const LocalMatrix* Aii_p = nullptr;
+
+  LocalVector ri = f.originalLinearPart()[i];
+  const auto& Ai = f.quadraticPart()[i];
+  sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
+      // TODO Here we must implement a wrapper to guarantee that this will work with proxy matrices!
+      H::ifElse(H::equals(j, i), [&](auto&& id){
+        Aii_p = id(&Aij);
+      });
+      Dune::TNNMG::Imp::mmv(Aij, f.origin()[j], ri, Dune::PriorityTag<1>());
+  });
+
+  LocalLowerObstacle dli = f.originalLowerObstacle()[i];
+  LocalUpperObstacle dui = f.originalUpperObstacle()[i];
+  dli -= f.origin()[i];
+  dui -= f.origin()[i];
+
+  auto&& phii = f.phi().restriction(i);
+
+  auto v = ri;
+  double const vnorm = v.two_norm();
+  if (vnorm > 1.0)
+        v /= vnorm;
+
+  return FirstOrderModelFunctional<LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, LocalVector, Range>(*Aii_p, std::move(ri), std::move(phii), std::move(dli), std::move(dui), std::move(f.origin()[i]), std::move(v));
+}
+
+
+/** \brief Box constrained quadratic functional with nonlinearity
+ *
+ *  \tparam M Matrix type
+ *  \tparam V Vector type
+ *  \tparam L Lower obstacle type
+ *  \tparam U Upper obstacle type
+ *  \tparam R Range type
+ */
+template<class M, class V, class N, class L, class U, class R>
+class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L, U, R>
+{
+private:
+  using Base = Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L, U, R>;
+
+
+public:
+  using Nonlinearity = std::decay_t<N>;
+
+  using Matrix = typename Base::Matrix;
+  using Vector = typename Base::Vector;
+  using Range = typename Base::Range;
+  using LowerObstacle = typename Base::LowerObstacle;
+  using UpperObstacle = typename Base::UpperObstacle;
+
+private:
+    Dune::Solvers::ConstCopyOrReference<N> phi_;
+
+public:
+
+  template <class MM, class VV, class NN, class LL, class UU>
+  Functional(
+          MM&& matrix,
+          VV&& linearPart,
+          NN&& phi,
+          LL& lower,
+          UU& upper) :
+    Base(std::forward<MM>(matrix), std::forward<VV>(linearPart), std::forward<LL>(lower), std::forward<UU>(upper)),
+    phi_(std::forward<NN>(phi))
+  {}
+
+  const Nonlinearity& phi() const {
+    return phi_.get();
+  }
+
+  Range operator()(const Vector& v) const
+  { 
+    if (Dune::TNNMG::checkInsideBox(v, this->lower_.get(), this->upper_.get()))
+      return Base::Base::operator()(v) + phi_.get().operator()(v);
+    return std::numeric_limits<Range>::max();
+  }
+
+  friend auto directionalRestriction(const Functional& f, const Vector& origin, const Vector& direction)
+    -> DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>
+  {
+    return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin, direction);
+  }
+
+  friend auto shift(const Functional& f, const Vector& origin)
+    -> ShiftedFunctional<Matrix, Vector, Nonlinearity, Range>
+  {
+    return ShiftedFunctional<Matrix, Vector, Nonlinearity, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
+  }
+};
+
+
+#endif

src/spatial-solving/tnnmg/levelglobalpreconditioner.hh

+#ifndef LEVEL_GLOBAL_PRECONDITIONER_HH
+#define LEVEL_GLOBAL_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+
+#include <dune/istl/umfpack.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+
+#include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
+#include <dune/faultnetworks/levelinterfacenetwork.hh>
+#include <dune/faultnetworks/utils/debugutils.hh>
+
+
+template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+class LevelGlobalPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
+
+public:
+    enum BoundaryMode {homogeneous, fromIterate};
+
+private:
+    typedef typename BasisType::GridView GridView;
+    typedef typename GridView::Grid GridType;
+
+    const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
+    const LocalAssembler& localAssembler_;
+    const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
+
+    const GridType& grid_;
+    const int level_;
+    const BasisType basis_;
+
+    Dune::BitSetVector<1> boundaryDofs_;
+    MatrixType matrix_;
+    VectorType rhs_;
+
+    bool requireBuild_;
+
+public:
+
+    // for each active fault in levelInterfaceNetwork: set up local problem, compute local correction
+    LevelGlobalPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+                             const LocalAssembler& localAssembler,
+                             const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers) :
+          levelInterfaceNetwork_(levelInterfaceNetwork),
+          localAssembler_(localAssembler),
+          localInterfaceAssemblers_(localInterfaceAssemblers),
+          grid_(levelInterfaceNetwork_.grid()),
+          level_(levelInterfaceNetwork_.level()),
+          basis_(levelInterfaceNetwork_)
+    {
+        assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
+    }
+
+    void build() {
+        //printBasisDofLocation(basis_);
+
+        GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
+        globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
+
+        typedef typename MatrixType::row_type RowType;
+        typedef typename RowType::ConstIterator ColumnIterator;
+
+        const GridView& gridView = levelInterfaceNetwork_.levelGridView();
+
+        // set boundary conditions
+        BoundaryPatch<GridView> boundaryPatch(gridView, true);
+        constructBoundaryDofs(boundaryPatch, basis_, boundaryDofs_);
+
+        for(size_t i=0; i<boundaryDofs_.size(); i++) {
+            if(!boundaryDofs_[i][0])
+                continue;
+
+            RowType& row = matrix_[i];
+
+            ColumnIterator cIt    = row.begin();
+            ColumnIterator cEndIt = row.end();
+
+            for(; cIt!=cEndIt; ++cIt) {
+                row[cIt.index()] = 0;
+            }
+
+            row[i] = 1;
+        }
+
+        requireBuild_ = false;
+    }
+
+
+
+    virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
+        this->x_ = &x;
+        updateRhs(rhs);
+        this->mat_ = Dune::stackobject_to_shared_ptr(mat); // mat will be set but not used
+    }
+
+    virtual void iterate() {      
+        //print(matrix_, "coarse matrix: ");
+        //print(rhs_, "coarse rhs: ");
+
+        // compute solution directly
+
+        if (requireBuild_)
+            DUNE_THROW(Dune::Exception, "LevelGlobalPreconditioner::iterate() Call build() before solving the global problem!");
+
+        Dune::Timer timer;
+        timer.reset();
+        timer.start();
+
+        this->x_->resize(matrix_.M());
+        *(this->x_) = 0;
+
+        #if HAVE_UMFPACK
+        Dune::InverseOperatorResult res;
+        VectorType rhsCopy(rhs_);
+
+        Dune::UMFPack<MatrixType> solver(matrix_);
+        solver.apply(*(this->x_), rhsCopy, res);
+
+        #else
+        #error No UMFPack!
+        #endif
+
+        //std::cout << "LevelGlobalPreconditioner::iterate() Solving global problem took: " << timer.elapsed() << " seconds" << std::endl;
+
+    }
+
+    const BasisType& basis() const {
+        return basis_;
+    }
+
+    const GridView& gridView() const {
+        return basis_.getGridView();
+    }
+
+    const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
+        return levelInterfaceNetwork_;
+    }
+
+private:
+    void updateRhs(const VectorType& rhs){
+        rhs_.resize(matrix_.M());
+        rhs_ = 0;
+
+        for (size_t i=0; i<rhs.size(); i++) {
+            if (!boundaryDofs_[i][0])
+                rhs_[i] = rhs[i];
+        }
+    }
+};
+
+#endif
+

src/spatial-solving/tnnmg/levelpatchpreconditioner.hh

+#ifndef LEVEL_PATCH_PRECONDITIONER_HH
+#define LEVEL_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/common/timer.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+
+#include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
+#include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
+#include <dune/faultnetworks/localproblem.hh>
+#include <dune/faultnetworks/levelinterfacenetwork.hh>
+#include <dune/faultnetworks/utils/debugutils.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+
+
+template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
+
+public:
+    enum Mode {additive, multiplicative};
+    enum BoundaryMode {homogeneous, fromIterate};
+
+private:
+    typedef typename BasisType::GridView GridView;
+    typedef typename GridView::Grid GridType;
+
+    const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
+    const BasisType& patchLevelBasis_;
+    const LocalAssembler& localAssembler_;
+    const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
+    const Mode mode_;
+
+    const GridType& grid_;
+    const int level_;
+    const BasisType basis_;
+
+    size_t patchDepth_;
+    BoundaryMode boundaryMode_;
+
+    MatrixType matrix_;
+    std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
+
+public:
+
+    // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
+    LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+                             const BasisType& patchLevelBasis,
+                             const LocalAssembler& localAssembler,
+                             const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+                             const Mode mode = LevelPatchPreconditioner::Mode::additive) :
+          levelInterfaceNetwork_(levelInterfaceNetwork),
+          patchLevelBasis_(patchLevelBasis),
+          localAssembler_(localAssembler),
+          localInterfaceAssemblers_(localInterfaceAssemblers),
+          mode_(mode),
+          grid_(levelInterfaceNetwork_.grid()),
+          level_(levelInterfaceNetwork_.level()),
+          basis_(levelInterfaceNetwork_)
+    {
+        setPatchDepth();
+        setBoundaryMode();
+
+        assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
+    }
+
+    void build() {
+        // assemble stiffness matrix for entire level including all faults
+        GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
+        globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
+
+        // set boundary conditions
+        Dune::BitSetVector<1> globalBoundaryDofs;
+        BoundaryPatch<GridView> boundaryPatch(levelInterfaceNetwork_.levelGridView(), true);
+        constructBoundaryDofs(boundaryPatch, basis_, globalBoundaryDofs);
+
+        typedef typename MatrixType::row_type RowType;
+        typedef typename RowType::ConstIterator ColumnIterator;
+
+        for(size_t i=0; i<globalBoundaryDofs.size(); i++) {
+            if(!globalBoundaryDofs[i][0])
+                continue;
+
+            RowType& row = matrix_[i];
+
+            ColumnIterator cIt    = row.begin();
+            ColumnIterator cEndIt = row.end();
+
+            for(; cIt!=cEndIt; ++cIt) {
+                row[cIt.index()] = 0;
+            }
+            row[i] = 1;
+        }
+
+        // init vertexInElements
+        const int dim = GridType::dimension;
+        typedef typename GridType::template Codim<0>::Entity EntityType;
+        std::vector<std::vector<EntityType>> vertexInElements;
+
+        const GridView& patchLevelGridView = patchLevelBasis_.getGridView();
+        vertexInElements.resize(patchLevelGridView.size(dim));
+        for (size_t i=0; i<vertexInElements.size(); i++) {
+            vertexInElements[i].resize(0);
+        }
+
+        typedef typename BasisType::LocalFiniteElement FE;
+        typedef  typename GridView::template Codim <0>::Iterator  ElementLevelIterator;
+        ElementLevelIterator endElemIt = patchLevelGridView.template end <0>();
+        for (ElementLevelIterator  it = patchLevelGridView. template begin <0>(); it!=endElemIt; ++it) {
+
+            // compute coarseGrid vertexInElements
+            const FE& coarseFE = patchLevelBasis_.getLocalFiniteElement(*it);
+
+            for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
+                int dofIndex = patchLevelBasis_.indexInGridView(*it, i);
+                vertexInElements[dofIndex].push_back(*it);
+            }
+        }
+
+        localProblems_.resize(vertexInElements.size());
+
+        std::cout << std::endl;
+        std::cout << "---------------------------------------------" << std::endl;
+        std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;
+
+        // init local fine level corrections
+        Dune::Timer timer;
+        timer.reset();
+        timer.start();
+
+        const int patchLevel = patchLevelBasis_.faultNetwork().level();
+        for (size_t i=0; i<vertexInElements.size(); i++) {
+            //std::cout << i << std::endl;
+            //std::cout << "---------------" << std::endl;
+
+            SupportPatchFactory<BasisType> patchFactory(patchLevelBasis_, basis_, patchLevel, level_, vertexInElements, i, patchDepth_);
+            std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
+            Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
+            VectorType rhs;
+            rhs.resize(matrix_.N());
+            rhs = 0;
+
+            //print(localToGlobal, "localToGlobal: ");
+            //print(boundaryDofs, "boundaryDofs: ");
+
+            localProblems_[i] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs, localToGlobal, boundaryDofs);
+
+            /*
+            if ((i+1) % 10 == 0) {
+                std::cout << '\r' << std::floor((i+1.0)/patchLevelBasis_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*patchLevelBasis_.size() << " seconds." << std::flush;
+            }*/
+        }
+
+        std::cout << std::endl;
+        std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
+        std::cout << "---------------------------------------------" << std::endl;
+        timer.stop();
+    }
+
+    void setPatchDepth(const size_t patchDepth = 0) {
+        patchDepth_ = patchDepth;
+    }
+
+    void setBoundaryMode(const BoundaryMode boundaryMode = LevelPatchPreconditioner::BoundaryMode::homogeneous) {
+        boundaryMode_ = boundaryMode;
+    }
+
+    virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
+        this->x_ = &x;
+        this->rhs_ = &rhs;
+        this->mat_ = Dune::stackobject_to_shared_ptr(mat);
+
+        for (size_t i=0; i<localProblems_.size(); i++) {
+            if (boundaryMode_ == BoundaryMode::homogeneous)
+                localProblems_[i]->updateRhs(rhs);
+            else
+                localProblems_[i]->updateRhsAndBoundary(rhs, x);
+        }
+    }
+
+    virtual void iterate() {
+        if (mode_ == additive)
+            iterateAdd();
+        else
+            iterateMult();
+    }
+
+    void iterateAdd() {
+        *(this->x_) = 0;	
+
+        VectorType it, x;
+        for (size_t i=0; i<localProblems_.size(); i++) {
+            localProblems_[i]->solve(it);
+            localProblems_[i]->prolong(it, x);
+
+            /*if (i==5) {
+                writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
+            }*/
+
+            *(this->x_) += x;
+        }
+    }
+
+    void iterateMult() {
+        *(this->x_) = 0;
+
+        VectorType it, x;
+        for (size_t i=0; i<localProblems_.size(); i++) {
+            VectorType updatedRhs(*(this->rhs_));
+            matrix_.mmv(*(this->x_), updatedRhs);
+
+            //step(i);
+            //print(updatedRhs, "localRhs: ");
+            //writeToVTK(basis_, updatedRhs, "/storage/mi/podlesny/data/faultnetworks/rhs/local", "exactvertexdata_step_"+std::to_string(i));
+
+            if (boundaryMode_ == BoundaryMode::homogeneous)
+                localProblems_[i]->updateRhs(updatedRhs);
+            else
+                localProblems_[i]->updateRhsAndBoundary(updatedRhs, *(this->x_));
+
+            localProblems_[i]->solve(it);
+            localProblems_[i]->prolong(it, x);
+
+
+            //print(it, "it: ");
+            //print(x, "x: ");
+
+            //writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/sol/local", "exactvertexdata_step_"+std::to_string(i));
+
+            /*if (i==5) {
+                writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
+            }*/
+
+            *(this->x_) += x;
+        }
+    }
+
+    const BasisType& basis() const {
+        return basis_;
+    }
+
+    const GridView& gridView() const {
+        return basis_.getGridView();
+    }
+
+    const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
+        return levelInterfaceNetwork_;
+    }
+
+    size_t size() const {
+        return localProblems_.size();
+    }
+};
+
+#endif
+

src/spatial-solving/tnnmg/linearization.hh

+#ifndef DUNE_TECTONIC_LINEARIZATION_HH
+#define DUNE_TECTONIC_LINEARIZATION_HH
+
+#include <cstddef>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/hybridutilities.hh>
+
+#include <dune/matrix-vector/algorithm.hh>
+
+#include <dune/istl/matrixindexset.hh>
+
+//#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
+
+#include "../../utils/debugutils.hh"
+
+/**
+ * derived from Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<F,BV>
+ */
+template<class F, class BV>
+class Linearization {
+public:
+  using Matrix = typename F::Matrix;
+  using Vector = typename F::Vector;
+  using BitVector = BV;
+  using ConstrainedMatrix = Matrix;
+  using ConstrainedVector = Vector;
+  using ConstrainedBitVector = BitVector;
+
+private:
+  using This = Linearization<F,BV>;
+
+  template <class T>
+  void determineRegularityTruncation(const Vector& x, BitVector& truncationFlags, const T& truncationTolerance)
+  {
+    for (size_t i = 0; i < x.size(); ++i) {
+      if (f_.phi().regularity(i, x[i]) > truncationTolerance) {
+        truncationFlags[i] = true;
+      }
+    }
+  }
+
+  template<class NV, class NBV, class T>
+  static void determineTruncation(const NV& x, const NV& lower, const NV& upper, NBV&& truncationFlags, const T& truncationTolerance)
+  {
+    namespace H = Dune::Hybrid;
+    H::ifElse(Dune::IsNumber<NV>(), [&](auto id){
+      if ((id(x) <= id(lower)+truncationTolerance) || (id(x) >= id(upper) - truncationTolerance))
+        id(truncationFlags) = true;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(x))), [&](auto&& i) {
+        This::determineTruncation(x[i], lower[i], upper[i], truncationFlags[i], truncationTolerance);
+      });
+    });
+  }
+
+  template<class NV, class NBV>
+  static void truncateVector(NV& x, const NBV& truncationFlags)
+  {
+    namespace H = Dune::Hybrid;
+    H::ifElse(Dune::IsNumber<NV>(), [&](auto id){
+      if (id(truncationFlags))
+        id(x) = 0;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(x))), [&](auto&& i) {
+        This::truncateVector(x[i], truncationFlags[i]);
+      });
+    });
+  }
+
+  template<class NM, class RBV, class CBV>
+  static void truncateMatrix(NM& A, const RBV& rowTruncationFlags, const CBV& colTruncationFlags)
+  {
+    namespace H = Dune::Hybrid;
+    using namespace Dune::MatrixVector;
+    H::ifElse(Dune::IsNumber<NM>(), [&](auto id){
+      if(id(rowTruncationFlags) or id(colTruncationFlags))
+        A = 0;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(rowTruncationFlags))), [&](auto&& i) {
+        auto&& Ai = A[i];
+        sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
+            This::truncateMatrix(Aij, rowTruncationFlags[i], colTruncationFlags[j]);
+        });
+      });
+    });
+  }
+
+public:
+  Linearization(const F& f, const BitVector& ignore) :
+      f_(f),
+      ignore_(ignore),
+      truncationTolerance_(1e-10)
+  {}
+
+  void bind(const Vector& x) {
+      //std::cout << "Linearization::bind()" << std::endl;
+
+      auto&& A = f_.quadraticPart();
+      auto&& phi = f_.phi();
+
+      // determine which components to truncate
+      // ----------------
+      truncationFlags_ = ignore_;
+
+      determineTruncation(x, f_.lowerObstacle(), f_.upperObstacle(), truncationFlags_, truncationTolerance_); // obstacle truncation
+      determineRegularityTruncation(x, truncationFlags_, truncationTolerance_); // truncation due to regularity deficit of nonlinearity
+
+
+      // compute hessian
+      // ---------------
+
+      // construct sparsity pattern for linearization
+      Dune::MatrixIndexSet indices(A.N(), A.M());
+      indices.import(A);
+      phi.addHessianIndices(indices);
+
+      // construct matrix from pattern and initialize it
+      indices.exportIdx(hessian_);
+      hessian_ = 0.0;
+
+      // quadratic part
+      for (size_t i = 0; i < A.N(); ++i) {
+        auto const end = std::end(A[i]);
+        for (auto it = std::begin(A[i]); it != end; ++it)
+          hessian_[i][it.index()] += *it;
+      }
+
+      //print(hessian_, "Hessian:");
+      //std::cout << "--------" << (double) hessian_[21][21] << "------------" << std::endl;
+
+      // nonlinearity part
+      phi.addHessian(x, hessian_);
+
+      // compute gradient
+      // ----------------
+
+      // quadratic part
+      negativeGradient_ = derivative(f_)(x); // A*x - b
+
+      // nonlinearity part
+      phi.addGradient(x, negativeGradient_);
+
+      // -grad is needed for Newton step
+      negativeGradient_ *= -1;
+
+
+
+      // truncate gradient and hessian
+      truncateVector(negativeGradient_, truncationFlags_);
+      truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
+
+      //print(hessian_, "Hessian:");
+  }
+
+
+  void extendCorrection(ConstrainedVector& cv, Vector& v) const {
+    v = cv;
+    truncateVector(v, truncationFlags_);
+  }
+
+  const BitVector& truncated() const {
+    return truncationFlags_;
+  }
+
+  const auto& negativeGradient() const {
+    return negativeGradient_;
+  }
+
+  const auto& hessian() const {
+    return hessian_;
+  }
+
+private:
+  const F& f_;
+  const BitVector& ignore_;
+
+  double truncationTolerance_;
+
+  Vector negativeGradient_;
+  Matrix hessian_;
+  BitVector truncationFlags_;
+};
+
+#endif

src/spatial-solving/tnnmg/linesearchsolver.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_LINESEARCHSOLVER_HH
+#define DUNE_TECTONIC_LINESEARCHSOLVER_HH
+
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include "../../utils/almostequal.hh"
+
+/**
+ * \brief A local solver for scalar quadratic obstacle problems with nonlinearity
+ * using bisection
+ */
+class LineSearchSolver
+{
+public:
+  template<class Vector, class Functional, class BitVector>
+  void operator()(Vector& x, const Functional& f, const BitVector& ignore) const {
+    x = 0;
+
+    Dune::Solvers::Interval<double> D;
+    D = f.subDifferential(0);
+    if (D[1] > 0) // NOTE: Numerical instability can actually get us here
+        return;
+
+    if (almost_equal(f.domain()[0], f.domain()[1], 2)) {
+        x = f.domain()[0];
+        return;
+    }
+
+    if (not ignore) {
+        int bisectionsteps = 0;
+        const Bisection globalBisection;
+       // std::cout << "LineSearchSolver: starting bisection" << std::endl;
+        x = globalBisection.minimize(f, 1.0, 0.0, bisectionsteps);
+
+        x = f.domain().projectIn(x);
+
+
+        //std::cout << "LineSearchSolver: bisection terminated" << std::endl;
+        // x = globalBisection.minimize(f, vNorm, 0.0, bisectionsteps) / vNorm;  // used rescaled v in f?
+    }
+  }
+};
+
+#endif
+