DOFMatrix.hpp 9.39 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
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
185
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
186
          if (dirichletNodes[row(i)]) {
187
            // set identity row
188
189
190
            value(i, (row(i) == col(i) ? value_type(1) : value_type(0)) );
          }
          else if (dirichletNodes[col(i)]) {
191
            columns.push_back({row(i), col(i), value(i)});
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
242
243
244
            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)) );
          }
        }
      }
245

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

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

271
  private:
272
273
    /// The minimal number of nonzeros per row
    static constexpr int MIN_NNZ_PER_ROW = 20;
274

275
276
    /// The finite element space / basis of the row
    RowFeSpace const& rowFeSpace;
277

278
279
    /// The finite element space / basis of the column
    ColFeSpace const& colFeSpace;
280

281
282
    /// The name of the DOFMatrix
    std::string name;
283

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

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

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

293
294
    /// The size of the slots used to insert values per row
    int slot_size = 20;
295

296
297
298
299
300
    // friend class declarations
    template <class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS