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

Comparing trunk/OOPSE-4/src/integrators/SMIPDForceManager.cpp (file contents):
Revision 3458 by chuckv, Fri Oct 3 17:47:08 2008 UTC vs.
Revision 3473 by chuckv, Fri Nov 14 15:44:34 2008 UTC

# Line 58 | Line 58 | namespace oopse {
58      
59      // Create Hull, Convex Hull for now, other options later.
60      surfaceMesh_ = new ConvexHull();
61 +
62      
63      
64      /* Check that the simulation has target pressure and target
# Line 79 | Line 80 | namespace oopse {
80        painCave.isFatal = 1;
81        simError();
82      } else {
83 <      targetPressure_ = simParams->getTargetPressure();
83 >      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
84      }
85  
86    
# Line 91 | Line 92 | namespace oopse {
92        simError();
93      }
94  
95 +    if (!simParams->haveViscosity()) {
96 +      sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
97 +      painCave.isFatal = 1;
98 +      painCave.severity = OOPSE_ERROR;
99 +      simError();
100 +    }
101  
102 +
103 +    
104 +
105 +    //Compute initial hull
106 +    /*
107 +    surfaceMesh_->computeHull(localSites_);
108 +    Area0_ = surfaceMesh_->getArea();
109 +    int nTriangles = surfaceMesh_->getNMeshElements();
110 +    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
111 +    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
112 +      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
113 +    //RealType eta0 = gamma_0/
114 +    */
115 +
116      // Build the hydroProp map:
117      std::map<std::string, HydroProp*> hydroPropMap;
118  
# Line 121 | Line 142 | namespace oopse {
142          hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
143        } else {              
144          sprintf( painCave.errMsg,
145 <                 "HydroPropFile must be set to a file name if Langevin Dynamics\n"
145 >                 "HydroPropFile must be set to a file name if SMIPDynamics\n"
146                   "\tis specified for rigidBodies which contain more than one atom\n"
147 <                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n");
147 >                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
148 >                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
149 >                 "\tset to 1.0 because the friction and random forces will be\n"
150 >                 "\tdynamically re-set assuming this is true.\n");
151          painCave.severity = OOPSE_ERROR;
152          painCave.isFatal = 1;
153          simError();  
# Line 192 | Line 216 | namespace oopse {
216                    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
217                    if (ljData != NULL) {
218                      LJParam ljParam = ljData->getData();
219 <                    currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
219 >                    currShape = new Sphere(atom->getPos(), 2.0);
220                    } else {
221                      sprintf( painCave.errMsg,
222                               "Can not cast GenericData to LJParam\n");
# Line 204 | Line 228 | namespace oopse {
228                } else {
229                  int aNum = etab.GetAtomicNum((atom->getType()).c_str());
230                  if (aNum != 0) {
231 <                  currShape = new Sphere(atom->getPos(), etab.GetVdwRad(aNum));
231 >                  currShape = new Sphere(atom->getPos(), 2.0);
232                  } else {
233                    sprintf( painCave.errMsg,
234                             "Could not find atom type in default element.txt\n");
# Line 215 | Line 239 | namespace oopse {
239                }
240              }
241            }
242 <          HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
242 >          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
243            std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
244            if (iter != hydroPropMap.end())
245              hydroProps_.push_back(iter->second);
# Line 230 | Line 254 | namespace oopse {
254  
255      /* Compute hull first time through to get the area of t=0*/
256  
257 <    /* Build a vector of integrable objects to determine if the are surface atoms */
257 >    //Build a vector of integrable objects to determine if the are surface atoms
258      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
259        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
260             integrableObject = mol->nextIntegrableObject(j)) {  
# Line 238 | Line 262 | namespace oopse {
262        }
263      }
264  
265 <    surfaceMesh_->computeHull(localSites_);
242 <    Area0_ = surfaceMesh_->getArea();
243 <    //variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
244 <    
265 >
266    }  
267  
268    std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
# Line 282 | Line 303 | namespace oopse {
303      surfaceMesh_->computeHull(localSites_);
304  
305      /* Get area and number of surface stunt doubles and compute new variance */
306 <     //RealType area = surfaceMesh_->getArea();
307 <     //RealType nSurfaceSDs = surfaceMesh_->getNs();
306 >     RealType area = surfaceMesh_->getArea();
307 >     int nSurfaceSDs = surfaceMesh_->getNs();
308  
288 – //    std::cerr << "Surface Area is: " << area << " nSurfaceSDs is: " << nSurfaceSDs << std::endl;
289 –
290 –    /* Compute variance for random forces */
309      
310 < //    variance_ = sqrt(2.0*NumericConstant::PI)*(targetPressure_*area/nSurfaceSDs);
311 <    
312 < //    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
313 < //    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
314 <  
310 >    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
311 >    int nTriangles = sMesh.size();
312 >
313 >
314 >
315 >     /* Compute variance for random forces */
316 >  
317 >
318 >
319 >
320 >    std::vector<RealType>  randNums = genTriangleForces(nTriangles, 1.0);
321 >
322      /* Loop over the mesh faces and apply random force to each of the faces*/
323      
324      
325 < //    std::vector<Triangle*>::iterator face;
326 < //    std::vector<StuntDouble*>::iterator vertex;
327 < /*    
325 >    std::vector<Triangle>::iterator face;
326 >    std::vector<StuntDouble*>::iterator vertex;
327 >    int thisNumber = 0;
328      for (face = sMesh.begin(); face != sMesh.end(); ++face){
329      
330 <      Triangle* thisTriangle = *face;
331 <      std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
330 >      Triangle thisTriangle = *face;
331 >      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
332 >      RealType thisArea = thisTriangle.getArea();
333 >      // RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*thisArea) /OOPSEConstant::energyConvert;
334 >      // gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());
335  
336 +      /* Get Random Force */
337 +      Vector3d unitNormal = thisTriangle.getNormal();
338 +      unitNormal.normalize();
339 +      //Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
340 +      Vector3d centroid = thisTriangle.getCentroid();
341 +      Vector3d facetVel = thisTriangle.getFacetVelocity();
342 +      RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/3.14;
343 +
344 +      RealType f_normal = simParams->getViscosity()*hydroLength*1.439326479e4*OOPSEConstant::energyConvert;
345 +      RealType extPressure = -(targetPressure_ * thisArea);
346 +      RealType randomForce = randNums[thisNumber] * f_normal * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
347 +      RealType dragForce = -f_normal * dot(facetVel, unitNormal);      
348 +      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
349 +      
350 +      //      Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
351 +      
352 +      //std::cout << "randomForce " << randomForce << " dragForce " << dragForce <<  " hydro  " << hydroLength << std::endl;
353 +
354 +
355        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
356 <            std::cout << (*vertex)->getPos() << std::endl;
357 <         // mass = integrableObject->getMass();
356 >        if ((*vertex) != NULL){
357 >          // mass = integrableObject->getMass();
358 >          Vector3d vertexForce = langevinForce/3.0;
359 >          (*vertex)->addFrc(vertexForce);
360  
361 <         //      integrableObject->addFrc(randomForce);          
361 >          if ((*vertex)->isDirectional()){
362 >
363 >            Vector3d vertexPos = (*vertex)->getPos();
364 >            Vector3d vertexCentroidVector = vertexPos - centroid;
365 >            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
366 >          }
367 >        }  
368        }
369 +    }
370  
371 <
316 <    }
317 <  */
371 >    /* Now loop over all surface particles and apply the drag*/
372      /*
373 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
373 >    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
374 >    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
375 >      integrableObject = *vertex;
376 >      mass = integrableObject->getMass();
377 >      if (integrableObject->isDirectional()){
378 >        
379 >        // preliminaries for directional objects:
380 >        
381 >        A = integrableObject->getA();
382 >        Atrans = A.transpose();
383 >        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
384 >        //apply random force and torque at center of resistance
385 >        Mat3x3d I = integrableObject->getI();
386 >        Vector3d omegaBody;
387 >        
388 >        // What remains contains velocity explicitly, but the velocity required
389 >        // is at the full step: v(t + h), while we have initially the velocity
390 >        // at the half step: v(t + h/2).  We need to iterate to converge the
391 >        // friction force and friction torque vectors.
392 >        
393 >        // this is the velocity at the half-step:
394 >        
395 >        Vector3d vel =integrableObject->getVel();
396 >        Vector3d angMom = integrableObject->getJ();
397 >        
398 >        //estimate velocity at full-step using everything but friction forces:          
399 >        
400 >        frc = integrableObject->getFrc();
401 >        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
402 >        
403 >        Tb = integrableObject->lab2Body(integrableObject->getTrq());
404 >        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
405 >        
406 >        Vector3d omegaLab;
407 >        Vector3d vcdLab;
408 >        Vector3d vcdBody;
409 >        Vector3d frictionForceBody;
410 >        Vector3d frictionForceLab(0.0);
411 >        Vector3d oldFFL;  // used to test for convergence
412 >        Vector3d frictionTorqueBody(0.0);
413 >        Vector3d oldFTB;  // used to test for convergence
414 >        Vector3d frictionTorqueLab;
415 >        RealType fdot;
416 >        RealType tdot;
417  
418 <      
419 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
420 <           integrableObject = mol->nextIntegrableObject(j)) {
421 <          
422 <          mass = integrableObject->getMass();
423 <          if (integrableObject->isDirectional()){
424 <
425 <            // preliminaries for directional objects:
426 <
427 <            A = integrableObject->getA();
428 <            Atrans = A.transpose();
429 <            Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
430 <
431 <            //apply random force and torque at center of resistance
432 <
433 <            Vector3d randomForceBody;
434 <            Vector3d randomTorqueBody;
435 <            genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
436 <            Vector3d randomForceLab = Atrans * randomForceBody;
437 <            Vector3d randomTorqueLab = Atrans * randomTorqueBody;
438 <            integrableObject->addFrc(randomForceLab);            
439 <            integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
440 <
441 <            Mat3x3d I = integrableObject->getI();
442 <            Vector3d omegaBody;
443 <
444 <            // What remains contains velocity explicitly, but the velocity required
445 <            // is at the full step: v(t + h), while we have initially the velocity
446 <            // at the half step: v(t + h/2).  We need to iterate to converge the
447 <            // friction force and friction torque vectors.
448 <
352 <            // this is the velocity at the half-step:
353 <            
354 <            Vector3d vel =integrableObject->getVel();
355 <            Vector3d angMom = integrableObject->getJ();
356 <
357 <            //estimate velocity at full-step using everything but friction forces:          
358 <
359 <            frc = integrableObject->getFrc();
360 <            Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
361 <
362 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
363 <            Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
364 <
365 <            Vector3d omegaLab;
366 <            Vector3d vcdLab;
367 <            Vector3d vcdBody;
368 <            Vector3d frictionForceBody;
369 <            Vector3d frictionForceLab(0.0);
370 <            Vector3d oldFFL;  // used to test for convergence
371 <            Vector3d frictionTorqueBody(0.0);
372 <            Vector3d oldFTB;  // used to test for convergence
373 <            Vector3d frictionTorqueLab;
374 <            RealType fdot;
375 <            RealType tdot;
376 <
377 <            //iteration starts here:
378 <
379 <            for (int k = 0; k < maxIterNum_; k++) {
380 <                            
381 <              if (integrableObject->isLinear()) {
382 <                int linearAxis = integrableObject->linearAxis();
383 <                int l = (linearAxis +1 )%3;
384 <                int m = (linearAxis +2 )%3;
385 <                omegaBody[l] = angMomStep[l] /I(l, l);
386 <                omegaBody[m] = angMomStep[m] /I(m, m);
387 <                
388 <              } else {
389 <                omegaBody[0] = angMomStep[0] /I(0, 0);
390 <                omegaBody[1] = angMomStep[1] /I(1, 1);
391 <                omegaBody[2] = angMomStep[2] /I(2, 2);
392 <              }
418 >        //iteration starts here:
419 >        
420 >        for (int k = 0; k < maxIterNum_; k++) {
421 >          
422 >          if (integrableObject->isLinear()) {
423 >            int linearAxis = integrableObject->linearAxis();
424 >            int l = (linearAxis +1 )%3;
425 >            int m = (linearAxis +2 )%3;
426 >            omegaBody[l] = angMomStep[l] /I(l, l);
427 >            omegaBody[m] = angMomStep[m] /I(m, m);
428 >            
429 >          } else {
430 >            omegaBody[0] = angMomStep[0] /I(0, 0);
431 >            omegaBody[1] = angMomStep[1] /I(1, 1);
432 >            omegaBody[2] = angMomStep[2] /I(2, 2);
433 >          }
434 >          
435 >          omegaLab = Atrans * omegaBody;
436 >          
437 >          // apply friction force and torque at center of resistance
438 >          
439 >          vcdLab = velStep + cross(omegaLab, rcrLab);      
440 >          vcdBody = A * vcdLab;
441 >          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
442 >          oldFFL = frictionForceLab;
443 >          frictionForceLab = Atrans * frictionForceBody;
444 >          oldFTB = frictionTorqueBody;
445 >          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
446 >          frictionTorqueLab = Atrans * frictionTorqueBody;
447 >          
448 >          // re-estimate velocities at full-step using friction forces:
449                
450 <              omegaLab = Atrans * omegaBody;
451 <              
396 <              // apply friction force and torque at center of resistance
397 <              
398 <              vcdLab = velStep + cross(omegaLab, rcrLab);      
399 <              vcdBody = A * vcdLab;
400 <              frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
401 <              oldFFL = frictionForceLab;
402 <              frictionForceLab = Atrans * frictionForceBody;
403 <              oldFTB = frictionTorqueBody;
404 <              frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
405 <              frictionTorqueLab = Atrans * frictionTorqueBody;
406 <              
407 <              // re-estimate velocities at full-step using friction forces:
408 <              
409 <              velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
410 <              angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
450 >          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
451 >          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
452  
453 <              // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
453 >          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
454                
455 <              fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
456 <              tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
457 <              
458 <              if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
459 <                break; // iteration ends here
460 <            }
455 >          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
456 >          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
457 >          
458 >          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
459 >            break; // iteration ends here
460 >        }
461 >        
462 >        integrableObject->addFrc(frictionForceLab);
463 >        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
464  
421 –            integrableObject->addFrc(frictionForceLab);
422 –            integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
423 –
465              
466 <          } else {
467 <            //spherical atom
466 >      } else {
467 >        //spherical atom
468  
469 <            Vector3d randomForce;
470 <            Vector3d randomTorque;
471 <            genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
472 <            integrableObject->addFrc(randomForce);            
473 <
474 <            // What remains contains velocity explicitly, but the velocity required
475 <            // is at the full step: v(t + h), while we have initially the velocity
476 <            // at the half step: v(t + h/2).  We need to iterate to converge the
477 <            // friction force vector.
478 <
479 <            // this is the velocity at the half-step:
480 <            
481 <            Vector3d vel =integrableObject->getVel();
482 <
483 <            //estimate velocity at full-step using everything but friction forces:          
484 <
485 <            frc = integrableObject->getFrc();
486 <            Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
487 <
488 <            Vector3d frictionForce(0.0);
489 <            Vector3d oldFF;  // used to test for convergence
490 <            RealType fdot;
491 <
492 <            //iteration starts here:
493 <
494 <            for (int k = 0; k < maxIterNum_; k++) {
454 <
455 <              oldFF = frictionForce;                            
456 <              frictionForce = -hydroProps_[index]->getXitt() * velStep;
457 <
458 <              // re-estimate velocities at full-step using friction forces:
459 <              
460 <              velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
461 <
462 <              // check for convergence (if the vector has converged, fdot will be 1.0):
463 <              
464 <              fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
465 <              
466 <              if (fabs(1.0 - fdot) <= forceTolerance_)
467 <                break; // iteration ends here
468 <            }
469 <
470 <            integrableObject->addFrc(frictionForce);
471 <
469 >        // What remains contains velocity explicitly, but the velocity required
470 >        // is at the full step: v(t + h), while we have initially the velocity
471 >        // at the half step: v(t + h/2).  We need to iterate to converge the
472 >        // friction force vector.
473 >        
474 >        // this is the velocity at the half-step:
475 >        
476 >        Vector3d vel =integrableObject->getVel();
477 >        
478 >        //estimate velocity at full-step using everything but friction forces:          
479 >        
480 >        frc = integrableObject->getFrc();
481 >        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
482 >        
483 >        Vector3d frictionForce(0.0);
484 >        Vector3d oldFF;  // used to test for convergence
485 >        RealType fdot;
486 >        
487 >        //iteration starts here:
488 >        
489 >        for (int k = 0; k < maxIterNum_; k++) {
490 >          
491 >          oldFF = frictionForce;                            
492 >          frictionForce = -hydroProps_[index]->getXitt() * velStep;
493 >          //frictionForce = -gamma_t*velStep;
494 >          // re-estimate velocities at full-step using friction forces:
495            
496 <          }
496 >          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
497 >          
498 >          // check for convergence (if the vector has converged, fdot will be 1.0):
499            
500 <        ++index;
501 <    
500 >          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
501 >          
502 >          if (fabs(1.0 - fdot) <= forceTolerance_)
503 >            break; // iteration ends here
504 >        }
505 >        
506 >        integrableObject->addFrc(frictionForce);
507 >        
508 >        
509        }
510 <    }    
510 >  
511 >      
512 >  }
513      */
514 <    // info_->setFdf(fdf);
515 <    // veloMunge->removeComDrift();
516 <    // Remove angular drift if we are not using periodic boundary conditions.
483 <    //if(!simParams->getUsePeriodicBoundaryConditions())
484 <    //  veloMunge->removeAngularDrift();
485 <
486 <    //ForceManager::postCalculation(needStress);  
514 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
515 >    currSnapshot->setVolume(surfaceMesh_->getVolume());    
516 >    ForceManager::postCalculation(needStress);  
517    }
518  
519 < void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
519 >  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
520  
521 <
521 >    
522      Vector<RealType, 6> Z;
523      Vector<RealType, 6> generalForce;
524 <        
524 >    
525 >
526      Z[0] = randNumGen_.randNorm(0, variance);
527      Z[1] = randNumGen_.randNorm(0, variance);
528      Z[2] = randNumGen_.randNorm(0, variance);
# Line 521 | Line 552 | void SMIPDForceManager::genRandomForceAndTorque(Vector
552      if (worldRank == 0) {
553   #endif
554        for (int i = 0; i < nTriangles; i++) {
555 <        gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
555 >        //gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
556 >        gaussRand[i] = randNumGen_.randNorm(0.0, 1.0);    
557        }
558   #ifdef IS_MPI
559      }
# Line 536 | Line 568 | void SMIPDForceManager::genRandomForceAndTorque(Vector
568        MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
569      }
570   #endif
571 <
571 >  
572      for (int i = 0; i < nTriangles; i++) {
573        gaussRand[i] = gaussRand[i] * variance;
574      }
575 <
575 >  
576      return gaussRand;
577    }
578  
579 +
580 +
581 +
582 +
583   }

Diff Legend

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