Merge branch 'develop' into derivative_term

CMakeLists.txt

 cmake_minimum_required(VERSION 3.0)
 project(dune-amdis CXX)
 
-set(ALBERTA_ROOT /opt/software/alberta)
-set(ALBERTA_INCLUDE_DIR  /opt/software/alberta/include)
+# set(ALBERTA_ROOT /opt/software/alberta)
+# set(ALBERTA_INCLUDE_DIR  /opt/software/alberta/include)
 # set(CMAKE_PREFIX_PATH /opt/software/dune/lib/cmake)
 # set(UG_DIR /opt/software/dune/lib/cmake/ug)
 # set(Vc_DIR /opt/software/dune/lib/cmake/Vc)
 
-set(MTL_DIR ${CMAKE_SOURCE_DIR}/install/MTL/share/mtl)
-
 if(NOT (dune-common_DIR OR dune-common_ROOT OR
       "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
     string(REPLACE  ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR

doc/doxygen/Doxylocal

 # This file contains local changes to the doxygen configuration
 # please us '+=' to add file/directories to the lists
 
+FILE_PATTERNS          += *.hpp \
+                          *.cpp
+
+HIDE_SCOPE_NAMES       = YES
+HIDE_UNDOC_CLASSES     = NO
+INTERNAL_DOCS          = NO
+MARKDOWN_SUPPORT       = YES
+
+EXCLUDE_SYMBOLS        = AMDiS::Impl \
+                         AMDiS::traits::Impl \
+                         AMDiS::detail \
+                         itl::details
+                           
+PREDEFINED            += HAVE_UMFPACK \
+                         HAVE_ALBERTA \
+                         HAVE_UG \
+                         AMDIS_BACKEND_MTL
+                          
 # The INPUT tag can be used to specify the files and/or directories that contain
 # documented source files. You may enter file names like "myfile.cpp" or
 # directories like "/usr/src/myproject". Separate the files or directories
 # with spaces.
 
-INPUT                 += @top_srcdir@/dune/
+INPUT                 += @top_srcdir@/dune/amdis \
+                         @top_srcdir@/dune/amdis/common \
+                         @top_srcdir@/dune/amdis/utility \
+                         @top_srcdir@/dune/amdis/linear_algebra \
+                         @top_srcdir@/dune/amdis/linear_algebra/mtl
 # see e.g. dune-grid for the examples of mainpage and modules
-# INPUT                 += @srcdir@/mainpage \
-#                          @srcdir@/modules
+#INPUT                 += @srcdir@/mainpage \
+#                         @srcdir@/modules
 
 # The EXCLUDE tag can be used to specify files and/or directories that should
 # excluded from the INPUT source files. This way you can easily exclude a
 # subdirectory from a directory tree whose root is specified with the INPUT tag.
 
-# EXCLUDE               += @top_srcdir@/dune/amdis/test
+EXCLUDE               += @top_srcdir@/dune/amdis/test \
+                         @top_srcdir@/dune/amdis/linear_algebra/istl
 
 # The EXAMPLE_PATH tag can be used to specify one or more files or
 # directories that contain example code fragments that are included (see
 # the \include command).
 
-# EXAMPLE_PATH          += @top_srcdir@/src
+EXAMPLE_PATH          += @top_srcdir@/src
 
 # The IMAGE_PATH tag can be used to specify one or more files or
 # directories that contain image that are included in the documentation (see

dune-amdis.pc.in

@@ -9,7 +9,7 @@ DEPENDENCIES=@REQUIRES@
 Name: @PACKAGE_NAME@
 Version: @VERSION@
 Description: amdis module
-URL: http://dune-project.org/
+URL: http://www.github.com/spraetor
 Requires: dune-common dune-geometry dune-localfunctions dune-istl dune-typetree dune-grid dune-functions
 Libs: -L${libdir}
 Cflags: -I${includedir}

dune/amdis/BoundaryElementIterator.hpp

+#pragma once
+
+namespace AMDiS
+{
+  template <class MeshView>
+  class BoundaryFacetIterator;
+
+  // An Iterator over all elements and when element hasBoundaryIntersections
+  template <class MeshView>
+  class BoundaryElementIterator
+  {
+    friend class BoundaryFacetIterator<MeshView>;
+    
+    using Self = BoundaryElementIterator;
+    using Element = typename MeshView::template Codim<0>::Entity;
+    using ElementIterator = typename MeshView::template Codim<0>::Iterator;
+    
+    class Iterator
+    {
+    public:
+      Iterator(ElementIterator elementIt, ElementIterator endIt)
+        : elementIt(elementIt)
+        , endIt(endIt)
+      {}
+      
+      Iterator(Iterator const&) = default;
+      Iterator& operator=(Iterator const&) = default;
+      
+      // iterate to next element that has boundary intersections
+      Iterator& operator++()
+      {
+        ++elementIt;
+        while (elementIt != endIt && !elementIt->hasBoundaryIntersections())
+          ++elementIt;
+        return *this;
+      }
+      
+      Iterator operator++(int)
+      {
+        auto tmp = *this;
+        ++(*this);
+        return tmp;
+      }
+      
+      Element const& operator*() const
+      {
+        return *elementIt;
+      }
+      
+      Element const* operator->() const
+      {
+        return &(*elementIt);
+      }
+      
+      bool operator==(Iterator const& that) const
+      {
+        return elementIt == that.elementIt;
+      }
+      
+      bool operator!=(Iterator const& that) const
+      {
+        return !(*this == that);
+      }
+      
+    private:
+      ElementIterator elementIt;
+      ElementIterator endIt;
+    };
+
+  public:
+    /// Constructor.
+    BoundaryElementIterator(MeshView const& meshView)
+      : meshView(meshView)
+    {}
+    
+    Iterator begin() {
+      auto elementIt = elements(meshView).begin();
+      auto endIt = elements(meshView).end();
+      while (elementIt != endIt && !elementIt->hasBoundaryIntersections())
+        ++elementIt;
+      
+      return {elementIt, endIt};
+    }
+    
+    Iterator end() {
+      return {elements(meshView).end(), elements(meshView).end()};
+    }
+    
+  private:
+    MeshView const& meshView;
+  };
+  
+  
+  /// Generator function for the boundary-element iterator
+  template <class MeshView>
+  BoundaryElementIterator<MeshView> boundary_elements(MeshView const& meshView)
+  {
+    return {meshView};
+  }
+  
+
+  // An Iterator over all elements and when element hasBoundaryIntersections, then 
+  // iterate over all boundary-intersections with given Index, oder for thos with 
+  // predicate returns true
+  template <class MeshView>
+  class BoundaryFacetIterator
+  {
+    using Element = typename MeshView::template Codim<0>::Entity;
+    using Facet   = typename MeshView::template Codim<1>::Entity;
+    
+    using ElementIterator = typename BoundaryElementIterator<MeshView>::Iterator;
+    using FacetIterator = typename MeshView::IntersectionIterator;
+    
+    class Iterator
+    {
+    public:
+      Iterator(MeshView const& meshView, 
+               ElementIterator elementIt, 
+               FacetIterator facetIt, 
+               ElementIterator endIt)
+        : meshView(meshView)
+        , elementIt(elementIt)
+        , facetIt(facetIt)
+        , endIt(endIt)
+      {}
+      
+      Iterator(Iterator const&) = default;
+      Iterator& operator=(Iterator const&) = default;
+      
+      // iterate to next boundary face within current element or to the first
+      // boundary face in the next boundary element
+      Iterator& operator++()
+      {
+        ++facetIt;
+        do {
+          auto facetEndIt = intersections(meshView, *elementIt).end();
+          while (facetIt != facetEndIt && !facetIt->boundary())
+            ++facetIt;
+          if (facetIt != facetEndIt)
+            break;
+          
+          ++elementIt;
+          facetIt = intersections(meshView, *elementIt).begin();
+        } while (elementIt != endIt);
+        
+        return *this;
+      }
+      
+      Iterator operator++(int)
+      {
+        auto tmp = *this;
+        ++(*this);
+        return tmp;
+      }
+      
+      // returns Dune::Intersection
+      auto const& operator*() const
+      {
+        return *facetIt;
+      }
+      
+      auto const* operator->() const
+      {
+        return &(*facetIt);
+      }
+      
+      bool operator==(Iterator const& that) const
+      {
+        return elementIt == that.elementIt && (elementIt == endIt || facetIt == that.facetIt);
+      }
+      
+      bool operator!=(Iterator const& that) const
+      {
+        return !(*this == that);
+      }
+      
+    private:
+      MeshView const& meshView;
+      ElementIterator elementIt;
+      FacetIterator   facetIt;
+      ElementIterator endIt;
+    };
+
+  public:
+    /// Constructor.
+    BoundaryFacetIterator(MeshView const& meshView)
+      : meshView(meshView)
+    {}
+    
+    Iterator begin() {
+      auto elementIt = boundary_elements(meshView).begin();
+      auto endElementIt = boundary_elements(meshView).end();
+      auto facetIt = intersections(meshView, *elementIt).begin();
+      
+      while (!facetIt->boundary())
+        ++facetIt;
+      
+      return {meshView, elementIt, facetIt, endElementIt};
+    }
+    
+    Iterator end() {
+      auto elementIt = boundary_elements(meshView).begin();
+      auto endElementIt = boundary_elements(meshView).end();
+      auto facetIt = intersections(meshView, *elementIt).begin();
+      return {meshView, endElementIt, facetIt, endElementIt};
+    }
+    
+  private:
+    MeshView const& meshView;
+  };
+  
+  
+  /// Generator function for the boundary-element iterator
+  template <class MeshView>
+  BoundaryFacetIterator<MeshView> boundary_facets(MeshView const& meshView)
+  {
+    return {meshView};
+  }
+
+} // end namespace AMDiS

dune/amdis/CMakeLists.txt

@@ -39,7 +39,9 @@ if (Boost_FOUND)
 endif (Boost_FOUND)
 
 
-find_package(MTL REQUIRED)
+find_package(MTL REQUIRED
+             PATHS ${CMAKE_SOURCE_DIR}/install/MTL/share/mtl
+                   /iwr/modules/mtl4/share/mtl)
 if (MTL_FOUND)
   target_include_directories("duneamdis" PUBLIC ${MTL_INCLUDE_DIRS})
 #   target_link_libraries("duneamdis" PUBLIC ${MTL_LIBRARIES})
@@ -49,14 +51,14 @@ if (MTL_FOUND)
   foreach (feature ${CXX_ELEVEN_FEATURE_LIST})
       target_compile_definitions("duneamdis" PUBLIC MTL_WITH_${feature})
   endforeach ()
-  
+
   if (HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_FOUND)
       target_compile_definitions("duneamdis" PUBLIC MTL_HAS_UMFPACK)
   endif ()
 endif (MTL_FOUND)
 
 
-install(FILES 
+install(FILES
     AdaptBase.hpp
     AdaptInfo.hpp
     AdaptInstationary.hpp

dune/amdis/CreatorInterface.hpp

@@ -3,6 +3,7 @@
 #include <memory>
 #include <string>
 
+#include <dune/amdis/common/Utility.hpp>
 #include <dune/amdis/Log.hpp>
 
 namespace AMDiS
@@ -30,7 +31,7 @@ namespace AMDiS
      * Must be implemented by sub classes of CreatorInterface.
      * Creates a new instance of the sub class of BaseClass.
      */
-    virtual std::shared_ptr<BaseClass> create() = 0; 
+    virtual shared_ptr<BaseClass> create() = 0; 
   };
 
   /**
@@ -44,9 +45,9 @@ namespace AMDiS
   {
   public:
     
-    virtual std::shared_ptr<BaseClass> create() override
+    virtual shared_ptr<BaseClass> create() override
     { 
-      AMDIS_ERROR_EXIT("This create method should not be called. Call create(string) instead!"); 
+      AMDIS_ERROR_EXIT("Should not be called. Call create(string) instead!"); 
       return {}; 
     };
     
@@ -55,7 +56,7 @@ namespace AMDiS
      * Creates a new instance of the sub class of BaseClass by passing a 
      * string to the constructor.
      */
-    virtual std::shared_ptr<BaseClass> create(std::string) = 0;
+    virtual shared_ptr<BaseClass> create(std::string) = 0;
   };
   
   

dune/amdis/DirichletBC.hpp

@@ -10,13 +10,28 @@
 
 namespace AMDiS
 {
+  /// Implements a boundary condition of Dirichlet-type.
+  /**
+   * By calling the methods \ref init() and \ref finish before and after 
+   * assembling the system-matrix, respectively, dirichlet boundary conditions
+   * can be applied to the matrix and system vector. Therefore a predicate 
+   * functions indicates the DOFs where values should be enforced and a second
+   * functor provided in the constructor is responsible for determining the 
+   * values to be set at the DOFs.
+   * 
+   * In the \ref finish method the matrix is called with \ref applyDirichletBC 
+   * to erase the corresponding rows and columns for the DOF indices. This 
+   * application of boundary conditions can be symmetric if the matrix does 
+   * support this symmetric modification. As a result, this method returns a list
+   * of columns values, that should be subtracted from the rhs.
+   **/
   template <class WorldVector>
   class DirichletBC
   {
   public:
     template <class Predicate, class Values,
-      class = std::enable_if_t< concepts::Functor<Predicate, bool(WorldVector)> &&
-				concepts::Functor<Values,  double(WorldVector)> > >
+      class = std::enable_if_t< Concepts::Functor<Predicate, bool(WorldVector)> &&
+				Concepts::Functor<Values,  double(WorldVector)> > >
     DirichletBC(Predicate&& predicate, Values&& values)
       : predicate(std::forward<Predicate>(predicate))
       , values(std::forward<Values>(values))

dune/amdis/FileWriter.hpp

@@ -45,8 +45,8 @@ namespace AMDiS
       dir = "output";
       Parameters::get(base + "->output directory", dir);
       
-      vtkWriter = std::make_shared<Dune::VTKWriter<MeshView>>(meshView);
-      vtkSeqWriter = std::make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
+      vtkWriter = make_shared<Dune::VTKWriter<MeshView>>(meshView);
+      vtkSeqWriter = make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
     }
 
     
@@ -125,8 +125,8 @@ namespace AMDiS
 
   private:
     MeshView const& meshView;
-    std::shared_ptr<Dune::VTKWriter<MeshView>> vtkWriter;
-    std::shared_ptr<Dune::VTKSequenceWriter<MeshView>> vtkSeqWriter;
+    shared_ptr<Dune::VTKWriter<MeshView>> vtkWriter;
+    shared_ptr<Dune::VTKSequenceWriter<MeshView>> vtkSeqWriter;
     
     std::vector<std::string> names;
     std::array<std::vector<double>, nComponents> data_vectors;

dune/amdis/Math.hpp

-/** \defgroup Common Common
- */
-
 #pragma once
 
 #include <string>
@@ -8,26 +5,40 @@
 #include <cmath>
 #include <cfloat>
 
+#include <boost/math/special_functions/pow.hpp> 
+
+#include <dune/amdis/common/ScalarTypes.hpp>
+
 namespace AMDiS
 {
-  namespace math
+  namespace Math
   {
-    template <class T>
-    std::enable_if_t<std::is_arithmetic<T>::value, T>
-    constexpr abs(T a)
+    /// Implementation of the absolute value \f$|a|\f$ of arithmetic types.
+    template <class T,
+      class = std::enable_if_t<Concepts::Arithmetic<T>> >
+    constexpr T abs(T a)
     {
       return  a >= 0 ? a : -a;
     }
 
 
-    template <class T>
-    std::enable_if_t<std::is_arithmetic<T>::value, T>
-    constexpr sqr(T a)
+    /// Implementation of the square \f$ a^2 \f$ of arithmetic types.
+    template <class T,
+      class = std::enable_if_t<Concepts::Arithmetic<T>> >
+    constexpr auto sqr(T a)
     {
       return a*a;
     }
 
+    template <size_t p, class T,
+      class = std::enable_if_t<Concepts::Arithmetic<T>> >
+    constexpr auto pow(T v)
+    {
+      return boost::math::pow<p>(v);
+    }
 
+    /// Implementation of the minimum of two values \f$ min(a,b)\f$ of any type 
+    /// supporting the `>` relation.
     template <class T0, class T1>
     constexpr auto min(T0 a, T1 b)
     {
@@ -35,13 +46,15 @@ namespace AMDiS
     }
 
 
+    /// Implementation of the maximum of two values \f$ max(a,b)\f$ of any type 
+    /// supporting the `>` relation.
     template <class T0, class T1>
     constexpr auto max(T0 a, T1 b)
     {
       return a > b ? a : b;
     }
 
-  } // end namespace math
+  } // end namespace Math
 
 
   template <class T> inline void nullify(T& a)

dune/amdis/Mesh.hpp

@@ -7,61 +7,75 @@
 #include <array>
 #include <memory>
 
+#include <dune/common/array.hh>
 #include <dune/common/fvector.hh>
 
 #include <dune/grid/albertagrid.hh>
 #include <dune/grid/uggrid.hh>
 #include <dune/grid/yaspgrid.hh>
 
-#include "Initfile.hpp"
-#include "Log.hpp"
+#include <dune/grid/albertagrid/albertareader.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+
+#include <dune/amdis/Initfile.hpp>
+#include <dune/amdis/Log.hpp>
+#include <dune/amdis/common/Filesystem.hpp>
 
 namespace AMDiS
 {
-  struct _albertagrid {};
-  struct _uggrid {};
-  struct _yaspgrid {};
-  struct _unknowngrid {};
-
-  template <class Grid>
-  struct MeshTagImpl
+  namespace tag 
   {
-    using type = _unknowngrid;
-  };
-  
-  template <class Grid>
-  using MeshTag = typename MeshTagImpl<Grid>::type;
+    struct albertagrid {};
+    struct uggrid {};
+    struct yaspgrid {};
+    struct unknowngrid {};
+    
+  } // end namespace tag
   
-  // specialization for some grid types from DUNE
-#if HAVE_ALBERTA
-  template <int dim, int dimworld>
-  struct MeshTagImpl<Dune::AlbertaGrid<dim, dimworld>>
+  namespace Impl
   {
-    using type = _albertagrid;
-  };
+    template <class Grid>
+    struct MeshTag
+    {
+      using type = tag::unknowngrid;
+    };
+    
+    // specialization for some grid types from DUNE
+#if HAVE_ALBERTA
+    template <int dim, int dimworld>
+    struct MeshTag<Dune::AlbertaGrid<dim, dimworld>>
+    {
+      using type = tag::albertagrid;
+    };
 #endif
-  
+    
 #if HAVE_UG
-  template <int dim>
-  struct MeshTagImpl<Dune::UGGrid<dim>>
-  {
-    using type = _uggrid;
-  };
+    template <int dim>
+    struct MeshTag<Dune::UGGrid<dim>>
+    {
+      using type = tag::uggrid;
+    };
 #endif
+    
+    template <int dim, class Coordinates>
+    struct MeshTag<Dune::YaspGrid<dim, Coordinates>>
+    {
+      using type = tag::yaspgrid;
+    };
+    
+  } // end namespace Impl
   
-  template <int dim, class Coordinates>
-  struct MeshTagImpl<Dune::YaspGrid<dim, Coordinates>>
-  {
-    using type = _yaspgrid;
-  };
+  template <class Grid>
+  using MeshTag_t = typename Impl::MeshTag<Grid>::type;
   
   
+  /// A creator class for meshes. Each mesh needs different way of initialization
   template <class Grid>
   class MeshCreator
   {
-    static std::unique_ptr<Grid> create(std::string meshName)
+    static unique_ptr<Grid> create(std::string meshName)
     {
-      AMDIS_ERROR_EXIT("Not yet implemented");
+      AMDIS_ERROR_EXIT("Creator not yet implemented for this mesh type.");
     }
   };
     
@@ -71,7 +85,7 @@ namespace AMDiS
   {
     using Grid = Dune::AlbertaGrid<dim, dimworld>;
     
-    static std::unique_ptr<Grid> create(std::string meshName)
+    static unique_ptr<Grid> create(std::string meshName)
     {
       std::string macro_filename = "";
       Parameters::get(meshName + "->macro file name", macro_filename);
@@ -79,30 +93,65 @@ namespace AMDiS
       // TODO: if filename_extension is ".2d" or ".3d" read it directly from file
       // otherwise use a factory method
       
-      return std::make_unique<Grid>(macro_filename);
+      return make_unique<Grid>(macro_filename);
     }
   };
 #endif
   
   
+#if HAVE_UG
+  template <int dim>