DOFMatrix.hpp 2.95 KB
Newer Older
1
2
#pragma once

Praetorius, Simon's avatar
Praetorius, Simon committed
3
#include <list>
4
5
6
7
8
9
#include <string>
#include <memory>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>

10
#include <amdis/Output.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
11
#include <amdis/linear_algebra/Common.hpp>
12
#include <amdis/linear_algebra/DOFMatrixBase.hpp>
13
14
15

namespace AMDiS
{
16
17
18
19
20
21
22
23
24
25
26
27
  template <class T, class = void>
  struct BlockMatrixType
  {
    using type = Dune::FieldMatrix<T,1,1>;
  };

  template <class T>
  struct BlockMatrixType<T, typename T::field_type>
  {
    using type = T;
  };

28
29
  template <class ValueType>
  class IstlMatrix
30
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
31
  public:
32
33
    /// The type of the elements of the DOFMatrix
    using value_type = typename BlockMatrixType<ValueType>::type;
34

35
36
    /// The matrix type of the underlying base matrix
    using BaseMatrix = Dune::BCRSMatrix<value_type>;
37

38
39
40
    /// The index/size - type
    using size_type = typename BaseMatrix::size_type;

41
  public:
42
    /// Constructor. Constructs new BaseVector.
43
    IstlMatrix() = default;
44

45
    /// Return the data-vector \ref vector
46
    BaseMatrix const& matrix() const
47
    {
48
      return matrix_;
49
    }
50

51
    /// Return the data-vector \ref vector
52
    BaseMatrix& matrix()
53
    {
54
      return matrix_;
55
    }
56
57


58
59
    /// Insert a single value into the matrix (add to existing value)
    void insert(size_type r, size_type c, value_type const& value)
60
    {
61
      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
62
      test_exit_dbg( r < matrix_.N() && c < matrix_.M() ,
63
          "Indices out of range [0,{})x[0,{})", matrix_.N(), matrix_.M() );
64
      matrix_[r][c] += value;
65
    }
66

67
    /// create occupation pattern and apply it to the matrix
68
69
70
    template <class RowBasis, class ColBasis>
    void init(RowBasis const& rowBasis, ColBasis const& colBasis,
              bool prepareForInsertion)
71
    {
72
      auto occupationPattern = Dune::MatrixIndexSet{rowBasis.dimension(), colBasis.dimension()};
73

74
75
76
      auto rowLocalView = rowBasis.localView();
      auto colLocalView = colBasis.localView();
      for (const auto& element : elements(rowBasis.gridView())) {
77
78
        rowLocalView.bind(element);
        colLocalView.bind(element);
79

Praetorius, Simon's avatar
Praetorius, Simon committed
80
        for (std::size_t i = 0; i < rowLocalView.size(); ++i) {
81
          // The global index of the i-th vertex of the element
Praetorius, Simon's avatar
Praetorius, Simon committed
82
83
          auto row = rowLocalView.index(i);
          for (std::size_t j = 0; j < colLocalView.size(); ++j) {
84
            // The global index of the j-th vertex of the element
Praetorius, Simon's avatar
Praetorius, Simon committed
85
            auto col = colLocalView.index(j);
86
            occupationPattern.add(row, col);
87
88
89
          }
        }
      }
90
      occupationPattern.exportIdx(matrix_);
91

92
      initialized_ = true;
93
94
95
96
    }

    void finish()
    {
97
      initialized_ = false;
98
    }
99

100
    std::size_t nnz() const
101
    {
102
      return matrix_.nonzeroes();
103
    }
104

105
  private:
106
    BaseMatrix matrix_;
107

108
    bool initialized_ = false;
109
110
  };

111
112
113
  template <class RowBasisType, class ColBasisType, class ValueType = double>
  using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, IstlMatrix<ValueType>>;

114
} // end namespace AMDiS