DOFVector and SystemVector implemented as MultiType containers

.gitignore

+# ignore generated .vtu and .pvd files,
+*.vtu
+*.pvd
+*.pvtu
+
+# ignore ARH files
+*.arh
+*.parh
+*.tarh
+
+# ignore build directories
+build*/
\ No newline at end of file

CMakeLists.txt

-cmake_minimum_required(VERSION 2.8.12)
-project(amdis CXX)
+cmake_minimum_required(VERSION 3.0)
+project(dune-amdis CXX)
+
+set(UG_DIR /opt/software/dune/lib/cmake/ug)
+set(Vc_DIR /opt/software/dune/lib/cmake/Vc)
+set(ALBERTA_ROOT /opt/software/alberta)
 
 if(NOT (dune-common_DIR OR dune-common_ROOT OR
       "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))

config.h.cmake

-/* begin amdis
+/* begin dune-amdis
    put the definitions for config.h specific to
    your project here. Everything above will be
    overwritten

dune.module

@@ -3,7 +3,7 @@
 ################################
 
 #Name of the module
-Module: amdis
+Module: dune-amdis
 Version: 0.1
 Maintainer: simon.praetorius@tu-dresden.de
 #depending on

dune/amdis/AMDiS.cpp

+#include "AMDiS.hpp"
+
+// std c++ headers
+#include <string>
+
+#include <boost/program_options.hpp>
+
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+
+// AMDiS includes
+#include "Initfile.hpp"
+#include "Log.hpp"
+
+namespace AMDiS
+{
+//   using namespace std;
+
+  void init(int argc, char** argv, std::string initFileName)
+  {
+    // Maybe initialize MPI
+    Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
+    
+    Parameters::clearData();
+
+    // read commandline arguments
+    namespace po = boost::program_options;
+
+    // Declare the supported options.
+    po::options_description desc("Usage: " + std::string(argv[0]) + " init-file [options]\nAllowed options");
+    desc.add_options()
+    ("help", "produce help message")
+    ("init-file", po::value<std::string>(), "set init file")
+    ("parameters", po::value<std::string>(), "set parameter in init file\nsyntax: \"key1: value1; key2: value2...\"");
+
+    po::options_description hidden("Hidden options");
+    hidden.add_options()
+    ("unknown", po::value<std::vector<std::string>>(), "unknown options");
+
+    po::options_description cmdline_options;
+    cmdline_options.add(desc).add(hidden);
+
+    // first argument is init-filename
+    po::positional_options_description p;
+    p.add("init-file", 1);
+    p.add("unknown", -1);
+
+    // parse comandline
+    po::variables_map vm;
+    po::store(po::command_line_parser(argc, argv).options(cmdline_options).positional(p).allow_unregistered().run(), vm);
+    po::notify(vm);
+
+    // print help message
+    if (vm.count("help"))
+    {
+      std::cout << desc << "\n";
+      exit(1);
+    }
+
+    // set parameters before reading the initfile
+    // if (vm.count("parameters"))
+    //   Parameters::readArgv(vm["parameters"].as<std::string>());
+
+    if (initFileName == "")
+    {
+      if (vm.count("init-file"))
+        Parameters::init(vm["init-file"].as<std::string>());
+      else
+      {
+        AMDIS_ERROR_EXIT("No init file specified!\n");
+      }
+    }
+    else
+    {
+      Parameters::init(initFileName);
+    }
+
+    // reset parameters from command line
+    bool ignoreCommandline = false;
+    Parameters::get("ignore commandline options", ignoreCommandline);
+    // if (vm.count("parameters") && !ignoreCommandline)
+    //   Parameters::readArgv(vm["parameters"].as<std::string>(),0);
+    // TODO: add command-line arguments again.
+  }
+
+
+  void init(std::string initFileName)
+  {
+    Parameters::init(initFileName);
+  }
+
+
+  void finalize()
+  {}
+
+} // end namespace AMDiS

dune/amdis/AMDiS.hpp

+#pragma once
+
+// std c++ headers
+#include <string>
+
+#include <dune/common/exceptions.hh> // We use exceptions
+
+namespace AMDiS
+{
+  void init(int argc, char** argv, std::string initFileName = "");
+
+  void init(std::string initFileName);
+
+  void finalize();
+
+} // end namespace AMDiS

src/AdaptInfo.hpp → dune/amdis/AdaptInfo.hpp

@@ -6,8 +6,8 @@
 #include <vector>
 
 // AMDiS includes
-#include "Log.hpp"
-#include "Math.hpp"
+#include <dune/amdis/Log.hpp>
+#include <dune/amdis/Math.hpp>
 
 namespace AMDiS
 {
@@ -303,9 +303,9 @@ namespace AMDiS
     /// Returns \ref est_sum.
     double getEstSum(int index) const
     {
-      FUNCNAME_DBG("AdaptInfo::getEstSum()");
-      TEST_EXIT_DBG(static_cast<size_t>(index) < scalContents.size())
-      ("Wrong index for adaptInfo!\n");
+      AMDIS_FUNCNAME_DBG("AdaptInfo::getEstSum()");
+      AMDIS_TEST_EXIT_DBG(static_cast<size_t>(index) < scalContents.size(),
+	"Wrong index for adaptInfo!\n");
 
       return scalContents[index]->est_sum;
     }
@@ -319,9 +319,9 @@ namespace AMDiS
     /// Returns \ref est_max.
     double getEstMax(int index) const
     {
-      FUNCNAME_DBG("AdaptInfo::getEstSum()");
-      TEST_EXIT_DBG(static_cast<size_t>(index) < scalContents.size())
-      ("Wrong index for adaptInfo!\n");
+      AMDIS_FUNCNAME_DBG("AdaptInfo::getEstSum()");
+      AMDIS_TEST_EXIT_DBG(static_cast<size_t>(index) < scalContents.size(),
+	"Wrong index for adaptInfo!\n");
 
       return scalContents[index]->est_max;
     }

src/Basic.hpp → dune/amdis/Basic.hpp

@@ -2,6 +2,10 @@
 
 #include <cassert>
 #include <tuple>
+#include <map>
+#include <list>
+
+#include "IndexSeq.hpp"
 
 namespace AMDiS
 {  
@@ -11,6 +15,8 @@ namespace AMDiS
     template <bool B>
     using bool_ = std::integral_constant<bool, B>;
 
+    template <size_t I>
+    using index_ = std::integral_constant<size_t, I>;
     
     template <class T>
     using IdxPairList = std::map< std::pair<int, int>, std::list<std::shared_ptr<T> > >;
@@ -32,15 +38,15 @@ namespace AMDiS
 	template <class Tuple, class... Args>
 	Tuple construct_tuple_aux(int_<1>, Args&&... args)
 	{
-	    static_assert(std::tuples_size<Tuple>::value == sizeof...(args), 
+	    static_assert(std::tuple_size<Tuple>::value == sizeof...(args), 
 		"Nr. of argument != tuple-size");
-	    return {std::forward<Args>(args)...};
+	    return Tuple{std::forward<Args>(args)...};
 	}
 	
 	template <class Tuple, class... Args>
 	Tuple construct_tuple_aux(int_<0>, Args&&... args)
 	{
-	    static_assert(std::tuples_size<Tuple>::value == 0, 
+	    static_assert(std::tuple_size<Tuple>::value == 0, 
 		"Construction of empty tuples with empty argument list only!");
 	    return {};
 	}
@@ -50,8 +56,67 @@ namespace AMDiS
     template <class Tuple, class Arg>
     Tuple construct_tuple(Arg&& arg)
     {
-	return Impl::construct_tuple_aux<Tuple>(int_<std::tuples_size<Tuple>::value>(), 
+	return Impl::construct_tuple_aux<Tuple>(int_<std::tuple_size<Tuple>::value>(), 
 	std::forward<Arg>(arg));
     }
+    
+    // -----------
+    
+    template <template<class> class Base, class Tuple, class Indices> struct MakeTuple;
+    
+    template <template<class> class Base, class Tuple, size_t... I> 
+    struct MakeTuple<Base, Tuple, Seq<I...>>
+    {
+	using type = std::tuple<Base<std::tuple_element_t<I, Tuple>>...>;
+    };
+    
+    template <template<class> class Base, class Tuple> 
+    using MakeTuple_t = 
+	typename MakeTuple<Base, Tuple, MakeSeq_t<std::tuple_size<Tuple>::value>>::type;
+    
+    // -----------
+    
+    template <template<class,class> class Base, class Tuple1, class Tuple2, class Indices> struct MakeTuple2;
+    
+    template <template<class,class> class Base, class Tuple1, class Tuple2, size_t... I> 
+    struct MakeTuple2<Base, Tuple1, Tuple2, Seq<I...>>
+    {
+	using type = std::tuple<Base<std::tuple_element_t<I, Tuple1>, std::tuple_element_t<I, Tuple2>>...>;
+    };
+    
+    template <template<class,class> class Base, class Tuple1, class Tuple2> 
+    using MakeTuple2_t = 
+	typename MakeTuple2<Base, Tuple1, Tuple2, MakeSeq_t<std::tuple_size<Tuple1>::value>>::type;
+    
+    // -----------
+    
+    template <template<class...> class Base, class Tuple, class Indices> struct ExpandTuple;
+    
+    template <template<class...> class Base, class Tuple, size_t... I> 
+    struct ExpandTuple<Base, Tuple, Seq<I...>>
+    {
+	using type = Base<std::tuple_element_t<I, Tuple>...>;
+    };
+    
+    template <template<class...> class Base, class Tuple> 
+    using ExpandTuple_t = 
+	typename ExpandTuple<Base, Tuple, MakeSeq_t<std::tuple_size<Tuple>::value>>::type;
+    
+    
+    // -----------
+    
+    template <class T, class Indices> struct RepeatedTuple;
+    
+    template <class T, size_t... I> 
+    struct RepeatedTuple<T, Seq<I...>>
+    {
+	template <size_t, class U> using Id = U;	
+	using type = std::tuple<Id<I, T>...>;
+    };
+    
+    template <size_t N, class T> 
+    using Repeat_t = 
+	typename RepeatedTuple<T, MakeSeq_t<N>>::type;
+    
   
 } // end namespace AMDiS

dune/amdis/CMakeLists.txt

 #install headers
-install(FILES amdis.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis)
+
+dune_add_library("duneamdis" NO_EXPORT
+    AdaptInfo.cpp
+    AMDiS.cpp
+    Initfile.cpp
+    ProblemStat.cpp
+    SystemVector.cpp
+    )
+add_dune_alberta_flags("duneamdis" OBJECT USE_GENERIC)
+
+set(BOOST_VERSION "1.54")
+set(BOOST_LIBS_REQUIRED system program_options)
+if (NOT BUILD_SHARED_LIBS)
+    set(Boost_USE_STATIC_LIBS ON)
+endif (NOT BUILD_SHARED_LIBS)
+set(Boost_NO_SYSTEM_PATHS ON)
+find_package(Boost ${BOOST_VERSION} REQUIRED ${BOOST_LIBS_REQUIRED})
+if (Boost_FOUND)
+    add_library(boost INTERFACE)
+    target_include_directories(boost INTERFACE ${Boost_INCLUDE_DIR})
+    target_link_libraries(boost INTERFACE  ${Boost_LIBRARIES})
+    target_link_libraries("duneamdis" INTERFACE boost)
+
+    if (MSVC_SHARED_LIBS)
+        link_directories(${Boost_LIBRARY_DIRS})
+        target_compile_definitions("duneamdis" INTERFACE ${Boost_LIB_DIAGNOSTIC_DEFINITIONS})
+    endif (MSVC_SHARED_LIBS)
+endif (Boost_FOUND)
+
+
+install(FILES 
+    AdaptInfo.hpp
+    AMDiS.hpp
+    Basic.hpp
+    DirichletBC.hpp
+    DirichletBC.inc.hpp
+    DOFVector.hpp
+    Flag.hpp
+    IndexSeq.hpp
+    Initfile.hpp
+    Log.hpp
+    Loop.hpp
+    Math.hpp
+    Operator.hpp
+    Operator.inc.hpp
+    OperatorTerm.hpp
+    ProblemStat.hpp
+    ProblemStat.inc.hpp
+    ProblemStatBase.hpp
+    SystemVector.hpp
+    Timer.hpp
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis)

dune/amdis/DOFVector.hpp

+#pragma once
+
+#include <string>
+#include <memory>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/istl/bvector.hh>
+
+#include "Log.hpp"
+
+namespace AMDiS
+{
+    template <class FeSpaceType, class ValueType = Dune::FieldVector<double,1>>
+    class DOFVector
+    {
+    public:
+	using FeSpace    = FeSpaceType;
+	using BaseVector = Dune::BlockVector<ValueType>;
+	
+	using size_type  = typename FeSpace::size_type;
+	using field_type = typename ValueType::field_type;
+	
+	using value_type = ValueType;
+	
+	/// Constructor.
+	DOFVector(FeSpace const& feSpace, std::string name)
+	  : feSpace(feSpace)
+	  , name(name)
+	  , vector(new BaseVector())
+	  , allocated(true)
+	{
+	    compress();
+	}
+	
+	/// Constructor. Takes pointer of data-vector.
+	DOFVector(FeSpace const& feSpace, std::string name, 
+		  BaseVector& vector)
+	  : feSpace(feSpace)
+	  , name(name)
+	  , vector(&vector)
+	{}
+	
+	~DOFVector()
+	{
+	    if (allocated)
+		delete vector;
+	}
+	
+	/// Return the basis \ref feSpace of the vector
+	FeSpace const& getFeSpace() const
+	{
+	    return feSpace;
+	}
+	
+	/// Return the data-vector \ref vector
+	BaseVector const& getVector() const
+	{
+	    return *vector;
+	}
+	
+	/// Return the data-vector \ref vector
+	BaseVector& getVector()
+	{
+	    return *vector;
+	}
+	
+	/// Return the size of the \ref feSpace
+	size_type getSize() const
+	{
+	    return feSpace.size();
+	}
+	
+	/// Return the \ref name of this vector
+	std::string getName() const
+	{
+	    return name;
+	}
+	
+	/// Resize the \ref vector to the size of the \ref feSpace.
+	void compress()
+	{
+	    vector->resize(getSize());
+	}
+	
+	
+	/// Access the entry \p i of the \ref vector with read-access.
+	value_type const& operator[](size_type i) const
+	{
+	    AMDIS_TEST_EXIT_DBG( i < vector->size() , 
+		"index " << i << " out of range [" << 0 << "," << vector->size() << ")" );
+	    return (*vector)[i];
+	}
+	
+	/// Access the entry \p i of the \ref vector with write-access.
+	value_type& operator[](size_type i)
+	{
+	    AMDIS_TEST_EXIT_DBG( i < vector->size() , 
+		"index " << i << " out of range [" << 0 << "," << vector->size() << ")" );
+	    return (*vector)[i];
+	}
+	
+	/// interpolate a function \p f to the basis \ref feSpace and store the
+	/// coefficients in \ref vector.
+	template <class F>
+	void interpol(F const& f)
+	{
+	    Dune::Functions::interpolate(feSpace, *vector, f);
+	}
+	
+    private:
+	FeSpace const&	feSpace;
+	std::string	name;
+	
+	BaseVector* vector;
+	bool allocated = false;
+	
+    // friend class declarations
+    template <class, class>
+    friend class SystemVector;
+    };
+
+} // end namespace AMDiS

src/DirichletBC.hpp → dune/amdis/DirichletBC.hpp

@@ -13,11 +13,11 @@ namespace AMDiS
   {
   public:
     template <class Predicate, class Values,
-      class = std::enable_if_t< Functions::Concept::isFunction<Predicate, bool(WorldVector)>() &&
-				Functions::Concept::isFunction<Values,  double(WorldVector)>() > >
+      class = std::enable_if_t< Dune::Functions::Concept::isFunction<Predicate, bool(WorldVector)>() &&
+				Dune::Functions::Concept::isFunction<Values,  double(WorldVector)>() > >
     DirichletBC(Predicate&& predicate, Values&& values)
       : predicate(std::forward<Predicate>(predicate))
-      , values(std::forward<Values>(values)
+      , values(std::forward<Values>(values))
     {}
     
     

src/DirichletBC.inc.hpp → dune/amdis/DirichletBC.inc.hpp

@@ -2,9 +2,12 @@
 
 #include <dune/functions/functionspacebases/interpolate.hh>
 
+#include "Log.hpp"
+
 namespace AMDiS
 {
-  template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
+  template <class WorldVector>
+    template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
   void DirichletBC<WorldVector>::init(bool apply, 
 				      RowFeSpace const& rowFeSpace, 
 				      ColFeSpace const& colFeSpace, 
@@ -14,12 +17,15 @@ namespace AMDiS
   {
     using Dune::Functions::interpolate;
     
-    if (!initialized)
+    if (!initialized) {
       interpolate(rowFeSpace, dirichletNodes, predicate);
+      initialized = true;
+    }
   }
     
     
-  template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
+  template <class WorldVector>
+    template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
   void DirichletBC<WorldVector>::finish(bool apply, 
 					RowFeSpace const& rowFeSpace, 
 					ColFeSpace const& colFeSpace, 
@@ -29,10 +35,10 @@ namespace AMDiS
   {
     using Dune::Functions::interpolate;
     
-    assert( initialized );
+    AMDIS_TEST_EXIT( initialized, "Boundary condition not initialized!" );
     
     // loop over the matrix rows
-    for (size_t i = 0; i < stiffnessMatrix.N(); ++i) {
+    for (size_t i = 0; i < matrix->N(); ++i) {
 	if (dirichletNodes[i]) {
 	    auto cIt = (*matrix)[i].begin();
 	    auto cEndIt = (*matrix)[i].end();

dune/amdis/IndexSeq.hpp

+#pragma once
+
+// from http://stackoverflow.com/questions/17424477/implementation-c14-make-integer-sequence/17426611#17426611
+namespace AMDiS
+{
+  template <size_t...> struct Seq {};
+
+  namespace detail
+  {
+    template <size_t s, class S>
+    struct Concat;
+
+    template <size_t s, size_t ... i>
+    struct Concat<s, Seq<i... >>
+    {
+      using type = Seq<i..., (s + i)... >;
+    };
+
+    template <bool, class S>
+    struct IncSeq_if 
+    { 
+      using type = S; 
+    };
+
+    template <size_t... Is>
+    struct IncSeq_if<true, Seq<Is...>>
+    {
+      using type = Seq<Is..., sizeof...(Is)>;
+    };
+    
+  } // end namespace detail
+  
+  template <size_t N>
+  struct MakeSeq 
+  {
+    using type = typename detail::IncSeq_if< (N % 2 != 0), 
+      typename detail::Concat<N/2, typename MakeSeq<N/2>::type>::type >::type;
+  };
+
+  // break condition
+  template <> struct MakeSeq<0> { using type = Seq<>; };
+
+  // alias template
+  template <size_t N>
+  using MakeSeq_t = typename MakeSeq<N>::type;
+  
+} // end namespace AMDiS

src/Initfile.cpp → dune/amdis/Initfile.cpp

@@ -12,7 +12,7 @@
 #include <boost/property_tree/json_parser.hpp>
 
 // a parser for arithmetic expressions
-#include <muParser.h>
+// #include <muParser.h>
 
 #include "Log.hpp"
 
@@ -29,18 +29,18 @@ namespace AMDiS
   }
 
 
-  namespace detail
-  {
-    double mu_parser_eval(std::string const& valStr)
-    {
-      mu::Parser parser;
-      parser.DefineConst(_T("M_PI"), m_pi);
-      parser.DefineConst(_T("M_E"), m_e);
-
-      parser.SetExpr(valStr);
-      return parser.Eval();
-    }
-  }
+//   namespace detail
+//   {
+//     double mu_parser_eval(std::string const& valStr)
+//     {
+//       mu::Parser parser;