ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/RigidBody.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/RigidBody.cpp (file contents):
Revision 1100 by gezelter, Mon Apr 12 21:02:01 2004 UTC vs.
Revision 1254 by tim, Wed Jun 9 16:16:33 2004 UTC

# Line 6 | Line 6 | RigidBody::RigidBody() : StuntDouble() {
6  
7   RigidBody::RigidBody() : StuntDouble() {
8    objType = OT_RIGIDBODY;
9 <  com_good = false;
10 <  precalc_done = false;
9 >  is_linear = false;
10 >  linear_axis =  -1;
11 >  momIntTol = 1e-6;
12   }
13  
14   RigidBody::~RigidBody() {
# Line 97 | Line 98 | void RigidBody::zeroForces() {
98      trq[i] = 0.0;
99    }
100  
100 –  forces_good = false;
101 –
101   }
102  
103 < void RigidBody::setEulerAngles( double phi, double theta, double psi ){
103 > void RigidBody::setEuler( double phi, double theta, double psi ){
104    
105      A[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
106      A[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
# Line 191 | Line 190 | void RigidBody::getA( double the_A[3][3] ){
190    
191    for (int i = 0; i < 3; i++)
192      for (int j = 0; j < 3; j++)
193 <      the_A[i][j] = the_A[i][j];
193 >      the_A[i][j] = A[i][j];
194  
195   }
196  
# Line 228 | Line 227 | void RigidBody::getI( double the_I[3][3] ){  
227   }      
228  
229   void RigidBody::getI( double the_I[3][3] ){  
231 –
232 –  if (precalc_done) {
230  
231      for (int i = 0; i < 3; i++)
232        for (int j = 0; j < 3; j++)
233          the_I[i][j] = I[i][j];
234  
238 –  } else {
239 –
240 –  }
235   }
236  
237   void RigidBody::lab2Body( double r[3] ){
# Line 265 | Line 259 | void RigidBody::body2Lab( double r[3] ){
259    r[0] = (A[0][0] * rb[0]) + (A[1][0] * rb[1]) + (A[2][0] * rb[2]);
260    r[1] = (A[0][1] * rb[0]) + (A[1][1] * rb[1]) + (A[2][1] * rb[2]);
261    r[2] = (A[0][2] * rb[0]) + (A[1][2] * rb[1]) + (A[2][2] * rb[2]);
262 +
263 + }
264 +
265 + double RigidBody::getZangle( ){
266 +    return zAngle;
267 + }
268 +
269 + void RigidBody::setZangle( double zAng ){
270 +    zAngle = zAng;
271 + }
272  
273 + void RigidBody::addZangle( double zAng ){
274 +    zAngle += zAng;
275   }
276  
277   void RigidBody::calcRefCoords( ) {
278  
279 <  int i,j,k;
279 >  int i,j,k, it, n_linear_coords;
280    double mtmp;
281    vec3 apos;
282    double refCOM[3];
283 +  vec3 ptmp;
284 +  double Itmp[3][3];
285 +  double evals[3];
286 +  double evects[3][3];
287 +  double r, r2, len;
288  
289 +  // First, find the center of mass:
290 +  
291    mass = 0.0;
292    for (j=0; j<3; j++)
293      refCOM[j] = 0.0;
# Line 293 | Line 306 | void RigidBody::calcRefCoords( ) {
306    for(j = 0; j < 3; j++)
307      refCOM[j] /= mass;
308  
309 + // Next, move the origin of the reference coordinate system to the COM:
310 +
311    for (i = 0; i < myAtoms.size(); i++) {
312      apos = refCoords[i];
313      for (j=0; j < 3; j++) {
314        apos[j] = apos[j] - refCOM[j];
315      }
316      refCoords[i] = apos;
317 +  }
318 +
319 + // Moment of Inertia calculation
320 +
321 +  for (i = 0; i < 3; i++)
322 +    for (j = 0; j < 3; j++)
323 +      Itmp[i][j] = 0.0;  
324 +  
325 +  for (it = 0; it < myAtoms.size(); it++) {
326 +
327 +    mtmp = myAtoms[it]->getMass();
328 +    ptmp = refCoords[it];
329 +    r= norm3(ptmp.vec);
330 +    r2 = r*r;
331 +    
332 +    for (i = 0; i < 3; i++) {
333 +      for (j = 0; j < 3; j++) {
334 +        
335 +        if (i==j) Itmp[i][j] += mtmp * r2;
336 +
337 +        Itmp[i][j] -= mtmp * ptmp.vec[i]*ptmp.vec[j];
338 +      }
339 +    }
340    }
341 +  
342 +  diagonalize3x3(Itmp, evals, sU);
343 +  
344 +  // zero out I and then fill the diagonals with the moments of inertia:
345  
346 +  n_linear_coords = 0;
347 +
348 +  for (i = 0; i < 3; i++) {
349 +    for (j = 0; j < 3; j++) {
350 +      I[i][j] = 0.0;  
351 +    }
352 +    I[i][i] = evals[i];
353 +
354 +    if (fabs(evals[i]) < momIntTol) {
355 +      is_linear = true;
356 +      n_linear_coords++;
357 +      linear_axis = i;
358 +    }
359 +  }
360 +
361 +  if (n_linear_coords > 1) {
362 +          sprintf( painCave.errMsg,
363 +               "RigidBody error.\n"
364 +               "\tOOPSE found more than one axis in this rigid body with a vanishing \n"
365 +               "\tmoment of inertia.  This can happen in one of three ways:\n"
366 +               "\t 1) Only one atom was specified, or \n"
367 +               "\t 2) All atoms were specified at the same location, or\n"
368 +               "\t 3) The programmers did something stupid.\n"
369 +               "\tIt is silly to use a rigid body to describe this situation.  Be smarter.\n"
370 +               );
371 +      painCave.isFatal = 1;
372 +      simError();
373 +  }
374 +  
375 +  // renormalize column vectors:
376 +  
377 +  for (i=0; i < 3; i++) {
378 +    len = 0.0;
379 +    for (j = 0; j < 3; j++) {
380 +      len += sU[i][j]*sU[i][j];
381 +    }
382 +    len = sqrt(len);
383 +    for (j = 0; j < 3; j++) {
384 +      sU[i][j] /= len;
385 +    }
386 +  }
387   }
388  
389   void RigidBody::doEulerToRotMat(vec3 &euler, mat3x3 &myA ){
# Line 366 | Line 449 | void RigidBody::calcForcesAndTorques() {
449    // (Actually, on second thought, don't.  Integrator does this now.)
450    // lab2Body(trq);
451  
369 –  forces_good = true;
370 –
452   }
453  
454   void RigidBody::updateAtoms() {
# Line 556 | Line 637 | void RigidBody::findCOM() {
637      vel[j] /= mass;
638    }
639  
559 –  com_good = true;
640   }
641 <  
642 < void RigidBody::findOrient() {
643 <  
644 <  size_t it;
565 <  int i, j;
566 <  double ptmp[3];
567 <  double Itmp[3][3];
568 <  double evals[3];
569 <  double evects[3][3];
570 <  double r2, mtmp, len;
641 >
642 > void RigidBody::accept(BaseVisitor* v){
643 >  vector<Atom*>::iterator atomIter;
644 >  v->visit(this);
645  
646 <  if (!com_good) findCOM();
646 >  //for(atomIter = myAtoms.begin(); atomIter != myAtoms.end(); ++atomIter)
647 >  //  (*atomIter)->accept(v);
648 > }
649  
650 <  // Calculate inertial tensor matrix elements:
651 <
576 <  for (i = 0; i < 3; i++)
577 <    for (j = 0; j < 3; j++)
578 <      Itmp[i][j] = 0.0;  
579 <  
580 <  for (it = 0; it < myAtoms.size(); it++) {
581 <    
582 <    mtmp = myAtoms[it]->getMass();    
583 <    myAtoms[it]->getPos(ptmp);
650 > void RigidBody::getAtomPos(double theP[3], int index){
651 >  vec3 ref;
652  
653 <    for (j = 0; j < 3; j++)
654 <      ptmp[j] = pos[j] - ptmp[j];
653 >  if (index >= myAtoms.size())
654 >    cerr << index << " is an invalid index, current rigid body contains " << myAtoms.size() << "atoms" << endl;
655  
656 <    r2 = norm3(ptmp);
657 <    
658 <    for (i = 0; i < 3; i++) {
659 <      for (j = 0; j < 3; j++) {
660 <        
661 <        if (i==j) Itmp[i][j] = mtmp * r2;
656 >  ref = refCoords[index];
657 >  body2Lab(ref.vec);
658 >  
659 >  theP[0] = pos[0] + ref[0];
660 >  theP[1] = pos[1] + ref[1];
661 >  theP[2] = pos[2] + ref[2];
662 > }
663  
595 –        Itmp[i][j] -= mtmp * ptmp[i]*ptmp[j];
596 –      }
597 –    }
598 –  }
599 –  
600 –  diagonalize3x3(Itmp, evals, sU);
601 –  
602 –  // zero out I and then fill the diagonals with the moments of inertia:
664  
665 <  for (i = 0; i < 3; i++) {
666 <    for (j = 0; j < 3; j++) {
606 <      I[i][j] = 0.0;  
607 <    }
608 <    I[i][i] = evals[i];
609 <  }
610 <  
611 <  // renormalize column vectors:
612 <  
613 <  for (i=0; i < 3; i++) {
614 <    len = 0.0;
615 <    for (j = 0; j < 3; j++) {
616 <      len += sU[i][j]*sU[i][j];
617 <    }
618 <    len = sqrt(len);
619 <    for (j = 0; j < 3; j++) {
620 <      sU[i][j] /= len;
621 <    }
622 <  }
623 <  
624 <  // sU now contains the coordinates of the 'special' frame;
625 <  
626 <  orient_good = true;
665 > void RigidBody::getAtomRefCoor(double pos[3], int index){
666 >  vec3 ref;
667  
668 +  ref = refCoords[index];
669 +  pos[0] = ref[0];
670 +  pos[1] = ref[1];
671 +  pos[2] = ref[2];
672 +  
673   }

Diff Legend

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