ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/integrators/LDForceManager.cpp
(Generate patch)

Comparing branches/development/src/integrators/LDForceManager.cpp (file contents):
Revision 1710 by gezelter, Fri May 18 21:44:02 2012 UTC vs.
Revision 1769 by gezelter, Mon Jul 9 14:15:52 2012 UTC

# Line 95 | Line 95 | namespace OpenMD {
95      std::map<std::string, HydroProp*> hydroPropMap;
96  
97      Molecule* mol;
98 <    StuntDouble* integrableObject;
98 >    StuntDouble* sd;
99      SimInfo::MoleculeIterator i;
100      Molecule::IntegrableObjectIterator  j;              
101      bool needHydroPropFile = false;
102      
103      for (mol = info->beginMolecule(i); mol != NULL;
104           mol = info->nextMolecule(i)) {
105 <      for (integrableObject = mol->beginIntegrableObject(j);
106 <           integrableObject != NULL;
107 <           integrableObject = mol->nextIntegrableObject(j)) {
105 >
106 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
107 >           sd = mol->nextIntegrableObject(j)) {
108          
109 <        if (integrableObject->isRigidBody()) {
110 <          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
109 >        if (sd->isRigidBody()) {
110 >          RigidBody* rb = static_cast<RigidBody*>(sd);
111            if (rb->getNumAtoms() > 1) needHydroPropFile = true;
112          }
113          
# Line 130 | Line 130 | namespace OpenMD {
130  
131        for (mol = info->beginMolecule(i); mol != NULL;
132             mol = info->nextMolecule(i)) {
133 –        for (integrableObject = mol->beginIntegrableObject(j);
134 –             integrableObject != NULL;
135 –             integrableObject = mol->nextIntegrableObject(j)) {
133  
134 <          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
134 >        for (sd = mol->beginIntegrableObject(j);  sd != NULL;
135 >             sd = mol->nextIntegrableObject(j)) {
136 >
137 >          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(sd->getType());
138            if (iter != hydroPropMap.end()) {
139              hydroProps_.push_back(iter->second);
140            } else {
141              sprintf( painCave.errMsg,
142 <                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
142 >                     "Can not find resistance tensor for atom [%s]\n", sd->getType().c_str());
143              painCave.severity = OPENMD_ERROR;
144              painCave.isFatal = 1;
145              simError();  
# Line 151 | Line 151 | namespace OpenMD {
151        std::map<std::string, HydroProp*> hydroPropMap;
152        for (mol = info->beginMolecule(i); mol != NULL;
153             mol = info->nextMolecule(i)) {
154 <        for (integrableObject = mol->beginIntegrableObject(j);
155 <             integrableObject != NULL;
156 <             integrableObject = mol->nextIntegrableObject(j)) {
154 >
155 >        for (sd = mol->beginIntegrableObject(j); sd != NULL;
156 >             sd = mol->nextIntegrableObject(j)) {
157 >
158            Shape* currShape = NULL;
159  
160 <          if (integrableObject->isAtom()){
161 <            Atom* atom = static_cast<Atom*>(integrableObject);
160 >          if (sd->isAtom()){
161 >            Atom* atom = static_cast<Atom*>(sd);
162              AtomType* atomType = atom->getAtomType();
163              GayBerneAdapter gba = GayBerneAdapter(atomType);
164              if (gba.isGayBerne()) {
# Line 199 | Line 200 | namespace OpenMD {
200  
201  
202            HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
203 <          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
203 >          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(sd->getType());
204            if (iter != hydroPropMap.end())
205              hydroProps_.push_back(iter->second);
206            else {
207              currHydroProp->complete();
208 <            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
208 >            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(sd->getType(), currHydroProp));
209              hydroProps_.push_back(currHydroProp);
210            }
211          }
# Line 234 | Line 235 | namespace OpenMD {
235      SimInfo::MoleculeIterator i;
236      Molecule::IntegrableObjectIterator  j;
237      Molecule* mol;
238 <    StuntDouble* integrableObject;
238 >    StuntDouble* sd;
239      RealType mass;
240      Vector3d pos;
241      Vector3d frc;
# Line 271 | Line 272 | namespace OpenMD {
272          }
273        }
274        
275 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
276 <           integrableObject = mol->nextIntegrableObject(j)) {
275 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
276 >           sd = mol->nextIntegrableObject(j)) {
277            
278          if (freezeMolecule)
279 <          fdf += integrableObject->freeze();
279 >          fdf += sd->freeze();
280          
281          if (doLangevinForces) {  
282 <          mass = integrableObject->getMass();
283 <          if (integrableObject->isDirectional()){
282 >          mass = sd->getMass();
283 >          if (sd->isDirectional()){
284  
285              // preliminaries for directional objects:
286  
287 <            A = integrableObject->getA();
287 >            A = sd->getA();
288              Atrans = A.transpose();
289              Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
290  
# Line 294 | Line 295 | namespace OpenMD {
295              genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
296              Vector3d randomForceLab = Atrans * randomForceBody;
297              Vector3d randomTorqueLab = Atrans * randomTorqueBody;
298 <            integrableObject->addFrc(randomForceLab);            
299 <            integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
298 >            sd->addFrc(randomForceLab);            
299 >            sd->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
300  
301 <            Mat3x3d I = integrableObject->getI();
301 >            Mat3x3d I = sd->getI();
302              Vector3d omegaBody;
303  
304              // What remains contains velocity explicitly, but the velocity required
# Line 307 | Line 308 | namespace OpenMD {
308  
309              // this is the velocity at the half-step:
310              
311 <            Vector3d vel =integrableObject->getVel();
312 <            Vector3d angMom = integrableObject->getJ();
311 >            Vector3d vel =sd->getVel();
312 >            Vector3d angMom = sd->getJ();
313  
314              //estimate velocity at full-step using everything but friction forces:          
315  
316 <            frc = integrableObject->getFrc();
316 >            frc = sd->getFrc();
317              Vector3d velStep = vel + (dt2_ /mass * PhysicalConstants::energyConvert) * frc;
318  
319 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
319 >            Tb = sd->lab2Body(sd->getTrq());
320              Vector3d angMomStep = angMom + (dt2_ * PhysicalConstants::energyConvert) * Tb;                            
321  
322              Vector3d omegaLab;
# Line 334 | Line 335 | namespace OpenMD {
335  
336              for (int k = 0; k < maxIterNum_; k++) {
337                              
338 <              if (integrableObject->isLinear()) {
339 <                int linearAxis = integrableObject->linearAxis();
338 >              if (sd->isLinear()) {
339 >                int linearAxis = sd->linearAxis();
340                  int l = (linearAxis +1 )%3;
341                  int m = (linearAxis +2 )%3;
342                  omegaBody[l] = angMomStep[l] /I(l, l);
# Line 374 | Line 375 | namespace OpenMD {
375                  break; // iteration ends here
376              }
377  
378 <            integrableObject->addFrc(frictionForceLab);
379 <            integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
378 >            sd->addFrc(frictionForceLab);
379 >            sd->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
380  
381              
382            } else {
# Line 384 | Line 385 | namespace OpenMD {
385              Vector3d randomForce;
386              Vector3d randomTorque;
387              genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
388 <            integrableObject->addFrc(randomForce);            
388 >            sd->addFrc(randomForce);            
389  
390              // What remains contains velocity explicitly, but the velocity required
391              // is at the full step: v(t + h), while we have initially the velocity
# Line 393 | Line 394 | namespace OpenMD {
394  
395              // this is the velocity at the half-step:
396              
397 <            Vector3d vel =integrableObject->getVel();
397 >            Vector3d vel =sd->getVel();
398  
399              //estimate velocity at full-step using everything but friction forces:          
400  
401 <            frc = integrableObject->getFrc();
401 >            frc = sd->getFrc();
402              Vector3d velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * frc;
403  
404              Vector3d frictionForce(0.0);
# Line 423 | Line 424 | namespace OpenMD {
424                  break; // iteration ends here
425              }
426  
427 <            integrableObject->addFrc(frictionForce);
427 >            sd->addFrc(frictionForce);
428  
429            }
430          }

Diff Legend

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