DOFMatrix.hpp 4.81 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
#include <dune/amdis/Output.hpp>
10
#include <dune/amdis/common/ClonablePtr.hpp>
11
12
13

namespace AMDiS
{
14
  template <class RowFeSpaceType, class ColFeSpaceType,
15
16
17
18
19
20
21
            class ValueType = Dune::FieldMatrix<double,1,1>>
  class DOFMatrix
  {
  public:
    using RowFeSpace = RowFeSpaceType;
    using ColFeSpace = ColFeSpaceType;
    using BaseMatrix = Dune::BCRSMatrix<ValueType>;
22

23
24
    using size_type  = typename RowFeSpaceType::size_type;
    using field_type = typename ValueType::field_type;
25

26
    using value_type = ValueType;
27

28
29
30
31
32
33
34
    /// Constructor. Constructs new BaseVector.
    DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name)
      : rowFeSpace(rowFeSpace)
      , colFeSpace(colFeSpace)
      , name(name)
      , matrix(ClonablePtr<BaseMatrix>::make())
    {}
35

36
    /// Constructor. Takes pointer of data-vector.
37
    DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name,
38
39
40
41
42
43
              BaseMatrix& matrix_ref)
      : rowFeSpace(rowFeSpace)
      , colFeSpace(colFeSpace)
      , name(name)
      , matrix(matrix_ref)
    {}
44

45
46
47
48
49
    /// Return the row-basis \ref rowFeSpace of the matrix
    RowFeSpace const& getRowFeSpace() const
    {
      return rowFeSpace;
    }
50

51
52
53
54
55
    /// Return the col-basis \ref colFeSpace of the matrix
    ColFeSpace const& getColFeSpace() const
    {
      return colFeSpace;
    }
56

57
58
59
60
61
    /// Return the data-vector \ref vector
    BaseMatrix const& getMatrix() const
    {
      return *matrix;
    }
62

63
64
65
66
67
    /// Return the data-vector \ref vector
    BaseMatrix& getMatrix()
    {
      return *matrix;
    }
68

69
70
71
72
73
    /// Return the size of the \ref feSpace
    size_type N() const
    {
      return rowFeSpace.size();
    }
74

75
76
77
78
    size_type M() const
    {
      return colFeSpace.size();
    }
79

80
81
82
83
84
    /// Return the \ref name of this vector
    std::string getName() const
    {
      return name;
    }
85
86


87
88
89
    /// Access the entry \p i of the \ref vector with read-access.
    value_type const& operator()(size_type r, size_type c) const
    {
90
91
92
      test_exit_dbg( initialized, "Occupation pattern not initialized!");
      test_exit_dbg( r < N() && c < M() ,
          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
93
94
      return (*matrix)[r][c];
    }
95

96
97
98
    /// Access the entry \p i of the \ref vector with write-access.
    value_type& operator()(size_type r, size_type c)
    {
99
100
101
      test_exit_dbg( initialized, "Occupation pattern not initialized!");
      test_exit_dbg( r < N() && c < M() ,
          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
102
103
      return (*matrix)[r][c];
    }
104

105
106
107
108
109
110
111
112
113
    /// create occupation pattern and apply it to the matrix
    void init()
    {
      occupationPattern.resize(rowFeSpace.size(), colFeSpace.size());
      auto meshView = rowFeSpace.gridView();

      // A loop over all elements of the grid
      auto rowLocalView = rowFeSpace.localView();
      auto rowIndexSet  = rowFeSpace.localIndexSet();
114

115
116
117
118
119
120
      auto colLocalView = colFeSpace.localView();
      auto colIndexSet  = colFeSpace.localIndexSet();

      for (const auto& element : elements(meshView)) {
        rowLocalView.bind(element);
        rowIndexSet.bind(rowLocalView); // TODO: can this be moved out of the loop?
121

122
123
        colLocalView.bind(element);
        colIndexSet.bind(colLocalView); // TODO: can this be moved out of the loop?
124

125
126
127
128
129
130
131
132
133
134
135
        for (size_t i = 0; i < rowIndexSet.size(); ++i) {
          // The global index of the i-th vertex of the element
          auto row = rowIndexSet.index(i);
          for (size_t j = 0; j < colIndexSet.size(); ++j) {
            // The global index of the j-th vertex of the element
            auto col = colIndexSet.index(j);
            occupationPattern.add(row, col);
          }
        }
      }
      occupationPattern.exportIdx(*matrix);
136

137
138
139
140
141
142
143
      initialized = true;
    }

    void finish()
    {
      initialized = false;
    }
144

145
146
147
148
    size_t getNnz()
    {
      return matrix->nonzeroes();
    }
149

150
    void clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
151
    {
152
153
154
155
156
157
158
159
160
161
162
      // loop over the matrix rows
      for (size_t i = 0; i < matrix->N(); ++i) {
        if (dirichletNodes[i]) {
            auto cIt = (*matrix)[i].begin();
            auto cEndIt = (*matrix)[i].end();
            // loop over nonzero matrix entries in current row
            for (; cIt != cEndIt; ++cIt)
                *cIt = (apply && i == cIt.index() ? 1 : 0);
        }
      }
    }
163

164
165
166
167
  private:
    RowFeSpace const&	rowFeSpace;
    ColFeSpace const&	colFeSpace;
    std::string	        name;
168

169
170
    ClonablePtr<BaseMatrix> matrix;
    Dune::MatrixIndexSet    occupationPattern;
171

172
    bool initialized = false;
173

174
175
176
177
178
    // friend class declarations
    template <class, class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS