cleanup for updated dune-* modules

dune/amdis/AMDiS.cpp

@@ -5,8 +5,6 @@
 
 #include <boost/program_options.hpp>
 
-#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
-
 // AMDiS includes
 #include "Initfile.hpp"
 #include "Log.hpp"
@@ -15,12 +13,10 @@ namespace AMDiS
 {
 //   using namespace std;
 
-  void init(int& argc, char**& argv, std::string initFileName)
+  Dune::MPIHelper& init(int& argc, char**& argv, std::string initFileName)
   {
     // Maybe initialize MPI
-    std::cout << "call MPIHelper::instance...";
-    Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
-    std::cout << " [ ok ]\n";
+    Dune::MPIHelper& mpiHelper = Dune::MPIHelper::instance(argc, argv);
     
     Parameters::clearData();
 
@@ -82,6 +78,8 @@ namespace AMDiS
     // if (vm.count("parameters") && !ignoreCommandline)
     //   Parameters::readArgv(vm["parameters"].as<std::string>(),0);
     // TODO: add command-line arguments again.
+    
+    return mpiHelper;
   }
 
 

dune/amdis/AMDiS.hpp

@@ -7,10 +7,11 @@
 #include <string>
 
 #include <dune/common/exceptions.hh> // We use exceptions
+#include <dune/common/parallel/mpihelper.hh>
 
 namespace AMDiS
 {
-  void init(int& argc, char**& argv, std::string initFileName = "");
+  Dune::MPIHelper& init(int& argc, char**& argv, std::string initFileName = "");
 
   void init(std::string initFileName);
 

dune/amdis/DirichletBC.inc.hpp

 #pragma once
 
 #include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/amdis/LinearAlgebra.hpp>
 
 namespace AMDiS
 {
@@ -14,7 +15,7 @@ namespace AMDiS
     using Dune::Functions::interpolate;
     
     if (!initialized) {
-//       interpolate(matrix.getRowFeSpace(), dirichletNodes, predicate);
+      interpolate(matrix.getRowFeSpace(), dirichletNodes, predicate);
       initialized = true;
     }
   }
@@ -33,8 +34,8 @@ namespace AMDiS
     matrix.clearDirichletRows(dirichletNodes, apply);
 
     if (apply) {
-//       interpolate(matrix.getRowFeSpace(), rhs.getVector(), values, dirichletNodes);
-//       interpolate(matrix.getColFeSpace(), solution.getVector(), values, dirichletNodes);
+      interpolate(matrix.getRowFeSpace(), wrapper(rhs.getVector()), values, dirichletNodes);
+      interpolate(matrix.getColFeSpace(), wrapper(solution.getVector()), values, dirichletNodes);
     }
   }
   

dune/amdis/Log.hpp

@@ -39,7 +39,7 @@
 #ifdef NDEBUG
   #define AMDIS_FUNCNAME_DBG(nn)
 #else
-  #define AMDIS_FUNCNAME_DBG(nn) const char *funcName; funcName = nn;
+  #define AMDIS_FUNCNAME_DBG(nn) AMDIS_UNUSED(const char *funcName); funcName = nn;
 #endif
 
 

dune/amdis/Operator.inc.hpp

@@ -388,7 +388,7 @@ namespace AMDiS
     auto geometry = element.geometry();
     
     const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
-    const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
+//     const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
     
     int order = getQuadratureDegree(2);
     const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);

dune/amdis/ProblemInstat.hpp

@@ -69,7 +69,7 @@ namespace AMDiS
     }
 
     /// Implementation of \ref ProblemTimeInterface::transferInitialSolution().
-    virtual void transferInitialSolution(AdaptInfo& adaptInfo);
+    virtual void transferInitialSolution(AdaptInfo& adaptInfo) override;
 
   protected:
     /// Space problem solved in each timestep.

dune/amdis/ProblemInstatBase.hpp

@@ -73,7 +73,7 @@ namespace AMDiS
     virtual void closeTimestep(AdaptInfo&) override { /* by default, do nothing */ }
 
     /// Implementation of \ref ProblemStatBase::getName().
-    std::string getName() const
+    virtual std::string getName() const override
     {
       return name;
     }

dune/amdis/ProblemStat.hpp

@@ -150,7 +150,7 @@ namespace AMDiS
     
     /// Return the i'th \ref solution component
     template <size_t I = 0>
-    decltype(auto) getSolution(const index_<I> _i = index_<I>())
+    decltype(auto) getSolution(const index_<I> _i = {})
     {
       return (*solution)[_i];
     }
@@ -176,7 +176,7 @@ namespace AMDiS
     
     /// Return the I'th \ref feSpaces component
     template <size_t I = 0>
-    FeSpace<I> const& getFeSpace(const index_<I> = index_<I>()) const
+    FeSpace<I> const& getFeSpace(const index_<I> = {}) const
     {
       return std::get<I>(*feSpaces);
     }
@@ -372,7 +372,7 @@ namespace AMDiS
        **/
       virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
       {
-        for (int i = 0; i < ProblemStatType::getNumComponents(); i++)
+        for (size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
           if (adaptInfo.spaceToleranceReached(i))
             adaptInfo.allowRefinement(false, i);
           else
@@ -391,13 +391,13 @@ namespace AMDiS
       virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
       
       /// Implementation of \ref ProblemStatBase::refineMesh.
-      virtual Flag refineMesh(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
+      virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return {0}; }
 
       /// Implementation of \ref ProblemStatBase::coarsenMesh.
-      virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
+      virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return {0}; }
       
       /// Implementation of \ref ProblemStatBase::markElements.
-      virtual Flag markElements(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
+      virtual Flag markElements(AdaptInfo& adaptInfo) override { return {0}; }
     };
     
   } // end namespace Impl

dune/amdis/ProblemStat.inc.hpp

@@ -222,9 +222,6 @@ namespace AMDiS
 
       For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector, _r](auto const _c) 
       {
-	auto const& rowFeSpace = this->getFeSpace(_r);
-        auto const& colFeSpace = this->getFeSpace(_c);
-	
         // The DOFMatrix which should be assembled
         auto& dofmatrix    = (*systemMatrix)(_r, _c);
         auto& solution_vec = (*solution)[_c];

dune/amdis/linear_algebra/LinearAlgebraBase.hpp

@@ -12,4 +12,41 @@ namespace AMDiS
   template <class Vector>
   using BaseVector = typename Impl::BaseVector<Vector>::type;
     
+  
+  namespace Impl
+  {
+    template <class Vector>
+    class BaseWrapper
+    {
+    public:
+      using value_type = typename Vector::value_type;
+      using size_type = typename Vector::size_type;
+      
+      BaseWrapper(Vector& vec) : vec(vec) {}
+      
+      void resize(size_type s)
+      {
+        vec.resize(s);
+      }
+      
+      size_type size() const
+      {
+        return vec.size();
+      }
+    
+      value_type&       operator[](size_type i)       { return vec[i]; }
+      value_type const& operator[](size_type i) const { return vec[i]; }
+      
+    protected:
+      Vector& vec;
+    };
+    
+  } // end namespace Impl
+  
+  template <class Vector>
+  Impl::BaseWrapper<std::remove_reference_t<Vector>> wrapper(Vector&& vec)
+  {
+    return {std::forward<Vector>(vec)};
+  }
+  
 } // end namespace AMDiS

dune/amdis/linear_algebra/LinearSolverInterface.hpp

@@ -48,7 +48,7 @@ namespace AMDiS
     }
     
     // return the runner/worker corresponding to this solver (optional)
-    virtual std::shared_ptr<RunnerBase> getRunner() { /* Does nothing */ };
+    virtual std::shared_ptr<RunnerBase> getRunner() { return {}; };
 
   private:
     /// main methods that all solvers must implement

dune/amdis/linear_algebra/RunnerInterface.hpp

@@ -33,8 +33,8 @@ namespace AMDiS
       return 0;
     }
                       
-    virtual std::shared_ptr<PreconBase> getLeftPrecon() {}
-    virtual std::shared_ptr<PreconBase> getRightPrecon() {}
+    virtual std::shared_ptr<PreconBase> getLeftPrecon() { return {}; }
+    virtual std::shared_ptr<PreconBase> getRightPrecon() { return {}; }
   };
 
 } // end namespace AMDiS

dune/amdis/linear_algebra/SolverInfo.hpp

@@ -2,6 +2,8 @@
 
 #include <string>
 
+#include <dune/amdis/Initfile.hpp>
+
 namespace AMDiS
 {
   /// Class that stores information about the solution process, like tolerances 

dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp

 #pragma once
 
+#include <boost/numeric/mtl/utility/irange.hpp>
+
 #include <dune/amdis/linear_algebra/LinearAlgebraBase.hpp>
+#include <dune/amdis/linear_algebra/mtl/MTLDenseVector.hpp>
 
 namespace AMDiS
 {
@@ -134,7 +137,7 @@ namespace AMDiS
   };
 
   
-  template <class BlockVector, class Vector = mtl::dense_vector<typename BlockVector::value_type>>
+  template <class BlockVector, class Vector = /**/MTLDenseVector<typename BlockVector::value_type>>
   BlockVectorWrapper<BlockVector, Vector>
   blockWrapper(BlockVector& bvec)
   {

dune/amdis/linear_algebra/mtl/DOFVector.hpp

@@ -3,13 +3,11 @@
 #include <string>
 #include <memory>
 
-
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
-
 #include <dune/functions/functionspacebases/interpolate.hh>
 
 #include <dune/amdis/ClonablePtr.hpp>
 #include <dune/amdis/Log.hpp>
+#include <dune/amdis/linear_algebra/mtl/MTLDenseVector.hpp>
 
 namespace AMDiS
 {
@@ -24,7 +22,7 @@ namespace AMDiS
     using FeSpace    = FeSpaceType;
     
     /// The type of the base vector
-    using BaseVector = mtl::dense_vector<ValueType>;
+    using BaseVector = MTLDenseVector<ValueType>;
     
     /// The index/size - type
     using size_type  = typename FeSpace::size_type;

dune/amdis/linear_algebra/mtl/LinearSolver.hpp

@@ -34,7 +34,7 @@ namespace AMDiS
     {}
     
     /// Implements \ref LinearSolverInterface::getRunner()
-    std::shared_ptr<RunnerBase> getRunner()
+    virtual std::shared_ptr<RunnerBase> getRunner() override
     {
       return runner;
     }

dune/amdis/linear_algebra/mtl/MTLDenseVector.hpp

+#pragma once
+
+
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+
+#include <dune/amdis/linear_algebra/LinearAlgebraBase.hpp>
+
+namespace AMDiS
+{
+  template <class Value>
+  using MTLDenseVector = mtl::vec::dense_vector<Value>;
+  
+  namespace Impl
+  {
+    template <class Vector>
+    class MTLWrapper
+      : public BaseWrapper<Vector>
+    {
+      using Super = BaseWrapper<Vector>;
+      
+    public:
+      using value_type = typename Super::value_type;
+      using size_type = typename Super::size_type;
+      
+      template <class Vector_>
+      MTLWrapper(Vector_&& vec) 
+        : Super(std::forward<Vector_>(vec)) 
+      {}
+      
+      void resize(size_type s)
+      {
+        this->vec.change_dim(s);
+      }
+      
+      size_type size() const
+      {
+        return num_rows(this->vec);
+      }
+    };
+    
+  } // end namespace Impl
+  
+  template <class Value>
+  Impl::MTLWrapper<MTLDenseVector<Value>> wrapper(MTLDenseVector<Value>& vec)
+  {
+    return {vec};
+  }
+  
+}

dune/amdis/linear_algebra/mtl/SystemVector.hpp

@@ -5,8 +5,6 @@
 #include <memory>
 #include <tuple>
 
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
-
 #include <dune/amdis/Basic.hpp>
 #include <dune/amdis/Log.hpp>
 #include <dune/amdis/Loops.hpp>
@@ -80,7 +78,8 @@ namespace AMDiS
     AMDIS_STATIC_ASSERT( nComponents > 0 );
     
     /// The type of the block-vector
-    using MultiVector = BlockMTLVector<mtl::dense_vector<ValueType>, nComponents>;
+    using BaseVector = typename DOFVectorGenerator<std::tuple_element_t<0, FeSpaces>>::BaseVector;
+    using MultiVector = BlockMTLVector<BaseVector, nComponents>;
     
     /// The type of the vector of DOFVectors
     using DOFVectors  = MakeTuple_t<DOFVectorGenerator, FeSpaces>;

install/MTL/include/boost/numeric/itl/iteration/basic_iteration.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_BASIC_ITERATION_INCLUDE
+#define ITL_BASIC_ITERATION_INCLUDE
+
+#include <iostream>
+#include <complex>
+#include <string>
+
+namespace itl {
+
+
+template <class Real>
+class basic_iteration
+{
+  public:
+    typedef basic_iteration self;
+    typedef Real            real;
+
+    template <class Vector>
+    basic_iteration(const Vector& r0, int max_iter_, Real t, Real a = Real(0))
+      : error(0), i(0), norm_r0(std::abs(two_norm(r0))),
+	max_iter(max_iter_), rtol_(t), atol_(a), is_finished(false), my_quite(false), my_suppress(false) { }
+
+    basic_iteration(Real nb, int max_iter_, Real t, Real a = Real(0))
+      : error(0), i(0), norm_r0(nb), max_iter(max_iter_), rtol_(t), atol_(a), is_finished(false), 
+	my_quite(false), my_suppress(false) {}
+
+    virtual ~basic_iteration() {}
+
+    bool check_max()
+    {
+	if (i >= max_iter) 
+	    error= 1, is_finished= true, err_msg= "Too many iterations.";
+	return is_finished;
+    }
+
+    template <class Vector>
+    bool finished(const Vector& r) 
+    {
+	if (converged(two_norm(r)))
+	    return is_finished= true;
+	return check_max();
+    }
+
+    bool finished(const Real& r) 
+    {
+	if (converged(r))
+	    return is_finished= true;
+	return check_max();
+    }
+
+    template <typename T>
+    bool finished(const std::complex<T>& r) 
+    {
+	if (converged(std::abs(r))) 
+	    return is_finished= true;
+	return check_max();
+    }
+
+    bool finished() const { return is_finished; }
+
+    template <class T>
+    int terminate(const T& r) { finished(r); return error; }
+
+    bool converged(const Real& r) { resid_= r; return converged(); }
+
+    bool converged() const 
+    { 
+	if (norm_r0 == 0)
+	    return resid_ <= atol_;  // ignore relative tolerance if |r0| is zero
+	return resid_ <= rtol_ * norm_r0 || resid_ <= atol_;
+    }
+
+    self& operator++() { ++i; return *this; }
+
+    self& operator+=(int n) { i+= n; return *this; }
+
+    bool first() const { return i <= 1; }
+
+    virtual operator int() const { return error; }
+
+    virtual int error_code() const { return error; }
+
+    bool is_converged() const { return is_finished && error == 0; }
+
+    int iterations() const { return i; }
+    
+    int max_iterations() const { return max_iter; }
+
+    void set_max_iterations(int m) { max_iter= m; }
+
+    Real resid() const { return resid_; }
+
+    Real relresid() const { return resid_ / norm_r0; }
+
+    Real normb() const { return norm_r0; }
+
+    Real tol() const { return rtol_; }
+    Real atol() const { return atol_; }
+
+    int fail(int err_code) { error = err_code; return error_code(); }
+
+    int fail(int err_code, const std::string& msg)
+    { error = err_code; err_msg = msg; return error_code(); }
+
+    void set(Real v) { norm_r0 = v; }
+
+    void set_quite(bool q) { my_quite= q; }
+
+    bool is_quite() const { return my_quite; }
+
+    void suppress_resume(bool s) { my_suppress= s; }
+
+    bool resume_suppressed() const { return my_suppress; }
+
+    void update_progress(const basic_iteration& that)
+    {
+	i= that.i;
+	resid_= that.resid_;
+	if (that.error > 1) { // copy error except too many iterations
+	    error= that.error;
+	    err_msg= that.err_msg;
+	    is_finished= true;
+	} else 
+	    finished(resid_);
+    }
+
+  protected:
+    int          error, i;
+    Real         norm_r0;
+    int          max_iter;
+    Real         rtol_, atol_, resid_;
+    std::string  err_msg;
+    bool         is_finished, my_quite, my_suppress;
+};
+
+
+} // namespace itl
+
+#endif // ITL_BASIC_ITERATION_INCLUDE

install/MTL/include/boost/numeric/itl/iteration/cyclic_iteration.hpp

+// Software License for MTL
+// 
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+// 
+// This file is part of the Matrix Template Library2
+// 
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_CYCLIC_ITERATION_INCLUDE
+#define ITL_CYCLIC_ITERATION_INCLUDE
+
+#include <iostream>
+#include <boost/numeric/itl/iteration/basic_iteration.hpp>
+
+namespace itl {
+
+  /// Class for iteration control that cyclically prints residual
+  template <class Real, class OStream = std::ostream>
+  class cyclic_iteration : public basic_iteration<Real> 
+  {
+      typedef basic_iteration<Real> super;
+      typedef cyclic_iteration self;
+
+      void print_resid()
+      {
+	  if (!this->my_quite && this->i % cycle == 0)
+	      if (multi_print || this->i != last_print) { // Avoid multiple print-outs in same iteration
+		  out << "iteration " << this->i << ": resid " << this->resid() 
+		      // << " / " << this->norm_r0 << " = " << this->resid() / this->norm_r0 << " (rel. error)"
+		      << std::endl;
+		  last_print= this->i;
+	      }
+      }
+
+    public:
+  
+      template <class Vector>
+      cyclic_iteration(const Vector& r0, int max_iter_, Real tol_, Real atol_ = Real(0), int cycle_ = 100,
+		       OStream& out = std::cout)
+	: super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out)
+      {}
+
+      cyclic_iteration(Real r0, int max_iter_, Real tol_, Real atol_ = Real(0), int cycle_ = 100,
+		       OStream& out = std::cout)
+	: super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out)
+      {}
+      
+
+      bool finished() { return super::finished(); }
+
+      template <typename T>
+      bool finished(const T& r) 
+      {
+	  bool ret= super::finished(r);
+	  print_resid();
+	  return ret;
+      }
+
+      inline self& operator++() { ++this->i; return *this; }
+      
+      inline self& operator+=(int n) { this->i+= n; return *this; }
+
+      operator int() const { return error_code(); }
+
+      /// Whether the residual is printed multiple times in iteration
+      bool is_multi_print() const { return multi_print; }
+
+      /// Set whether the residual is printed multiple times in iteration
+      void set_multi_print(bool m) { multi_print= m; }
+
+      int error_code() const 
+      {
+	  if (!this->my_suppress)
+	      out << "finished! error code = " << this->error << '\n'
+		  << this->iterations() << " iterations\n"
+		  << this->resid() << " is actual final residual. \n"