added a boundary-element/facet-iterator implementation, to be used in the...

added a boundary-element/facet-iterator implementation, to be used in the problemStat->assemble method

dune/amdis/BoundaryElementIterator.hpp

+#pragma once
+
+namespace AMDiS
+{
+  template <MeshView>
+  class BoundaryFacetIterator;
+
+  // An Iterator over all elements and when element hasBoundaryIntersections
+  template <MeshView>
+  class BoundaryElementIterator
+  {
+    friend template <class> class BoundaryFacetIterator;
+    
+    using Element = typename MeshView::Codim<0>::Entity;
+    using ElementIterator = typename MeshView::Codim<1>::Iterator;
+    
+    class Iterator
+    {
+    public:
+      Iterator(ElementIterator elementIt, ElementIterator endIt)
+        : elementIt(elementIt)
+        , endIt(endIt)
+      {}
+      
+      Iterator(Iterator const&) = default;
+      Iterator& operator=(Iterator const&) = default;
+      
+      Iterator& operator++()
+      {
+        ++elementIt;
+        while (!elementIt->hasBoundaryIntersections() && elementIt != endIt)
+          ++elementIt;
+        return *this;
+      }
+      
+      Iterator operator++(int)
+      {
+        auto tmp = *this;
+        ++(*this);
+        return tmp;
+      }
+      
+      Element& operator*() const
+      {
+        return *elementIt;
+      }
+      
+      Element* operator->() const
+      {
+        return &(*elementIt);
+      }
+      
+      bool operator==(Iterator const& that)
+      {
+        return elementIt == that.elementIt;
+      }
+      
+      bool operator!=(Iterator const& that)
+      {
+        return !(*this == that);
+      }
+      
+    private:
+      ElementIterator elementIt;
+      ElementIterator endIt;
+    };
+
+  public:
+    /// Constructor.
+    BoundaryElementIterator(MeshView& meshView)
+      : meshView(meshView)
+    {}
+    
+    Iterator begin() {
+      auto elementIt = elements(meshView).begin();
+      auto endIt = elements(meshView).end();
+      while (!elementIt->hasBoundaryIntersections() && elementIt != endIt)
+        ++elementIt;
+      
+      return {elementIt, endIt};
+    }
+    
+    Iterator end() {
+      return {elements(meshView).end(), elements(meshView).end()};
+    }
+    
+  private:
+    MeshView& meshView;
+  };
+  
+  
+  /// Generator function for the boundary-element iterator
+  template <class MeshView>
+  BoundaryElementIterator<MeshView> boundary_elements(MeshView& 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 <MeshView>
+  class BoundaryFacetIterator
+  {
+    using Element = typename MeshView::Codim<0>::Entity;
+    using Facet   = typename MeshView::Codim<1>::Entity;
+    
+    using ElementIterator = typename BoundaryElementIterator<MeshView>::Iterator;
+    using FacetIterator = typename MeshView::IntersectionIterator;
+    
+    class Iterator
+    {
+    public:
+      Iterator(MeshView& meshView, 
+               ElementIterator elementIt, 
+               FacetIterator facetIt, 
+               ElementIterator endIt)
+        : meshView(meshView)
+        , elementIt(elementIt)
+        , facetIt(facetIt)
+        , endIt(endIt)
+      {}
+      
+      Iterator(Iterator const&) = default;
+      Iterator& operator=(Iterator const&) = default;
+      
+      Iterator& operator++()
+      {
+        ++facetIt;
+        do {
+          auto facetEndIt = intersections(meshView, *elementIt).end();
+          while (!facetIt->boundary() && facetIt != facetEndIt)
+            ++facetIt;
+          if (facetIt == facetEndIt)
+            ++elementIt;
+        } while (elementIt != endIt);
+        
+        return *this;
+      }
+      
+      Iterator operator++(int)
+      {
+        auto tmp = *this;
+        ++(*this);
+        return tmp;
+      }
+      
+      Facet& operator*() const
+      {
+        return *facetIt;
+      }
+      
+      Facet* operator->() const
+      {
+        return &(*facetIt);
+      }
+      
+      bool operator==(Iterator const& that)
+      {
+        return elementIt == that.elementIt && (elementIt == endIt || facetIt == that.facetIt);
+      }
+      
+      bool operator!=(Iterator const& that)
+      {
+        return !(*this == that);
+      }
+      
+    private:
+      MeshView&       meshView;
+      ElementIterator elementIt;
+      ElementIterator endIt;
+      FacetIterator   facetIt;
+    };
+
+  public:
+    /// Constructor.
+    BoundaryFacetIterator(MeshView& meshView)
+      : meshView(meshView)
+    {}
+    
+    Iterator begin() {
+      auto elementIt = boundary_elements(meshView).begin();
+      auto endElementIt = boundary_elements(meshView).end();
+      auto facetIt = intersections(meshView, *elementIt).begin();
+      auto endFacetIt = intersections(meshView, *elementIt).end();
+      
+      while (!facetIt->boundary() && facetIt != facetEndIt)
+        ++facetIt;
+      
+      return {elementIt, facetIt, endElementIt};
+    }
+    
+    Iterator end() {
+      auto elementIt = boundary_elements(meshView).begin();
+      auto endElementIt = boundary_elements(meshView).end();
+      auto facetIt = intersections(meshView, *elementIt).begin(); // TODO: what is the correct end_facet_iterator
+      return {endElementIt, facetIt, endElementIt};
+    }
+    
+  private:
+    MeshView& meshView;
+  };
+  
+  
+  /// Generator function for the boundary-element iterator
+  template <class MeshView>
+  BoundaryFacetIterator<MeshView> boundary_facets(MeshView& meshView)
+  {
+    return {meshView};
+  }
+
+} // end namespace AMDiS

dune/amdis/Mesh.hpp

@@ -147,7 +147,7 @@ namespace AMDiS
       Dune::FieldVector<double, dim> L; L = 1.0;  // extension of the domain
       Parameters::get(meshName + "->dimension", L);
       
-      auto s = Dune::fill_array<dim>(int{2}); // 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)      
@@ -168,7 +168,7 @@ namespace AMDiS
       Parameters::get(meshName + "->min corner", lowerleft);
       Parameters::get(meshName + "->max corner", upperright);
       
-      auto s = Dune::fill_array<dim>(int{2}); // 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)      

dune/amdis/ProblemStat.hpp

@@ -278,10 +278,8 @@ namespace AMDiS
     
   protected: // sub-methods to assemble DOFMatrix
     
-    template <class Matrix, class Vector>
-    void assemble(std::pair<int, int> row_col,
-                  Matrix& matrix, bool asmMatrix,
-                  Vector& rhs,    bool asmVector);
+    template <class LhsData, class RhsData, class Elements>
+    void assemble(LhsData lhs, RhsData rhs, Elements const& elements);
     
     template <class RowView, class ColView>
     bool getElementMatrix(RowView const& rowView,
@@ -382,6 +380,34 @@ namespace AMDiS
     VectorEntries<double*>                  vectorFactors;
     std::map< int, bool >                   vectorAssembled; // if false, do reassemble
     std::map< int, bool >                   vectorChanging;  // if true, or vectorAssembled false, do reassemble
+    
+    template <int r, int c>
+    struct MatrixData 
+    {
+      using DOFMatrixType = 
+        std::tuple_element_t<c, std::tuple_element_t<r, typename SystemMatrixType::DOFMatrices>>;
+        
+      DOFMatrixType&                       matrix;
+      std::list<shared_ptr<OperatorType>>& operators;
+      std::list<double*> const&            factors;
+      bool                                 assemble;
+      
+      std::pair<int,int> row_col = {r, c};
+    };
+    
+    template <int r>
+    struct VectorData
+    {
+      using DOFVectorType = 
+        std::tuple_element_t<r, typename SystemVectorType::DOFVectors>;
+      
+      DOFVectorType&                       vector;
+      std::list<shared_ptr<OperatorType>>& operators;
+      std::list<double*> const&            factors;
+      bool                                 assemble;
+      
+      int row = r;
+    };
   };
   
   

dune/amdis/ProblemStat.inc.hpp

@@ -277,16 +277,20 @@ namespace AMDiS
 
       For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector, _r](auto const _c) 
       {
+        using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<_r, _c>;
+        using VectorData = typename ProblemStatSeq<Traits>::template VectorData<_r>;
+        
         // The DOFMatrix which should be assembled
         auto& dofmatrix    = (*systemMatrix)(_r, _c);
         auto& solution_vec = (*solution)[_c];
         auto& rhs_vec      = (*rhs)[_r];
 	  
         auto row_col = std::make_pair(int(_r), int(_c));
-        bool asmMatrix = asmMatrix_ && (!matrixAssembled[row_col] || matrixChanging[row_col]);
+        bool assembleMatrix = asmMatrix_ && (!matrixAssembled[row_col] || matrixChanging[row_col]);
+        bool assembleVector = asmVector  && _r == _c;
         
         int r = 0, c = 0;
-	if (asmMatrix) {
+	if (assembleMatrix) {
 	  // init boundary condition
 	  for (auto bc_list : dirichletBc) {
 	    std::tie(r, c) = bc_list.first;
@@ -296,11 +300,22 @@ namespace AMDiS
 	    }
 	  }
         }
+        
+        auto mat = MatrixData{dofmatrix, matrixOperators[row_col], matrixFactors[row_col], assembleMatrix};
+        auto vec = VectorData{rhs_vec,   vectorOperators[int(_r)], vectorFactors[int(_r)], assembleVector};
 	  
         // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
-        this->assemble(row_col, dofmatrix, asmMatrix, rhs_vec, (_r == _c && asmVector));
+    
+        dofmatrix.init(assembleMatrix);
+        this->assemble(mat, vec, elements(*meshView));
+        dofmatrix.finish();
+      
+        if (assembleMatrix)
+          matrixAssembled[row_col] = true;
+        if (assembleVector)
+          vectorAssembled[int(_r)] = true;
 
-        if (asmMatrix) {
+        if (assembleMatrix) {
 	  // finish boundary condition
 	  for (auto bc_list : dirichletBc) {
 	    std::tie(r, c) = bc_list.first;
@@ -331,71 +346,51 @@ namespace AMDiS
 
 
   template <class Traits>
-    template <class Matrix, class Vector>
-  void ProblemStatSeq<Traits>::assemble(std::pair<int, int> row_col,
-                                        Matrix& dofmatrix, bool asmMatrix,
-                                        Vector& rhs,       bool asmVector)
+    template <class LhsData, class RhsData, class Elements>
+  void ProblemStatSeq<Traits>::assemble(LhsData lhs, RhsData rhs,
+                                        Elements const& elements)
   {    
-    auto const& rowFeSpace = dofmatrix.getRowFeSpace();
-    auto const& colFeSpace = dofmatrix.getColFeSpace();
-    
-    auto matrixOp = matrixOperators[row_col];
-    auto matrixOpFac = matrixFactors[row_col];
-    auto vectorOp = vectorOperators[row_col.first];
-    auto vectorOpFac = vectorFactors[row_col.first];
-    
-    dofmatrix.init(asmMatrix);
+    auto const& rowFeSpace = lhs.matrix.getRowFeSpace();
+    auto const& colFeSpace = lhs.matrix.getColFeSpace();
     
-    // return if nothing to to.
-    if ((matrixOp.empty() || !asmMatrix) && 
-        (vectorOp.empty()  || !asmVector))
-    {
-      dofmatrix.finish();
-      return;
-    }
+    if ((lhs.operators.empty() || !lhs.assemble) && 
+        (rhs.operators.empty() || !rhs.assemble))
+      return; // nothing to do
     
-    for (auto op : matrixOp)
+    for (auto op : lhs.operators)
       op->init(rowFeSpace, colFeSpace);
-    for (auto op : vectorOp)
+    for (auto op : rhs.operators)
       op->init(rowFeSpace, colFeSpace);
     
     auto rowLocalView = rowFeSpace.localView();
-    auto rowIndexSet = rowFeSpace.localIndexSet();
+    auto rowIndexSet  = rowFeSpace.localIndexSet();
     
     auto colLocalView = colFeSpace.localView();
-    auto colIndexSet = colFeSpace.localIndexSet();
+    auto colIndexSet  = colFeSpace.localIndexSet();
     
-    for (auto const& element : elements(*meshView)) {
+    for (auto const& element : elements) {
       rowLocalView.bind(element);
       rowIndexSet.bind(rowLocalView);
       
       colLocalView.bind(element);
       colIndexSet.bind(colLocalView);
       
-      if (asmMatrix) {
+      if (lhs.assemble) {
         ElementMatrix elementMatrix;
         bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix, 
-                                    matrixOp, matrixOpFac);
+                                    lhs.operators, lhs.factors);
         if (add)
-          addElementMatrix(dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
+          addElementMatrix(lhs.matrix, rowIndexSet, colIndexSet, elementMatrix);
       }
       
-      if (asmVector) {
+      if (rhs.assemble) {
         ElementVector elementVector;
         bool add = getElementVector(rowLocalView, elementVector, 
-                                    vectorOp, vectorOpFac);
+                                    rhs.operators, rhs.factors);
         if (add)
-          addElementVector(rhs, rowIndexSet, elementVector);
+          addElementVector(rhs.vector, rowIndexSet, elementVector);
       }
     }
-    
-    if (asmMatrix) {
-      dofmatrix.finish();
-      matrixAssembled[row_col] = true;
-    }
-    if (asmVector) {      
-      vectorAssembled[row_col.first] = true;
-    }
   }