DOFMatrix.hpp 6.19 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#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>

#include <dune/amdis/ClonablePtr.hpp>
#include <dune/amdis/Log.hpp>

namespace AMDiS
{
16
17
  /// \brief The basic container that stores a base matrix and a corresponding 
  /// row/column feSpace
18
19
20
21
22
  template <class RowFeSpaceType, class ColFeSpaceType, 
            class ValueType = double>
  class DOFMatrix
  {
  public:
23
    /// The type of the finite element space / basis of the row
24
    using RowFeSpace = RowFeSpaceType;
25
26
    
    /// The type of the finite element space / basis of the column
27
    using ColFeSpace = ColFeSpaceType;
28
29
    
    /// The vector type of the underlying base matrix
30
31
    using BaseMatrix = mtl::compressed2D<ValueType>;
    
32
    /// The index/size - type
33
34
    using size_type  = typename BaseMatrix::size_type;
    
35
    /// The type of the elements of the DOFMatrix
36
    using value_type = ValueType;
37
38
    
    /// The underlying field-type (typically the same as \ref value_type)
39
40
    using field_type = typename BaseMatrix::value_type;
    
41
    /// The type of the inserter used when filling the matrix. NOTE: need not be public
42
43
    using Inserter = mtl::matrix::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;
    
44
45
46
47
48
  public:
    /// Constructor. Constructs new BaseMatrix.
    DOFMatrix(RowFeSpace const& rowFeSpace, 
              ColFeSpace const& colFeSpace, 
              std::string name)
49
50
51
52
53
54
      : rowFeSpace(rowFeSpace)
      , colFeSpace(colFeSpace)
      , name(name)
      , matrix(ClonablePtr<BaseMatrix>::make())
    {}
    
55
56
57
58
    /// Constructor. Takes pointer of data-matrix.
    DOFMatrix(RowFeSpace const& rowFeSpace, 
              ColFeSpace const& colFeSpace, 
              std::string name, 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
              BaseMatrix& matrix_ref)
      : rowFeSpace(rowFeSpace)
      , colFeSpace(colFeSpace)
      , name(name)
      , matrix(matrix_ref)
    {}
    
    /// Return the row-basis \ref rowFeSpace of the matrix
    RowFeSpace const& getRowFeSpace() const
    {
      return rowFeSpace;
    }
    
    /// Return the col-basis \ref colFeSpace of the matrix
    ColFeSpace const& getColFeSpace() const
    {
      return colFeSpace;
    }
    
78
79
    /// Return a reference to the data-matrix \ref matrix
    BaseMatrix& getMatrix()
80
81
82
83
    {
      return *matrix;
    }
    
84
85
    /// Return a reference to the data-matrix \ref matrix
    BaseMatrix const& getMatrix() const
86
87
88
89
    {
      return *matrix;
    }
    
90
    /// Return the size of the row \ref feSpace
91
92
93
94
95
    size_type N() const
    {
      return rowFeSpace.size();
    }
    
96
    /// Return the size of the column \ref feSpace
97
98
99
100
101
    size_type M() const
    {
      return colFeSpace.size();
    }
    
102
    /// Return the \ref name of this matrix
103
104
105
106
107
    std::string getName() const
    {
      return name;
    }
    
108
109
110
111
    /// \brief Returns an update-proxy of the inserter, to inserte/update a value at 
    /// 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.
    auto operator()(size_type r, size_type c)
112
113
114
115
116
117
118
    {
      AMDIS_TEST_EXIT_DBG( initialized, "Inserter not initialized!");
      AMDIS_TEST_EXIT_DBG( r < N() && c < M() , 
          "Indices out of range [0," << N() << ")x[0," << M() << ")" );
      return (*inserter)[r][c];
    }
    
119
    /// Create inserter. Assumes that no inserter is currently active on this matrix.
120
121
122
123
    void init()
    {
      AMDIS_TEST_EXIT(!initialized && !inserter, "Matrix already in insertion mode!");
      
124
      calculateSlotSize();
125
126
127
128
129
130
      matrix->change_dim(rowFeSpace.size(), colFeSpace.size());
      set_to_zero(*matrix);
      inserter = new Inserter(*matrix, slot_size);
      initialized = true;
    }

131
132
    /// Delete inserter -> finish insertion. Must be called in order to fill the 
    /// final construction of the matrix.
133
134
135
136
137
138
139
    void finish()
    {
      delete inserter;
      inserter = NULL;
      initialized = false;
    }
    
140
    /// Return the number of nonzeros in the matrix
141
142
143
144
145
    size_t getNnz() const
    {
      return matrix->nnz();
    }
    
146
147
    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
    /// a one on the diagonal of the matrix.
148
149
150
    void clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
    {    
      // Define the property maps
151
152
153
      auto row   = mtl::matrix::row_map(*matrix); 
      auto col   = mtl::matrix::col_map(*matrix);
      auto value = mtl::matrix::value_map(*matrix); 
154

155
      // iterate over the matrix    
156
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
157
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
158
159
160
161
162
163
164
165
166
          if (!dirichletNodes[row(i)])
            break;
          
          // set identity row
          value(i, (apply && row(i) == col(i) ? value_type(1) : value_type(0)) );
        }
      }
    }
    
167
168
169
  protected:
    // Estimates the slot-size used in the inserter to reserve enough space per row.
    void calculateSlotSize()
170
171
172
173
174
175
176
177
178
179
    {
      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;
    }
      
  private:
180
181
182
183
184
    /// The minimal number of nonzeros per row
    static constexpr int MIN_NNZ_PER_ROW = 20;
    
    /// The finite element space / basis of the row
    RowFeSpace const& rowFeSpace;
185
    
186
187
    /// The finite element space / basis of the column
    ColFeSpace const& colFeSpace;
188
    
189
190
191
192
193
194
195
196
197
198
    /// The name of the DOFMatrix
    std::string name;
    
    /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
    ClonablePtr<BaseMatrix> matrix;
    
    /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
    Inserter* inserter = NULL;
    
    /// A flag that indicates whether the matrix is in insertion mode
199
    bool initialized = false;
200
201
202
    
    /// The size of the slots used to insert values per row
    int slot_size = 20;
203
204
205
206
207
208
    
    // friend class declarations
    template <class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS