| 1 |
|
/* |
| 2 |
< |
* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project |
| 2 |
> |
* Copyright (C) 2000-2009 The Open Molecular Dynamics Engine (OpenMD) project |
| 3 |
|
* |
| 4 |
< |
* Contact: oopse@oopse.org |
| 4 |
> |
* Contact: gezelter@openscience.org |
| 5 |
|
* |
| 6 |
|
* This program is free software; you can redistribute it and/or |
| 7 |
|
* modify it under the terms of the GNU Lesser General Public License |
| 25 |
|
|
| 26 |
|
#include "primitives/RigidBody.hpp" |
| 27 |
|
|
| 28 |
< |
namespace oopse { |
| 28 |
> |
namespace OpenMD { |
| 29 |
|
|
| 30 |
< |
RigidBody::RigidBody() : objType_(otRigidBody), storage_(&Snapshot::rigidbodyData){ |
| 30 |
> |
RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData){ |
| 31 |
|
|
| 32 |
|
} |
| 33 |
|
|
| 34 |
|
void RigidBody::setPrevA(const RotMat3x3d& a) { |
| 35 |
< |
(snapshotMan_->getPrevSnapshot())->storage_->aMat[localIndex_] = a; |
| 36 |
< |
(snapshotMan_->getPrevSnapshot())->storage_->unitVector[localIndex_] = a.inverse() * sU_.getColum(2); |
| 35 |
> |
((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a; |
| 36 |
> |
((snapshotMan_->getPrevSnapshot())->*storage_).unitVector[localIndex_] = a.inverse() * sU_.getColumn(2); |
| 37 |
|
|
| 38 |
|
std::vector<Atom*>::iterator i; |
| 39 |
|
for (i = atoms_.begin(); i != atoms_.end(); ++i) { |
| 46 |
|
|
| 47 |
|
|
| 48 |
|
void RigidBody::setA(const RotMat3x3d& a) { |
| 49 |
< |
(snapshotMan_->getCurrentSnapshot())->storage_->aMat[localIndex_] = a; |
| 50 |
< |
(snapshotMan_->getCurrentSnapshot())->storage_->unitVector[localIndex_] = a.inverse() * sU_.getColum(2); |
| 49 |
> |
((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a; |
| 50 |
> |
((snapshotMan_->getCurrentSnapshot())->*storage_).unitVector[localIndex_] = a.inverse() * sU_.getColumn(2); |
| 51 |
|
|
| 52 |
|
std::vector<Atom*>::iterator i; |
| 53 |
|
for (i = atoms_.begin(); i != atoms_.end(); ++i) { |
| 58 |
|
} |
| 59 |
|
|
| 60 |
|
void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) { |
| 61 |
< |
(snapshotMan_->getSnapshot(snapshotNo))->storage_->aMat[localIndex_] = a; |
| 62 |
< |
(snapshotMan_->getSnapshot(snapshotNo))->storage_->unitVector[localIndex_] = a.inverse() * sU_.getColum(2); |
| 61 |
> |
((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a; |
| 62 |
> |
((snapshotMan_->getSnapshot(snapshotNo))->*storage_).unitVector[localIndex_] = a.inverse() * sU_.getColumn(2); |
| 63 |
|
|
| 64 |
|
std::vector<Atom*>::iterator i; |
| 65 |
|
for (i = atoms_.begin(); i != atoms_.end(); ++i) { |
| 78 |
|
return inertiaTensor_; |
| 79 |
|
} |
| 80 |
|
|
| 81 |
– |
void RigidBody::setI(Mat3x3d& I) { |
| 82 |
– |
inertiaTensor_ = I; |
| 83 |
– |
} |
| 84 |
– |
|
| 81 |
|
std::vector<double> RigidBody::getGrad() { |
| 82 |
|
vector<double> grad(6, 0.0); |
| 83 |
|
Vector3d force; |
| 221 |
|
if (n_linear_coords > 1) { |
| 222 |
|
sprintf( painCave.errMsg, |
| 223 |
|
"RigidBody error.\n" |
| 224 |
< |
"\tOOPSE found more than one axis in this rigid body with a vanishing \n" |
| 224 |
> |
"\tOpenMD found more than one axis in this rigid body with a vanishing \n" |
| 225 |
|
"\tmoment of inertia. This can happen in one of three ways:\n" |
| 226 |
|
"\t 1) Only one atom was specified, or \n" |
| 227 |
|
"\t 2) All atoms were specified at the same location, or\n" |
| 322 |
|
|
| 323 |
|
|
| 324 |
|
bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) { |
| 325 |
< |
if (index < atoms_.size() { |
| 325 |
> |
if (index < atoms_.size()) { |
| 326 |
|
|
| 327 |
|
Vector3d ref = body2Lab(refCoords_[index]); |
| 328 |
|
pos = getPos() + ref; |
| 329 |
< |
return true |
| 329 |
> |
return true; |
| 330 |
|
} else { |
| 331 |
|
std::cerr << index << " is an invalid index, current rigid body contains " |
| 332 |
|
<< atoms_.size() << "atoms" << std::endl; |
| 345 |
|
} else { |
| 346 |
|
std::cerr << "Atom " << atom->getGlobalIndex() |
| 347 |
|
<<" does not belong to Rigid body "<< getGlobalIndex() << std::endl; |
| 348 |
+ |
return false; |
| 349 |
|
} |
| 350 |
|
} |
| 351 |
|
bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) { |
| 352 |
|
|
| 353 |
|
//velRot = $(A\cdot skew(I^{-1}j))^{T}refCoor$ |
| 354 |
|
|
| 355 |
< |
if (index < atoms_.size() { |
| 355 |
> |
if (index < atoms_.size()) { |
| 356 |
|
|
| 360 |
– |
Vector3d ref; |
| 357 |
|
Vector3d velRot; |
| 358 |
|
Mat3x3d skewMat;; |
| 359 |
|
Vector3d ref = refCoords_[index]; |
| 375 |
|
velRot = (getA() * skewMat).transpose() * ref; |
| 376 |
|
|
| 377 |
|
vel =getVel() + velRot; |
| 378 |
+ |
return true; |
| 379 |
|
|
| 380 |
|
} else { |
| 381 |
< |
std::cerr << "Atom " << atom->getGlobalIndex() |
| 382 |
< |
<<" does not belong to Rigid body "<< getGlobalIndex() << std::endl; |
| 381 |
> |
std::cerr << index << " is an invalid index, current rigid body contains " |
| 382 |
> |
<< atoms_.size() << "atoms" << std::endl; |
| 383 |
|
return false; |
| 384 |
|
} |
| 385 |
|
} |
| 398 |
|
} |
| 399 |
|
|
| 400 |
|
bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) { |
| 401 |
< |
if (index < atoms_.size() { |
| 401 |
> |
if (index < atoms_.size()) { |
| 402 |
|
|
| 403 |
|
coor = refCoords_[index]; |
| 404 |
< |
return true |
| 404 |
> |
return true; |
| 405 |
|
} else { |
| 406 |
|
std::cerr << index << " is an invalid index, current rigid body contains " |
| 407 |
|
<< atoms_.size() << "atoms" << std::endl; |
