merged derivative term and implemented coords and functor expression

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.*"))

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

@@ -38,7 +38,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})

dune/amdis/Mesh.hpp

@@ -7,6 +7,7 @@
 #include <array>
 #include <memory>
 
+#include <dune/common/array.hh>
 #include <dune/common/fvector.hh>
 
 #include <dune/grid/albertagrid.hh>
@@ -148,10 +149,10 @@ namespace AMDiS
 
       msg("L = ", L);
 
-      std::array<int, dim> s; // number of cells on coarse mesh in each direction
+      auto s = Dune::fill_array<int,dim>(2); // number of cells on coarse mesh in each direction
       Parameters::get(meshName + "->num cells", s);
 
-      // TODO: add more parameters for yasp-grid (see constructor)
+      // TODO: add more parameters for yasp-grid (see constructor)      
 
       return make_unique<Grid>(L, s);
     }
@@ -165,15 +166,15 @@ namespace AMDiS
 
     static unique_ptr<Grid> create(std::string meshName)
     {
-      Dune::FieldVector<double, dim> lowerleft;  // Lower left corner of the domain
-      Dune::FieldVector<double, dim> upperright; // Upper right corner of the domain
+      Dune::FieldVector<double, dim> lowerleft;  lowerleft = 0.0; // Lower left corner of the domain
+      Dune::FieldVector<double, dim> upperright; upperright = 1.0; // Upper right corner of the domain
       Parameters::get(meshName + "->min corner", lowerleft);
       Parameters::get(meshName + "->max corner", upperright);
 
-      std::array<int, dim> s; // number of cells on coarse mesh in each direction
+      auto s = Dune::fill_array<int,dim>(2); // number of cells on coarse mesh in each direction
       Parameters::get(meshName + "->num cells", s);
 
-      // TODO: add more parameters for yasp-grid (see constructor)
+      // TODO: add more parameters for yasp-grid (see constructor)      
 
       return make_unique<Grid>(lowerleft, upperright, s);
     }

dune/amdis/Operator.hpp

@@ -3,6 +3,8 @@
 #include <list>
 #include <vector>
 
+#include <dune/geometry/type.hh>
+
 #include <dune/amdis/OperatorTermBase.hpp>
 #include <dune/amdis/terms/TermGenerator.hpp>
 
@@ -16,8 +18,8 @@ namespace AMDiS
     GRD_PSI,
     GRD_PHI
   };
-    
-    
+
+
   template <class MeshView>
   class Operator
   {
@@ -26,10 +28,10 @@ namespace AMDiS
 
     using ElementVector = Impl::ElementVector;
     using ElementMatrix = Impl::ElementMatrix;
-    
+
     using IdType = typename MeshView::Grid::LocalIdSet::IdType;
-    
-  public:    
+
+  public:
     /// Add coefficients for zero-order operator < c * u, v >.
     /// The coefficient \p c must be a scalar expression.
     template <class Term>
@@ -37,7 +39,7 @@ namespace AMDiS
     {
       return addZOTImpl(toTerm(c));
     }
-    
+
     template <class... Args>
     static shared_ptr<Self> zot(Args&&... args)
     {
@@ -45,10 +47,10 @@ namespace AMDiS
       op->addZOT(std::forward<Args>(args)...);
       return op;
     }
-    
-    
-    /// Add coefficients for first-order operator < psi, term * grad(phi) > or 
-    /// < grad(psi), term * phi >, where the first corresponds to 
+
+
+    /// Add coefficients for first-order operator < psi, term * grad(phi) > or
+    /// < grad(psi), term * phi >, where the first corresponds to
     /// \p firstOrderType == GRD_PHI and the second one to \p firstOrderType == GRD_PSI.
     /// The coefficient \p b must be a vector expression.
     template <class Term>
@@ -56,9 +58,9 @@ namespace AMDiS
     {
       return addFOTImpl(toTerm(b), firstOrderType);
     }
-    
-    /// Add coefficients for first-order operator < psi, b * d_i(phi) > or 
-    /// < d_i(psi), b * phi >, where the first corresponds to 
+
+    /// Add coefficients for first-order operator < psi, b * d_i(phi) > or
+    /// < d_i(psi), b * phi >, where the first corresponds to
     /// \p firstOrderType == GRD_PHI and the second one to \p firstOrderType == GRD_PSI.
     /// The coefficient \p b must be a scalar expression.
     template <class Term, size_t I>
@@ -66,7 +68,7 @@ namespace AMDiS
     {
       return addFOTImpl(toTerm(b), _i, firstOrderType);
     }
-    
+
     template <class... Args>
     static shared_ptr<Self> fot(Args&&... args)
     {
@@ -74,16 +76,16 @@ namespace AMDiS
       op->addFOT(std::forward<Args>(args)...);
       return op;
     }
-    
-    
+
+
     /// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
-    /// where \p A can be a matrix or a scalar expression. 
+    /// where \p A can be a matrix or a scalar expression.
     template <class Term>
     Self& addSOT(Term const& A)
     {
       return addSOTImpl(toTerm(A));
     }
-    
+
     /// Add coefficients for second-order operator < d_i(psi), A * d_j(phi) >,
     /// where \p A can be a matrix or a scalar expression. The first index \p _i
     /// corresponds to the derivative component of the test-function and the
@@ -93,7 +95,7 @@ namespace AMDiS
     {
       return addSOTImpl(toTerm(A), _i, _j);
     }
-    
+
     template <class... Args>
     static shared_ptr<Self> sot(Args&&... args)
     {
@@ -101,110 +103,111 @@ namespace AMDiS
       op->addSOT(std::forward<Args>(args)...);
       return op;
     }
-    
-    
-    /// Calls \ref zot(), \ref for() or \ref sot(), depending on template 
+
+
+    /// Calls \ref zot(), \ref for() or \ref sot(), depending on template
     /// parameter \p Order.
     template <size_t Order, class... Args>
     static shared_ptr<Self> create(Args&&... args)
     {
       return create(index_<Order>{}, std::forward<Args>(args)...);
     }
-    
-    
+
+
     /// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
     template <class RowFeSpace, class ColFeSpace>
     void init(RowFeSpace const& rowFeSpace,
               ColFeSpace const& colFeSpace);
-    
+
 
     /// Calculates the needed quadrature degree for the given term-order \p order.
     /// This takes into account the polynomial degree of trial and test functions
     /// and the polynomial degree of the coefficient functions
-    int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);
-      
-    
+    template <class Element>
+    int getQuadratureDegree(Element const& element, int order, FirstOrderType firstOrderType = GRD_PHI);
+
+
     template <class RowView, class ColView>
     bool getElementMatrix(RowView const& rowView,
                           ColView const& colView,
                           ElementMatrix& elementMatrix,
                           double* factor = NULL);
-    
-    
+
+
     template <class RowView>
     bool getElementVector(RowView const& rowView,
                           ElementVector& elementVector,
                           double* factor = NULL);
-    
+
   protected:
     template <class RowView, class ColView>
     void assembleZeroOrderTerms(RowView const& rowView,
 				ColView const& colView,
 				ElementMatrix& elementMatrix);
-    
+
     template <class RowView>
     void assembleZeroOrderTerms(RowView const& rowView,
 				ElementVector& elementVector);
-    
+
     template <class RowView, class ColView>
     void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
 				       ColView const& colView,
 				       ElementMatrix& elementMatrix);
-    
+
     template <class RowView, class ColView>
     void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
 				       ColView const& colView,
 				       ElementMatrix& elementMatrix);
-    
+
     template <class RowView>
     void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
                                        ElementVector& elementVector);
-    
+
     template <class RowView, class ColView>
     void assembleSecondOrderTerms(RowView const& rowView,
 				  ColView const& colView,
 				  ElementMatrix& elementMatrix);
-    
+
   private:
     template <class Term>
     Self& addZOTImpl(Term const& term);
-    
+
     template <class Term>
     Self& addFOTImpl(Term const& term, FirstOrderType firstOrderType);
-    
+
     template <class Term, size_t I>
-    Self& addFOTImpl(Term const& term, const index_<I>, 
+    Self& addFOTImpl(Term const& term, const index_<I>,
                      FirstOrderType firstOrderType);
-    
+
     template <class Term>
     Self& addSOTImpl(Term const& term);
-    
+
     template <class Term, size_t I, size_t J>
     Self& addSOTImpl(Term const& term, const index_<I>, const index_<J>);
-    
-    
+
+
     template <class... Args>
     static shared_ptr<Self> create(index_<0>, Args&&... args)
     {
       return zot(std::forward<Args>(args)...);
     }
-    
+
     template <class... Args>
     static shared_ptr<Self> create(index_<1>, Args&&... args)
     {
       return fot(std::forward<Args>(args)...);
     }
-    
+
     template <class... Args>
     static shared_ptr<Self> create(index_<2>, Args&&... args)
     {
       return sot(std::forward<Args>(args)...);
     }
-      
+
   private:
     template <class LocalView>
     IdType getElementId(LocalView const& localView);
-      
+
   private:
     /// List of all second order terms
     std::list<OperatorTermType*> secondOrder;
@@ -217,17 +220,17 @@ namespace AMDiS
 
     /// List of all zero order terms
     std::list<OperatorTermType*> zeroOrder;
-    
+
     int psiDegree = 1;
     int phiDegree = 1;
-    
+
     IdType lastMatrixId = 0;
     IdType lastVectorId = 0;
-    
+
     ElementMatrix cachedElementMatrix;
     ElementVector cachedElementVector;
   };
-  
+
 } // end namespace AMDiS
 
 #include "Operator.inc.hpp"

dune/amdis/Operator.inc.hpp

@@ -27,10 +27,14 @@ namespace AMDiS
 
 
   template <class MeshView>
-  int Operator<MeshView>::getQuadratureDegree(int order, FirstOrderType firstOrderType)
+    template <class Element>
+  int Operator<MeshView>::getQuadratureDegree(Element const& element, int order, FirstOrderType firstOrderType)
   {
     std::list<OperatorTermType*>* terms = NULL;
 
+    auto const& t = element.type();
+    auto geometry = element.geometry();
+
     switch(order)
     {
     case 0:
@@ -49,9 +53,16 @@ namespace AMDiS
     int maxTermDegree = 0;
 
     for (OperatorTermType* term : *terms)
-      maxTermDegree = std::max(maxTermDegree, term->getDegree());
+      maxTermDegree = std::max(maxTermDegree, term->getDegree(t));
+
+    int degree = psiDegree + phiDegree + maxTermDegree;
+    if (t.isSimplex())
+      degree -= order;
+
+    if (!geometry.affine())
+      degree += 1; // oder += (order+1)
 
-    return psiDegree + phiDegree - order + maxTermDegree;
+    return degree;
   }
 
 
@@ -177,7 +188,7 @@ namespace AMDiS
     auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
     auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
 
-    int order = getQuadratureDegree(0);
+    int order = getQuadratureDegree(element, 0);
     auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
 
     for (auto* operatorTerm : zeroOrder)
@@ -226,7 +237,7 @@ namespace AMDiS
 
     auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
 
-    int order = getQuadratureDegree(0);
+    int order = getQuadratureDegree(element, 0);
     auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
 
     for (auto* operatorTerm : zeroOrder)
@@ -269,7 +280,7 @@ namespace AMDiS
     auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
     auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
 
-    int order = getQuadratureDegree(2);
+    int order = getQuadratureDegree(element, 1, GRD_PHI);
     auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
 
     for (auto* operatorTerm : firstOrderGrdPhi)
@@ -329,7 +340,7 @@ namespace AMDiS
     auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
     auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
 
-    int order = getQuadratureDegree(