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

#include <string>
#include <memory>

#include <boost/numeric/mtl/matrix/compressed2D.hpp>
#include <boost/numeric/mtl/matrix/inserter.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_wrapper.hpp>

11
12
#include <dune/common/std/optional.hh>

13
#include <dune/amdis/Output.hpp>
14
#include <dune/amdis/common/ClonablePtr.hpp>
15
16
17

namespace AMDiS
{
18
  /// \brief The basic container that stores a base matrix and a corresponding
19
  /// row/column feSpace.
20
  template <class RowFeSpaceType, class ColFeSpaceType,
21
22
23
24
            class ValueType = double>
  class DOFMatrix
  {
  public:
25
    /// The type of the finite element space / basis of the row
26
    using RowFeSpace = RowFeSpaceType;
27

28
    /// The type of the finite element space / basis of the column
29
    using ColFeSpace = ColFeSpaceType;
30

31
    /// The matrix type of the underlying base matrix
32
    using BaseMatrix = mtl::compressed2D<ValueType>;
33

34
    /// The index/size - type
35
    using size_type  = typename BaseMatrix::size_type;
36

37
    /// The type of the elements of the DOFMatrix
38
    using value_type = ValueType;
39

40
    /// The underlying field-type (typically the same as \ref value_type)
41
    using field_type = typename BaseMatrix::value_type;
42

43
    /// The type of the inserter used when filling the matrix. NOTE: need not be public
44
    using Inserter = mtl::mat::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;
45

46
47
  public:
    /// Constructor. Constructs new BaseMatrix.
48
49
    DOFMatrix(RowFeSpace const& rowFeSpace,
              ColFeSpace const& colFeSpace,
50
              std::string name)
51
52
53
54
55
      : rowFeSpace(rowFeSpace)
      , colFeSpace(colFeSpace)
      , name(name)
      , matrix(ClonablePtr<BaseMatrix>::make())
    {}
56

57
    /// Constructor. Takes pointer of data-matrix.
58
59
60
    DOFMatrix(RowFeSpace const& rowFeSpace,
              ColFeSpace const& colFeSpace,
              std::string name,
61
62
63
64
65
66
              BaseMatrix& matrix_ref)
      : rowFeSpace(rowFeSpace)
      , colFeSpace(colFeSpace)
      , name(name)
      , matrix(matrix_ref)
    {}
67

68
69
70
71
72
    /// Return the row-basis \ref rowFeSpace of the matrix
    RowFeSpace const& getRowFeSpace() const
    {
      return rowFeSpace;
    }
73

74
75
76
77
78
    /// Return the col-basis \ref colFeSpace of the matrix
    ColFeSpace const& getColFeSpace() const
    {
      return colFeSpace;
    }
79

80
81
    /// Return a reference to the data-matrix \ref matrix
    BaseMatrix& getMatrix()
82
83
84
    {
      return *matrix;
    }
85

86
87
    /// Return a reference to the data-matrix \ref matrix
    BaseMatrix const& getMatrix() const
88
89
90
    {
      return *matrix;
    }
91

92
    /// Return the size of the row \ref feSpace
93
94
95
96
    size_type N() const
    {
      return rowFeSpace.size();
    }
97

98
    /// Return the size of the column \ref feSpace
99
100
101
102
    size_type M() const
    {
      return colFeSpace.size();
    }
103

104
    /// Return the \ref name of this matrix
105
106
107
108
    std::string getName() const
    {
      return name;
    }
109
110

    /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
111
112
113
    /// position (\p r, \p c) in the matrix. Need an initialized inserter, that can
    /// be created by \ref init and must be closed by \ref finish after insertion.
    auto operator()(size_type r, size_type c)
114
    {
115
116
117
      test_exit_dbg( initialized, "Inserter not initialized!");
      test_exit_dbg( r < N() && c < M() ,
          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
118
119
      return (*inserter)[r][c];
    }
120

121
    /// Create inserter. Assumes that no inserter is currently active on this matrix.
122
    void init(bool prepareForInsertion)
123
    {
124
125
      test_exit(!initialized && !inserter, "Matrix already in insertion mode!");

126
      calculateSlotSize();
127
      matrix->change_dim(rowFeSpace.size(), colFeSpace.size());
128
129
130
131
132
      if (prepareForInsertion) {
        set_to_zero(*matrix);
        inserter = new Inserter(*matrix, slot_size);
        initialized = true;
      }
133
134
    }

135
    /// Delete inserter -> finish insertion. Must be called in order to fill the
136
    /// final construction of the matrix.
137
138
139
140
141
142
    void finish()
    {
      delete inserter;
      inserter = NULL;
      initialized = false;
    }
143

144
    /// Return the number of nonzeros in the matrix
145
    std::size_t getNnz() const
146
147
148
    {
      return matrix->nnz();
    }
149

150
151
    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
    /// a one on the diagonal of the matrix.
152
    auto applyDirichletBC(std::vector<char> const& dirichletNodes, bool apply)
153
    {
154
155
      std::vector<std::map<size_t,value_type>> columns(num_rows(*matrix));

156
      // Define the property maps
157
      auto row   = mtl::mat::row_map(*matrix);
158
      auto col   = mtl::mat::col_map(*matrix);
159
      auto value = mtl::mat::value_map(*matrix);
160

161
      // iterate over the matrix
162
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
163
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
          if (dirichletNodes[row(i)]) {
            // set identity row
            value(i, (apply && row(i) == col(i) ? value_type(1) : value_type(0)) );
          } else if (apply && dirichletNodes[col(i)]) {
            columns[row(i)][col(i)] = value(i);
            value(i, value_type(0));
          }
        }
      }

      return columns;
    }

#if 0
    void applyPeriodicBC(std::vector<char> const& periodicNodes,
                         std::map<size_t, size_t> const& association, bool applyRow, bool applyCol)
    {
      bool apply = applyRow && applyCol;

      // Define the property maps
      auto row   = mtl::mat::row_map(*matrix);
      auto col   = mtl::mat::col_map(*matrix);
      auto value = mtl::mat::value_map(*matrix);

      std::vector<std::map<size_t, std::list<value_type>>> row_values(num_cols(*matrix));
      std::vector<std::map<size_t, std::list<value_type>>> col_values(num_rows(*matrix));
      std::vector<char> dualNodes(periodicNodes.size(), 0);

      // iterate over the matrix to collect the values and erase the source-row/col
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
          if (applyRow && periodicNodes[row(i)]) {
            row_values[col(i)][association[row(i)]].push_back(value(i));
            value(i, value_type(0));
            dualNodes[association[row(i)]] = 1;
          } else if (applyCol && periodicNodes[col(i)]) {
            col_values[row(i)][association[col(i)]].push_back(value(i));
            value(i, value_type(0));
          }
        }
      }

      // TODO: use inserter for assignment of values!!!
      // iterate over the matrix to assign the values
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
          if (applyRow && dualNodes[row(i)]) {
            value(i, std::accumulate(row_values[col(i)][row(i)].begin(),
                                     row_values[col(i)][row(i)].end(),
                                     value(i)) );
          }
          if (applyCol && dualNodes[col(i)]) {
            value(i, std::accumulate(col_values[row(i)][col(i)].begin(),
                                     col_values[row(i)][col(i)].end(),
                                     value(i)) );
          }
        }
      }
222

223
224
225
226
227
228
229
230
      // assign values 1, -1 to diagonal and associated entries
      if (apply) {
        Inserter ins(*matrix);
        for (auto const& a : association) {
          ins[a.first][a.first]   << value_type( 1);
          ins[a.second][a.second] << value_type( 1);
          ins[a.first][a.second]  << value_type(-1);
          ins[a.second][a.first]  << value_type(-1);
231
232
233
        }
      }
    }
234
#endif
235

236
237
238
  protected:
    // Estimates the slot-size used in the inserter to reserve enough space per row.
    void calculateSlotSize()
239
240
241
242
243
244
245
246
    {
      slot_size = 0;

      if (num_rows(*matrix) != 0)
        slot_size = int(double(matrix->nnz()) / num_rows(*matrix) * 1.2);
      if (slot_size < MIN_NNZ_PER_ROW)
        slot_size = MIN_NNZ_PER_ROW;
    }
247

248
  private:
249
250
    /// The minimal number of nonzeros per row
    static constexpr int MIN_NNZ_PER_ROW = 20;
251

252
253
    /// The finite element space / basis of the row
    RowFeSpace const& rowFeSpace;
254

255
256
    /// The finite element space / basis of the column
    ColFeSpace const& colFeSpace;
257

258
259
    /// The name of the DOFMatrix
    std::string name;
260

261
262
    /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
    ClonablePtr<BaseMatrix> matrix;
263

264
265
    /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
    Inserter* inserter = NULL;
266

267
    /// A flag that indicates whether the matrix is in insertion mode
268
    bool initialized = false;
269

270
271
    /// The size of the slots used to insert values per row
    int slot_size = 20;
272

273
274
275
276
277
    // friend class declarations
    template <class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS