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 1126 by gezelter, Fri Apr 6 21:53:43 2007 UTC vs.
Revision 1211 by gezelter, Wed Jan 23 16:38:22 2008 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 <
54 >    
55      for (int i =0 ; i < atoms_.size(); ++i){
56        if (atoms_[i]->isDirectional()) {
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;
66  
# Line 70 | Line 70 | namespace oopse {
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(refOrients_[i].transpose() * a, snapshotNo);
82        }
83      }
84 <
84 >    
85    }  
86 <
86 >  
87    Mat3x3d RigidBody::getI() {
88      return inertiaTensor_;
89    }    
90 <
90 >  
91    std::vector<RealType> RigidBody::getGrad() {
92      std::vector<RealType> grad(6, 0.0);
93      Vector3d force;
# Line 97 | Line 98 | namespace oopse {
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 <
120 >    
121      etheta[0] = cphi;
122      etheta[1] = sphi;
123      etheta[2] = 0.0;
124 <
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    }    
# Line 155 | 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 172 | Line 173 | namespace oopse {
173        IAtom(1, 1) += mtmp * r2;
174        IAtom(2, 2) += mtmp * r2;
175        Itmp += IAtom;
176 <
176 >      
177        //project the inertial moment of directional atoms into this rigid body
178        if (atoms_[i]->isDirectional()) {
179          Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i];
# Line 252 | Line 253 | namespace oopse {
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 <
263 >      
264        afrc = atoms_[i]->getFrc();
265        apos = atoms_[i]->getPos();
266        rpos = apos - pos;
# Line 286 | Line 288 | namespace oopse {
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 <    }        
291 >
292 >    }
293      addFrc(frc);
294      addTrq(trq);
295      return tau_;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines