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