ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/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 1360 by cli2, Mon Sep 7 16:31:51 2009 UTC

# Line 44 | Line 44 | namespace oopse {
44   #include "utils/simError.h"
45   #include "utils/NumericConstant.hpp"
46   namespace oopse {
47 <
48 <  RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){
49 <
47 >  
48 >  RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData),
49 >                           inertiaTensor_(0.0){    
50    }
51 <
51 >  
52    void RigidBody::setPrevA(const RotMat3x3d& a) {
53      ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
54 <    //((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
55 <
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 <
60 >    
61    }
62 <
63 <      
62 >  
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    }    
73 <    
73 >  
74    void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
75      ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
76 +    
77      //((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;    
78 <
78 >    
79      for (int i =0 ; i < atoms_.size(); ++i){
80        if (atoms_[i]->isDirectional()) {
81 <        atoms_[i]->setA(a * refOrients_[i], snapshotNo);
81 >        atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo);
82        }
83      }
84 <
84 >    
85    }  
86 <
86 >  
87    Mat3x3d RigidBody::getI() {
88      return inertiaTensor_;
89    }    
90 <
91 <  std::vector<double> RigidBody::getGrad() {
92 <    std::vector<double> grad(6, 0.0);
90 >  
91 >  std::vector<RealType> RigidBody::getGrad() {
92 >    std::vector<RealType> grad(6, 0.0);
93      Vector3d force;
94      Vector3d torque;
95      Vector3d myEuler;
96 <    double phi, theta, psi;
97 <    double cphi, sphi, ctheta, stheta;
96 >    RealType phi, theta, psi;
97 >    RealType cphi, sphi, ctheta, stheta;
98      Vector3d ephi;
99      Vector3d etheta;
100      Vector3d epsi;
101 <
101 >    
102      force = getFrc();
103      torque =getTrq();
104      myEuler = getA().toEulerAngles();
105 <
105 >    
106      phi = myEuler[0];
107      theta = myEuler[1];
108      psi = myEuler[2];
109 <
109 >    
110      cphi = cos(phi);
111      sphi = sin(phi);
112      ctheta = cos(theta);
113      stheta = sin(theta);
114 <
114 >    
115      // get unit vectors along the phi, theta and psi rotation axes
116 <
116 >    
117      ephi[0] = 0.0;
118      ephi[1] = 0.0;
119      ephi[2] = 1.0;
120 <
121 <    etheta[0] = cphi;
122 <    etheta[1] = sphi;
123 <    etheta[2] = 0.0;
124 <
120 >    
121 >    etheta[0] = -sphi;
122 >    etheta[1] =  cphi;
123 >    etheta[2] =  0.0;
124 >    
125      epsi[0] = stheta * cphi;
126      epsi[1] = stheta * sphi;
127      epsi[2] = ctheta;
128 <
128 >    
129      //gradient is equal to -force
130      for (int j = 0 ; j<3; j++)
131        grad[j] = -force[j];
132 <
132 >    
133      for (int j = 0; j < 3; j++ ) {
134 <
134 >      
135        grad[3] += torque[j]*ephi[j];
136        grad[4] += torque[j]*etheta[j];
137        grad[5] += torque[j]*epsi[j];
138 <
138 >      
139      }
140      
141      return grad;
142    }    
143 <
143 >  
144    void RigidBody::accept(BaseVisitor* v) {
145      v->visit(this);
146    }    
147  
148    /**@todo need modification */
149    void  RigidBody::calcRefCoords() {
150 <    double mtmp;
150 >    RealType mtmp;
151      Vector3d refCOM(0.0);
152      mass_ = 0.0;
153      for (std::size_t i = 0; i < atoms_.size(); ++i) {
# Line 157 | Line 156 | namespace oopse {
156        refCOM += refCoords_[i]*mtmp;
157      }
158      refCOM /= mass_;
159 <
159 >    
160      // Next, move the origin of the reference coordinate system to the COM:
161      for (std::size_t i = 0; i < atoms_.size(); ++i) {
162        refCoords_[i] -= refCOM;
# Line 169 | Line 168 | namespace oopse {
168        Mat3x3d IAtom(0.0);  
169        mtmp = atoms_[i]->getMass();
170        IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
171 <      double r2 = refCoords_[i].lengthSquare();
171 >      RealType r2 = refCoords_[i].lengthSquare();
172        IAtom(0, 0) += mtmp * r2;
173        IAtom(1, 1) += mtmp * r2;
174        IAtom(2, 2) += mtmp * r2;
175 <
175 >      Itmp += IAtom;
176 >      
177        //project the inertial moment of directional atoms into this rigid body
178        if (atoms_[i]->isDirectional()) {
179 <        IAtom += atoms_[i]->getI();
180 <        Itmp += refOrients_[i].transpose() * IAtom * refOrients_[i];
181 <      } else {
182 <        Itmp += IAtom;
183 <      }
179 >        Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i];
180 >      }
181      }
182  
183 +    //    std::cout << Itmp << std::endl;
184 +
185      //diagonalize
186      Vector3d evals;
187      Mat3x3d::diagonalize(Itmp, evals, sU_);
# Line 223 | Line 222 | namespace oopse {
222      Vector3d apos;
223      Vector3d rpos;
224      Vector3d frc(0.0);
225 <    Vector3d trq(0.0);
225 >    Vector3d trq(0.0);    
226      Vector3d pos = this->getPos();
227      for (int i = 0; i < atoms_.size(); i++) {
228  
# Line 243 | Line 242 | namespace oopse {
242        if (atoms_[i]->isDirectional()) {
243          atrq = atoms_[i]->getTrq();
244          trq += atrq;
245 <      }
245 >      }      
246 >    }        
247 >    addFrc(frc);
248 >    addTrq(trq);    
249 >  }
250 >
251 >  Mat3x3d RigidBody::calcForcesAndTorquesAndVirial() {
252 >    Vector3d afrc;
253 >    Vector3d atrq;
254 >    Vector3d apos;
255 >    Vector3d rpos;
256 >    Vector3d dfrc;
257 >    Vector3d frc(0.0);
258 >    Vector3d trq(0.0);    
259 >    Vector3d pos = this->getPos();
260 >    Mat3x3d tau_(0.0);
261 >
262 >    for (int i = 0; i < atoms_.size(); i++) {
263 >      
264 >      afrc = atoms_[i]->getFrc();
265 >      apos = atoms_[i]->getPos();
266 >      rpos = apos - pos;
267          
268 +      frc += afrc;
269 +
270 +      trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
271 +      trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
272 +      trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
273 +
274 +      // If the atom has a torque associated with it, then we also need to
275 +      // migrate the torques onto the center of mass:
276 +
277 +      if (atoms_[i]->isDirectional()) {
278 +        atrq = atoms_[i]->getTrq();
279 +        trq += atrq;
280 +      }
281 +      
282 +      tau_(0,0) -= rpos[0]*afrc[0];
283 +      tau_(0,1) -= rpos[0]*afrc[1];
284 +      tau_(0,2) -= rpos[0]*afrc[2];
285 +      tau_(1,0) -= rpos[1]*afrc[0];
286 +      tau_(1,1) -= rpos[1]*afrc[1];
287 +      tau_(1,2) -= rpos[1]*afrc[2];
288 +      tau_(2,0) -= rpos[2]*afrc[0];
289 +      tau_(2,1) -= rpos[2]*afrc[1];
290 +      tau_(2,2) -= rpos[2]*afrc[2];
291 +
292      }
293 <    
294 <    setFrc(frc);
295 <    setTrq(trq);
252 <    
293 >    addFrc(frc);
294 >    addTrq(trq);
295 >    return tau_;
296    }
297  
298    void  RigidBody::updateAtoms() {
# Line 271 | Line 314 | namespace oopse {
314        if (atoms_[i]->isDirectional()) {
315            
316          dAtom = (DirectionalAtom *) atoms_[i];
317 <        dAtom->setA(refOrients_[i] * a);
317 >        dAtom->setA(refOrients_[i].transpose() * a);
318        }
319  
320      }
# Line 298 | Line 341 | namespace oopse {
341        if (atoms_[i]->isDirectional()) {
342            
343          dAtom = (DirectionalAtom *) atoms_[i];
344 <        dAtom->setA(refOrients_[i] * a, frame);
344 >        dAtom->setA(refOrients_[i].transpose() * a, frame);
345        }
346  
347      }
# Line 483 | Line 526 | namespace oopse {
526                 "RigidBody error.\n"
527                 "\tAtom %s does not have a position specified.\n"
528                 "\tThis means RigidBody cannot set up reference coordinates.\n",
529 <               ats->getType() );
529 >               ats->getType().c_str() );
530        painCave.isFatal = 1;
531        simError();
532      }
# Line 503 | Line 546 | namespace oopse {
546                   "RigidBody error.\n"
547                   "\tAtom %s does not have an orientation specified.\n"
548                   "\tThis means RigidBody cannot set up reference orientations.\n",
549 <                 ats->getType() );
549 >                 ats->getType().c_str() );
550          painCave.isFatal = 1;
551          simError();
552        }    

Diff Legend

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