.

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

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

src/multi-body-problem.cfg

@@ -11,7 +11,7 @@ vtk.write       = true
 
 [problem]
 finalTime       = 100  # [s] #1000
-bodyCount       = 4
+bodyCount       = 2
 
 [body]
 bulkModulus     = 0.5e5 # [Pa]

src/spatial-solving/tnnmg/functional.hh

@@ -16,10 +16,12 @@
 
 #include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
 
+#include "../../utils/debugutils.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>
+template<class M, class E, class V, class N, class R>
 class ShiftedFunctional;
 
 template <class M, class V, class N, class L, class U, class R>
@@ -33,21 +35,23 @@ class Functional;
  *  \tparam V nonlinearity type
  *  \tparam R Range type
  */
-template<class M, class V, class N, class R>
+template<class M, class E, 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 Eigenvalues = E;
   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) :
+  ShiftedFunctional(const Matrix& quadraticPart, const Eigenvalues& maxEigenvalues, const Vector& linearPart, const Nonlinearity& phi, const LowerObstacle& lowerObstacle, const UpperObstacle& upperObstacle, const Vector& origin) :
     Base(quadraticPart, linearPart, lowerObstacle, upperObstacle, origin),
-    phi_(phi)
+    phi_(phi),
+    maxEigenvalues_(maxEigenvalues)
   {}
 
   Range operator()(const Vector& v) const
@@ -63,14 +67,13 @@ 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);
+  const auto& maxEigenvalues() const {
+    return maxEigenvalues_;
   }
 
 protected:
   const Nonlinearity& phi_;
+  const Eigenvalues& maxEigenvalues_;
 };
 
 
@@ -167,80 +170,56 @@ 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>
-class FirstOrderModelFunctional {
+template<class N, class V, class R>
+class FirstOrderFunctional {
     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 Vector = V;
+    using LowerObstacle = V;
+    using UpperObstacle = V;
     using Range = R;
 
-    using field_type = typename Vector::field_type;
-
+    using field_type = typename V::field_type;
 public:
-    template <class VV1, class NN, class LL, class UU, class VV2>
-    FirstOrderModelFunctional(
+    template <class NN, class LL, class UU, class OO, class DD>
+    FirstOrderFunctional(
             const Range& maxEigenvalue,
-            VV1&& linearTerm,
+            const Range& linearPart,
             NN&& phi,
             LL&& lower,
             UU&& upper,
-            VV2&& origin) :
-        maxEigenvalue_(maxEigenvalue),
-        linearTerm_(std::forward<VV1>(linearTerm)),
+            OO&& origin,
+            DD&& direction) :
+        quadraticPart_(maxEigenvalue),
+        linearPart_(linearPart),
         lower_(std::forward<LL>(lower)),
         upper_(std::forward<UU>(upper)),
-        origin_(std::forward<VV2>(origin)),
-        //direction_(linearTerm_),
+        origin_(std::forward<OO>(origin)),
+        direction_(std::forward<DD>(direction)),
         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();
         defectUpper_ = std::numeric_limits<field_type>::max();
-        Dune::TNNMG::directionalObstacles(origin, direction, lower_.get(), upper_.get(), defectLower_, defectUpper_);
+        Dune::TNNMG::directionalObstacles(origin_.get(), direction_.get(), lower_.get(), upper_.get(), defectLower_, defectUpper_);
 
         // set domain
-        phi_.get().directionalDomain(origin, direction, domain_);
+        phi_.get().directionalDomain(origin_.get(), direction_.get(), domain_);
 
-        if (domain_[0] < defectLower_) {
+        /*if (domain_[0] < defectLower_) {
            domain_[0] = defectLower_;
         }
 
         if (domain_[1] > defectUpper_) {
            domain_[1] = defectUpper_;
-        }
+        }*/
     }
 
-  template <class BitVector>
-  auto setIgnore(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
-
-      // set quadratic and linear parts of functional
-      quadraticPart_ = maxEigenvalue_ * direction.two_norm2();
-      linearPart_ = linearTerm_.get() * direction;
-  }
-
   Range operator()(const Vector& v) const
   {
-    DUNE_THROW(Dune::NotImplemented, "Evaluation of FirstOrderModelFunctional not implemented");
+    DUNE_THROW(Dune::NotImplemented, "Evaluation of FirstOrderFunctional not implemented");
   }
 
   auto subDifferential(double x) const {
@@ -261,49 +240,53 @@ class FirstOrderModelFunctional {
     return domain_;
   }
 
-  const Vector& origin() const {
+  const auto& origin() const {
       return origin_.get();
   }
 
-  const Vector& direction() const {
+  const auto& direction() const {
       return direction_.get();
   }
 
-  const field_type& defectLower() const {
+  const auto& lowerObstacle() const {
       return defectLower_;
   }
 
-  const field_type& defectUpper() const {
+  const auto& upperObstacle() const {
       return defectUpper_;
   }
 
+  const auto& quadraticPart() const {
+      return quadraticPart_;
+  }
+
+  const auto& linearPart() const {
+      return linearPart_;
+  }
 private:
-  Dune::Solvers::ConstCopyOrReference<V> linearTerm_;
-  Dune::Solvers::ConstCopyOrReference<L> lower_;
-  Dune::Solvers::ConstCopyOrReference<U> upper_;
+  const Range quadraticPart_;
+  const Range linearPart_;
+
+  Dune::Solvers::ConstCopyOrReference<LowerObstacle> lower_;
+  Dune::Solvers::ConstCopyOrReference<UpperObstacle> upper_;
 
-  Dune::Solvers::ConstCopyOrReference<O> origin_;
-  Dune::Solvers::ConstCopyOrReference<D> direction_;
+  Dune::Solvers::ConstCopyOrReference<Vector> origin_;
+  Dune::Solvers::ConstCopyOrReference<Vector> 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)
+template<class Index, class M, class E, class V, class Nonlinearity, class R>
+auto coordinateRestriction(const ShiftedFunctional<M, E, V, Nonlinearity, R>& f, const Index& i)
 {
   using Range = R;
   using LocalMatrix = std::decay_t<decltype(f.quadraticPart()[i][i])>;
@@ -316,6 +299,8 @@ auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, co
 
   const LocalMatrix* Aii_p = nullptr;
 
+  //print(f.origin(), "f.origin: ");
+
   LocalVector ri = f.originalLinearPart()[i];
   const auto& Ai = f.quadraticPart()[i];
   sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
@@ -326,35 +311,9 @@ 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];
-  dui -= f.origin()[i];
-
   auto&& phii = f.phi().restriction(i);
 
-  /*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>(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]);
+  return ShiftedFunctional<LocalMatrix, Range, LocalVector, decltype(phii), Range>(*Aii_p, f.maxEigenvalues()[i], std::move(ri), std::move(phii), f.originalLowerObstacle()[i], f.originalUpperObstacle()[i], f.origin()[i]);
 }
 
 
@@ -382,8 +341,11 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
   using LowerObstacle = typename Base::LowerObstacle;
   using UpperObstacle = typename Base::UpperObstacle;
 
+  using Eigenvalues = std::vector<Range>;
+
 private:
     Dune::Solvers::ConstCopyOrReference<N> phi_;
+    Eigenvalues maxEigenvalues_;
 
 public:
 
@@ -395,13 +357,25 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
           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))
-  {}
+    phi_(std::forward<NN>(phi)),
+    maxEigenvalues_(this->linearPart().size())
+  {
+     for (size_t i=0; i<this->quadraticPart().N(); ++i) {
+       typename Vector::block_type eigenvalues;
+       Dune::FMatrixHelp::eigenValues(this->quadraticPart()[i][i], eigenvalues);
+       maxEigenvalues_[i] =
+           *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+     }
+  }
 
   const Nonlinearity& phi() const {
     return phi_.get();
   }
 
+  const auto& maxEigenvalues() const {
+    return maxEigenvalues_;
+  }
+
   Range operator()(const Vector& v) const
   { 
     if (Dune::TNNMG::checkInsideBox(v, this->lower_.get(), this->upper_.get()))
@@ -416,9 +390,9 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
   }
 
   friend auto shift(const Functional& f, const Vector& origin)
-    -> ShiftedFunctional<Matrix, Vector, Nonlinearity, Range>
+    -> ShiftedFunctional<Matrix, Eigenvalues, Vector, Nonlinearity, Range>
   {
-    return ShiftedFunctional<Matrix, Vector, Nonlinearity, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
+    return ShiftedFunctional<Matrix, Eigenvalues, Vector, Nonlinearity, Range>(f.quadraticPart(), f.maxEigenvalues(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
   }
 };
 

src/spatial-solving/tnnmg/localbisectionsolver.hh

@@ -3,8 +3,13 @@
 #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"
 
 /**
@@ -21,26 +26,33 @@ class LocalBisectionSolver
           return;
       }
 
-      Vector origin = f.origin();
-      Vector direction = f.originalLinearPart();
+      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])
-          direction[j] = 0; // makes sure result remains in subspace after correction
+          linearPart[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)
+      double const linearPartNorm = linearPart.two_norm();
+      if (linearPartNorm <= 0.0)
         return;
-      direction /= directionNorm;
+      linearPart /= linearPartNorm;
 
-      auto&& directionalF = directionalRestriction(f, origin, direction);
+      using Nonlinearity = std::decay_t<decltype(f.phi())>;
+      using Range = std::decay_t<decltype(f.maxEigenvalues())>;
+      using FirstOrderFunctional = FirstOrderFunctional<Nonlinearity, Vector, Range>;
+
+      FirstOrderFunctional firstOrderFunctional(maxEigenvalue, linearPartNorm, f.phi(), f.originalLowerObstacle(), f.originalUpperObstacle(), origin, linearPart);
 
       // scalar obstacle solver without nonlinearity
       typename Functional::Range alpha;
-      Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
-      obstacleSolver(alpha, directionalF, false);
+      /*Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
+      obstacleSolver(alpha, firstOrderFunctional, false);*/
 
       //direction *= alpha;
 
@@ -59,23 +71,31 @@ class LocalBisectionSolver
       auto D = directionalF.subDifferential(0);
       std::cout << "subDiff at x: " << D[0] << " " << D[1] << std::endl;*/
 
-      const auto& domain = directionalF.domain();
+      const auto& domain = firstOrderFunctional.domain();
       if (std::abs(domain[0]-domain[1])>safety) {
         int bisectionsteps = 0;
-        const Bisection globalBisection;
+        const Bisection globalBisection(0.0, 1.0, 0.0, 0.0);
+
+        alpha = globalBisection.minimize(firstOrderFunctional, 0.0, 0.0, bisectionsteps);
+
+        if (alpha < firstOrderFunctional.lowerObstacle() ) {
+                   alpha = firstOrderFunctional.lowerObstacle();
+                }
 
-        alpha = globalBisection.minimize(directionalF, alpha, 0.0, bisectionsteps);
+         if (alpha > firstOrderFunctional.upperObstacle() ) {
+                   alpha = firstOrderFunctional.upperObstacle();
+         }
       }
 
-      direction *= alpha;
+      linearPart *= alpha;
 
 #ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
       x = origin;
-      x += direction;
+      x += linearPart;
 #else
       //x = origin;
-      x -= f.origin();
-      x += direction;
+      //x -= f.origin();
+      x = linearPart;
 #endif
   }
 };