DOFMatrix.hpp 3.48 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
13
#include <amdis/linear_algebra/DOFMatrixBase.hpp>
#include <amdis/utility/MultiIndex.hpp>
14
15
16

namespace AMDiS
{
17
18
19
20
21
22
23
24
25
26
27
28
  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;
  };

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

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

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

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

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

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


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

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

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

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

93
      initialized_ = true;
94
95
96
97
    }

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

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

106
    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
107
    {
108
      // loop over the matrix rows
109
      for (std::size_t i = 0; i < matrix_.N(); ++i) {
110
        if (dirichletNodes[i]) {
111
112
            auto cIt = matrix_[i].begin();
            auto cEndIt = matrix_[i].end();
113
114
            // loop over nonzero matrix entries in current row
            for (; cIt != cEndIt; ++cIt)
115
                *cIt = (i == cIt.index() ? 1 : 0);
116
117
        }
      }
118

Praetorius, Simon's avatar
Praetorius, Simon committed
119
      return std::list<Triplet<value_type>>{};
120
    }
121

122
  private:
123
    BaseMatrix matrix_;
124

125
    bool initialized_ = false;
126
127
  };

128
129
130
  template <class RowBasisType, class ColBasisType, class ValueType = double>
  using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, IstlMatrix<ValueType>>;

131
} // end namespace AMDiS