git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6026 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/fix_rigid.cpp

@@ -648,8 +648,7 @@ void FixRigid::init()
 	sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
 	sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
 	sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
-      }
-      if (eflags[i] & INERTIA_ELLIPSOID) {
+      } else if (eflags[i] & INERTIA_ELLIPSOID) {
 	shape = ebonus[ellipsoid[i]].shape;
 	quatatom = ebonus[ellipsoid[i]].quat;
 	MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
@@ -665,11 +664,12 @@ void FixRigid::init()
   
   MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
   
+  // diagonalize inertia tensor for each body via Jacobi rotations
   // inertia = 3 eigenvalues = principal moments of inertia
-  // ex_space,ey_space,ez_space = 3 eigenvectors = principal axes of rigid body
-  
+  // evectors and exzy_space = 3 evectors = principal axes of rigid body
+
   int ierror;
-  double ez0,ez1,ez2;
+  double cross[3];
   double tensor[3][3],evectors[3][3];
 
   for (ibody = 0; ibody < nbody; ibody++) {
@@ -686,11 +686,9 @@ void FixRigid::init()
     ex_space[ibody][0] = evectors[0][0];
     ex_space[ibody][1] = evectors[1][0];
     ex_space[ibody][2] = evectors[2][0];
-    
     ey_space[ibody][0] = evectors[0][1];
     ey_space[ibody][1] = evectors[1][1];
     ey_space[ibody][2] = evectors[2][1];
-    
     ez_space[ibody][0] = evectors[0][2];
     ez_space[ibody][1] = evectors[1][2];
     ez_space[ibody][2] = evectors[2][2];
@@ -706,21 +704,11 @@ void FixRigid::init()
     if (inertia[ibody][2] < EPSILON*max) inertia[ibody][2] = 0.0;
 
     // enforce 3 evectors as a right-handed coordinate system
-    // flip 3rd evector if needed
-  
-    ez0 = ex_space[ibody][1]*ey_space[ibody][2] -
-      ex_space[ibody][2]*ey_space[ibody][1];
-    ez1 = ex_space[ibody][2]*ey_space[ibody][0] -
-      ex_space[ibody][0]*ey_space[ibody][2];
-    ez2 = ex_space[ibody][0]*ey_space[ibody][1] -
-      ex_space[ibody][1]*ey_space[ibody][0];
-  
-    if (ez0*ez_space[ibody][0] + ez1*ez_space[ibody][1] + 
-	ez2*ez_space[ibody][2] < 0.0) {
-      ez_space[ibody][0] = -ez_space[ibody][0];
-      ez_space[ibody][1] = -ez_space[ibody][1];
-      ez_space[ibody][2] = -ez_space[ibody][2];
-    }
+    // flip 3rd vector if needed
+
+    MathExtra::cross3(ex_space[ibody],ey_space[ibody],cross);
+    if (MathExtra::dot3(cross,ez_space[ibody]) < 0.0)
+      MathExtra::negate3(ez_space[ibody]);
 
     // create initial quaternion
   
@@ -823,8 +811,7 @@ void FixRigid::init()
 	sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
 	sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
 	sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
-      }
-      if (eflags[i] & INERTIA_ELLIPSOID) {
+      } else if (eflags[i] & INERTIA_ELLIPSOID) {
 	shape = ebonus[ellipsoid[i]].shape;
 	MathExtra::inertia_ellipsoid(shape,qorient[i],massone,ivec);
 	sum[ibody][0] += ivec[0];

src/math_extra.cpp

@@ -487,7 +487,7 @@ void inertia_line(double length, double theta, double mass, double *inertia)
 /* ----------------------------------------------------------------------
    compute space-frame inertia tensor of a triangle
    v0,v1,v2 = 3 vertices of triangle
-   from http://en.wikipedia.org/wiki/Inertia_tensor_of_triangle:
+   from http://en.wikipedia.org/wiki/Inertia_tensor_of_triangle
    inertia tensor = a/24 (v0^2 + v1^2 + v2^2 + (v0+v1+v2)^2) I - a Vt S V
    a = 2*area of tri = |(v1-v0) x (v2-v0)|
    I = 3x3 identity matrix
@@ -523,9 +523,7 @@ void inertia_triangle(double *v0, double *v1, double *v2,
   sub3(v2,v0,v2mv0);
   cross3(v1mv0,v2mv0,normal);
   double a = len3(normal);
-  double inv24 = 1.0/24.0;
-
-  // NOTE: use mass
+  double inv24 = mass/24.0;
 
   inertia[0] = inv24*a*(sum-vtsv[0][0]);
   inertia[1] = inv24*a*(sum-vtsv[1][1]);

src/math_extra.h

@@ -30,6 +30,7 @@ namespace MathExtra {
   inline void norm3(double *v);
   inline void normalize3(const double *v, double *ans);
   inline void snormalize3(const double, const double *v, double *ans);
+  inline void negate3(double *v);
   inline void add3(const double *v1, const double *v2, double *ans);
   inline void sub3(const double *v1, const double *v2, double *ans);
   inline double len3(const double *v);
@@ -156,6 +157,17 @@ void MathExtra::snormalize3(const double length, const double *v, double *ans)
   ans[2] = v[2]*scale;
 }
 
+/* ----------------------------------------------------------------------
+   negate vector v
+------------------------------------------------------------------------- */
+
+void MathExtra::negate3(double *v)
+{
+  v[0] = -v[0];
+  v[1] = -v[1];
+  v[2] = -v[2];
+}
+
 /* ----------------------------------------------------------------------
    ans = v1 + v2
 ------------------------------------------------------------------------- */