# | 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 45 | Line 45 | namespace oopse { | |
45 | #include "utils/NumericConstant.hpp" | |
46 | namespace oopse { | |
47 | ||
48 | < | RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){ |
48 | > | RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){ |
49 | ||
50 | < | } |
50 | > | } |
51 | ||
52 | < | void RigidBody::setPrevA(const RotMat3x3d& a) { |
52 | > | void RigidBody::setPrevA(const RotMat3x3d& a) { |
53 | ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a; | |
54 | //((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_; | |
55 | ||
56 | for (int i =0 ; i < atoms_.size(); ++i){ | |
57 | < | if (atoms_[i]->isDirectional()) { |
58 | < | atoms_[i]->setPrevA(a * refOrients_[i]); |
59 | < | } |
57 | > | if (atoms_[i]->isDirectional()) { |
58 | > | atoms_[i]->setPrevA(a * refOrients_[i]); |
59 | > | } |
60 | } | |
61 | ||
62 | < | } |
62 | > | } |
63 | ||
64 | ||
65 | < | void RigidBody::setA(const RotMat3x3d& a) { |
65 | > | void RigidBody::setA(const RotMat3x3d& a) { |
66 | ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a; | |
67 | //((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_; | |
68 | ||
69 | for (int i =0 ; i < atoms_.size(); ++i){ | |
70 | < | if (atoms_[i]->isDirectional()) { |
71 | < | atoms_[i]->setA(a * refOrients_[i]); |
72 | < | } |
70 | > | if (atoms_[i]->isDirectional()) { |
71 | > | atoms_[i]->setA(a * refOrients_[i]); |
72 | > | } |
73 | } | |
74 | < | } |
74 | > | } |
75 | ||
76 | < | void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) { |
76 | > | void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) { |
77 | ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a; | |
78 | //((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_; | |
79 | ||
80 | for (int i =0 ; i < atoms_.size(); ++i){ | |
81 | < | if (atoms_[i]->isDirectional()) { |
82 | < | atoms_[i]->setA(a * refOrients_[i], snapshotNo); |
83 | < | } |
81 | > | if (atoms_[i]->isDirectional()) { |
82 | > | atoms_[i]->setA(a * refOrients_[i], snapshotNo); |
83 | > | } |
84 | } | |
85 | ||
86 | < | } |
86 | > | } |
87 | ||
88 | < | Mat3x3d RigidBody::getI() { |
88 | > | Mat3x3d RigidBody::getI() { |
89 | return inertiaTensor_; | |
90 | < | } |
90 | > | } |
91 | ||
92 | < | std::vector<double> RigidBody::getGrad() { |
93 | < | std::vector<double> grad(6, 0.0); |
92 | > | std::vector<double> RigidBody::getGrad() { |
93 | > | std::vector<double> grad(6, 0.0); |
94 | Vector3d force; | |
95 | Vector3d torque; | |
96 | Vector3d myEuler; | |
# | Line 129 | Line 129 | std::vector<double> RigidBody::getGrad() { | |
129 | ||
130 | //gradient is equal to -force | |
131 | for (int j = 0 ; j<3; j++) | |
132 | < | grad[j] = -force[j]; |
132 | > | grad[j] = -force[j]; |
133 | ||
134 | for (int j = 0; j < 3; j++ ) { | |
135 | ||
136 | < | grad[3] += torque[j]*ephi[j]; |
137 | < | grad[4] += torque[j]*etheta[j]; |
138 | < | grad[5] += torque[j]*epsi[j]; |
136 | > | grad[3] += torque[j]*ephi[j]; |
137 | > | grad[4] += torque[j]*etheta[j]; |
138 | > | grad[5] += torque[j]*epsi[j]; |
139 | ||
140 | } | |
141 | ||
142 | return grad; | |
143 | < | } |
143 | > | } |
144 | ||
145 | < | void RigidBody::accept(BaseVisitor* v) { |
145 | > | void RigidBody::accept(BaseVisitor* v) { |
146 | v->visit(this); | |
147 | < | } |
147 | > | } |
148 | ||
149 | < | /**@todo need modification */ |
150 | < | void RigidBody::calcRefCoords() { |
149 | > | /**@todo need modification */ |
150 | > | void RigidBody::calcRefCoords() { |
151 | double mtmp; | |
152 | Vector3d refCOM(0.0); | |
153 | mass_ = 0.0; | |
154 | for (std::size_t i = 0; i < atoms_.size(); ++i) { | |
155 | < | mtmp = atoms_[i]->getMass(); |
156 | < | mass_ += mtmp; |
157 | < | refCOM += refCoords_[i]*mtmp; |
155 | > | mtmp = atoms_[i]->getMass(); |
156 | > | mass_ += mtmp; |
157 | > | refCOM += refCoords_[i]*mtmp; |
158 | } | |
159 | refCOM /= mass_; | |
160 | ||
161 | // Next, move the origin of the reference coordinate system to the COM: | |
162 | for (std::size_t i = 0; i < atoms_.size(); ++i) { | |
163 | < | refCoords_[i] -= refCOM; |
163 | > | refCoords_[i] -= refCOM; |
164 | } | |
165 | ||
166 | < | // Moment of Inertia calculation |
167 | < | Mat3x3d Itmp(0.0); |
168 | < | |
166 | > | // Moment of Inertia calculation |
167 | > | Mat3x3d Itmp(0.0); |
168 | for (std::size_t i = 0; i < atoms_.size(); i++) { | |
169 | < | mtmp = atoms_[i]->getMass(); |
170 | < | Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp; |
171 | < | double r2 = refCoords_[i].lengthSquare(); |
172 | < | Itmp(0, 0) += mtmp * r2; |
173 | < | Itmp(1, 1) += mtmp * r2; |
174 | < | Itmp(2, 2) += mtmp * r2; |
175 | < | } |
169 | > | Mat3x3d IAtom(0.0); |
170 | > | mtmp = atoms_[i]->getMass(); |
171 | > | IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp; |
172 | > | double r2 = refCoords_[i].lengthSquare(); |
173 | > | IAtom(0, 0) += mtmp * r2; |
174 | > | IAtom(1, 1) += mtmp * r2; |
175 | > | IAtom(2, 2) += mtmp * r2; |
176 | ||
177 | < | //project the inertial moment of directional atoms into this rigid body |
178 | < | for (std::size_t i = 0; i < atoms_.size(); i++) { |
179 | < | if (atoms_[i]->isDirectional()) { |
180 | < | RectMatrix<double, 3, 3> Iproject = refOrients_[i].transpose() * atoms_[i]->getI(); |
181 | < | Itmp(0, 0) += Iproject(0, 0); |
182 | < | Itmp(1, 1) += Iproject(1, 1); |
183 | < | Itmp(2, 2) += Iproject(2, 2); |
185 | < | } |
177 | > | //project the inertial moment of directional atoms into this rigid body |
178 | > | if (atoms_[i]->isDirectional()) { |
179 | > | IAtom += atoms_[i]->getI(); |
180 | > | Itmp += refOrients_[i].transpose() * IAtom * refOrients_[i]; |
181 | > | } else { |
182 | > | Itmp += IAtom; |
183 | > | } |
184 | } | |
185 | ||
186 | //diagonalize | |
# | Line 196 | Line 194 | void RigidBody::calcRefCoords() { | |
194 | ||
195 | int nLinearAxis = 0; | |
196 | for (int i = 0; i < 3; i++) { | |
197 | < | if (fabs(evals[i]) < oopse::epsilon) { |
198 | < | linear_ = true; |
199 | < | linearAxis_ = i; |
200 | < | ++ nLinearAxis; |
201 | < | } |
197 | > | if (fabs(evals[i]) < oopse::epsilon) { |
198 | > | linear_ = true; |
199 | > | linearAxis_ = i; |
200 | > | ++ nLinearAxis; |
201 | > | } |
202 | } | |
203 | ||
204 | if (nLinearAxis > 1) { | |
205 | < | sprintf( painCave.errMsg, |
206 | < | "RigidBody error.\n" |
207 | < | "\tOOPSE found more than one axis in this rigid body with a vanishing \n" |
208 | < | "\tmoment of inertia. This can happen in one of three ways:\n" |
209 | < | "\t 1) Only one atom was specified, or \n" |
210 | < | "\t 2) All atoms were specified at the same location, or\n" |
211 | < | "\t 3) The programmers did something stupid.\n" |
212 | < | "\tIt is silly to use a rigid body to describe this situation. Be smarter.\n" |
213 | < | ); |
214 | < | painCave.isFatal = 1; |
215 | < | simError(); |
205 | > | sprintf( painCave.errMsg, |
206 | > | "RigidBody error.\n" |
207 | > | "\tOOPSE found more than one axis in this rigid body with a vanishing \n" |
208 | > | "\tmoment of inertia. This can happen in one of three ways:\n" |
209 | > | "\t 1) Only one atom was specified, or \n" |
210 | > | "\t 2) All atoms were specified at the same location, or\n" |
211 | > | "\t 3) The programmers did something stupid.\n" |
212 | > | "\tIt is silly to use a rigid body to describe this situation. Be smarter.\n" |
213 | > | ); |
214 | > | painCave.isFatal = 1; |
215 | > | simError(); |
216 | } | |
217 | ||
218 | < | } |
218 | > | } |
219 | ||
220 | < | void RigidBody::calcForcesAndTorques() { |
220 | > | void RigidBody::calcForcesAndTorques() { |
221 | Vector3d afrc; | |
222 | Vector3d atrq; | |
223 | Vector3d apos; | |
# | Line 229 | Line 227 | void RigidBody::calcForcesAndTorques() { | |
227 | Vector3d pos = this->getPos(); | |
228 | for (int i = 0; i < atoms_.size(); i++) { | |
229 | ||
230 | < | afrc = atoms_[i]->getFrc(); |
231 | < | apos = atoms_[i]->getPos(); |
232 | < | rpos = apos - pos; |
230 | > | afrc = atoms_[i]->getFrc(); |
231 | > | apos = atoms_[i]->getPos(); |
232 | > | rpos = apos - pos; |
233 | ||
234 | < | frc += afrc; |
234 | > | frc += afrc; |
235 | ||
236 | < | trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1]; |
237 | < | trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2]; |
238 | < | trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0]; |
236 | > | trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1]; |
237 | > | trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2]; |
238 | > | trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0]; |
239 | ||
240 | < | // If the atom has a torque associated with it, then we also need to |
241 | < | // migrate the torques onto the center of mass: |
240 | > | // If the atom has a torque associated with it, then we also need to |
241 | > | // migrate the torques onto the center of mass: |
242 | ||
243 | < | if (atoms_[i]->isDirectional()) { |
244 | < | atrq = atoms_[i]->getTrq(); |
245 | < | trq += atrq; |
246 | < | } |
243 | > | if (atoms_[i]->isDirectional()) { |
244 | > | atrq = atoms_[i]->getTrq(); |
245 | > | trq += atrq; |
246 | > | } |
247 | ||
248 | } | |
249 | ||
250 | setFrc(frc); | |
251 | setTrq(trq); | |
252 | ||
253 | < | } |
253 | > | } |
254 | ||
255 | < | void RigidBody::updateAtoms() { |
255 | > | void RigidBody::updateAtoms() { |
256 | unsigned int i; | |
257 | Vector3d ref; | |
258 | Vector3d apos; | |
# | Line 264 | Line 262 | void RigidBody::updateAtoms() { | |
262 | ||
263 | for (i = 0; i < atoms_.size(); i++) { | |
264 | ||
265 | < | ref = body2Lab(refCoords_[i]); |
265 | > | ref = body2Lab(refCoords_[i]); |
266 | ||
267 | < | apos = pos + ref; |
267 | > | apos = pos + ref; |
268 | ||
269 | < | atoms_[i]->setPos(apos); |
269 | > | atoms_[i]->setPos(apos); |
270 | ||
271 | < | if (atoms_[i]->isDirectional()) { |
271 | > | if (atoms_[i]->isDirectional()) { |
272 | ||
273 | < | dAtom = (DirectionalAtom *) atoms_[i]; |
274 | < | dAtom->setA(a * refOrients_[i]); |
275 | < | //dAtom->rotateBy( A ); |
278 | < | } |
273 | > | dAtom = (DirectionalAtom *) atoms_[i]; |
274 | > | dAtom->setA(refOrients_[i] * a); |
275 | > | } |
276 | ||
277 | } | |
278 | ||
279 | < | } |
279 | > | } |
280 | ||
281 | ||
282 | < | void RigidBody::updateAtoms(int frame) { |
282 | > | void RigidBody::updateAtoms(int frame) { |
283 | unsigned int i; | |
284 | Vector3d ref; | |
285 | Vector3d apos; | |
# | Line 292 | Line 289 | void RigidBody::updateAtoms(int frame) { | |
289 | ||
290 | for (i = 0; i < atoms_.size(); i++) { | |
291 | ||
292 | < | ref = body2Lab(refCoords_[i], frame); |
292 | > | ref = body2Lab(refCoords_[i], frame); |
293 | ||
294 | < | apos = pos + ref; |
294 | > | apos = pos + ref; |
295 | ||
296 | < | atoms_[i]->setPos(apos, frame); |
296 | > | atoms_[i]->setPos(apos, frame); |
297 | ||
298 | < | if (atoms_[i]->isDirectional()) { |
298 | > | if (atoms_[i]->isDirectional()) { |
299 | ||
300 | < | dAtom = (DirectionalAtom *) atoms_[i]; |
301 | < | dAtom->setA(a * refOrients_[i], frame); |
302 | < | } |
300 | > | dAtom = (DirectionalAtom *) atoms_[i]; |
301 | > | dAtom->setA(refOrients_[i] * a, frame); |
302 | > | } |
303 | ||
304 | } | |
305 | ||
306 | < | } |
306 | > | } |
307 | ||
308 | < | void RigidBody::updateAtomVel() { |
308 | > | void RigidBody::updateAtomVel() { |
309 | Mat3x3d skewMat;; | |
310 | ||
311 | Vector3d ji = getJ(); | |
# | Line 332 | Line 329 | void RigidBody::updateAtomVel() { | |
329 | ||
330 | Vector3d velRot; | |
331 | for (int i =0 ; i < refCoords_.size(); ++i) { | |
332 | < | atoms_[i]->setVel(rbVel + mat * refCoords_[i]); |
332 | > | atoms_[i]->setVel(rbVel + mat * refCoords_[i]); |
333 | } | |
334 | ||
335 | < | } |
335 | > | } |
336 | ||
337 | < | void RigidBody::updateAtomVel(int frame) { |
337 | > | void RigidBody::updateAtomVel(int frame) { |
338 | Mat3x3d skewMat;; | |
339 | ||
340 | Vector3d ji = getJ(frame); | |
# | Line 361 | Line 358 | void RigidBody::updateAtomVel(int frame) { | |
358 | ||
359 | Vector3d velRot; | |
360 | for (int i =0 ; i < refCoords_.size(); ++i) { | |
361 | < | atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame); |
361 | > | atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame); |
362 | } | |
363 | ||
364 | < | } |
364 | > | } |
365 | ||
366 | ||
367 | ||
368 | < | bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) { |
368 | > | bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) { |
369 | if (index < atoms_.size()) { | |
370 | ||
371 | < | Vector3d ref = body2Lab(refCoords_[index]); |
372 | < | pos = getPos() + ref; |
373 | < | return true; |
371 | > | Vector3d ref = body2Lab(refCoords_[index]); |
372 | > | pos = getPos() + ref; |
373 | > | return true; |
374 | } else { | |
375 | < | std::cerr << index << " is an invalid index, current rigid body contains " |
376 | < | << atoms_.size() << "atoms" << std::endl; |
377 | < | return false; |
375 | > | std::cerr << index << " is an invalid index, current rigid body contains " |
376 | > | << atoms_.size() << "atoms" << std::endl; |
377 | > | return false; |
378 | } | |
379 | < | } |
379 | > | } |
380 | ||
381 | < | bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) { |
381 | > | bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) { |
382 | std::vector<Atom*>::iterator i; | |
383 | i = std::find(atoms_.begin(), atoms_.end(), atom); | |
384 | if (i != atoms_.end()) { | |
385 | < | //RigidBody class makes sure refCoords_ and atoms_ match each other |
386 | < | Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]); |
387 | < | pos = getPos() + ref; |
388 | < | return true; |
385 | > | //RigidBody class makes sure refCoords_ and atoms_ match each other |
386 | > | Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]); |
387 | > | pos = getPos() + ref; |
388 | > | return true; |
389 | } else { | |
390 | < | std::cerr << "Atom " << atom->getGlobalIndex() |
391 | < | <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl; |
392 | < | return false; |
390 | > | std::cerr << "Atom " << atom->getGlobalIndex() |
391 | > | <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl; |
392 | > | return false; |
393 | } | |
394 | < | } |
395 | < | bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) { |
394 | > | } |
395 | > | bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) { |
396 | ||
397 | //velRot = $(A\cdot skew(I^{-1}j))^{T}refCoor$ | |
398 | ||
399 | if (index < atoms_.size()) { | |
400 | ||
401 | < | Vector3d velRot; |
402 | < | Mat3x3d skewMat;; |
403 | < | Vector3d ref = refCoords_[index]; |
404 | < | Vector3d ji = getJ(); |
405 | < | Mat3x3d I = getI(); |
401 | > | Vector3d velRot; |
402 | > | Mat3x3d skewMat;; |
403 | > | Vector3d ref = refCoords_[index]; |
404 | > | Vector3d ji = getJ(); |
405 | > | Mat3x3d I = getI(); |
406 | ||
407 | < | skewMat(0, 0) =0; |
408 | < | skewMat(0, 1) = ji[2] /I(2, 2); |
409 | < | skewMat(0, 2) = -ji[1] /I(1, 1); |
407 | > | skewMat(0, 0) =0; |
408 | > | skewMat(0, 1) = ji[2] /I(2, 2); |
409 | > | skewMat(0, 2) = -ji[1] /I(1, 1); |
410 | ||
411 | < | skewMat(1, 0) = -ji[2] /I(2, 2); |
412 | < | skewMat(1, 1) = 0; |
413 | < | skewMat(1, 2) = ji[0]/I(0, 0); |
411 | > | skewMat(1, 0) = -ji[2] /I(2, 2); |
412 | > | skewMat(1, 1) = 0; |
413 | > | skewMat(1, 2) = ji[0]/I(0, 0); |
414 | ||
415 | < | skewMat(2, 0) =ji[1] /I(1, 1); |
416 | < | skewMat(2, 1) = -ji[0]/I(0, 0); |
417 | < | skewMat(2, 2) = 0; |
415 | > | skewMat(2, 0) =ji[1] /I(1, 1); |
416 | > | skewMat(2, 1) = -ji[0]/I(0, 0); |
417 | > | skewMat(2, 2) = 0; |
418 | ||
419 | < | velRot = (getA() * skewMat).transpose() * ref; |
419 | > | velRot = (getA() * skewMat).transpose() * ref; |
420 | ||
421 | < | vel =getVel() + velRot; |
422 | < | return true; |
421 | > | vel =getVel() + velRot; |
422 | > | return true; |
423 | ||
424 | } else { | |
425 | < | std::cerr << index << " is an invalid index, current rigid body contains " |
426 | < | << atoms_.size() << "atoms" << std::endl; |
427 | < | return false; |
425 | > | std::cerr << index << " is an invalid index, current rigid body contains " |
426 | > | << atoms_.size() << "atoms" << std::endl; |
427 | > | return false; |
428 | } | |
429 | < | } |
429 | > | } |
430 | ||
431 | < | bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) { |
431 | > | bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) { |
432 | ||
433 | std::vector<Atom*>::iterator i; | |
434 | i = std::find(atoms_.begin(), atoms_.end(), atom); | |
435 | if (i != atoms_.end()) { | |
436 | < | return getAtomVel(vel, i - atoms_.begin()); |
436 | > | return getAtomVel(vel, i - atoms_.begin()); |
437 | } else { | |
438 | < | std::cerr << "Atom " << atom->getGlobalIndex() |
439 | < | <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl; |
440 | < | return false; |
438 | > | std::cerr << "Atom " << atom->getGlobalIndex() |
439 | > | <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl; |
440 | > | return false; |
441 | } | |
442 | < | } |
442 | > | } |
443 | ||
444 | < | bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) { |
444 | > | bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) { |
445 | if (index < atoms_.size()) { | |
446 | ||
447 | < | coor = refCoords_[index]; |
448 | < | return true; |
447 | > | coor = refCoords_[index]; |
448 | > | return true; |
449 | } else { | |
450 | < | std::cerr << index << " is an invalid index, current rigid body contains " |
451 | < | << atoms_.size() << "atoms" << std::endl; |
452 | < | return false; |
450 | > | std::cerr << index << " is an invalid index, current rigid body contains " |
451 | > | << atoms_.size() << "atoms" << std::endl; |
452 | > | return false; |
453 | } | |
454 | ||
455 | < | } |
455 | > | } |
456 | ||
457 | < | bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) { |
457 | > | bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) { |
458 | std::vector<Atom*>::iterator i; | |
459 | i = std::find(atoms_.begin(), atoms_.end(), atom); | |
460 | if (i != atoms_.end()) { | |
461 | < | //RigidBody class makes sure refCoords_ and atoms_ match each other |
462 | < | coor = refCoords_[i - atoms_.begin()]; |
463 | < | return true; |
461 | > | //RigidBody class makes sure refCoords_ and atoms_ match each other |
462 | > | coor = refCoords_[i - atoms_.begin()]; |
463 | > | return true; |
464 | } else { | |
465 | < | std::cerr << "Atom " << atom->getGlobalIndex() |
466 | < | <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl; |
467 | < | return false; |
465 | > | std::cerr << "Atom " << atom->getGlobalIndex() |
466 | > | <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl; |
467 | > | return false; |
468 | } | |
469 | ||
470 | < | } |
470 | > | } |
471 | ||
472 | ||
473 | < | void RigidBody::addAtom(Atom* at, AtomStamp* ats) { |
473 | > | void RigidBody::addAtom(Atom* at, AtomStamp* ats) { |
474 | ||
475 | < | Vector3d coords; |
476 | < | Vector3d euler; |
475 | > | Vector3d coords; |
476 | > | Vector3d euler; |
477 | ||
478 | ||
479 | < | atoms_.push_back(at); |
479 | > | atoms_.push_back(at); |
480 | ||
481 | < | if( !ats->havePosition() ){ |
482 | < | sprintf( painCave.errMsg, |
483 | < | "RigidBody error.\n" |
484 | < | "\tAtom %s does not have a position specified.\n" |
485 | < | "\tThis means RigidBody cannot set up reference coordinates.\n", |
486 | < | ats->getType() ); |
487 | < | painCave.isFatal = 1; |
488 | < | simError(); |
489 | < | } |
481 | > | if( !ats->havePosition() ){ |
482 | > | sprintf( painCave.errMsg, |
483 | > | "RigidBody error.\n" |
484 | > | "\tAtom %s does not have a position specified.\n" |
485 | > | "\tThis means RigidBody cannot set up reference coordinates.\n", |
486 | > | ats->getType() ); |
487 | > | painCave.isFatal = 1; |
488 | > | simError(); |
489 | > | } |
490 | ||
491 | < | coords[0] = ats->getPosX(); |
492 | < | coords[1] = ats->getPosY(); |
493 | < | coords[2] = ats->getPosZ(); |
491 | > | coords[0] = ats->getPosX(); |
492 | > | coords[1] = ats->getPosY(); |
493 | > | coords[2] = ats->getPosZ(); |
494 | ||
495 | < | refCoords_.push_back(coords); |
495 | > | refCoords_.push_back(coords); |
496 | ||
497 | < | RotMat3x3d identMat = RotMat3x3d::identity(); |
497 | > | RotMat3x3d identMat = RotMat3x3d::identity(); |
498 | ||
499 | < | if (at->isDirectional()) { |
499 | > | if (at->isDirectional()) { |
500 | ||
501 | < | if( !ats->haveOrientation() ){ |
502 | < | sprintf( painCave.errMsg, |
503 | < | "RigidBody error.\n" |
504 | < | "\tAtom %s does not have an orientation specified.\n" |
505 | < | "\tThis means RigidBody cannot set up reference orientations.\n", |
506 | < | ats->getType() ); |
507 | < | painCave.isFatal = 1; |
508 | < | simError(); |
509 | < | } |
501 | > | if( !ats->haveOrientation() ){ |
502 | > | sprintf( painCave.errMsg, |
503 | > | "RigidBody error.\n" |
504 | > | "\tAtom %s does not have an orientation specified.\n" |
505 | > | "\tThis means RigidBody cannot set up reference orientations.\n", |
506 | > | ats->getType() ); |
507 | > | painCave.isFatal = 1; |
508 | > | simError(); |
509 | > | } |
510 | ||
511 | < | euler[0] = ats->getEulerPhi() * NumericConstant::PI /180.0; |
512 | < | euler[1] = ats->getEulerTheta() * NumericConstant::PI /180.0; |
513 | < | euler[2] = ats->getEulerPsi() * NumericConstant::PI /180.0; |
511 | > | euler[0] = ats->getEulerPhi() * NumericConstant::PI /180.0; |
512 | > | euler[1] = ats->getEulerTheta() * NumericConstant::PI /180.0; |
513 | > | euler[2] = ats->getEulerPsi() * NumericConstant::PI /180.0; |
514 | ||
515 | < | RotMat3x3d Atmp(euler); |
516 | < | refOrients_.push_back(Atmp); |
515 | > | RotMat3x3d Atmp(euler); |
516 | > | refOrients_.push_back(Atmp); |
517 | ||
518 | < | }else { |
519 | < | refOrients_.push_back(identMat); |
520 | < | } |
518 | > | }else { |
519 | > | refOrients_.push_back(identMat); |
520 | > | } |
521 | ||
522 | ||
523 | < | } |
523 | > | } |
524 | ||
525 | } | |
526 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |