Mapper.hpp 6.95 KB
Newer Older
1
2
3
4
5
6
#pragma once

#include <vector>

#include <boost/numeric/mtl/matrix/mapped_inserter.hpp>

7
8
#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
9

10
namespace AMDiS
11
{
12

13
14
  /**
   * \brief BaseClass for all mapper types
15
   *
16
17
18
19
20
   * A mapper assigned a local row/column index of a cell of
   * a block matrix to a global matrix index.
   * Call \ref setRow and \ref setCol first, to select a block
   * in the block matrix. Then get with \ref row, respective \ref col,
   * the global matrix index to the assigned local index.
21
   **/
22
23
24
25
  template <class Derived>
  struct MapperBase
  {
    using size_type = std::size_t;
26

27
28
29
30
    template <class BaseInserter>
    struct Inserter {
      using type = mtl::mat::mapped_inserter<BaseInserter, Derived>;
    };
31

32
33
    /// set the current block row
    void setRow(size_type rb) { self().setRow(rb); }
34

35
36
    /// set the current block columns
    void setCol(size_type cb) { self().setCol(cb); }
37

38
39
    /// return global matrix row, for local row in the current matrix block
    size_type row(size_type r) const { return self().row(r); }
40

41
42
    /// return global matrix column, for local column in the current matrix block
    size_type col(size_type c) const { return self().col(c); }
43

44
45
    /// return overall number of rows
    size_type getNumRows() const { return self().getNumRows(); }
46

47
48
49
50
51
    /// return overall number of columns
    size_type getNumCols() const { return self().getNumCols(); }

    size_type getNumRows(size_type comp) const { return self().getNumRows(comp); }
    size_type getNumCols(size_type comp) const { return self().getNumCols(comp); }
52

53
54
    /// return number of components/blocks
    size_type getNumComponents() const { return self().getNumComponents(); }
55

56
57
58
59
60
61
    size_type row(size_type rb, size_type cb, size_type r)
    {
      setRow(rb);
      setCol(cb);
      return row(r);
    }
62

63
64
65
66
67
68
    size_type col(size_type rb, size_type cb, size_type c)
    {
      setRow(rb);
      setCol(cb);
      return col(c);
    }
69

70
71
72
73
    Derived&       self()       { return static_cast<Derived&>(*this); }
    Derived const& self() const { return static_cast<Derived const&>(*this); }
  };

74

75
76
77
78
79
80
81
  /// Mapper implementation for non-parallel block matrices
  struct BlockMapper : public MapperBase<BlockMapper>
  {
    /// Default constructor
    BlockMapper() = default;
    BlockMapper(BlockMapper const& other) = default;
    BlockMapper& operator=(BlockMapper const& other) = default;
82

83
    /// Constructor for BlockMTLMatrix
84
    template <class Value, std::size_t _N>
85
    explicit BlockMapper(BlockMTLMatrix<Value, _N, _N> const& mat)
86
      : nComp(_N)
87
88
89
      , rowOffset(0), colOffset(0)
      , nrow(num_rows(mat)), ncol(nrow)
      , sizes(nComp)
90
    {
91
      for (std::size_t i = 0; i < _N; ++i)
92
93
        sizes[i] = num_rows(mat[i][i]);
    }
94

95
96
    /// Constructor for SystemMatrix
    template <class FeSpaces, class Value>
97
    explicit BlockMapper(SystemMatrix<FeSpaces, Value> const& mat)
98
99
      : BlockMapper(mat.getMatrix())
    {}
100

101
102
    /// Constructor for single DOFMatrix
    template <class FeSpace, class Value>
103
    explicit BlockMapper(DOFMatrix<FeSpace, Value> const& mat)
104
      : nComp(1)
105
106
107
108
109
110
111
112
113
114
115
      , rowOffset(0), colOffset(0)
      , nrow(0), ncol(0)
      , sizes(nComp)
    {
      sizes[0] = mat.N();
      nrow += sizes[0];
      ncol = nrow;
    }

    /// Constructor for system with equal components
    BlockMapper(size_type nComp, size_type nDOFperComp)
116
117
    : nComp(nComp),
      rowOffset(0), colOffset(0),
118
119
120
121
122
123
      nrow(nComp*nDOFperComp), ncol(nrow),
      sizes(nComp)
    {
       for (size_type i = 0; i < nComp; ++i)
         sizes[i] = nDOFperComp;
    }
124

125
126
    /// calculate row offset for row component \param r
    void setRow(size_type rb)
127
128
129
    {
      test_exit_dbg(rb <= sizes.size(), "row nr out of range!");
      rowOffset = sum(rb);
130
131
132
133
134
    }

    /// calculate column offset for col component \param c
    void setCol(size_type cb)
    {
135
      test_exit_dbg(cb <= sizes.size(), "column nr out of range!");
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
      colOffset = sum(cb);
    }

    size_type row(size_type r) const { return r + rowOffset; }
    size_type col(size_type c) const { return c + colOffset; }

    size_type getNumRows() const { return nrow; }
    size_type getNumCols() const { return ncol; }

    size_type getNumRows(size_type comp) const { return sizes[comp]; }
    size_type getNumCols(size_type comp) const { return sizes[comp]; }

    size_type getNumComponents() const { return nComp; }

  private: // methods

    ///compute the sum of sizes from [0, end)
153
    size_type sum(size_type end) const
154
155
156
157
158
159
    {
      size_type ret = 0;
      for (size_type i = 0; i < end; ++i)
	ret += sizes[i];
      return ret;
    }
160

161
162
163
164
165
166
167
  private: // variables

    size_type nComp = 0;
    size_type rowOffset = 0;
    size_type colOffset = 0;
    size_type nrow = 0;
    size_type ncol = 0;
168

169
170
171
    std::vector<size_type> sizes;
  };

172

173
174
175
176
  /// Mapper implementation for non-parallel rectangular block matrices
  struct RectangularMapper : public MapperBase<RectangularMapper>
  {
    /// Constructor for block-matrices
177
    template <class Value, std::size_t _N, std::size_t _M>
178
    explicit RectangularMapper(BlockMTLMatrix<Value, _N, _M> const& mat)
179
180
181
182
183
184
185
186
      : nRowComp(_N), nColComp(_M)
      , rowOffset(0), colOffset(0)
      , nrow(num_rows(mat)), ncol(num_cols(mat))
      , sizes_rows(nRowComp), sizes_cols(nColComp)
    {
      for (size_type i = 0; i < _N; ++i)
        sizes_rows[i] = num_rows(mat[i][0]);
      for (size_type j = 0; j < _M; ++j)
187
        sizes_cols[j] = num_cols(mat[0][j]);
188
    }
189

190
191
    /// calculate row offset for row component \param r
    void setRow(size_type rb)
192
193
194
    {
      test_exit_dbg(rb <= sizes_rows.size(), "row nr out of range!");
      rowOffset = sum(rb, sizes_rows);
195
196
197
198
199
    }

    /// calculate column offset for col component \param c
    void setCol(size_type cb)
    {
200
      test_exit_dbg(cb <= sizes_cols.size(), "column nr out of range!");
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
      colOffset = sum(cb, sizes_cols);
    }

    size_type row(size_type r) const { return r + rowOffset; }
    size_type col(size_type c) const { return c + colOffset; }

    size_type getNumRows() const { return nrow; }
    size_type getNumCols() const { return ncol; }

    size_type getNumRows(size_type comp) const { return sizes_rows[comp]; }
    size_type getNumCols(size_type comp) const { return sizes_cols[comp]; }

    size_type getNumComponents() const { return nRowComp; }
    size_type getNumRowComponents() const { return nRowComp; }
    size_type getNumColComponents() const { return nColComp; }

  private: // methods

    ///compute the sum of sizes from [0, end)
220
    size_type sum(size_type end, std::vector<size_type> const& sizes) const
221
222
223
224
225
226
    {
      size_type ret = 0;
      for (size_type i = 0; i < end; ++i)
	ret += sizes[i];
      return ret;
    }
227

228
229
230
231
232
233
234
235
  private: // variables

    size_type nRowComp = 0;
    size_type nColComp = 0;
    size_type rowOffset = 0;
    size_type colOffset = 0;
    size_type nrow = 0;
    size_type ncol = 0;
236

237
238
239
    std::vector<size_type> sizes_rows;
    std::vector<size_type> sizes_cols;
  };
240

241
} // end namespace AMDiS