SparsityPattern.hpp

#pragma once

#include <memory>

#include <dune/istl/matrixindexset.hh>

#include <amdis/common/Index.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>

namespace AMDiS
{
  class SparsityPattern
  {
  public:
    template <class Basis>
    SparsityPattern(Basis const& basis, SymmetryStructure symmetry = SymmetryStructure::unknown)
      : rows_(basis.dimension())
      , cols_(rows_)
      , pattern_(rows_, cols_)
    {
      init(basis);
    }

    template <class RowBasis, class ColBasis>
    SparsityPattern(RowBasis const& rowBasis, ColBasis const& colBasis,
                    SymmetryStructure /*symmetry*/ = SymmetryStructure::unknown)
      : rows_(rowBasis.dimension())
      , cols_(colBasis.dimension())
      , pattern_(rows_, cols_)
    {
      if (uintptr_t(&rowBasis) == uintptr_t(&colBasis))
        init(rowBasis);
      else
        init(rowBasis, colBasis);
    }

    std::size_t rows() const
    {
      return rows_;
    }

    std::size_t cols() const
    {
      return cols_;
    }

    std::size_t avgRowSize() const
    {
      assert(rows_ == pattern_.rows());
      std::size_t sum = 0;
      for (std::size_t r = 0; r < rows_; ++r)
        sum += rowSize(r);
      return (sum + rows_ - 1) / rows_;
    }

    std::size_t rowSize(std::size_t r) const
    {
      return pattern_.rowsize(r);
    }

    std::size_t nnz() const
    {
      return pattern_.size();
    }

    template <class Matrix>
    void applyTo(Matrix& matrix) const
    {
      pattern_.exportIdx(matrix);
    }


  protected:
    // Create pattern from basis. This method is called if rowBasis == colBasis.
    template <class Basis>
    void init(Basis const& basis)
    {
      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) {
          // The global index of the i-th vertex of the element
          auto row = localView.index(i);
          for (std::size_t j = 0; j < size; ++j) {
            // The global index of the j-th vertex of the element
            auto col = localView.index(j);
            pattern_.add(row, col);
          }
        }
      }
    }

    // Create pattern from basis. This method is called if rowBasis != colBasis.
    template <class RowBasis, class ColBasis>
    void init(RowBasis const& rowBasis, ColBasis const& colBasis)
    {
      auto rowLocalView = rowBasis.localView();
      auto colLocalView = colBasis.localView();
      for (const auto& element : elements(rowBasis.gridView())) {
        rowLocalView.bind(element);
        colLocalView.bind(element);

        for (std::size_t i = 0; i < rowLocalView.size(); ++i) {
          // The global index of the i-th vertex of the element
          auto row = rowLocalView.index(i);
          for (std::size_t j = 0; j < colLocalView.size(); ++j) {
            // The global index of the j-th vertex of the element
            auto col = colLocalView.index(j);
            pattern_.add(row, col);
          }
        }
      }
    }

  private:
    std::size_t rows_;
    std::size_t cols_;
    Dune::MatrixIndexSet pattern_;
  };

} // end namespace AMDiS