DOFMatrix.hpp 9.35 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
#include <dune/common/std/optional.hh>
12
13
#include <dune/common/reservedvector.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
14

15
#include <dune/amdis/Output.hpp>
16
#include <dune/amdis/common/ClonablePtr.hpp>
17
18
19

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

31
    /// The type of the finite element space / basis of the column
32
    using ColFeSpace = ColFeSpaceType;
33

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

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

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

43
    /// The underlying field-type (typically the same as \ref value_type)
44
    using field_type = typename BaseMatrix::value_type;
45

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

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

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

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

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

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

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

97
    /// Return the size of the row \ref feSpace
98
99
100
101
    size_type N() const
    {
      return rowFeSpace.size();
    }
102

103
    /// Return the size of the column \ref feSpace
104
105
106
107
    size_type M() const
    {
      return colFeSpace.size();
    }
108

109
    /// Return the \ref name of this matrix
110
111
112
113
    std::string getName() const
    {
      return name;
    }
114

115
116
117
118

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

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

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

145
      calculateSlotSize();
146
      matrix->change_dim(rowFeSpace.size(), colFeSpace.size());
147
148
149
      if (prepareForInsertion) {
        set_to_zero(*matrix);
        inserter = new Inserter(*matrix, slot_size);
150
        insertionMode = true;
151
      }
152
153
    }

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

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

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

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

180
      // iterate over the matrix
181
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
182
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
183
          if (dirichletNodes[row(i)]) {
184
            // set identity row
185
186
187
            value(i, (row(i) == col(i) ? value_type(1) : value_type(0)) );
          }
          else if (dirichletNodes[col(i)]) {
188
            columns.push_back({row(i), col(i), value(i)});
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
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
            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)) );
          }
        }
      }
242

243
244
245
246
247
248
249
250
      // 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);
251
252
253
        }
      }
    }
254
#endif
255

256
257
258
  protected:
    // Estimates the slot-size used in the inserter to reserve enough space per row.
    void calculateSlotSize()
259
260
261
262
263
264
265
266
    {
      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;
    }
267

268
  private:
269
270
    /// The minimal number of nonzeros per row
    static constexpr int MIN_NNZ_PER_ROW = 20;
271

272
273
    /// The finite element space / basis of the row
    RowFeSpace const& rowFeSpace;
274

275
276
    /// The finite element space / basis of the column
    ColFeSpace const& colFeSpace;
277

278
279
    /// The name of the DOFMatrix
    std::string name;
280

281
282
    /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
    ClonablePtr<BaseMatrix> matrix;
283

284
285
    /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
    Inserter* inserter = NULL;
286

287
    /// A flag that indicates whether the matrix is in insertion mode
288
    bool insertionMode = false;
289

290
291
    /// The size of the slots used to insert values per row
    int slot_size = 20;
292

293
294
295
296
297
    // friend class declarations
    template <class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS