ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/primitives/RigidBody.cpp
(Generate patch)

Comparing trunk/src/primitives/RigidBody.cpp (file contents):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 770 by tim, Fri Dec 2 15:38:03 2005 UTC

# Line 164 | Line 164 | namespace oopse {
164      }
165  
166      // Moment of Inertia calculation
167 <    Mat3x3d Itmp(0.0);
168 <  
167 >    Mat3x3d Itmp(0.0);    
168      for (std::size_t i = 0; i < atoms_.size(); i++) {
169 +      Mat3x3d IAtom(0.0);  
170        mtmp = atoms_[i]->getMass();
171 <      Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
171 >      IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
172        double r2 = refCoords_[i].lengthSquare();
173 <      Itmp(0, 0) += mtmp * r2;
174 <      Itmp(1, 1) += mtmp * r2;
175 <      Itmp(2, 2) += mtmp * r2;
176 <    }
173 >      IAtom(0, 0) += mtmp * r2;
174 >      IAtom(1, 1) += mtmp * r2;
175 >      IAtom(2, 2) += mtmp * r2;
176 >      Itmp += IAtom;
177  
178 <    //project the inertial moment of directional atoms into this rigid body
179 <    for (std::size_t i = 0; i < atoms_.size(); i++) {
178 >      //project the inertial moment of directional atoms into this rigid body
179        if (atoms_[i]->isDirectional()) {
180 <        RectMatrix<double, 3, 3> Iproject = refOrients_[i].transpose() * atoms_[i]->getI();
181 <        Itmp(0, 0) += Iproject(0, 0);
183 <        Itmp(1, 1) += Iproject(1, 1);
184 <        Itmp(2, 2) += Iproject(2, 2);
185 <      }
180 >        Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i];
181 >      }
182      }
183  
184 +    //    std::cout << Itmp << std::endl;
185 +
186      //diagonalize
187      Vector3d evals;
188      Mat3x3d::diagonalize(Itmp, evals, sU_);
# Line 273 | Line 271 | namespace oopse {
271        if (atoms_[i]->isDirectional()) {
272            
273          dAtom = (DirectionalAtom *) atoms_[i];
274 <        dAtom->setA(a * refOrients_[i]);
277 <        //dAtom->rotateBy( A );      
274 >        dAtom->setA(refOrients_[i] * a);
275        }
276  
277      }
# Line 301 | Line 298 | namespace oopse {
298        if (atoms_[i]->isDirectional()) {
299            
300          dAtom = (DirectionalAtom *) atoms_[i];
301 <        dAtom->setA(a * refOrients_[i], frame);
301 >        dAtom->setA(refOrients_[i] * a, frame);
302        }
303  
304      }
# Line 486 | Line 483 | namespace oopse {
483                 "RigidBody error.\n"
484                 "\tAtom %s does not have a position specified.\n"
485                 "\tThis means RigidBody cannot set up reference coordinates.\n",
486 <               ats->getType() );
486 >               ats->getType().c_str() );
487        painCave.isFatal = 1;
488        simError();
489      }
# Line 506 | Line 503 | namespace oopse {
503                   "RigidBody error.\n"
504                   "\tAtom %s does not have an orientation specified.\n"
505                   "\tThis means RigidBody cannot set up reference orientations.\n",
506 <                 ats->getType() );
506 >                 ats->getType().c_str() );
507          painCave.isFatal = 1;
508          simError();
509        }    

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines