Dune Core Modules (2.4.1)

algebra.hh
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
8
9namespace 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 ( 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 K sum = matrix[ 0 ][ 0 ] * matrix[ 0 ][ 0 ];
43 for( int i = 1; i < m; ++i )
44 sum += matrix[ 0 ][ i ] * matrix[ 0 ][ i ];
45 return sqrt( sum );
46 }
47
48 template< class K >
49 inline static K determinant ( const FieldMatrix< K, 2, 2 > &matrix )
50 {
51 return matrix[ 0 ][ 0 ] * matrix[ 1 ][ 1 ] - matrix[ 0 ][ 1 ] * matrix[ 1 ][ 0 ];
52 }
53
54 template< class K >
55 inline static K determinant ( const FieldMatrix< K, 2, 3 > &matrix )
56 {
57 return vectorProduct( matrix[ 0 ], matrix[ 1 ] ).two_norm();
58 }
59
60 template< class K, int m >
61 inline static K determinant ( const FieldMatrix< K, 2, m > &matrix )
62 {
63 const K tmpA = matrix[ 0 ].two_norm2();
64 const K tmpB = matrix[ 1 ].two_norm2();
65 const K tmpC = matrix[ 0 ] * matrix[ 1 ];
66 return sqrt( tmpA * tmpB - tmpC * tmpC );
67 }
68
69 template< class K >
70 inline static K determinant ( const FieldMatrix< K, 3, 3 > &matrix )
71 {
72 return matrix[ 0 ] * vectorProduct( matrix[ 1 ], matrix[ 2 ] );
73 }
74
75
76 template< class K, int m >
77 inline static K invert ( const FieldMatrix< K, 0, m > &matrix,
78 FieldMatrix< K, m, 0 > &inverse )
79 {
80 return K( 1 );
81 }
82
83 template< class K >
84 inline static K invert ( const FieldMatrix< K, 1, 1 > &matrix,
85 FieldMatrix< K, 1, 1 > &inverse )
86 {
87 inverse[ 0 ][ 0 ] = K( 1 ) / matrix[ 0 ][ 0 ];
88 return matrix[ 0 ][ 0 ];
89 }
90
91 template< class K, int m >
92 inline static K invert ( const FieldMatrix< K, 1, m > &matrix,
93 FieldMatrix< K, m, 1 > &inverse )
94 {
95 K detSqr = matrix[ 0 ].two_norm2();
96 K invDetSqr = K( 1 ) / detSqr;
97 for( int i = 0; i < m; ++i )
98 inverse[ i ][ 0 ] = invDetSqr * matrix[ 0 ][ i ];
99 return sqrt( detSqr );
100 }
101
102 template< class K >
103 inline static K invert ( const FieldMatrix< K, 2, 2 > &matrix,
104 FieldMatrix< K, 2, 2 > &inverse )
105 {
106 K det = determinant( matrix );
107 K invDet = K( 1 ) / det;
108 inverse[ 0 ][ 0 ] = invDet * matrix[ 1 ][ 1 ];
109 inverse[ 0 ][ 1 ] = - invDet * matrix[ 0 ][ 1 ];
110 inverse[ 1 ][ 0 ] = - invDet * matrix[ 1 ][ 0 ];
111 inverse[ 1 ][ 1 ] = invDet * matrix[ 0 ][ 0 ];
112 return det;
113 }
114
115 template< class K, int m >
116 inline static K invert ( const FieldMatrix< K, 2, m > &matrix,
117 FieldMatrix< K, m, 2 > &inverse )
118 {
119 const K tmpA = matrix[ 0 ].two_norm2();
120 const K tmpB = matrix[ 1 ].two_norm2();
121 const K tmpC = matrix[ 0 ] * matrix[ 1 ];
122 const K detSqr = tmpA * tmpB - tmpC * tmpC;
123 const K invDetSqr = K( 1 ) / detSqr;
124 for( int i = 0; i < m; ++i )
125 {
126 inverse[ i ][ 0 ] = invDetSqr * (tmpB * matrix[ 0 ][ i ] - tmpC * matrix[ 1 ][ i ]);
127 inverse[ i ][ 1 ] = invDetSqr * (tmpA * matrix[ 1 ][ i ] - tmpC * matrix[ 0 ][ i ]);
128 }
129 return sqrt( detSqr );
130 }
131
132 template< class K >
133 inline static K invert ( const FieldMatrix< K, 3, 3 > &matrix,
134 FieldMatrix< K, 3, 3 > &inverse )
135 {
136 return FMatrixHelp::invertMatrix( matrix, inverse );
137 }
138 }
139
140}
141
142#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: alignment.hh:10
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)