DOFMatrix.hpp 4.66 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
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
114
115
116
117
    /// 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 colLocalView = colFeSpace.localView();

      for (const auto& element : elements(meshView)) {
        rowLocalView.bind(element);
        colLocalView.bind(element);
118

Praetorius, Simon's avatar
Praetorius, Simon committed
119
        for (std::size_t i = 0; i < rowLocalView.size(); ++i) {
120
          // The global index of the i-th vertex of the element
Praetorius, Simon's avatar
Praetorius, Simon committed
121
122
          auto row = rowLocalView.index(i);
          for (std::size_t j = 0; j < colLocalView.size(); ++j) {
123
            // The global index of the j-th vertex of the element
Praetorius, Simon's avatar
Praetorius, Simon committed
124
            auto col = colLocalView.index(j);
125
126
127
128
129
            occupationPattern.add(row, col);
          }
        }
      }
      occupationPattern.exportIdx(*matrix);
130

131
132
133
134
135
136
137
      initialized = true;
    }

    void finish()
    {
      initialized = false;
    }
138

139
    std::size_t nnz()
140
141
142
    {
      return matrix->nonzeroes();
    }
143

144
    auto clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
145
    {
146
      // loop over the matrix rows
147
      for (std::size_t i = 0; i < matrix->N(); ++i) {
148
149
150
151
152
153
154
155
        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);
        }
      }
156
157
158

      std::vector<std::map<size_t,value_type>> columns; // symmetric dbc not yet implemented
      return columns;
159
    }
160

161
162
163
164
  private:
    RowFeSpace const&	rowFeSpace;
    ColFeSpace const&	colFeSpace;
    std::string	        name;
165

166
167
    ClonablePtr<BaseMatrix> matrix;
    Dune::MatrixIndexSet    occupationPattern;
168

169
    bool initialized = false;
170

171
172
173
174
175
    // friend class declarations
    template <class, class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS