6 #include <dune/functions/functionspacebases/flatvectorview.hh> 8 #include <amdis/common/Math.hpp> 9 #include <amdis/operations/Arithmetic.hpp> 10 #include <amdis/operations/MaxMin.hpp> 19 auto negate(A
const& a)
21 return multiplies(a, -1);
25 template <
class A,
class B>
26 auto plus(A
const& a, B
const& b)
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};
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)
41 template <
class A,
class B>
42 auto minus(A
const& a, B
const& b)
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};
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)
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)
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)
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};
72 typename MakeMatVec<A,T>::type a_{a};
77 template <
class T,
int N,
class S>
78 auto multiplies(FieldVector<T,N>
const& a, FieldVector<S,N>
const& b)
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)
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;
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)
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;
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)
111 FieldMatrix<std::common_type_t<T,S>,L,N> C;
113 for (
int i = 0; i < L; ++i) {
114 for (
int j = 0; j < N; ++j) {
116 for (
int k = 0; k < M; ++k)
117 C[i][j] += a[i][k]*b[k][j];
123 template <
class T,
int SIZE>
124 class MatrixView<DiagonalMatrix<T,SIZE>>
126 using Matrix = DiagonalMatrix<T,SIZE>;
127 using size_type =
typename Matrix::size_type;
131 T operator[](size_type c)
const 133 assert(0 <= c && c < mat_->M());
134 return c == r_ ? mat_->diagonal(r_) : T(0);
142 MatrixView(Matrix
const& mat)
146 RowView operator[](size_type r)
const 148 assert(0 <= r && r < mat_.N());
172 FieldVector<T, 2> cross(FieldVector<T, 2>
const& a)
174 return {{ a[1], -a[0] }};
179 FieldVector<T, 3> cross(FieldVector<T, 3>
const& a, FieldVector<T, 3>
const& b)
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] }};
187 template <
class T,
class S,
int N>
188 auto dot(FieldVector<T,N>
const& vec1, FieldVector<S,N>
const& vec2)
190 return vec1.dot(vec2);
193 template <
class T,
class S,
int N>
194 auto dot(FieldMatrix<T,1,N>
const& vec1, FieldMatrix<S,1,N>
const& vec2)
196 return vec1[0].dot(vec2[0]);
203 template <
class T,
int N,
class Operation>
204 T accumulate(FieldVector<T, N>
const& x, T init, Operation op)
206 for (
int i = 0; i < N; ++i)
207 init = op(init, x[i]);
211 template <
class T,
int N,
class Operation>
212 T accumulate(FieldMatrix<T, 1, N>
const& x, T init, Operation op)
214 for (
int i = 0; i < N; ++i)
215 init = op(init, x[0][i]);
222 template <
class T,
int N>
223 T sum(FieldVector<T, N>
const& x)
228 template <
class T,
int N>
229 T sum(FieldMatrix<T, 1, N>
const& x)
237 std::enable_if_t<Dune::IsNumber<T>::value,
int> >
238 auto unary_dot(T
const& x)
241 return AMDiS::Math::sqr(abs(x));
244 template <
class T,
int N>
245 auto unary_dot(FieldVector<T, N>
const& x)
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);
251 template <
class T,
int N>
252 auto unary_dot(FieldMatrix<T, 1, N>
const& x)
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);
259 template <
class T,
int N>
260 auto max(FieldVector<T, N>
const& x)
265 template <
class T,
int N>
266 auto max(FieldMatrix<T, 1, N>
const& x)
272 template <
class T,
int N>
273 auto min(FieldVector<T, N>
const& x)
278 template <
class T,
int N>
279 auto min(FieldMatrix<T, 1, N>
const& x)
285 template <
class T,
int N>
286 auto abs_max(FieldVector<T, N>
const& x)
291 template <
class T,
int N>
292 auto abs_max(FieldMatrix<T, 1, N>
const& x)
298 template <
class T,
int N>
299 auto abs_min(FieldVector<T, N>
const& x)
304 template <
class T,
int N>
305 auto abs_min(FieldMatrix<T, 1, N>
const& x)
315 template <
class T,
int N>
316 auto one_norm(FieldVector<T, N>
const& x)
318 auto op = [](
auto const& a,
auto const& b) {
using std::abs;
return a + abs(b); };
319 return Impl::accumulate(x, T(0), op);
322 template <
class T,
int N>
323 auto one_norm(FieldMatrix<T, 1, N>
const& x)
325 auto op = [](
auto const& a,
auto const& b) {
using std::abs;
return a + abs(b); };
326 return Impl::accumulate(x, T(0), op);
333 std::enable_if_t<Dune::IsNumber<T>::value,
int> >
334 auto two_norm(T
const& x)
340 template <
class T,
int N>
341 auto two_norm(FieldVector<T, N>
const& x)
344 return sqrt(unary_dot(x));
347 template <
class T,
int N>
348 auto two_norm(FieldMatrix<T, 1, N>
const& x)
351 return sqrt(unary_dot(x));
357 template <
int p,
class T,
int N>
358 auto p_norm(FieldVector<T, N>
const& x)
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 );
364 template <
int p,
class T,
int N>
365 auto p_norm(FieldMatrix<T, 1, N>
const& x)
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 );
374 template <
class T,
int N>
375 auto infty_norm(FieldVector<T, N>
const& x)
380 template <
class T,
int N>
381 auto infty_norm(FieldMatrix<T, 1, N>
const& x)
390 std::enable_if_t<Dune::IsNumber<T>::value,
int> >
391 T distance(T
const& lhs, T
const& rhs)
394 return abs(lhs - rhs);
397 template <
class T,
int N>
398 T distance(FieldVector<T, N>
const& lhs, FieldVector<T, N>
const& rhs)
402 for (
int i = 0; i < N; ++i)
403 result += AMDiS::Math::sqr(lhs[i] - rhs[i]);
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)
413 using result_type = FieldMatrix<TYPEOF( std::declval<T>() * std::declval<S>() ), N, M>;
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]);
422 template <
class T,
class S,
int N,
int M>
423 auto outer(FieldVector<T,N>
const& vec1, FieldVector<S,M>
const& vec2)
425 using result_type = FieldMatrix<TYPEOF( std::declval<T>() * std::declval<S>() ), N, M>;
427 for (
int i = 0; i < N; ++i)
428 for (
int j = 0; j < M; ++j)
429 mat[i][j] = vec1[i] * vec2[j];
436 T det(FieldMatrix<T, 0, 0>
const& )
443 T det(FieldMatrix<T, 1, 1>
const& mat)
450 T det(FieldMatrix<T, 2, 2>
const& mat)
452 return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
457 T det(FieldMatrix<T, 3, 3>
const& mat)
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]);
464 template <
class T,
int N>
465 T det(FieldMatrix<T, N, N>
const& mat)
467 return mat.determinant();
471 template <
class T,
int N>
472 auto inv(FieldMatrix<T, N, N> mat)
479 template <
class T,
int N>
480 void solve(FieldMatrix<T, N, N>
const& A, FieldVector<T, N>& x, FieldVector<T, N>
const& b)
487 template <
class T,
int N,
int M>
488 T gramian(FieldMatrix<T,N,M>
const& DF)
491 return sqrt( det(outer(DF, DF)) );
495 template <
class T,
int M>
496 T gramian(FieldMatrix<T, 1, M>
const& DF)
499 return sqrt(dot(DF[0], DF[0]));
505 template <
class T,
int M,
int N>
506 FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N>
const& A)
508 FieldMatrix<T,N,M> At;
509 for (
int i = 0; i < M; ++i)
510 for (
int j = 0; j < N; ++j)
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)
520 FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
522 for (
int m = 0; m < M; ++m) {
523 for (
int n = 0; n < N; ++n) {
525 for (
int l = 0; l < L; ++l)
526 C[m][n] += A[l][m] * B[n][l];
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)
535 FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
536 return multiplies_ABt(A,B,C);
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)
542 for (
int m = 0; m < M; ++m) {
543 for (
int n = 0; n < N; ++n) {
545 for (
int l = 0; l < L; ++l)
546 C[m][n] += A[m][l] * B[n][l];
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)
555 for (
int n = 0; n < N; ++n) {
557 for (
int l = 0; l < L; ++l)
558 C[n] += A[0][l] * B[n][l];
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)
566 for (
int n = 0; n < N; ++n) {
568 for (
int l = 0; l < L; ++l)
569 C[n] += A[l] * B[n][l];
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)
577 for (
int n = 0; n < N; ++n) {
579 for (
int l = 0; l < L; ++l)
580 C[0][n] += A[l] * B[n][l];
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)
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);
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)
601 for (
int n = 0; n < N; ++n) {
602 C[n] = A[0][n] * B.diagonal(n);
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)
610 for (
int n = 0; n < N; ++n) {
611 C[n] = A[n] * B.diagonal(n);
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)
619 for (
int n = 0; n < N; ++n) {
620 C[0][n] = A[n] * B.diagonal(n);
626 template <
class T,
int N>
627 T
const& at(FieldMatrix<T,N,1>
const& vec, std::size_t i)
632 template <
class T,
int M>
633 T
const& at(FieldMatrix<T,1,M>
const& vec, std::size_t i)
639 T
const& at(FieldMatrix<T,1,1>
const& vec, std::size_t i)
644 template <
class T,
int N>
645 T
const& at(FieldVector<T,N>
const& vec, std::size_t i)
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