AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
FieldMatVec.hpp
1#pragma once
2
3#include <type_traits>
4
5#include <dune/common/diagonalmatrix.hh>
6#include <dune/common/fmatrix.hh>
7#include <dune/common/fvector.hh>
8#include <dune/common/typetraits.hh>
9
10#include <amdis/common/TypeTraits.hpp>
11
12namespace std
13{
14 template <class T, int N>
15 struct common_type<Dune::FieldVector<T,N>, T>
16 {
17 using type = T;
18 };
19
20 template <class T, int N, int M>
21 struct common_type<Dune::FieldMatrix<T,N,M>, T>
22 {
23 using type = T;
24 };
25}
26
27namespace Dune
28{
29 namespace MatVec
30 {
33 template <class T>
34 struct IsMatrix : std::false_type {};
35
36 template <class T, int M, int N>
37 struct IsMatrix<FieldMatrix<T,M,N>> : std::true_type {};
38
39 template <class T, int N>
40 struct IsMatrix<DiagonalMatrix<T,N>> : std::true_type {};
41
42 template <class T>
43 struct IsScalarMatrix : std::false_type {};
44
45 template <class T>
46 struct IsScalarMatrix<FieldMatrix<T,1,1>> : std::true_type {};
47
48 template <class T>
49 struct IsScalarMatrix<DiagonalMatrix<T,1>> : std::true_type {};
50
51
52 template <class T>
53 struct IsVector : std::false_type {};
54
55 template <class T, int N>
56 struct IsVector<FieldVector<T,N>> : std::true_type {};
57
58 template <class T>
59 struct IsScalarVector : std::false_type {};
60
61 template <class T>
62 struct IsScalarVector<FieldVector<T,1>> : std::true_type {};
63
64
65
66 template <class T>
67 struct IsMatVec
68 : std::integral_constant<bool, IsMatrix<T>::value || IsVector<T>::value> {};
70
73 template <class A, class S>
74 struct MakeMatVec
75 {
76 using type = S;
77 };
78
79 template <class T, int M, int N, class S>
80 struct MakeMatVec<FieldMatrix<T,M,N>,S>
81 {
82 using type = FieldMatrix<S,M,N>;
83 };
84
85 template <class T, int N, class S>
86 struct MakeMatVec<DiagonalMatrix<T,N>,S>
87 {
88 using type = DiagonalMatrix<S,N>;
89 };
90
91 template <class T, int N, class S>
92 struct MakeMatVec<FieldVector<T,N>,S>
93 {
94 using type = FieldVector<S,N>;
95 };
97
100 template <class T>
101 decltype(auto) as_scalar(T&& t)
102 {
103 return FWD(t);
104 }
105
106 template <class T>
107 T as_scalar(FieldVector<T,1> const& t)
108 {
109 return t[0];
110 }
111
112 template <class T>
113 T as_scalar(FieldMatrix<T,1,1> const& t)
114 {
115 return t[0][0];
116 }
117
118 template <class T>
119 T as_scalar(DiagonalMatrix<T,1> const& t)
120 {
121 return t.diagonal(0);
122 }
124
127 template <class T>
128 decltype(auto) as_vector(T&& t)
129 {
130 return FWD(t);
131 }
132
133 template <class T, int N>
134 FieldVector<T,N> const& as_vector(FieldMatrix<T,1,N> const& t)
135 {
136 return t[0];
137 }
138
139 template <class T, int N>
140 FieldVector<T,N>& as_vector(FieldMatrix<T,1,N>& t)
141 {
142 return t[0];
143 }
145
146
149 template <class T>
150 decltype(auto) as_matrix(T&& t)
151 {
152 return FWD(t);
153 }
154
155 template <class Mat>
156 class MatrixView;
157
158 template <class T, int N>
159 MatrixView<DiagonalMatrix<T,N>> as_matrix(DiagonalMatrix<T,N> const& mat)
160 {
161 return {mat};
162 }
164
165 // returns -a
166 template <class A>
167 auto negate(A const& a);
168
169 // returns a+b
170 template <class A, class B>
171 auto plus(A const& a, B const& b);
172
173 // returns a-b
174 template <class A, class B>
175 auto minus(A const& a, B const& b);
176
177 // returns a*b
178 template <class A, class B,
179 std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value, int> = 0>
180 auto multiplies(A const& a, B const& b);
181
182 template <class A, class B,
183 std::enable_if_t<IsNumber<A>::value != IsNumber<B>::value, int> = 0>
184 auto multiplies(A const& a, B const& b);
185
186 template <class T, int N, class S>
187 auto multiplies(FieldVector<T,N> const& a, FieldVector<S,N> const& b);
188
189 template <class Mat, class Vec,
190 std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value, int> = 0>
191 auto multiplies(Mat const& mat, Vec const& vec);
192
193 template <class Vec, class Mat,
194 std::enable_if_t<IsVector<Vec>::value && IsMatrix<Mat>::value, int> = 0>
195 auto multiplies(Vec const& vec, Mat const& mat);
196
197 template <class T, int L, int M, int N, class S>
198 auto multiplies(FieldMatrix<T,L,M> const& a, FieldMatrix<S,M,N> const& b);
199
200 // return a/b
201 template <class A, class B>
202 auto divides(A a, B const& b)
203 {
204 return a /= b;
205 }
206
207 } // end namespace MatVec
208
209 // some arithmetic operations with FieldVector and FieldMatrix
210
211 template <class A,
212 std::enable_if_t<MatVec::IsMatVec<A>::value, int> = 0>
213 auto operator-(A const& a)
214 {
215 return MatVec::negate(MatVec::as_scalar(a));
216 }
217
218 template <class A, class B,
219 std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
220 auto operator+(A const& a, B const& b)
221 {
222 return MatVec::plus(MatVec::as_scalar(a), MatVec::as_scalar(b));
223 }
224
225 template <class A, class B,
226 std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
227 auto operator-(A const& a, B const& b)
228 {
229 return MatVec::minus(MatVec::as_scalar(a), MatVec::as_scalar(b));
230 }
231
232 template <class A, class B,
233 std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
234 auto operator*(A const& a, B const& b)
235 {
236 return MatVec::multiplies(MatVec::as_scalar(a), MatVec::as_scalar(b));
237 }
238
239 template <class A, class B,
240 std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
241 auto operator/(A const& a, B const& b)
242 {
243 return MatVec::divides(MatVec::as_scalar(a), MatVec::as_scalar(b));
244 }
245
246 // ----------------------------------------------------------------------------
247
249 template <class T>
250 FieldVector<T, 2> cross(FieldVector<T, 2> const& a);
251
253 template <class T>
254 FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b);
255
257 template <class T, class S, int N>
258 auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2);
259
260 template <class T, class S, int N>
261 auto dot(FieldMatrix<T,1,N> const& vec1, FieldMatrix<S,1,N> const& vec2);
262
263 // ----------------------------------------------------------------------------
264
266 template <class T, int N>
267 T sum(FieldVector<T, N> const& x);
268
269 template <class T, int N>
270 T sum(FieldMatrix<T, 1, N> const& x);
271
272
274 template <class T,
275 std::enable_if_t<Dune::IsNumber<T>::value, int> = 0>
276 auto unary_dot(T const& x);
277
278 template <class T, int N>
279 auto unary_dot(FieldVector<T, N> const& x);
280
281 template <class T, int N>
282 auto unary_dot(FieldMatrix<T, 1, N> const& x);
283
285 template <class T, int N>
286 auto max(FieldVector<T, N> const& x);
287
288 template <class T, int N>
289 auto max(FieldMatrix<T, 1, N> const& x);
290
292 template <class T, int N>
293 auto min(FieldVector<T, N> const& x);
294
295 template <class T, int N>
296 auto min(FieldMatrix<T, 1, N> const& x);
297
299 template <class T, int N>
300 auto abs_max(FieldVector<T, N> const& x);
301
302 template <class T, int N>
303 auto abs_max(FieldMatrix<T, 1, N> const& x);
304
306 template <class T, int N>
307 auto abs_min(FieldVector<T, N> const& x);
308
309 template <class T, int N>
310 auto abs_min(FieldMatrix<T, 1, N> const& x);
311
312 // ----------------------------------------------------------------------------
313
317 template <class T, int N>
318 auto one_norm(FieldVector<T, N> const& x);
319
320 template <class T, int N>
321 auto one_norm(FieldMatrix<T, 1, N> const& x);
322
326 template <class T,
327 std::enable_if_t<Dune::IsNumber<T>::value, int> = 0>
328 auto two_norm(T const& x);
329
330 template <class T, int N>
331 auto two_norm(FieldVector<T, N> const& x);
332
333 template <class T, int N>
334 auto two_norm(FieldMatrix<T, 1, N> const& x);
335
339 template <int p, class T, int N>
340 auto p_norm(FieldVector<T, N> const& x);
341
342 template <int p, class T, int N>
343 auto p_norm(FieldMatrix<T, 1, N> const& x);
344
348 template <class T, int N>
349 auto infty_norm(FieldVector<T, N> const& x);
350
351 template <class T, int N>
352 auto infty_norm(FieldMatrix<T, 1, N> const& x);
353
354 // ----------------------------------------------------------------------------
355
357 template <class T,
358 std::enable_if_t<Dune::IsNumber<T>::value, int> = 0>
359 T distance(T const& lhs, T const& rhs);
360
361 template <class T, int N>
362 T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs);
363
364 // ----------------------------------------------------------------------------
365
367 template <class T, class S, int N, int M, int K>
368 auto outer(FieldMatrix<T,N,K> const& vec1, FieldMatrix<S,M,K> const& vec2);
369
370 template <class T, class S, int N, int M>
371 auto outer(FieldVector<T,N> const& vec1, FieldVector<S,M> const& vec2);
372
373 // ----------------------------------------------------------------------------
374
375 template <class T>
376 T det(FieldMatrix<T, 0, 0> const& /*mat*/);
377
379 template <class T>
380 T det(FieldMatrix<T, 1, 1> const& mat);
381
383 template <class T>
384 T det(FieldMatrix<T, 2, 2> const& mat);
385
387 template <class T>
388 T det(FieldMatrix<T, 3, 3> const& mat);
389
391 template <class T, int N>
392 T det(FieldMatrix<T, N, N> const& mat);
393
394
396 template <class T, int N>
397 auto inv(FieldMatrix<T, N, N> mat);
398
400 template <class T, int N>
401 void solve(FieldMatrix<T, N, N> const& A, FieldVector<T, N>& x, FieldVector<T, N> const& b);
402
403
405 template <class T, int N, int M>
406 T gramian(FieldMatrix<T,N,M> const& DF);
407
409 template <class T, int M>
410 T gramian(FieldMatrix<T, 1, M> const& DF);
411
412 // ----------------------------------------------------------------------------
413 // some arithmetic operations with FieldMatrix
414
415 template <class T, int M, int N>
416 FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A);
417
418 template <class T, int N>
419 DiagonalMatrix<T,N> const& trans(DiagonalMatrix<T,N> const& A)
420 {
421 return A;
422 }
423
424 // -----------------------------------------------------------------------------
425
426 template <class T1, class T2, int M, int N, int L>
427 FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M> const& A, FieldMatrix<T2, N, L> const& B);
428
429 template <class T1, class T2, int M, int N, int L>
430 FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B);
431
432
433 template <class T1, class T2, class T3, int M, int N, int L>
434 FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C);
435
436 template <class T1, class T2, class T3, int N, int L>
437 FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& C);
438
439 template <class T1, class T2, class T3, int N, int L>
440 FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& C);
441
442 template <class T1, class T2, class T3, int N, int L>
443 FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,1,N>& C);
444
445
446 template <class T1, class T2, class T3, int M, int N>
447 FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C);
448
449 template <class T1, class T2, class T3, int N>
450 FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C);
451
452 template <class T1, class T2, class T3, int N>
453 FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C);
454
455 template <class T1, class T2, class T3, int N>
456 FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,1,N>& C);
457
458 // -----------------------------------------------------------------------------
459
460 template <class T, int N>
461 T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i);
462
463 template <class T, int M>
464 T const& at(FieldMatrix<T,1,M> const& vec, std::size_t i);
465
466 // necessary specialization to resolve ambiguous calls
467 template <class T>
468 T const& at(FieldMatrix<T,1,1> const& vec, std::size_t i);
469
470 template <class T, int N>
471 T const& at(FieldVector<T,N> const& vec, std::size_t i);
472
473} // end namespace Dune
474
475namespace AMDiS
476{
477 using Dune::FieldMatrix;
478 using Dune::FieldVector;
479}
480
481#include "FieldMatVec.inc.hpp"