Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2007 Sandia National Laboratories. Developed at the
5 : : University of Wisconsin--Madison under SNL contract number
6 : : 624796. The U.S. Government and the University of Wisconsin
7 : : retain certain rights to this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : (2008) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file SymMatrix3D.hpp
28 : : * \brief Symetric 3x3 Matrix
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #ifndef MSQ_SYM_MATRIX_3D_HPP
33 : : #define MSQ_SYM_MATRIX_3D_HPP
34 : :
35 : : #include "Mesquite.hpp"
36 : : #include "Vector3D.hpp"
37 : :
38 : : namespace MBMesquite
39 : : {
40 : :
41 : : class MESQUITE_EXPORT SymMatrix3D
42 : : {
43 : : private:
44 : : double d_[6];
45 : :
46 : : public:
47 : : enum Term
48 : : {
49 : : T00 = 0,
50 : : T01 = 1,
51 : : T02 = 2,
52 : : T10 = T01,
53 : : T11 = 3,
54 : : T12 = 4,
55 : : T20 = T02,
56 : : T21 = T12,
57 : : T22 = 5
58 : : };
59 : :
60 : 0 : inline static Term term( unsigned r, unsigned c )
61 : : {
62 [ # # ]: 0 : return ( Term )( r <= c ? 3 * r - r * ( r + 1 ) / 2 + c : 3 * c - c * ( c + 1 ) / 2 + r );
63 : : }
64 : :
65 : 0 : SymMatrix3D() {}
66 : :
67 : 0 : SymMatrix3D( double diagonal_value )
68 : : {
69 : 0 : d_[T00] = d_[T11] = d_[T22] = diagonal_value;
70 : 0 : d_[T01] = d_[T02] = d_[T12] = 0.0;
71 : 0 : }
72 : :
73 : 0 : SymMatrix3D( double t00, double t01, double t02, double t11, double t12, double t22 )
74 : : {
75 : 0 : d_[T00] = t00;
76 : 0 : d_[T01] = t01;
77 : 0 : d_[T02] = t02;
78 : 0 : d_[T11] = t11;
79 : 0 : d_[T12] = t12;
80 : 0 : d_[T22] = t22;
81 : 0 : }
82 : :
83 : : /**\brief Outer product */
84 : 0 : SymMatrix3D( const Vector3D& u )
85 : : {
86 : 0 : d_[T00] = u[0] * u[0];
87 : 0 : d_[T01] = u[0] * u[1];
88 : 0 : d_[T02] = u[0] * u[2];
89 : 0 : d_[T11] = u[1] * u[1];
90 : 0 : d_[T12] = u[1] * u[2];
91 : 0 : d_[T22] = u[2] * u[2];
92 : 0 : }
93 : :
94 : 0 : double& operator[]( unsigned t )
95 : : {
96 : 0 : return d_[t];
97 : : }
98 : 0 : double operator[]( unsigned t ) const
99 : : {
100 : 0 : return d_[t];
101 : : }
102 : :
103 : 0 : double& operator()( unsigned short r, unsigned short c )
104 : : {
105 : 0 : return d_[term( r, c )];
106 : : }
107 : 0 : double operator()( unsigned short r, unsigned short c ) const
108 : : {
109 : 0 : return d_[term( r, c )];
110 : : }
111 : :
112 : : inline SymMatrix3D& operator+=( const SymMatrix3D& other );
113 : : inline SymMatrix3D& operator-=( const SymMatrix3D& other );
114 : : inline SymMatrix3D& operator*=( double scalar );
115 : : inline SymMatrix3D& operator/=( double scalar );
116 : : };
117 : :
118 : : inline SymMatrix3D operator-( const SymMatrix3D& m )
119 : : {
120 : : return SymMatrix3D( -m[SymMatrix3D::T00], -m[SymMatrix3D::T01], -m[SymMatrix3D::T02], -m[SymMatrix3D::T11],
121 : : -m[SymMatrix3D::T12], -m[SymMatrix3D::T22] );
122 : : }
123 : :
124 : 0 : inline SymMatrix3D& SymMatrix3D::operator+=( const SymMatrix3D& other )
125 : : {
126 : 0 : d_[0] += other.d_[0];
127 : 0 : d_[1] += other.d_[1];
128 : 0 : d_[2] += other.d_[2];
129 : 0 : d_[3] += other.d_[3];
130 : 0 : d_[4] += other.d_[4];
131 : 0 : d_[5] += other.d_[5];
132 : 0 : return *this;
133 : : }
134 : :
135 : 0 : inline SymMatrix3D& SymMatrix3D::operator-=( const SymMatrix3D& other )
136 : : {
137 : 0 : d_[0] -= other.d_[0];
138 : 0 : d_[1] -= other.d_[1];
139 : 0 : d_[2] -= other.d_[2];
140 : 0 : d_[3] -= other.d_[3];
141 : 0 : d_[4] -= other.d_[4];
142 : 0 : d_[5] -= other.d_[5];
143 : 0 : return *this;
144 : : }
145 : :
146 : 0 : inline SymMatrix3D& SymMatrix3D::operator*=( double s )
147 : : {
148 : 0 : d_[0] *= s;
149 : 0 : d_[1] *= s;
150 : 0 : d_[2] *= s;
151 : 0 : d_[3] *= s;
152 : 0 : d_[4] *= s;
153 : 0 : d_[5] *= s;
154 : 0 : return *this;
155 : : }
156 : :
157 : : inline SymMatrix3D& SymMatrix3D::operator/=( double s )
158 : : {
159 : : d_[0] /= s;
160 : : d_[1] /= s;
161 : : d_[2] /= s;
162 : : d_[3] /= s;
163 : : d_[4] /= s;
164 : : d_[5] /= s;
165 : : return *this;
166 : : }
167 : :
168 : : inline SymMatrix3D operator+( const SymMatrix3D& a, const SymMatrix3D& b )
169 : : {
170 : : SymMatrix3D r( a );
171 : : r += b;
172 : : return r;
173 : : }
174 : 0 : inline SymMatrix3D operator-( const SymMatrix3D& a, const SymMatrix3D& b )
175 : : {
176 : 0 : SymMatrix3D r( a );
177 : 0 : r -= b;
178 : 0 : return r;
179 : : }
180 : : inline SymMatrix3D operator*( const SymMatrix3D& a, double s )
181 : : {
182 : : SymMatrix3D r( a );
183 : : r *= s;
184 : : return r;
185 : : }
186 : 0 : inline SymMatrix3D operator*( double s, const SymMatrix3D& a )
187 : : {
188 : 0 : SymMatrix3D r( a );
189 : 0 : r *= s;
190 : 0 : return r;
191 : : }
192 : : inline SymMatrix3D operator/( const SymMatrix3D& a, double s )
193 : : {
194 : : SymMatrix3D r( a );
195 : : r /= s;
196 : : return r;
197 : : }
198 : : inline SymMatrix3D operator/( double s, const SymMatrix3D& a )
199 : : {
200 : : SymMatrix3D r( a );
201 : : r /= s;
202 : : return r;
203 : : }
204 : :
205 : : inline Vector3D operator*( const Vector3D& v, const SymMatrix3D& m )
206 : : {
207 : : return Vector3D( v[0] * m[0] + v[1] * m[1] + v[2] * m[2], v[0] * m[1] + v[1] * m[3] + v[2] * m[4],
208 : : v[0] * m[2] + v[1] * m[4] + v[2] * m[5] );
209 : : }
210 : : inline Vector3D operator*( const SymMatrix3D& m, const Vector3D& v )
211 : : {
212 : : return v * m;
213 : : }
214 : :
215 : : /** Calculate the outer product of a vector with itself */
216 : 0 : inline SymMatrix3D outer( const Vector3D& v )
217 : : {
218 : 0 : return SymMatrix3D( v[0] * v[0], v[0] * v[1], v[0] * v[2], v[1] * v[1], v[1] * v[2], v[2] * v[2] );
219 : : }
220 : :
221 : : /** Given to vectors u and v, calculate the symmetric matrix
222 : : * equal to outer(u,v) + transpose(outer(u,v))
223 : : * equal to outer(v,u) + transpose(outer(v,u))
224 : : */
225 : : inline SymMatrix3D outer_plus_transpose( const Vector3D& u, const Vector3D& v )
226 : : {
227 : : return SymMatrix3D( 2 * u[0] * v[0], u[0] * v[1] + u[1] * v[0], u[0] * v[2] + u[2] * v[0], 2 * u[1] * v[1],
228 : : u[1] * v[2] + u[2] * v[1], 2 * u[2] * v[2] );
229 : : }
230 : :
231 : : inline const SymMatrix3D& transpose( const SymMatrix3D& a )
232 : : {
233 : : return a;
234 : : }
235 : :
236 : : inline double det( const SymMatrix3D& a )
237 : : {
238 : : return a[0] * a[3] * a[5] + 2.0 * a[1] * a[2] * a[4] - a[0] * a[4] * a[4] - a[3] * a[2] * a[2] - a[5] * a[1] * a[1];
239 : : }
240 : :
241 : : inline SymMatrix3D inverse( const SymMatrix3D& a )
242 : : {
243 : : SymMatrix3D result( a[3] * a[5] - a[4] * a[4], a[2] * a[4] - a[1] * a[5], a[1] * a[4] - a[2] * a[3],
244 : : a[0] * a[5] - a[2] * a[2], a[1] * a[2] - a[0] * a[4], a[0] * a[3] - a[1] * a[1] );
245 : : result /= det( a );
246 : : return result;
247 : : }
248 : :
249 : 0 : inline double Frobenius_2( const SymMatrix3D& a )
250 : : {
251 : 0 : return a[0] * a[0] + 2 * a[1] * a[1] + 2 * a[2] * a[2] + a[3] * a[3] + 2 * a[4] * a[5] + a[5] * a[5];
252 : : }
253 : :
254 : : inline double Frobenius( const SymMatrix3D& a )
255 : : {
256 : : return std::sqrt( Frobenius_2( a ) );
257 : : }
258 : :
259 : : } // namespace MBMesquite
260 : :
261 : : #endif
|
