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

3
4
#include <algorithm>

5
6
7
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>

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

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

  // some arithmetic operations with FieldVector
20
21
22

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

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

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

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

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

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

60
61
62
63
  // ----------------------------------------------------------------------------

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

  /// Cross-product of two vectors (in 3d only)
  template <class T>
Praetorius, Simon's avatar
Praetorius, Simon committed
71
  FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b)
72
73
74
75
76
77
78
79
  {
    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
80
  auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2)
81
82
83
84
  {
    return vec1.dot(vec2);
  }

85
86
87
88
89
90
91
92
  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);
  }

93
94
95
96
97
  // ----------------------------------------------------------------------------

  namespace Impl
  {
    template <class T, int N, class Operation>
Praetorius, Simon's avatar
Praetorius, Simon committed
98
    T accumulate(FieldVector<T, N> const& x, Operation op)
99
100
101
102
103
104
105
    {
      T result = 0;
      for (int i = 0; i < N; ++i)
        result = op(result, x[i]);
      return result;
    }

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

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, Operation::Plus{});
122
123
  }

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

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
137
138
    return Impl::accumulate(x, op);
  }

139
140
141
142
143
144
145
  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)); };
    return Impl::accumulate(x, op);
  }

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, Operation::Max{});
151
152
  }

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

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, Operation::Min{});
164
165
  }

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

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, Operation::AbsMax{});
177
178
  }

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

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, Operation::AbsMin{});
190
191
  }

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

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
207
208
209
  {
    auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
    return Impl::accumulate(x, op);
  }

210
211
212
213
214
215
216
  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); };
    return Impl::accumulate(x, op);
  }

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
239
240
241
  {
    auto op = [](auto const& a, auto const& b) { return a + Math::pow<p>(std::abs(b)); };
    return std::pow( Impl::accumulate(x, op), 1.0/p );
  }

242
243
244
245
246
247
248
  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)); };
    return std::pow( Impl::accumulate(x, op), 1.0/p );
  }

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
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
  template <class T, int M, int N>
  FieldMatrix<T,M,N> operator*(T scalar, FieldMatrix<T, M, N> A)
  {
    return A *= scalar;
  }

  template <class T, int M, int N>
  FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, T scalar)
  {
    return A *= scalar;
  }

  template <class T, int M, int N>
  FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, T scalar)
  {
    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);
  }



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



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

    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
437
  FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> const& A,  FieldMatrix<T, N, L> const& B)
438
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
439
    FieldMatrix<T,M,N> C;
440
441
442
443
444
445
446
447
448
449
450
451

    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
452
  FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A,  FieldMatrix<T, N, L> const& B, FieldMatrix<T,M,N>& C)
453
454
455
456
457
458
459
460
461
462
463
  {
    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;
  }

464
  template <class T, int M, int N>
465
  FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A,  Dune::DiagonalMatrix<T, N> const& B, FieldMatrix<T,M,N>& C)
466
467
468
469
470
471
472
473
474
  {
    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;
  }

475
} // end namespace AMDiS