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;
    if (dAtomType->isMultipole()) {
        electroBodyFrame_ = dAtomType->getElectroBodyFrame();
    }
}

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

void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
    ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
    if (atomType_->isMultipole()) {
        ((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
    }
}

      
void DirectionalAtom::setA(const RotMat3x3d& a) {
    ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;

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

    if (atomType_->isMultipole()) {
        ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;    
    }
}    

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

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:	1815
Committed:	Wed Dec 1 18:42:45 2004 UTC (19 years, 7 months ago) by tim
File size:	3711 byte(s)
Log Message:	except integrator and constraint, other directories get built
#	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	tim	1813	if (dAtomType->isMultipole()) {
34	tim	1815	electroBodyFrame_ = dAtomType->getElectroBodyFrame();
35	tim	1813	}
36	gezelter	1490	}
37
38	tim	1692	Mat3x3d DirectionalAtom::getI() {
39			return static_cast<DirectionalAtomType*>(getAtomType())->getI();
40			}
41	gezelter	1490
42	tim	1692	void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
43			((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
44	tim	1813	if (atomType_->isMultipole()) {
45			((snapshotMan_->getPrevSnapshot())->storage_).electroFrame[localIndex_] = a.transpose() electroBodyFrame_;
46			}
47	gezelter	1490	}
48
49	tim	1692
50			void DirectionalAtom::setA(const RotMat3x3d& a) {
51			((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
52	tim	1813
53			if (atomType_->isMultipole()) {
54			((snapshotMan_->getCurrentSnapshot())->storage_).electroFrame[localIndex_] = a.transpose() electroBodyFrame_;
55			}
56	tim	1692	}
57	gezelter	1490
58	tim	1692	void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
59			((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
60	tim	1813
61			if (atomType_->isMultipole()) {
62			((snapshotMan_->getSnapshot(snapshotNo))->storage_).electroFrame[localIndex_] = a.transpose() electroBodyFrame_;
63			}
64	tim	1692	}
65	gezelter	1490
66	tim	1692	void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
67			setA(m *getA());
68	gezelter	1490	}
69
70	tim	1692	std::vector<double> DirectionalAtom::getGrad() {
71			vector<double> grad(6, 0.0);
72			Vector3d force;
73			Vector3d torque;
74			Vector3d myEuler;
75			double phi, theta, psi;
76			double cphi, sphi, ctheta, stheta;
77			Vector3d ephi;
78			Vector3d etheta;
79			Vector3d epsi;
80	gezelter	1490
81	tim	1692	force = getFrc();
82			torque =getTrq();
83			myEuler = getA().toEulerAngles();
84	gezelter	1490
85	tim	1692	phi = myEuler[0];
86			theta = myEuler[1];
87			psi = myEuler[2];
88	gezelter	1490
89	tim	1692	cphi = cos(phi);
90			sphi = sin(phi);
91			ctheta = cos(theta);
92			stheta = sin(theta);
93	gezelter	1490
94	tim	1692	// get unit vectors along the phi, theta and psi rotation axes
95	gezelter	1490
96	tim	1692	ephi[0] = 0.0;
97			ephi[1] = 0.0;
98			ephi[2] = 1.0;
99	gezelter	1490
100	tim	1692	etheta[0] = cphi;
101			etheta[1] = sphi;
102			etheta[2] = 0.0;
103	gezelter	1490
104	tim	1692	epsi[0] = stheta * cphi;
105			epsi[1] = stheta * sphi;
106			epsi[2] = ctheta;
107	gezelter	1490
108	tim	1692	//gradient is equal to -force
109			for (int j = 0 ; j<3; j++)
110			grad[j] = -force[j];
111	gezelter	1490
112	tim	1692	for (int j = 0; j < 3; j++ ) {
113	gezelter	1490
114	tim	1692	grad[3] += torque[j]*ephi[j];
115			grad[4] += torque[j]*etheta[j];
116			grad[5] += torque[j]*epsi[j];
117	gezelter	1490
118			}
119
120	tim	1692	return grad;
121			}
122	gezelter	1490
123	tim	1692	void DirectionalAtom::accept(BaseVisitor* v) {
124			v->visit(this);
125			}
126	gezelter	1490
127			}
128