[ViewVC] Diff of: OpenMD/branches/development/src/primitives/RigidBody.cpp

Comparing trunk/src/primitives/RigidBody.cpp (file contents):
Revision 253 by tim, Thu Jan 13 19:40:37 2005 UTC vs.
Revision 374 by tim, Tue Feb 22 18:56:25 2005 UTC

+#include <math.h>
+#include "primitives/RigidBody.hpp"
+#include "utils/simError.h"
-+
+#include "utils/NumericConstant.hpp"
+namespace oopse {
+RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){
+        Itmp(0, 0) += mtmp * r2;
+        Itmp(1, 1) += mtmp * r2;
+        Itmp(2, 2) += mtmp * r2;
-+
+    }
-+
-+
+    //project the inertial moment of directional atoms into this rigid body
-+
+    for (std::size_t i = 0; i < atoms_.size(); i++) {
-+
+        if (atoms_[i]->isDirectional()) {
-+
+            RectMatrix<double, 3, 3> 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);
-+
+        }
+    //diagonalize
-+
+void  RigidBody::updateAtoms(int frame) {
-+
+    unsigned int i;
-+
+    Vector3d ref;
-+
+    Vector3d apos;
-+
+    DirectionalAtom* dAtom;
-+
+    Vector3d pos = getPos(frame);
-+
+    RotMat3x3d a = getA(frame);
-+
-+
+    for (i = 0; i < atoms_.size(); i++) {
-+
-+
+        ref = body2Lab(refCoords_[i], frame);
-+
-+
+        apos = pos + ref;
-+
-+
+        atoms_[i]->setPos(apos, frame);
-+
-+
+        if (atoms_[i]->isDirectional()) {
-+
-+
+          dAtom = (DirectionalAtom *) atoms_[i];
-+
+          dAtom->setA(a * refOrients_[i], frame);
-+
+        }
-+
-+
+    }
-+
-+
+}
-+
-+
+void RigidBody::updateAtomVel() {
-+
+    Mat3x3d skewMat;;
-+
-+
+    Vector3d ji = getJ();
-+
+    Mat3x3d I =  getI();
-+
-+
+    skewMat(0, 0) =0;
-+
+    skewMat(0, 1) = ji[2] /I(2, 2);
-+
+    skewMat(0, 2) = -ji[1] /I(1, 1);
-+
-+
+    skewMat(1, 0) = -ji[2] /I(2, 2);
-+
+    skewMat(1, 1) = 0;
-+
+    skewMat(1, 2) = ji[0]/I(0, 0);
-+
-+
+    skewMat(2, 0) =ji[1] /I(1, 1);
-+
+    skewMat(2, 1) = -ji[0]/I(0, 0);
-+
+    skewMat(2, 2) = 0;
-+
-+
+    Mat3x3d mat = (getA() * skewMat).transpose();
-+
+    Vector3d rbVel = getVel();
-+
-+
-+
+    Vector3d velRot;
-+
+    for (int i =0 ; i < refCoords_.size(); ++i) {
-+
+        atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
-+
+    }
-+
-+
+}
-+
-+
+void RigidBody::updateAtomVel(int frame) {
-+
+    Mat3x3d skewMat;;
-+
-+
+    Vector3d ji = getJ(frame);
-+
+    Mat3x3d I =  getI();
-+
-+
+    skewMat(0, 0) =0;
-+
+    skewMat(0, 1) = ji[2] /I(2, 2);
-+
+    skewMat(0, 2) = -ji[1] /I(1, 1);
-+
-+
+    skewMat(1, 0) = -ji[2] /I(2, 2);
-+
+    skewMat(1, 1) = 0;
-+
+    skewMat(1, 2) = ji[0]/I(0, 0);
-+
-+
+    skewMat(2, 0) =ji[1] /I(1, 1);
-+
+    skewMat(2, 1) = -ji[0]/I(0, 0);
-+
+    skewMat(2, 2) = 0;
-+
-+
+    Mat3x3d mat = (getA(frame) * skewMat).transpose();
-+
+    Vector3d rbVel = getVel(frame);
-+
-+
-+
+    Vector3d velRot;
-+
+    for (int i =0 ; i < refCoords_.size(); ++i) {
-+
+        atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
-+
+    }
-+
-+
+}
-+
-+
-+
+bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
+    if (index < atoms_.size()) {
+      simError();
-<
+    euler[0] = ats->getEulerPhi();
-<
+    euler[1] = ats->getEulerTheta();
-<
+    euler[2] = ats->getEulerPsi();
->
+    euler[0] = ats->getEulerPhi() * NumericConstant::PI /180.0;
->
+    euler[1] = ats->getEulerTheta() * NumericConstant::PI /180.0;
->
+    euler[2] = ats->getEulerPsi() * NumericConstant::PI /180.0;
+    RotMat3x3d Atmp(euler);
+    refOrients_.push_back(Atmp);

Diff Legend

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

Comparing trunk/src/primitives/RigidBody.cpp (file contents): Revision 253 by tim, Thu Jan 13 19:40:37 2005 UTC vs. Revision 374 by tim, Tue Feb 22 18:56:25 2005 UTC

Diff Legend

Comparing trunk/src/primitives/RigidBody.cpp (file contents):
Revision 253 by tim, Thu Jan 13 19:40:37 2005 UTC vs.
Revision 374 by tim, Tue Feb 22 18:56:25 2005 UTC