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

Comparing trunk/OOPSE/libmdtools/Thermo.cpp (file contents):
Revision 1253 by gezelter, Tue Jun 8 16:49:46 2004 UTC vs.
Revision 1452 by tim, Mon Aug 23 15:11:36 2004 UTC

# Line 11 | Line 11 | using namespace std;
11   #include "Integrator.hpp"
12   #include "simError.h"
13   #include "MatVec3.h"
14 + #include "ConstraintManager.hpp"
15 + #include "Mat3x3d.hpp"
16  
17   #ifdef IS_MPI
18   #define __C
# Line 26 | Line 28 | Thermo::Thermo( SimInfo* the_info ) {
28    int baseSeed = the_info->getSeed();
29    
30    gaussStream = new gaussianSPRNG( baseSeed );
31 +
32 +  cpIter = info->consMan->createPairIterator();
33   }
34  
35   Thermo::~Thermo(){
36    delete gaussStream;
37 +  delete cpIter;
38   }
39  
40   double Thermo::getKinetic(){
# Line 438 | Line 443 | void Thermo::removeCOMdrift() {
443      info->integrableObjects[vd]->setVel( aVel );
444    }
445   }
446 +
447 + void Thermo::removeAngularMomentum(){
448 +  Vector3d vcom;
449 +  Vector3d qcom;
450 +  Vector3d pos;
451 +  Vector3d vel;
452 +  double mass;  
453 +  double xx;
454 +  double yy;
455 +  double zz;
456 +  double xy;
457 +  double xz;
458 +  double yz;
459 +  Vector3d localAngMom;
460 +  Vector3d angMom;
461 +  Vector3d omega;
462 +  vector<StuntDouble *> integrableObjects;
463 +  double localInertiaVec[9];
464 +  double inertiaVec[9];
465 +  vector<Vector3d> qMinusQCom;
466 +  vector<Vector3d> vMinusVCom;
467 +  Mat3x3d inertiaMat;
468 +  Mat3x3d inverseInertiaMat;
469 +  
470 +  integrableObjects = info->integrableObjects;
471 +  qMinusQCom.resize(integrableObjects.size());
472 +  vMinusVCom.resize(integrableObjects.size());
473 +  
474 +  getCOM(qcom.vec);
475 +  getCOMVel(vcom.vec);
476 +        
477 +  //initialize components for inertia tensor
478 +  xx = 0.0;
479 +  yy = 0.0;
480 +  zz = 0.0;
481 +  xy = 0.0;
482 +  xz = 0.0;
483 +  yz = 0.0;
484 +  
485 +   //build components of Inertia tensor
486 +  //
487 +  //       [  Ixx -Ixy  -Ixz ]
488 +  //   J = | -Iyx  Iyy  -Iyz |
489 +  //       [ -Izx -Iyz   Izz ]
490 +  //See Fowles and Cassidy Chapter 9 or Goldstein Chapter 5
491 +  for(size_t i = 0; i < integrableObjects.size(); i++){
492 +    integrableObjects[i]->getPos(pos.vec);
493 +    integrableObjects[i]->getVel(vel.vec);
494 +    mass = integrableObjects[i]->getMass();
495 +    
496 +    qMinusQCom[i] = pos - qcom;
497 +    info->wrapVector(qMinusQCom[i].vec);
498 +    
499 +    vMinusVCom[i] = vel - vcom;
500 +
501 +    //compute moment of inertia coefficents
502 +    xx += qMinusQCom[i].x * qMinusQCom[i].x * mass;
503 +    yy += qMinusQCom[i].y * qMinusQCom[i].y * mass;
504 +    zz += qMinusQCom[i].z * qMinusQCom[i].z * mass;
505 +
506 +    // compute products of inertia
507 +    xy += qMinusQCom[i].x * qMinusQCom[i].y * mass;
508 +    xz += qMinusQCom[i].x * qMinusQCom[i].z * mass;
509 +    yz += qMinusQCom[i].y * qMinusQCom[i].z * mass;
510 +
511 +    localAngMom += crossProduct(qMinusQCom[i] , vMinusVCom[i] ) * mass;
512 +    
513 +  }
514 +
515 +  localInertiaVec[0] =yy+zz;
516 +  localInertiaVec[1] = -xy;
517 +  localInertiaVec[2] = -xz;
518 +  localInertiaVec[3] = -xy;
519 +  localInertiaVec[4] = xx+zz;
520 +  localInertiaVec[5] = -yz;
521 +  localInertiaVec[6] = -xz;
522 +  localInertiaVec[7] = -yz;
523 +  localInertiaVec[8] = xx+yy;
524 +
525 +  //Sum and distribute inertia and angmom arrays
526 + #ifdef MPI
527 +
528 +  MPI_Allreduce(localInertiaVec, inertiaVec, 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
529 +
530 +  MPI_Allreduce(localAngMom.vec, angMom.vec, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
531 +
532 +  inertiaMat.element[0][0] = inertiaVec[0];
533 +  inertiaMat.element[0][1] = inertiaVec[1];
534 +  inertiaMat.element[0][2] = inertiaVec[2];
535 +
536 +  inertiaMat.element[1][0] = inertiaVec[3];
537 +  inertiaMat.element[1][1] = inertiaVec[4];
538 +  inertiaMat.element[1][2] = inertiaVec[5];
539 +
540 +  inertiaMat.element[2][0] = inertiaVec[6];
541 +  inertiaMat.element[2][1] = inertiaVec[7];
542 +  inertiaMat.element[2][2] = inertiaVec[8];
543 +
544 + #else
545 +
546 +    inertiaMat.element[0][0] = localInertiaVec[0];
547 +    inertiaMat.element[0][1] = localInertiaVec[1];
548 +    inertiaMat.element[0][2] = localInertiaVec[2];
549 +
550 +    inertiaMat.element[1][0] = localInertiaVec[3];
551 +    inertiaMat.element[1][1] = localInertiaVec[4];
552 +    inertiaMat.element[1][2] = localInertiaVec[5];
553 +
554 +    inertiaMat.element[2][0] = localInertiaVec[6];
555 +    inertiaMat.element[2][1] = localInertiaVec[7];
556 +    inertiaMat.element[2][2] = localInertiaVec[8];
557 +  
558 +    angMom     = localAngMom;
559 + #endif
560 +
561 +    //invert the moment of inertia tensor by LU-decomposition / backsolving:
562 +
563 +    inverseInertiaMat = inertiaMat.inverse();
564 +
565 +    //calculate the angular velocities: omega = I^-1 . L
566 +
567 +    omega = inverseInertiaMat * angMom;
568 +
569 +    //subtract out center of mass velocity and angular momentum from
570 +    //particle velocities
571 +
572 +    for(size_t i = 0; i < integrableObjects.size(); i++){
573 +      vel = vMinusVCom[i] - crossProduct(omega, qMinusQCom[i]);
574 +      integrableObjects[i]->setVel(vel.vec);      
575 +    }
576 + }
577 +
578 + double Thermo::getConsEnergy(){
579 +  ConstraintPair* consPair;
580 +  double totConsEnergy;
581 +  double bondLen2;
582 +  double dist;
583 +  double lamda;
584 +  
585 +  totConsEnergy = 0;
586 +  
587 +  for(cpIter->first(); !cpIter->isEnd(); cpIter->next()){
588 +    consPair =  cpIter->currentItem();
589 +    bondLen2 = consPair->getBondLength2();
590 +    lamda = consPair->getLamda();
591 +    //dist = consPair->getDistance();
592 +
593 +    //totConsEnergy += lamda * (dist*dist - bondLen2);
594 +  }
595 +
596 +  return totConsEnergy;
597 + }
598 +
599 +

Diff Legend

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