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

3
4
#include <type_traits>

5
#include <dune/common/diagonalmatrix.hh>
6
7
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
8
#include <dune/common/typetraits.hh>
9

Praetorius, Simon's avatar
Praetorius, Simon committed
10
11
#include <amdis/common/TypeTraits.hpp>

12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
namespace std
{
  template <class T, int N>
  struct common_type<Dune::FieldVector<T,N>, T>
  {
    using type = T;
  };

  template <class T, int N, int M>
  struct common_type<Dune::FieldMatrix<T,N,M>, T>
  {
    using type = T;
  };
}

27
namespace Dune
28
{
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
  namespace MatVec
  {
    /// Traits to detect fixed size containers like FieldVector and FieldMatrix
    /// @{
    template <class T>
    struct IsMatrix : std::false_type {};

    template <class T, int M, int N>
    struct IsMatrix<FieldMatrix<T,M,N>> : std::true_type {};

    template <class T, int N>
    struct IsMatrix<DiagonalMatrix<T,N>> : std::true_type {};


    template <class T>
    struct IsVector : std::false_type {};

    template <class T, int N>
    struct IsVector<FieldVector<T,N>> : std::true_type {};


    template <class T>
    struct IsMatVec
      : std::integral_constant<bool, IsMatrix<T>::value || IsVector<T>::value> {};
    /// @}

    /// Convert the field_Type of container to type S
    /// @{
    template <class A, class S>
    struct MakeMatVec
    {
      using type = A;
    };

    template <class T, int M, int N, class S>
    struct MakeMatVec<FieldMatrix<T,M,N>,S>
    {
      using type = FieldMatrix<S,M,N>;
    };

    template <class T, int N, class S>
    struct MakeMatVec<DiagonalMatrix<T,N>,S>
    {
      using type = DiagonalMatrix<S,N>;
    };

    template <class T, int N, class S>
    struct MakeMatVec<FieldVector<T,N>,S>
    {
      using type = FieldVector<S,N>;
    };
    /// @}

    /// Convert pseudo-scalar to real scalar type
    /// @{
    template <class T>
    decltype(auto) simplify(T&& t)
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
87
      return FWD(t);
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
    }

    template <class T>
    T simplify(FieldVector<T,1> const& t)
    {
      return t[0];
    }

    template <class T>
    T simplify(FieldMatrix<T,1,1> const& t)
    {
      return t[0][0];
    }

    template <class T>
    T simplify(DiagonalMatrix<T,1> const& t)
    {
      return t.diagonal(0);
    }
    /// @}

    // returns -a
    template <class A>
    auto negate(A const& a);

    // returns a+b
    template <class A, class B>
    auto plus(A a, B const& b)
    {
      return a += b;
    }

    // returns a-b
    template <class A, class B>
    auto minus(A a, B const& b)
    {
      return a -= b;
    }

    // returns a*b
    template <class A, class B,
      std::enable_if_t<IsNumber<A>::value || IsNumber<B>::value, int> = 0>
    auto multiplies(A const& a, B const& b);

    template <class T, int N, class S>
    auto multiplies(FieldVector<T,N> const& a, FieldVector<S,N> const& b);

    template <class Mat, class Vec,
      std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value, int> = 0>
    auto multiplies(Mat const& mat, Vec const& vec);

    template <class Vec, class Mat,
      std::enable_if_t<IsVector<Vec>::value && IsMatrix<Mat>::value, int> = 0>
    auto multiplies(Vec const& vec, Mat const& mat);

    template <class T, int L, int M, int N, class S>
    auto multiplies(FieldMatrix<T,L,M> const& a, FieldMatrix<S,M,N> const& b);

    // return a/b
    template <class A, class B>
    auto divides(A a, B const& b)
    {
      return a /= b;
    }

  } // end namespace MatVec

  // some arithmetic operations with FieldVector and FieldMatrix

  template <class A,
    std::enable_if_t<MatVec::IsMatVec<A>::value, int> = 0>
  auto operator-(A const& a)
  {
    return MatVec::negate(MatVec::simplify(a));
  }
163

164
165
166
167
168
169
  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator+(A const& a, B const& b)
  {
    return MatVec::plus(MatVec::simplify(a), MatVec::simplify(b));
  }
170

171
172
173
174
175
176
  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator-(A const& a, B const& b)
  {
    return MatVec::minus(MatVec::simplify(a), MatVec::simplify(b));
  }
177

178
179
180
181
182
183
  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator*(A const& a, B const& b)
  {
    return MatVec::multiplies(MatVec::simplify(a), MatVec::simplify(b));
  }
184

185
186
187
188
189
190
  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator/(A const& a, B const& b)
  {
    return MatVec::divides(MatVec::simplify(a), MatVec::simplify(b));
  }
191
192
193
194
195

  // ----------------------------------------------------------------------------

  /// Cross-product a 2d-vector = orthogonal vector
  template <class T>
196
  FieldVector<T, 2> cross(FieldVector<T, 2> const& a);
197

198
  /// Cross-product of two 3d-vectors
199
  template <class T>
200
  FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b);
201
202
203

  /// Dot product (vec1^T * vec2)
  template <class T, class S, int N>
204
  auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2);
205

206
207
208
  template <class T, class S, int N>
  auto dot(FieldMatrix<T,1,N> const& vec1, FieldMatrix<S,1,N> const& vec2);

209
210
211
212
  // ----------------------------------------------------------------------------

  /// Sum of vector entires.
  template <class T, int N>
213
  T sum(FieldVector<T, N> const& x);
214

215
  template <class T, int N>
216
  T sum(FieldMatrix<T, 1, N> const& x);
217

218
219
220

  /// Dot-product with the vector itself
  template <class T, int N>
221
  auto unary_dot(FieldVector<T, N> const& x);
222

223
  template <class T, int N>
224
  auto unary_dot(FieldMatrix<T, 1, N> const& x);
225

226
227
  /// Maximum over all vector entries
  template <class T, int N>
228
  auto max(FieldVector<T, N> const& x);
229

230
  template <class T, int N>
231
  auto max(FieldMatrix<T, 1, N> const& x);
232

233
234
  /// Minimum over all vector entries
  template <class T, int N>
235
  auto min(FieldVector<T, N> const& x);
236

237
  template <class T, int N>
238
  auto min(FieldMatrix<T, 1, N> const& x);
239

240
241
  /// Maximum of the absolute values of vector entries
  template <class T, int N>
242
  auto abs_max(FieldVector<T, N> const& x);
243

244
  template <class T, int N>
245
  auto abs_max(FieldMatrix<T, 1, N> const& x);
246

247
248
  /// Minimum of the absolute values of vector entries
  template <class T, int N>
249
  auto abs_min(FieldVector<T, N> const& x);
250

251
  template <class T, int N>
252
  auto abs_min(FieldMatrix<T, 1, N> const& x);
253

254
255
256
257
258
259
  // ----------------------------------------------------------------------------

  /** \ingroup vector_norms
   *  \brief The 1-norm of a vector = sum_i |x_i|
   **/
  template <class T, int N>
260
  auto one_norm(FieldVector<T, N> const& x);
261

262
  template <class T, int N>
263
  auto one_norm(FieldMatrix<T, 1, N> const& x);
264

265
266
267
268
  /** \ingroup vector_norms
   *  \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 )
   **/
  template <class T, int N>
269
  auto two_norm(FieldVector<T, N> const& x);
270

271
  template <class T, int N>
272
  auto two_norm(FieldMatrix<T, 1, N> const& x);
273

274
275
276
277
  /** \ingroup vector_norms
   *  \brief The p-norm of a vector = ( sum_i |x_i|^p )^(1/p)
   **/
  template <int p, class T, int N>
278
  auto p_norm(FieldVector<T, N> const& x);
279

280
  template <int p, class T, int N>
281
  auto p_norm(FieldMatrix<T, 1, N> const& x);
282

283
284
285
286
  /** \ingroup vector_norms
   *  \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max
   **/
  template <class T, int N>
287
  auto infty_norm(FieldVector<T, N> const& x);
288

289
  template <class T, int N>
290
  auto infty_norm(FieldMatrix<T, 1, N> const& x);
291

292
293
294
295
  // ----------------------------------------------------------------------------

  /// The euklidean distance between two vectors = |lhs-rhs|_2
  template <class T, int N>
296
  T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs);
297
298
299
300
301

  // ----------------------------------------------------------------------------

  /// Outer product (vec1 * vec2^T)
  template <class T, class S, int N, int M, int K>
302
  auto outer(FieldMatrix<T,N,K> const& vec1, FieldMatrix<S,M,K> const& vec2);
303
304
305
306

  // ----------------------------------------------------------------------------

  template <class T>
307
  T det(FieldMatrix<T, 0, 0> const& /*mat*/);
308
309
310

  /// Determinant of a 1x1 matrix
  template <class T>
311
  T det(FieldMatrix<T, 1, 1> const& mat);
312
313
314

  /// Determinant of a 2x2 matrix
  template <class T>
315
  T det(FieldMatrix<T, 2, 2> const& mat);
316
317
318

  /// Determinant of a 3x3 matrix
  template <class T>
319
  T det(FieldMatrix<T, 3, 3> const& mat);
320
321
322

  /// Determinant of a NxN matrix
  template <class T,  int N>
323
324
  T det(FieldMatrix<T, N, N> const& mat);

325
326
327

  /// Return the inverse of the matrix `mat`
  template <class T, int N>
328
  auto inv(FieldMatrix<T, N, N> mat);
329
330
331

  /// Solve the linear system A*x = b
  template <class T, int N>
332
  void solve(FieldMatrix<T, N, N> const& A,  FieldVector<T, N>& x,  FieldVector<T, N> const& b);
333
334
335
336


  /// Gramian determinant = sqrt( det( DT^T * DF ) )
  template <class T, int N, int M>
337
  T gramian(FieldMatrix<T,N,M> const& DF);
338
339
340

  /// Gramian determinant, specialization for 1 column matrices
  template <class T, int M>
341
  T gramian(FieldMatrix<T, 1, M> const& DF);
342

343
344
345
  // ----------------------------------------------------------------------------
  // some arithmetic operations with FieldMatrix

346
  template <class T, int M, int N>
347
  FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A);
348

349
  // -----------------------------------------------------------------------------
350

351
  template <class T, int M, int N, int L>
352
  FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A,  FieldMatrix<T, N, L> const& B);
353
354

  template <class T, int M, int N, int L>
355
  FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> const& A,  FieldMatrix<T, N, L> const& B);
356
357

  template <class T, int M, int N, int L>
358
  FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A,  FieldMatrix<T, N, L> const& B, FieldMatrix<T,M,N>& C);
359

360
  template <class T, int M, int N>
361
  FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A,  DiagonalMatrix<T, N> const& B, FieldMatrix<T,M,N>& C);
362

363
  // -----------------------------------------------------------------------------
364
365
366
367
368
369
370

  template <class T, int N>
  T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i);

  template <class T, int M>
  T const& at(FieldMatrix<T,1,M> const& vec, std::size_t i);

371
  // necessary specialization to resolve ambiguous calls
372
373
374
  template <class T>
  T const& at(FieldMatrix<T,1,1> const& vec, std::size_t i);

375
376
377
  template <class T, int N>
  T const& at(FieldVector<T,N> const& vec, std::size_t i);

378
379
380
381
382
383
384
} // end namespace Dune

namespace AMDiS
{
  using Dune::FieldMatrix;
  using Dune::FieldVector;
}
385
386

#include "FieldMatVec.inc.hpp"