--- trunk/src/primitives/RigidBody.cpp 2005/04/15 22:04:00 507 +++ trunk/src/primitives/RigidBody.cpp 2006/05/17 21:51:42 963 @@ -51,11 +51,10 @@ namespace oopse { void RigidBody::setPrevA(const RotMat3x3d& a) { ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a; - //((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_; for (int i =0 ; i < atoms_.size(); ++i){ if (atoms_[i]->isDirectional()) { - atoms_[i]->setPrevA(a * refOrients_[i]); + atoms_[i]->setPrevA(refOrients_[i].transpose() * a); } } @@ -64,11 +63,10 @@ namespace oopse { void RigidBody::setA(const RotMat3x3d& a) { ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a; - //((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_; for (int i =0 ; i < atoms_.size(); ++i){ if (atoms_[i]->isDirectional()) { - atoms_[i]->setA(a * refOrients_[i]); + atoms_[i]->setA(refOrients_[i].transpose() * a); } } } @@ -79,7 +77,7 @@ namespace oopse { for (int i =0 ; i < atoms_.size(); ++i){ if (atoms_[i]->isDirectional()) { - atoms_[i]->setA(a * refOrients_[i], snapshotNo); + atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo); } } @@ -89,13 +87,13 @@ namespace oopse { return inertiaTensor_; } - std::vector RigidBody::getGrad() { - std::vector grad(6, 0.0); + std::vector RigidBody::getGrad() { + std::vector grad(6, 0.0); Vector3d force; Vector3d torque; Vector3d myEuler; - double phi, theta, psi; - double cphi, sphi, ctheta, stheta; + RealType phi, theta, psi; + RealType cphi, sphi, ctheta, stheta; Vector3d ephi; Vector3d etheta; Vector3d epsi; @@ -148,7 +146,7 @@ namespace oopse { /**@todo need modification */ void RigidBody::calcRefCoords() { - double mtmp; + RealType mtmp; Vector3d refCOM(0.0); mass_ = 0.0; for (std::size_t i = 0; i < atoms_.size(); ++i) { @@ -164,26 +162,24 @@ namespace oopse { } // Moment of Inertia calculation - Mat3x3d Itmp(0.0); - + Mat3x3d Itmp(0.0); for (std::size_t i = 0; i < atoms_.size(); i++) { + Mat3x3d IAtom(0.0); mtmp = atoms_[i]->getMass(); - Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp; - double r2 = refCoords_[i].lengthSquare(); - Itmp(0, 0) += mtmp * r2; - Itmp(1, 1) += mtmp * r2; - Itmp(2, 2) += mtmp * r2; - } + IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp; + RealType r2 = refCoords_[i].lengthSquare(); + IAtom(0, 0) += mtmp * r2; + IAtom(1, 1) += mtmp * r2; + IAtom(2, 2) += mtmp * r2; + Itmp += IAtom; - //project the inertial moment of directional atoms into this rigid body - for (std::size_t i = 0; i < atoms_.size(); i++) { + //project the inertial moment of directional atoms into this rigid body if (atoms_[i]->isDirectional()) { - RectMatrix Iproject = refOrients_[i].transpose() * atoms_[i]->getI(); - Itmp(0, 0) += Iproject(0, 0); - Itmp(1, 1) += Iproject(1, 1); - Itmp(2, 2) += Iproject(2, 2); - } + Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i]; + } } + + // std::cout << Itmp << std::endl; //diagonalize Vector3d evals; @@ -249,8 +245,8 @@ namespace oopse { } - setFrc(frc); - setTrq(trq); + addFrc(frc); + addTrq(trq); } @@ -273,8 +269,7 @@ namespace oopse { if (atoms_[i]->isDirectional()) { dAtom = (DirectionalAtom *) atoms_[i]; - dAtom->setA(a * refOrients_[i]); - //dAtom->rotateBy( A ); + dAtom->setA(refOrients_[i].transpose() * a); } } @@ -301,7 +296,7 @@ namespace oopse { if (atoms_[i]->isDirectional()) { dAtom = (DirectionalAtom *) atoms_[i]; - dAtom->setA(a * refOrients_[i], frame); + dAtom->setA(refOrients_[i].transpose() * a, frame); } } @@ -486,7 +481,7 @@ namespace oopse { "RigidBody error.\n" "\tAtom %s does not have a position specified.\n" "\tThis means RigidBody cannot set up reference coordinates.\n", - ats->getType() ); + ats->getType().c_str() ); painCave.isFatal = 1; simError(); } @@ -506,7 +501,7 @@ namespace oopse { "RigidBody error.\n" "\tAtom %s does not have an orientation specified.\n" "\tThis means RigidBody cannot set up reference orientations.\n", - ats->getType() ); + ats->getType().c_str() ); painCave.isFatal = 1; simError(); }