DOFMatrix.hpp 9.28 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
87
    {
      return *matrix;
    }
88

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

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

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

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

113
114
115
116

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

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

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

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

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

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

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

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

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

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

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

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

270
271
    /// The finite element space / basis of the row
    RowFeSpace const& rowFeSpace;
272

273
274
    /// The finite element space / basis of the column
    ColFeSpace const& colFeSpace;
275

276
277
    /// The name of the DOFMatrix
    std::string name;
278

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

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

285
    /// A flag that indicates whether the matrix is in insertion mode
286
    bool initialized = false;
287

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

291
292
293
294
295
    // friend class declarations
    template <class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS