add solverfactorytest

src/CMakeLists.txt

@@ -64,17 +64,42 @@ set(UGW_SOURCE_FILES
   one-body-problem-data/mygrid.cc
 )
 
+set(SFT_SOURCE_FILES
+  assemblers.cc
+  data-structures/body.cc
+  data-structures/levelcontactnetwork.cc
+  data-structures/contactnetwork.cc
+  data-structures/enumparser.cc
+  factories/stackedblocksfactory.cc
+  io/vtk.cc
+  multi-body-problem-data/grid/cuboidgeometry.cc
+  multi-body-problem-data/grid/mygrids.cc
+  multi-body-problem-data/grid/simplexmanager.cc
+  #spatial-solving/solverfactory.cc
+  spatial-solving/fixedpointiterator.cc
+  #spatial-solving/solverfactory_old.cc
+  time-stepping/adaptivetimestepper.cc
+  time-stepping/coupledtimestepper.cc
+  time-stepping/rate.cc
+  time-stepping/rate/rateupdater.cc
+  time-stepping/state.cc
+  solverfactorytest.cc
+)
+
 foreach(_dim 2 3)
   set(_sw_target one-body-problem-${_dim}D)
   set(_msw_target multi-body-problem-${_dim}D)
   set(_ugw_target uniform-grid-writer-${_dim}D)
+  set(_sft_target solverfactorytest-${_dim}D)
 
   add_executable(${_sw_target} ${SW_SOURCE_FILES})
   add_executable(${_msw_target} ${MSW_SOURCE_FILES})
   add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
+  add_executable(${_sft_target} ${SFT_SOURCE_FILES})
   add_dune_ug_flags(${_sw_target})
   add_dune_ug_flags(${_msw_target})
   add_dune_ug_flags(${_ugw_target})
+  add_dune_ug_flags(${_sft_target})
 
   add_dune_hdf5_flags(${_sw_target})
   add_dune_hdf5_flags(${_msw_target})
@@ -82,4 +107,5 @@ 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 ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  set_property(TARGET ${_sft_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
 endforeach()

src/data-structures/program_state.hh

@@ -243,6 +243,8 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       solver.preprocess();
       solver.solve();
 
+      std::cout << "ProgramState: Energy of TNNMG solution: " << J(tnnmgStep->getSol()) << std::endl;
+
       nBodyAssembler.postprocess(tnnmgStep->getSol(), _x);
 
       Vector res = totalRhs;

src/factories/stackedblocksfactory.cc

@@ -73,10 +73,13 @@ void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {
 
         for (size_t i=1; i<this->bodyCount_-1; i++) {
             origins[i] = cuboidGeometries_[i-1]->upperLeft();
+            auto height_i = height/3.0;
+
             GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
-            GlobalCoords upperWeakOrigin_i = upperWeakOrigin + origins[i];
+            GlobalCoords upperWeakOrigin_i = {upperWeakOrigin[0], height_i};
+            upperWeakOrigin_i += origins[i];
 
-            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height);
+            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height_i);
 
 
             cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));

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

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

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

@@ -7,7 +7,6 @@ class VelocityDirichletCondition
     : public Dune::VirtualFunction<double, double> {
   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;
@@ -20,8 +19,6 @@ class VelocityDirichletCondition
     y = (relativeTime <= 0.1)
             ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
             : finalVelocity;
-            */
-      y = relativeTime;
   }
 };
 

src/obstaclesolver.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_OBSTACLESOLVER_HH
+#define DUNE_TECTONIC_OBSTACLESOLVER_HH
+
+#include "spatial-solving/tnnmg/localproblem.hh"
+
+
+/**
+ * \brief A local solver for quadratic obstacle problems
+ *
+ * \todo Add concept check for the function interface
+ */
+class ObstacleSolver
+{
+public:
+  template<class Vector, class Functional, class BitVector>
+  constexpr void operator()(Vector& x, const Functional& f, const BitVector& ignore) const
+  {
+      auto& A = f.quadraticPart();
+      auto& r = f.linearPart();
+
+      LocalProblem<std::decay_t<decltype(A)> , std::decay_t<decltype(r)>> localProblem(A, r, ignore);
+      Vector newR;
+      localProblem.getLocalRhs(x, newR);
+
+      auto directSolver = Dune::Solvers::UMFPackSolver<std::decay_t<decltype(A)>, std::decay_t<decltype(r)>>();
+
+      directSolver->setProblem(localProblem.getMat(), x, newR);
+      directSolver->preprocess();
+      directSolver->solve();
+
+      const auto& lower = f.lowerObstacle();
+      const auto& upper = f.upperObstacle();
+
+      for (size_t i=0; i<x.size(); i++) {
+        if (ignore[i][0])
+            continue;
+
+        if (x[i] < lower[i])
+            x[i] = lower[i];
+        if (x[i] > upper[i])
+            x[i] = upper[i];
+      }
+  }
+};
+
+#endif

src/spatial-solving/fixedpointiterator.cc

@@ -210,6 +210,9 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
     std::cout << "-- Solved" << std::endl;
 
     const auto& tnnmgSol = step->getSol();
+
+    std::cout << "FixPointIterator: Energy of TNNMG solution: " << J(tnnmgSol) << std::endl;
+
     nBodyAssembler.postprocess(tnnmgSol, velocityIterates);
     //nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
 

src/spatial-solving/solverfactory.cc

@@ -19,7 +19,14 @@ SolverFactory<Functional, BitVector>::SolverFactory(
     const BitVector& ignoreNodes) :
         J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))) {
 
-    nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
+    auto localSolver = Dune::TNNMG::gaussSeidelLocalSolver(LocalSolver());
+    //nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, localSolver);
+
+    //using MyNonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, decltype(localSolver), BitVector>;
+    nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, localSolver);
+
+    //nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
+
 
     auto linearSolver_ptr = Dune::Solvers::wrap_own_share<std::decay_t<LinearSolver>>(std::forward<LinearSolver>(linearSolver));
 

src/spatial-solving/solverfactory.hh

@@ -10,11 +10,12 @@
 #include <dune/solvers/solvers/loopsolver.hh>
 #include <dune/solvers/iterationsteps/cgstep.hh>
 
-//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
 #include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
 #include <dune/tnnmg/projections/obstacledefectprojection.hh>
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
 
-#include "tnnmg/tnnmgstep.hh"
+//#include "tnnmg/tnnmgstep.hh"
 #include "tnnmg/linearization.hh"
 #include "tnnmg/linesearchsolver.hh"
 #include "tnnmg/localbisectionsolver.hh"
@@ -27,13 +28,13 @@ class SolverFactory {
     using Vector = typename Functional::Vector;
     using BitVector = BitVectorType;
 
-    using LocalSolver = LocalBisectionSolver;
-    using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, LocalSolver>;
+    using LocalSolver = Dune::TNNMG::ScalarObstacleSolver;//LocalBisectionSolver;
+    using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, Dune::TNNMG::GaussSeidelLocalSolver<LocalSolver>, BitVector>;
     using Linearization = Linearization<Functional, BitVector>;
     using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
 
-    using Step = typename Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
-
+    //using Step = TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
+    using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
 
   template <class LinearSolver>
   SolverFactory(const Dune::ParameterTree&,

src/spatial-solving/tnnmg/functional.hh

@@ -16,6 +16,12 @@
 
 #include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
 
+template<class M, class V, class N, 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
  *
@@ -346,12 +352,13 @@ auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, co
 
   auto&& phii = f.phi().restriction(i);
 
-  auto v = ri;
+  /*auto v = ri;
   double const vnorm = v.two_norm();
   if (vnorm > 1.0)
-        v /= vnorm;
+        v /= vnorm;*/
 
-  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 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));
 }
 
 
@@ -389,8 +396,8 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
           MM&& matrix,
           VV&& linearPart,
           NN&& phi,
-          LL& lower,
-          UU& upper) :
+          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))
   {}

src/spatial-solving/tnnmg/linearcorrection.hh

-#ifndef DUNE_TNNMG_ITERATIONSTEPS_LINEARCORRECTION_HH
-#define DUNE_TNNMG_ITERATIONSTEPS_LINEARCORRECTION_HH 1
+#ifndef DUNE_TECTONIC_ITERATIONSTEPS_LINEARCORRECTION_HH
+#define DUNE_TECTONIC_ITERATIONSTEPS_LINEARCORRECTION_HH
 
 #include <functional>
 #include <memory>
@@ -38,7 +38,7 @@ emptyIgnore(const Vector& v)
 
 template<typename Matrix, typename Vector>
 LinearCorrection<Matrix, Vector>
-makeLinearCorrection(std::shared_ptr< Dune::Solvers::LinearSolver<Matrix, Vector> > linearSolver)
+buildLinearCorrection(std::shared_ptr< Dune::Solvers::LinearSolver<Matrix, Vector> > linearSolver)
 {
   return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
 
@@ -56,7 +56,7 @@ makeLinearCorrection(std::shared_ptr< Dune::Solvers::LinearSolver<Matrix, Vector
 
 template<typename Matrix, typename Vector>
 LinearCorrection<Matrix, Vector>
-makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > iterativeSolver)
+buildLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > iterativeSolver)
 {
   return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
 
@@ -76,16 +76,16 @@ makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > i
       directSolver->preprocess();
       directSolver->solve();
 
-      x = directX;
-    /*using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
+      //x = directX;
+    using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
 
     auto linearIterationStep = dynamic_cast<LinearIterationStep*>(&iterativeSolver->getIterationStep());
     if (not linearIterationStep)
       DUNE_THROW(Dune::Exception, "iterative solver must use a linear iteration step");
 
     auto empty = emptyIgnore(x);
-
     linearIterationStep->setIgnore(empty);
+    //linearIterationStep->setIgnore(ignore);
     linearIterationStep->setProblem(A, x, b);
     iterativeSolver->preprocess();
     iterativeSolver->solve();
@@ -94,13 +94,13 @@ makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > i
     const auto& norm = iterativeSolver->getErrorNorm();
     auto error = norm.diff(linearIterationStep->getSol(), directX);
 
-    std::cout << "Linear solver iterations: " << iterativeSolver->getResult().iterations << " Error: " << error << std::endl;*/
+    std::cout << "Linear solver iterations: " << iterativeSolver->getResult().iterations << " Error: " << error << std::endl;
   };
 }
 
 template<typename Matrix, typename Vector>
 LinearCorrection<Matrix, Vector>
-makeLinearCorrection(std::shared_ptr< Dune::Solvers::LinearIterationStep<Matrix, Vector> > linearIterationStep, int nIterationSteps = 1)
+buildLinearCorrection(std::shared_ptr< Dune::Solvers::LinearIterationStep<Matrix, Vector> > linearIterationStep, int nIterationSteps = 1)
 {
   return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
     linearIterationStep->setIgnore(ignore);

src/spatial-solving/tnnmg/tnnmgstep.hh

-#ifndef DUNE_TNNMG_ITERATIONSTEPS_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
-#define DUNE_TNNMG_ITERATIONSTEPS_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
+#ifndef DUNE_TECTONIC_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
+#define DUNE_TECTONIC_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
 
 #include <string>
 #include <sstream>
@@ -24,8 +24,6 @@
 
 #include <dune/tectonic/../../src/utils/debugutils.hh>
 
-namespace Dune {
-namespace TNNMG {
 
 /**
  * \brief One iteration of the TNNMG method
@@ -49,8 +47,8 @@ class TNNMGStep :
   using BitVector = typename Base::BitVector;
   using ConstrainedBitVector = typename Linearization::ConstrainedBitVector;
   using Functional = F;
-  using IterativeSolver = Solvers::IterativeSolver< ConstrainedVector, Solvers::DefaultBitVector_t<ConstrainedVector> >;
-  using LinearSolver = Solvers::LinearSolver< ConstrainedMatrix,  ConstrainedVector >;
+  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
@@ -67,7 +65,7 @@ class TNNMGStep :
   : Base(x),
     f_(&f),
     nonlinearSmoother_(nonlinearSmoother),
-    linearCorrection_(makeLinearCorrection<ConstrainedMatrix>(iterativeSolver)),
+    linearCorrection_(buildLinearCorrection<ConstrainedMatrix>(iterativeSolver)),
     projection_(projection),
     lineSolver_(lineSolver)
   {}
@@ -87,7 +85,7 @@ class TNNMGStep :
   : Base(x),
     f_(&f),
     nonlinearSmoother_(nonlinearSmoother),
-    linearCorrection_(makeLinearCorrection(linearSolver)),
+    linearCorrection_(buildLinearCorrection(linearSolver)),
     projection_(projection),
     lineSolver_(lineSolver)
   {}
@@ -100,15 +98,15 @@ class TNNMGStep :
    */
   TNNMGStep(const Functional& f,
             Vector& x,
-            std::shared_ptr<Solvers::IterationStep<Vector,BitVector> > nonlinearSmoother,
-            std::shared_ptr<Solvers::LinearIterationStep<ConstrainedMatrix,ConstrainedVector> > linearIterationStep,
+            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),
-    linearCorrection_(makeLinearCorrection(linearIterationStep, noOfLinearIterationSteps)),
+    linearCorrection_(buildLinearCorrection(linearIterationStep, noOfLinearIterationSteps)),
     projection_(projection),
     lineSolver_(lineSolver)
   {}
@@ -145,24 +143,24 @@ class TNNMGStep :
     const auto& ignore = (*this->ignoreNodes_);
     auto& x = *getIterate();
 
-    std::cout << "TNNMGStep::iterate " << std::endl;
+    //std::cout << "TNNMGStep::iterate " << std::endl;
 
     //print(f.quadraticPart(), "f.quadraticPart():");
 
     //print(f.linearPart(), "f.linearPart():");
 
-    std::cout << "-- energy before smoothing: " << f(x) << std::endl;
+    //std::cout << "-- energy before smoothing: " << f(x) << std::endl;
 
-    print(x, "TNNMG iterate after smoothing:");
+    //print(x, "TNNMG iterate after smoothing:");
 
     // Nonlinear presmoothing
     for (std::size_t i=0; i<preSmoothingSteps_; ++i)
         nonlinearSmoother_->iterate();
 
     //std::cout << "- nonlinear presmoothing: success" << std::endl;
-    print(x, "TNNMG iterate after smoothing:");
+    //print(x, "TNNMG iterate after smoothing:");
 
-    std::cout << "-- energy after presmoothing: " << f(x) << std::endl;
+    //std::cout << "-- energy after presmoothing: " << f(x) << std::endl;
 
     // Compute constraint/truncated linearization
     if (not linearization_)
@@ -173,20 +171,20 @@ class TNNMGStep :
     auto&& A = linearization_->hessian();
     auto&& r = linearization_->negativeGradient();
 
-    print(A, "TNNMG Linear problem matrix:");
+    /*print(A, "TNNMG Linear problem matrix:");
     print(r, "TNNMG Linear problem rhs:");
     print(ignore, "ignore:");
-    print(linearization_->truncated(), "truncation:");
+    print(linearization_->truncated(), "truncation:");*/
 
     // Compute inexact solution of the linearized problem
-    Solvers::resizeInitializeZero(correction_, x);
-    Solvers::resizeInitializeZero(constrainedCorrection_, r);
+    Dune::Solvers::resizeInitializeZero(correction_, x);
+    Dune::Solvers::resizeInitializeZero(constrainedCorrection_, r);
 
 
 
     //print(constrainedCorrection_compare, "direct solver solution:");
 
-    std::cout << "- computing linear correction..." << std::endl;
+    //std::cout << "- computing linear correction..." << std::endl;
 
     //linearCorrection_(A, constrainedCorrection_, r, ignore);
 
@@ -207,9 +205,13 @@ class TNNMGStep :
     //std::cout << "- extended correction: success" << std::endl;
     //print(correction_, "correction:");
 
+    //std::cout << "-- energy before projection: " << f(x) << std::endl;
+
     // Project onto admissible set
     projection_(f, x, correction_);
 
+    //std::cout << "-- energy after projection: " << f(x) << std::endl;
+
     //std::cout << "- projection onto admissible set: success" << std::endl;
 
     //print(correction_, "correction:");
@@ -231,6 +233,8 @@ class TNNMGStep :
     dampingFactor_ = 0;
     lineSolver_(dampingFactor_, fv, false);
 
+    //std::cout << "damping factor: " << dampingFactor_ << std::endl;
+
    // std::cout << "- computing damping factor: success" << std::endl;
     if (std::isnan(dampingFactor_))
       dampingFactor_ = 0;
@@ -240,7 +244,7 @@ class TNNMGStep :
 
     x += correction_;
 
-    std::cout << "-- energy after correction: " << energy(x) << std::endl;
+    //std::cout << "-- energy after correction: " << energy(x) << std::endl;
   }
 
   /**
@@ -279,7 +283,4 @@ class TNNMGStep :
 };
 
 
-} // end namespace TNNMG
-} // end namespace Dune
-
 #endif