MatrixFacade.hpp 3.44 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#pragma once

#include <type_traits>

#include <dune/common/shared_ptr.hh>
#include <dune/functions/functionspacebases/concepts.hh>

#include <amdis/common/Concepts.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/functions/NodeIndices.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
#include <amdis/typetree/MultiIndex.hpp>

namespace AMDiS
{
  /**
   * Basis implementation of DOFMatrix, i.e. a sparse matrix storing all the
   * assembled Operators indexed with DOF indices. The matrix data is associated
   * to a row and column global basis.
   *
   * \tparam T          The coefficient type of the matrix
   * \tparam Pattern    The type of the sparsity pattern
   * \tparam MatrixImpl A linear-algebra backend for the matrix storage
   **/
25
  template <class T, template <class> class MatrixImpl>
Praetorius, Simon's avatar
Praetorius, Simon committed
26
27
28
29
30
31
  class MatrixFacade
  {
    using Self = MatrixFacade;
    using Impl = MatrixImpl<T>;

  public:
32
    /// Constructor. Forwards the bases to the implementation class and
Praetorius, Simon's avatar
Praetorius, Simon committed
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
    /// constructs a matrix sparsity pattern.
    template <class RowBasis, class ColBasis>
    MatrixFacade(RowBasis const& rowBasis, ColBasis const& colBasis)
      : impl_(rowBasis, colBasis)
    {}

    /// Return the underlying matrix backend
    Impl const& impl() const { return impl_; }
    Impl&       impl()       { return impl_; }

    /// \brief Initialize the matrix for insertion, i.e. allocate the non-zero pattern
    /**
     * With the optional parameter \p symmetry some additional information about the
     * structure of the values or the sparsity pattern can be provided.
     * See \ref SymmetryStructure.
     **/
49
50
    template <class SparsityPattern>
    void init(SparsityPattern const& pattern, SymmetryStructure symmetry = SymmetryStructure::unknown)
Praetorius, Simon's avatar
Praetorius, Simon committed
51
    {
52
      impl_.init(pattern, symmetry);
Praetorius, Simon's avatar
Praetorius, Simon committed
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
    }

    /// Finish the matrix insertion, e.g. cleanup or final insertion
    void finish()
    {
      impl_.finish();
    }

    /// Insert a single value into the matrix (add to existing value)
    template <class RowIndex, class ColIndex,
      REQUIRES(Concepts::MultiIndex<RowIndex>),
      REQUIRES(Concepts::MultiIndex<ColIndex>)>
    void insert(RowIndex const& row, ColIndex const& col, typename Impl::value_type const& value)
    {
      impl_.insert(row, col, value);
    }

    /// Insert a block of values into the sparse matrix (add to existing values)
    /// The global matrix indices are determined by the corresponding localviews.
    template <class RowLocalView, class ColLocalView, class LocalMatrix,
      REQUIRES(Dune::models<Dune::Functions::Concept::LocalView<typename RowLocalView::GlobalBasis>, RowLocalView>()),
      REQUIRES(Dune::models<Dune::Functions::Concept::LocalView<typename ColLocalView::GlobalBasis>, ColLocalView>())>
    void scatter(RowLocalView const& r, ColLocalView const& c, LocalMatrix const& localMatrix)
    {
      assert(r.size() * c.size() == localMatrix.size());
      assert(r.size() == localMatrix.rows());
      assert(c.size() == localMatrix.cols());

      const bool optimized = std::is_same<RowLocalView,ColLocalView>::value
        && r.tree().treeIndex() == c.tree().treeIndex();

      if (optimized)
        impl_.scatter(nodeIndices(r), localMatrix);
      else
        impl_.scatter(nodeIndices(r), nodeIndices(c), localMatrix);
    }

    /// Number of nonzeros in the matrix
    std::size_t nnz() const
    {
      return impl_.nnz();
    }

  protected:
    /// The matrix backend
    Impl impl_;
  };

} // end namespace AMDiS