diff --git a/dune/contact/solvers/primaldualsemismoothnewtonsolver.cc b/dune/contact/solvers/primaldualsemismoothnewtonsolver.cc
new file mode 100644
index 0000000000000000000000000000000000000000..01436e76b2dc3b3d5c09a7079e2f5b87a39ca821
--- /dev/null
+++ b/dune/contact/solvers/primaldualsemismoothnewtonsolver.cc
@@ -0,0 +1,198 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=4 sw=2 et sts=2:
+
+#include <limits>
+#include <dune/common/timer.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/solvers/solvers/quadraticipopt.hh>
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/umfpack.hh>
+#include <dune/istl/solvers.hh>
+
+namespace Dune {
+namespace Contact {
+
+template<class ProblemType>
+void PrimalDualSemiSmoothNewtonSolver<ProblemType>::solve()
+{
+ auto& iterates = *iterates_;
+
+ size_t totalSize(0);
+ for (size_t i=0; i<iterates.size(); i++)
+ totalSize += iterates[i].size();
+
+ // If no Dirichlet conditions are prescribed we have to setup the initial correction
+ if (correction_.size() != dim*totalSize) {
+ correction_.resize(dim*totalSize);
+ correction_ = 0;
+ }
+
+ problem_->contactAssembler()->getContactCouplings()[0]->path_ = debugPath_;
+ problem_->contactAssembler()->getContactCouplings()[0]->it_ = -1;
+
+ // compute energy infeasibility of initial iterate and build the contact mapping
+ energy_ = problem_->energy(iterates);
+ infeasibility_ = problem_->infeasibility(iterates);
+
+ // if there are non-homogeneous Dirichlet conditions,
+ // then we have to adjust the initial infeasibility
+ infeasibility_ = std::max(infeasibility_, correction_.infinity_norm());
+
+ if (this->verbosity_ == Solver::REDUCED or this->verbosity_== Solver::FULL)
+ std::cout<< " Initial energy: " << energy_
+ << " and infeasibility: " << infeasibility_ << std::endl;
+
+ // initialise the iteration counter and error
+ int it(0);
+ field_type error = std::numeric_limits<field_type>::max();
+
+ // assemble linearisations of initial deformation
+ problem_->assembleLinearisations(iterates);
+ // initialise Lagrange multiplier with zero
+ lagrangeMultiplier_.resize(problem_->getConstraints().size());
+ // lagrangeMultiplier_ = -1000;
+ lagrangeMultiplier_ = problem_->leastSquaresLagrangeMultiplier();
+
+ problem_->initialiseEmptyActiveSet();
+
+ // create IpOpt solver
+ QuadraticIPOptSolver<ScalarMatrix, ScalarVector> ipoptSolver(localTol_, localMaxIterations_,
+ Solver::QUIET, linearSolverType_);
+
+ bool activeSetChanged(true);
+ // Loop until desired tolerance or maximum number of iterations is reached
+ for (; activeSetChanged or (it<this->maxIterations_ and (error>this->tolerance_ or std::isnan(error))); it++)
+ {
+
+ if (this->verbosity_ == Solver::REDUCED or this->verbosity_== Solver::FULL)
+ std::cout << "************************\n"
+ << "* PriDual iteration " << it << "\n"
+ << "************************\n";
+
+ problem_->contactAssembler()->getContactCouplings()[0]->path_ = debugPath_;
+ problem_->contactAssembler()->getContactCouplings()[0]->it_ = it;
+
+ // Measure norm of the old iterate
+ field_type oldNorm(0);
+ for (size_t i=0; i<iterates_->size(); i++)
+ oldNorm += (*errorNorms_[i])(iterates[i]);
+
+ // assemble the local quadratic problem
+ problem_->assemblePrimalDualProblem(lagrangeMultiplier_);
+
+ ScalarVector priDualIterate = combineVector(correction_, ScalarVector(lagrangeMultiplier_.size(),0));
+
+ field_type newEnergy;
+ field_type newInfeasibility;
+ std::vector<BlockVector> newIterates;
+
+ // setup the linearised problem
+ problem_->addHardDirichlet(correction_, dirichletNodes_);
+ ipoptSolver.setProblem(problem_->priDualMatrix(), priDualIterate, problem_->priDualRhs());
+
+ Timer timer;
+
+#if HAVE_SUITESPARSE_UMFPACK
+ UMFPack<ScalarMatrix> umfPack(problem_->priDualMatrix(), 0);
+ std::cout<<"took "<<timer.elapsed()<<" to decompose matrix\n";
+
+ // solve
+ InverseOperatorResult res;
+ auto copyRhs = problem_->priDualRhs();
+ umfPack.apply(priDualIterate, copyRhs, res);
+ if (!res.converged)
+ std::cout<<"Warning: Not converged! Reduction "<<res.reduction<<" iterations "<<res.iterations<<std::endl;
+#endif
+
+ auto residual = problem_->priDualRhs();
+ problem_->priDualMatrix().mmv(priDualIterate, residual);
+ std::cout<<"Residual " << residual.infinity_norm()<<std::endl;
+
+
+
+ if (this->verbosity_ == Solver::FULL)
+ std::cout << "Solved local problem in " << timer.stop() << " secs\n";
+
+ std::tie(correction_, lagrangeMultiplier_) = decomposeVector(priDualIterate);
+
+ // new iterate
+ int counter(0);
+ newIterates = iterates;
+ for (size_t i=0; i<newIterates.size(); i++)
+ for (size_t j=0; j<newIterates[i].size(); j++)
+ for (int k = 0; k < dim; k++)
+ newIterates[i][j][k] += correction_[counter++][0];
+
+ std::cout<<"iterates inf "<<newIterates[0].infinity_norm() << " and " << newIterates[1].infinity_norm()<<std::endl;
+ // compute new energy and infeasibility
+ newEnergy = problem_->energy(newIterates);
+ newInfeasibility = problem_->infeasibility(newIterates);
+
+ // update the active set, needs to be called after the deformation is updated
+ // which is done in "infeasibility"
+ problem_->updateActiveSetAndLagrangeMultiplier(lagrangeMultiplier_);
+ activeSetChanged = (problem_->changedActiveNodes() > 0);
+
+ if (this->verbosity_ == Solver::FULL) {
+ std::cout << "Active nodes changed " << problem_->changedActiveNodes()
+ << " active now " << problem_->getActiveSet().count() << " \n";
+ std::cout << "New energy of potential iterate " << newEnergy
+ << " and infeasibility " << newInfeasibility << std::endl;
+ }
+
+ // update the Newton problem
+ problem_->assembleLinearisations(newIterates);
+
+ // if the correction is very small stop the loop
+ field_type infNorm = correction_.infinity_norm();
+ if(this->verbosity_==Solver::FULL)
+ std::cout << "Correction infinite norm: " << infNorm << std::endl;
+
+ // compute error estimate
+ std::vector<BlockVector> corrections(iterates.size());
+ counter = 0;
+ for (size_t i=0; i<iterates.size();i++) {
+ corrections[i].resize(iterates[i].size());
+ for (size_t j=0; j<corrections[i].size(); j++)
+ for (int k = 0; k < dim; k++)
+ corrections[i][j][k] = correction_[counter++][0];
+ }
+
+ auto residuals = problem_->residuals(lagrangeMultiplier_);
+ for (int k=0; k<2; k++)
+ for (size_t j=0; j<residuals[k].size();j++)
+ if (dN_[k][j].any())
+ residuals[k][j] = 0;
+
+ field_type absError = 0;
+ for (size_t i=0; i < corrections.size(); ++i)
+ absError += (*errorNorms_[i])(residuals[i]);
+ absError += newInfeasibility;
+
+ error = absError;
+ if (this->verbosity_ == Solver::REDUCED or this->verbosity_== Solver::FULL)
+ std::cout << "In filter iteration "<< it
+ << ", absolut norm of correction : " << absError
+ << " relative error " << error
+ << "\nEnergy " << newEnergy
+ << " and infeasibility " << newInfeasibility << std::endl;
+
+ // update iterates, energy and infeasibility
+ iterates = newIterates;
+ energy_ = newEnergy;
+ infeasibility_ = newInfeasibility;
+
+ // set correction to zero to remove any Dirichlet values
+ correction_ = 0;
+
+ } // end of filter loop
+
+ if (this->verbosity_ == Solver::REDUCED or this->verbosity_== Solver::FULL)
+ std::cout << "Solved problem in " << it << " iterations with final energy " << energy_
+ << " and infeasibility " << infeasibility_ << std::endl;
+}
+
+} /* namespace Contact */
+} /* namespace Dune */
diff --git a/dune/contact/solvers/primaldualsemismoothnewtonsolver.hh b/dune/contact/solvers/primaldualsemismoothnewtonsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..741e6020f3a0b693c6ff9242780b887ba7313e32
--- /dev/null
+++ b/dune/contact/solvers/primaldualsemismoothnewtonsolver.hh
@@ -0,0 +1,168 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=4 sw=2 et sts=2:
+#ifndef DUNE_CONTACT_SOLVERS_PRIMAL_DUAL_SEMI_SMOOTH_NEWTON_SOLVER_HH
+#define DUNE_CONTACT_SOLVERS_PRIMAL_DUAL_SEMI_SMOOTH_NEWTON_SOLVER_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/norms/norm.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+
+namespace Dune {
+namespace Contact {
+
+/** \brief Primal-Dual Semi-Smooth Newton method
+ *
+ * See Popp, Wall et al. 2011
+*/
+template <class Problem>
+class PrimalDualSemiSmoothNewtonSolver
+ : public Dune::Solvers::IterativeSolver<typename Problem::BlockVector> {
+
+ using BlockVector = typename Problem::BlockVector;
+ using field_type = typename Dune::FieldTraits<BlockVector>::field_type;
+ using ScalarVector = typename Problem::ScalarVector;
+ using ScalarMatrix = typename Problem::ScalarMatrix;
+ using Matrix = typename Problem::MatrixType;
+ using JacobianType = typename Problem::JacobianType;
+
+ static constexpr int dim = BlockVector::block_type::dimension;
+ using Base = Dune::Solvers::IterativeSolver<BlockVector>;
+
+public:
+ PrimalDualSemiSmoothNewtonSolver(
+ const Dune::ParameterTree& config, std::vector<BlockVector>& x,
+ Problem& problem,
+ const std::vector<std::shared_ptr<Norm<BlockVector> > >& errorNorms)
+ : Base(config.get<field_type>("tolerance"),
+ config.get<int>("maxIterations"),
+ config.get("verbosity", NumProc::FULL)) {
+ setupParameter(config);
+ setProblem(x, problem);
+ setErrorNorm(errorNorms);
+
+ size_t totalSize(0);
+ for (size_t i = 0; i < x.size(); i++)
+ totalSize += x[i].size();
+
+ dirichletNodes_.resize(totalSize*dim, false);
+ }
+
+ //! Setup the trust region parameter from a config file
+ void setupParameter(const Dune::ParameterTree& config) {
+ localTol_ = config.get<field_type>("localTolerance");
+ localMaxIterations_ = config.get<int>("localMaxIterations");
+ linearSolverType_ = config.get("linearSolverType", "");
+ }
+
+ /** \brief Set the Dirichlet values.
+ *
+ * The Dirichlet conditions are enforced during the solution
+ * of the first linearised problem rather than modifying the
+ * initial iterate of the filter scheme. This avoids possible problems
+ *in the energy evaluation due to inverted elements.
+ * Hence Dirichlet values for the correction must be handed over!
+ */
+ void setDirichletValues(
+ const std::vector<BlockVector>& dirichletValues,
+ const std::vector<Dune::BitSetVector<dim> >& dirichletNodes) {
+ // only adjust correction Dirichlet conditions are prescribed
+ correction_.resize(dirichletNodes_.size());
+ correction_ = 0;
+ int counter(0);
+ for (size_t i = 0; i < dirichletNodes.size(); i++)
+ for (size_t j = 0; j < dirichletNodes[i].size(); counter++, j++)
+ if (dirichletNodes[i][j].any())
+ for (int k = 0; k < dim; k++) {
+ dirichletNodes_[dim*counter + k] = dirichletNodes[i][j][k];
+ correction_[dim*counter + k][0] = dirichletValues[i][j][k];
+ }
+
+ dN_ = dirichletNodes;
+ }
+
+ //! Set initial iterate and the problem class
+ void setProblem(std::vector<BlockVector>& x, Problem& problem) {
+ iterates_ = Dune::stackobject_to_shared_ptr<std::vector<BlockVector> >(x);
+ problem_ = Dune::stackobject_to_shared_ptr<Problem>(problem);
+ }
+
+ //! Set the error norms to be used within the filter method
+ void
+ setErrorNorm(const std::vector<std::shared_ptr<Norm<BlockVector> > >& errorNorms) {
+ errorNorms_ = errorNorms;
+ }
+
+ //! Solve the problem
+ void solve();
+
+ //! Return solution object
+ const std::vector<BlockVector>& getSol() const { return *iterates_; }
+
+ std::string debugPath_;
+ std::vector<Dune::BitSetVector<dim> > dN_;
+
+private:
+ ScalarVector combineVector(const ScalarVector& primal, const ScalarVector& dual) const {
+ size_t index(0);
+ ScalarVector combinedVector(primal.size() + dual.size());
+ for (size_t i = 0; i < primal.size(); i++)
+ combinedVector[index++] = primal[i];
+
+ for (size_t i = 0; i < dual.size(); i++)
+ combinedVector[index++] = dual[i];
+ return combinedVector;
+ }
+
+ std::tuple<ScalarVector, ScalarVector> decomposeVector(const ScalarVector& primalDual) const {
+ size_t index(0);
+ ScalarVector primal(correction_.size());
+ for (size_t i = 0; i < primal.size(); i++)
+ primal[i] = primalDual[index++];
+
+ ScalarVector dual(lagrangeMultiplier_.size());
+ for (size_t i = 0; i < dual.size(); i++)
+ dual[i] = primalDual[index++];
+ return std::tie(primal, dual);
+ }
+
+ //! The iterates of the global scheme
+ std::shared_ptr<std::vector<BlockVector> > iterates_;
+
+ ScalarVector lagrangeMultiplier_;
+
+ //! Dirichlet degrees of freedom
+ Dune::BitSetVector<1> dirichletNodes_;
+
+ //! The iterate for the local problems, containing possible Dirichlet values
+ ScalarVector correction_;
+
+ //! The problem type, providing methods to assemble the nonlinear and
+ //! linearised
+ //! energy and constraints
+ std::shared_ptr<Problem> problem_;
+
+ //! The norms to estimate the errors
+ std::vector<std::shared_ptr<Norm<BlockVector> > > errorNorms_;
+
+ //! Tolerance for the solution of the local problem
+ field_type localTol_;
+ //! Max iterations for the solution of the local problem
+ int localMaxIterations_;
+
+ //! Store the old energy
+ field_type energy_;
+ //! Store the infeasibility of the old iterate
+ field_type infeasibility_;
+
+ //! The linear solver to be used within IpOpt
+ std::string linearSolverType_;
+};
+
+} /* namespace Contact */
+} /* namespace Dune */
+
+#include "primaldualsemismoothnewtonsolver.cc"
+
+#endif