diff --git a/dune/tectonic/spatial-solving/tnnmg/functional_old.hh b/dune/tectonic/spatial-solving/tnnmg/functional_old.hh
new file mode 100755
index 0000000000000000000000000000000000000000..15674fb378df4337709f4f3244f04a69cf76765e
--- /dev/null
+++ b/dune/tectonic/spatial-solving/tnnmg/functional_old.hh
@@ -0,0 +1,505 @@
+// -*- 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>
+
+#include "../../utils/maxeigenvalue.hh"
+#include "../../utils/debugutils.hh"
+
+template<class M, class V, class N, class L, class U, class R>
+class DirectionalRestriction;
+
+template<class M, class E, class V, class N, class L, class U, class O, class R>
+class ShiftedFunctional;
+
+template <class M, class V, class N, class L, class U, class R>
+class Functional;
+
+/** \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 E, class V, class N, class L, class U, class O, class R>
+class ShiftedFunctional {
+
+public:
+ using Matrix = std::decay_t<M>;
+ using Eigenvalues = E;
+ 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 Origin = std::decay_t<O>;
+ using Range = R;
+
+ template <class MM, class VV, class LL, class UU, class OO>
+ ShiftedFunctional(MM&& matrix, const Eigenvalues& maxEigenvalues, VV&& linearPart, const Nonlinearity& phi,
+ LL&& lower, UU&& upper, OO&& origin) :
+ quadraticPart_(std::forward<MM>(matrix)),
+ maxEigenvalues_(maxEigenvalues),
+ originalLinearPart_(std::forward<VV>(linearPart)),
+ phi_(phi),
+ originalLowerObstacle_(std::forward<LL>(lower)),
+ originalUpperObstacle_(std::forward<UU>(upper)),
+ origin_(std::forward<OO>(origin))
+ {}
+
+ Range operator()(const Vector& v) const
+ {
+ auto w = origin();
+ w += v;
+ if (Dune::TNNMG::checkInsideBox(w, originalLowerObstacle(), originalUpperObstacle())) {
+ Vector temp;
+ Dune::Solvers::resizeInitializeZero(temp,v);
+ quadraticPart().umv(w, temp);
+ temp *= 0.5;
+ temp -= originalLinearPart();
+ return temp*w + phi()(w);
+ }
+ return std::numeric_limits<Range>::max();
+
+ /*Vector temp;
+ Dune::Solvers::resizeInitializeZero(temp,v);
+ quadraticPart().umv(w, temp);
+ temp *= 0.5;
+ temp -= originalLinearPart();
+ return temp*w + phi()(w);*/
+ }
+
+ const Matrix& quadraticPart() const
+ {
+ return quadraticPart_.get();
+ }
+
+ const Vector& originalLinearPart() const
+ {
+ return originalLinearPart_.get();
+ }
+
+ const Origin& origin() const
+ {
+ return origin_.get();
+ }
+
+ void updateOrigin()
+ {}
+
+ void updateOrigin(std::size_t i)
+ {}
+
+ const LowerObstacle& originalLowerObstacle() const
+ {
+ return originalLowerObstacle_.get();
+ }
+
+ const UpperObstacle& originalUpperObstacle() const
+ {
+ return originalUpperObstacle_.get();
+ }
+
+ const auto& phi() const {
+ return phi_;
+ }
+
+ const auto& maxEigenvalues() const {
+ return maxEigenvalues_;
+ }
+
+protected:
+ Dune::Solvers::ConstCopyOrReference<M> quadraticPart_;
+ const Eigenvalues& maxEigenvalues_;
+
+ Dune::Solvers::ConstCopyOrReference<V> originalLinearPart_;
+ const Nonlinearity& phi_;
+
+ Dune::Solvers::ConstCopyOrReference<L> originalLowerObstacle_;
+ Dune::Solvers::ConstCopyOrReference<U> originalUpperObstacle_;
+
+ Dune::Solvers::ConstCopyOrReference<O> origin_;
+};
+
+/** \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, const double scaling) :
+ Base(matrix, linearTerm, lower, upper, origin, direction),
+ origin_(origin),
+ direction_(direction),
+ phi_(phi),
+ scaling_(scaling)
+ {
+ 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_;
+ }
+
+ double scaling() const {
+ return scaling_;
+ }
+protected:
+ const GlobalVector& origin_;
+ GlobalVector direction_;
+
+ const Nonlinearity& phi_;
+ Interval domain_;
+
+ const double scaling_;
+};
+
+
+/** \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 N, class V, class R>
+class FirstOrderFunctional {
+ using Interval = typename Dune::Solvers::Interval<R>;
+
+public:
+ using Nonlinearity = std::decay_t<N>;
+ using Vector = V;
+ using LowerObstacle = V;
+ using UpperObstacle = V;
+ using Range = R;
+
+ using field_type = typename V::field_type;
+public:
+ template <class LL, class UU, class OO, class DD>
+ FirstOrderFunctional(
+ const Range& maxEigenvalue,
+ const Range& linearPart,
+ const Nonlinearity& phi,
+ LL&& lower,
+ UU&& upper,
+ OO&& origin,
+ DD&& direction) :
+ quadraticPart_(maxEigenvalue),
+ linearPart_(linearPart),
+ lower_(std::forward<LL>(lower)),
+ upper_(std::forward<UU>(upper)),
+ origin_(std::forward<OO>(origin)),
+ direction_(std::forward<DD>(direction)),
+ phi_(phi) {
+
+ // set defect obstacles
+ defectLower_ = -std::numeric_limits<field_type>::max();
+ defectUpper_ = std::numeric_limits<field_type>::max();
+ Dune::TNNMG::directionalObstacles(origin_.get(), direction_.get(), lower_.get(), upper_.get(), defectLower_, defectUpper_);
+
+ // set domain
+ phi_.directionalDomain(origin_.get(), direction_.get(), domain_);
+
+ if (domain_[0] < defectLower_) {
+ domain_[0] = defectLower_;
+ }
+
+ if (domain_[1] > defectUpper_) {
+ domain_[1] = defectUpper_;
+ }
+ }
+
+ Range operator()(const Vector& v) const
+ {
+ DUNE_THROW(Dune::NotImplemented, "Evaluation of FirstOrderFunctional not implemented");
+ }
+
+ auto subDifferential(double x) const {
+ Interval Di;
+ Vector uxv = origin_.get();
+
+ Dune::MatrixVector::addProduct(uxv, x, direction_.get());
+ phi_.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 auto& origin() const {
+ return origin_.get();
+ }
+
+ const auto& direction() const {
+ return direction_.get();
+ }
+
+ const auto& lowerObstacle() const {
+ return defectLower_;
+ }
+
+ const auto& upperObstacle() const {
+ return defectUpper_;
+ }
+
+ const auto& quadraticPart() const {
+ return quadraticPart_;
+ }
+
+ const auto& linearPart() const {
+ return linearPart_;
+ }
+private:
+ const Range quadraticPart_;
+ const Range linearPart_;
+
+ Dune::Solvers::ConstCopyOrReference<LowerObstacle> lower_;
+ Dune::Solvers::ConstCopyOrReference<UpperObstacle> upper_;
+
+ Dune::Solvers::ConstCopyOrReference<Vector> origin_;
+ Dune::Solvers::ConstCopyOrReference<Vector> direction_;
+
+ const Nonlinearity& phi_;
+
+ Interval domain_;
+
+ 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 GlobalShiftedFunctional, class Index>
+auto coordinateRestriction(const GlobalShiftedFunctional& f, const Index& i)
+{
+ using Range = typename GlobalShiftedFunctional::Range;
+ 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;
+
+ //print(f.originalLinearPart(), "f.linearPart: ");
+ //print(f.quadraticPart(), "f.quadraticPart: ");
+
+ 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>());
+ });
+
+ //print(*Aii_p, "Aii_p:");
+ //print(ri, "ri:");
+
+ //print(f.originalLowerObstacle()[i], "lower:");
+ //print(f.originalUpperObstacle()[i], "upper:");
+
+ auto& phii = f.phi().restriction(i);
+
+ return ShiftedFunctional<LocalMatrix&, Range, LocalVector, std::decay_t<decltype(phii)>, LocalLowerObstacle&, LocalUpperObstacle&, LocalVector&, Range>(*Aii_p, f.maxEigenvalues()[i], std::move(ri), phii, f.originalLowerObstacle()[i], f.originalUpperObstacle()[i], f.origin()[i]);
+}
+
+
+/** \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;
+
+ using Eigenvalues = std::vector<Range>;
+
+private:
+ Dune::Solvers::ConstCopyOrReference<N> phi_;
+ Eigenvalues maxEigenvalues_;
+
+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)),
+ maxEigenvalues_(this->linearPart().size())
+ {
+ for (size_t i=0; i<this->quadraticPart().N(); ++i) {
+ const auto& Aii = this->quadraticPart()[i][i];
+ maxEigenvalues_[i] = maxEigenvalue(Aii);
+
+ // due to numerical imprecision, Aii might not be symmetric leading to complex eigenvalues
+ // consider Toeplitz decomposition of Aii and take only symmetric part
+
+ /*auto symBlock = Aii;
+
+ for (size_t j=0; j<symBlock.N(); j++) {
+ for (size_t k=0; k<symBlock.M(); k++) {
+ if (symBlock[j][k]/symBlock[j][j] < 1e-14)
+ symBlock[j][k] = 0;
+ }
+ }
+
+ try {
+ typename Vector::block_type eigenvalues;
+ Dune::FMatrixHelp::eigenValues(symBlock, eigenvalues);
+ maxEigenvalues_[i] =
+ *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+ } catch (Dune::Exception &e) {
+ print(symBlock, "symBlock");
+ typename Vector::block_type eigenvalues;
+ Dune::FMatrixHelp::eigenValues(symBlock, eigenvalues);
+ std::cout << "complex eig in dof: " << i << std::endl;
+ } */
+ }
+ }
+
+ const Nonlinearity& phi() const {
+ return phi_.get();
+ }
+
+ const auto& maxEigenvalues() const {
+ return maxEigenvalues_;
+ }
+
+ Range operator()(const Vector& v) const
+ {
+ //print(v, "v:");
+ //print(this->lowerObstacle(), "lower: ");
+ //print(this->upperObstacle(), "upper: ");
+ //std::cout << Dune::TNNMG::checkInsideBox(v, this->lowerObstacle(), this->upperObstacle()) << " " << Dune::TNNMG::QuadraticFunctional<M,V,R>::operator()(v) << " " << phi_.get().operator()(v) << std::endl;
+ if (Dune::TNNMG::checkInsideBox(v, this->lowerObstacle(), this->upperObstacle()))
+ return Dune::TNNMG::QuadraticFunctional<M,V,R>::operator()(v) + phi_.get()(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>
+ {
+ // rescale direction for numerical stability
+ auto dir = direction;
+ auto dirNorm = dir.two_norm();
+ if (dirNorm > 0.0)
+ dir /= dirNorm;
+ else {
+ dirNorm = 0.0;
+ dir = 0.0;
+ }
+
+
+ return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin, dir, dirNorm);
+ }
+
+ friend auto shift(const Functional& f, const Vector& origin)
+ -> ShiftedFunctional<Matrix&, Eigenvalues, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>
+ {
+ return ShiftedFunctional<Matrix&, Eigenvalues, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>(f.quadraticPart(), f.maxEigenvalues(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
+ }
+};
+
+
+#endif
diff --git a/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver_old.hh b/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver_old.hh
new file mode 100644
index 0000000000000000000000000000000000000000..454d898edba5f8f4de6f618e6899e50d449ed57a
--- /dev/null
+++ b/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver_old.hh
@@ -0,0 +1,130 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
+#define DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
+
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include "functional.hh"
+
+#include "../../utils/debugutils.hh"
+
+/**
+ * \brief A local solver for quadratic obstacle problems with nonlinearity
+ * using bisection
+ */
+class LocalBisectionSolver
+{
+public:
+ template<class Vector, class Functional, class BitVector>
+ void operator()(Vector& x, const Functional& f, const BitVector& ignore) const {
+ double safety = 1e-14; //or (f.upperObstacle()-f.lowerObstacle()<safety)
+ if (ignore.all()) {
+ return;
+ }
+
+ auto maxEigenvalue = f.maxEigenvalues();
+
+ auto origin = f.origin();
+ auto linearPart = f.originalLinearPart();
+
+ Dune::MatrixVector::addProduct(linearPart, maxEigenvalue, origin);
+
+ for (size_t j = 0; j < ignore.size(); ++j)
+ if (ignore[j])
+ linearPart[j] = 0; // makes sure result remains in subspace after correction
+ else
+ origin[j] = 0; // shift affine offset
+
+ double const linearPartNorm = linearPart.two_norm();
+ if (linearPartNorm <= 0.0)
+ return;
+ linearPart /= linearPartNorm;
+
+ std::cout << "maxEigenvalue " << maxEigenvalue << std::endl;
+ print(linearPart, "linearPart:");
+ std::cout << "linearPartNorm " << linearPartNorm << std::endl;
+ print(origin, "origin:");
+
+ using Nonlinearity = std::decay_t<decltype(f.phi())>;
+ using Range = std::decay_t<decltype(f.maxEigenvalues())>;
+ using FirstOrderFunctional = FirstOrderFunctional<Nonlinearity, Vector, Range>;
+
+ auto lower = f.originalLowerObstacle();
+ auto upper = f.originalUpperObstacle();
+
+ //print(lower, "lower:");
+ //print(upper, "upper:");
+
+ FirstOrderFunctional firstOrderFunctional(maxEigenvalue, linearPartNorm, f.phi(),
+ lower, upper,
+ origin, linearPart);
+
+ //std::cout << "lower: " << firstOrderFunctional.lowerObstacle() << " upper: " << firstOrderFunctional.upperObstacle() << std::endl;
+
+ // scalar obstacle solver without nonlinearity
+ typename Functional::Range alpha;
+ //Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
+ //obstacleSolver(alpha, firstOrderFunctional, false);
+
+ //direction *= alpha;
+
+ /*const auto& A = f.quadraticPart();
+ std::cout << "f.quadratic(): " << std::endl;
+ for (size_t i=0; i<A.N(); i++) {
+ for (size_t j=0; j<A.M(); j++) {
+ std::cout << A[i][j] << " ";
+ }
+ std::cout << std::endl;
+ }*/
+ //std::cout << f.quadraticPart() << std::endl;
+ //std::cout << "A: " << directionalF.quadraticPart() << " " << (directionalF.quadraticPart()==f.quadraticPart()[0][0]) << std::endl;
+ /*std::cout << "b: " << directionalF.linearPart() << std::endl;
+ std::cout << "domain: " << directionalF.domain()[0] << " " << directionalF.domain()[1] << std::endl;
+ auto D = directionalF.subDifferential(0);
+ std::cout << "subDiff at x: " << D[0] << " " << D[1] << std::endl;*/
+
+ //std::cout << "domain: " << firstOrderFunctional.domain()[0] << " " << firstOrderFunctional.domain()[1] << std::endl;
+
+
+ const auto& domain = firstOrderFunctional.domain();
+ if (std::abs(domain[0]-domain[1])>safety) {
+ int bisectionsteps = 0;
+ const Bisection globalBisection(0.0, 1.0, 0.0, 0.0);
+
+ alpha = globalBisection.minimize(firstOrderFunctional, 0.0, 0.0, bisectionsteps);
+ } else {
+ alpha = domain[0];
+ }
+
+ std::cout << "alpha: " << alpha << std::endl;
+ linearPart *= alpha;
+
+ std::cout << linearPart << std::endl;
+
+#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
+ if (std::abs(alpha)> safety) {
+ x = origin;
+ x += linearPart;
+ } else {
+ x = f.origin();
+ }
+#else
+ if (std::abs(alpha)> safety) {
+ x = origin;
+ x += linearPart;
+ x -= f.origin();
+ }
+#endif
+ std::cout << "new x: " << x << std::endl;
+ }
+};
+
+#endif
+
+
+
+