DOFMatrix.hpp 4.99 KB
Newer Older
1
2
3
4
5
6
7
8
#pragma once

#include <string>
#include <memory>

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

9
10
#include <amdis/Output.hpp>
#include <amdis/common/ClonablePtr.hpp>
11
12
#include <amdis/linear_algebra/DOFMatrixBase.hpp>
#include <amdis/utility/MultiIndex.hpp>
13
14
15

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

  template <class RowBasisType, class ColBasisType,
            class ValueType = double>
30
31
32
  class DOFMatrix
  {
  public:
33
34
35
36
37
    /// The type of the finite element space / basis of the row
    using RowBasis = RowBasisType;

    /// The type of the finite element space / basis of the column
    using ColBasis = ColBasisType;
38

39
40
    /// The type of the elements of the DOFMatrix
    using value_type = typename BlockMatrixType<ValueType>::type;
41

42
43
    /// The matrix type of the underlying base matrix
    using BaseMatrix = Dune::BCRSMatrix<value_type>;
44

45
46
47
48
    /// The index/size - type
    using size_type = typename BaseMatrix::size_type;

  public:
49
    /// Constructor. Constructs new BaseVector.
50
51
52
53
    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis)
      : rowBasis_(&rowBasis)
      , colBasis_(&colBasis)
      , matrix_(ClonablePtr<BaseMatrix>::make())
54
    {}
55

56
57
58
59
60
    /// Constructor. Takes reference to a base matrix.
    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis, BaseMatrix& matrix)
      : rowBasis_(&rowBasis)
      , colBasis_(&colBasis)
      , matrix_(matrix)
61
    {}
62

63
64
    /// Return the row-basis \ref rowBasis of the matrix
    RowBasis const& rowBasis() const
65
    {
66
      return *rowBasis_;
67
    }
68

69
70
    /// Return the col-basis \ref colBasis of the matrix
    ColBasis const& colBasis() const
71
    {
72
      return *colBasis_;
73
    }
74

75
    /// Return the data-vector \ref vector
76
    BaseMatrix const& matrix() const
77
    {
78
      return *matrix_;
79
    }
80

81
    /// Return the data-vector \ref vector
82
    BaseMatrix& matrix()
83
    {
84
      return *matrix_;
85
    }
86

87
    /// Return the size of the \ref feSpace
88
    size_type rows() const
89
    {
90
      return rowBasis_->size();
91
    }
92

93
    size_type cols() const
94
    {
95
      return colBasis_->size();
96
    }
97
98


99
    /// Access the entry \p i of the \ref vector with read-access.
100
101
    template <class Index>
    value_type const& operator()(Index row, Index col) const
102
    {
103
104
105
106
107
      size_type r = flatMultiIndex(row), c = flatMultiIndex(col);
      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
      test_exit_dbg( r < rows() && c < cols() ,
          "Indices out of range [0,", rows(), ")x[0,", cols(), ")" );
      return (*matrix_)[r][c];
108
    }
109

110
    /// Access the entry \p i of the \ref vector with write-access.
111
112
    template <class Index>
    value_type& operator()(Index row, Index col)
113
    {
114
115
116
117
118
      size_type r = flatMultiIndex(row), c = flatMultiIndex(col);
      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
      test_exit_dbg( r < rows() && c < cols() ,
          "Indices out of range [0,", rows(), ")x[0,", cols(), ")" );
      return (*matrix_)[r][c];
119
    }
120

121
    /// create occupation pattern and apply it to the matrix
122
    void init(bool prepareForInsertion)
123
    {
124
      occupationPattern_.resize(rowBasis_->size(), colBasis_->size());
125
126

      // A loop over all elements of the grid
127
128
      auto rowLocalView = rowBasis_->localView();
      auto colLocalView = colBasis_->localView();
129

130
      for (const auto& element : elements(rowBasis_->gridView())) {
131
132
        rowLocalView.bind(element);
        colLocalView.bind(element);
133

Praetorius, Simon's avatar
Praetorius, Simon committed
134
        for (std::size_t i = 0; i < rowLocalView.size(); ++i) {
135
          // The global index of the i-th vertex of the element
Praetorius, Simon's avatar
Praetorius, Simon committed
136
137
          auto row = rowLocalView.index(i);
          for (std::size_t j = 0; j < colLocalView.size(); ++j) {
138
            // The global index of the j-th vertex of the element
Praetorius, Simon's avatar
Praetorius, Simon committed
139
            auto col = colLocalView.index(j);
140
            occupationPattern_.add(row, col);
141
142
143
          }
        }
      }
144
      occupationPattern_.exportIdx(*matrix_);
145

146
      initialized_ = true;
147
148
149
150
    }

    void finish()
    {
151
      initialized_ = false;
152
    }
153

154
    std::size_t nnz() const
155
    {
156
      return matrix_->nonzeroes();
157
    }
158

159
    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
160
    {
161
      // loop over the matrix rows
162
      for (std::size_t i = 0; i < matrix_->N(); ++i) {
163
        if (dirichletNodes[i]) {
164
165
            auto cIt = (*matrix_)[i].begin();
            auto cEndIt = (*matrix_)[i].end();
166
167
            // loop over nonzero matrix entries in current row
            for (; cIt != cEndIt; ++cIt)
168
                *cIt = (i == cIt.index() ? 1 : 0);
169
170
        }
      }
171

172
      std::list<Triplet<value_type>> columns; // symmetric dbc not yet implemented
173
      return columns;
174
    }
175

176
  private:
177
178
    RowBasis const*	rowBasis_;
    ColBasis const*	colBasis_;
179

180
181
    ClonablePtr<BaseMatrix> matrix_;
    Dune::MatrixIndexSet occupationPattern_;
182

183
    bool initialized_ = false;
184
185
186
  };

} // end namespace AMDiS