DUNE-FEM (unstable)

algebra.hh
1// SPDX-FileCopyrightText: Copyright © 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
10
11namespace 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
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)