dune-grid  2.8.0
algebra.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ALBERTA_ALGEBRA_HH
4 #define DUNE_ALBERTA_ALGEBRA_HH
5 
6 #include <dune/common/fvector.hh>
7 #include <dune/common/fmatrix.hh>
8 
9 namespace Dune
10 {
11 
12  namespace Alberta
13  {
14 
15  template< class K >
16  inline static FieldVector< K, 3 >
17  vectorProduct ( const FieldVector< K, 3 > &u, const FieldVector< K, 3 > &v )
18  {
19  FieldVector< K, 3 > w;
20  w[ 0 ] = u[ 1 ] * v[ 2 ] - u[ 2 ] * v[ 1 ];
21  w[ 1 ] = u[ 2 ] * v[ 0 ] - u[ 0 ] * v[ 2 ];
22  w[ 2 ] = u[ 0 ] * v[ 1 ] - u[ 1 ] * v[ 0 ];
23  return w;
24  }
25 
26 
27  template< class K, int m >
28  inline static K determinant ( [[maybe_unused]] const FieldMatrix< K, 0, m > &matrix )
29  {
30  return K( 1 );
31  }
32 
33  template< class K >
34  inline static K determinant ( const FieldMatrix< K, 1, 1 > &matrix )
35  {
36  return matrix[ 0 ][ 0 ];
37  }
38 
39  template< class K, int m >
40  inline static K determinant ( const FieldMatrix< K, 1, m > &matrix )
41  {
42  using std::sqrt;
43  K sum = matrix[ 0 ][ 0 ] * matrix[ 0 ][ 0 ];
44  for( int i = 1; i < m; ++i )
45  sum += matrix[ 0 ][ i ] * matrix[ 0 ][ i ];
46  return sqrt( sum );
47  }
48 
49  template< class K >
50  inline static K determinant ( const FieldMatrix< K, 2, 2 > &matrix )
51  {
52  return matrix[ 0 ][ 0 ] * matrix[ 1 ][ 1 ] - matrix[ 0 ][ 1 ] * matrix[ 1 ][ 0 ];
53  }
54 
55  template< class K >
56  inline static K determinant ( const FieldMatrix< K, 2, 3 > &matrix )
57  {
58  return vectorProduct( matrix[ 0 ], matrix[ 1 ] ).two_norm();
59  }
60 
61  template< class K, int m >
62  inline static K determinant ( const FieldMatrix< K, 2, m > &matrix )
63  {
64  using std::sqrt;
65  const K tmpA = matrix[ 0 ].two_norm2();
66  const K tmpB = matrix[ 1 ].two_norm2();
67  const K tmpC = matrix[ 0 ] * matrix[ 1 ];
68  return sqrt( tmpA * tmpB - tmpC * tmpC );
69  }
70 
71  template< class K >
72  inline static K determinant ( const FieldMatrix< K, 3, 3 > &matrix )
73  {
74  return matrix[ 0 ] * vectorProduct( matrix[ 1 ], matrix[ 2 ] );
75  }
76 
77 
78  template< class K, int m >
79  inline static K invert ( [[maybe_unused]] const FieldMatrix< K, 0, m > &matrix,
80  [[maybe_unused]] FieldMatrix< K, m, 0 > &inverse )
81  {
82  return K( 1 );
83  }
84 
85  template< class K >
86  inline static K invert ( const FieldMatrix< K, 1, 1 > &matrix,
87  FieldMatrix< K, 1, 1 > &inverse )
88  {
89  inverse[ 0 ][ 0 ] = K( 1 ) / matrix[ 0 ][ 0 ];
90  return matrix[ 0 ][ 0 ];
91  }
92 
93  template< class K, int m >
94  inline static K invert ( const FieldMatrix< K, 1, m > &matrix,
95  FieldMatrix< K, m, 1 > &inverse )
96  {
97  using std::sqrt;
98  K detSqr = matrix[ 0 ].two_norm2();
99  K invDetSqr = K( 1 ) / detSqr;
100  for( int i = 0; i < m; ++i )
101  inverse[ i ][ 0 ] = invDetSqr * matrix[ 0 ][ i ];
102  return sqrt( detSqr );
103  }
104 
105  template< class K >
106  inline static K invert ( const FieldMatrix< K, 2, 2 > &matrix,
107  FieldMatrix< K, 2, 2 > &inverse )
108  {
109  K det = determinant( matrix );
110  K invDet = K( 1 ) / det;
111  inverse[ 0 ][ 0 ] = invDet * matrix[ 1 ][ 1 ];
112  inverse[ 0 ][ 1 ] = - invDet * matrix[ 0 ][ 1 ];
113  inverse[ 1 ][ 0 ] = - invDet * matrix[ 1 ][ 0 ];
114  inverse[ 1 ][ 1 ] = invDet * matrix[ 0 ][ 0 ];
115  return det;
116  }
117 
118  template< class K, int m >
119  inline static K invert ( const FieldMatrix< K, 2, m > &matrix,
120  FieldMatrix< K, m, 2 > &inverse )
121  {
122  using std::sqrt;
123  const K tmpA = matrix[ 0 ].two_norm2();
124  const K tmpB = matrix[ 1 ].two_norm2();
125  const K tmpC = matrix[ 0 ] * matrix[ 1 ];
126  const K detSqr = tmpA * tmpB - tmpC * tmpC;
127  const K invDetSqr = K( 1 ) / detSqr;
128  for( int i = 0; i < m; ++i )
129  {
130  inverse[ i ][ 0 ] = invDetSqr * (tmpB * matrix[ 0 ][ i ] - tmpC * matrix[ 1 ][ i ]);
131  inverse[ i ][ 1 ] = invDetSqr * (tmpA * matrix[ 1 ][ i ] - tmpC * matrix[ 0 ][ i ]);
132  }
133  return sqrt( detSqr );
134  }
135 
136  template< class K >
137  inline static K invert ( const FieldMatrix< K, 3, 3 > &matrix,
138  FieldMatrix< K, 3, 3 > &inverse )
139  {
140  return FMatrixHelp::invertMatrix( matrix, inverse );
141  }
142  }
143 
144 }
145 
146 #endif // #ifndef DUNE_ALBERTA_ALGEBRA_HH
Include standard header files.
Definition: agrid.hh:58
static K invert([[maybe_unused]] const FieldMatrix< K, 0, m > &matrix, [[maybe_unused]] FieldMatrix< K, m, 0 > &inverse)
Definition: algebra.hh:79
static K determinant([[maybe_unused]] const FieldMatrix< K, 0, m > &matrix)
Definition: algebra.hh:28
static FieldVector< K, 3 > vectorProduct(const FieldVector< K, 3 > &u, const FieldVector< K, 3 > &v)
Definition: algebra.hh:17