DOFMatrix.hpp 4.96 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
172
173
174
175
176
177
178
#pragma once

#include <string>
#include <memory>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>

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

namespace AMDiS
{
  template <class RowFeSpaceType, class ColFeSpaceType, 
            class ValueType = Dune::FieldMatrix<double,1,1>>
  class DOFMatrix
  {
  public:
    using RowFeSpace = RowFeSpaceType;
    using ColFeSpace = ColFeSpaceType;
    using BaseMatrix = Dune::BCRSMatrix<ValueType>;
    
    using size_type  = typename RowFeSpaceType::size_type;
    using field_type = typename ValueType::field_type;
    
    using value_type = ValueType;
    
    /// 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 read-access.
    value_type const& operator()(size_type r, size_type c) const
    {
      AMDIS_TEST_EXIT_DBG( initialized, "Occupation pattern not initialized!");
      AMDIS_TEST_EXIT_DBG( r < N() && c < M() , 
          "Indices out of range [0," << N() << ")x[0," << M() << ")" );
      return (*matrix)[r][c];
    }
    
    /// Access the entry \p i of the \ref vector with write-access.
    value_type& operator()(size_type r, size_type c)
    {
      AMDIS_TEST_EXIT_DBG( initialized, "Occupation pattern not initialized!");
      AMDIS_TEST_EXIT_DBG( r < N() && c < M() , 
          "Indices out of range [0," << N() << ")x[0," << M() << ")" );
      return (*matrix)[r][c];
    }
    
    /// create occupation pattern and apply it to the matrix
    void init()
    {
      occupationPattern.resize(rowFeSpace.size(), colFeSpace.size());
      auto meshView = rowFeSpace.gridView();

      // A loop over all elements of the grid
      auto rowLocalView = rowFeSpace.localView();
      auto rowIndexSet  = rowFeSpace.localIndexSet();
      
      auto colLocalView = colFeSpace.localView();
      auto colIndexSet  = colFeSpace.localIndexSet();

      for (const auto& element : elements(meshView)) {
        rowLocalView.bind(element);
        rowIndexSet.bind(rowLocalView); // TODO: can this be moved out of the loop?
        
        colLocalView.bind(element);
        colIndexSet.bind(colLocalView); // TODO: can this be moved out of the loop?
        
        for (size_t i = 0; i < rowIndexSet.size(); ++i) {
          // The global index of the i-th vertex of the element
          auto row = rowIndexSet.index(i);
          for (size_t j = 0; j < colIndexSet.size(); ++j) {
            // The global index of the j-th vertex of the element
            auto col = colIndexSet.index(j);
            occupationPattern.add(row, col);
          }
        }
      }
      occupationPattern.exportIdx(*matrix);
      
      initialized = true;
    }

    void finish()
    {
      initialized = false;
    }
    
    size_t getNnz()
    {
      return matrix->nonzeroes();
    }
    
    void clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
    {    
      // loop over the matrix rows
      for (size_t i = 0; i < matrix->N(); ++i) {
        if (dirichletNodes[i]) {
            auto cIt = (*matrix)[i].begin();
            auto cEndIt = (*matrix)[i].end();
            // loop over nonzero matrix entries in current row
            for (; cIt != cEndIt; ++cIt)
                *cIt = (apply && i == cIt.index() ? 1 : 0);
        }
      }
    }
      
  private:
    RowFeSpace const&	rowFeSpace;
    ColFeSpace const&	colFeSpace;
    std::string	        name;
    
    ClonablePtr<BaseMatrix> matrix;
    Dune::MatrixIndexSet    occupationPattern;
    
    bool initialized = false;
    
    // friend class declarations
    template <class, class, class> friend class SystemMatrix;
  };

} // end namespace AMDiS