unit_filter.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/lafem/sparse_matrix_csr.hpp>
#include <kernel/lafem/sparse_matrix_bcsr.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/sparse_vector.hpp>
#include <kernel/lafem/arch/unit_filter.hpp>
 
namespace FEAT
{
  namespace LAFEM
  {
    template<
      typename DT_,
      typename IT_ = Index>
    class UnitFilter
    {
    public:
      typedef DT_ DataType;
      typedef IT_ IndexType;
 
      typedef DenseVector<DataType, IndexType> VectorType;
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      using FilterType = UnitFilter<DT2_, IT2_>;
 
      template <typename DataType2_, typename IndexType2_>
      using FilterTypeByDI = FilterType<DataType2_, IndexType2_>;
 
      static constexpr bool is_global = false;
      static constexpr bool is_local = true;
 
    private:
      SparseVector<DT_, IT_> _sv;
 
    public:
      UnitFilter() :
        _sv()
      {
      }
 
      explicit UnitFilter(Index size_in) :
        _sv(size_in)
      {
      }
 
      explicit UnitFilter(Index size_in, DenseVector<DT_, IT_> & values, DenseVector<IT_, IT_> & indices) :
        _sv(size_in, values, indices)
      {
        XASSERTM(values.size() == indices.size(), "Vector size mismatch!");
      }
 
      UnitFilter(UnitFilter && other) :
        _sv(std::move(other._sv))
      {
      }
 
      UnitFilter & operator=(UnitFilter && other)
      {
        if(this != &other)
        {
          _sv = std::forward<decltype(other._sv)>(other._sv);
        }
        return *this;
      }
 
      virtual ~UnitFilter()
      {
      }
 
      UnitFilter clone(CloneMode clone_mode = CloneMode::Deep) const
      {
        UnitFilter other;
        other.clone(*this, clone_mode);
        return other;
      }
 
      void clone(const UnitFilter & other, CloneMode clone_mode = CloneMode::Deep)
      {
        _sv.clone(other.get_filter_vector(), clone_mode);
      }
 
      template<typename DT2_, typename IT2_>
      void convert(const UnitFilter<DT2_, IT2_>& other)
      {
        _sv.convert(other.get_filter_vector());
      }
 
      void clear()
      {
        _sv.clear();
      }
 
      std::size_t bytes() const
      {
        return _sv.bytes();
      }
 
      SparseVector<DT_, IT_>& get_filter_vector()
      {
        return _sv;
      }
 
      const SparseVector<DT_, IT_>& get_filter_vector() const
      {
        return _sv;
      }
 
      void add(IndexType idx, DataType val)
      {
        _sv(idx, val);
      }
 
      Index size() const
      {
        return _sv.size();
      }
 
      Index used_elements() const
      {
        return _sv.used_elements();
      }
 
      IT_* get_indices()
      {
        return _sv.indices();
      }
 
      const IT_* get_indices() const
      {
        return _sv.indices();
      }
 
      DT_* get_values()
      {
        return _sv.elements();
      }
 
      const DT_* get_values() const
      {
        return _sv.elements();
      }
 
      void permute(Adjacency::Permutation & perm)
      {
        _sv.permute(perm);
      }
 
#ifdef DOXYGEN
      // The following documentation block is visible to Doxygen only. The actual implementation is matrix type
      // specific and provided below.
 
      void filter_mat(MatrixType & matrix) const
      {
      }
 
      void filter_offdiag_row_mat(MatrixType & matrix) const
      {
      }
 
      void filter_offdiag_col_mat(MatrixType & matrix) const
      {
      }
 
#endif
      void filter_mat(SparseMatrixCSR<DT_, IT_> & matrix) const
      {
        if(_sv.used_elements() == Index(0))
          return;
 
        XASSERTM(_sv.size() == matrix.rows(), "Matrix size does not match!");
 
        const IndexType* row_ptr(matrix.row_ptr());
        const IndexType* col_idx(matrix.col_ind());
        Arch::UnitFilter::filter_unit_mat(matrix.val(), row_ptr, col_idx, _sv.indices(), _sv.used_elements());
      }
 
      void filter_offdiag_row_mat(SparseMatrixCSR<DT_, IT_> & matrix) const
      {
        if(_sv.used_elements() == Index(0))
          return;
 
        XASSERTM(_sv.size() == matrix.rows(), "Matrix size does not match!");
 
        const IndexType* row_ptr(matrix.row_ptr());
        Arch::UnitFilter::filter_offdiag_row_mat(matrix.val(), row_ptr, 1, _sv.indices(), _sv.used_elements());
      }
 
      void filter_offdiag_col_mat(SparseMatrixCSR<DT_, IT_> &) const
      {
        // nothing to do here
      }
 
      template<int block_width_>
      void filter_offdiag_row_mat(SparseMatrixBCSR<DT_, IT_, 1, block_width_> & matrix) const
      {
        if(_sv.used_elements() == Index(0))
          return;
 
        XASSERTM(_sv.size() == matrix.rows(), "Matrix size does not match!");
 
        const IndexType* row_ptr(matrix.row_ptr());
        Arch::UnitFilter::filter_offdiag_row_mat(matrix.template val<Perspective::pod>(), row_ptr, block_width_, _sv.indices(), _sv.used_elements());
      }
 
      template<int block_height_>
      void filter_offdiag_row_mat(SparseMatrixBCSR<DT_, IT_, block_height_, 1> &) const
      {
        // nothing to do here
      }
 
      void filter_weak_matrix_rows(SparseMatrixCSR<DT_, IT_> & matrix_a, const SparseMatrixCSR<DT_, IT_>& matrix_m) const
      {
        if(_sv.used_elements() == Index(0))
          return;
 
        XASSERTM(_sv.size() == matrix_a.rows(), "Matrix size does not match!");
        XASSERTM(_sv.size() == matrix_m.rows(), "Matrix size does not match!");
 
        const IndexType* row_ptr(matrix_a.row_ptr());
        const IndexType* col_idx(matrix_m.col_ind());
        XASSERTM(row_ptr == matrix_m.row_ptr(), "matrix A and M must share their layout");
        XASSERTM(col_idx == matrix_m.col_ind(), "matrix A and M must share their layout");
 
        DT_* val_a(matrix_a.val());
        const DT_* val_m(matrix_m.val());
 
        Arch::UnitFilter::filter_weak_matrix_rows(val_a, val_m, row_ptr, get_values(), get_indices(), _sv.used_elements());
      }
 
 
      void filter_rhs(DenseVector<DT_, IT_> & vector) const
      {
        if(_sv.used_elements() == Index(0))
          return;
        XASSERTM(_sv.size() == vector.size(), "Vector size does not match!");
        Arch::UnitFilter::filter_rhs(vector.elements(), _sv.elements(), _sv.indices(), _sv.used_elements());
      }
 
      void filter_sol(DenseVector<DT_, IT_> & vector) const
      {
        // same as rhs
        filter_rhs(vector);
      }
 
      void filter_def(DenseVector<DT_, IT_> & vector) const
      {
        if(_sv.used_elements() == Index(0))
          return;
        XASSERTM(_sv.size() == vector.size(), "Vector size does not match!");
        Arch::UnitFilter::filter_def(vector.elements(), _sv.indices(), _sv.used_elements());
      }
 
      void filter_cor(DenseVector<DT_, IT_> & vector) const
      {
        // same as def
        filter_def(vector);
      }
    }; // class UnitFilter<...>
  } // namespace LAFEM
} // namespace FEAT