.

src/CMakeLists.txt

@@ -106,6 +106,8 @@ foreach(_dim 2 3)
 
   set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
   set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  #set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
   set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
   set_property(TARGET ${_sft_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  #set_property(TARGET ${_sft_target} APPEND PROPERTY COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
 endforeach()

src/multi-body-problem-2D.cfg

 # -*- mode:conf -*-
 [boundary.friction]
-smallestDiameter = 0.05  # 2e-3 [m]
+smallestDiameter = 0.01  # 2e-3 [m]
 
 [timeSteps]
 refinementTolerance = 1e-5 # 1e-5

src/multi-body-problem.cfg

@@ -47,10 +47,10 @@ relativeTau = 1e-4 # 1e-6
 
 [timeSteps]
 scheme = newmark
-timeSteps = 5
+timeSteps = 1
 
 [u0.solver]
-maximumIterations = 100
+maximumIterations = 20
 verbosity         = full
 
 [a0.solver]

src/solverfactorytest.cc

@@ -74,6 +74,31 @@ size_t timeSteps = 0;
 const std::string path = "";
 const std::string outputFile = "solverfactorytest.log";
 
+std::vector<std::vector<double>> allReductionFactors;
+
+template<class IterationStepType, class NormType, class ReductionFactorContainer>
+Dune::Solvers::Criterion reductionFactorCriterion(
+      const IterationStepType& iterationStep,
+      const NormType& norm,
+      ReductionFactorContainer& reductionFactors)
+{
+  double normOfOldCorrection = 1;
+  auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+  return Dune::Solvers::Criterion(
+      [&, lastIterate, normOfOldCorrection] () mutable {
+        double normOfCorrection = norm.diff(*lastIterate, *iterationStep.getIterate());
+        double convRate = (normOfOldCorrection > 0) ? normOfCorrection / normOfOldCorrection : 0.0;
+
+        normOfOldCorrection = normOfCorrection;
+        *lastIterate = *iterationStep.getIterate();
+
+        reductionFactors.push_back(convRate);
+        return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+      },
+      " reductionFactor");
+}
+
 template <class ContactNetwork>
 void solveProblem(const ContactNetwork& contactNetwork,
                   const Matrix& mat, const Vector& rhs, Vector& x,
@@ -141,15 +166,15 @@ void solveProblem(const ContactNetwork& contactNetwork,
 
     refSolver.addCriterion(
             [&](){
-            return Dune::formatString("   % 12.5e", I(x));
+            return Dune::formatString("   % 12.5e", I(refX));
             },
             "   energy      ");
 
-    double initialEnergy = I(x);
+    double initialEnergy = I(refX);
     refSolver.addCriterion(
             [&](){
             static double oldEnergy=initialEnergy;
-            double currentEnergy = I(x);
+            double currentEnergy = I(refX);
             double decrease = currentEnergy - oldEnergy;
             oldEnergy = currentEnergy;
             return Dune::formatString("   % 12.5e", decrease);
@@ -169,13 +194,13 @@ void solveProblem(const ContactNetwork& contactNetwork,
             },
             "   truncated   ");
 
-
-    std::vector<double> correctionNorms;
-    refSolver.addCriterion(Dune::Solvers::correctionNormCriterion(step, norm, 1e-10, correctionNorms));
+    if (timeStep>0 and initial) {
+        allReductionFactors[timeStep].resize(0);
+        refSolver.addCriterion(reductionFactorCriterion(step, norm, allReductionFactors[timeStep]));
+    }
 
     refSolver.preprocess();
     refSolver.solve();
-    std::cout << correctionNorms.size() << std::endl;
 
     /*if (initial) {
         x = refX;
@@ -216,6 +241,35 @@ void solveProblem(const ContactNetwork& contactNetwork,
                 parset.get<double>("u0.solver.tolerance"), &norm,
         Solver::FULL, true, &refX); // absolute error
 
+    solver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12.5e", J(x));
+            },
+            "   energy      ");
+
+    initialEnergy = J(x);
+    solver.addCriterion(
+            [&](){
+            static double oldEnergy=initialEnergy;
+            double currentEnergy = J(x);
+            double decrease = currentEnergy - oldEnergy;
+            oldEnergy = currentEnergy;
+            return Dune::formatString("   % 12.5e", decrease);
+            },
+            "   decrease    ");
+
+    solver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12.5e", step.lastDampingFactor());
+            },
+            "   damping     ");
+
+    solver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12d", tnnmgStep->linearization().truncated().count());
+            },
+            "   truncated   ");
+
     solver.preprocess();
     solver.solve();
 
@@ -479,7 +533,7 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
     // contribution from nonlinearity
     globalFriction.updateAlpha(alpha);
 
-    solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, lower, upper);
+    solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, lower, upper, false);
 
     nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
 
@@ -619,6 +673,8 @@ int main(int argc, char *argv[]) { try {
     // -----------------
     // init input/output
     // -----------------
+    timeSteps = parset.get<size_t>("timeSteps.timeSteps");
+    allReductionFactors.resize(timeSteps+1);
 
     setupInitialConditions(contactNetwork);
 
@@ -739,7 +795,6 @@ int main(int argc, char *argv[]) { try {
 
     //const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
 
-    timeSteps = parset.get<size_t>("timeSteps.timeSteps");
     for (timeStep=1; timeStep<=timeSteps; timeStep++) {
 
       advanceStep(current, contactNetwork);
@@ -756,6 +811,29 @@ int main(int argc, char *argv[]) { try {
       report();
     }
 
+    // output reduction factors
+    size_t count = 0;
+    for (size_t i=0; i<allReductionFactors.size(); i++) {
+        count = std::max(count, allReductionFactors[i].size());
+    }
+    std::vector<double> avgReductionFactors(count);
+    for (size_t i=0; i<count; i++) {
+        avgReductionFactors[i] = 1;
+        size_t c = 0;
+
+        for (size_t j=0; j<allReductionFactors.size(); j++) {
+            if (!(i<allReductionFactors[j].size()))
+                continue;
+
+            avgReductionFactors[i] *= allReductionFactors[j][i];
+            c++;
+        }
+
+        avgReductionFactors[i] = std::pow(avgReductionFactors[i], 1.0/((double)c));
+    }
+
+    print(avgReductionFactors, "average reduction factors: ");
+
     bool passed = true;
 
     std::cout << "Overall, the test " << (passed ? "was successful!" : "failed!") << std::endl;

src/spatial-solving/solverfactory.cc

@@ -26,9 +26,9 @@ SolverFactory<Functional, BitVector>::SolverFactory(
 
     auto linearSolver_ptr = Dune::Solvers::wrap_own_share<std::decay_t<LinearSolver>>(std::forward<LinearSolver>(linearSolver));
 
-    //tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, 10, DefectProjection(), LineSearchSolver());
+    tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, DefectProjection(), Dune::TNNMG::ScalarObstacleSolver());
 
-    tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, DefectProjection(), LineSearchSolver());
+    //tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, DefectProjection(), LineSearchSolver());
     tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
     tnnmgStep_->setIgnore(ignoreNodes);
 }

src/spatial-solving/solverfactory.hh

 #ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
 #define SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
 
-#define NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
+//#define NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
 
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/parametertree.hh>
@@ -34,7 +34,7 @@ class SolverFactory {
     using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
 
     //using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
-    using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
+    using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, Dune::TNNMG::ScalarObstacleSolver>;
 
   template <class LinearSolver>
   SolverFactory(const Dune::ParameterTree&,

src/spatial-solving/tnnmg/functional.hh

@@ -16,6 +16,9 @@
 
 #include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
 
+template<class M, class V, class N, class L, class U, class R>
+class DirectionalRestriction;
+
 template<class M, class V, class N, class R>
 class ShiftedFunctional;
 
@@ -60,6 +63,12 @@ class ShiftedFunctional : public Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunc
       return phi_;
   }
 
+  friend auto directionalRestriction(const ShiftedFunctional& 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.originalLinearPart(), f.phi(), f.originalLowerObstacle(), f.originalUpperObstacle(), origin, direction);
+  }
+
 protected:
   const Nonlinearity& phi_;
 };
@@ -148,6 +157,7 @@ class DirectionalRestriction :
   Interval domain_;
 };
 
+
 /** \brief Box constrained quadratic functional with nonlinearity
  *         Note: call setIgnore() to set up functional in affine subspace with ignore information
  *
@@ -157,7 +167,7 @@ class DirectionalRestriction :
  *  \tparam U Upper obstacle type
  *  \tparam R Range type
  */
-template<class V, class N, class L, class U, class O, class D, class R>
+/*template<class V, class N, class L, class U, class O, class D, class R>
 class FirstOrderModelFunctional {
     using Interval = typename Dune::Solvers::Interval<R>;
 
@@ -170,10 +180,30 @@ class FirstOrderModelFunctional {
 
     using field_type = typename Vector::field_type;
 
-private:
-    void init() {
-        auto origin = origin_.get();
-        auto direction = direction_.get();
+public:
+    template <class VV1, class NN, class LL, class UU, class VV2>
+    FirstOrderModelFunctional(
+            const Range& maxEigenvalue,
+            VV1&& linearTerm,
+            NN&& phi,
+            LL&& lower,
+            UU&& upper,
+            VV2&& origin) :
+        maxEigenvalue_(maxEigenvalue),
+        linearTerm_(std::forward<VV1>(linearTerm)),
+        lower_(std::forward<LL>(lower)),
+        upper_(std::forward<UU>(upper)),
+        origin_(std::forward<VV2>(origin)),
+        //direction_(linearTerm_),
+        phi_(std::forward<NN>(phi)) {
+
+        auto& origin = origin_.get();
+        auto& direction = direction_.get();
+
+        auto v = ri;
+        double const vnorm = v.two_norm();
+        if (vnorm > 1.0)
+              v /= vnorm;
 
         // set defect obstacles
         defectLower_ = -std::numeric_limits<field_type>::max();
@@ -190,61 +220,10 @@ class FirstOrderModelFunctional {
         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 {
+  auto setIgnore(const BitVector& ignore) const {
       Vector direction = direction_.get();
       Vector origin = origin_.get();
      // Dune::Solvers::resizeInitializeZero(direction, linearPart());
@@ -254,7 +233,9 @@ class FirstOrderModelFunctional {
         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));
+      // set quadratic and linear parts of functional
+      quadraticPart_ = maxEigenvalue_ * direction.two_norm2();
+      linearPart_ = linearTerm_.get() * direction;
   }
 
   Range operator()(const Vector& v) const
@@ -315,7 +296,7 @@ class FirstOrderModelFunctional {
 
   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
@@ -345,6 +326,18 @@ auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, co
       Dune::TNNMG::Imp::mmv(Aij, f.origin()[j], ri, Dune::PriorityTag<1>());
   });
 
+  // set maxEigenvalue from matrix
+  LocalVector eigenvalues;
+  Dune::FMatrixHelp::eigenValues(*Aii_p, eigenvalues);
+  auto maxEigenvalue = *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+
+  Dune::MatrixVector::addProduct(ri, maxEigenvalue, f.origin()[i]);
+
+  LocalMatrix Aii = 0;
+  for (size_t d=0; d<Aii.N(); d++) {
+      Aii[d][d] = maxEigenvalue;
+  }
+
   LocalLowerObstacle dli = f.originalLowerObstacle()[i];
   LocalUpperObstacle dui = f.originalUpperObstacle()[i];
   dli -= f.origin()[i];
@@ -352,13 +345,16 @@ auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, co
 
   auto&& phii = f.phi().restriction(i);
 
-  auto v = ri;
-  double const vnorm = v.two_norm();
-  if (vnorm > 1.0)
-        v /= vnorm;
+  /*std::cout << "maxEigenvalue matrix: " << std::endl;
+  for (size_t i=0; i<Aii.N(); i++) {
+        for (size_t j=0; j<Aii.M(); j++) {
+            std::cout << Aii[i][j] << " ";
+        }
+        std::cout << std::endl;
+  }*/
 
-  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));
-  //return Functional<LocalMatrix&, LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, Range>(*Aii_p, std::move(ri), std::move(phii), std::move(dli), std::move(dui));
+  //return FirstOrderModelFunctional<LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, LocalVector, Range>(maxEigenvalue, std::move(ri), std::move(phii), std::move(dli), std::move(dui), std::move(f.origin()[i]));
+  return ShiftedFunctional<LocalMatrix, LocalVector, decltype(phii), Range>(std::move(Aii), std::move(ri), std::move(phii), std::move(dli), std::move(dui), f.origin()[i]);
 }
 
 

src/spatial-solving/tnnmg/localbisectionsolver.hh

@@ -3,7 +3,9 @@
 #ifndef DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
 #define DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
 
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
 #include <dune/tnnmg/problem-classes/bisection.hh>
+#include "../../utils/debugutils.hh"
 
 /**
  * \brief A local solver for quadratic obstacle problems with nonlinearity
@@ -14,48 +16,67 @@ 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;
       }
 
-      double safety = 1e-14;
-      if (f.defectUpper() - f.defectLower() < safety) {
-          x = f.origin();
+      Vector origin = f.origin();
+      Vector direction = f.originalLinearPart();
+
+      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
+
+      double const directionNorm = direction.two_norm();
+      if (directionNorm <= 0.0)
         return;
+      direction /= directionNorm;
+
+      auto&& directionalF = directionalRestriction(f, origin, direction);
+
+      // scalar obstacle solver without nonlinearity
+      typename Functional::Range alpha;
+      Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
+      obstacleSolver(alpha, directionalF, 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;*/
 
+      const auto& domain = directionalF.domain();
+      if (std::abs(domain[0]-domain[1])>safety) {
         int bisectionsteps = 0;
         const Bisection globalBisection;
 
-      if (ignore.any()) {
-            auto&& restrictedf = f.getIgnoreFunctional(ignore);
-            const auto alpha = globalBisection.minimize(restrictedf, 0.0, 0.0, bisectionsteps);
-
-#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
-            x = restrictedf.origin();
-            auto correction = restrictedf.direction();
-            correction *= alpha;
-            x += correction;
-#else
-            auto x = restrictedf.direction();
-            x *= alpha;
-#endif
-
+        alpha = globalBisection.minimize(directionalF, alpha, 0.0, bisectionsteps);
+      }
 
-      } else { // ignore.none()
-          const auto alpha = globalBisection.minimize(f, 0.0, 0.0, bisectionsteps);
+      direction *= alpha;
 
 #ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
-          x = f.origin();
-          auto correction = f.direction();
-          correction *= alpha;
-          x += correction;
+      x = origin;
+      x += direction;
 #else
-          auto x = f.direction();
-          x *= alpha;
+      //x = origin;
+      x -= f.origin();
+      x += direction;
 #endif
-
-      }
   }
 };