test/brains/DirectionalAtom.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */

#include "primitives/DirectionalAtom.hpp"

namespace oopse {

DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType) 
                         : Atom(dAtomType){
    objType_= otDAtom;
}

Mat3x3d DirectionalAtom::getI() {
    return static_cast<DirectionalAtomType*>(getAtomType())->getI();
}    

void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
    ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
    ((snapshotMan_->getPrevSnapshot())->*storage_).unitVector[localIndex_] = a.inverse() * sU_.getColum(2);
}

      
void DirectionalAtom::setA(const RotMat3x3d& a) {
    ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
    ((snapshotMan_->getCurrentSnapshot())->*storage_).unitVector[localIndex_] = a.inverse() * sU_.getColum(2);
}    
    
void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
    ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
    ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).unitVector[localIndex_] = a.inverse() * sU_.getColum(2);    
}    

void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
    setA(m *getA());
}

void  DirectionalAtom::setUnitFrameFromEuler(double phi, double theta, double psi) {
    sU_.setupRotMat(phi,theta,psi);
}

std::vector<double> DirectionalAtom::getGrad() {
    vector<double> grad(6, 0.0);
    Vector3d force;
    Vector3d torque;
    Vector3d myEuler;
    double phi, theta, psi;
    double cphi, sphi, ctheta, stheta;
    Vector3d ephi;
    Vector3d etheta;
    Vector3d epsi;

    force = getFrc();
    torque =getTrq();
    myEuler = getA().toEulerAngles();

    phi = myEuler[0];
    theta = myEuler[1];
    psi = myEuler[2];

    cphi = cos(phi);
    sphi = sin(phi);
    ctheta = cos(theta);
    stheta = sin(theta);

    // get unit vectors along the phi, theta and psi rotation axes

    ephi[0] = 0.0;
    ephi[1] = 0.0;
    ephi[2] = 1.0;

    etheta[0] = cphi;
    etheta[1] = sphi;
    etheta[2] = 0.0;

    epsi[0] = stheta * cphi;
    epsi[1] = stheta * sphi;
    epsi[2] = ctheta;

    //gradient is equal to -force
    for (int j = 0 ; j<3; j++)
        grad[j] = -force[j];

    for (int j = 0; j < 3; j++ ) {

        grad[3] += torque[j]*ephi[j];
        grad[4] += torque[j]*etheta[j];
        grad[5] += torque[j]*epsi[j];

    }
    
    return grad;
}    

void DirectionalAtom::accept(BaseVisitor* v) {
    v->visit(this);
}    

}

Revision:	197
Committed:	Mon Nov 1 19:57:07 2004 UTC (21 years ago) by tim
File size:	3573 byte(s)
Log Message:	these core classes get compiled
#	User	Rev	Content
1	tim	197	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.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
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
22			* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23			*
24			*/
25
26			#include "primitives/DirectionalAtom.hpp"
27
28			namespace oopse {
29
30			DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
31			: Atom(dAtomType){
32			objType_= otDAtom;
33			}
34
35			Mat3x3d DirectionalAtom::getI() {
36			return static_cast<DirectionalAtomType*>(getAtomType())->getI();
37			}
38
39			void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
40			((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
41			((snapshotMan_->getPrevSnapshot())->storage_).unitVector[localIndex_] = a.inverse() sU_.getColum(2);
42			}
43
44
45			void DirectionalAtom::setA(const RotMat3x3d& a) {
46			((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
47			((snapshotMan_->getCurrentSnapshot())->storage_).unitVector[localIndex_] = a.inverse() sU_.getColum(2);
48			}
49
50			void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
51			((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
52			((snapshotMan_->getSnapshot(snapshotNo))->storage_).unitVector[localIndex_] = a.inverse() sU_.getColum(2);
53			}
54
55			void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
56			setA(m *getA());
57			}
58
59			void DirectionalAtom::setUnitFrameFromEuler(double phi, double theta, double psi) {
60			sU_.setupRotMat(phi,theta,psi);
61			}
62
63			std::vector<double> DirectionalAtom::getGrad() {
64			vector<double> grad(6, 0.0);
65			Vector3d force;
66			Vector3d torque;
67			Vector3d myEuler;
68			double phi, theta, psi;
69			double cphi, sphi, ctheta, stheta;
70			Vector3d ephi;
71			Vector3d etheta;
72			Vector3d epsi;
73
74			force = getFrc();
75			torque =getTrq();
76			myEuler = getA().toEulerAngles();
77
78			phi = myEuler[0];
79			theta = myEuler[1];
80			psi = myEuler[2];
81
82			cphi = cos(phi);
83			sphi = sin(phi);
84			ctheta = cos(theta);
85			stheta = sin(theta);
86
87			// get unit vectors along the phi, theta and psi rotation axes
88
89			ephi[0] = 0.0;
90			ephi[1] = 0.0;
91			ephi[2] = 1.0;
92
93			etheta[0] = cphi;
94			etheta[1] = sphi;
95			etheta[2] = 0.0;
96
97			epsi[0] = stheta * cphi;
98			epsi[1] = stheta * sphi;
99			epsi[2] = ctheta;
100
101			//gradient is equal to -force
102			for (int j = 0 ; j<3; j++)
103			grad[j] = -force[j];
104
105			for (int j = 0; j < 3; j++ ) {
106
107			grad[3] += torque[j]*ephi[j];
108			grad[4] += torque[j]*etheta[j];
109			grad[5] += torque[j]*epsi[j];
110
111			}
112
113			return grad;
114			}
115
116			void DirectionalAtom::accept(BaseVisitor* v) {
117			v->visit(this);
118			}
119
120			}
121