50 |
|
#include "math/ParallelRandNumGen.hpp" |
51 |
|
#endif |
52 |
|
|
53 |
+ |
/* Remove me after testing*/ |
54 |
+ |
#include <cstdio> |
55 |
+ |
#include <iostream> |
56 |
+ |
/*End remove me*/ |
57 |
+ |
|
58 |
|
namespace oopse { |
59 |
|
|
60 |
|
Velocitizer::Velocitizer(SimInfo* info) : info_(info) { |
97 |
|
double av2; |
98 |
|
double kebar; |
99 |
|
|
100 |
+ |
Globals * simParams = info_->getSimParams(); |
101 |
+ |
|
102 |
|
SimInfo::MoleculeIterator i; |
103 |
|
Molecule::IntegrableObjectIterator j; |
104 |
|
Molecule * mol; |
156 |
|
|
157 |
|
|
158 |
|
removeComDrift(); |
159 |
+ |
// Remove angular drift if we are not using periodic boundary conditions. |
160 |
+ |
if(simParams->getPBC()) removeAngularDrift(); |
161 |
|
|
162 |
|
} |
163 |
|
|
185 |
|
|
186 |
|
} |
187 |
|
|
188 |
+ |
|
189 |
+ |
void Velocitizer::removeAngularDrift() { |
190 |
+ |
// Get the Center of Mass drift velocity. |
191 |
+ |
|
192 |
+ |
Vector3d vdrift; |
193 |
+ |
Vector3d com; |
194 |
+ |
|
195 |
+ |
info_->getComAll(com,vdrift); |
196 |
+ |
|
197 |
+ |
Mat3x3d inertiaTensor; |
198 |
+ |
Vector3d angularMomentum; |
199 |
+ |
Vector3d omega; |
200 |
+ |
|
201 |
+ |
|
202 |
+ |
|
203 |
+ |
info_->getInertiaTensor(inertiaTensor,angularMomentum); |
204 |
+ |
// We now need the inverse of the inertia tensor. |
205 |
+ |
|
206 |
+ |
std::cerr << "Angular Momentum before is " |
207 |
+ |
<< angularMomentum << std::endl; |
208 |
+ |
|
209 |
+ |
|
210 |
+ |
inertiaTensor.inverse(); |
211 |
+ |
|
212 |
+ |
|
213 |
+ |
omega = inertiaTensor*angularMomentum; |
214 |
+ |
|
215 |
+ |
SimInfo::MoleculeIterator i; |
216 |
+ |
Molecule::IntegrableObjectIterator j; |
217 |
+ |
Molecule * mol; |
218 |
+ |
StuntDouble * integrableObject; |
219 |
+ |
Vector3d tempComPos; |
220 |
+ |
|
221 |
+ |
// Corrects for the center of mass angular drift. |
222 |
+ |
// sums all the angular momentum and divides by total mass. |
223 |
+ |
for( mol = info_->beginMolecule(i); mol != NULL; |
224 |
+ |
mol = info_->nextMolecule(i) ) { |
225 |
+ |
for( integrableObject = mol->beginIntegrableObject(j); |
226 |
+ |
integrableObject != NULL; |
227 |
+ |
integrableObject = mol->nextIntegrableObject(j) ) { |
228 |
+ |
tempComPos = integrableObject->getPos()-com; |
229 |
+ |
integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos)); |
230 |
+ |
} |
231 |
+ |
} |
232 |
+ |
|
233 |
+ |
angularMomentum = info_->getAngularMomentum(); |
234 |
+ |
std::cerr << "Angular Momentum after is " |
235 |
+ |
<< angularMomentum << std::endl; |
236 |
+ |
|
237 |
+ |
|
238 |
+ |
} |
239 |
+ |
|
240 |
+ |
|
241 |
+ |
|
242 |
+ |
|
243 |
|
} |