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

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

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

13
14
#include <amdis/Output.hpp>
#include <amdis/common/ClonablePtr.hpp>
15
#include <amdis/linear_algebra/DOFMatrixBase.hpp>
16
#include <amdis/utility/MultiIndex.hpp>
17
18
19

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

30
    /// The type of the finite element space / basis of the column
31
    using ColBasis = ColBasisType;
32

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

36
    /// The index/size - type
37
    using size_type = typename BaseMatrix::size_type;
38

39
    /// The type of the elements of the DOFMatrix
40
    using value_type = ValueType;
41

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

45
46
  public:
    /// Constructor. Constructs new BaseMatrix.
47
    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis)
48
49
50
      : rowBasis_(&rowBasis)
      , colBasis_(&colBasis)
      , matrix_(ClonablePtr<BaseMatrix>::make())
51
    {}
52

53
54
    /// Constructor. Takes a reference to a base matrix
    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis, BaseMatrix& matrix)
55
56
      : rowBasis_(&rowBasis)
      , colBasis_(&colBasis)
57
      , matrix_(matrix)
58
    {}
59

60
61
    /// Return the row-basis \ref rowBasis_ of the matrix
    RowBasis const& rowBasis() const
62
    {
63
      return *rowBasis_;
64
    }
65

66
67
    /// Return the col-basis \ref colBasis_ of the matrix
    ColBasis const& colBasis() const
68
    {
69
      return *colBasis_;
70
    }
71

72
    /// Return a reference to the data-matrix \ref matrix
73
    BaseMatrix& matrix()
74
    {
75
      assert( !inserter_ );
76
      return *matrix_;
77
    }
78

79
    /// Return a reference to the data-matrix \ref matrix
80
    BaseMatrix const& matrix() const
81
    {
82
      assert( !inserter_ );
83
      return *matrix_;
84
    }
85

86
    /// Return the size of the row \ref feSpace
87
    size_type rows() const
88
    {
89
      return rowBasis_->size();
90
    }
91

92
    /// Return the size of the column \ref feSpace
93
    size_type cols() const
94
    {
95
      return colBasis_->size();
96
    }
97

98

99
    /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
100
    /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
101
    /// be created by \ref init and must be closed by \ref finish after insertion.
102
103
    template <class Index>
    auto operator()(Index row, Index col)
104
    {
105
106
107
108
      size_type r = flatMultiIndex(row), c = flatMultiIndex(col);
      test_exit_dbg(inserter_, "Inserter not initilized!");
      test_exit_dbg(r < rows() && c < cols(),
          "Indices out of range [0,", rows(), ")x[0,", cols(), ")" );
109
      return (*inserter_)[r][c];
110
111
    }

112

113
    /// Create inserter. Assumes that no inserter is currently active on this matrix.
114
    void init(bool prepareForInsertion)
115
    {
116
      test_exit(!inserter_, "Matrix already in insertion mode!");
117

118
      calculateSlotSize();
119
120
      matrix_->change_dim(rowBasis_->dimension(), colBasis_->dimension());
      test_exit(num_rows(*matrix_) == rowBasis_->dimension() && num_cols(*matrix_) == colBasis_->dimension(),
121
        "Wrong dimensions in matrix");
122
      if (prepareForInsertion) {
123
124
        set_to_zero(*matrix_);
        inserter_ = new Inserter(*matrix_, slotSize_);
125
      }
126
127
    }

128
    /// Delete inserter -> finish insertion. Must be called in order to fill the
129
    /// final construction of the matrix.
130
131
    void finish()
    {
132
      delete inserter_;
133
      inserter_ = nullptr;
134
    }
135

136
    /// Return the number of nonzeros in the matrix
137
    std::size_t nnz() const
138
    {
139
      return matrix_->nnz();
140
    }
141

142
143
    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
    /// a one on the diagonal of the matrix.
144
    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
145
    {
146
      std::list<Triplet<value_type>> columns;
147

148
      // Define the property maps
149
150
151
      auto row   = mtl::mat::row_map(*matrix_);
      auto col   = mtl::mat::col_map(*matrix_);
      auto value = mtl::mat::value_map(*matrix_);
152

153
      // iterate over the matrix
154
#if 0
155
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
156
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
157
          if (dirichletNodes[row(i)]) {
158
            // set identity row
159
160
161
            value(i, (row(i) == col(i) ? value_type(1) : value_type(0)) );
          }
          else if (dirichletNodes[col(i)]) {
162
            columns.push_back({row(i), col(i), value(i)});
163
164
165
166
            value(i, value_type(0));
          }
        }
      }
167
168
#endif

169
      for (auto r : mtl::rows_of(*matrix_)) {   // rows or columns
170
171
172
173
174
175
176
        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)) );
          }
        }
      }
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
222
223
224
225

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

227
228
229
230
231
232
233
234
      // 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);
235
236
237
        }
      }
    }
238
#endif
239

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

246
247
248
249
      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;
250
    }
251

252
  private:
253
    /// The minimal number of nonzeros per row
254
    static constexpr int MIN_NNZ_PER_ROW = 50;
255

256
    /// The finite element space / basis of the row
257
    RowBasis const* rowBasis_;
258

259
    /// The finite element space / basis of the column
260
    ColBasis const* colBasis_;
261

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

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

268
    /// The size of the slots used to insert values per row
269
    int slotSize_ = MIN_NNZ_PER_ROW;
270
271
272
  };

} // end namespace AMDiS