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 374 by tim, Tue Feb 22 18:56:25 2005 UTC vs.
Revision 663 by gezelter, Wed Oct 12 21:00:59 2005 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 45 | Line 45 | namespace oopse {
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;
54      //((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
55  
56      for (int i =0 ; i < atoms_.size(); ++i){
57 <        if (atoms_[i]->isDirectional()) {
58 <            atoms_[i]->setPrevA(a * refOrients_[i]);
59 <        }
57 >      if (atoms_[i]->isDirectional()) {
58 >        atoms_[i]->setPrevA(a * refOrients_[i]);
59 >      }
60      }
61  
62 < }
62 >  }
63  
64        
65 < void RigidBody::setA(const RotMat3x3d& a) {
65 >  void RigidBody::setA(const RotMat3x3d& a) {
66      ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
67      //((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
68  
69      for (int i =0 ; i < atoms_.size(); ++i){
70 <        if (atoms_[i]->isDirectional()) {
71 <            atoms_[i]->setA(a * refOrients_[i]);
72 <        }
70 >      if (atoms_[i]->isDirectional()) {
71 >        atoms_[i]->setA(a * refOrients_[i]);
72 >      }
73      }
74 < }    
74 >  }    
75      
76 < void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
76 >  void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
77      ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
78      //((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;    
79  
80      for (int i =0 ; i < atoms_.size(); ++i){
81 <        if (atoms_[i]->isDirectional()) {
82 <            atoms_[i]->setA(a * refOrients_[i], snapshotNo);
83 <        }
81 >      if (atoms_[i]->isDirectional()) {
82 >        atoms_[i]->setA(a * refOrients_[i], snapshotNo);
83 >      }
84      }
85  
86 < }  
86 >  }  
87  
88 < Mat3x3d RigidBody::getI() {
88 >  Mat3x3d RigidBody::getI() {
89      return inertiaTensor_;
90 < }    
90 >  }    
91  
92 < std::vector<double> RigidBody::getGrad() {
93 <     std::vector<double> grad(6, 0.0);
92 >  std::vector<double> RigidBody::getGrad() {
93 >    std::vector<double> grad(6, 0.0);
94      Vector3d force;
95      Vector3d torque;
96      Vector3d myEuler;
# Line 129 | Line 129 | std::vector<double> RigidBody::getGrad() {
129  
130      //gradient is equal to -force
131      for (int j = 0 ; j<3; j++)
132 <        grad[j] = -force[j];
132 >      grad[j] = -force[j];
133  
134      for (int j = 0; j < 3; j++ ) {
135  
136 <        grad[3] += torque[j]*ephi[j];
137 <        grad[4] += torque[j]*etheta[j];
138 <        grad[5] += torque[j]*epsi[j];
136 >      grad[3] += torque[j]*ephi[j];
137 >      grad[4] += torque[j]*etheta[j];
138 >      grad[5] += torque[j]*epsi[j];
139  
140      }
141      
142      return grad;
143 < }    
143 >  }    
144  
145 < void RigidBody::accept(BaseVisitor* v) {
145 >  void RigidBody::accept(BaseVisitor* v) {
146      v->visit(this);
147 < }    
147 >  }    
148  
149 < /**@todo need modification */
150 < void  RigidBody::calcRefCoords() {
149 >  /**@todo need modification */
150 >  void  RigidBody::calcRefCoords() {
151      double mtmp;
152      Vector3d refCOM(0.0);
153      mass_ = 0.0;
154      for (std::size_t i = 0; i < atoms_.size(); ++i) {
155 <        mtmp = atoms_[i]->getMass();
156 <        mass_ += mtmp;
157 <        refCOM += refCoords_[i]*mtmp;
155 >      mtmp = atoms_[i]->getMass();
156 >      mass_ += mtmp;
157 >      refCOM += refCoords_[i]*mtmp;
158      }
159      refCOM /= mass_;
160  
161      // Next, move the origin of the reference coordinate system to the COM:
162      for (std::size_t i = 0; i < atoms_.size(); ++i) {
163 <        refCoords_[i] -= refCOM;
163 >      refCoords_[i] -= refCOM;
164      }
165  
166 < // Moment of Inertia calculation
167 <    Mat3x3d Itmp(0.0);
168 <  
166 >    // Moment of Inertia calculation
167 >    Mat3x3d Itmp(0.0);    
168      for (std::size_t i = 0; i < atoms_.size(); i++) {
169 <        mtmp = atoms_[i]->getMass();
170 <        Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
171 <        double r2 = refCoords_[i].lengthSquare();
172 <        Itmp(0, 0) += mtmp * r2;
173 <        Itmp(1, 1) += mtmp * r2;
174 <        Itmp(2, 2) += mtmp * r2;
175 <    }
169 >      Mat3x3d IAtom(0.0);  
170 >      mtmp = atoms_[i]->getMass();
171 >      IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
172 >      double r2 = refCoords_[i].lengthSquare();
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++) {
180 <        if (atoms_[i]->isDirectional()) {
181 <            RectMatrix<double, 3, 3> Iproject = refOrients_[i].transpose() * atoms_[i]->getI();
182 <            Itmp(0, 0) += Iproject(0, 0);
183 <            Itmp(1, 1) += Iproject(1, 1);
184 <            Itmp(2, 2) += Iproject(2, 2);
185 <        }
178 >      //project the inertial moment of directional atoms into this rigid body
179 >      if (atoms_[i]->isDirectional()) {
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 196 | Line 194 | void  RigidBody::calcRefCoords() {
194          
195      int nLinearAxis = 0;
196      for (int i = 0; i < 3; i++) {    
197 <        if (fabs(evals[i]) < oopse::epsilon) {
198 <            linear_ = true;
199 <            linearAxis_ = i;
200 <            ++ nLinearAxis;
201 <        }
197 >      if (fabs(evals[i]) < oopse::epsilon) {
198 >        linear_ = true;
199 >        linearAxis_ = i;
200 >        ++ nLinearAxis;
201 >      }
202      }
203  
204      if (nLinearAxis > 1) {
205 <        sprintf( painCave.errMsg,
206 <            "RigidBody error.\n"
207 <            "\tOOPSE found more than one axis in this rigid body with a vanishing \n"
208 <            "\tmoment of inertia.  This can happen in one of three ways:\n"
209 <            "\t 1) Only one atom was specified, or \n"
210 <            "\t 2) All atoms were specified at the same location, or\n"
211 <            "\t 3) The programmers did something stupid.\n"
212 <            "\tIt is silly to use a rigid body to describe this situation.  Be smarter.\n"
213 <            );
214 <        painCave.isFatal = 1;
215 <        simError();
205 >      sprintf( painCave.errMsg,
206 >               "RigidBody error.\n"
207 >               "\tOOPSE found more than one axis in this rigid body with a vanishing \n"
208 >               "\tmoment of inertia.  This can happen in one of three ways:\n"
209 >               "\t 1) Only one atom was specified, or \n"
210 >               "\t 2) All atoms were specified at the same location, or\n"
211 >               "\t 3) The programmers did something stupid.\n"
212 >               "\tIt is silly to use a rigid body to describe this situation.  Be smarter.\n"
213 >               );
214 >      painCave.isFatal = 1;
215 >      simError();
216      }
217    
218 < }
218 >  }
219  
220 < void  RigidBody::calcForcesAndTorques() {
220 >  void  RigidBody::calcForcesAndTorques() {
221      Vector3d afrc;
222      Vector3d atrq;
223      Vector3d apos;
# Line 229 | Line 227 | void  RigidBody::calcForcesAndTorques() {
227      Vector3d pos = this->getPos();
228      for (int i = 0; i < atoms_.size(); i++) {
229  
230 <        afrc = atoms_[i]->getFrc();
231 <        apos = atoms_[i]->getPos();
232 <        rpos = apos - pos;
230 >      afrc = atoms_[i]->getFrc();
231 >      apos = atoms_[i]->getPos();
232 >      rpos = apos - pos;
233          
234 <        frc += afrc;
234 >      frc += afrc;
235  
236 <        trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
237 <        trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
238 <        trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
236 >      trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
237 >      trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
238 >      trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
239  
240 <        // If the atom has a torque associated with it, then we also need to
241 <        // migrate the torques onto the center of mass:
240 >      // If the atom has a torque associated with it, then we also need to
241 >      // migrate the torques onto the center of mass:
242  
243 <        if (atoms_[i]->isDirectional()) {
244 <            atrq = atoms_[i]->getTrq();
245 <            trq += atrq;
246 <        }
243 >      if (atoms_[i]->isDirectional()) {
244 >        atrq = atoms_[i]->getTrq();
245 >        trq += atrq;
246 >      }
247          
248      }
249      
250      setFrc(frc);
251      setTrq(trq);
252      
253 < }
253 >  }
254  
255 < void  RigidBody::updateAtoms() {
255 >  void  RigidBody::updateAtoms() {
256      unsigned int i;
257      Vector3d ref;
258      Vector3d apos;
# Line 264 | Line 262 | void  RigidBody::updateAtoms() {
262      
263      for (i = 0; i < atoms_.size(); i++) {
264      
265 <        ref = body2Lab(refCoords_[i]);
265 >      ref = body2Lab(refCoords_[i]);
266  
267 <        apos = pos + ref;
267 >      apos = pos + ref;
268  
269 <        atoms_[i]->setPos(apos);
269 >      atoms_[i]->setPos(apos);
270  
271 <        if (atoms_[i]->isDirectional()) {
271 >      if (atoms_[i]->isDirectional()) {
272            
273 <          dAtom = (DirectionalAtom *) atoms_[i];
274 <          dAtom->setA(a * refOrients_[i]);
275 <          //dAtom->rotateBy( A );      
278 <        }
273 >        dAtom = (DirectionalAtom *) atoms_[i];
274 >        dAtom->setA(refOrients_[i] * a);
275 >      }
276  
277      }
278    
279 < }
279 >  }
280  
281  
282 < void  RigidBody::updateAtoms(int frame) {
282 >  void  RigidBody::updateAtoms(int frame) {
283      unsigned int i;
284      Vector3d ref;
285      Vector3d apos;
# Line 292 | Line 289 | void  RigidBody::updateAtoms(int frame) {
289      
290      for (i = 0; i < atoms_.size(); i++) {
291      
292 <        ref = body2Lab(refCoords_[i], frame);
292 >      ref = body2Lab(refCoords_[i], frame);
293  
294 <        apos = pos + ref;
294 >      apos = pos + ref;
295  
296 <        atoms_[i]->setPos(apos, frame);
296 >      atoms_[i]->setPos(apos, frame);
297  
298 <        if (atoms_[i]->isDirectional()) {
298 >      if (atoms_[i]->isDirectional()) {
299            
300 <          dAtom = (DirectionalAtom *) atoms_[i];
301 <          dAtom->setA(a * refOrients_[i], frame);
302 <        }
300 >        dAtom = (DirectionalAtom *) atoms_[i];
301 >        dAtom->setA(refOrients_[i] * a, frame);
302 >      }
303  
304      }
305    
306 < }
306 >  }
307  
308 < void RigidBody::updateAtomVel() {
308 >  void RigidBody::updateAtomVel() {
309      Mat3x3d skewMat;;
310  
311      Vector3d ji = getJ();
# Line 332 | Line 329 | void RigidBody::updateAtomVel() {
329  
330      Vector3d velRot;        
331      for (int i =0 ; i < refCoords_.size(); ++i) {
332 <        atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
332 >      atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
333      }
334  
335 < }
335 >  }
336  
337 < void RigidBody::updateAtomVel(int frame) {
337 >  void RigidBody::updateAtomVel(int frame) {
338      Mat3x3d skewMat;;
339  
340      Vector3d ji = getJ(frame);
# Line 361 | Line 358 | void RigidBody::updateAtomVel(int frame) {
358  
359      Vector3d velRot;        
360      for (int i =0 ; i < refCoords_.size(); ++i) {
361 <        atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
361 >      atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
362      }
363  
364 < }
364 >  }
365  
366          
367  
368 < bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
368 >  bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
369      if (index < atoms_.size()) {
370  
371 <        Vector3d ref = body2Lab(refCoords_[index]);
372 <        pos = getPos() + ref;
373 <        return true;
371 >      Vector3d ref = body2Lab(refCoords_[index]);
372 >      pos = getPos() + ref;
373 >      return true;
374      } else {
375 <        std::cerr << index << " is an invalid index, current rigid body contains "
376 <                      << atoms_.size() << "atoms" << std::endl;
377 <        return false;
375 >      std::cerr << index << " is an invalid index, current rigid body contains "
376 >                << atoms_.size() << "atoms" << std::endl;
377 >      return false;
378      }    
379 < }
379 >  }
380  
381 < bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) {
381 >  bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) {
382      std::vector<Atom*>::iterator i;
383      i = std::find(atoms_.begin(), atoms_.end(), atom);
384      if (i != atoms_.end()) {
385 <        //RigidBody class makes sure refCoords_ and atoms_ match each other
386 <        Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]);
387 <        pos = getPos() + ref;
388 <        return true;
385 >      //RigidBody class makes sure refCoords_ and atoms_ match each other
386 >      Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]);
387 >      pos = getPos() + ref;
388 >      return true;
389      } else {
390 <        std::cerr << "Atom " << atom->getGlobalIndex()
391 <                      <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
392 <        return false;
390 >      std::cerr << "Atom " << atom->getGlobalIndex()
391 >                <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
392 >      return false;
393      }
394 < }
395 < bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) {
394 >  }
395 >  bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) {
396  
397      //velRot = $(A\cdot skew(I^{-1}j))^{T}refCoor$
398  
399      if (index < atoms_.size()) {
400  
401 <        Vector3d velRot;
402 <        Mat3x3d skewMat;;
403 <        Vector3d ref = refCoords_[index];
404 <        Vector3d ji = getJ();
405 <        Mat3x3d I =  getI();
401 >      Vector3d velRot;
402 >      Mat3x3d skewMat;;
403 >      Vector3d ref = refCoords_[index];
404 >      Vector3d ji = getJ();
405 >      Mat3x3d I =  getI();
406  
407 <        skewMat(0, 0) =0;
408 <        skewMat(0, 1) = ji[2] /I(2, 2);
409 <        skewMat(0, 2) = -ji[1] /I(1, 1);
407 >      skewMat(0, 0) =0;
408 >      skewMat(0, 1) = ji[2] /I(2, 2);
409 >      skewMat(0, 2) = -ji[1] /I(1, 1);
410  
411 <        skewMat(1, 0) = -ji[2] /I(2, 2);
412 <        skewMat(1, 1) = 0;
413 <        skewMat(1, 2) = ji[0]/I(0, 0);
411 >      skewMat(1, 0) = -ji[2] /I(2, 2);
412 >      skewMat(1, 1) = 0;
413 >      skewMat(1, 2) = ji[0]/I(0, 0);
414  
415 <        skewMat(2, 0) =ji[1] /I(1, 1);
416 <        skewMat(2, 1) = -ji[0]/I(0, 0);
417 <        skewMat(2, 2) = 0;
415 >      skewMat(2, 0) =ji[1] /I(1, 1);
416 >      skewMat(2, 1) = -ji[0]/I(0, 0);
417 >      skewMat(2, 2) = 0;
418  
419 <        velRot = (getA() * skewMat).transpose() * ref;
419 >      velRot = (getA() * skewMat).transpose() * ref;
420  
421 <        vel =getVel() + velRot;
422 <        return true;
421 >      vel =getVel() + velRot;
422 >      return true;
423          
424      } else {
425 <        std::cerr << index << " is an invalid index, current rigid body contains "
426 <                      << atoms_.size() << "atoms" << std::endl;
427 <        return false;
425 >      std::cerr << index << " is an invalid index, current rigid body contains "
426 >                << atoms_.size() << "atoms" << std::endl;
427 >      return false;
428      }
429 < }
429 >  }
430  
431 < bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) {
431 >  bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) {
432  
433      std::vector<Atom*>::iterator i;
434      i = std::find(atoms_.begin(), atoms_.end(), atom);
435      if (i != atoms_.end()) {
436 <        return getAtomVel(vel, i - atoms_.begin());
436 >      return getAtomVel(vel, i - atoms_.begin());
437      } else {
438 <        std::cerr << "Atom " << atom->getGlobalIndex()
439 <                      <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
440 <        return false;
438 >      std::cerr << "Atom " << atom->getGlobalIndex()
439 >                <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
440 >      return false;
441      }    
442 < }
442 >  }
443  
444 < bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) {
444 >  bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) {
445      if (index < atoms_.size()) {
446  
447 <        coor = refCoords_[index];
448 <        return true;
447 >      coor = refCoords_[index];
448 >      return true;
449      } else {
450 <        std::cerr << index << " is an invalid index, current rigid body contains "
451 <                      << atoms_.size() << "atoms" << std::endl;
452 <        return false;
450 >      std::cerr << index << " is an invalid index, current rigid body contains "
451 >                << atoms_.size() << "atoms" << std::endl;
452 >      return false;
453      }
454  
455 < }
455 >  }
456  
457 < bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) {
457 >  bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) {
458      std::vector<Atom*>::iterator i;
459      i = std::find(atoms_.begin(), atoms_.end(), atom);
460      if (i != atoms_.end()) {
461 <        //RigidBody class makes sure refCoords_ and atoms_ match each other
462 <        coor = refCoords_[i - atoms_.begin()];
463 <        return true;
461 >      //RigidBody class makes sure refCoords_ and atoms_ match each other
462 >      coor = refCoords_[i - atoms_.begin()];
463 >      return true;
464      } else {
465 <        std::cerr << "Atom " << atom->getGlobalIndex()
466 <                      <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
467 <        return false;
465 >      std::cerr << "Atom " << atom->getGlobalIndex()
466 >                <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
467 >      return false;
468      }
469  
470 < }
470 >  }
471  
472  
473 < void RigidBody::addAtom(Atom* at, AtomStamp* ats) {
473 >  void RigidBody::addAtom(Atom* at, AtomStamp* ats) {
474  
475 <  Vector3d coords;
476 <  Vector3d euler;
475 >    Vector3d coords;
476 >    Vector3d euler;
477    
478  
479 <  atoms_.push_back(at);
479 >    atoms_.push_back(at);
480  
481 <  if( !ats->havePosition() ){
482 <    sprintf( painCave.errMsg,
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() );
487 <    painCave.isFatal = 1;
488 <    simError();
489 <  }
481 >    if( !ats->havePosition() ){
482 >      sprintf( painCave.errMsg,
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() );
487 >      painCave.isFatal = 1;
488 >      simError();
489 >    }
490    
491 <  coords[0] = ats->getPosX();
492 <  coords[1] = ats->getPosY();
493 <  coords[2] = ats->getPosZ();
491 >    coords[0] = ats->getPosX();
492 >    coords[1] = ats->getPosY();
493 >    coords[2] = ats->getPosZ();
494  
495 <  refCoords_.push_back(coords);
495 >    refCoords_.push_back(coords);
496  
497 <  RotMat3x3d identMat = RotMat3x3d::identity();
497 >    RotMat3x3d identMat = RotMat3x3d::identity();
498    
499 <  if (at->isDirectional()) {  
499 >    if (at->isDirectional()) {  
500  
501 <    if( !ats->haveOrientation() ){
502 <      sprintf( painCave.errMsg,
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() );
507 <      painCave.isFatal = 1;
508 <      simError();
509 <    }    
501 >      if( !ats->haveOrientation() ){
502 >        sprintf( painCave.errMsg,
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() );
507 >        painCave.isFatal = 1;
508 >        simError();
509 >      }    
510      
511 <    euler[0] = ats->getEulerPhi() * NumericConstant::PI /180.0;
512 <    euler[1] = ats->getEulerTheta() * NumericConstant::PI /180.0;
513 <    euler[2] = ats->getEulerPsi() * NumericConstant::PI /180.0;
511 >      euler[0] = ats->getEulerPhi() * NumericConstant::PI /180.0;
512 >      euler[1] = ats->getEulerTheta() * NumericConstant::PI /180.0;
513 >      euler[2] = ats->getEulerPsi() * NumericConstant::PI /180.0;
514  
515 <    RotMat3x3d Atmp(euler);
516 <    refOrients_.push_back(Atmp);
515 >      RotMat3x3d Atmp(euler);
516 >      refOrients_.push_back(Atmp);
517      
518 <  }else {
519 <    refOrients_.push_back(identMat);
520 <  }
518 >    }else {
519 >      refOrients_.push_back(identMat);
520 >    }
521    
522    
523 < }
523 >  }
524  
525   }
526  

Diff Legend

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