ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/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.
branches/development/src/primitives/RigidBody.cpp (file contents), Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC

#	Line 1 | Line 1
1	<	/*
1	>	/*
2		* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3		*
4		* The University of Notre Dame grants you ("Licensee") a
#	Line 6 | Line 6
6		* redistribute this software in source and binary code form, provided
7		* that the following conditions are met:
8		*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13		*    notice, this list of conditions and the following disclaimer in the
14		*    documentation and/or other materials provided with the
15		*    distribution.
#	Line 37 | Line 28
28		* arising out of the use of or inability to use software, even if the
29		* University of Notre Dame has been advised of the possibility of
30		* such damages.
31	+	*
32	+	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	+	* research, please cite the appropriate papers when you publish your
34	+	* work.  Good starting points are:
35	+	*
36	+	* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	+	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	+	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	+	* [4]  Vardeman & Gezelter, in progress (2009).
40		*/
41		#include <algorithm>
42		#include <math.h>
43		#include "primitives/RigidBody.hpp"
44		#include "utils/simError.h"
45	<	namespace oopse {
46	<
47	<	RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){
48	<
49	<	}
50	<
51	<	void RigidBody::setPrevA(const RotMat3x3d& a) {
45	>	#include "utils/NumericConstant.hpp"
46	>	namespace OpenMD {
47	>
48	>	RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData),
49	>	inertiaTensor_(0.0){
50	>	}
51	>
52	>	void RigidBody::setPrevA(const RotMat3x3d& a) {
53		((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
54	<	//((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
54	<
54	>
55		for (int i =0 ; i < atoms_.size(); ++i){
56	<	if (atoms_[i]->isDirectional()) {
57	<	atoms_[i]->setPrevA(a * refOrients_[i]);
58	<	}
56	>	if (atoms_[i]->isDirectional()) {
57	>	atoms_[i]->setPrevA(refOrients_[i].transpose() * a);
58	>	}
59		}
60	<
61	<	}
62	<
63	<
64	<	void RigidBody::setA(const RotMat3x3d& a) {
60	>
61	>	}
62	>
63	>
64	>	void RigidBody::setA(const RotMat3x3d& a) {
65		((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
66	–	//((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
66
67		for (int i =0 ; i < atoms_.size(); ++i){
68	<	if (atoms_[i]->isDirectional()) {
69	<	atoms_[i]->setA(a * refOrients_[i]);
70	<	}
68	>	if (atoms_[i]->isDirectional()) {
69	>	atoms_[i]->setA(refOrients_[i].transpose() * a);
70	>	}
71		}
72	<	}
73	<
74	<	void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
72	>	}
73	>
74	>	void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
75		((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
76	+
77		//((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
78	<
78	>
79		for (int i =0 ; i < atoms_.size(); ++i){
80	<	if (atoms_[i]->isDirectional()) {
81	<	atoms_[i]->setA(a * refOrients_[i], snapshotNo);
82	<	}
80	>	if (atoms_[i]->isDirectional()) {
81	>	atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo);
82	>	}
83		}
84	<
85	<	}
86	<
87	<	Mat3x3d RigidBody::getI() {
84	>
85	>	}
86	>
87	>	Mat3x3d RigidBody::getI() {
88		return inertiaTensor_;
89	<	}
90	<
91	<	std::vector<double> RigidBody::getGrad() {
92	<	std::vector<double> grad(6, 0.0);
89	>	}
90	>
91	>	std::vector<RealType> RigidBody::getGrad() {
92	>	std::vector<RealType> grad(6, 0.0);
93		Vector3d force;
94		Vector3d torque;
95		Vector3d myEuler;
96	<	double phi, theta, psi;
97	<	double cphi, sphi, ctheta, stheta;
96	>	RealType phi, theta, psi;
97	>	RealType cphi, sphi, ctheta, stheta;
98		Vector3d ephi;
99		Vector3d etheta;
100		Vector3d epsi;
101	<
101	>
102		force = getFrc();
103		torque =getTrq();
104		myEuler = getA().toEulerAngles();
105	<
105	>
106		phi = myEuler[0];
107		theta = myEuler[1];
108		psi = myEuler[2];
109	<
109	>
110		cphi = cos(phi);
111		sphi = sin(phi);
112		ctheta = cos(theta);
113		stheta = sin(theta);
114	<
114	>
115		// get unit vectors along the phi, theta and psi rotation axes
116	<
116	>
117		ephi[0] = 0.0;
118		ephi[1] = 0.0;
119		ephi[2] = 1.0;
120	<
120	>
121	>	//etheta[0] = -sphi;
122	>	//etheta[1] =  cphi;
123	>	//etheta[2] =  0.0;
124	>
125		etheta[0] = cphi;
126		etheta[1] = sphi;
127	<	etheta[2] = 0.0;
128	<
127	>	etheta[2] =  0.0;
128	>
129		epsi[0] = stheta * cphi;
130		epsi[1] = stheta * sphi;
131		epsi[2] = ctheta;
132	<
132	>
133		//gradient is equal to -force
134		for (int j = 0 ; j<3; j++)
135	<	grad[j] = -force[j];
136	<
135	>	grad[j] = -force[j];
136	>
137		for (int j = 0; j < 3; j++ ) {
138	<
139	<	grad[3] += torque[j]*ephi[j];
140	<	grad[4] += torque[j]*etheta[j];
141	<	grad[5] += torque[j]*epsi[j];
142	<
138	>
139	>	grad[3] += torque[j]*ephi[j];
140	>	grad[4] += torque[j]*etheta[j];
141	>	grad[5] += torque[j]*epsi[j];
142	>
143		}
144
145		return grad;
146	<	}
147	<
148	<	void RigidBody::accept(BaseVisitor* v) {
146	>	}
147	>
148	>	void RigidBody::accept(BaseVisitor* v) {
149		v->visit(this);
150	<	}
150	>	}
151
152	<	/**@todo need modification */
153	<	void  RigidBody::calcRefCoords() {
154	<	double mtmp;
152	>	/**@todo need modification */
153	>	void  RigidBody::calcRefCoords() {
154	>	RealType mtmp;
155		Vector3d refCOM(0.0);
156		mass_ = 0.0;
157		for (std::size_t i = 0; i < atoms_.size(); ++i) {
158	<	mtmp = atoms_[i]->getMass();
159	<	mass_ += mtmp;
160	<	refCOM += refCoords_[i]*mtmp;
158	>	mtmp = atoms_[i]->getMass();
159	>	mass_ += mtmp;
160	>	refCOM += refCoords_[i]*mtmp;
161		}
162		refCOM /= mass_;
163	<
163	>
164		// Next, move the origin of the reference coordinate system to the COM:
165		for (std::size_t i = 0; i < atoms_.size(); ++i) {
166	<	refCoords_[i] -= refCOM;
166	>	refCoords_[i] -= refCOM;
167		}
168
169	<	// Moment of Inertia calculation
170	<	Mat3x3d Itmp(0.0);
167	<
169	>	// Moment of Inertia calculation
170	>	Mat3x3d Itmp(0.0);
171		for (std::size_t i = 0; i < atoms_.size(); i++) {
172	<	mtmp = atoms_[i]->getMass();
173	<	Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
174	<	double r2 = refCoords_[i].lengthSquare();
175	<	Itmp(0, 0) += mtmp * r2;
176	<	Itmp(1, 1) += mtmp * r2;
177	<	Itmp(2, 2) += mtmp * r2;
172	>	Mat3x3d IAtom(0.0);
173	>	mtmp = atoms_[i]->getMass();
174	>	IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
175	>	RealType r2 = refCoords_[i].lengthSquare();
176	>	IAtom(0, 0) += mtmp * r2;
177	>	IAtom(1, 1) += mtmp * r2;
178	>	IAtom(2, 2) += mtmp * r2;
179	>	Itmp += IAtom;
180	>
181	>	//project the inertial moment of directional atoms into this rigid body
182	>	if (atoms_[i]->isDirectional()) {
183	>	Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i];
184	>	}
185		}
186
187	+	//    std::cout << Itmp << std::endl;
188	+
189		//diagonalize
190		Vector3d evals;
191		Mat3x3d::diagonalize(Itmp, evals, sU_);
#	Line 185 | Line 197 | void  RigidBody::calcRefCoords() {
197
198		int nLinearAxis = 0;
199		for (int i = 0; i < 3; i++) {
200	<	if (fabs(evals[i]) < oopse::epsilon) {
201	<	linear_ = true;
202	<	linearAxis_ = i;
203	<	++ nLinearAxis;
204	<	}
200	>	if (fabs(evals[i]) < OpenMD::epsilon) {
201	>	linear_ = true;
202	>	linearAxis_ = i;
203	>	++ nLinearAxis;
204	>	}
205		}
206
207		if (nLinearAxis > 1) {
208	<	sprintf( painCave.errMsg,
209	<	"RigidBody error.\n"
210	<	"\tOOPSE found more than one axis in this rigid body with a vanishing \n"
211	<	"\tmoment of inertia.  This can happen in one of three ways:\n"
212	<	"\t 1) Only one atom was specified, or \n"
213	<	"\t 2) All atoms were specified at the same location, or\n"
214	<	"\t 3) The programmers did something stupid.\n"
215	<	"\tIt is silly to use a rigid body to describe this situation.  Be smarter.\n"
216	<	);
217	<	painCave.isFatal = 1;
218	<	simError();
208	>	sprintf( painCave.errMsg,
209	>	"RigidBody error.\n"
210	>	"\tOpenMD found more than one axis in this rigid body with a vanishing \n"
211	>	"\tmoment of inertia.  This can happen in one of three ways:\n"
212	>	"\t 1) Only one atom was specified, or \n"
213	>	"\t 2) All atoms were specified at the same location, or\n"
214	>	"\t 3) The programmers did something stupid.\n"
215	>	"\tIt is silly to use a rigid body to describe this situation.  Be smarter.\n"
216	>	);
217	>	painCave.isFatal = 1;
218	>	simError();
219		}
220
221	<	}
221	>	}
222
223	<	void  RigidBody::calcForcesAndTorques() {
223	>	void  RigidBody::calcForcesAndTorques() {
224		Vector3d afrc;
225		Vector3d atrq;
226		Vector3d apos;
227		Vector3d rpos;
228		Vector3d frc(0.0);
229	<	Vector3d trq(0.0);
229	>	Vector3d trq(0.0);
230		Vector3d pos = this->getPos();
231		for (int i = 0; i < atoms_.size(); i++) {
232
233	<	afrc = atoms_[i]->getFrc();
234	<	apos = atoms_[i]->getPos();
235	<	rpos = apos - pos;
233	>	afrc = atoms_[i]->getFrc();
234	>	apos = atoms_[i]->getPos();
235	>	rpos = apos - pos;
236
237	<	frc += afrc;
237	>	frc += afrc;
238
239	<	trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
240	<	trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
241	<	trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
239	>	trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
240	>	trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
241	>	trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
242
243	<	// If the atom has a torque associated with it, then we also need to
244	<	// migrate the torques onto the center of mass:
243	>	// If the atom has a torque associated with it, then we also need to
244	>	// migrate the torques onto the center of mass:
245
246	<	if (atoms_[i]->isDirectional()) {
247	<	atrq = atoms_[i]->getTrq();
248	<	trq += atrq;
249	<	}
246	>	if (atoms_[i]->isDirectional()) {
247	>	atrq = atoms_[i]->getTrq();
248	>	trq += atrq;
249	>	}
250	>	}
251	>	addFrc(frc);
252	>	addTrq(trq);
253	>	}
254	>
255	>	Mat3x3d RigidBody::calcForcesAndTorquesAndVirial() {
256	>	Vector3d afrc;
257	>	Vector3d atrq;
258	>	Vector3d apos;
259	>	Vector3d rpos;
260	>	Vector3d dfrc;
261	>	Vector3d frc(0.0);
262	>	Vector3d trq(0.0);
263	>	Vector3d pos = this->getPos();
264	>	Mat3x3d tau_(0.0);
265	>
266	>	for (int i = 0; i < atoms_.size(); i++) {
267	>
268	>	afrc = atoms_[i]->getFrc();
269	>	apos = atoms_[i]->getPos();
270	>	rpos = apos - pos;
271
272	+	frc += afrc;
273	+
274	+	trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
275	+	trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
276	+	trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
277	+
278	+	// If the atom has a torque associated with it, then we also need to
279	+	// migrate the torques onto the center of mass:
280	+
281	+	if (atoms_[i]->isDirectional()) {
282	+	atrq = atoms_[i]->getTrq();
283	+	trq += atrq;
284	+	}
285	+
286	+	tau_(0,0) -= rpos[0]*afrc[0];
287	+	tau_(0,1) -= rpos[0]*afrc[1];
288	+	tau_(0,2) -= rpos[0]*afrc[2];
289	+	tau_(1,0) -= rpos[1]*afrc[0];
290	+	tau_(1,1) -= rpos[1]*afrc[1];
291	+	tau_(1,2) -= rpos[1]*afrc[2];
292	+	tau_(2,0) -= rpos[2]*afrc[0];
293	+	tau_(2,1) -= rpos[2]*afrc[1];
294	+	tau_(2,2) -= rpos[2]*afrc[2];
295	+
296		}
297	<
298	<	setFrc(frc);
299	<	setTrq(trq);
300	<
244	<	}
297	>	addFrc(frc);
298	>	addTrq(trq);
299	>	return tau_;
300	>	}
301
302	<	void  RigidBody::updateAtoms() {
302	>	void  RigidBody::updateAtoms() {
303		unsigned int i;
304		Vector3d ref;
305		Vector3d apos;
#	Line 253 | Line 309 | void  RigidBody::updateAtoms() {
309
310		for (i = 0; i < atoms_.size(); i++) {
311
312	<	ref = body2Lab(refCoords_[i]);
312	>	ref = body2Lab(refCoords_[i]);
313
314	<	apos = pos + ref;
314	>	apos = pos + ref;
315
316	<	atoms_[i]->setPos(apos);
316	>	atoms_[i]->setPos(apos);
317
318	<	if (atoms_[i]->isDirectional()) {
318	>	if (atoms_[i]->isDirectional()) {
319
320	<	dAtom = (DirectionalAtom *) atoms_[i];
321	<	dAtom->setA(a * refOrients_[i]);
322	<	//dAtom->rotateBy( A );
267	<	}
320	>	dAtom = (DirectionalAtom *) atoms_[i];
321	>	dAtom->setA(refOrients_[i].transpose() * a);
322	>	}
323
324		}
325
326	<	}
326	>	}
327
328
329	<	bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
329	>	void  RigidBody::updateAtoms(int frame) {
330	>	unsigned int i;
331	>	Vector3d ref;
332	>	Vector3d apos;
333	>	DirectionalAtom* dAtom;
334	>	Vector3d pos = getPos(frame);
335	>	RotMat3x3d a = getA(frame);
336	>
337	>	for (i = 0; i < atoms_.size(); i++) {
338	>
339	>	ref = body2Lab(refCoords_[i], frame);
340	>
341	>	apos = pos + ref;
342	>
343	>	atoms_[i]->setPos(apos, frame);
344	>
345	>	if (atoms_[i]->isDirectional()) {
346	>
347	>	dAtom = (DirectionalAtom *) atoms_[i];
348	>	dAtom->setA(refOrients_[i].transpose() * a, frame);
349	>	}
350	>
351	>	}
352	>
353	>	}
354	>
355	>	void RigidBody::updateAtomVel() {
356	>	Mat3x3d skewMat;;
357	>
358	>	Vector3d ji = getJ();
359	>	Mat3x3d I =  getI();
360	>
361	>	skewMat(0, 0) =0;
362	>	skewMat(0, 1) = ji[2] /I(2, 2);
363	>	skewMat(0, 2) = -ji[1] /I(1, 1);
364	>
365	>	skewMat(1, 0) = -ji[2] /I(2, 2);
366	>	skewMat(1, 1) = 0;
367	>	skewMat(1, 2) = ji[0]/I(0, 0);
368	>
369	>	skewMat(2, 0) =ji[1] /I(1, 1);
370	>	skewMat(2, 1) = -ji[0]/I(0, 0);
371	>	skewMat(2, 2) = 0;
372	>
373	>	Mat3x3d mat = (getA() * skewMat).transpose();
374	>	Vector3d rbVel = getVel();
375	>
376	>
377	>	Vector3d velRot;
378	>	for (int i =0 ; i < refCoords_.size(); ++i) {
379	>	atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
380	>	}
381	>
382	>	}
383	>
384	>	void RigidBody::updateAtomVel(int frame) {
385	>	Mat3x3d skewMat;;
386	>
387	>	Vector3d ji = getJ(frame);
388	>	Mat3x3d I =  getI();
389	>
390	>	skewMat(0, 0) =0;
391	>	skewMat(0, 1) = ji[2] /I(2, 2);
392	>	skewMat(0, 2) = -ji[1] /I(1, 1);
393	>
394	>	skewMat(1, 0) = -ji[2] /I(2, 2);
395	>	skewMat(1, 1) = 0;
396	>	skewMat(1, 2) = ji[0]/I(0, 0);
397	>
398	>	skewMat(2, 0) =ji[1] /I(1, 1);
399	>	skewMat(2, 1) = -ji[0]/I(0, 0);
400	>	skewMat(2, 2) = 0;
401	>
402	>	Mat3x3d mat = (getA(frame) * skewMat).transpose();
403	>	Vector3d rbVel = getVel(frame);
404	>
405	>
406	>	Vector3d velRot;
407	>	for (int i =0 ; i < refCoords_.size(); ++i) {
408	>	atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
409	>	}
410	>
411	>	}
412	>
413	>
414	>
415	>	bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
416		if (index < atoms_.size()) {
417
418	<	Vector3d ref = body2Lab(refCoords_[index]);
419	<	pos = getPos() + ref;
420	<	return true;
418	>	Vector3d ref = body2Lab(refCoords_[index]);
419	>	pos = getPos() + ref;
420	>	return true;
421		} else {
422	<	std::cerr << index << " is an invalid index, current rigid body contains "
423	<	<< atoms_.size() << "atoms" << std::endl;
424	<	return false;
422	>	std::cerr << index << " is an invalid index, current rigid body contains "
423	>	<< atoms_.size() << "atoms" << std::endl;
424	>	return false;
425		}
426	<	}
426	>	}
427
428	<	bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) {
428	>	bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) {
429		std::vector<Atom*>::iterator i;
430		i = std::find(atoms_.begin(), atoms_.end(), atom);
431		if (i != atoms_.end()) {
432	<	//RigidBody class makes sure refCoords_ and atoms_ match each other
433	<	Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]);
434	<	pos = getPos() + ref;
435	<	return true;
432	>	//RigidBody class makes sure refCoords_ and atoms_ match each other
433	>	Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]);
434	>	pos = getPos() + ref;
435	>	return true;
436		} else {
437	<	std::cerr << "Atom " << atom->getGlobalIndex()
438	<	<<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
439	<	return false;
437	>	std::cerr << "Atom " << atom->getGlobalIndex()
438	>	<<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
439	>	return false;
440		}
441	<	}
442	<	bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) {
441	>	}
442	>	bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) {
443
444		//velRot = $(A\cdot skew(I^{-1}j))^{T}refCoor$
445
446		if (index < atoms_.size()) {
447
448	<	Vector3d velRot;
449	<	Mat3x3d skewMat;;
450	<	Vector3d ref = refCoords_[index];
451	<	Vector3d ji = getJ();
452	<	Mat3x3d I =  getI();
448	>	Vector3d velRot;
449	>	Mat3x3d skewMat;;
450	>	Vector3d ref = refCoords_[index];
451	>	Vector3d ji = getJ();
452	>	Mat3x3d I =  getI();
453
454	<	skewMat(0, 0) =0;
455	<	skewMat(0, 1) = ji[2] /I(2, 2);
456	<	skewMat(0, 2) = -ji[1] /I(1, 1);
454	>	skewMat(0, 0) =0;
455	>	skewMat(0, 1) = ji[2] /I(2, 2);
456	>	skewMat(0, 2) = -ji[1] /I(1, 1);
457
458	<	skewMat(1, 0) = -ji[2] /I(2, 2);
459	<	skewMat(1, 1) = 0;
460	<	skewMat(1, 2) = ji[0]/I(0, 0);
458	>	skewMat(1, 0) = -ji[2] /I(2, 2);
459	>	skewMat(1, 1) = 0;
460	>	skewMat(1, 2) = ji[0]/I(0, 0);
461
462	<	skewMat(2, 0) =ji[1] /I(1, 1);
463	<	skewMat(2, 1) = -ji[0]/I(0, 0);
464	<	skewMat(2, 2) = 0;
462	>	skewMat(2, 0) =ji[1] /I(1, 1);
463	>	skewMat(2, 1) = -ji[0]/I(0, 0);
464	>	skewMat(2, 2) = 0;
465
466	<	velRot = (getA() * skewMat).transpose() * ref;
466	>	velRot = (getA() * skewMat).transpose() * ref;
467
468	<	vel =getVel() + velRot;
469	<	return true;
468	>	vel =getVel() + velRot;
469	>	return true;
470
471		} else {
472	<	std::cerr << index << " is an invalid index, current rigid body contains "
473	<	<< atoms_.size() << "atoms" << std::endl;
474	<	return false;
472	>	std::cerr << index << " is an invalid index, current rigid body contains "
473	>	<< atoms_.size() << "atoms" << std::endl;
474	>	return false;
475		}
476	<	}
476	>	}
477
478	<	bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) {
478	>	bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) {
479
480		std::vector<Atom*>::iterator i;
481		i = std::find(atoms_.begin(), atoms_.end(), atom);
482		if (i != atoms_.end()) {
483	<	return getAtomVel(vel, i - atoms_.begin());
483	>	return getAtomVel(vel, i - atoms_.begin());
484		} else {
485	<	std::cerr << "Atom " << atom->getGlobalIndex()
486	<	<<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
487	<	return false;
485	>	std::cerr << "Atom " << atom->getGlobalIndex()
486	>	<<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
487	>	return false;
488		}
489	<	}
489	>	}
490
491	<	bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) {
491	>	bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) {
492		if (index < atoms_.size()) {
493
494	<	coor = refCoords_[index];
495	<	return true;
494	>	coor = refCoords_[index];
495	>	return true;
496		} else {
497	<	std::cerr << index << " is an invalid index, current rigid body contains "
498	<	<< atoms_.size() << "atoms" << std::endl;
499	<	return false;
497	>	std::cerr << index << " is an invalid index, current rigid body contains "
498	>	<< atoms_.size() << "atoms" << std::endl;
499	>	return false;
500		}
501
502	<	}
502	>	}
503
504	<	bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) {
504	>	bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) {
505		std::vector<Atom*>::iterator i;
506		i = std::find(atoms_.begin(), atoms_.end(), atom);
507		if (i != atoms_.end()) {
508	<	//RigidBody class makes sure refCoords_ and atoms_ match each other
509	<	coor = refCoords_[i - atoms_.begin()];
510	<	return true;
508	>	//RigidBody class makes sure refCoords_ and atoms_ match each other
509	>	coor = refCoords_[i - atoms_.begin()];
510	>	return true;
511		} else {
512	<	std::cerr << "Atom " << atom->getGlobalIndex()
513	<	<<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
514	<	return false;
512	>	std::cerr << "Atom " << atom->getGlobalIndex()
513	>	<<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
514	>	return false;
515		}
516
517	<	}
517	>	}
518
519
520	<	void RigidBody::addAtom(Atom* at, AtomStamp* ats) {
520	>	void RigidBody::addAtom(Atom* at, AtomStamp* ats) {
521
522	<	Vector3d coords;
523	<	Vector3d euler;
522	>	Vector3d coords;
523	>	Vector3d euler;
524
525
526	<	atoms_.push_back(at);
526	>	atoms_.push_back(at);
527
528	<	if( !ats->havePosition() ){
529	<	sprintf( painCave.errMsg,
530	<	"RigidBody error.\n"
531	<	"\tAtom %s does not have a position specified.\n"
532	<	"\tThis means RigidBody cannot set up reference coordinates.\n",
533	<	ats->getType() );
534	<	painCave.isFatal = 1;
535	<	simError();
536	<	}
528	>	if( !ats->havePosition() ){
529	>	sprintf( painCave.errMsg,
530	>	"RigidBody error.\n"
531	>	"\tAtom %s does not have a position specified.\n"
532	>	"\tThis means RigidBody cannot set up reference coordinates.\n",
533	>	ats->getType().c_str() );
534	>	painCave.isFatal = 1;
535	>	simError();
536	>	}
537
538	<	coords[0] = ats->getPosX();
539	<	coords[1] = ats->getPosY();
540	<	coords[2] = ats->getPosZ();
538	>	coords[0] = ats->getPosX();
539	>	coords[1] = ats->getPosY();
540	>	coords[2] = ats->getPosZ();
541
542	<	refCoords_.push_back(coords);
542	>	refCoords_.push_back(coords);
543
544	<	RotMat3x3d identMat = RotMat3x3d::identity();
544	>	RotMat3x3d identMat = RotMat3x3d::identity();
545
546	<	if (at->isDirectional()) {
546	>	if (at->isDirectional()) {
547
548	<	if( !ats->haveOrientation() ){
549	<	sprintf( painCave.errMsg,
550	<	"RigidBody error.\n"
551	<	"\tAtom %s does not have an orientation specified.\n"
552	<	"\tThis means RigidBody cannot set up reference orientations.\n",
553	<	ats->getType() );
554	<	painCave.isFatal = 1;
555	<	simError();
556	<	}
548	>	if( !ats->haveOrientation() ){
549	>	sprintf( painCave.errMsg,
550	>	"RigidBody error.\n"
551	>	"\tAtom %s does not have an orientation specified.\n"
552	>	"\tThis means RigidBody cannot set up reference orientations.\n",
553	>	ats->getType().c_str() );
554	>	painCave.isFatal = 1;
555	>	simError();
556	>	}
557
558	<	euler[0] = ats->getEulerPhi();
559	<	euler[1] = ats->getEulerTheta();
560	<	euler[2] = ats->getEulerPsi();
558	>	euler[0] = ats->getEulerPhi() * NumericConstant::PI /180.0;
559	>	euler[1] = ats->getEulerTheta() * NumericConstant::PI /180.0;
560	>	euler[2] = ats->getEulerPsi() * NumericConstant::PI /180.0;
561
562	<	RotMat3x3d Atmp(euler);
563	<	refOrients_.push_back(Atmp);
562	>	RotMat3x3d Atmp(euler);
563	>	refOrients_.push_back(Atmp);
564
565	<	}else {
566	<	refOrients_.push_back(identMat);
567	<	}
565	>	}else {
566	>	refOrients_.push_back(identMat);
567	>	}
568
569
570	<	}
570	>	}
571
572		}
573

Comparing:
trunk/src/primitives/RigidBody.cpp (property svn:keywords), Revision 253 by tim, Thu Jan 13 19:40:37 2005 UTC vs.
branches/development/src/primitives/RigidBody.cpp (property svn:keywords), Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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