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 334 by tim, Mon Feb 14 17:57:01 2005 UTC vs.
Revision 1126 by gezelter, Fri Apr 6 21:53:43 2007 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 42 | Line 42
42   #include <math.h>
43   #include "primitives/RigidBody.hpp"
44   #include "utils/simError.h"
45 + #include "utils/NumericConstant.hpp"
46   namespace oopse {
47  
48 < RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){
48 >  RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){
49  
50 < }
50 >  }
51  
52 < void RigidBody::setPrevA(const RotMat3x3d& a) {
52 >  void RigidBody::setPrevA(const RotMat3x3d& a) {
53      ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
53 –    //((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]);
58 <        }
56 >      if (atoms_[i]->isDirectional()) {
57 >        atoms_[i]->setPrevA(refOrients_[i].transpose() * a);
58 >      }
59      }
60  
61 < }
61 >  }
62  
63        
64 < void RigidBody::setA(const RotMat3x3d& a) {
64 >  void RigidBody::setA(const RotMat3x3d& a) {
65      ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
66 –    //((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]);
70 <        }
68 >      if (atoms_[i]->isDirectional()) {
69 >        atoms_[i]->setA(refOrients_[i].transpose() * a);
70 >      }
71      }
72 < }    
72 >  }    
73      
74 < void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
74 >  void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
75      ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
76      //((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;    
77  
78      for (int i =0 ; i < atoms_.size(); ++i){
79 <        if (atoms_[i]->isDirectional()) {
80 <            atoms_[i]->setA(a * refOrients_[i], snapshotNo);
81 <        }
79 >      if (atoms_[i]->isDirectional()) {
80 >        atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo);
81 >      }
82      }
83  
84 < }  
84 >  }  
85  
86 < Mat3x3d RigidBody::getI() {
86 >  Mat3x3d RigidBody::getI() {
87      return inertiaTensor_;
88 < }    
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 128 | Line 127 | std::vector<double> RigidBody::getGrad() {
127  
128      //gradient is equal to -force
129      for (int j = 0 ; j<3; j++)
130 <        grad[j] = -force[j];
130 >      grad[j] = -force[j];
131  
132      for (int j = 0; j < 3; j++ ) {
133  
134 <        grad[3] += torque[j]*ephi[j];
135 <        grad[4] += torque[j]*etheta[j];
136 <        grad[5] += torque[j]*epsi[j];
134 >      grad[3] += torque[j]*ephi[j];
135 >      grad[4] += torque[j]*etheta[j];
136 >      grad[5] += torque[j]*epsi[j];
137  
138      }
139      
140      return grad;
141 < }    
141 >  }    
142  
143 < void RigidBody::accept(BaseVisitor* v) {
143 >  void RigidBody::accept(BaseVisitor* v) {
144      v->visit(this);
145 < }    
145 >  }    
146  
147 < /**@todo need modification */
148 < void  RigidBody::calcRefCoords() {
149 <    double mtmp;
147 >  /**@todo need modification */
148 >  void  RigidBody::calcRefCoords() {
149 >    RealType mtmp;
150      Vector3d refCOM(0.0);
151      mass_ = 0.0;
152      for (std::size_t i = 0; i < atoms_.size(); ++i) {
153 <        mtmp = atoms_[i]->getMass();
154 <        mass_ += mtmp;
155 <        refCOM += refCoords_[i]*mtmp;
153 >      mtmp = atoms_[i]->getMass();
154 >      mass_ += mtmp;
155 >      refCOM += refCoords_[i]*mtmp;
156      }
157      refCOM /= mass_;
158  
159      // Next, move the origin of the reference coordinate system to the COM:
160      for (std::size_t i = 0; i < atoms_.size(); ++i) {
161 <        refCoords_[i] -= refCOM;
161 >      refCoords_[i] -= refCOM;
162      }
163  
164 < // Moment of Inertia calculation
165 <    Mat3x3d Itmp(0.0);
167 <  
164 >    // Moment of Inertia calculation
165 >    Mat3x3d Itmp(0.0);    
166      for (std::size_t i = 0; i < atoms_.size(); i++) {
167 <        mtmp = atoms_[i]->getMass();
168 <        Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
169 <        double r2 = refCoords_[i].lengthSquare();
170 <        Itmp(0, 0) += mtmp * r2;
171 <        Itmp(1, 1) += mtmp * r2;
172 <        Itmp(2, 2) += mtmp * r2;
173 <    }
167 >      Mat3x3d IAtom(0.0);  
168 >      mtmp = atoms_[i]->getMass();
169 >      IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
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 <    for (std::size_t i = 0; i < atoms_.size(); i++) {
178 <        if (atoms_[i]->isDirectional()) {
179 <            RectMatrix<double, 3, 3> Iproject = refOrients_[i].transpose() * atoms_[i]->getI();
181 <            Itmp(0, 0) += Iproject(0, 0);
182 <            Itmp(1, 1) += Iproject(1, 1);
183 <            Itmp(2, 2) += Iproject(2, 2);
184 <        }
176 >      //project the inertial moment of directional atoms into this rigid body
177 >      if (atoms_[i]->isDirectional()) {
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 195 | Line 192 | void  RigidBody::calcRefCoords() {
192          
193      int nLinearAxis = 0;
194      for (int i = 0; i < 3; i++) {    
195 <        if (fabs(evals[i]) < oopse::epsilon) {
196 <            linear_ = true;
197 <            linearAxis_ = i;
198 <            ++ nLinearAxis;
199 <        }
195 >      if (fabs(evals[i]) < oopse::epsilon) {
196 >        linear_ = true;
197 >        linearAxis_ = i;
198 >        ++ nLinearAxis;
199 >      }
200      }
201  
202      if (nLinearAxis > 1) {
203 <        sprintf( painCave.errMsg,
204 <            "RigidBody error.\n"
205 <            "\tOOPSE found more than one axis in this rigid body with a vanishing \n"
206 <            "\tmoment of inertia.  This can happen in one of three ways:\n"
207 <            "\t 1) Only one atom was specified, or \n"
208 <            "\t 2) All atoms were specified at the same location, or\n"
209 <            "\t 3) The programmers did something stupid.\n"
210 <            "\tIt is silly to use a rigid body to describe this situation.  Be smarter.\n"
211 <            );
212 <        painCave.isFatal = 1;
213 <        simError();
203 >      sprintf( painCave.errMsg,
204 >               "RigidBody error.\n"
205 >               "\tOOPSE found more than one axis in this rigid body with a vanishing \n"
206 >               "\tmoment of inertia.  This can happen in one of three ways:\n"
207 >               "\t 1) Only one atom was specified, or \n"
208 >               "\t 2) All atoms were specified at the same location, or\n"
209 >               "\t 3) The programmers did something stupid.\n"
210 >               "\tIt is silly to use a rigid body to describe this situation.  Be smarter.\n"
211 >               );
212 >      painCave.isFatal = 1;
213 >      simError();
214      }
215    
216 < }
216 >  }
217  
218 < void  RigidBody::calcForcesAndTorques() {
218 >  void  RigidBody::calcForcesAndTorques() {
219      Vector3d afrc;
220      Vector3d atrq;
221      Vector3d apos;
222      Vector3d rpos;
223      Vector3d frc(0.0);
224 <    Vector3d trq(0.0);
224 >    Vector3d trq(0.0);    
225      Vector3d pos = this->getPos();
226      for (int i = 0; i < atoms_.size(); i++) {
227  
228 <        afrc = atoms_[i]->getFrc();
229 <        apos = atoms_[i]->getPos();
230 <        rpos = apos - pos;
228 >      afrc = atoms_[i]->getFrc();
229 >      apos = atoms_[i]->getPos();
230 >      rpos = apos - pos;
231          
232 <        frc += afrc;
232 >      frc += afrc;
233  
234 <        trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
235 <        trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
236 <        trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
234 >      trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
235 >      trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
236 >      trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
237  
238 <        // If the atom has a torque associated with it, then we also need to
239 <        // migrate the torques onto the center of mass:
238 >      // If the atom has a torque associated with it, then we also need to
239 >      // migrate the torques onto the center of mass:
240  
241 <        if (atoms_[i]->isDirectional()) {
242 <            atrq = atoms_[i]->getTrq();
243 <            trq += atrq;
244 <        }
241 >      if (atoms_[i]->isDirectional()) {
242 >        atrq = atoms_[i]->getTrq();
243 >        trq += atrq;
244 >      }      
245 >    }        
246 >    addFrc(frc);
247 >    addTrq(trq);    
248 >  }
249 >
250 >  Mat3x3d RigidBody::calcForcesAndTorquesAndVirial() {
251 >    Vector3d afrc;
252 >    Vector3d atrq;
253 >    Vector3d apos;
254 >    Vector3d rpos;
255 >    Vector3d frc(0.0);
256 >    Vector3d trq(0.0);    
257 >    Vector3d pos = this->getPos();
258 >    Mat3x3d tau_(0.0);
259 >
260 >    for (int i = 0; i < atoms_.size(); i++) {
261 >
262 >      afrc = atoms_[i]->getFrc();
263 >      apos = atoms_[i]->getPos();
264 >      rpos = apos - pos;
265          
266 <    }
250 <    
251 <    setFrc(frc);
252 <    setTrq(trq);
253 <    
254 < }
266 >      frc += afrc;
267  
268 < void  RigidBody::updateAtoms() {
268 >      trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
269 >      trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
270 >      trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
271 >
272 >      // If the atom has a torque associated with it, then we also need to
273 >      // migrate the torques onto the center of mass:
274 >
275 >      if (atoms_[i]->isDirectional()) {
276 >        atrq = atoms_[i]->getTrq();
277 >        trq += atrq;
278 >      }
279 >      
280 >      tau_(0,0) -= rpos[0]*afrc[0];
281 >      tau_(0,1) -= rpos[0]*afrc[1];
282 >      tau_(0,2) -= rpos[0]*afrc[2];
283 >      tau_(1,0) -= rpos[1]*afrc[0];
284 >      tau_(1,1) -= rpos[1]*afrc[1];
285 >      tau_(1,2) -= rpos[1]*afrc[2];
286 >      tau_(2,0) -= rpos[2]*afrc[0];
287 >      tau_(2,1) -= rpos[2]*afrc[1];
288 >      tau_(2,2) -= rpos[2]*afrc[2];
289 >      
290 >    }        
291 >    addFrc(frc);
292 >    addTrq(trq);
293 >    return tau_;
294 >  }
295 >
296 >  void  RigidBody::updateAtoms() {
297      unsigned int i;
298      Vector3d ref;
299      Vector3d apos;
# Line 263 | Line 303 | void  RigidBody::updateAtoms() {
303      
304      for (i = 0; i < atoms_.size(); i++) {
305      
306 <        ref = body2Lab(refCoords_[i]);
306 >      ref = body2Lab(refCoords_[i]);
307  
308 <        apos = pos + ref;
308 >      apos = pos + ref;
309  
310 <        atoms_[i]->setPos(apos);
310 >      atoms_[i]->setPos(apos);
311  
312 <        if (atoms_[i]->isDirectional()) {
312 >      if (atoms_[i]->isDirectional()) {
313            
314 <          dAtom = (DirectionalAtom *) atoms_[i];
315 <          dAtom->setA(a * refOrients_[i]);
316 <          //dAtom->rotateBy( A );      
277 <        }
314 >        dAtom = (DirectionalAtom *) atoms_[i];
315 >        dAtom->setA(refOrients_[i].transpose() * a);
316 >      }
317  
318      }
319    
320 < }
320 >  }
321  
322  
323 < void  RigidBody::updateAtoms(int frame) {
323 >  void  RigidBody::updateAtoms(int frame) {
324      unsigned int i;
325      Vector3d ref;
326      Vector3d apos;
# Line 291 | Line 330 | void  RigidBody::updateAtoms(int frame) {
330      
331      for (i = 0; i < atoms_.size(); i++) {
332      
333 <        ref = body2Lab(refCoords_[i], frame);
333 >      ref = body2Lab(refCoords_[i], frame);
334  
335 <        apos = pos + ref;
335 >      apos = pos + ref;
336  
337 <        atoms_[i]->setPos(apos, frame);
337 >      atoms_[i]->setPos(apos, frame);
338  
339 <        if (atoms_[i]->isDirectional()) {
339 >      if (atoms_[i]->isDirectional()) {
340            
341 <          dAtom = (DirectionalAtom *) atoms_[i];
342 <          dAtom->setA(a * refOrients_[i], frame);
343 <        }
341 >        dAtom = (DirectionalAtom *) atoms_[i];
342 >        dAtom->setA(refOrients_[i].transpose() * a, frame);
343 >      }
344  
345      }
346    
347 < }
347 >  }
348  
349 < void RigidBody::updateAtomVel() {
349 >  void RigidBody::updateAtomVel() {
350      Mat3x3d skewMat;;
351  
352      Vector3d ji = getJ();
# Line 331 | Line 370 | void RigidBody::updateAtomVel() {
370  
371      Vector3d velRot;        
372      for (int i =0 ; i < refCoords_.size(); ++i) {
373 <        atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
373 >      atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
374      }
375  
376 < }
376 >  }
377  
378 < void RigidBody::updateAtomVel(int frame) {
378 >  void RigidBody::updateAtomVel(int frame) {
379      Mat3x3d skewMat;;
380  
381      Vector3d ji = getJ(frame);
# Line 360 | Line 399 | void RigidBody::updateAtomVel(int frame) {
399  
400      Vector3d velRot;        
401      for (int i =0 ; i < refCoords_.size(); ++i) {
402 <        atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
402 >      atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
403      }
404  
405 < }
405 >  }
406  
407          
408  
409 < bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
409 >  bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
410      if (index < atoms_.size()) {
411  
412 <        Vector3d ref = body2Lab(refCoords_[index]);
413 <        pos = getPos() + ref;
414 <        return true;
412 >      Vector3d ref = body2Lab(refCoords_[index]);
413 >      pos = getPos() + ref;
414 >      return true;
415      } else {
416 <        std::cerr << index << " is an invalid index, current rigid body contains "
417 <                      << atoms_.size() << "atoms" << std::endl;
418 <        return false;
416 >      std::cerr << index << " is an invalid index, current rigid body contains "
417 >                << atoms_.size() << "atoms" << std::endl;
418 >      return false;
419      }    
420 < }
420 >  }
421  
422 < bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) {
422 >  bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) {
423      std::vector<Atom*>::iterator i;
424      i = std::find(atoms_.begin(), atoms_.end(), atom);
425      if (i != atoms_.end()) {
426 <        //RigidBody class makes sure refCoords_ and atoms_ match each other
427 <        Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]);
428 <        pos = getPos() + ref;
429 <        return true;
426 >      //RigidBody class makes sure refCoords_ and atoms_ match each other
427 >      Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]);
428 >      pos = getPos() + ref;
429 >      return true;
430      } else {
431 <        std::cerr << "Atom " << atom->getGlobalIndex()
432 <                      <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
433 <        return false;
431 >      std::cerr << "Atom " << atom->getGlobalIndex()
432 >                <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
433 >      return false;
434      }
435 < }
436 < bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) {
435 >  }
436 >  bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) {
437  
438      //velRot = $(A\cdot skew(I^{-1}j))^{T}refCoor$
439  
440      if (index < atoms_.size()) {
441  
442 <        Vector3d velRot;
443 <        Mat3x3d skewMat;;
444 <        Vector3d ref = refCoords_[index];
445 <        Vector3d ji = getJ();
446 <        Mat3x3d I =  getI();
442 >      Vector3d velRot;
443 >      Mat3x3d skewMat;;
444 >      Vector3d ref = refCoords_[index];
445 >      Vector3d ji = getJ();
446 >      Mat3x3d I =  getI();
447  
448 <        skewMat(0, 0) =0;
449 <        skewMat(0, 1) = ji[2] /I(2, 2);
450 <        skewMat(0, 2) = -ji[1] /I(1, 1);
448 >      skewMat(0, 0) =0;
449 >      skewMat(0, 1) = ji[2] /I(2, 2);
450 >      skewMat(0, 2) = -ji[1] /I(1, 1);
451  
452 <        skewMat(1, 0) = -ji[2] /I(2, 2);
453 <        skewMat(1, 1) = 0;
454 <        skewMat(1, 2) = ji[0]/I(0, 0);
452 >      skewMat(1, 0) = -ji[2] /I(2, 2);
453 >      skewMat(1, 1) = 0;
454 >      skewMat(1, 2) = ji[0]/I(0, 0);
455  
456 <        skewMat(2, 0) =ji[1] /I(1, 1);
457 <        skewMat(2, 1) = -ji[0]/I(0, 0);
458 <        skewMat(2, 2) = 0;
456 >      skewMat(2, 0) =ji[1] /I(1, 1);
457 >      skewMat(2, 1) = -ji[0]/I(0, 0);
458 >      skewMat(2, 2) = 0;
459  
460 <        velRot = (getA() * skewMat).transpose() * ref;
460 >      velRot = (getA() * skewMat).transpose() * ref;
461  
462 <        vel =getVel() + velRot;
463 <        return true;
462 >      vel =getVel() + velRot;
463 >      return true;
464          
465      } else {
466 <        std::cerr << index << " is an invalid index, current rigid body contains "
467 <                      << atoms_.size() << "atoms" << std::endl;
468 <        return false;
466 >      std::cerr << index << " is an invalid index, current rigid body contains "
467 >                << atoms_.size() << "atoms" << std::endl;
468 >      return false;
469      }
470 < }
470 >  }
471  
472 < bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) {
472 >  bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) {
473  
474      std::vector<Atom*>::iterator i;
475      i = std::find(atoms_.begin(), atoms_.end(), atom);
476      if (i != atoms_.end()) {
477 <        return getAtomVel(vel, i - atoms_.begin());
477 >      return getAtomVel(vel, i - atoms_.begin());
478      } else {
479 <        std::cerr << "Atom " << atom->getGlobalIndex()
480 <                      <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
481 <        return false;
479 >      std::cerr << "Atom " << atom->getGlobalIndex()
480 >                <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
481 >      return false;
482      }    
483 < }
483 >  }
484  
485 < bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) {
485 >  bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) {
486      if (index < atoms_.size()) {
487  
488 <        coor = refCoords_[index];
489 <        return true;
488 >      coor = refCoords_[index];
489 >      return true;
490      } else {
491 <        std::cerr << index << " is an invalid index, current rigid body contains "
492 <                      << atoms_.size() << "atoms" << std::endl;
493 <        return false;
491 >      std::cerr << index << " is an invalid index, current rigid body contains "
492 >                << atoms_.size() << "atoms" << std::endl;
493 >      return false;
494      }
495  
496 < }
496 >  }
497  
498 < bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) {
498 >  bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) {
499      std::vector<Atom*>::iterator i;
500      i = std::find(atoms_.begin(), atoms_.end(), atom);
501      if (i != atoms_.end()) {
502 <        //RigidBody class makes sure refCoords_ and atoms_ match each other
503 <        coor = refCoords_[i - atoms_.begin()];
504 <        return true;
502 >      //RigidBody class makes sure refCoords_ and atoms_ match each other
503 >      coor = refCoords_[i - atoms_.begin()];
504 >      return true;
505      } else {
506 <        std::cerr << "Atom " << atom->getGlobalIndex()
507 <                      <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
508 <        return false;
506 >      std::cerr << "Atom " << atom->getGlobalIndex()
507 >                <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
508 >      return false;
509      }
510  
511 < }
511 >  }
512  
513  
514 < void RigidBody::addAtom(Atom* at, AtomStamp* ats) {
514 >  void RigidBody::addAtom(Atom* at, AtomStamp* ats) {
515  
516 <  Vector3d coords;
517 <  Vector3d euler;
516 >    Vector3d coords;
517 >    Vector3d euler;
518    
519  
520 <  atoms_.push_back(at);
520 >    atoms_.push_back(at);
521  
522 <  if( !ats->havePosition() ){
523 <    sprintf( painCave.errMsg,
524 <             "RigidBody error.\n"
525 <             "\tAtom %s does not have a position specified.\n"
526 <             "\tThis means RigidBody cannot set up reference coordinates.\n",
527 <             ats->getType() );
528 <    painCave.isFatal = 1;
529 <    simError();
530 <  }
522 >    if( !ats->havePosition() ){
523 >      sprintf( painCave.errMsg,
524 >               "RigidBody error.\n"
525 >               "\tAtom %s does not have a position specified.\n"
526 >               "\tThis means RigidBody cannot set up reference coordinates.\n",
527 >               ats->getType().c_str() );
528 >      painCave.isFatal = 1;
529 >      simError();
530 >    }
531    
532 <  coords[0] = ats->getPosX();
533 <  coords[1] = ats->getPosY();
534 <  coords[2] = ats->getPosZ();
532 >    coords[0] = ats->getPosX();
533 >    coords[1] = ats->getPosY();
534 >    coords[2] = ats->getPosZ();
535  
536 <  refCoords_.push_back(coords);
536 >    refCoords_.push_back(coords);
537  
538 <  RotMat3x3d identMat = RotMat3x3d::identity();
538 >    RotMat3x3d identMat = RotMat3x3d::identity();
539    
540 <  if (at->isDirectional()) {  
540 >    if (at->isDirectional()) {  
541  
542 <    if( !ats->haveOrientation() ){
543 <      sprintf( painCave.errMsg,
544 <               "RigidBody error.\n"
545 <               "\tAtom %s does not have an orientation specified.\n"
546 <               "\tThis means RigidBody cannot set up reference orientations.\n",
547 <               ats->getType() );
548 <      painCave.isFatal = 1;
549 <      simError();
550 <    }    
542 >      if( !ats->haveOrientation() ){
543 >        sprintf( painCave.errMsg,
544 >                 "RigidBody error.\n"
545 >                 "\tAtom %s does not have an orientation specified.\n"
546 >                 "\tThis means RigidBody cannot set up reference orientations.\n",
547 >                 ats->getType().c_str() );
548 >        painCave.isFatal = 1;
549 >        simError();
550 >      }    
551      
552 <    euler[0] = ats->getEulerPhi();
553 <    euler[1] = ats->getEulerTheta();
554 <    euler[2] = ats->getEulerPsi();
552 >      euler[0] = ats->getEulerPhi() * NumericConstant::PI /180.0;
553 >      euler[1] = ats->getEulerTheta() * NumericConstant::PI /180.0;
554 >      euler[2] = ats->getEulerPsi() * NumericConstant::PI /180.0;
555  
556 <    RotMat3x3d Atmp(euler);
557 <    refOrients_.push_back(Atmp);
556 >      RotMat3x3d Atmp(euler);
557 >      refOrients_.push_back(Atmp);
558      
559 <  }else {
560 <    refOrients_.push_back(identMat);
561 <  }
559 >    }else {
560 >      refOrients_.push_back(identMat);
561 >    }
562    
563    
564 < }
564 >  }
565  
566   }
567  

Diff Legend

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