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

Comparing trunk/src/math/SquareMatrix3.hpp (file contents):
Revision 633 by tim, Tue Sep 27 20:02:57 2005 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 6 | Line 6
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
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   /**
# Line 47 | Line 47
47   */
48   #ifndef MATH_SQUAREMATRIX3_HPP
49   #define  MATH_SQUAREMATRIX3_HPP
50 <
50 > #include <vector>
51   #include "Quaternion.hpp"
52   #include "SquareMatrix.hpp"
53   #include "Vector3.hpp"
54   #include "utils/NumericConstant.hpp"
55 < namespace oopse {
55 > namespace OpenMD {
56  
57    template<typename Real>
58    class SquareMatrix3 : public SquareMatrix<Real, 3> {
# Line 166 | Line 166 | namespace oopse {
166      void setupRotMat(Real w, Real x, Real y, Real z) {
167        Quaternion<Real> q(w, x, y, z);
168        *this = q.toRotationMatrix3();
169 +    }
170 +
171 +    void setupSkewMat(Vector3<Real> v) {
172 +        setupSkewMat(v[0], v[1], v[2]);
173 +    }
174 +
175 +    void setupSkewMat(Real v1, Real v2, Real v3) {
176 +        this->data_[0][0] = 0;
177 +        this->data_[0][1] = -v3;
178 +        this->data_[0][2] = v2;
179 +        this->data_[1][0] = v3;
180 +        this->data_[1][1] = 0;
181 +        this->data_[1][2] = -v1;
182 +        this->data_[2][0] = -v2;
183 +        this->data_[2][1] = v1;
184 +        this->data_[2][2] = 0;
185 +        
186 +        
187      }
188  
189 +
190 +
191      /**
192       * Returns the quaternion from this rotation matrix
193       * @return the quaternion from this rotation matrix
# Line 179 | Line 199 | namespace oopse {
199        Real ad1, ad2, ad3;    
200        t = this->data_[0][0] + this->data_[1][1] + this->data_[2][2] + 1.0;
201  
202 <      if( t > 0.0 ){
202 >      if( t > NumericConstant::epsilon ){
203  
204          s = 0.5 / sqrt( t );
205          q[0] = 0.25 / s;
# Line 224 | Line 244 | namespace oopse {
244       * @return the euler angles in a vector
245       * @exception invalid rotation matrix
246       * We use so-called "x-convention", which is the most common definition.
247 <     * In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
248 <     * rotation is by an angle phi about the z-axis, the second is by an angle  
249 <     * theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
250 <     * z-axis (again).
247 >     * In this convention, the rotation given by Euler angles (phi, theta,
248 >     * psi), where the first rotation is by an angle phi about the z-axis,
249 >     * the second is by an angle theta (0 <= theta <= 180) about the x-axis,
250 >     * and the third is by an angle psi about the z-axis (again).
251       */            
252      Vector3<Real> toEulerAngles() {
253        Vector3<Real> myEuler;
# Line 239 | Line 259 | namespace oopse {
259                  
260        // set the tolerance for Euler angles and rotation elements
261  
262 <      theta = acos(std::min(1.0, std::max(-1.0,this->data_[2][2])));
262 >      theta = acos(std::min((RealType)1.0, std::max((RealType)-1.0,this->data_[2][2])));
263        ctheta = this->data_[2][2];
264        stheta = sqrt(1.0 - ctheta * ctheta);
265  
266 <      // when sin(theta) is close to 0, we need to consider singularity
267 <      // In this case, we can assign an arbitary value to phi (or psi), and then determine
268 <      // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
269 <      // in cases of singularity.  
266 >      // when sin(theta) is close to 0, we need to consider
267 >      // singularity In this case, we can assign an arbitary value to
268 >      // phi (or psi), and then determine the psi (or phi) or
269 >      // vice-versa. We'll assume that phi always gets the rotation,
270 >      // and psi is 0 in cases of singularity.
271        // we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
272 <      // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
273 <      // change the sign of both of the parameters passed to atan2.
272 >      // Since 0 <= theta <= 180, sin(theta) will be always
273 >      // non-negative. Therefore, it will never change the sign of both of
274 >      // the parameters passed to atan2.
275  
276 <      if (fabs(stheta) <= oopse::epsilon){
276 >      if (fabs(stheta) < 1e-6){
277          psi = 0.0;
278          phi = atan2(-this->data_[1][0], this->data_[0][0]);  
279        }
# Line 263 | Line 285 | namespace oopse {
285  
286        //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
287        if (phi < 0)
288 <        phi += M_PI;
288 >        phi += 2.0 * M_PI;
289  
290        if (psi < 0)
291 <        psi += M_PI;
291 >        psi += 2.0 * M_PI;
292  
293        myEuler[0] = phi;
294        myEuler[1] = theta;
# Line 298 | Line 320 | namespace oopse {
320       */
321      SquareMatrix3<Real>  inverse() const {
322        SquareMatrix3<Real> m;
323 <      double det = determinant();
324 <      if (fabs(det) <= oopse::epsilon) {
323 >      RealType det = determinant();
324 >      if (fabs(det) <= OpenMD::epsilon) {
325          //"The method was called on a matrix with |determinant| <= 1e-6.",
326          //"This is a runtime or a programming error in your application.");
327 <      }
327 >        std::vector<int> zeroDiagElementIndex;
328 >        for (int i =0; i < 3; ++i) {
329 >            if (fabs(this->data_[i][i]) <= OpenMD::epsilon) {
330 >                zeroDiagElementIndex.push_back(i);
331 >            }
332 >        }
333  
334 <      m(0, 0) = this->data_[1][1]*this->data_[2][2] - this->data_[1][2]*this->data_[2][1];
335 <      m(1, 0) = this->data_[1][2]*this->data_[2][0] - this->data_[1][0]*this->data_[2][2];
336 <      m(2, 0) = this->data_[1][0]*this->data_[2][1] - this->data_[1][1]*this->data_[2][0];
337 <      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];
334 >        if (zeroDiagElementIndex.size() == 2) {
335 >            int index = zeroDiagElementIndex[0];
336 >            m(index, index) = 1.0 / this->data_[index][index];
337 >        }else if (zeroDiagElementIndex.size() == 1) {
338  
339 <      m /= det;
339 >            int a = (zeroDiagElementIndex[0] + 1) % 3;
340 >            int b = (zeroDiagElementIndex[0] + 2) %3;
341 >            RealType denom = this->data_[a][a] * this->data_[b][b] - this->data_[b][a]*this->data_[a][b];
342 >            m(a, a) = this->data_[b][b] /denom;
343 >            m(b, a) = -this->data_[b][a]/denom;
344 >
345 >            m(a,b) = -this->data_[a][b]/denom;
346 >            m(b, b) = this->data_[a][a]/denom;
347 >                
348 >        }
349 >      
350 > /*
351 >        for(std::vector<int>::iterator iter = zeroDiagElementIndex.begin(); iter != zeroDiagElementIndex.end() ++iter) {
352 >            if (this->data_[*iter][0] > OpenMD::epsilon || this->data_[*iter][1] ||this->data_[*iter][2] ||
353 >                this->data_[0][*iter] > OpenMD::epsilon || this->data_[1][*iter] ||this->data_[2][*iter] ) {
354 >                std::cout << "can not inverse matrix" << std::endl;
355 >            }
356 >        }
357 > */
358 >      } else {
359 >
360 >          m(0, 0) = this->data_[1][1]*this->data_[2][2] - this->data_[1][2]*this->data_[2][1];
361 >          m(1, 0) = this->data_[1][2]*this->data_[2][0] - this->data_[1][0]*this->data_[2][2];
362 >          m(2, 0) = this->data_[1][0]*this->data_[2][1] - this->data_[1][1]*this->data_[2][0];
363 >          m(0, 1) = this->data_[2][1]*this->data_[0][2] - this->data_[2][2]*this->data_[0][1];
364 >          m(1, 1) = this->data_[2][2]*this->data_[0][0] - this->data_[2][0]*this->data_[0][2];
365 >          m(2, 1) = this->data_[2][0]*this->data_[0][1] - this->data_[2][1]*this->data_[0][0];
366 >          m(0, 2) = this->data_[0][1]*this->data_[1][2] - this->data_[0][2]*this->data_[1][1];
367 >          m(1, 2) = this->data_[0][2]*this->data_[1][0] - this->data_[0][0]*this->data_[1][2];
368 >          m(2, 2) = this->data_[0][0]*this->data_[1][1] - this->data_[0][1]*this->data_[1][0];
369 >
370 >          m /= det;
371 >        }
372        return m;
373      }
374 +
375 +    SquareMatrix3<Real> transpose() const{
376 +      SquareMatrix3<Real> result;
377 +                
378 +      for (unsigned int i = 0; i < 3; i++)
379 +        for (unsigned int j = 0; j < 3; j++)              
380 +          result(j, i) = this->data_[i][j];
381 +
382 +      return result;
383 +    }
384      /**
385       * Extract the eigenvalues and eigenvectors from a 3x3 matrix.
386       * The eigenvectors (the columns of V) will be normalized.
# Line 499 | Line 563 | namespace oopse {
563    }
564  
565      
566 <  typedef SquareMatrix3<double> Mat3x3d;
567 <  typedef SquareMatrix3<double> RotMat3x3d;
566 >  typedef SquareMatrix3<RealType> Mat3x3d;
567 >  typedef SquareMatrix3<RealType> RotMat3x3d;
568  
569 < } //namespace oopse
569 > } //namespace OpenMD
570   #endif // MATH_SQUAREMATRIX_HPP
571  

Diff Legend

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