AMDiS  0.1
The Adaptive Multi-Dimensional Simulation Toolbox
FieldMatVec.inc.hpp
1 #pragma once
2 
3 #include <algorithm>
4 #include <limits>
5 
6 #include <dune/functions/functionspacebases/flatvectorview.hh>
7 
8 #include <amdis/common/Math.hpp>
9 #include <amdis/operations/Arithmetic.hpp>
10 #include <amdis/operations/MaxMin.hpp>
11 
12 #ifndef DOXYGEN
13 
14 namespace Dune {
15 
16 namespace MatVec {
17 
18  template <class A>
19  auto negate(A const& a)
20  {
21  return multiplies(a, -1);
22  }
23 
24  // returns a+b
25  template <class A, class B>
26  auto plus(A const& a, B const& b)
27  {
28  using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
29  typename MakeMatVec<A,T>::type c{a};
30 
31  auto b_ = Dune::Functions::flatVectorView(b);
32  auto c_ = Dune::Functions::flatVectorView(c);
33  assert(int(b_.size()) == int(c_.size()));
34  for(int i = 0; i < int(b_.size()); ++i)
35  c_[i] += b_[i];
36 
37  return c;
38  }
39 
40  // returns a-b
41  template <class A, class B>
42  auto minus(A const& a, B const& b)
43  {
44  using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
45  typename MakeMatVec<A,T>::type c{a};
46 
47  auto b_ = Dune::Functions::flatVectorView(b);
48  auto c_ = Dune::Functions::flatVectorView(c);
49  assert(int(b_.size()) == int(c_.size()));
50  for(int i = 0; i < int(b_.size()); ++i)
51  c_[i] -= b_[i];
52 
53  return c;
54  }
55 
56  template <class A, class B,
57  std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value, int>>
58  auto multiplies(A const& a, B const& b)
59  {
60  return a * b;
61  }
62 
63  template <class A, class B,
64  std::enable_if_t<IsNumber<A>::value != IsNumber<B>::value, int>>
65  auto multiplies(A const& a, B const& b)
66  {
67  using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
68  if constexpr(IsNumber<A>::value) {
69  typename MakeMatVec<B,T>::type b_{b};
70  return b_ *= a;
71  } else {
72  typename MakeMatVec<A,T>::type a_{a};
73  return a_ *= b;
74  }
75  }
76 
77  template <class T, int N, class S>
78  auto multiplies(FieldVector<T,N> const& a, FieldVector<S,N> const& b)
79  {
80  return a.dot(b);
81  }
82 
83 
84  template <class Mat, class Vec,
85  std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value, int>>
86  auto multiplies(Mat const& mat, Vec const& vec)
87  {
88  static_assert(int(Mat::cols) == int(Vec::dimension), "");
89  using T = std::common_type_t<typename FieldTraits<Vec>::field_type, typename FieldTraits<Mat>::field_type>;
90  FieldVector<T,Mat::rows> y;
91  mat.mv(vec, y);
92  return y;
93  }
94 
95 
96  template <class Vec, class Mat,
97  std::enable_if_t<IsVector<Vec>::value && IsMatrix<Mat>::value, int>>
98  auto multiplies(Vec const& vec, Mat const& mat)
99  {
100  static_assert(int(Mat::rows) == int(Vec::dimension), "");
101  using T = std::common_type_t<typename FieldTraits<Vec>::field_type, typename FieldTraits<Mat>::field_type>;
102  FieldVector<T,Mat::cols> y;
103  mat.mtv(vec, y);
104  return y;
105  }
106 
107 
108  template <class T, int L, int M, int N, class S>
109  auto multiplies(FieldMatrix<T,L,M> const& a, FieldMatrix<S,M,N> const& b)
110  {
111  FieldMatrix<std::common_type_t<T,S>,L,N> C;
112 
113  for (int i = 0; i < L; ++i) {
114  for (int j = 0; j < N; ++j) {
115  C[i][j] = 0;
116  for (int k = 0; k < M; ++k)
117  C[i][j] += a[i][k]*b[k][j];
118  }
119  }
120  return C;
121  }
122 
123  template <class T, int SIZE>
124  class MatrixView<DiagonalMatrix<T,SIZE>>
125  {
126  using Matrix = DiagonalMatrix<T,SIZE>;
127  using size_type = typename Matrix::size_type;
128 
129  struct RowView
130  {
131  T operator[](size_type c) const
132  {
133  assert(0 <= c && c < mat_->M());
134  return c == r_ ? mat_->diagonal(r_) : T(0);
135  }
136 
137  Matrix const* mat_;
138  size_type r_;
139  };
140 
141  public:
142  MatrixView(Matrix const& mat)
143  : mat_(mat)
144  {}
145 
146  RowView operator[](size_type r) const
147  {
148  assert(0 <= r && r < mat_.N());
149  return {&mat_,r};
150  }
151 
152  size_type N() const
153  {
154  return mat_.N();
155  }
156 
157  size_type M() const
158  {
159  return mat_.M();
160  }
161 
162  private:
163  Matrix const& mat_;
164  };
165 
166 } // end namespace MatVec
167 
168 // ----------------------------------------------------------------------------
169 
171 template <class T>
172 FieldVector<T, 2> cross(FieldVector<T, 2> const& a)
173 {
174  return {{ a[1], -a[0] }};
175 }
176 
178 template <class T>
179 FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b)
180 {
181  return {{ a[1]*b[2] - a[2]*b[1],
182  a[2]*b[0] - a[0]*b[2],
183  a[0]*b[1] - a[1]*b[0] }};
184 }
185 
187 template <class T, class S, int N>
188 auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2)
189 {
190  return vec1.dot(vec2);
191 }
192 
193 template <class T, class S, int N>
194 auto dot(FieldMatrix<T,1,N> const& vec1, FieldMatrix<S,1,N> const& vec2)
195 {
196  return vec1[0].dot(vec2[0]);
197 }
198 
199 // ----------------------------------------------------------------------------
200 
201 namespace Impl
202 {
203  template <class T, int N, class Operation>
204  T accumulate(FieldVector<T, N> const& x, T init, Operation op)
205  {
206  for (int i = 0; i < N; ++i)
207  init = op(init, x[i]);
208  return init;
209  }
210 
211  template <class T, int N, class Operation>
212  T accumulate(FieldMatrix<T, 1, N> const& x, T init, Operation op)
213  {
214  for (int i = 0; i < N; ++i)
215  init = op(init, x[0][i]);
216  return init;
217  }
218 
219 } // end namespace Impl
220 
222 template <class T, int N>
223 T sum(FieldVector<T, N> const& x)
224 {
225  return Impl::accumulate(x, T(0), AMDiS::Operation::Plus{});
226 }
227 
228 template <class T, int N>
229 T sum(FieldMatrix<T, 1, N> const& x)
230 {
231  return Impl::accumulate(x, T(0), AMDiS::Operation::Plus{});
232 }
233 
234 
236 template <class T,
237  std::enable_if_t<Dune::IsNumber<T>::value, int> >
238 auto unary_dot(T const& x)
239 {
240  using std::abs;
241  return AMDiS::Math::sqr(abs(x));
242 }
243 
244 template <class T, int N>
245 auto unary_dot(FieldVector<T, N> const& x)
246 {
247  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::sqr(abs(b)); };
248  return Impl::accumulate(x, T(0), op);
249 }
250 
251 template <class T, int N>
252 auto unary_dot(FieldMatrix<T, 1, N> const& x)
253 {
254  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::sqr(abs(b)); };
255  return Impl::accumulate(x, T(0), op);
256 }
257 
259 template <class T, int N>
260 auto max(FieldVector<T, N> const& x)
261 {
262  return Impl::accumulate(x, std::numeric_limits<T>::lowest(), AMDiS::Operation::Max{});
263 }
264 
265 template <class T, int N>
266 auto max(FieldMatrix<T, 1, N> const& x)
267 {
268  return Impl::accumulate(x, std::numeric_limits<T>::lowest(), AMDiS::Operation::Max{});
269 }
270 
272 template <class T, int N>
273 auto min(FieldVector<T, N> const& x)
274 {
275  return Impl::accumulate(x, std::numeric_limits<T>::max(), AMDiS::Operation::Min{});
276 }
277 
278 template <class T, int N>
279 auto min(FieldMatrix<T, 1, N> const& x)
280 {
281  return Impl::accumulate(x, std::numeric_limits<T>::max(), AMDiS::Operation::Min{});
282 }
283 
285 template <class T, int N>
286 auto abs_max(FieldVector<T, N> const& x)
287 {
288  return Impl::accumulate(x, T(0), AMDiS::Operation::AbsMax{});
289 }
290 
291 template <class T, int N>
292 auto abs_max(FieldMatrix<T, 1, N> const& x)
293 {
294  return Impl::accumulate(x, T(0), AMDiS::Operation::AbsMax{});
295 }
296 
298 template <class T, int N>
299 auto abs_min(FieldVector<T, N> const& x)
300 {
301  return Impl::accumulate(x, std::numeric_limits<T>::max(), AMDiS::Operation::AbsMin{});
302 }
303 
304 template <class T, int N>
305 auto abs_min(FieldMatrix<T, 1, N> const& x)
306 {
307  return Impl::accumulate(x, std::numeric_limits<T>::max(), AMDiS::Operation::AbsMin{});
308 }
309 
310 // ----------------------------------------------------------------------------
311 
315 template <class T, int N>
316 auto one_norm(FieldVector<T, N> const& x)
317 {
318  auto op = [](auto const& a, auto const& b) { using std::abs; return a + abs(b); };
319  return Impl::accumulate(x, T(0), op);
320 }
321 
322 template <class T, int N>
323 auto one_norm(FieldMatrix<T, 1, N> const& x)
324 {
325  auto op = [](auto const& a, auto const& b) { using std::abs; return a + abs(b); };
326  return Impl::accumulate(x, T(0), op);
327 }
328 
332 template <class T,
333  std::enable_if_t<Dune::IsNumber<T>::value, int> >
334 auto two_norm(T const& x)
335 {
336  using std::abs;
337  return abs(x);
338 }
339 
340 template <class T, int N>
341 auto two_norm(FieldVector<T, N> const& x)
342 {
343  using std::sqrt;
344  return sqrt(unary_dot(x));
345 }
346 
347 template <class T, int N>
348 auto two_norm(FieldMatrix<T, 1, N> const& x)
349 {
350  using std::sqrt;
351  return sqrt(unary_dot(x));
352 }
353 
357 template <int p, class T, int N>
358 auto p_norm(FieldVector<T, N> const& x)
359 {
360  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::pow<p>(abs(b)); };
361  return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
362 }
363 
364 template <int p, class T, int N>
365 auto p_norm(FieldMatrix<T, 1, N> const& x)
366 {
367  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::pow<p>(abs(b)); };
368  return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
369 }
370 
374 template <class T, int N>
375 auto infty_norm(FieldVector<T, N> const& x)
376 {
377  return abs_max(x);
378 }
379 
380 template <class T, int N>
381 auto infty_norm(FieldMatrix<T, 1, N> const& x)
382 {
383  return abs_max(x);
384 }
385 
386 // ----------------------------------------------------------------------------
387 
389 template <class T,
390  std::enable_if_t<Dune::IsNumber<T>::value, int> >
391 T distance(T const& lhs, T const& rhs)
392 {
393  using std::abs;
394  return abs(lhs - rhs);
395 }
396 
397 template <class T, int N>
398 T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs)
399 {
400  using std::sqrt;
401  T result = 0;
402  for (int i = 0; i < N; ++i)
403  result += AMDiS::Math::sqr(lhs[i] - rhs[i]);
404  return sqrt(result);
405 }
406 
407 // ----------------------------------------------------------------------------
408 
410 template <class T, class S, int N, int M, int K>
411 auto outer(FieldMatrix<T,N,K> const& vec1, FieldMatrix<S,M,K> const& vec2)
412 {
413  using result_type = FieldMatrix<TYPEOF( std::declval<T>() * std::declval<S>() ), N, M>;
414  result_type mat;
415  for (int i = 0; i < N; ++i)
416  for (int j = 0; j < M; ++j)
417  mat[i][j] = vec1[i].dot(vec2[j]);
418  return mat;
419 }
420 
422 template <class T, class S, int N, int M>
423 auto outer(FieldVector<T,N> const& vec1, FieldVector<S,M> const& vec2)
424 {
425  using result_type = FieldMatrix<TYPEOF( std::declval<T>() * std::declval<S>() ), N, M>;
426  result_type mat;
427  for (int i = 0; i < N; ++i)
428  for (int j = 0; j < M; ++j)
429  mat[i][j] = vec1[i] * vec2[j];
430  return mat;
431 }
432 
433 // ----------------------------------------------------------------------------
434 
435 template <class T>
436 T det(FieldMatrix<T, 0, 0> const& /*mat*/)
437 {
438  return 0;
439 }
440 
442 template <class T>
443 T det(FieldMatrix<T, 1, 1> const& mat)
444 {
445  return mat[0][0];
446 }
447 
449 template <class T>
450 T det(FieldMatrix<T, 2, 2> const& mat)
451 {
452  return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
453 }
454 
456 template <class T>
457 T det(FieldMatrix<T, 3, 3> const& mat)
458 {
459  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]
460  - (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]);
461 }
462 
464 template <class T, int N>
465 T det(FieldMatrix<T, N, N> const& mat)
466 {
467  return mat.determinant();
468 }
469 
471 template <class T, int N>
472 auto inv(FieldMatrix<T, N, N> mat)
473 {
474  mat.invert();
475  return mat;
476 }
477 
479 template <class T, int N>
480 void solve(FieldMatrix<T, N, N> const& A, FieldVector<T, N>& x, FieldVector<T, N> const& b)
481 {
482  A.solve(x, b);
483 }
484 
485 
487 template <class T, int N, int M>
488 T gramian(FieldMatrix<T,N,M> const& DF)
489 {
490  using std::sqrt;
491  return sqrt( det(outer(DF, DF)) );
492 }
493 
495 template <class T, int M>
496 T gramian(FieldMatrix<T, 1, M> const& DF)
497 {
498  using std::sqrt;
499  return sqrt(dot(DF[0], DF[0]));
500 }
501 
502 // ----------------------------------------------------------------------------
503 // some arithmetic operations with FieldMatrix
504 
505 template <class T, int M, int N>
506 FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A)
507 {
508  FieldMatrix<T,N,M> At;
509  for (int i = 0; i < M; ++i)
510  for (int j = 0; j < N; ++j)
511  At[j][i] = A[i][j];
512 
513  return At;
514 }
515 
516 
517 template <class T1, class T2, int M, int N, int L>
518 FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M> const& A, FieldMatrix<T2, N, L> const& B)
519 {
520  FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
521 
522  for (int m = 0; m < M; ++m) {
523  for (int n = 0; n < N; ++n) {
524  C[m][n] = 0;
525  for (int l = 0; l < L; ++l)
526  C[m][n] += A[l][m] * B[n][l];
527  }
528  }
529  return C;
530 }
531 
532 template <class T1, class T2, int M, int N, int L>
533 FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B)
534 {
535  FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
536  return multiplies_ABt(A,B,C);
537 }
538 
539 template <class T1, class T2, class T3, int M, int N, int L>
540 FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C)
541 {
542  for (int m = 0; m < M; ++m) {
543  for (int n = 0; n < N; ++n) {
544  C[m][n] = 0;
545  for (int l = 0; l < L; ++l)
546  C[m][n] += A[m][l] * B[n][l];
547  }
548  }
549  return C;
550 }
551 
552 template <class T1, class T2, class T3, int N, int L>
553 FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& C)
554 {
555  for (int n = 0; n < N; ++n) {
556  C[n] = 0;
557  for (int l = 0; l < L; ++l)
558  C[n] += A[0][l] * B[n][l];
559  }
560  return C;
561 }
562 
563 template <class T1, class T2, class T3, int N, int L>
564 FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& C)
565 {
566  for (int n = 0; n < N; ++n) {
567  C[n] = 0;
568  for (int l = 0; l < L; ++l)
569  C[n] += A[l] * B[n][l];
570  }
571  return C;
572 }
573 
574 template <class T1, class T2, class T3, int N, int L>
575 FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,1,N>& C)
576 {
577  for (int n = 0; n < N; ++n) {
578  C[0][n] = 0;
579  for (int l = 0; l < L; ++l)
580  C[0][n] += A[l] * B[n][l];
581  }
582  return C;
583 }
584 
585 
586 
587 template <class T1, class T2, class T3, int M, int N>
588 FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C)
589 {
590  for (int m = 0; m < M; ++m) {
591  for (int n = 0; n < N; ++n) {
592  C[m][n] = A[m][n] * B.diagonal(n);
593  }
594  }
595  return C;
596 }
597 
598 template <class T1, class T2, class T3, int N>
599 FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C)
600 {
601  for (int n = 0; n < N; ++n) {
602  C[n] = A[0][n] * B.diagonal(n);
603  }
604  return C;
605 }
606 
607 template <class T1, class T2, class T3, int N>
608 FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C)
609 {
610  for (int n = 0; n < N; ++n) {
611  C[n] = A[n] * B.diagonal(n);
612  }
613  return C;
614 }
615 
616 template <class T1, class T2, class T3, int N>
617 FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,1,N>& C)
618 {
619  for (int n = 0; n < N; ++n) {
620  C[0][n] = A[n] * B.diagonal(n);
621  }
622  return C;
623 }
624 
625 
626 template <class T, int N>
627 T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i)
628 {
629  return vec[i][0];
630 }
631 
632 template <class T, int M>
633 T const& at(FieldMatrix<T,1,M> const& vec, std::size_t i)
634 {
635  return vec[0][i];
636 }
637 
638 template <class T>
639 T const& at(FieldMatrix<T,1,1> const& vec, std::size_t i)
640 {
641  return vec[0][i];
642 }
643 
644 template <class T, int N>
645 T const& at(FieldVector<T,N> const& vec, std::size_t i)
646 {
647  return vec[i];
648 }
649 
650 } // end namespace AMDiS
651 
652 #endif
Operation that represents min(A,B)
Definition: MaxMin.hpp:20
Definition: AdaptiveGrid.hpp:373
Operation that represents max(A,B)
Definition: MaxMin.hpp:8
Operation that represents max(|A|,|B|)
Definition: MaxMin.hpp:32
Functor that represents A+B.
Definition: Arithmetic.hpp:19
Operation that represents min(|A|,|B|)
Definition: MaxMin.hpp:45