.

CMakeLists.txt

@@ -17,7 +17,7 @@ include(DuneMacros)
 
 # start a dune project with information from dune.module
 dune_project()
-
+dune_enable_all_packages()
 find_package(HDF5 COMPONENTS C REQUIRED)
 
 add_subdirectory("src")

src/factories/stackedblocksfactory.cc

@@ -51,7 +51,7 @@ void StackedBlocksFactory<GridType, dims>::setBodies() {
 #elif MY_DIM == 2
         origins[0] = {0,0};
         lowerWeakOrigins[0] = {0.2,0};
-        upperWeakOrigins[0] = {0.2,width};
+        upperWeakOrigins[0] = {0.6,width};
 
         cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width);
         for (size_t i=1; i<this->bodyCount_-1; i++) {

src/multi-body-problem-data/bc.hh

@@ -8,6 +8,14 @@ class VelocityDirichletCondition
   void evaluate(double const &relativeTime, double &y) const {
     // Assumed to vanish at time zero
     double const finalVelocity = -5e-5;
+    
+    std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;
+    
+    if (relativeTime <= 0.1)
+        std::cout << "- loading phase" << std::endl;
+    else
+        std::cout << "- final velocity reached" << std::endl;
+    
     y = (relativeTime <= 0.1)
             ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
             : finalVelocity;

src/multi-body-problem.cc

@@ -138,13 +138,14 @@ int main(int argc, char *argv[]) {
      /* Vector def(levelContactNetwork.deformedGrids()[i]->size(dims));
       def = 1;
       deformedGridComplex.setDeformation(def, i);*/
-        auto& levelViews = levelContactNetwork.levelViews(i);
+
+        /*auto& levelViews = levelContactNetwork.levelViews(i);
 
         for (size_t j=0; j<levelViews.size(); j++) {
             writeToVTK(*levelViews[j], "", "body_" + std::to_string(i) + "_level_" + std::to_string(j));
         }
 
-        writeToVTK(levelContactNetwork.leafView(i), "", "body_" + std::to_string(i) + "_leaf");
+        writeToVTK(levelContactNetwork.leafView(i), "", "body_" + std::to_string(i) + "_leaf"); */
     }
 
     // ----------------------------
@@ -226,13 +227,12 @@ int main(int argc, char *argv[]) {
                         frictionBoundaries) //, weakPatches)
                   : nullptr;*/
 
-    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "body");
+    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
 
     IterationRegister iterationCount;
 
     auto const report = [&](bool initial = false) {
-      /*
-      if (writeData) {
+      /*if (writeData) {
         dataWriter->reportSolution(programState, levelContactNetwork.globalFriction());
         if (!initial)
           dataWriter->reportIterations(programState, iterationCount);
@@ -243,14 +243,14 @@ int main(int argc, char *argv[]) {
           programState.timeStep % restartSpacing == 0) {
         restartIO->write(programState);
         restartFile->flush();
-      }
+      }*/
 
       if (parset.get<bool>("io.printProgress"))
         std::cout << "timeStep = " << std::setw(6) << programState.timeStep
                   << ", time = " << std::setw(12) << programState.relativeTime
                   << ", tau = " << std::setw(12) << programState.relativeTau
                   << std::endl;
-      */
+
       if (parset.get<bool>("io.vtk.write")) {
         std::vector<ScalarVector> stress(bodyCount);
 
@@ -276,6 +276,20 @@ int main(int argc, char *argv[]) {
     BitVector totalDirichletNodes;
     levelContactNetwork.totalNodes("dirichlet", totalDirichletNodes);
 
+    /*for (size_t i=0; i<totalDirichletNodes.size(); i++) {
+        bool val = false;
+        for (size_t d=0; d<dims; d++) {
+            val = val || totalDirichletNodes[i][d];
+        }
+
+        totalDirichletNodes[i] = val;
+        for (size_t d=0; d<dims; d++) {
+            totalDirichletNodes[i][d] = val;
+        }
+    }*/
+
+    //print(totalDirichletNodes, "totalDirichletNodes:");
+
     using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
     using NonlinearFactory = SolverFactory<Functional, BitVector>;
 
@@ -316,14 +330,28 @@ int main(int argc, char *argv[]) {
 
     auto const mustRefine = [&](Updaters &coarseUpdater,
                                 Updaters &fineUpdater) {
+
+        //return false;
+      std::cout << "Step 1" << std::endl;
+
       std::vector<ScalarVector> coarseAlpha;
       coarseUpdater.state_->extractAlpha(coarseAlpha);
 
+      std::cout << "Step 2" << std::endl;
+
       std::vector<ScalarVector> fineAlpha;
       fineUpdater.state_->extractAlpha(fineAlpha);
 
+      std::cout << "Step 3" << std::endl;
+
       ScalarVector::field_type energyNorm = 0;
       for (size_t i=0; i<stateEnergyNorms.size(); i++) {
+          std::cout << "for " << i << std::endl;
+
+          std::cout << not stateEnergyNorms[i] << std::endl;
+
+          if (coarseAlpha[i].size()==0 || fineAlpha[i].size()==0)
+              continue;
           energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
       }
       return energyNorm > refinementTolerance;

src/multi-body-problem.cfg

@@ -3,11 +3,11 @@ gravity         = 9.81  # [m/s^2]
 
 [io]
 data.write      = false #true
-printProgress   = false
+printProgress   = true
 restarts.first  = 0
 restarts.spacing= 20
 restarts.write  = false #true
-vtk.write       = false
+vtk.write       = true
 
 [problem]
 finalTime       = 1000  # [s]
@@ -52,8 +52,8 @@ maximumIterations = 10000
 verbosity         = full
 
 [v.solver]
-maximumIterations = 100000
-verbosity         = quiet
+maximumIterations = 100
+verbosity         = full
 
 [v.fpi]
 maximumIterations = 10000

src/spatial-solving/fixedpointiterator.cc

@@ -64,6 +64,8 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
     Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
     std::vector<Vector>& velocityIterates) {
 
+  std::cout << "FixedPointIterator::run()" << std::endl;
+
   const auto nBodies = nBodyAssembler_.nGrids();
 
   std::vector<const Matrix*> matrices_ptr(nBodies);
@@ -75,9 +77,13 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
   Matrix bilinearForm;
   nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
 
+  print(bilinearForm, "bilinearForm:");
+
   Vector totalRhs;
   nBodyAssembler_.assembleRightHandSide(velocityRHSs, totalRhs);
 
+  print(totalRhs, "totalRhs:");
+
   Vector totalX;
   nBodyAssembler_.nodalToTransformed(velocityIterates, totalX);
 
@@ -96,11 +102,15 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
     }
   }
 
+  std::cout << "- Problem assembled: success" << std::endl;
+
   // set up functional and TNMMG solver
   Functional J(bilinearForm, totalRhs, globalFriction_, lower, upper);
   Factory solverFactory(parset_.sub("solver.tnnmg"), J, ignoreNodes_);
   auto step = solverFactory.step();
 
+  std::cout << "- Functional and TNNMG step setup: success" << std::endl;
+
   EnergyNorm<Matrix, Vector> energyNorm(bilinearForm);
   LoopSolver<Vector> velocityProblemSolver(step.get(), velocityMaxIterations_,
                                            velocityTolerance_, &energyNorm,
@@ -119,13 +129,22 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
     Vector totalVelocityIterate;
     nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
 
+    std::cout << "- FixedPointIteration iterate" << std::endl;
+
     // solve a velocity problem
     step->setProblem(totalVelocityIterate);
 
+    std::cout << "- Velocity problem set" << std::endl;
+
     velocityProblemSolver.preprocess();
+    std::cout << "-- Preprocessed" << std::endl;
     velocityProblemSolver.solve();
+    std::cout << "-- Solved" << std::endl;
     nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
 
+    //print(totalVelocityIterate, "totalVelocityIterate:");
+    //print(velocityIterates, "velocityIterates:");
+
     multigridIterations += velocityProblemSolver.getResult().iterations;
 
     std::vector<Vector> v_m;
@@ -140,15 +159,28 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
     std::vector<Vector> v_rel;
     relativeVelocities(totalVelocityIterate, v_rel);
 
+    //print(v_rel, "v_rel");
+
+    std::cout << "- State problem set" << std::endl;
+
     // solve a state problem
     updaters.state_->solve(v_rel);
+
+    std::cout << "-- Solved" << std::endl;
+
     std::vector<ScalarVector> newAlpha(nBodies);
     updaters.state_->extractAlpha(newAlpha);
 
     bool breakCriterion = true;
     for (size_t i=0; i<nBodies; i++) {
+        if (alpha[i].size()==0 || newAlpha[i].size()==0)
+            continue;
+
+        //print(alpha[i], "alpha i:");
+        //print(newAlpha[i], "new alpha i:");
         if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
             breakCriterion = false;
+            std::cout << "fixedPoint error: " << errorNorms_[i]->diff(alpha[i], newAlpha[i]) << std::endl;
             break;
         }
     }
@@ -160,6 +192,8 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
     alpha = newAlpha;
   }
 
+  std::cout << "-FixedPointIteration finished!" << std::endl;
+
   if (fixedPointIteration == fixedPointMaxIterations_)
     DUNE_THROW(Dune::Exception, "FPI failed to converge");
 

src/spatial-solving/solverfactory.cc

@@ -3,25 +3,29 @@
 #endif
 
 #include <dune/solvers/common/wrapownshare.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
 
 #include "solverfactory.hh"
 
+#include "../utils/debugutils.hh"
 
 template <class Functional, class BitVector>
 SolverFactory<Functional, BitVector>::SolverFactory(
     const Dune::ParameterTree& parset,
     Functional& J,
     const BitVector& ignoreNodes) :
-      J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))),
+      J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J)))
       //linear solver
       //preconditioner_(),
-      linearSolverStep_(),
-      energyNorm_(linearSolverStep_)
+      //linearSolverStep_(),
     {
     nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
 
     // linearSolverStep_.setPreconditioner(preconditioner_);
-    linearSolver_ = std::make_shared<LinearSolver>(linearSolverStep_, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), energyNorm_, Solver::QUIET);
+    linearSolverStep_ = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::gs(0.0, 0.0));
+    energyNorm_ = std::make_shared<EnergyNorm<Matrix, Vector>>(*linearSolverStep_);
+    linearSolver_ = std::make_shared<LinearSolver>(*linearSolverStep_, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), *energyNorm_, Solver::QUIET);
+    //linearSolver_ = std::make_shared<LinearSolver>();
 
     tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_, DefectProjection(), LineSearchSolver());
     tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));

src/spatial-solving/solverfactory.hh

@@ -8,13 +8,18 @@
 
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+//#include <dune/solvers/iterationsteps/cgstep.hh>
+//#include <dune/solvers/iterationsteps/amgstep.hh>
+//#include <dune/solvers/solvers/umfpacksolver.hh>
 
 #include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
 #include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
 #include <dune/tnnmg/projections/obstacledefectprojection.hh>
+#include <dune/tnnmg/linearsolvers/fastamgsolver.hh>
 
 //#include "tnnmg/functional.hh"
+//#include "tnnmg/tnnmgstep.hh"
 #include "tnnmg/linearization.hh"
 #include "tnnmg/linesearchsolver.hh"
 #include "tnnmg/localbisectionsolver.hh"
@@ -30,10 +35,14 @@ class SolverFactory {
     using LocalSolver = LocalBisectionSolver;
     using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, LocalSolver>;
 
-    using Linearization = typename Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<Functional, BitVector>;
+    using Linearization = Linearization<Functional, BitVector>;
     //using Preconditioner = Dune::Solvers::Preconditioner; //TODO Platzhalter für PatchPreconditioner
-    using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
+    //using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
+    //using LinearSolverStep = typename Dune::Solvers::AMGStep<Matrix, Vector>;
+    using LinearSolverStep = LinearIterationStep<Matrix, Vector>;
     using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    //using LinearSolver = typename Dune::Solvers::UMFPackSolver<Matrix, Vector>;
+    //using LinearSolver = typename Dune::TNNMG::FastAMGSolver<Matrix, Vector>;
 
     using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
 
@@ -57,8 +66,8 @@ class SolverFactory {
 
   // linear solver
   //Preconditioner preconditioner_;
-  LinearSolverStep linearSolverStep_;
-  EnergyNorm<Matrix, Vector> energyNorm_;
+  std::shared_ptr<LinearSolverStep> linearSolverStep_;
+  std::shared_ptr<EnergyNorm<Matrix, Vector>> energyNorm_;
   std::shared_ptr<LinearSolver> linearSolver_;
 
   // TNNMGStep

src/spatial-solving/tnnmg/functional.hh

@@ -97,6 +97,8 @@ class DirectionalRestriction :
   {
     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_;
     }

src/spatial-solving/tnnmg/linearization.hh

 #ifndef DUNE_TECTONIC_LINEARIZATION_HH
 #define DUNE_TECTONIC_LINEARIZATION_HH
 
-#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
-#include <dune/istl/matrixindexset.hh>
+#include <cstddef>
 
-//#include <>
+#include <dune/common/typetraits.hh>
+#include <dune/common/hybridutilities.hh>
 
-/*#include <dune/common/bitsetvector.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/fmatrixev.hh>
+#include <dune/matrix-vector/algorithm.hh>
 
-#include <dune/fufem/arithmetic.hh>
-#include <dune/solvers/common/interval.hh>
-#include <dune/solvers/computeenergy.hh>
-#include <dune/tnnmg/problem-classes/bisection.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineargsproblem.hh>
+#include <dune/istl/matrixindexset.hh>
 
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/minimisation.hh>
-#include <dune/tectonic/mydirectionalconvexfunction.hh>
-#include <dune/tectonic/quadraticenergy.hh>*/
+//#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
 
+#include "../../utils/debugutils.hh"
 
+/**
+ * derived from Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<F,BV>
+ */
 template<class F, class BV>
-class Linearization : public Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<F,BV> {
+class Linearization {
+public:
+  using Matrix = typename F::Matrix;
+  using Vector = typename F::Vector;
+  using BitVector = BV;
+  using ConstrainedMatrix = Matrix;
+  using ConstrainedVector = Vector;
+  using ConstrainedBitVector = BitVector;
+
 private:
-  using Base = Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<F,BV>;
-  using Vector = typename Base::Vector;
-  using BitVector = typename Base::BitVector;
+  using This = Linearization<F,BV>;
 
   template <class T>
-  void determineRegularityTruncation(const Vector& x, BitVector&& truncationFlags, const T& truncationTolerance)
+  void determineRegularityTruncation(const Vector& x, BitVector& truncationFlags, const T& truncationTolerance)
   {
     for (size_t i = 0; i < x.size(); ++i) {
-      if (this->f_.phi().regularity(i, x[i]) > truncationTolerance) {
+      if (f_.phi().regularity(i, x[i]) > truncationTolerance) {
         truncationFlags[i] = true;
       }
     }
   }
 
+  template<class NV, class NBV, class T>
+  static void determineTruncation(const NV& x, const NV& lower, const NV& upper, NBV&& truncationFlags, const T& truncationTolerance)
+  {
+    namespace H = Dune::Hybrid;
+    H::ifElse(Dune::IsNumber<NV>(), [&](auto id){
+      if ((id(x) <= id(lower)+truncationTolerance) || (id(x) >= id(upper) - truncationTolerance))
+        id(truncationFlags) = true;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(x))), [&](auto&& i) {
+        This::determineTruncation(x[i], lower[i], upper[i], truncationFlags[i], truncationTolerance);
+      });
+    });
+  }
+
+  template<class NV, class NBV>
+  static void truncateVector(NV& x, const NBV& truncationFlags)
+  {
+    namespace H = Dune::Hybrid;
+    H::ifElse(Dune::IsNumber<NV>(), [&](auto id){
+      if (id(truncationFlags))
+        id(x) = 0;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(x))), [&](auto&& i) {
+        This::truncateVector(x[i], truncationFlags[i]);
+      });
+    });
+  }
+
+  template<class NM, class RBV, class CBV>
+  static void truncateMatrix(NM& A, const RBV& rowTruncationFlags, const CBV& colTruncationFlags)
+  {
+    namespace H = Dune::Hybrid;
+    using namespace Dune::MatrixVector;
+    H::ifElse(Dune::IsNumber<NM>(), [&](auto id){
+      if(id(rowTruncationFlags) or id(colTruncationFlags))
+        A = 0;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(rowTruncationFlags))), [&](auto&& i) {
+        auto&& Ai = A[i];
+        sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
+            This::truncateMatrix(Aij, rowTruncationFlags[i], colTruncationFlags[j]);
+        });
+      });
+    });
+  }
+
 public:
   Linearization(const F& f, const BitVector& ignore) :
-    Base(f, ignore)
+      f_(f),
+      ignore_(ignore),
+      truncationTolerance_(1e-10)
   {}
 
   void bind(const Vector& x) {
-      auto&& A = this->f_.quadraticPart();
-      auto&& phi = this->f_.phi();
+      //std::cout << "Linearization::bind()" << std::endl;
+
+      auto&& A = f_.quadraticPart();
+      auto&& phi = f_.phi();
+
+      // determine which components to truncate
+      // ----------------
+      truncationFlags_ = ignore_;
+
+      determineTruncation(x, f_.lowerObstacle(), f_.upperObstacle(), truncationFlags_, truncationTolerance_); // obstacle truncation
+      determineRegularityTruncation(x, truncationFlags_, truncationTolerance_); // truncation due to regularity deficit of nonlinearity
+
 
       // compute hessian
       // ---------------
@@ -57,43 +116,70 @@ class Linearization : public Dune::TNNMG::BoxConstrainedQuadraticFunctionalConst
       phi.addHessianIndices(indices);
 
       // construct matrix from pattern and initialize it
-      indices.exportIdx(this->hessian_);
-      this->hessian_ = 0.0;
+      indices.exportIdx(hessian_);
+      hessian_ = 0.0;
 
       // quadratic part
       for (size_t i = 0; i < A.N(); ++i) {
         auto const end = std::end(A[i]);
         for (auto it = std::begin(A[i]); it != end; ++it)
-          this->hessian_[i][it.index()] += *it;
+          hessian_[i][it.index()] += *it;
       }
 
+      //print(hessian_, "Hessian:");
+      //std::cout << "--------" << (double) hessian_[21][21] << "------------" << std::endl;
+
       // nonlinearity part
-      phi.addHessian(x, this->hessian_);
+      phi.addHessian(x, hessian_);
 
       // compute gradient
       // ----------------
 
       // quadratic part
-      this->negativeGradient_ = derivative(this->f_)(x); // A*x - b
+      negativeGradient_ = derivative(f_)(x); // A*x - b
 
       // nonlinearity part
-      phi.addGradient(x, this->negativeGradient_);
+      phi.addGradient(x, negativeGradient_);
 
       // -grad is needed for Newton step
-      this->negativeGradient_ *= -1;
+      negativeGradient_ *= -1;
 
 
-      // determine which components to truncate
-      // ----------------
-      this->truncationFlags_ = this->ignore_;
-
-      determineTruncation(x, this->f_.lowerObstacle(), this->f_.upperObstacle(), this->truncationFlags_, this->truncationTolerance_); // obstacle truncation
-      determineRegularityTruncation(x, this->truncationFlags_, this->truncationTolerance_); // truncation due to regularity deficit of nonlinearity
 
       // truncate gradient and hessian
-      truncateVector(this->negativeGradient_, this->truncationFlags_);
-      truncateMatrix(this->hessian_, this->truncationFlags_, this->truncationFlags_);
+      truncateVector(negativeGradient_, truncationFlags_);
+      truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
+
+      //print(hessian_, "Hessian:");
+  }
+
+
+  void extendCorrection(ConstrainedVector& cv, Vector& v) const {
+    v = cv;
+    truncateVector(v, truncationFlags_);
+  }
+
+  const BitVector& truncated() const {
+    return truncationFlags_;
+  }
+
+  const auto& negativeGradient() const {
+    return negativeGradient_;
+  }
+
+  const auto& hessian() const {
+    return hessian_;
   }
+
+private:
+  const F& f_;
+  const BitVector& ignore_;
+
+  double truncationTolerance_;
+
+  Vector negativeGradient_;
+  Matrix hessian_;
+  BitVector truncationFlags_;
 };
 
 #endif

src/spatial-solving/tnnmg/linesearchsolver.hh

@@ -5,6 +5,7 @@
 
 #include <dune/tnnmg/problem-classes/bisection.hh>
 
+#include "../../utils/almostequal.hh"
 
 /**
  * \brief A local solver for scalar quadratic obstacle problems with nonlinearity
@@ -22,10 +23,17 @@ class LineSearchSolver
     if (D[1] > 0) // NOTE: Numerical instability can actually get us here
         return;
 
+    if (almost_equal(f.domain()[0], f.domain()[1], 2)) {
+        x = f.domain()[0];
+        return;
+    }
+
     if (not ignore) {
         int bisectionsteps = 0;
         const Bisection globalBisection;
+        std::cout << "LineSearchSolver: starting bisection" << std::endl;
         x = globalBisection.minimize(f, 1.0, 0.0, bisectionsteps);
+        std::cout << "LineSearchSolver: bisection terminated" << std::endl;
         // x = globalBisection.minimize(f, vNorm, 0.0, bisectionsteps) / vNorm;  // used rescaled v in f?
     }
   }

src/spatial-solving/tnnmg/tnnmgstep.hh

+#ifndef DUNE_TNNMG_ITERATIONSTEPS_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
+#define DUNE_TNNMG_ITERATIONSTEPS_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
+
+#include <string>
+#include <sstream>
+#include <vector>
+#include <iomanip>
+
+#include <dune/common/timer.hh>
+
+#include <dune/solvers/common/resize.hh>
+#include "dune/solvers/iterationsteps/iterationstep.hh"
+#include "dune/solvers/iterationsteps/lineariterationstep.hh"
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/solvers/linearsolver.hh>
+
+#include <dune/tectonic/../../src/utils/debugutils.hh>
+
+
+/**
+ * \brief One iteration of the TNNMG method
+ *
+ * \tparam F Functional to minimize
+ * \tparam BV Bit-vector type for marking ignored components
+ */
+template<class F, class BV, class Linearization,
+                                  class DefectProjection,
+                                  class LineSearchSolver>
+class TNNMGStep :
+  public Dune::Solvers::IterationStep<typename F::Vector, BV>
+{
+  using Base = Dune::Solvers::IterationStep<typename F::Vector, BV>;
+
+public:
+
+  using Vector = typename F::Vector;
+  using ConstrainedVector = typename Linearization::ConstrainedVector;
+  using ConstrainedMatrix = typename Linearization::ConstrainedMatrix;
+  using BitVector = typename Base::BitVector;
+  using ConstrainedBitVector = typename Linearization::ConstrainedBitVector;
+  using Functional = F;
+  using IterativeSolver = Dune::Solvers::IterativeSolver< ConstrainedVector, Dune::Solvers::DefaultBitVector_t<ConstrainedVector> >;
+  using LinearSolver = Dune::Solvers::LinearSolver< ConstrainedMatrix,  ConstrainedVector >;
+
+  /** \brief Constructor with an iterative solver object for the linear correction
+   * \param iterativeSolver This is a callback used to solve the constrained linearized system
+   * \param projection This is a callback used to compute a projection into a defect-admissible set
+   * \param lineSolver This is a callback used to minimize a directional restriction of the functional
+   *        for computing a damping parameter
+   */
+  TNNMGStep(const Functional& f,
+            Vector& x,
+            std::shared_ptr<IterationStep<Vector,BitVector> > nonlinearSmoother,
+            std::shared_ptr<IterativeSolver> iterativeSolver,
+            const DefectProjection& projection,
+            const LineSearchSolver& lineSolver)
+  : Base(x),
+    f_(&f),
+    nonlinearSmoother_(nonlinearSmoother),
+    projection_(projection),
+    lineSolver_(lineSolver)
+  {}
+
+  /** \brief Constructor with a linear solver object for the linear correction
+   * \param linearSolver This is a callback used to solve the constrained linearized system
+   * \param projection This is a callback used to compute a projection into a defect-admissible set
+   * \param lineSolver This is a callback used to minimize a directional restriction of the functional
+   *        for computing a damping parameter
+   */
+  TNNMGStep(const Functional& f,
+            Vector& x,
+            std::shared_ptr<IterationStep<Vector,BitVector> > nonlinearSmoother,
+            std::shared_ptr<LinearSolver> linearSolver,
+            const DefectProjection& projection,
+            const LineSearchSolver& lineSolver)
+  : Base(x),
+    f_(&f),
+    nonlinearSmoother_(nonlinearSmoother),
+    projection_(projection),
+    lineSolver_(lineSolver)
+  {}
+
+  /** \brief Constructor with a LinearIterationStep object for the linear correction
+   * \param linearIterationStep This is a callback used to solve the constrained linearized system
+   * \param projection This is a callback used to compute a projection into a defect-admissible set
+   * \param lineSolver This is a callback used to minimize a directional restriction of the functional
+   *        for computing a damping parameter
+   */
+  TNNMGStep(const Functional& f,
+            Vector& x,
+            std::shared_ptr<Dune::Solvers::IterationStep<Vector,BitVector> > nonlinearSmoother,
+            std::shared_ptr<Dune::Solvers::LinearIterationStep<ConstrainedMatrix,ConstrainedVector> > linearIterationStep,
+            unsigned int noOfLinearIterationSteps,
+            const DefectProjection& projection,
+            const LineSearchSolver& lineSolver)
+  : Base(x),
+    f_(&f),
+    nonlinearSmoother_(nonlinearSmoother),
+    projection_(projection),
+    lineSolver_(lineSolver)
+  {}
+
+  using Base::getIterate;
+
+  void preprocess() override
+  {
+    nonlinearSmoother_->setIgnore(this->ignore());
+    nonlinearSmoother_->preprocess();
+  }
+
+  void setPreSmoothingSteps(std::size_t i)
+  {
+    preSmoothingSteps_ = i;
+  }
+
+  /**
+   * \brief Do one TNNMG step
+   */
+  void iterate() override
+  {
+
+    const auto& f = *f_;
+    const auto& ignore = (*this->ignoreNodes_);
+    auto& x = *getIterate();
+
+    print(f.quadraticPart(), "f.quadraticPart():");
+
+    print(f.linearPart(), "f.linearPart():");
+
+    // Nonlinear presmoothing
+    for (std::size_t i=0; i<preSmoothingSteps_; ++i)
+        nonlinearSmoother_->iterate();
+
+    print(x, "TNNMG iterate after smoothing:");
+
+    // Compute constraint/truncated linearization
+    if (not linearization_)
+      linearization_ = std::make_shared<Linearization>(f, ignore);
+
+    linearization_->bind(x);
+
+    auto&& A = linearization_->hessian();
+    auto&& r = linearization_->negativeGradient();
+
+    print(A, "TNNMG Linear problem matrix:");
+    print(r, "TNNMG Linear problem rhs:");
+
+    // Compute inexact solution of the linearized problem
+    Solvers::resizeInitializeZero(correction_, x);
+    Solvers::resizeInitializeZero(constrainedCorrection_, r);
+
+    linearCorrection_(A, constrainedCorrection_, r);
+
+    print(constrainedCorrection_, "contrained correction:");
+
+    linearization_->extendCorrection(constrainedCorrection_, correction_);
+
+    print(correction_, "correction:");
+
+    // Project onto admissible set
+    projection_(f, x, correction_);
+
+    // Line search
+    auto fv = directionalRestriction(f, x, correction_);
+    dampingFactor_ = 0;
+    lineSolver_(dampingFactor_, fv, false);
+    if (std::isnan(dampingFactor_))
+      dampingFactor_ = 0;
+    correction_ *= dampingFactor_;
+
+    x += correction_;
+  }
+
+  /**
+   * \brief Export the last computed damping factor
+   */
+  double lastDampingFactor() const
+  {
+    return dampingFactor_;
+  }
+
+  /**
+   * \brief Export the last used linearization
+   */
+  const Linearization& linearization() const
+  {
+    return *linearization_;
+  }
+
+private:
+
+  const Functional* f_;
+
+  std::shared_ptr<IterationStep<Vector,BitVector> > nonlinearSmoother_;
+  std::size_t preSmoothingSteps_ = 1;
+
+  std::shared_ptr<Linearization> linearization_;
+
+  typename Linearization::ConstrainedVector constrainedCorrection_;
+  Vector correction_;
+  DefectProjection projection_;
+  LineSearchSolver lineSolver_;
+  double dampingFactor_;
+};
+
+#endif

src/spatial-solving/tnnmg/zerononlinearity.hh

@@ -34,7 +34,7 @@ class ZeroNonlinearity
     void addHessianIndices(IndexSet& indices) const {}
 
     /** \brief Returns the interval \f$ [0,0]\f$ */
-    void subDiff(int i, double x, Dune::Solvers::Interval<double>& D, int j) const
+    void subDiff(int i, double x, Dune::Solvers::Interval<double>& D) const
     {
         D[0] = 0.0;
         D[1] = 0.0;
@@ -48,12 +48,13 @@ class ZeroNonlinearity
     }
 
     /** \brief Returns 0 */
-    double regularity(int i, double x, int j) const
+    template <class VectorType>
+    double regularity(int i, const VectorType& x) const
     {
         return 0.0;
     }
 
-    void domain(int i, Dune::Solvers::Interval<double>& dom, int j) const
+    void domain(int i, Dune::Solvers::Interval<double>& dom) const
     {
         dom[0] = -std::numeric_limits<double>::max();
         dom[1] = std::numeric_limits<double>::max();

src/time-stepping/adaptivetimestepper.cc

@@ -59,16 +59,24 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
     refine, we compare the result of (F1+F2) with R1, i.e. two steps
     of size tau/2 with one of size tau. The method makes multiple
     coarsening/refining attempts, with coarsening coming first. */
+
+  std::cout << "AdaptiveTimeStepper::advance()" << std::endl;
+
   if (R1_.updaters == Updaters())
     R1_ = step(current_, relativeTime_, relativeTau_);
 
+  std::cout << "AdaptiveTimeStepper Step 1" << std::endl;
+
+
   bool didCoarsen = false;
   iterationRegister_.reset();
   UpdatersWithCount R2;
   UpdatersWithCount C;
   while (relativeTime_ + relativeTau_ <= 1.0) {
     R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
+    std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
     C = step(current_, relativeTime_, 2 * relativeTau_);
+    std::cout << "AdaptiveTimeStepper C computed!" << std::endl << std::endl;
     if (mustRefine_(R2.updaters, C.updaters))
       break;
 
@@ -76,15 +84,21 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
     relativeTau_ *= 2;
     R1_ = C;
   }
+
+  std::cout << "AdaptiveTimeStepper Step 1" << std::endl;
   UpdatersWithCount F1;
   UpdatersWithCount F2;
   if (!didCoarsen) {
     while (true) {
       F1 = step(current_, relativeTime_, relativeTau_ / 2.0);
+      std::cout << "AdaptiveTimeStepper F1 computed!" << std::endl << std::endl;
       F2 = step(F1.updaters, relativeTime_ + relativeTau_ / 2.0,
                 relativeTau_ / 2.0);
-      if (!mustRefine_(F2.updaters, R1_.updaters))
+      std::cout << "AdaptiveTimeStepper F2 computed!" << std::endl << std::endl;
+      if (!mustRefine_(F2.updaters, R1_.updaters)) {
+        std::cout << "Sufficiently refined!" << std::endl;
         break;
+      }
 
       relativeTau_ /= 2.0;
       R1_ = F1;
@@ -92,11 +106,18 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
     }
   }
 
+  std::cout << "AdaptiveTimeStepper::advance() ...";
+
   iterationRegister_.registerFinalCount(R1_.count);
   relativeTime_ += relativeTau_;
   current_ = R1_.updaters;
+
+  //UpdatersWithCount emptyR1;
+  //R1_ = emptyR1;
   R1_ = R2;
 
+  std::cout << " done" << std::endl;
+
   return iterationRegister_;
 }
 

src/time-stepping/coupledtimestepper.cc

@@ -28,6 +28,9 @@ template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
 FixedPointIterationCounter
 CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(double relativeTime,
                                                        double relativeTau) {
+
+  std::cout << "CoupledTimeStepper::step()" << std::endl;
+
   updaters_.state_->nextTimeStep();
   updaters_.rate_->nextTimeStep();