DOFMatrix.hpp 4.47 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
#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
{
  template <class RowFeSpaceType, class ColFeSpaceType, 
            class ValueType = double>
  class DOFMatrix
  {
  public:
    using RowFeSpace = RowFeSpaceType;
    using ColFeSpace = ColFeSpaceType;
    using BaseMatrix = mtl::compressed2D<ValueType>;
    
    using size_type  = typename BaseMatrix::size_type;
    
    using value_type = ValueType;
    using field_type = typename BaseMatrix::value_type;
    
    using Inserter = mtl::matrix::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;
    
    static constexpr int MIN_NNZ_PER_ROW = 20;
    
    /// Constructor. Constructs new BaseVector.
    DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name)
      : rowFeSpace(rowFeSpace)
      , colFeSpace(colFeSpace)
      , name(name)
      , matrix(ClonablePtr<BaseMatrix>::make())
    {}
    
    /// Constructor. Takes pointer of data-vector.
    DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name, 
              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;
    }
    
    /// Return the data-vector \ref vector
    BaseMatrix const& getMatrix() const
    {
      return *matrix;
    }
    
    /// Return the data-vector \ref vector
    BaseMatrix& getMatrix()
    {
      return *matrix;
    }
    
    /// Return the size of the \ref feSpace
    size_type N() const
    {
      return rowFeSpace.size();
    }
    
    size_type M() const
    {
      return colFeSpace.size();
    }
    
    /// Return the \ref name of this vector
    std::string getName() const
    {
      return name;
    }
    
    /// Access the entry \p i of the \ref vector with write-access.
    auto operator()(size_type r, size_type c) // returns update_proxy
    {
      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];
    }
    
    /// create inserter
    void init()
    {
      AMDIS_TEST_EXIT(!initialized && !inserter, "Matrix already in insertion mode!");
      
      calculateNnz();
      matrix->change_dim(rowFeSpace.size(), colFeSpace.size());
      set_to_zero(*matrix);
      inserter = new Inserter(*matrix, slot_size);
      initialized = true;
    }

    // delete inserter -> finish insertion
    void finish()
    {
      delete inserter;
      inserter = NULL;
      initialized = false;
    }
    
    size_t getNnz() const
    {
      return matrix->nnz();
    }
    
    void clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
    {    
      // Define the property maps
      auto row=   mtl::matrix::row_map(*matrix); 
      auto col=   mtl::matrix::col_map(*matrix);
      auto value= mtl::matrix::value_map(*matrix); 

      // Now iterate over the matrix    
      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
        for (auto i : mtl::nz_of(r)) {    // non-zeros within
          if (!dirichletNodes[row(i)])
            break;
          
          // set identity row
          value(i, (apply && row(i) == col(i) ? value_type(1) : value_type(0)) );
        }
      }
    }
    
  protected:    
    void calculateNnz()
    {
      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:
    RowFeSpace const&	rowFeSpace;
    ColFeSpace const&	colFeSpace;
    std::string	        name;
    
    ClonablePtr<BaseMatrix>   matrix;
    Inserter*                 inserter = NULL;
    
    bool initialized = false;
    int slot_size    = 20;
    
    // friend class declarations
    template <class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS