diff --git a/dune/tectonic/spatial-solving/functionalfactory.hh b/dune/tectonic/spatial-solving/functionalfactory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c96557a53142d37612e715854d1ca72f92fbf7fd
--- /dev/null
+++ b/dune/tectonic/spatial-solving/functionalfactory.hh
@@ -0,0 +1,102 @@
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_GLOBALPROBLEM_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_GLOBALPROBLEM_HH
+
+#include <vector>
+#include <memory>
+
+#include "tnnmg/functional.hh"
+
+template <class NBodyAssembler, class Nonlinearity, class MatrixType, class VectorType>
+class FunctionalFactory {
+ using Vector = VectorType;
+ using Matrix = MatrixType;
+
+ const static int dims = Vector::block_type::dimension;
+
+public:
+ using Functional = Functional<Matrix&, Vector&, Nonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
+
+ FunctionalFactory(const NBodyAssembler& nBodyAssembler, const Nonlinearity& nonlinearity, const std::vector<Matrix>& matrices, const std::vector<Vector>& rhs) :
+ nBodyAssembler_(nBodyAssembler),
+ nonlinearity_(nonlinearity),
+ rhs_(rhs) {
+
+ const auto nBodies = nBodyAssembler_.nGrids();
+
+ matrices_.resize(nBodies);
+ for (int i=0; i<nBodies; i++) {
+ matrices_[i] = &matrices[i];
+ }
+ }
+
+ FunctionalFactory(const NBodyAssembler& nBodyAssembler, const Nonlinearity& nonlinearity, const std::vector<std::shared_ptr<Matrix>>& matrices, const std::vector<Vector>& rhs) :
+ nBodyAssembler_(nBodyAssembler),
+ nonlinearity_(nonlinearity),
+ rhs_(rhs) {
+
+ const auto nBodies = nBodyAssembler_.nGrids();
+
+ matrices_.resize(nBodies);
+ for (size_t i=0; i<nBodies; i++) {
+ matrices_[i] = matrices[i].get();
+ }
+ }
+
+ void build() {
+ // assemble full global contact problem
+ nBodyAssembler_.assembleJacobian(matrices_, bilinearForm_);
+ nBodyAssembler_.assembleRightHandSide(rhs_, totalRhs_);
+
+ // get lower and upper obstacles
+ const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
+ lower_.resize(totalObstacles.size());
+ upper_.resize(totalObstacles.size());
+
+ for (size_t j=0; j<totalObstacles.size(); ++j) {
+ const auto& totalObstaclesj = totalObstacles[j];
+ auto& lowerj = lower_[j];
+ auto& upperj = upper_[j];
+ for (size_t d=0; d<dims; ++d) {
+ lowerj[d] = totalObstaclesj[d][0];
+ upperj[d] = totalObstaclesj[d][1];
+ }
+ }
+
+ // set up functional
+ functional_ = std::make_shared<Functional>(bilinearForm_, totalRhs_, nonlinearity_, lower_, upper_);
+ }
+
+ std::shared_ptr<Functional> functional() const {
+ return functional_;
+ }
+
+ const Matrix& bilinearForm() const {
+ return bilinearForm_;
+ }
+
+ const Vector& rhs() const {
+ return totalRhs_;
+ }
+
+ const Vector& lower() const {
+ return lower_;
+ }
+
+ const Vector& upper() const {
+ return upper_;
+ }
+
+private:
+ const NBodyAssembler nBodyAssembler_;
+ const Nonlinearity& nonlinearity_;
+ std::vector<const Matrix*> matrices_;
+ const std::vector<Vector>& rhs_;
+
+ Matrix bilinearForm_;
+ Vector totalRhs_;
+ Vector lower_;
+ Vector upper_;
+
+ std::shared_ptr<Functional> functional_;
+};
+#endif
diff --git a/dune/tectonic/spatial-solving/makelinearsolver.hh b/dune/tectonic/spatial-solving/makelinearsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ada0083548dbc7f95c8b3f69092d72fd301f8ae6
--- /dev/null
+++ b/dune/tectonic/spatial-solving/makelinearsolver.hh
@@ -0,0 +1,86 @@
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH
+
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "preconditioners/multilevelpatchpreconditioner.hh"
+
+template <class ContactNetwork, class VectorType>
+auto makeLinearSolver(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) {
+ // make linear solver for linear correction in TNNMGStep
+ using Vector = VectorType;
+ using Matrix = typename ContactNetwork::Matrix;
+
+ /* old, pre multi threading, was unused
+ using Norm = EnergyNorm<Matrix, Vector>;
+ using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+ using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+ const auto& preconditionerParset = parset_.sub("solver.tnnmg.linear.preconditioner");
+
+ Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
+ Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
+ preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
+ preconditioner.build();
+
+ auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
+ cgStep->setPreconditioner(preconditioner);
+
+ Norm norm(*cgStep);
+
+ return std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.linear.tolerance"), norm, Solver::QUIET);
+ */
+
+ // patch preconditioner only needs to be computed once per advance()
+ // make linear solver for linear correction in TNNMGStep
+ using Norm = EnergyNorm<Matrix, Vector>;
+ using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+ using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+ /*const auto& preconditionerParset = parset_.sub("solver.tnnmg.preconditioner");
+
+ Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
+ Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
+ preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
+ preconditioner.build();
+
+ auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
+ cgStep->setPreconditioner(preconditioner);
+
+ Norm norm(*cgStep);
+
+ auto linearSolver = std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+ */
+
+ // set multigrid solver
+ auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
+
+ // transfer operators need to be recomputed on change due to setDeformation()
+ using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+ using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+ TransferOperators transfer(contactNetwork.nLevels()-1);
+ for (size_t i=0; i<transfer.size(); i++) {
+ transfer[i] = std::make_shared<TransferOperator>();
+ transfer[i]->setup(contactNetwork, i, i+1);
+ }
+
+ // Remove any recompute filed so that initially the full transferoperator is assembled
+ for (size_t i=0; i<transfer.size(); i++)
+ std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+ auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+ linearMultigridStep->setMGType(1, 3, 3);
+ linearMultigridStep->setSmoother(smoother);
+ linearMultigridStep->setTransferOperators(transfer);
+
+ Norm norm(*linearMultigridStep);
+
+ return std::make_shared<LinearSolver>(linearMultigridStep, parset.template get<int>("solver.tnnmg.main.multi"), parset.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+ }
+#endif
diff --git a/dune/tectonic/spatial-solving/nonlinearsolver.hh b/dune/tectonic/spatial-solving/nonlinearsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..81245b0f7654b796f7b8f13a349a4fb438d9fc7b
--- /dev/null
+++ b/dune/tectonic/spatial-solving/nonlinearsolver.hh
@@ -0,0 +1,107 @@
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_SOLVENONLINEARPROBLEM_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_SOLVENONLINEARPROBLEM_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/contact/common/dualbasisadapter.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/functions/gridfunctions/gridfunction.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
+#include "../data-structures/enums.hh"
+#include "../data-structures/enumparser.hh"
+
+
+
+#include "../utils/tobool.hh"
+#include "../utils/debugutils.hh"
+
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include "../spatial-solving/contact/nbodycontacttransfer.hh"
+
+#include "solverfactory.hh"
+#include "solverfactory.cc"
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+
+#include "fixedpointiterator.hh"
+
+
+template <class Functional, class BitVector>
+class NonlinearSolver {
+protected:
+ using Matrix = typename Functional::Matrix;
+ using Vector = typename Functional::Vector;
+ using Factory = SolverFactory<Functional, BitVector>;
+
+ using Norm = EnergyNorm<Matrix, Vector>;
+ using SolverType = Dune::Solvers::LoopSolver<Vector>;
+
+ const static int dims = Vector::block_type::dimension;
+
+public:
+ template <class LinearSolver>
+ NonlinearSolver(const Dune::ParameterTree& tnnmgParset, std::shared_ptr<LinearSolver> linearSolver, std::shared_ptr<Functional> functional, const BitVector& dirichletNodes) :
+ functional_(functional),
+ norm_(functional_->quadraticPart()) {
+
+ // set up TNMMG step
+ solverFactory_ = std::make_shared<Factory>(tnnmgParset, functional_, dirichletNodes);
+ solverFactory_->build(linearSolver);
+ }
+
+ auto solve(const Dune::ParameterTree& solverParset, Vector& x) {
+ auto tnnmgStep = solverFactory_->step();
+
+ SolverType solver(*tnnmgStep.get(), solverParset.get<size_t>("maximumIterations"),
+ solverParset.get<double>("tolerance"), norm_,
+ solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+
+ const auto& lower = functional_->lowerObstacle();
+ const auto& upper = functional_->upperObstacle();
+
+ // project in onto admissible set
+ for (size_t i=0; i<x.size(); i++) {
+ for (size_t j=0; j<dims; j++) {
+ if (x[i][j] < lower[i][j]) {
+ x[i][j] = lower[i][j];
+ }
+
+ if (x[i][j] > upper[i][j]) {
+ x[i][j] = upper[i][j];
+ }
+ }
+ }
+
+ solverFactory_->setProblem(x);
+
+ solver.preprocess();
+ solver.solve();
+
+ return solver.getResult();
+ }
+
+private:
+ std::shared_ptr<Functional> functional_;
+ std::shared_ptr<Factory> solverFactory_;
+
+ Norm norm_;
+};
+
+#endif