SystemMatrix.hpp 10.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#pragma once

#include <vector>
#include <string>
#include <memory>
#include <tuple>

#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/istl/multitypeblockvector.hh>

11
12
13
14
#include <amdis/Output.hpp>
#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/Utility.hpp>
#include <amdis/linear_algebra/istl/DOFMatrix.hpp>
15
16

// for explicit instantiation
17
#include <amdis/ProblemStatTraits.hpp>
18
19
20
21
22
23
24
25
26

namespace AMDiS
{
  namespace Impl
  {
    // DOFMatrices = std::tuple< std::tuple< DOFMatrix<RowFeSpace, ColFeSpace, ValueType>... >... >
    template <class DOFMatrices>
    struct BuildDOFMatrices
    {
27
      template <std::size_t R>
28
      using DOFMatrixRow = std::tuple_element_t<R, DOFMatrices>;
29

30
      template <std::size_t R, std::size_t C>
31
      using DOFMatrixRowCol = std::tuple_element_t<C, DOFMatrixRow<R>>;
32

33
      // add arg to repeated constructor argument list
34
      template <std::size_t... R, class RowFeSpaces, class ColFeSpaces, class MultiMatrix>
35
36
      static DOFMatrices make(Indices<R...>,
                              RowFeSpaces&& rowFeSpaces,
37
38
39
40
                              ColFeSpaces&& colFeSpace,
                              MultiMatrix&& multiMatrix)
      {
        return DOFMatrices{make_row(
41
          index_<R>,
42
43
44
          MakeSeq_t<std::tuple_size<std::decay_t<ColFeSpaces>>::value>(),
          std::get<R>(std::forward<RowFeSpaces>(rowFeSpaces)),
          std::forward<ColFeSpaces>(colFeSpace),
45
          std::get<R>(std::forward<MultiMatrix>(multiMatrix))
46
47
          )...};
      }
48

49
50
      template <std::size_t R, std::size_t... C, class RowFeSpace, class ColFeSpaces, class MultiVector>
      static DOFMatrixRow<R> make_row(index_t<R>,
51
52
                                      Indices<C...>,
                                      RowFeSpace&& rowFeSpace,
53
54
55
56
                                      ColFeSpaces&& colFeSpace,
                                      MultiVector&& multiVector)
      {
        return DOFMatrixRow<R>{DOFMatrixRowCol<R,C>(
57
58
          std::forward<RowFeSpace>(rowFeSpace),
          std::get<C>(std::forward<ColFeSpaces>(colFeSpace)),
59
          "matrix[" + std::to_string(R) + "][" + std::to_string(C) + "]",
60
          std::get<C>(std::forward<MultiVector>(multiVector))
61
62
63
          )...};
      }
    };
64

65
    template <class DOFMatrices, class RowFeSpaces, class ColFeSpaces, class MultiMatrix>
66
    DOFMatrices buildDOFMatrices(RowFeSpaces&& rowFeSpaces,
67
68
69
70
71
72
73
74
75
                                  ColFeSpaces&& colFeSpace,
                                  MultiMatrix&& multiMatrix)
    {
        return BuildDOFMatrices<DOFMatrices>::make(
            MakeSeq_t<std::tuple_size<std::decay_t<RowFeSpaces>>::value>(), // traverse rows first
            std::forward<RowFeSpaces>(rowFeSpaces),
            std::forward<ColFeSpaces>(colFeSpace),
            std::forward<MultiMatrix>(multiMatrix));
    }
76

77
  } // end namespace Impl
78
79
80
81




82
  /// \brief Container that repesents multiple data-Vectors of different value types.
83
84
85
  /**
    *  Represents a \ref Dune::MultiTypeBlockVector + a tuple of corresponding
    *  feSpaces. This I'th \ref Dune::BlockVector compbined with the I'th FeSpace
86
    *  builds a \ref DOFVector  and can be return by the \ref operator[].
87
    *
88
89
    *  By default, the ValueTypes are all \ref Dune::FieldMatrix of 1x1 double entry.
    **/
90
  template <class FeSpaces,
91
            class FieldType  = double,
92
            class Dimensions = Repeat_t<std::tuple_size<FeSpaces>::value, index_t<1>> >
93
94
95
96
97
98
99
100
  class SystemMatrix
  {
    template <class RowFeSpace, class RowDim>
    struct MultiMatrixRowGenerator
    {
      template <class ColDim>
      using ValueTypeGenerator = Dune::FieldMatrix<FieldType, RowDim::value, ColDim::value>;
      using ValueTypes = MakeTuple_t<ValueTypeGenerator, Dimensions>;
101

102
103
104
      template <class ColFeSpace, class VT>
      using DOFMatrixGenerator = DOFMatrix<RowFeSpace, ColFeSpace, VT>;
      using DOFMatrices  = MakeTuple2_t<DOFMatrixGenerator, FeSpaces, ValueTypes>;
105
106

      template <class VT>
107
108
      using BaseMatrixGenerator = Dune::BCRSMatrix<VT, std::allocator<VT>>;
      using BaseMatrices = MakeTuple_t<BaseMatrixGenerator, ValueTypes>;
109

110
111
      using MultiVector = ExpandTuple_t<Dune::MultiTypeBlockVector, BaseMatrices>;
    };
112

113
114
    template <class RowFeSpace, class RowDim>
    using MultiVectorGenerator = typename MultiMatrixRowGenerator<RowFeSpace, RowDim>::MultiVector;
115

116
117
    template <class RowFeSpace, class RowDim>
    using ValueTypeGenerator = typename MultiMatrixRowGenerator<RowFeSpace, RowDim>::ValueTypes;
118

119
120
    template <class RowFeSpace, class RowDim>
    using DOFMatrixGenerator = typename MultiMatrixRowGenerator<RowFeSpace, RowDim>::DOFMatrices;
121

122
123
124
125
126
  public:
    using MultiVectors = MakeTuple2_t<MultiVectorGenerator, FeSpaces, Dimensions>;
    using ValueTypes   = MakeTuple2_t<ValueTypeGenerator, FeSpaces, Dimensions>;
    using DOFMatrices  = MakeTuple2_t<DOFMatrixGenerator, FeSpaces, Dimensions>;
    using MultiMatrix  = ExpandTuple_t<Dune::MultiTypeBlockMatrix, MultiVectors>;
127

128
    static constexpr std::size_t nComponents = std::tuple_size<FeSpaces>::value;
129

130
131
    AMDIS_STATIC_ASSERT( nComponents > 0 );
    AMDIS_STATIC_ASSERT( nComponents == std::tuple_size<Dimensions>::value );
132

133
134
135
136
    /// Constructor that constructs new base-vectors
    SystemMatrix(FeSpaces const& feSpaces)
      : feSpaces(feSpaces)
      , matrix{}
137
      , dofmatrices(Impl::buildDOFMatrices<DOFMatrices>(feSpaces, feSpaces, matrix))
138
    {}
139

140
141

    /// Return the number of vector entries
142
    static constexpr std::size_t N()
143
144
145
    {
      return std::tuple_size<FeSpaces>::value;
    }
146

147
    static constexpr std::size_t M()
148
149
150
    {
      return std::tuple_size<FeSpaces>::value;
    }
151

152
    /// Return a shared pointer to the i'th underlaying base vector
153
154
    template <std::size_t R, std::size_t C>
    auto& getMatrix(const index_t<R> = {}, const index_t<C> = {})
155
156
157
158
    {
      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
      return std::get<C>(std::get<R>(matrix));
    }
159

160
    /// Return a shared pointer to the i'th underlaying base vector
161
162
    template <std::size_t R, std::size_t C>
    auto const& getMatrix(const index_t<R> = {}, const index_t<C> = {}) const
163
164
165
166
    {
      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
      return std::get<C>(std::get<R>(matrix));
    }
167

168
169
170
171
172
    /// Return the underlying multi vector
    MultiMatrix& getMatrix()
    {
        return matrix;
    }
173

174
175
176
177
    MultiMatrix const& getMatrix() const
    {
      return matrix;
    }
178

179
    /// Return the I'th finite element space
180
181
    template <std::size_t I = 0>
    auto& getRowFeSpace(const index_t<I> = {}) const
182
183
184
185
    {
      static_assert(I < N(), "Index out of range [0,N)");
      return std::get<I>(feSpaces);
    }
186

187
188
    template <std::size_t I = 0>
    auto const& getColFeSpace(const index_t<I> = {}) const
189
190
191
192
    {
      static_assert(I < M(), "Index out of range [0,M)");
      return std::get<I>(feSpaces);
    }
193

194
195
    template <std::size_t R = 0, std::size_t C = 0>
    auto& operator()(const index_t<R> = {}, const index_t<C> = {})
196
197
198
199
    {
      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
      return std::get<C>(std::get<R>(dofmatrices));
    }
200

201
    /// Create a DOFVector with references to the unerlaying data
202
203
    template <std::size_t R = 0, std::size_t C = 0>
    auto const& operator()(const index_t<R> = {}, const index_t<C> = {}) const
204
205
206
207
    {
      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
      return std::get<C>(std::get<R>(dofmatrices));
    }
208

209
210
211
212
213
  private:
    FeSpaces const& feSpaces;	///< a tuple of feSpaces
    MultiMatrix 	matrix;		///< a tuple of data-vectors
    DOFMatrices     dofmatrices;
  };
214
215


216
217
  // specialization for 1 component
  // -------------------------------------------------------------------------
218
219

  template <class FeSpace,
220
221
            class FieldType,
            int dim>
222
  class SystemMatrix<std::tuple<FeSpace>, FieldType, std::tuple<index_t<dim>>>
223
224
225
  {
    using ValueType     = Dune::FieldMatrix<FieldType, dim, dim>;
    using DOFMatrixType = DOFMatrix<FeSpace, FeSpace, ValueType>;
226

227
228
  public:
    using MultiMatrix   = typename DOFMatrixType::BaseMatrix;
229

230
    static constexpr std::size_t nComponent = 1;
231
232


233
234
235
236
237
    /// Constructor that constructs new base-vectors
    template <class FeSpaces>
    SystemMatrix(FeSpaces const& feSpaces)
      : dofmatrix(std::get<0>(feSpaces), std::get<0>(feSpaces), "")
    {}
238

239
240

    /// Return the number of vector entries
241
    static constexpr std::size_t N()
242
243
244
    {
      return 1;
    }
245

246
    static constexpr std::size_t M()
247
248
249
    {
      return 1;
    }
250

251
    /// Return a shared pointer to the i'th underlaying base vector
252
253
    template <std::size_t R = 0, std::size_t C = 0>
    MultiMatrix& getMatrix(const index_t<R> = {}, const index_t<C> = {})
254
255
256
257
    {
      static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
      return dofmatrix.getMatrix();
    }
258

259
    /// Return a shared pointer to the i'th underlaying base vector
260
261
    template <std::size_t R = 0, std::size_t C = 0>
    MultiMatrix const& getMatrix(const index_t<R> = {}, const index_t<C> = {}) const
262
263
264
265
    {
      static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
      return dofmatrix.getMatrix();
    }
266

267
268
269
270
271
    /// Return the underlying multi vector
    MultiMatrix const& getMatrix() const
    {
      return dofmatrix.getMatrix();
    }
272

273
274
275
276
    MultiMatrix& getMatrix()
    {
      return dofmatrix.getMatrix();
    }
277

278
    /// Return the I'th finite element space
279
280
    template <std::size_t I = 0>
    FeSpace const& getRowFeSpace(const index_t<I> = {}) const
281
282
283
284
    {
      static_assert(I == 0, "Index out of range [0,1)");
      return dofmatrix.getRowFeSpace();
    }
285

286
287
    template <std::size_t I = 0>
    FeSpace const& getColFeSpace(const index_t<I> = {}) const
288
289
290
291
    {
      static_assert(I == 0, "Index out of range [0,1)");
      return dofmatrix.getColFeSpace();
    }
292

293
294
    template <std::size_t R = 0, std::size_t C = 0>
    DOFMatrixType& operator()(const index_t<R> = {}, const index_t<C> = {})
295
296
297
298
    {
      static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
      return dofmatrix;
    }
299

300
    /// Create a DOFVector with references to the unerlaying data
301
302
    template <std::size_t R = 0, std::size_t C = 0>
    DOFMatrixType const& operator()(const index_t<R> = {}, const index_t<C> = {}) const
303
304
305
306
    {
      static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
      return dofmatrix;
    }
307
308

  private:
309
310
311
312
313
    DOFMatrixType  dofmatrix;
  };

#ifndef AMDIS_NO_EXTERN_SYSTEMMATRIX
    // explicit instantiation in SystemMatrix.cpp
314
315
    extern template class SystemMatrix<typename TestTraits<2,1>::FeSpaces>;
    extern template class SystemMatrix<typename TestTraits<2,1,2>::FeSpaces>;
316
#endif
317

318
} // end namespace AMDiS