DOFMatrix.hpp 9.5 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 RowBasisType,
            class ColBasisType,
27
28
29
30
            class ValueType = double>
  class DOFMatrix
  {
  public:
31
    /// The type of the finite element space / basis of the row
32
    using RowBasis = RowBasisType;
33

34
    /// The type of the finite element space / basis of the column
35
    using ColBasis = ColBasisType;
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
56
57
58
    DOFMatrix(RowBasis const& rowBasis,
              ColBasis const& colBasis)
      : rowBasis_(&rowBasis)
      , colBasis_(&colBasis)
      , matrix_(ClonablePtr<BaseMatrix>::make())
59
    {}
60

61
    /// Constructor. Takes pointer of data-matrix.
62
63
    DOFMatrix(RowBasis const& rowBasis,
              ColBasis const& colBasis,
64
              BaseMatrix& matrix_ref)
65
66
67
      : rowBasis_(&rowBasis)
      , colBasis_(&colBasis)
      , matrix_(matrix_ref)
68
    {}
69

70
71
    /// Return the row-basis \ref rowBasis_ of the matrix
    RowBasis const& rowBasis() const
72
    {
73
      return *rowBasis_;
74
    }
75

76
77
    /// Return the col-basis \ref colBasis_ of the matrix
    ColBasis const& colBasis() const
78
    {
79
      return *colBasis_;
80
    }
81

82
    /// Return a reference to the data-matrix \ref matrix
83
    BaseMatrix& matrix()
84
    {
85
86
      assert( !insertionMode_ );
      return *matrix_;
87
    }
88

89
    /// Return a reference to the data-matrix \ref matrix
90
    BaseMatrix const& matrix() const
91
    {
92
93
      assert( !insertionMode_ );
      return *matrix_;
94
    }
95

96
    /// Return the size of the row \ref feSpace
97
98
    size_type N() const
    {
99
      return rowBasis_->size();
100
    }
101

102
    /// Return the size of the column \ref feSpace
103
104
    size_type M() const
    {
105
      return colBasis_->size();
106
    }
107

108
109
110
111

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

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

    auto operator()(BlockIndex row, BlockIndex col)
125
    {
126
      size_type r = row[0], c = col[0];
127
      test_exit_dbg( insertionMode_, "Inserter not initilized!");
128
129
      test_exit_dbg( r < N() && c < M() ,
          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
130
      return (*inserter_)[r][c];
131
    }
132

133
    /// Create inserter. Assumes that no inserter is currently active on this matrix.
134
    void init(bool prepareForInsertion)
135
    {
136
      test_exit(!insertionMode_ && !inserter_, "Matrix already in insertion mode!");
137

138
      calculateSlotSize();
139
140
      matrix_->change_dim(rowBasis_->dimension(), colBasis_->dimension());
      test_exit(num_rows(*matrix_) == rowBasis_->dimension() && num_cols(*matrix_) == colBasis_->dimension(),
141
        "Wrong dimensions in matrix");
142
      if (prepareForInsertion) {
143
144
145
        set_to_zero(*matrix_);
        inserter_ = new Inserter(*matrix_, slotSize_);
        insertionMode_ = true;
146
      }
147
148
    }

149
    /// Delete inserter -> finish insertion. Must be called in order to fill the
150
    /// final construction of the matrix.
151
152
    void finish()
    {
153
154
155
      delete inserter_;
      inserter_ = NULL;
      insertionMode_ = false;
156
    }
157

158
    /// Return the number of nonzeros in the matrix
159
    std::size_t nnz() const
160
    {
161
      return matrix_->nnz();
162
    }
163

164
165
    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
    /// a one on the diagonal of the matrix.
166
    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
167
    {
168
      std::list<Triplet<value_type>> columns;
169

170
      // Define the property maps
171
172
173
      auto row   = mtl::mat::row_map(*matrix_);
      auto col   = mtl::mat::col_map(*matrix_);
      auto value = mtl::mat::value_map(*matrix_);
174

175
      // iterate over the matrix
176
#if 0
177
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
178
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
179
          if (dirichletNodes[row(i)]) {
180
            // set identity row
181
182
183
            value(i, (row(i) == col(i) ? value_type(1) : value_type(0)) );
          }
          else if (dirichletNodes[col(i)]) {
184
            columns.push_back({row(i), col(i), value(i)});
185
186
187
188
            value(i, value_type(0));
          }
        }
      }
189
190
#endif

191
      for (auto r : mtl::rows_of(*matrix_)) {   // rows or columns
192
193
194
195
196
197
198
        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)) );
          }
        }
      }
199
200
201
202
203
204
205
206
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

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

249
250
251
252
253
254
255
256
      // 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);
257
258
259
        }
      }
    }
260
#endif
261

262
263
264
  protected:
    // Estimates the slot-size used in the inserter to reserve enough space per row.
    void calculateSlotSize()
265
    {
266
      slotSize_ = 0;
267

268
269
270
271
      if (num_rows(*matrix_) != 0)
        slotSize_ = int(double(matrix_->nnz()) / num_rows(*matrix_) * 1.2);
      if (slotSize_ < MIN_NNZ_PER_ROW)
        slotSize_ = MIN_NNZ_PER_ROW;
272
    }
273

274
  private:
275
    /// The minimal number of nonzeros per row
276
    static constexpr int MIN_NNZ_PER_ROW = 50;
277

278
    /// The finite element space / basis of the row
279
    RowBasis const* rowBasis_;
280

281
    /// The finite element space / basis of the column
282
    ColBasis const* colBasis_;
283

284
    /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
285
    ClonablePtr<BaseMatrix> matrix_;
286

287
    /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
288
    Inserter* inserter_ = nullptr;
289

290
    /// A flag that indicates whether the matrix is in insertion mode
291
    bool insertionMode_ = false;
292

293
    /// The size of the slots used to insert values per row
294
    int slotSize_ = MIN_NNZ_PER_ROW;
295
296
297
  };

} // end namespace AMDiS