moved the sparsity pattern update method to the bilinear form

CMakeLists.txt

@@ -23,12 +23,13 @@ add_subdirectory("src")
 add_subdirectory("test")
 
 target_link_libraries(amdis fmt)
+target_compile_options(amdis PUBLIC "-Wfatal-errors" "-Wall" "-Wpedantic")
 
 option(ENABLE_ALL_WARNINGS "enable all meaningful warnings" OFF)
 if (ENABLE_ALL_WARNINGS)
-  target_compile_options(amdis PUBLIC "-Wall" "-Wextra" "-pedantic" "-Wnon-virtual-dtor"
-                                      "-Wold-style-cast" "-Wcast-align" "-Woverloaded-virtual"
-                                      "-Wpedantic" "-Wconversion")
+  target_compile_options(amdis PUBLIC "-Wextra" "-Wnon-virtual-dtor"
+                                      "-Wold-style-cast" "-Wcast-align"
+                                      "-Woverloaded-virtual" "-Wconversion")
 endif (ENABLE_ALL_WARNINGS)
 
 # finalize the dune project, e.g. generating config.h etc.

src/amdis/BiLinearForm.hpp

@@ -3,6 +3,7 @@
 #include <cmath>
 
 #include <amdis/LinearAlgebra.hpp>
+#include <amdis/Observer.hpp>
 #include <amdis/OperatorList.hpp>
 #include <amdis/common/FlatMatrix.hpp>
 #include <amdis/common/TypeTraits.hpp>
@@ -22,10 +23,11 @@ namespace AMDiS
    **/
   template <class RB, class CB, class T = double, class Traits = BackendTraits<RB,T>>
   class BiLinearForm
-      : public MatrixFacade<T, typename Traits::SparsityPattern, Traits::template MatrixImpl>
+      : public MatrixFacade<T, Traits::template MatrixImpl>
+      , private ObserverSequence<event::adapt,2>
   {
     using Self = BiLinearForm;
-    using Super = MatrixFacade<T, typename Traits::SparsityPattern, Traits::template MatrixImpl>;
+    using Super = MatrixFacade<T, Traits::template MatrixImpl>;
 
   public:
     /// The type of the finite element space / basis of the row
@@ -42,13 +44,18 @@ namespace AMDiS
     /// The type of the matrix filled on an element with local contributions
     using ElementMatrix = FlatMatrix<CoefficientType>;
 
+    /// Type of the sparsity pattern of the backend
+    using SparsityPattern = typename Traits::SparsityPattern;
+
   public:
     /// Constructor. Stores the row and column basis in a local `shared_ptr` to const
     BiLinearForm(std::shared_ptr<RB> const& rowBasis, std::shared_ptr<CB> const& colBasis)
       : Super(*rowBasis, *colBasis)
+      , ObserverSequence<event::adapt,2>(*rowBasis, *colBasis)
       , rowBasis_(rowBasis)
       , colBasis_(colBasis)
     {
+      pattern_.init(*rowBasis_, *colBasis_);
       operators_.init(*rowBasis_, *colBasis_);
 
       auto const rowSize = rowBasis_->localView().maxSize();
@@ -128,7 +135,34 @@ namespace AMDiS
     /// Assemble all matrix operators, TODO: incorporate boundary conditions
     void assemble(SymmetryStructure symmetry = SymmetryStructure::unknown);
 
+    void init(SymmetryStructure symmetry = SymmetryStructure::unknown)
+    {
+      Super::init(pattern_, symmetry);
+    }
+
+    using Super::init;
+
+    void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(e,i); }
+    void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(e,i); }
+
+  private:
+    // Track for each basis wether updated and if both are, update the sparsity pattern
+    template <std::size_t I>
+    void updateImpl2(event::adapt, index_t<I>)
+    {
+      assert(!updateCounter_.test(I));
+      updateCounter_.set(I);
+
+      if (updateCounter_.all()) {
+        pattern_.init(*rowBasis_, *colBasis_);
+        updateCounter_.reset();
+      }
+    }
+
   protected:
+    /// The structure of the non-zeros in the matrix
+    SparsityPattern pattern_;
+
     /// Dense matrix to store coefficients during \ref assemble()
     ElementMatrix elementMatrix_;
 
@@ -137,6 +171,9 @@ namespace AMDiS
 
     std::shared_ptr<RowBasis const> rowBasis_;
     std::shared_ptr<ColBasis const> colBasis_;
+
+  private:
+    std::bitset<2> updateCounter_ = 0;
   };
 
 

src/amdis/DOFVector.hpp

@@ -145,6 +145,19 @@ namespace AMDiS
     /// Read backup data from file
     void restore(std::string const& filename);
 
+    void resize()
+    {
+      Coefficients::resize(sizeInfo(*basis_));
+    }
+
+    using Coefficients::resize;
+
+    void resizeZero()
+    {
+      Coefficients::resizeZero(sizeInfo(*basis_));
+    }
+
+    using Coefficients::resizeZero;
 
     /// Return the associated DataTransfer object
     DataTransfer const& dataTransfer() const
@@ -189,7 +202,7 @@ namespace AMDiS
     void updateImpl(event::adapt e) override
     {
       assert(e.value);
-      this->resize();
+      resize();
       dataTransfer_.adapt(*this);
     }
 

src/amdis/DOFVector.inc.hpp

@@ -48,7 +48,7 @@ restore(std::string const& filename)
   // assume the order of element traversal is not changed
   auto localView = this->basis()->localView();
   std::vector<value_type> data;
-  this->init(true);
+  this->init(sizeInfo(*this->basis()), true);
   for (auto const& element : elements(this->basis()->gridView()))
   {
     std::uint64_t len = 0;

src/amdis/DataTransfer.inc.hpp

@@ -148,7 +148,7 @@ preAdapt(C const& coeff, bool mightCoarsen)
 template <class C, class B>
 void DataTransfer<C,B>::adapt(C& coeff)
 {
-  coeff.init(false);
+  coeff.resize();
 
   test_exit(!persistentContainer_.empty(), "No data was saved before adapting the grid, make "
     "sure to call DataTransfer::preAdapt before calling adapt() on the grid");

src/amdis/Initfile.hpp

@@ -20,17 +20,15 @@ namespace AMDiS
      * Looks for the `key` in the parameter-tree and returns
      * the stored and parsed value if found and parsable.
      *
-     * Does not thrown an exception if something goes wrong!
+     * May throw an exception if the value can not be parsed into type T
      **/
     template <class T>
     static Dune::Std::optional<T> get(std::string const& key)
     {
-      try {
+      if (pt().hasKey(key))
         return pt().get<T>(key);
-      }
-      catch (...) {
+      else
         return {};
-      }
     }
 
     /// \brief Get parameter-values from parameter-tree with default value

src/amdis/LinearForm.inc.hpp

@@ -59,7 +59,7 @@ assemble()
 {
   auto localView = this->basis()->localView();
 
-  this->init(true);
+  this->init(sizeInfo(*this->basis()), true);
   for (auto const& element : elements(this->basis()->gridView(), typename Traits::PartitionSet{})) {
     localView.bind(element);
     this->assemble(localView);

src/amdis/Observer.hpp

@@ -132,7 +132,7 @@ namespace AMDiS
       return *this;
     }
 
-    /// Set the Notifier* to nullptr. Used by the Notifer to avoid segfaults when destruction occurs
+    /// Set the Notifier* to nullptr. Used by the Notifier to avoid segfaults when destruction occurs
     /// out of order.
     void unset() final
     {

src/amdis/ProblemStat.inc.hpp

@@ -167,7 +167,7 @@ restore(Flag initFlag)
   if (initFlag.isSet(INIT_FILEWRITER))
     createFileWriter();
 
-  solution_->resize();
+  solution_->resize(sizeInfo(*globalBasis_));
   solution_->restore(solution_filename);
 }
 
@@ -516,7 +516,7 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as
   Parameters::get(name_ + "->symmetry", symmetryStr);
 
   systemMatrix_->init(symmetryStructure(symmetryStr));
-  rhs_->init(asmVector);
+  rhs_->init(sizeInfo(*globalBasis_), asmVector);
 
   // statistic about system size
   if (Environment::mpiSize() > 1)
@@ -543,7 +543,7 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as
   info(2,"  assemble operators needed {} seconds", t2.elapsed());
   t2.reset();
 
-  solution_->resize();
+  solution_->resize(sizeInfo(*globalBasis_));
   for_each_node(localView.tree(), [&,this](auto const& rowNode, auto row_tp) -> void {
     for_each_node(localView.tree(), [&,this](auto const& colNode, auto col_tp) -> void {
       // finish boundary condition

src/amdis/common/FakeContainer.hpp

@@ -60,7 +60,8 @@ namespace AMDiS
     template <class... Args>
     constexpr explicit FakeContainer(Args&&...) noexcept {}
 
-    constexpr void init(bool) noexcept { /* do nothing */ }
+    template <class S>
+    constexpr void init(S const&, bool) noexcept { /* do nothing */ }
     constexpr void finish() noexcept { /* do nothing */ }
 
 

src/amdis/functions/Interpolate.hpp

@@ -53,8 +53,8 @@ namespace AMDiS
       // Obtain a local view of f
       auto lf = localFunction(gf);
 
-      vector.init(false);
-      counter.init(true); // set to zero
+      vector.init(sizeInfo(basis), false);
+      counter.init(sizeInfo(basis), true); // set to zero
       for (const auto& e : elements(basis.gridView(), typename BackendTraits<B>::PartitionSet{}))
       {
         localView.bind(e);

src/amdis/linearalgebra/MatrixFacade.hpp

@@ -22,22 +22,18 @@ namespace AMDiS
    * \tparam Pattern    The type of the sparsity pattern
    * \tparam MatrixImpl A linear-algebra backend for the matrix storage
    **/
-  template <class T, class Pattern, template <class> class MatrixImpl>
+  template <class T, template <class> class MatrixImpl>
   class MatrixFacade
   {
     using Self = MatrixFacade;
     using Impl = MatrixImpl<T>;
 
-    /// Type of the sparsity pattern of the backend
-    using SparsityPattern = Pattern;
-
   public:
-    /// Constructor. Forwards the bases to the implementation class and 
+    /// Constructor. Forwards the bases to the implementation class and
     /// constructs a matrix sparsity pattern.
     template <class RowBasis, class ColBasis>
     MatrixFacade(RowBasis const& rowBasis, ColBasis const& colBasis)
       : impl_(rowBasis, colBasis)
-      , pattern_(rowBasis, colBasis)
     {}
 
     /// Return the underlying matrix backend
@@ -50,9 +46,10 @@ namespace AMDiS
      * structure of the values or the sparsity pattern can be provided.
      * See \ref SymmetryStructure.
      **/
-    void init(SymmetryStructure symmetry = SymmetryStructure::unknown)
+    template <class SparsityPattern>
+    void init(SparsityPattern const& pattern, SymmetryStructure symmetry = SymmetryStructure::unknown)
     {
-      impl_.init(pattern_, symmetry);
+      impl_.init(pattern, symmetry);
     }
 
     /// Finish the matrix insertion, e.g. cleanup or final insertion
@@ -99,9 +96,6 @@ namespace AMDiS
   protected:
     /// The matrix backend
     Impl impl_;
-
-    /// The structure of the non-zeros in the matrix
-    SparsityPattern pattern_;
   };
 
 } // end namespace AMDiS

src/amdis/linearalgebra/SparsityPattern.hpp

@@ -4,27 +4,15 @@
 
 #include <dune/istl/matrixindexset.hh>
 
-#include <amdis/Observer.hpp>
 #include <amdis/common/Index.hpp>
 
 namespace AMDiS
 {
   /// \brief A general sparsity pattern implementation using the full pattern of the
   /// basis by adding all local indices
-  template <class RowBasis, class ColBasis>
   class SparsityPattern
-      : private ObserverSequence<event::adapt,2>
   {
   public:
-    /// Constructor. Stores pointers to the passed bases.
-    SparsityPattern(RowBasis const& rowBasis, ColBasis const& colBasis)
-      : ObserverSequence<event::adapt,2>(rowBasis, colBasis)
-      , rowBasis_(&rowBasis)
-      , colBasis_(&colBasis)
-    {
-      updateImpl3();
-    }
-
     /// Number of rows in the matrix
     std::size_t rows() const
     {
@@ -66,42 +54,17 @@ namespace AMDiS
       pattern_.exportIdx(matrix);
     }
 
-  protected:
-
-    void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(e,i); }
-    void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(e,i); }
-
-  private:
-
-    // Track for each basis wether updated and if both are, call updateImpl3()
-    template <std::size_t I>
-    void updateImpl2(event::adapt, index_t<I>)
-    {
-      assert(!updateCounter_.test(I));
-      updateCounter_.set(I);
-
-      if (updateCounter_.all()) {
-        updateImpl3();
-        updateCounter_.reset();
-      }
-    }
-
-    void updateImpl3()
-    {
-      rowBasis_ == colBasis_
-          ? updateSameBasis()
-          : updateDifferentBasis();
-    }
-
     // Update pattern when basis is updated. This method is called if rowBasis == colBasis.
-    void updateSameBasis()
+    template <class Basis>
+    void init(Basis const& basis)
     {
-      rows_ = rowBasis_->dimension();
+      rows_ = basis.dimension();
       cols_ = rows_;
+      pattern_.resize(0, 0); // clear the old pattern
       pattern_.resize(rows_, cols_);
 
-      auto localView = rowBasis_->localView();
-      for (const auto& element : elements(rowBasis_->gridView())) {
+      auto localView = basis.localView();
+      for (const auto& element : elements(basis.gridView())) {
         localView.bind(element);
 
         for (std::size_t i = 0, size = localView.size(); i < size; ++i) {
@@ -117,15 +80,20 @@ namespace AMDiS
     }
 
     // Update pattern when basis is updated. This method is called if rowBasis != colBasis.
-    void updateDifferentBasis()
+    template <class RowBasis, class ColBasis>
+    void init(RowBasis const& rowBasis, ColBasis const& colBasis)
     {
-      rows_ = rowBasis_->dimension();
-      cols_ = colBasis_->dimension();
+      if (uintptr_t(&rowBasis) == uintptr_t(&colBasis))
+        return init(rowBasis);
+
+      rows_ = rowBasis.dimension();
+      cols_ = colBasis.dimension();
+      pattern_.resize(0, 0); // clear the old pattern
       pattern_.resize(rows_, cols_);
 
-      auto rowLocalView = rowBasis_->localView();
-      auto colLocalView = colBasis_->localView();
-      for (const auto& element : elements(rowBasis_->gridView())) {
+      auto rowLocalView = rowBasis.localView();
+      auto colLocalView = colBasis.localView();
+      for (const auto& element : elements(rowBasis.gridView())) {
         rowLocalView.bind(element);
         colLocalView.bind(element);
 
@@ -142,11 +110,6 @@ namespace AMDiS
     }
 
   private:
-    RowBasis const* rowBasis_;
-    ColBasis const* colBasis_;
-
-    std::bitset<2> updateCounter_ = 0;
-
     std::size_t rows_;
     std::size_t cols_;
     Dune::MatrixIndexSet pattern_;

src/amdis/linearalgebra/VectorFacade.hpp

@@ -18,7 +18,6 @@
 #include <amdis/common/FakeContainer.hpp>
 #include <amdis/common/TypeTraits.hpp>
 #include <amdis/functions/NodeIndices.hpp>
-#include <amdis/functions/SizeInfo.hpp>
 #include <amdis/operations/Assigner.hpp>
 #include <amdis/typetree/MultiIndex.hpp>
 
@@ -57,15 +56,14 @@ namespace AMDiS
     using size_type = typename Impl::size_type;
     using value_type = typename Impl::value_type;
 
-    /// Constructor. Forwards the basis to the implementation class and 
+    /// Constructor. Forwards the basis to the implementation class and
     /// constructs a (type-erased) size-info object.
     template <class GlobalBasis,
       class = void_t<decltype(std::declval<GlobalBasis>().dimension())> >
     VectorFacade(GlobalBasis const& basis)
       : impl_(basis)
-      , sizeInfo_(sizeInfo(basis))
     {
-      resizeZero();
+      resizeZero(sizeInfo(basis));
     }
 
     /// Return the underlying linear algebra backend
@@ -96,22 +94,25 @@ namespace AMDiS
     }
 
     /// Resize the \ref vector to the size of the \ref basis
-    void resize()
+    template <class SizeInfo>
+    void resize(SizeInfo const& sizeInfo)
     {
-      init(false);
+      init(sizeInfo, false);
     }
 
     /// Resize the \ref vector to the size of the \ref basis and set to zero
-    void resizeZero()
+    template <class SizeInfo>
+    void resizeZero(SizeInfo const& sizeInfo)
     {
-      init(true);
+      init(sizeInfo, true);
     }
 
     /// Prepare the vector for insertion of values, finish the insertion with
     /// \ref finish().
-    void init(bool clear)
+    template <class SizeInfo>
+    void init(SizeInfo const& sizeInfo, bool clear)
     {
-      impl_.init(sizeInfo_, clear);
+      impl_.init(sizeInfo, clear);
       state_ = clear ? VectorState::synchronized : VectorState::unknown;
     }
 
@@ -341,9 +342,6 @@ namespace AMDiS
     /// Data backend
     Impl impl_;
 
-    /// Functor returning the basis dimension
-    SizeInfo sizeInfo_;
-
     /// The current state of the vector, one of {synchronized, insert_values,
     /// add_values, unknown}
     VectorState state_ = VectorState::unknown;

src/amdis/linearalgebra/eigen/MatrixSize.hpp

@@ -3,35 +3,13 @@
 #include <functional>
 #include <memory>
 
-#include <amdis/Observer.hpp>
 #include <amdis/common/Index.hpp>
 
 namespace AMDiS
 {
   class MatrixSize
-      : private ObserverSequence<event::adapt,2>
   {
-    struct BasisCallback
-    {
-      template <class Basis>
-      explicit BasisCallback(Basis const* basis)
-        : dimension([basis](){ return basis->dimension(); })
-      {}
-
-      std::function<std::size_t()> dimension;
-    };
-
   public:
-    template <class RowBasis, class ColBasis>
-    MatrixSize(RowBasis const& rowBasis, ColBasis const& colBasis)
-      : ObserverSequence<event::adapt,2>(rowBasis, colBasis)
-      , rowBasis_(&rowBasis)
-      , colBasis_(&colBasis)
-    {
-      updateImpl(event::adapt{true}, index_t<0>{});
-      updateImpl(event::adapt{true}, index_t<1>{});
-    }
-
     /// Number of rows in the matrix
     std::size_t rows() const
     {
@@ -44,17 +22,16 @@ namespace AMDiS
       return cols_;
     }
 
-  protected:
-
-    void updateImpl(event::adapt e, index_t<0> i) final { rows_ = rowBasis_.dimension(); }
-    void updateImpl(event::adapt e, index_t<1> i) final { cols_ = colBasis_.dimension(); }
+    template <class RowBasis, class ColBasis>
+    void init(RowBasis const& rowBasis, ColBasis const& colBasis)
+    {
+      rows_ = rowBasis.dimension();
+      cols_ = colBasis.dimension();
+    }
 
   private:
-    BasisCallback rowBasis_;
-    BasisCallback colBasis_;
-
-    std::size_t rows_;
-    std::size_t cols_;
+    std::size_t rows_ = 0;
+    std::size_t cols_ = 0;
   };
 
 } // end namespace AMDiS

src/amdis/linearalgebra/istl/Traits.hpp

@@ -24,7 +24,7 @@ namespace AMDiS
   {
     template <class Value>
     using MatrixImpl = ISTLBCRSMatrix<Value>;
-    
+
     template <class Value>
     using VectorImpl = ISTLBlockVector<Value>;
 
@@ -34,7 +34,7 @@ namespace AMDiS
 
     using CoefficientType = T;
     using PartitionSet    = Dune::Partitions::All;
-    using SparsityPattern = AMDiS::SparsityPattern<Basis, Basis>;
+    using SparsityPattern = AMDiS::SparsityPattern;
 
     using ScalProd        = Dune::ScalarProduct<Vec>;
     using LinOp           = Dune::AssembledLinearOperator<Mat, Vec, Vec>;

src/amdis/linearalgebra/mtl/SlotSize.hpp

@@ -12,32 +12,8 @@
 namespace AMDiS
 {
   class SlotSize
-      : private ObserverSequence<event::adapt,2>
   {
-    struct BasisCallback
-    {
-      template <class Basis>
-      explicit BasisCallback(Basis const* basis)
-        : dimension([basis](){ return basis->dimension(); })
-        , maxSize([basis](){ return basis->localView().maxSize(); })
-      {}
-
-      std::function<std::size_t()> dimension;
-      std::function<std::size_t()> maxSize;
-    };
-
   public:
-    template <class RowBasis, class ColBasis>
-    SlotSize(RowBasis const& rowBasis, ColBasis const& colBasis)
-      : ObserverSequence<event::adapt,2>(rowBasis, colBasis)
-      , rowBasis_(&rowBasis)
-      , colBasis_(&colBasis)
-      , surrounding_(Math::pow<RowBasis::GridView::dimension>(2) * Dune::factorial(int(RowBasis::GridView::dimension)))
-    {
-      updateImpl(event::adapt{true}, index_t<0>{});
-      updateImpl(event::adapt{true}, index_t<1>{});
-    }
-
     /// Number of rows in the matrix
     std::size_t rows() const
     {
@@ -52,7 +28,7 @@ namespace AMDiS
 
     /// Average number of non-zeros per row
     /// In the first time this function is called, use the estRowSize,
-    /// otherwise take `1.5 * nnz/rows`
+    /// otherwise take `1.2 * nnz/rows`
     template <class Matrix>
     std::size_t avgRowSize(Matrix const& matrix) const
     {
@@ -60,31 +36,25 @@ namespace AMDiS
       return matrix.nnz() > 0 ? 6*matrix.nnz() / (5*rows_) : estRowSize_;
     }
 
-  protected:
-
-    // update of the row basis, just update the dimension
-    void updateImpl(event::adapt e, index_t<0> i) final
-    {
-      rows_ = rowBasis_.dimension();
-    }
-
     // estimate the number of columns by multiplying the maximal node size with the
     // number of element surrounding a vertex. This number is approximated by the
-    // number of simplices surrounding a vertex in a kuhn tringulation
-    void updateImpl(event::adapt e, index_t<1> i) final