From e2707ca96f37c10fc982ea99fdde9db68b7c0784 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 28 Apr 2011 16:16:54 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6026 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/fix_rigid.cpp | 35 +++++++++++------------------------ src/math_extra.cpp | 6 ++---- src/math_extra.h | 12 ++++++++++++ 3 files changed, 25 insertions(+), 28 deletions(-) diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp index 1099f09d7a..4d0686409d 100644 --- a/src/fix_rigid.cpp +++ b/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]; diff --git a/src/math_extra.cpp b/src/math_extra.cpp index c4318e8bbe..5160262aff 100644 --- a/src/math_extra.cpp +++ b/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]); diff --git a/src/math_extra.h b/src/math_extra.h index 3ca98f8f12..44af2e9a8a 100755 --- a/src/math_extra.h +++ b/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 ------------------------------------------------------------------------- */