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

10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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;
  };
}

25
namespace Dune
26
{
27
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>
    decltype(auto) simplify(T&& t)
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
85
      return FWD(t);
86
87
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
    }

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

162
163
164
165
166
167
  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));
  }
168

169
170
171
172
173
174
  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));
  }
175

176
177
178
179
180
181
  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));
  }
182

183
184
185
186
187
188
  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));
  }
189
190
191
192
193

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

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

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

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

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

207
208
209
210
  // ----------------------------------------------------------------------------

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

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

216
217
218

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

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

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

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

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

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

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

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

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

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

252
253
254
255
256
257
  // ----------------------------------------------------------------------------

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

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

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

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

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

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

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

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

290
291
292
293
  // ----------------------------------------------------------------------------

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

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

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

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

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

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

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

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

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

323
324
325

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

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


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

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

341
342
343
  // ----------------------------------------------------------------------------
  // some arithmetic operations with FieldMatrix

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

347
  // -----------------------------------------------------------------------------
348

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

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

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

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

361
  // -----------------------------------------------------------------------------
362
363
364
365
366
367
368

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

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

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

376
377
378
379
380
381
382
} // end namespace Dune

namespace AMDiS
{
  using Dune::FieldMatrix;
  using Dune::FieldVector;
}
383
384

#include "FieldMatVec.inc.hpp"