src/primitives/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_).unitFrame[localIndex_] = a.transpose() * sU_;
}

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

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:	1775
Committed:	Wed Nov 24 17:42:37 2004 UTC (19 years, 7 months ago) by tim
File size:	3540 byte(s)
Log Message:	add SetForceManager method to integrator
#	User	Rev	Content
1	tim	1692	/*
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	gezelter	1490
26	tim	1492	#include "primitives/DirectionalAtom.hpp"
27	gezelter	1490
28	tim	1692	namespace oopse {
29	gezelter	1490
30	tim	1692	DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
31			: Atom(dAtomType){
32			objType_= otDAtom;
33	gezelter	1490	}
34
35	tim	1692	Mat3x3d DirectionalAtom::getI() {
36			return static_cast<DirectionalAtomType*>(getAtomType())->getI();
37			}
38	gezelter	1490
39	tim	1692	void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
40			((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
41	tim	1775	((snapshotMan_->getPrevSnapshot())->storage_).unitFrame[localIndex_] = a.transpose() sU_;
42	gezelter	1490	}
43
44	tim	1692
45			void DirectionalAtom::setA(const RotMat3x3d& a) {
46			((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
47	tim	1775	((snapshotMan_->getCurrentSnapshot())->storage_).unitFrame[localIndex_] = a.transpose() sU_;
48	tim	1692	}
49	gezelter	1490
50	tim	1692	void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
51			((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
52	tim	1775	((snapshotMan_->getSnapshot(snapshotNo))->storage_).unitFrame[localIndex_] = a.transpose() sU_;
53	tim	1692	}
54	gezelter	1490
55	tim	1692	void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
56			setA(m *getA());
57	gezelter	1490	}
58
59	tim	1692	void DirectionalAtom::setUnitFrameFromEuler(double phi, double theta, double psi) {
60			sU_.setupRotMat(phi,theta,psi);
61	gezelter	1490	}
62
63	tim	1692	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	gezelter	1490
74	tim	1692	force = getFrc();
75			torque =getTrq();
76			myEuler = getA().toEulerAngles();
77	gezelter	1490
78	tim	1692	phi = myEuler[0];
79			theta = myEuler[1];
80			psi = myEuler[2];
81	gezelter	1490
82	tim	1692	cphi = cos(phi);
83			sphi = sin(phi);
84			ctheta = cos(theta);
85			stheta = sin(theta);
86	gezelter	1490
87	tim	1692	// get unit vectors along the phi, theta and psi rotation axes
88	gezelter	1490
89	tim	1692	ephi[0] = 0.0;
90			ephi[1] = 0.0;
91			ephi[2] = 1.0;
92	gezelter	1490
93	tim	1692	etheta[0] = cphi;
94			etheta[1] = sphi;
95			etheta[2] = 0.0;
96	gezelter	1490
97	tim	1692	epsi[0] = stheta * cphi;
98			epsi[1] = stheta * sphi;
99			epsi[2] = ctheta;
100	gezelter	1490
101	tim	1692	//gradient is equal to -force
102			for (int j = 0 ; j<3; j++)
103			grad[j] = -force[j];
104	gezelter	1490
105	tim	1692	for (int j = 0; j < 3; j++ ) {
106	gezelter	1490
107	tim	1692	grad[3] += torque[j]*ephi[j];
108			grad[4] += torque[j]*etheta[j];
109			grad[5] += torque[j]*epsi[j];
110	gezelter	1490
111			}
112
113	tim	1692	return grad;
114			}
115	gezelter	1490
116	tim	1692	void DirectionalAtom::accept(BaseVisitor* v) {
117			v->visit(this);
118			}
119	gezelter	1490
120			}
121