FieldMatVec.hpp 11.7 KB
Newer Older
1
2
#pragma once

3
4
#include <type_traits>

5
#include <dune/common/diagonalmatrix.hh>
6
7
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
8
#include <dune/common/typetraits.hh>
9

Praetorius, Simon's avatar
Praetorius, Simon committed
10
11
#include <amdis/common/TypeTraits.hpp>

12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
namespace std
{
  template <class T, int N>
  struct common_type<Dune::FieldVector<T,N>, T>
  {
    using type = T;
  };

  template <class T, int N, int M>
  struct common_type<Dune::FieldMatrix<T,N,M>, T>
  {
    using type = T;
  };
}

27
namespace Dune
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
  namespace MatVec
  {
    /// Traits to detect fixed size containers like FieldVector and FieldMatrix
    /// @{
    template <class T>
    struct IsMatrix : std::false_type {};

    template <class T, int M, int N>
    struct IsMatrix<FieldMatrix<T,M,N>> : std::true_type {};

    template <class T, int N>
    struct IsMatrix<DiagonalMatrix<T,N>> : std::true_type {};


    template <class T>
    struct IsVector : std::false_type {};

    template <class T, int N>
    struct IsVector<FieldVector<T,N>> : std::true_type {};


    template <class T>
    struct IsMatVec
      : std::integral_constant<bool, IsMatrix<T>::value || IsVector<T>::value> {};
    /// @}

    /// Convert the field_Type of container to type S
    /// @{
    template <class A, class S>
    struct MakeMatVec
    {
      using type = A;
    };

    template <class T, int M, int N, class S>
    struct MakeMatVec<FieldMatrix<T,M,N>,S>
    {
      using type = FieldMatrix<S,M,N>;
    };

    template <class T, int N, class S>
    struct MakeMatVec<DiagonalMatrix<T,N>,S>
    {
      using type = DiagonalMatrix<S,N>;
    };

    template <class T, int N, class S>
    struct MakeMatVec<FieldVector<T,N>,S>
    {
      using type = FieldVector<S,N>;
    };
    /// @}

    /// Convert pseudo-scalar to real scalar type
    /// @{
    template <class T>
85
    decltype(auto) as_scalar(T&& t)
86
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
87
      return FWD(t);
88
89
90
    }

    template <class T>
91
    T as_scalar(FieldVector<T,1> const& t)
92
93
94
95
96
    {
      return t[0];
    }

    template <class T>
97
    T as_scalar(FieldMatrix<T,1,1> const& t)
98
99
100
101
102
    {
      return t[0][0];
    }

    template <class T>
103
    T as_scalar(DiagonalMatrix<T,1> const& t)
104
105
106
107
108
    {
      return t.diagonal(0);
    }
    /// @}

109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
    /// Convert pseudo-vector to real vector type
    /// @{
    template <class T>
    decltype(auto) as_vector(T&& t)
    {
      return FWD(t);
    }

    template <class T, int N>
    FieldVector<T,N> const& as_vector(FieldMatrix<T,1,N> const& t)
    {
      return t[0];
    }

    template <class T, int N>
    FieldVector<T,N>& as_vector(FieldMatrix<T,1,N>& t)
    {
      return t[0];
    }
    /// @}

130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147

    /// Convert pseudo-matrix to real matrix type with proper operator[][]
    /// @{
    template <class T>
    decltype(auto) as_matrix(T&& t)
    {
      return FWD(t);
    }

    template <class Mat>
    class MatrixView;

    template <class T, int N>
    MatrixView<DiagonalMatrix<T,N>> as_matrix(DiagonalMatrix<T,N> const& mat)
    {
      return {mat};
    }

148
149
150
151
152
153
    // returns -a
    template <class A>
    auto negate(A const& a);

    // returns a+b
    template <class A, class B>
154
    auto plus(A const& a, B const& b);
155
156
157

    // returns a-b
    template <class A, class B>
158
    auto minus(A const& a, B const& b);
159
160
161

    // returns a*b
    template <class A, class B,
162
163
164
165
166
      std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value, int> = 0>
    auto multiplies(A const& a, B const& b);

    template <class A, class B,
      std::enable_if_t<IsNumber<A>::value != IsNumber<B>::value, int> = 0>
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
    auto multiplies(A const& a, B const& b);

    template <class T, int N, class S>
    auto multiplies(FieldVector<T,N> const& a, FieldVector<S,N> const& b);

    template <class Mat, class Vec,
      std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value, int> = 0>
    auto multiplies(Mat const& mat, Vec const& vec);

    template <class Vec, class Mat,
      std::enable_if_t<IsVector<Vec>::value && IsMatrix<Mat>::value, int> = 0>
    auto multiplies(Vec const& vec, Mat const& mat);

    template <class T, int L, int M, int N, class S>
    auto multiplies(FieldMatrix<T,L,M> const& a, FieldMatrix<S,M,N> const& b);

    // return a/b
    template <class A, class B>
    auto divides(A a, B const& b)
    {
      return a /= b;
    }

  } // end namespace MatVec

  // some arithmetic operations with FieldVector and FieldMatrix

  template <class A,
    std::enable_if_t<MatVec::IsMatVec<A>::value, int> = 0>
  auto operator-(A const& a)
  {
198
    return MatVec::negate(MatVec::as_scalar(a));
199
  }
200

201
202
203
204
  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator+(A const& a, B const& b)
  {
205
    return MatVec::plus(MatVec::as_scalar(a), MatVec::as_scalar(b));
206
  }
207

208
209
210
211
  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator-(A const& a, B const& b)
  {
212
    return MatVec::minus(MatVec::as_scalar(a), MatVec::as_scalar(b));
213
  }
214

215
216
217
218
  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator*(A const& a, B const& b)
  {
219
    return MatVec::multiplies(MatVec::as_scalar(a), MatVec::as_scalar(b));
220
  }
221

222
223
224
225
  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator/(A const& a, B const& b)
  {
226
    return MatVec::divides(MatVec::as_scalar(a), MatVec::as_scalar(b));
227
  }
228
229
230
231
232

  // ----------------------------------------------------------------------------

  /// Cross-product a 2d-vector = orthogonal vector
  template <class T>
233
  FieldVector<T, 2> cross(FieldVector<T, 2> const& a);
234

235
  /// Cross-product of two 3d-vectors
236
  template <class T>
237
  FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b);
238
239
240

  /// Dot product (vec1^T * vec2)
  template <class T, class S, int N>
241
  auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2);
242

243
244
245
  template <class T, class S, int N>
  auto dot(FieldMatrix<T,1,N> const& vec1, FieldMatrix<S,1,N> const& vec2);

246
247
248
249
  // ----------------------------------------------------------------------------

  /// Sum of vector entires.
  template <class T, int N>
250
  T sum(FieldVector<T, N> const& x);
251

252
  template <class T, int N>
253
  T sum(FieldMatrix<T, 1, N> const& x);
254

255
256
257

  /// Dot-product with the vector itself
  template <class T, int N>
258
  auto unary_dot(FieldVector<T, N> const& x);
259

260
  template <class T, int N>
261
  auto unary_dot(FieldMatrix<T, 1, N> const& x);
262

263
264
  /// Maximum over all vector entries
  template <class T, int N>
265
  auto max(FieldVector<T, N> const& x);
266

267
  template <class T, int N>
268
  auto max(FieldMatrix<T, 1, N> const& x);
269

270
271
  /// Minimum over all vector entries
  template <class T, int N>
272
  auto min(FieldVector<T, N> const& x);
273

274
  template <class T, int N>
275
  auto min(FieldMatrix<T, 1, N> const& x);
276

277
278
  /// Maximum of the absolute values of vector entries
  template <class T, int N>
279
  auto abs_max(FieldVector<T, N> const& x);
280

281
  template <class T, int N>
282
  auto abs_max(FieldMatrix<T, 1, N> const& x);
283

284
285
  /// Minimum of the absolute values of vector entries
  template <class T, int N>
286
  auto abs_min(FieldVector<T, N> const& x);
287

288
  template <class T, int N>
289
  auto abs_min(FieldMatrix<T, 1, N> const& x);
290

291
292
293
294
295
296
  // ----------------------------------------------------------------------------

  /** \ingroup vector_norms
   *  \brief The 1-norm of a vector = sum_i |x_i|
   **/
  template <class T, int N>
297
  auto one_norm(FieldVector<T, N> const& x);
298

299
  template <class T, int N>
300
  auto one_norm(FieldMatrix<T, 1, N> const& x);
301

302
303
304
305
  /** \ingroup vector_norms
   *  \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 )
   **/
  template <class T, int N>
306
  auto two_norm(FieldVector<T, N> const& x);
307

308
  template <class T, int N>
309
  auto two_norm(FieldMatrix<T, 1, N> const& x);
310

311
312
313
314
  /** \ingroup vector_norms
   *  \brief The p-norm of a vector = ( sum_i |x_i|^p )^(1/p)
   **/
  template <int p, class T, int N>
315
  auto p_norm(FieldVector<T, N> const& x);
316

317
  template <int p, class T, int N>
318
  auto p_norm(FieldMatrix<T, 1, N> const& x);
319

320
321
322
323
  /** \ingroup vector_norms
   *  \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max
   **/
  template <class T, int N>
324
  auto infty_norm(FieldVector<T, N> const& x);
325

326
  template <class T, int N>
327
  auto infty_norm(FieldMatrix<T, 1, N> const& x);
328

329
330
331
332
  // ----------------------------------------------------------------------------

  /// The euklidean distance between two vectors = |lhs-rhs|_2
  template <class T, int N>
333
  T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs);
334
335
336
337
338

  // ----------------------------------------------------------------------------

  /// Outer product (vec1 * vec2^T)
  template <class T, class S, int N, int M, int K>
339
  auto outer(FieldMatrix<T,N,K> const& vec1, FieldMatrix<S,M,K> const& vec2);
340
341
342
343

  // ----------------------------------------------------------------------------

  template <class T>
344
  T det(FieldMatrix<T, 0, 0> const& /*mat*/);
345
346
347

  /// Determinant of a 1x1 matrix
  template <class T>
348
  T det(FieldMatrix<T, 1, 1> const& mat);
349
350
351

  /// Determinant of a 2x2 matrix
  template <class T>
352
  T det(FieldMatrix<T, 2, 2> const& mat);
353
354
355

  /// Determinant of a 3x3 matrix
  template <class T>
356
  T det(FieldMatrix<T, 3, 3> const& mat);
357
358
359

  /// Determinant of a NxN matrix
  template <class T,  int N>
360
361
  T det(FieldMatrix<T, N, N> const& mat);

362
363
364

  /// Return the inverse of the matrix `mat`
  template <class T, int N>
365
  auto inv(FieldMatrix<T, N, N> mat);
366
367
368

  /// Solve the linear system A*x = b
  template <class T, int N>
369
  void solve(FieldMatrix<T, N, N> const& A,  FieldVector<T, N>& x,  FieldVector<T, N> const& b);
370
371
372
373


  /// Gramian determinant = sqrt( det( DT^T * DF ) )
  template <class T, int N, int M>
374
  T gramian(FieldMatrix<T,N,M> const& DF);
375
376
377

  /// Gramian determinant, specialization for 1 column matrices
  template <class T, int M>
378
  T gramian(FieldMatrix<T, 1, M> const& DF);
379

380
381
382
  // ----------------------------------------------------------------------------
  // some arithmetic operations with FieldMatrix

383
  template <class T, int M, int N>
384
  FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A);
385

386
387
388
389
390
391
  template <class T, int N>
  DiagonalMatrix<T,N> const& trans(DiagonalMatrix<T,N> const& A)
  {
    return A;
  }

392
  // -----------------------------------------------------------------------------
393

394
395
  template <class T1, class T2, int M, int N, int L>
  FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M> const& A,  FieldMatrix<T2, N, L> const& B);
396

397
398
  template <class T1, class T2, int M, int N, int L>
  FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B);
399

400
401
  template <class T1, class T2, class T3, int M, int N, int L>
  FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C);
402

403
404
405
406
407
  template <class T1, class T2, class T3, int M, int N>
  FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C);

  template <class T1, class T2, class T3, int N>
  FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C);
408

409
  // -----------------------------------------------------------------------------
410
411
412
413
414
415
416

  template <class T, int N>
  T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i);

  template <class T, int M>
  T const& at(FieldMatrix<T,1,M> const& vec, std::size_t i);

417
  // necessary specialization to resolve ambiguous calls
418
419
420
  template <class T>
  T const& at(FieldMatrix<T,1,1> const& vec, std::size_t i);

421
422
423
  template <class T, int N>
  T const& at(FieldVector<T,N> const& vec, std::size_t i);

424
425
426
427
428
429
430
} // end namespace Dune

namespace AMDiS
{
  using Dune::FieldMatrix;
  using Dune::FieldVector;
}
431
432

#include "FieldMatVec.inc.hpp"