DOFMatrix.hpp 9.71 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
      matrix->change_dim(rowFeSpace.size(), colFeSpace.size());
150
151
152
      if (prepareForInsertion) {
        set_to_zero(*matrix);
        inserter = new Inserter(*matrix, slot_size);
153
        insertionMode = true;
154
      }
155
156
    }

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

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

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

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

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

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

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

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

282
  private:
283
284
    /// The minimal number of nonzeros per row
    static constexpr int MIN_NNZ_PER_ROW = 20;
285

286
287
    /// The finite element space / basis of the row
    RowFeSpace const& rowFeSpace;
288

289
290
    /// The finite element space / basis of the column
    ColFeSpace const& colFeSpace;
291

292
293
    /// The name of the DOFMatrix
    std::string name;
294

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

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

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

304
305
    /// The size of the slots used to insert values per row
    int slot_size = 20;
306

307
308
309
310
311
    // friend class declarations
    template <class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS