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 642 by tim, Mon Oct 3 15:54:23 2005 UTC vs.
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC

# Line 51 | Line 51 | namespace oopse {
51  
52    void RigidBody::setPrevA(const RotMat3x3d& a) {
53      ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
54 –    //((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
54  
55      for (int i =0 ; i < atoms_.size(); ++i){
56        if (atoms_[i]->isDirectional()) {
57 <        atoms_[i]->setPrevA(a * refOrients_[i]);
57 >        atoms_[i]->setPrevA(refOrients_[i].transpose() * a);
58        }
59      }
60  
# Line 64 | Line 63 | namespace oopse {
63        
64    void RigidBody::setA(const RotMat3x3d& a) {
65      ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
67 –    //((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
66  
67      for (int i =0 ; i < atoms_.size(); ++i){
68        if (atoms_[i]->isDirectional()) {
69 <        atoms_[i]->setA(a * refOrients_[i]);
69 >        atoms_[i]->setA(refOrients_[i].transpose() * a);
70        }
71      }
72    }    
# Line 79 | Line 77 | namespace oopse {
77  
78      for (int i =0 ; i < atoms_.size(); ++i){
79        if (atoms_[i]->isDirectional()) {
80 <        atoms_[i]->setA(a * refOrients_[i], snapshotNo);
80 >        atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo);
81        }
82      }
83  
# Line 89 | Line 87 | namespace oopse {
87      return inertiaTensor_;
88    }    
89  
90 <  std::vector<double> RigidBody::getGrad() {
91 <    std::vector<double> grad(6, 0.0);
90 >  std::vector<RealType> RigidBody::getGrad() {
91 >    std::vector<RealType> grad(6, 0.0);
92      Vector3d force;
93      Vector3d torque;
94      Vector3d myEuler;
95 <    double phi, theta, psi;
96 <    double cphi, sphi, ctheta, stheta;
95 >    RealType phi, theta, psi;
96 >    RealType cphi, sphi, ctheta, stheta;
97      Vector3d ephi;
98      Vector3d etheta;
99      Vector3d epsi;
# Line 148 | Line 146 | namespace oopse {
146  
147    /**@todo need modification */
148    void  RigidBody::calcRefCoords() {
149 <    double mtmp;
149 >    RealType mtmp;
150      Vector3d refCOM(0.0);
151      mass_ = 0.0;
152      for (std::size_t i = 0; i < atoms_.size(); ++i) {
# Line 169 | Line 167 | namespace oopse {
167        Mat3x3d IAtom(0.0);  
168        mtmp = atoms_[i]->getMass();
169        IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
170 <      double r2 = refCoords_[i].lengthSquare();
170 >      RealType r2 = refCoords_[i].lengthSquare();
171        IAtom(0, 0) += mtmp * r2;
172        IAtom(1, 1) += mtmp * r2;
173        IAtom(2, 2) += mtmp * r2;
174 +      Itmp += IAtom;
175  
176        //project the inertial moment of directional atoms into this rigid body
177        if (atoms_[i]->isDirectional()) {
178 <        IAtom += atoms_[i]->getI();
179 <        Itmp += refOrients_[i].transpose() * IAtom * refOrients_[i];
181 <      } else {
182 <        Itmp += IAtom;
183 <      }
178 >        Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i];
179 >      }
180      }
181  
182 +    //    std::cout << Itmp << std::endl;
183 +
184      //diagonalize
185      Vector3d evals;
186      Mat3x3d::diagonalize(Itmp, evals, sU_);
# Line 247 | Line 245 | namespace oopse {
245          
246      }
247      
248 <    setFrc(frc);
249 <    setTrq(trq);
248 >    addFrc(frc);
249 >    addTrq(trq);
250      
251    }
252  
# Line 271 | Line 269 | namespace oopse {
269        if (atoms_[i]->isDirectional()) {
270            
271          dAtom = (DirectionalAtom *) atoms_[i];
272 <        dAtom->setA(refOrients_[i] * a);
272 >        dAtom->setA(refOrients_[i].transpose() * a);
273        }
274  
275      }
# Line 298 | Line 296 | namespace oopse {
296        if (atoms_[i]->isDirectional()) {
297            
298          dAtom = (DirectionalAtom *) atoms_[i];
299 <        dAtom->setA(refOrients_[i] * a, frame);
299 >        dAtom->setA(refOrients_[i].transpose() * a, frame);
300        }
301  
302      }
# Line 483 | Line 481 | namespace oopse {
481                 "RigidBody error.\n"
482                 "\tAtom %s does not have a position specified.\n"
483                 "\tThis means RigidBody cannot set up reference coordinates.\n",
484 <               ats->getType() );
484 >               ats->getType().c_str() );
485        painCave.isFatal = 1;
486        simError();
487      }
# Line 503 | Line 501 | namespace oopse {
501                   "RigidBody error.\n"
502                   "\tAtom %s does not have an orientation specified.\n"
503                   "\tThis means RigidBody cannot set up reference orientations.\n",
504 <                 ats->getType() );
504 >                 ats->getType().c_str() );
505          painCave.isFatal = 1;
506          simError();
507        }    

Diff Legend

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