ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/integrators/Velocitizer.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/integrators/Velocitizer.cpp (file contents):
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2364 by tim, Thu Oct 13 22:26:47 2005 UTC

# Line 50 | Line 50 | namespace oopse {
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) {
# Line 92 | Line 97 | namespace oopse {
97      double av2;
98      double kebar;
99      
100 +    Globals * simParams = info_->getSimParams();
101 +    
102      SimInfo::MoleculeIterator i;
103      Molecule::IntegrableObjectIterator j;
104      Molecule * mol;
# Line 149 | Line 156 | namespace oopse {
156      
157      
158      removeComDrift();
159 +    // Remove angular drift if we are not using periodic boundary conditions.
160 +    if(!simParams->getUsePeriodicBoundaryConditions()) removeAngularDrift();
161      
162    }
163    
# Line 176 | Line 185 | namespace oopse {
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 +      std::cerr << "Inertia Tensor before is "
209 +                << inertiaTensor <<  std::endl;
210 +      */
211 +      
212 +      inertiaTensor =inertiaTensor.inverse();
213 +      /*
214 +       std::cerr << "Inertia Tensor after inverse is "
215 +                << inertiaTensor <<  std::endl;
216 +      */
217 +      omega = inertiaTensor*angularMomentum;
218 +      
219 +      SimInfo::MoleculeIterator i;
220 +      Molecule::IntegrableObjectIterator j;
221 +      Molecule * mol;
222 +      StuntDouble * integrableObject;
223 +      Vector3d tempComPos;
224 +      
225 +      //  Corrects for the center of mass angular drift.
226 +      // sums all the angular momentum and divides by total mass.
227 +      for( mol = info_->beginMolecule(i); mol != NULL;
228 +           mol = info_->nextMolecule(i) ) {
229 +         for( integrableObject = mol->beginIntegrableObject(j);
230 +              integrableObject != NULL;
231 +              integrableObject = mol->nextIntegrableObject(j) ) {
232 +            tempComPos = integrableObject->getPos()-com;
233 +            integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos));
234 +         }
235 +      }
236 +      
237 +      angularMomentum = info_->getAngularMomentum();
238 +      /*
239 +      std::cerr << "Angular Momentum after is "
240 +         << angularMomentum <<  std::endl;
241 +      */
242 +      
243 +   }
244 +  
245 +  
246 +  
247 +  
248   }

Diff Legend

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