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

3
#include <algorithm>
4
#include <limits>
5

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

10
11
12
13
14
#include <amdis/common/Math.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ScalarTypes.hpp>
#include <amdis/operations/Arithmetic.hpp>
#include <amdis/operations/MaxMin.hpp>
15
16
17

namespace AMDiS
{
Praetorius, Simon's avatar
Praetorius, Simon committed
18
19
20
21
  using Dune::FieldVector;
  using Dune::FieldMatrix;

  // some arithmetic operations with FieldVector
22
23
24

  template <class T, int N, class S,
    REQUIRES(Concepts::Arithmetic<S>) >
Praetorius, Simon's avatar
Praetorius, Simon committed
25
  FieldVector<T,N> operator*(FieldVector<T,N> v, S factor)
26
27
28
29
30
31
  {
    return v *= factor;
  }

  template <class S, class T, int N,
    REQUIRES(Concepts::Arithmetic<S>) >
Praetorius, Simon's avatar
Praetorius, Simon committed
32
  FieldVector<T,N> operator*(S factor, FieldVector<T,N> v)
33
34
35
36
37
38
  {
    return v *= factor;
  }

  template <class T, int N, class S,
    REQUIRES(Concepts::Arithmetic<S>) >
Praetorius, Simon's avatar
Praetorius, Simon committed
39
  FieldVector<T,N> operator/(FieldVector<T,N> v, S factor)
40
41
42
43
  {
    return v /= factor;
  }

44
45
  template <class T>
  FieldVector<T,1> operator*(FieldVector<T,1> const& v, FieldVector<T,1> const& w)
46
  {
47
    return {v[0] * w[0]};
48
49
  }

50
51
  template <class T, int N>
  FieldVector<T,N> operator*(FieldVector<T,1> const& factor, FieldVector<T,N> v)
52
  {
53
    return v *= factor[0];
54
55
  }

56
57
  template <class T, int N>
  FieldVector<T,N> operator*(FieldVector<T,N> v, FieldVector<T,1> const& factor)
58
  {
59
    return v *= factor[0];
60
61
  }

62
63
64
65
  // ----------------------------------------------------------------------------

  /// Cross-product a 2d-vector = orthogonal vector
  template <class T>
Praetorius, Simon's avatar
Praetorius, Simon committed
66
  FieldVector<T, 2> cross(FieldVector<T, 2> const& a)
67
68
69
70
71
72
  {
    return {{ a[1], -a[0] }};
  }

  /// Cross-product of two vectors (in 3d only)
  template <class T>
Praetorius, Simon's avatar
Praetorius, Simon committed
73
  FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b)
74
75
76
77
78
79
80
81
  {
    return {{ a[1]*b[2] - a[2]*b[1],
              a[2]*b[0] - a[0]*b[2],
              a[0]*b[1] - a[1]*b[0] }};
  }

  /// Dot product (vec1^T * vec2)
  template <class T, class S, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
82
  auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2)
83
84
85
86
  {
    return vec1.dot(vec2);
  }

87
88
89
90
91
92
93
94
  template <class T, int N, int M,
    REQUIRES( N!=1 && M!=1 )>
  auto operator*(FieldVector<T,N> v, FieldVector<T,M> const& w)
  {
    static_assert(M == N, "Requires vectors of the same type!");
    return v.dot(w);
  }

95
96
97
98
99
  // ----------------------------------------------------------------------------

  namespace Impl
  {
    template <class T, int N, class Operation>
100
    T accumulate(FieldVector<T, N> const& x, T init, Operation op)
101
102
    {
      for (int i = 0; i < N; ++i)
103
104
        init = op(init, x[i]);
      return init;
105
106
    }

107
    template <class T, int N, class Operation>
108
    T accumulate(FieldMatrix<T, 1, N> const& x, T init, Operation op)
109
110
    {
      for (int i = 0; i < N; ++i)
111
112
        init = op(init, x[0][i]);
      return init;
113
114
    }

115
116
117
118
  } // end namespace Impl

  /// Sum of vector entires.
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
119
  T sum(FieldVector<T, N> const& x)
120
  {
121
    return Impl::accumulate(x, T(0), Operation::Plus{});
122
123
  }

124
125
126
  template <class T, int N>
  T sum(FieldMatrix<T, 1, N> const& x)
  {
127
    return Impl::accumulate(x, T(0), Operation::Plus{});
128
129
  }

130
131
132

  /// Dot-product with the vector itself
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
133
  auto unary_dot(FieldVector<T, N> const& x)
134
  {
135
    auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); };
136
    return Impl::accumulate(x, T(0), op);
137
138
  }

139
140
141
142
  template <class T, int N>
  auto unary_dot(FieldMatrix<T, 1, N> const& x)
  {
    auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); };
143
    return Impl::accumulate(x, T(0), op);
144
145
  }

146
147
  /// Maximum over all vector entries
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
148
  auto max(FieldVector<T, N> const& x)
149
  {
150
    return Impl::accumulate(x, std::numeric_limits<T>::lowest(), Operation::Max{});
151
152
  }

153
154
155
  template <class T, int N>
  auto max(FieldMatrix<T, 1, N> const& x)
  {
156
    return Impl::accumulate(x, std::numeric_limits<T>::lowest(), Operation::Max{});
157
158
  }

159
160
  /// Minimum over all vector entries
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
161
  auto min(FieldVector<T, N> const& x)
162
  {
163
    return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::Min{});
164
165
  }

166
167
168
  template <class T, int N>
  auto min(FieldMatrix<T, 1, N> const& x)
  {
169
    return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::Min{});
170
171
  }

172
173
  /// Maximum of the absolute values of vector entries
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
174
  auto abs_max(FieldVector<T, N> const& x)
175
  {
176
    return Impl::accumulate(x, T(0), Operation::AbsMax{});
177
178
  }

179
180
181
  template <class T, int N>
  auto abs_max(FieldMatrix<T, 1, N> const& x)
  {
182
    return Impl::accumulate(x, T(0), Operation::AbsMax{});
183
184
  }

185
186
  /// Minimum of the absolute values of vector entries
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
187
  auto abs_min(FieldVector<T, N> const& x)
188
  {
189
    return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::AbsMin{});
190
191
  }

192
193
194
  template <class T, int N>
  auto abs_min(FieldMatrix<T, 1, N> const& x)
  {
195
    return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::AbsMin{});
196
197
  }

198
199
200
201
202
203
  // ----------------------------------------------------------------------------

  /** \ingroup vector_norms
   *  \brief The 1-norm of a vector = sum_i |x_i|
   **/
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
204
  auto one_norm(FieldVector<T, N> const& x)
205
206
  {
    auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
207
    return Impl::accumulate(x, T(0), op);
208
209
  }

210
211
212
213
  template <class T, int N>
  auto one_norm(FieldMatrix<T, 1, N> const& x)
  {
    auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
214
    return Impl::accumulate(x, T(0), op);
215
216
  }

217
218
219
220
  /** \ingroup vector_norms
   *  \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 )
   **/
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
221
  auto two_norm(FieldVector<T, N> const& x)
222
223
224
225
  {
    return std::sqrt(unary_dot(x));
  }

226
227
228
229
230
231
  template <class T, int N>
  auto two_norm(FieldMatrix<T, 1, N> const& x)
  {
    return std::sqrt(unary_dot(x));
  }

232
233
234
235
  /** \ingroup vector_norms
   *  \brief The p-norm of a vector = ( sum_i |x_i|^p )^(1/p)
   **/
  template <int p, class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
236
  auto p_norm(FieldVector<T, N> const& x)
237
238
  {
    auto op = [](auto const& a, auto const& b) { return a + Math::pow<p>(std::abs(b)); };
239
    return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
240
241
  }

242
243
244
245
  template <int p, class T, int N>
  auto p_norm(FieldMatrix<T, 1, N> const& x)
  {
    auto op = [](auto const& a, auto const& b) { return a + Math::pow<p>(std::abs(b)); };
246
    return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
247
248
  }

249
250
251
252
  /** \ingroup vector_norms
   *  \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max
   **/
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
253
  auto infty_norm(FieldVector<T, N> const& x)
254
255
256
257
  {
    return abs_max(x);
  }

258
259
260
261
262
263
  template <class T, int N>
  auto infty_norm(FieldMatrix<T, 1, N> const& x)
  {
    return abs_max(x);
  }

264
265
266
267
  // ----------------------------------------------------------------------------

  /// The euklidean distance between two vectors = |lhs-rhs|_2
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
268
  T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs)
269
270
271
  {
    T result = 0;
    for (int i = 0; i < N; ++i)
272
      result += Math::sqr(lhs[i] - rhs[i]);
273
274
275
276
277
278
279
    return std::sqrt(result);
  }

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

  /// Outer product (vec1 * vec2^T)
  template <class T, class S, int N, int M, int K>
Praetorius, Simon's avatar
Praetorius, Simon committed
280
  auto outer(FieldMatrix<T,N,K> const& vec1, FieldMatrix<S,M,K> const& vec2)
281
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
282
    using result_type = FieldMatrix<decltype( std::declval<T>() * std::declval<S>() ), N, M>;
283
284
285
286
287
288
289
290
291
292
    result_type mat;
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < M; ++j)
        mat[i][j] = vec1[i].dot(vec2[j]);
    return mat;
  }

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

  template <class T>
Praetorius, Simon's avatar
Praetorius, Simon committed
293
  T det(FieldMatrix<T, 0, 0> const& /*mat*/)
294
295
296
297
298
299
  {
    return 0;
  }

  /// Determinant of a 1x1 matrix
  template <class T>
Praetorius, Simon's avatar
Praetorius, Simon committed
300
  T det(FieldMatrix<T, 1, 1> const& mat)
301
302
303
304
305
306
  {
    return mat[0][0];
  }

  /// Determinant of a 2x2 matrix
  template <class T>
Praetorius, Simon's avatar
Praetorius, Simon committed
307
  T det(FieldMatrix<T, 2, 2> const& mat)
308
309
310
311
312
313
  {
    return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
  }

  /// Determinant of a 3x3 matrix
  template <class T>
Praetorius, Simon's avatar
Praetorius, Simon committed
314
  T det(FieldMatrix<T, 3, 3> const& mat)
315
316
317
318
319
320
321
  {
    return mat[0][0]*mat[1][1]*mat[2][2] + mat[0][1]*mat[1][2]*mat[2][0] + mat[0][2]*mat[1][0]*mat[2][1]
        - (mat[0][2]*mat[1][1]*mat[2][0] + mat[0][1]*mat[1][0]*mat[2][2] + mat[0][0]*mat[1][2]*mat[2][1]);
  }

  /// Determinant of a NxN matrix
  template <class T,  int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
322
  T det(FieldMatrix<T, N, N> const& mat)
323
324
325
326
327
328
  {
    return mat.determinant();
  }

  /// Return the inverse of the matrix `mat`
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
329
  auto inv(FieldMatrix<T, N, N> mat)
330
331
332
333
334
335
336
  {
    mat.invert();
    return mat;
  }

  /// Solve the linear system A*x = b
  template <class T, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
337
  void solve(FieldMatrix<T, N, N> const& A,  FieldVector<T, N>& x,  FieldVector<T, N> const& b)
338
339
340
341
342
343
344
  {
    A.solve(x, b);
  }


  /// Gramian determinant = sqrt( det( DT^T * DF ) )
  template <class T, int N, int M>
Praetorius, Simon's avatar
Praetorius, Simon committed
345
  T gramian(FieldMatrix<T,N,M> const& DF)
346
347
348
349
350
351
352
  {
    using std::sqrt;
    return sqrt( det(outer(DF, DF)) );
  }

  /// Gramian determinant, specialization for 1 column matrices
  template <class T, int M>
Praetorius, Simon's avatar
Praetorius, Simon committed
353
  T gramian(FieldMatrix<T, 1, M> const& DF)
354
355
356
357
358
  {
    using std::sqrt;
    return sqrt(dot(DF[0], DF[0]));
  }

359
360
361
  // ----------------------------------------------------------------------------
  // some arithmetic operations with FieldMatrix

362
  template <class T, int M, int N>
Praetorius, Simon's avatar
Praetorius, Simon committed
363
  FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A)
364
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
365
    FieldMatrix<T,N,M> At;
366
367
368
369
370
371
372
373
    for (int i = 0; i < M; ++i)
      for (int j = 0; j < N; ++j)
        At[j][i] = A[i][j];

    return At;
  }


374
375
376
  template <class T, int M, int N, class S,
    REQUIRES(Concepts::Arithmetic<S>) >
  FieldMatrix<T,M,N> operator*(S scalar, FieldMatrix<T, M, N> A)
377
378
379
380
  {
    return A *= scalar;
  }

381
382
383
  template <class T, int M, int N, class S,
    REQUIRES(Concepts::Arithmetic<S>) >
  FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, S scalar)
384
385
386
387
  {
    return A *= scalar;
  }

388
389
390
  template <class T, int M, int N, class S,
    REQUIRES(Concepts::Arithmetic<S>) >
  FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, S scalar)
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
  {
    return A /= scalar;
  }

  template <class T, int M, int N>
  FieldMatrix<T,M,N> operator+(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B)
  {
    return A += B;
  }

  template <class T, int M, int N>
  FieldMatrix<T,M,N> operator-(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B)
  {
    return A -= B;
  }


  template <class T, int N, int M>
  FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const& vec)
  {
    return Dune::FMatrixHelp::mult(mat, vec);
  }



416
  template <class T, int M, int N, int L>
Praetorius, Simon's avatar
Praetorius, Simon committed
417
  FieldMatrix<T,M,N> multiplies(FieldMatrix<T, M, L> const& A,  FieldMatrix<T, L, N> const& B)
418
419
420
421
422
423
424
  {
    return A.rightmultiplyany(B);
  }



  template <class T, int M, int N, int L>
Praetorius, Simon's avatar
Praetorius, Simon committed
425
  FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A,  FieldMatrix<T, N, L> const& B)
426
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
427
    FieldMatrix<T,M,N> C;
428
429
430
431
432
433
434
435
436
437
438
439

    for (int m = 0; m < M; ++m) {
      for (int n = 0; n < N; ++n) {
        C[m][n] = 0;
        for (int l = 0; l < L; ++l)
          C[m][n] += A[l][m] * B[n][l];
      }
    }
    return C;
  }

  template <class T, int M, int N, int L>
Praetorius, Simon's avatar
Praetorius, Simon committed
440
  FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> const& A,  FieldMatrix<T, N, L> const& B)
441
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
442
    FieldMatrix<T,M,N> C;
443
444
445
446
447
448
449
450
451
452
453
454

    for (int m = 0; m < M; ++m) {
      for (int n = 0; n < N; ++n) {
        C[m][n] = 0;
        for (int l = 0; l < L; ++l)
          C[m][n] += A[m][l] * B[n][l];
      }
    }
    return C;
  }

  template <class T, int M, int N, int L>
Praetorius, Simon's avatar
Praetorius, Simon committed
455
  FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A,  FieldMatrix<T, N, L> const& B, FieldMatrix<T,M,N>& C)
456
457
458
459
460
461
462
463
464
465
466
  {
    for (int m = 0; m < M; ++m) {
      for (int n = 0; n < N; ++n) {
        C[m][n] = 0;
        for (int l = 0; l < L; ++l)
          C[m][n] += A[m][l] * B[n][l];
      }
    }
    return C;
  }

467
  template <class T, int M, int N>
468
  FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A,  Dune::DiagonalMatrix<T, N> const& B, FieldMatrix<T,M,N>& C)
469
470
471
472
473
474
475
476
477
  {
    for (int m = 0; m < M; ++m) {
      for (int n = 0; n < N; ++n) {
        C[m][n] = A[m][n] * B.diagonal(n);
      }
    }
    return C;
  }

478
} // end namespace AMDiS