ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/math/SquareMatrix3.hpp
(Generate patch)

Comparing trunk/src/math/SquareMatrix3.hpp (file contents):
Revision 123 by tim, Wed Oct 20 18:07:08 2004 UTC vs.
Revision 451 by tim, Tue Mar 29 21:00:54 2005 UTC

# Line 1 | Line 1
1 < /*
2 < * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
3 < *
4 < * Contact: oopse@oopse.org
5 < *
6 < * This program is free software; you can redistribute it and/or
7 < * modify it under the terms of the GNU Lesser General Public License
8 < * as published by the Free Software Foundation; either version 2.1
9 < * of the License, or (at your option) any later version.
10 < * All we ask is that proper credit is given for our work, which includes
11 < * - but is not limited to - adding the above copyright notice to the beginning
12 < * of your source code files, and to any copyright notice that you may distribute
13 < * with programs based on this work.
14 < *
15 < * This program is distributed in the hope that it will be useful,
16 < * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 < * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 < * GNU Lesser General Public License for more details.
19 < *
20 < * You should have received a copy of the GNU Lesser General Public License
21 < * along with this program; if not, write to the Free Software
22 < * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
1 > /*
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4 + * The University of Notre Dame grants you ("Licensee") a
5 + * non-exclusive, royalty free, license to use, modify and
6 + * redistribute this software in source and binary code form, provided
7 + * that the following conditions are met:
8 + *
9 + * 1. Acknowledgement of the program authors must be made in any
10 + *    publication of scientific results based in part on use of the
11 + *    program.  An acceptable form of acknowledgement is citation of
12 + *    the article in which the program was described (Matthew
13 + *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 + *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 + *    Parallel Simulation Engine for Molecular Dynamics,"
16 + *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 + *
18 + * 2. Redistributions of source code must retain the above copyright
19 + *    notice, this list of conditions and the following disclaimer.
20 + *
21 + * 3. Redistributions in binary form must reproduce the above copyright
22 + *    notice, this list of conditions and the following disclaimer in the
23 + *    documentation and/or other materials provided with the
24 + *    distribution.
25 + *
26 + * This software is provided "AS IS," without a warranty of any
27 + * kind. All express or implied conditions, representations and
28 + * warranties, including any implied warranty of merchantability,
29 + * fitness for a particular purpose or non-infringement, are hereby
30 + * excluded.  The University of Notre Dame and its licensors shall not
31 + * be liable for any damages suffered by licensee as a result of
32 + * using, modifying or distributing the software or its
33 + * derivatives. In no event will the University of Notre Dame or its
34 + * licensors be liable for any lost revenue, profit or data, or for
35 + * direct, indirect, special, consequential, incidental or punitive
36 + * damages, however caused and regardless of the theory of liability,
37 + * arising out of the use of or inability to use software, even if the
38 + * University of Notre Dame has been advised of the possibility of
39 + * such damages.
40   */
41 <
41 >
42   /**
43   * @file SquareMatrix3.hpp
44   * @author Teng Lin
# Line 35 | Line 51
51   #include "Quaternion.hpp"
52   #include "SquareMatrix.hpp"
53   #include "Vector3.hpp"
54 <
54 > #include "utils/NumericConstant.hpp"
55   namespace oopse {
56  
57      template<typename Real>
58      class SquareMatrix3 : public SquareMatrix<Real, 3> {
59          public:
60 +
61 +            typedef Real ElemType;
62 +            typedef Real* ElemPoinerType;
63              
64              /** default constructor */
65              SquareMatrix3() : SquareMatrix<Real, 3>() {
66              }
67  
68 +            /** Constructs and initializes every element of this matrix to a scalar */
69 +            SquareMatrix3(Real s) : SquareMatrix<Real,3>(s){
70 +            }
71 +
72 +            /** Constructs and initializes from an array */
73 +            SquareMatrix3(Real* array) : SquareMatrix<Real,3>(array){
74 +            }
75 +
76 +
77              /** copy  constructor */
78              SquareMatrix3(const SquareMatrix<Real, 3>& m)  : SquareMatrix<Real, 3>(m) {
79              }
80 <
80 >            
81              SquareMatrix3( const Vector3<Real>& eulerAngles) {
82                  setupRotMat(eulerAngles);
83              }
# Line 75 | Line 103 | namespace oopse {
103                   return *this;
104              }
105  
106 +
107 +            SquareMatrix3<Real>& operator =(const Quaternion<Real>& q) {
108 +                this->setupRotMat(q);
109 +                return *this;
110 +            }
111 +
112              /**
113               * Sets this matrix to a rotation matrix by three euler angles
114               * @ param euler
# Line 100 | Line 134 | namespace oopse {
134                  ctheta = cos(theta);
135                  cpsi = cos(psi);
136  
137 <                data_[0][0] = cpsi * cphi - ctheta * sphi * spsi;
138 <                data_[0][1] = cpsi * sphi + ctheta * cphi * spsi;
139 <                data_[0][2] = spsi * stheta;
137 >                this->data_[0][0] = cpsi * cphi - ctheta * sphi * spsi;
138 >                this->data_[0][1] = cpsi * sphi + ctheta * cphi * spsi;
139 >                this->data_[0][2] = spsi * stheta;
140                  
141 <                data_[1][0] = -spsi * ctheta - ctheta * sphi * cpsi;
142 <                data_[1][1] = -spsi * stheta + ctheta * cphi * cpsi;
143 <                data_[1][2] = cpsi * stheta;
141 >                this->data_[1][0] = -spsi * ctheta - ctheta * sphi * cpsi;
142 >                this->data_[1][1] = -spsi * stheta + ctheta * cphi * cpsi;
143 >                this->data_[1][2] = cpsi * stheta;
144  
145 <                data_[2][0] = stheta * sphi;
146 <                data_[2][1] = -stheta * cphi;
147 <                data_[2][2] = ctheta;
145 >                this->data_[2][0] = stheta * sphi;
146 >                this->data_[2][1] = -stheta * cphi;
147 >                this->data_[2][2] = ctheta;
148              }
149  
150  
# Line 143 | Line 177 | namespace oopse {
177                  Quaternion<Real> q;
178                  Real t, s;
179                  Real ad1, ad2, ad3;    
180 <                t = data_[0][0] + data_[1][1] + data_[2][2] + 1.0;
180 >                t = this->data_[0][0] + this->data_[1][1] + this->data_[2][2] + 1.0;
181  
182 <                if( t > 0.0 ){
182 >                if( t > NumericConstant::epsilon ){
183  
184                      s = 0.5 / sqrt( t );
185                      q[0] = 0.25 / s;
186 <                    q[1] = (data_[1][2] - data_[2][1]) * s;
187 <                    q[2] = (data_[2][0] - data_[0][2]) * s;
188 <                    q[3] = (data_[0][1] - data_[1][0]) * s;
186 >                    q[1] = (this->data_[1][2] - this->data_[2][1]) * s;
187 >                    q[2] = (this->data_[2][0] - this->data_[0][2]) * s;
188 >                    q[3] = (this->data_[0][1] - this->data_[1][0]) * s;
189                  } else {
190  
191 <                    ad1 = fabs( data_[0][0] );
192 <                    ad2 = fabs( data_[1][1] );
193 <                    ad3 = fabs( data_[2][2] );
191 >                    ad1 = fabs( this->data_[0][0] );
192 >                    ad2 = fabs( this->data_[1][1] );
193 >                    ad3 = fabs( this->data_[2][2] );
194  
195                      if( ad1 >= ad2 && ad1 >= ad3 ){
196  
197 <                        s = 2.0 * sqrt( 1.0 + data_[0][0] - data_[1][1] - data_[2][2] );
198 <                        q[0] = (data_[1][2] + data_[2][1]) / s;
199 <                        q[1] = 0.5 / s;
200 <                        q[2] = (data_[0][1] + data_[1][0]) / s;
201 <                        q[3] = (data_[0][2] + data_[2][0]) / s;
197 >                        s = 0.5 / sqrt( 1.0 + this->data_[0][0] - this->data_[1][1] - this->data_[2][2] );
198 >                        q[0] = (this->data_[1][2] - this->data_[2][1]) * s;
199 >                        q[1] = 0.25 / s;
200 >                        q[2] = (this->data_[0][1] + this->data_[1][0]) * s;
201 >                        q[3] = (this->data_[0][2] + this->data_[2][0]) * s;
202                      } else if ( ad2 >= ad1 && ad2 >= ad3 ) {
203 <                        s = sqrt( 1.0 + data_[1][1] - data_[0][0] - data_[2][2] ) * 2.0;
204 <                        q[0] = (data_[0][2] + data_[2][0]) / s;
205 <                        q[1] = (data_[0][1] + data_[1][0]) / s;
206 <                        q[2] = 0.5 / s;
207 <                        q[3] = (data_[1][2] + data_[2][1]) / s;
203 >                        s = 0.5 / sqrt( 1.0 + this->data_[1][1] - this->data_[0][0] - this->data_[2][2] );
204 >                        q[0] = (this->data_[2][0] - this->data_[0][2] ) * s;
205 >                        q[1] = (this->data_[0][1] + this->data_[1][0]) * s;
206 >                        q[2] = 0.25 / s;
207 >                        q[3] = (this->data_[1][2] + this->data_[2][1]) * s;
208                      } else {
209  
210 <                        s = sqrt( 1.0 + data_[2][2] - data_[0][0] - data_[1][1] ) * 2.0;
211 <                        q[0] = (data_[0][1] + data_[1][0]) / s;
212 <                        q[1] = (data_[0][2] + data_[2][0]) / s;
213 <                        q[2] = (data_[1][2] + data_[2][1]) / s;
214 <                        q[3] = 0.5 / s;
210 >                        s = 0.5 / sqrt( 1.0 + this->data_[2][2] - this->data_[0][0] - this->data_[1][1] );
211 >                        q[0] = (this->data_[0][1] - this->data_[1][0]) * s;
212 >                        q[1] = (this->data_[0][2] + this->data_[2][0]) * s;
213 >                        q[2] = (this->data_[1][2] + this->data_[2][1]) * s;
214 >                        q[3] = 0.25 / s;
215                      }
216                  }            
217  
# Line 197 | Line 231 | namespace oopse {
231              */            
232              Vector3<Real> toEulerAngles() {
233                  Vector3<Real> myEuler;
234 <                Real phi,theta,psi,eps;
235 <                Real ctheta,stheta;
234 >                Real phi;
235 >                Real theta;
236 >                Real psi;
237 >                Real ctheta;
238 >                Real stheta;
239                  
240                  // set the tolerance for Euler angles and rotation elements
241  
242 <                theta = acos(std::min(1.0, std::max(-1.0,data_[2][2])));
243 <                ctheta = data_[2][2];
242 >                theta = acos(std::min(1.0, std::max(-1.0,this->data_[2][2])));
243 >                ctheta = this->data_[2][2];
244                  stheta = sqrt(1.0 - ctheta * ctheta);
245  
246                  // when sin(theta) is close to 0, we need to consider singularity
# Line 216 | Line 253 | namespace oopse {
253  
254                  if (fabs(stheta) <= oopse::epsilon){
255                      psi = 0.0;
256 <                    phi = atan2(-data_[1][0], data_[0][0]);  
256 >                    phi = atan2(-this->data_[1][0], this->data_[0][0]);  
257                  }
258                  // we only have one unique solution
259                  else{    
260 <                    phi = atan2(data_[2][0], -data_[2][1]);
261 <                    psi = atan2(data_[0][2], data_[1][2]);
260 >                    phi = atan2(this->data_[2][0], -this->data_[2][1]);
261 >                    psi = atan2(this->data_[0][2], this->data_[1][2]);
262                  }
263  
264                  //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
# Line 242 | Line 279 | namespace oopse {
279              Real determinant() const {
280                  Real x,y,z;
281  
282 <                x = data_[0][0] * (data_[1][1] * data_[2][2] - data_[1][2] * data_[2][1]);
283 <                y = data_[0][1] * (data_[1][2] * data_[2][0] - data_[1][0] * data_[2][2]);
284 <                z = data_[0][2] * (data_[1][0] * data_[2][1] - data_[1][1] * data_[2][0]);
282 >                x = this->data_[0][0] * (this->data_[1][1] * this->data_[2][2] - this->data_[1][2] * this->data_[2][1]);
283 >                y = this->data_[0][1] * (this->data_[1][2] * this->data_[2][0] - this->data_[1][0] * this->data_[2][2]);
284 >                z = this->data_[0][2] * (this->data_[1][0] * this->data_[2][1] - this->data_[1][1] * this->data_[2][0]);
285  
286                  return(x + y + z);
287              }            
288 +
289 +            /** Returns the trace of this matrix. */
290 +            Real trace() const {
291 +                return this->data_[0][0] + this->data_[1][1] + this->data_[2][2];
292 +            }
293              
294              /**
295               * Sets the value of this matrix to  the inversion of itself.
296               * @note since simple algorithm can be applied to inverse the 3 by 3 matrix, we hide the
297               * implementation of inverse in SquareMatrix class
298               */
299 <            SquareMatrix3<Real>  inverse() {
299 >            SquareMatrix3<Real>  inverse() const {
300                  SquareMatrix3<Real> m;
301                  double det = determinant();
302                  if (fabs(det) <= oopse::epsilon) {
# Line 262 | Line 304 | namespace oopse {
304                  //"This is a runtime or a programming error in your application.");
305                  }
306  
307 <                m(0, 0) = data_[1][1]*data_[2][2] - data_[1][2]*data_[2][1];
308 <                m(1, 0) = data_[1][2]*data_[2][0] - data_[1][0]*data_[2][2];
309 <                m(2, 0) = data_[1][0]*data_[2][1] - data_[1][1]*data_[2][0];
310 <                m(0, 1) = data_[2][1]*data_[0][2] - data_[2][2]*data_[0][1];
311 <                m(1, 1) = data_[2][2]*data_[0][0] - data_[2][0]*data_[0][2];
312 <                m(2, 1) = data_[2][0]*data_[0][1] - data_[2][1]*data_[0][0];
313 <                m(0, 2) = data_[0][1]*data_[1][2] - data_[0][2]*data_[1][1];
314 <                m(1, 2) = data_[0][2]*data_[1][0] - data_[0][0]*data_[1][2];
315 <                m(2, 2) = data_[0][0]*data_[1][1] - data_[0][1]*data_[1][0];
307 >                m(0, 0) = this->data_[1][1]*this->data_[2][2] - this->data_[1][2]*this->data_[2][1];
308 >                m(1, 0) = this->data_[1][2]*this->data_[2][0] - this->data_[1][0]*this->data_[2][2];
309 >                m(2, 0) = this->data_[1][0]*this->data_[2][1] - this->data_[1][1]*this->data_[2][0];
310 >                m(0, 1) = this->data_[2][1]*this->data_[0][2] - this->data_[2][2]*this->data_[0][1];
311 >                m(1, 1) = this->data_[2][2]*this->data_[0][0] - this->data_[2][0]*this->data_[0][2];
312 >                m(2, 1) = this->data_[2][0]*this->data_[0][1] - this->data_[2][1]*this->data_[0][0];
313 >                m(0, 2) = this->data_[0][1]*this->data_[1][2] - this->data_[0][2]*this->data_[1][1];
314 >                m(1, 2) = this->data_[0][2]*this->data_[1][0] - this->data_[0][0]*this->data_[1][2];
315 >                m(2, 2) = this->data_[0][0]*this->data_[1][1] - this->data_[0][1]*this->data_[1][0];
316  
317                  m /= det;
318                  return m;
# Line 424 | Line 466 | namespace oopse {
466          v = v.transpose();
467          return ;
468      }
469 +
470 +    /**
471 +    * Return the multiplication of two matrixes  (m1 * m2).
472 +    * @return the multiplication of two matrixes
473 +    * @param m1 the first matrix
474 +    * @param m2 the second matrix
475 +    */
476 +    template<typename Real>
477 +    inline SquareMatrix3<Real> operator *(const SquareMatrix3<Real>& m1, const SquareMatrix3<Real>& m2) {
478 +        SquareMatrix3<Real> result;
479 +
480 +            for (unsigned int i = 0; i < 3; i++)
481 +                for (unsigned int j = 0; j < 3; j++)
482 +                    for (unsigned int k = 0; k < 3; k++)
483 +                        result(i, j)  += m1(i, k) * m2(k, j);                
484 +
485 +        return result;
486 +    }
487 +
488 +    template<typename Real>
489 +    inline SquareMatrix3<Real> outProduct(const Vector3<Real>& v1, const Vector3<Real>& v2) {
490 +        SquareMatrix3<Real> result;
491 +
492 +            for (unsigned int i = 0; i < 3; i++) {
493 +                for (unsigned int j = 0; j < 3; j++) {
494 +                        result(i, j)  = v1[i] * v2[j];                
495 +                }
496 +            }
497 +            
498 +        return result;        
499 +    }
500 +
501 +    
502      typedef SquareMatrix3<double> Mat3x3d;
503      typedef SquareMatrix3<double> RotMat3x3d;
504  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines