DOFMatrix.hpp 4.6 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
#include <amdis/Output.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
14
#include <amdis/linear_algebra/Common.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
23
  template <class ValueType>
  class MtlMatrix
24
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
25
  public:
26
    /// The matrix type of the underlying base matrix
27
    using BaseMatrix = mtl::compressed2D<ValueType>;
28

29
    /// The type of the elements of the DOFMatrix
30
    using value_type = ValueType;
31

32
33
34
    /// The index/size - type
    using size_type = typename BaseMatrix::size_type;

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

38
39
  public:
    /// Constructor. Constructs new BaseMatrix.
40
    MtlMatrix() = default;
41

42
    /// Return a reference to the data-matrix \ref matrix
43
    BaseMatrix& matrix()
44
    {
45
      assert( !inserter_ );
46
      return matrix_;
47
    }
48

49
    /// Return a reference to the data-matrix \ref matrix
50
    BaseMatrix const& matrix() const
51
    {
52
      assert( !inserter_ );
53
      return matrix_;
54
    }
55

56

57
    /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
58
    /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
59
    /// be created by \ref init and must be closed by \ref finish after insertion.
60
    void insert(size_type r, size_type c, value_type const& value)
61
    {
62
      test_exit_dbg(inserter_, "Inserter not initilized!");
63
64
65
      test_exit_dbg(r < num_rows(matrix_) && c < num_cols(matrix_),
          "Indices out of range [0,", num_rows(matrix_), ")x[0,", num_cols(matrix_), ")" );
      (*inserter_)[r][c] += value;
66
67
    }

68

69
    /// Create inserter. Assumes that no inserter is currently active on this matrix.
70
71
72
    template <class RowBasis, class ColBasis>
    void init(RowBasis const& rowBasis, ColBasis const& colBasis,
              bool prepareForInsertion)
73
    {
74
      test_exit(!inserter_, "Matrix already in insertion mode!");
75

76
      calculateSlotSize();
77
      matrix_.change_dim(rowBasis.dimension(), colBasis.dimension());
78
      if (prepareForInsertion) {
79
80
        set_to_zero(matrix_);
        inserter_ = new Inserter(matrix_, slotSize_);
81
      }
82
83
    }

84
    /// Delete inserter -> finish insertion. Must be called in order to fill the
85
    /// final construction of the matrix.
86
87
    void finish()
    {
88
      delete inserter_;
89
      inserter_ = nullptr;
90
    }
91

92
    /// Return the number of nonzeros in the matrix
93
    std::size_t nnz() const
94
    {
95
      return matrix_.nnz();
96
    }
97

98
99
    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
    /// a one on the diagonal of the matrix.
100
    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
101
    {
102
      // Define the property maps
103
104
105
      auto row   = mtl::mat::row_map(matrix_);
      auto col   = mtl::mat::col_map(matrix_);
      auto value = mtl::mat::value_map(matrix_);
106

107
      // iterate over the matrix
108
      for (auto r : mtl::rows_of(matrix_)) {   // rows of the matrix
109
110
111
112
113
114
115
        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)) );
          }
        }
      }
116

Praetorius, Simon's avatar
Praetorius, Simon committed
117
      return std::list<Triplet<value_type>>{};
118
    }
119

120
121
122
  protected:
    // Estimates the slot-size used in the inserter to reserve enough space per row.
    void calculateSlotSize()
123
    {
124
      slotSize_ = 0;
125

126
127
      if (num_rows(matrix_) != 0)
        slotSize_ = int(double(matrix_.nnz()) / num_rows(matrix_) * 1.2);
128
129
      if (slotSize_ < MIN_NNZ_PER_ROW)
        slotSize_ = MIN_NNZ_PER_ROW;
130
    }
131

132
  private:
133
    /// The minimal number of nonzeros per row
134
    static constexpr int MIN_NNZ_PER_ROW = 50;
135

136
    /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
137
    BaseMatrix matrix_;
138

139
    /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
140
    Inserter* inserter_ = nullptr;
141

142
    /// The size of the slots used to insert values per row
143
    int slotSize_ = MIN_NNZ_PER_ROW;
144
145
  };

146
147
148
  template <class RowBasisType, class ColBasisType, class ValueType = double>
  using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, MtlMatrix<ValueType>>;

149
} // end namespace AMDiS