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

3
4
#include <type_traits>

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

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

24
namespace Dune
25
{
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
85
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
  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)
    {
      return std::forward<T>(t);
    }

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

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

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

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

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

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

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

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

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

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

206
207
208
209
  // ----------------------------------------------------------------------------

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

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

215
216
217

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

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

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

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

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

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

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

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

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

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

251
252
253
254
255
256
  // ----------------------------------------------------------------------------

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

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

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

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

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

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

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

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

289
290
291
292
  // ----------------------------------------------------------------------------

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

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

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

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

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

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

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

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

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

322
323
324

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

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


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

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

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

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

346
  // -----------------------------------------------------------------------------
347

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

  template <class T, int M, int N, int L>
352
  FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> 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, FieldMatrix<T,M,N>& C);
356

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

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

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

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

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

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

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

#include "FieldMatVec.inc.hpp"