DOFMatrix.hpp 9.88 KB
Newer Older
1
2
#pragma once

3
4
#include <list>
#include <map>
5
#include <memory>
6
7
#include <string>
#include <vector>
8
9
10
11
12
13

#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>

14
#include <dune/common/std/optional.hh>
15
16
#include <dune/common/reservedvector.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
17

18
19
#include <amdis/Output.hpp>
#include <amdis/common/ClonablePtr.hpp>
20
21
22

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

34
    /// The type of the finite element space / basis of the column
35
    using ColFeSpace = ColFeSpaceType;
36

37
    /// The matrix type of the underlying base matrix
38
    using BaseMatrix = mtl::compressed2D<ValueType>;
39

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

43
    /// The type of the elements of the DOFMatrix
44
    using value_type = ValueType;
45

46
    /// The underlying field-type (typically the same as \ref value_type)
47
    using field_type = typename BaseMatrix::value_type;
48

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

52
53
  public:
    /// Constructor. Constructs new BaseMatrix.
54
55
    DOFMatrix(RowFeSpace const& rowFeSpace,
              ColFeSpace const& colFeSpace,
56
              std::string name)
57
58
59
60
61
      : rowFeSpace(rowFeSpace)
      , colFeSpace(colFeSpace)
      , name(name)
      , matrix(ClonablePtr<BaseMatrix>::make())
    {}
62

63
    /// Constructor. Takes pointer of data-matrix.
64
65
66
    DOFMatrix(RowFeSpace const& rowFeSpace,
              ColFeSpace const& colFeSpace,
              std::string name,
67
68
69
70
71
72
              BaseMatrix& matrix_ref)
      : rowFeSpace(rowFeSpace)
      , colFeSpace(colFeSpace)
      , name(name)
      , matrix(matrix_ref)
    {}
73

74
75
76
77
78
    /// Return the row-basis \ref rowFeSpace of the matrix
    RowFeSpace const& getRowFeSpace() const
    {
      return rowFeSpace;
    }
79

80
81
82
83
84
    /// Return the col-basis \ref colFeSpace of the matrix
    ColFeSpace const& getColFeSpace() const
    {
      return colFeSpace;
    }
85

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

93
94
    /// Return a reference to the data-matrix \ref matrix
    BaseMatrix const& getMatrix() const
95
    {
96
      assert( !insertionMode );
97
98
      return *matrix;
    }
99

100
    /// Return the size of the row \ref feSpace
101
102
103
104
    size_type N() const
    {
      return rowFeSpace.size();
    }
105

106
    /// Return the size of the column \ref feSpace
107
108
109
110
    size_type M() const
    {
      return colFeSpace.size();
    }
111

112
    /// Return the \ref name of this matrix
113
114
115
116
    std::string getName() const
    {
      return name;
    }
117

118
119
120
121

    using FlatIndex = Dune::Functions::FlatMultiIndex<size_type>;
    using BlockIndex = Dune::ReservedVector<size_type,1>;

122
    /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
123
    /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
124
    /// be created by \ref init and must be closed by \ref finish after insertion.
125
126
127
    auto operator()(FlatIndex row, FlatIndex col)
    {
      size_type r = row[0], c = col[0];
128
      test_exit_dbg( insertionMode, "Inserter not initilized!");
129
130
131
132
133
134
      test_exit_dbg( r < N() && c < M() ,
          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
      return (*inserter)[r][c];
    }

    auto operator()(BlockIndex row, BlockIndex col)
135
    {
136
      size_type r = row[0], c = col[0];
137
      test_exit_dbg( insertionMode, "Inserter not initilized!");
138
139
      test_exit_dbg( r < N() && c < M() ,
          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
140
141
      return (*inserter)[r][c];
    }
142

143
    /// Create inserter. Assumes that no inserter is currently active on this matrix.
144
    void init(bool prepareForInsertion)
145
    {
146
      test_exit(!insertionMode && !inserter, "Matrix already in insertion mode!");
147

148
      calculateSlotSize();
149
150
151
      matrix->change_dim(rowFeSpace.dimension(), colFeSpace.dimension());
      test_exit(num_rows(*matrix) == rowFeSpace.dimension() && num_cols(*matrix) == colFeSpace.dimension(),
        "Wrong dimensions in matrix");
152
153
154
      if (prepareForInsertion) {
        set_to_zero(*matrix);
        inserter = new Inserter(*matrix, slot_size);
155
        insertionMode = true;
156
      }
157
158
    }

159
    /// Delete inserter -> finish insertion. Must be called in order to fill the
160
    /// final construction of the matrix.
161
162
163
164
    void finish()
    {
      delete inserter;
      inserter = NULL;
165
      insertionMode = false;
166
    }
167

168
    /// Return the number of nonzeros in the matrix
169
    std::size_t getNnz() const
170
171
172
    {
      return matrix->nnz();
    }
173

174
175
    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
    /// a one on the diagonal of the matrix.
176
    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
177
    {
178
      std::list<Triplet<value_type>> columns;
179

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

185
      // iterate over the matrix
186
#if 0
187
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
188
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
189
          if (dirichletNodes[row(i)]) {
190
            // set identity row
191
192
193
            value(i, (row(i) == col(i) ? value_type(1) : value_type(0)) );
          }
          else if (dirichletNodes[col(i)]) {
194
            columns.push_back({row(i), col(i), value(i)});
195
196
197
198
            value(i, value_type(0));
          }
        }
      }
199
200
201
202
203
204
205
206
207
208
#endif

      for (auto r : mtl::rows_of(*matrix)) {   // rows or columns
        if (dirichletNodes[r.value()]) {
          for (auto i : mtl::nz_of(r)) {          // non-zeros within
            // set identity row
            value(i, (row(i) == col(i) ? value_type(1) : value_type(0)) );
          }
        }
      }
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257

      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)) );
          }
        }
      }
258

259
260
261
262
263
264
265
266
      // 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);
267
268
269
        }
      }
    }
270
#endif
271

272
273
274
  protected:
    // Estimates the slot-size used in the inserter to reserve enough space per row.
    void calculateSlotSize()
275
276
277
278
279
280
281
282
    {
      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;
    }
283

284
  private:
285
    /// The minimal number of nonzeros per row
286
    static constexpr int MIN_NNZ_PER_ROW = 50;
287

288
289
    /// The finite element space / basis of the row
    RowFeSpace const& rowFeSpace;
290

291
292
    /// The finite element space / basis of the column
    ColFeSpace const& colFeSpace;
293

294
295
    /// The name of the DOFMatrix
    std::string name;
296

297
298
    /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
    ClonablePtr<BaseMatrix> matrix;
299

300
301
    /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
    Inserter* inserter = NULL;
302

303
    /// A flag that indicates whether the matrix is in insertion mode
304
    bool insertionMode = false;
305

306
    /// The size of the slots used to insert values per row
307
    int slot_size = MIN_NNZ_PER_ROW;
308

309
310
311
312
313
    // friend class declarations
    template <class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS