DOFMatrix.hpp 4.99 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>
Praetorius, Simon's avatar
Praetorius, Simon committed
15
#include <amdis/linear_algebra/Common.hpp>
16
#include <amdis/linear_algebra/DOFMatrixBase.hpp>
17
#include <amdis/utility/MultiIndex.hpp>
18
19
20

namespace AMDiS
{
21
  /// \brief The basic container that stores a base matrix and a corresponding
22
  /// row/column feSpace.
23
  template <class RowBasisType, class ColBasisType,
24
25
            class ValueType = double>
  class DOFMatrix
Praetorius, Simon's avatar
Praetorius, Simon committed
26
      : public DOFMatrixBase<RowBasisType, ColBasisType>
27
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
28
    using Super = DOFMatrixBase<RowBasisType, ColBasisType>;
29

Praetorius, Simon's avatar
Praetorius, Simon committed
30
  public:
31
    /// The matrix type of the underlying base matrix
32
    using BaseMatrix = mtl::compressed2D<ValueType>;
33

34
    /// The type of the elements of the DOFMatrix
35
    using value_type = ValueType;
36

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

40
41
  public:
    /// Constructor. Constructs new BaseMatrix.
Praetorius, Simon's avatar
Praetorius, Simon committed
42
43
    DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis)
      : Super(rowBasis, colBasis)
44
      , matrix_(ClonablePtr<BaseMatrix>::make())
45
    {}
46

47
    /// Constructor. Takes a reference to a base matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
48
49
    DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis, BaseMatrix& matrix)
      : Super(rowBasis, colBasis)
50
      , matrix_(matrix)
51
    {}
52

53
    /// Return a reference to the data-matrix \ref matrix
54
    BaseMatrix& matrix()
55
    {
56
      assert( !inserter_ );
57
      return *matrix_;
58
    }
59

60
    /// Return a reference to the data-matrix \ref matrix
61
    BaseMatrix const& matrix() const
62
    {
63
      assert( !inserter_ );
64
      return *matrix_;
65
    }
66

67

68
    /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
69
    /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
70
    /// be created by \ref init and must be closed by \ref finish after insertion.
71
72
    template <class Index>
    auto operator()(Index row, Index col)
73
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
74
      auto r = flatMultiIndex(row), c = flatMultiIndex(col);
75
      test_exit_dbg(inserter_, "Inserter not initilized!");
Praetorius, Simon's avatar
Praetorius, Simon committed
76
77
      test_exit_dbg(r < this->rows() && c < this->cols(),
          "Indices out of range [0,", this->rows(), ")x[0,", this->cols(), ")" );
78
      return (*inserter_)[r][c];
79
80
    }

81

82
    /// Create inserter. Assumes that no inserter is currently active on this matrix.
83
    void init(bool prepareForInsertion)
84
    {
85
      test_exit(!inserter_, "Matrix already in insertion mode!");
86

87
      calculateSlotSize();
Praetorius, Simon's avatar
Praetorius, Simon committed
88
89
      matrix_->change_dim(this->rows(), this->cols());
      test_exit(num_rows(*matrix_) == this->rows() && num_cols(*matrix_) == this->cols(),
90
        "Wrong dimensions in matrix");
91
      if (prepareForInsertion) {
92
93
        set_to_zero(*matrix_);
        inserter_ = new Inserter(*matrix_, slotSize_);
94
      }
95
96
    }

97
    /// Delete inserter -> finish insertion. Must be called in order to fill the
98
    /// final construction of the matrix.
99
100
    void finish()
    {
101
      delete inserter_;
102
      inserter_ = nullptr;
103
    }
104

105
    /// Return the number of nonzeros in the matrix
106
    std::size_t nnz() const
107
    {
108
      return matrix_->nnz();
109
    }
110

111
112
    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
    /// a one on the diagonal of the matrix.
113
    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
114
    {
115
      // Define the property maps
116
117
118
      auto row   = mtl::mat::row_map(*matrix_);
      auto col   = mtl::mat::col_map(*matrix_);
      auto value = mtl::mat::value_map(*matrix_);
119

120
      // iterate over the matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
121
      for (auto r : mtl::rows_of(*matrix_)) {   // rows of the matrix
122
123
124
125
126
127
128
        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)) );
          }
        }
      }
129

Praetorius, Simon's avatar
Praetorius, Simon committed
130
      return std::list<Triplet<value_type>>{};
131
    }
132

133
134
135
  protected:
    // Estimates the slot-size used in the inserter to reserve enough space per row.
    void calculateSlotSize()
136
    {
137
      slotSize_ = 0;
138

139
140
141
142
      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;
143
    }
144

145
  private:
146
    /// The minimal number of nonzeros per row
147
    static constexpr int MIN_NNZ_PER_ROW = 50;
148

149
    /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
150
    ClonablePtr<BaseMatrix> matrix_;
151

152
    /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
153
    Inserter* inserter_ = nullptr;
154

155
    /// The size of the slots used to insert values per row
156
    int slotSize_ = MIN_NNZ_PER_ROW;
157
158
159
  };

} // end namespace AMDiS