FieldMatVec.hpp 11.3 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
  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>
85
    decltype(auto) as_scalar(T&& t)
86
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
87
      return FWD(t);
88
89
90
    }

    template <class T>
91
    T as_scalar(FieldVector<T,1> const& t)
92
93
94
95
96
    {
      return t[0];
    }

    template <class T>
97
    T as_scalar(FieldMatrix<T,1,1> const& t)
98
99
100
101
102
    {
      return t[0][0];
    }

    template <class T>
103
    T as_scalar(DiagonalMatrix<T,1> const& t)
104
105
106
107
108
    {
      return t.diagonal(0);
    }
    /// @}

109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
    /// Convert pseudo-vector to real vector type
    /// @{
    template <class T>
    decltype(auto) as_vector(T&& t)
    {
      return FWD(t);
    }

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

    template <class T, int N>
    FieldVector<T,N>& as_vector(FieldMatrix<T,1,N>& t)
    {
      return t[0];
    }
    /// @}

130
131
132
133
134
135
    // returns -a
    template <class A>
    auto negate(A const& a);

    // returns a+b
    template <class A, class B>
136
    auto plus(A const& a, B const& b);
137
138
139

    // returns a-b
    template <class A, class B>
140
    auto minus(A const& a, B const& b);
141
142
143

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

    template <class A, class B,
      std::enable_if_t<IsNumber<A>::value != IsNumber<B>::value, int> = 0>
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
    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)
  {
180
    return MatVec::negate(MatVec::as_scalar(a));
181
  }
182

183
184
185
186
  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)
  {
187
    return MatVec::plus(MatVec::as_scalar(a), MatVec::as_scalar(b));
188
  }
189

190
191
192
193
  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)
  {
194
    return MatVec::minus(MatVec::as_scalar(a), MatVec::as_scalar(b));
195
  }
196

197
198
199
200
  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)
  {
201
    return MatVec::multiplies(MatVec::as_scalar(a), MatVec::as_scalar(b));
202
  }
203

204
205
206
207
  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)
  {
208
    return MatVec::divides(MatVec::as_scalar(a), MatVec::as_scalar(b));
209
  }
210
211
212
213
214

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

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

217
  /// Cross-product of two 3d-vectors
218
  template <class T>
219
  FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b);
220
221
222

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

225
226
227
  template <class T, class S, int N>
  auto dot(FieldMatrix<T,1,N> const& vec1, FieldMatrix<S,1,N> const& vec2);

228
229
230
231
  // ----------------------------------------------------------------------------

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

234
  template <class T, int N>
235
  T sum(FieldMatrix<T, 1, N> const& x);
236

237
238
239

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

242
  template <class T, int N>
243
  auto unary_dot(FieldMatrix<T, 1, N> const& x);
244

245
246
  /// Maximum over all vector entries
  template <class T, int N>
247
  auto max(FieldVector<T, N> const& x);
248

249
  template <class T, int N>
250
  auto max(FieldMatrix<T, 1, N> const& x);
251

252
253
  /// Minimum over all vector entries
  template <class T, int N>
254
  auto min(FieldVector<T, N> const& x);
255

256
  template <class T, int N>
257
  auto min(FieldMatrix<T, 1, N> const& x);
258

259
260
  /// Maximum of the absolute values of vector entries
  template <class T, int N>
261
  auto abs_max(FieldVector<T, N> const& x);
262

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

266
267
  /// Minimum of the absolute values of vector entries
  template <class T, int N>
268
  auto abs_min(FieldVector<T, N> const& x);
269

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

273
274
275
276
277
278
  // ----------------------------------------------------------------------------

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

281
  template <class T, int N>
282
  auto one_norm(FieldMatrix<T, 1, N> const& x);
283

284
285
286
287
  /** \ingroup vector_norms
   *  \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 )
   **/
  template <class T, int N>
288
  auto two_norm(FieldVector<T, N> const& x);
289

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

293
294
295
296
  /** \ingroup vector_norms
   *  \brief The p-norm of a vector = ( sum_i |x_i|^p )^(1/p)
   **/
  template <int p, class T, int N>
297
  auto p_norm(FieldVector<T, N> const& x);
298

299
  template <int p, class T, int N>
300
  auto p_norm(FieldMatrix<T, 1, N> const& x);
301

302
303
304
305
  /** \ingroup vector_norms
   *  \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max
   **/
  template <class T, int N>
306
  auto infty_norm(FieldVector<T, N> const& x);
307

308
  template <class T, int N>
309
  auto infty_norm(FieldMatrix<T, 1, N> const& x);
310

311
312
313
314
  // ----------------------------------------------------------------------------

  /// The euklidean distance between two vectors = |lhs-rhs|_2
  template <class T, int N>
315
  T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs);
316
317
318
319
320

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

  /// Outer product (vec1 * vec2^T)
  template <class T, class S, int N, int M, int K>
321
  auto outer(FieldMatrix<T,N,K> const& vec1, FieldMatrix<S,M,K> const& vec2);
322
323
324
325

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

  template <class T>
326
  T det(FieldMatrix<T, 0, 0> const& /*mat*/);
327
328
329

  /// Determinant of a 1x1 matrix
  template <class T>
330
  T det(FieldMatrix<T, 1, 1> const& mat);
331
332
333

  /// Determinant of a 2x2 matrix
  template <class T>
334
  T det(FieldMatrix<T, 2, 2> const& mat);
335
336
337

  /// Determinant of a 3x3 matrix
  template <class T>
338
  T det(FieldMatrix<T, 3, 3> const& mat);
339
340
341

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

344
345
346

  /// Return the inverse of the matrix `mat`
  template <class T, int N>
347
  auto inv(FieldMatrix<T, N, N> mat);
348
349
350

  /// Solve the linear system A*x = b
  template <class T, int N>
351
  void solve(FieldMatrix<T, N, N> const& A,  FieldVector<T, N>& x,  FieldVector<T, N> const& b);
352
353
354
355


  /// Gramian determinant = sqrt( det( DT^T * DF ) )
  template <class T, int N, int M>
356
  T gramian(FieldMatrix<T,N,M> const& DF);
357
358
359

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

362
363
364
  // ----------------------------------------------------------------------------
  // some arithmetic operations with FieldMatrix

365
  template <class T, int M, int N>
366
  FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A);
367

368
369
370
371
372
373
  template <class T, int N>
  DiagonalMatrix<T,N> const& trans(DiagonalMatrix<T,N> const& A)
  {
    return A;
  }

374
  // -----------------------------------------------------------------------------
375

376
377
  template <class T1, class T2, int M, int N, int L>
  FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M> const& A,  FieldMatrix<T2, N, L> const& B);
378

379
380
  template <class T1, class T2, int M, int N, int L>
  FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B);
381

382
383
  template <class T1, class T2, class T3, int M, int N, int L>
  FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C);
384

385
386
387
388
389
  template <class T1, class T2, class T3, int M, int N>
  FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C);

  template <class T1, class T2, class T3, int N>
  FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C);
390

391
  // -----------------------------------------------------------------------------
392
393
394
395
396
397
398

  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);

399
  // necessary specialization to resolve ambiguous calls
400
401
402
  template <class T>
  T const& at(FieldMatrix<T,1,1> const& vec, std::size_t i);

403
404
405
  template <class T, int N>
  T const& at(FieldVector<T,N> const& vec, std::size_t i);

406
407
408
409
410
411
412
} // end namespace Dune

namespace AMDiS
{
  using Dune::FieldMatrix;
  using Dune::FieldVector;
}
413
414

#include "FieldMatVec.inc.hpp"