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