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 3469 by chuckv, Wed Oct 22 17:52:40 2008 UTC vs.
Revision 3515 by gezelter, Mon Jul 20 17:50:40 2009 UTC

# Line 41 | Line 41
41   #include <fstream>
42   #include <iostream>
43   #include "integrators/SMIPDForceManager.hpp"
44 – #include "math/CholeskyDecomposition.hpp"
44   #include "utils/OOPSEConstant.hpp"
46 – #include "hydrodynamics/Sphere.hpp"
47 – #include "hydrodynamics/Ellipsoid.hpp"
48 – #include "utils/ElementsTable.hpp"
45   #include "math/ConvexHull.hpp"
46   #include "math/Triangle.hpp"
47  
52 –
48   namespace oopse {
49  
50 <  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
50 >  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
51 >
52      simParams = info->getSimParams();
53 +    thermo = new Thermo(info);
54      veloMunge = new Velocitizer(info);
55      
56      // Create Hull, Convex Hull for now, other options later.
57 +    
58      surfaceMesh_ = new ConvexHull();
61 –
59      
63 –    
60      /* Check that the simulation has target pressure and target
61 <       temperature set*/
62 <
61 >       temperature set */
62 >    
63      if (!simParams->haveTargetTemp()) {
64 <      sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
64 >      sprintf(painCave.errMsg,
65 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
66 >              "   without a targetTemp!\n");      
67        painCave.isFatal = 1;
68        painCave.severity = OOPSE_ERROR;
69        simError();
70      } else {
71        targetTemp_ = simParams->getTargetTemp();
72      }
73 <
73 >    
74      if (!simParams->haveTargetPressure()) {
75 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
76 <              "   without a targetPressure!\n");
77 <      
75 >      sprintf(painCave.errMsg,
76 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
77 >              "   without a targetPressure!\n");      
78        painCave.isFatal = 1;
79        simError();
80      } else {
81 <      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
81 >      // Convert pressure from atm -> amu/(fs^2*Ang)
82 >      targetPressure_ = simParams->getTargetPressure() /
83 >        OOPSEConstant::pressureConvert;
84      }
85 –
85    
86      if (simParams->getUsePeriodicBoundaryConditions()) {
87 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
88 <              "   with periodic boundary conditions !\n");
89 <      
87 >      sprintf(painCave.errMsg,
88 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
89 >              "   with periodic boundary conditions!\n");    
90        painCave.isFatal = 1;
91        simError();
92      }
93 +    
94 +    if (!simParams->haveThermalConductivity()) {
95 +      sprintf(painCave.errMsg,
96 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
97 +              "   without a thermalConductivity!\n");
98 +      painCave.isFatal = 1;
99 +      painCave.severity = OOPSE_ERROR;
100 +      simError();
101 +    }else{
102 +      thermalConductivity_ = simParams->getThermalConductivity();
103 +    }
104  
105 <
106 <
105 >    if (!simParams->haveThermalLength()) {
106 >      sprintf(painCave.errMsg,
107 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
108 >              "   without a thermalLength!\n");
109 >      painCave.isFatal = 1;
110 >      painCave.severity = OOPSE_ERROR;
111 >      simError();
112 >    }else{
113 >      thermalLength_ = simParams->getThermalLength();
114 >    }
115      
116 +    dt_ = simParams->getDt();
117 +  
118 +    variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;
119  
120 <    //Compute initial hull
121 <    /*
101 <    surfaceMesh_->computeHull(localSites_);
102 <    Area0_ = surfaceMesh_->getArea();
103 <    int nTriangles = surfaceMesh_->getNMeshElements();
104 <    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
105 <    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
106 <      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
107 <    //RealType eta0 = gamma_0/
108 <    */
109 <
110 <    // Build the hydroProp map:
111 <    std::map<std::string, HydroProp*> hydroPropMap;
112 <
120 >    // Build a vector of integrable objects to determine if the are
121 >    // surface atoms
122      Molecule* mol;
123      StuntDouble* integrableObject;
124      SimInfo::MoleculeIterator i;
125 <    Molecule::IntegrableObjectIterator  j;              
126 <    bool needHydroPropFile = false;
127 <    
128 <    for (mol = info->beginMolecule(i); mol != NULL;
120 <         mol = info->nextMolecule(i)) {
125 >    Molecule::IntegrableObjectIterator  j;
126 >
127 >    for (mol = info_->beginMolecule(i); mol != NULL;
128 >         mol = info_->nextMolecule(i)) {          
129        for (integrableObject = mol->beginIntegrableObject(j);
130             integrableObject != NULL;
123 –           integrableObject = mol->nextIntegrableObject(j)) {
124 –        
125 –        if (integrableObject->isRigidBody()) {
126 –          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
127 –          if (rb->getNumAtoms() > 1) needHydroPropFile = true;
128 –        }
129 –        
130 –      }
131 –    }
132 –        
133 –
134 –    if (needHydroPropFile) {              
135 –      if (simParams->haveHydroPropFile()) {
136 –        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
137 –      } else {              
138 –        sprintf( painCave.errMsg,
139 –                 "HydroPropFile must be set to a file name if SMIPDynamics\n"
140 –                 "\tis specified for rigidBodies which contain more than one atom\n"
141 –                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
142 –                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
143 –                 "\tset to 1.0 because the friction and random forces will be\n"
144 –                 "\tdynamically re-set assuming this is true.\n");
145 –        painCave.severity = OOPSE_ERROR;
146 –        painCave.isFatal = 1;
147 –        simError();  
148 –      }      
149 –
150 –      for (mol = info->beginMolecule(i); mol != NULL;
151 –           mol = info->nextMolecule(i)) {
152 –        for (integrableObject = mol->beginIntegrableObject(j);
153 –             integrableObject != NULL;
154 –             integrableObject = mol->nextIntegrableObject(j)) {
155 –
156 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
157 –          if (iter != hydroPropMap.end()) {
158 –            hydroProps_.push_back(iter->second);
159 –          } else {
160 –            sprintf( painCave.errMsg,
161 –                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
162 –            painCave.severity = OOPSE_ERROR;
163 –            painCave.isFatal = 1;
164 –            simError();  
165 –          }        
166 –        }
167 –      }
168 –    } else {
169 –      
170 –      std::map<std::string, HydroProp*> hydroPropMap;
171 –      for (mol = info->beginMolecule(i); mol != NULL;
172 –           mol = info->nextMolecule(i)) {
173 –        for (integrableObject = mol->beginIntegrableObject(j);
174 –             integrableObject != NULL;
175 –             integrableObject = mol->nextIntegrableObject(j)) {
176 –          Shape* currShape = NULL;
177 –
178 –          if (integrableObject->isAtom()){
179 –            Atom* atom = static_cast<Atom*>(integrableObject);
180 –            AtomType* atomType = atom->getAtomType();
181 –            if (atomType->isGayBerne()) {
182 –              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
183 –              GenericData* data = dAtomType->getPropertyByName("GayBerne");
184 –              if (data != NULL) {
185 –                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
186 –                
187 –                if (gayBerneData != NULL) {  
188 –                  GayBerneParam gayBerneParam = gayBerneData->getData();
189 –                  currShape = new Ellipsoid(V3Zero,
190 –                                            gayBerneParam.GB_l / 2.0,
191 –                                            gayBerneParam.GB_d / 2.0,
192 –                                            Mat3x3d::identity());
193 –                } else {
194 –                  sprintf( painCave.errMsg,
195 –                           "Can not cast GenericData to GayBerneParam\n");
196 –                  painCave.severity = OOPSE_ERROR;
197 –                  painCave.isFatal = 1;
198 –                  simError();  
199 –                }
200 –              } else {
201 –                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
202 –                painCave.severity = OOPSE_ERROR;
203 –                painCave.isFatal = 1;
204 –                simError();    
205 –              }
206 –            } else {
207 –              if (atomType->isLennardJones()){
208 –                GenericData* data = atomType->getPropertyByName("LennardJones");
209 –                if (data != NULL) {
210 –                  LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
211 –                  if (ljData != NULL) {
212 –                    LJParam ljParam = ljData->getData();
213 –                    currShape = new Sphere(atom->getPos(), 2.0);
214 –                  } else {
215 –                    sprintf( painCave.errMsg,
216 –                             "Can not cast GenericData to LJParam\n");
217 –                    painCave.severity = OOPSE_ERROR;
218 –                    painCave.isFatal = 1;
219 –                    simError();          
220 –                  }      
221 –                }
222 –              } else {
223 –                int aNum = etab.GetAtomicNum((atom->getType()).c_str());
224 –                if (aNum != 0) {
225 –                  currShape = new Sphere(atom->getPos(), 2.0);
226 –                } else {
227 –                  sprintf( painCave.errMsg,
228 –                           "Could not find atom type in default element.txt\n");
229 –                  painCave.severity = OOPSE_ERROR;
230 –                  painCave.isFatal = 1;
231 –                  simError();          
232 –                }
233 –              }
234 –            }
235 –          }
236 –          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
237 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
238 –          if (iter != hydroPropMap.end())
239 –            hydroProps_.push_back(iter->second);
240 –          else {
241 –            currHydroProp->complete();
242 –            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
243 –            hydroProps_.push_back(currHydroProp);
244 –          }
245 –        }
246 –      }
247 –    }
248 –
249 –    /* Compute hull first time through to get the area of t=0*/
250 –
251 –    //Build a vector of integrable objects to determine if the are surface atoms
252 –    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
253 –      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
131             integrableObject = mol->nextIntegrableObject(j)) {  
132          localSites_.push_back(integrableObject);
133        }
134 <    }
258 <
259 <
134 >    }  
135    }  
261 –
262 –  std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
263 –    std::map<std::string, HydroProp*> props;
264 –    std::ifstream ifs(filename.c_str());
265 –    if (ifs.is_open()) {
266 –      
267 –    }
268 –    
269 –    const unsigned int BufferSize = 65535;
270 –    char buffer[BufferSize];  
271 –    while (ifs.getline(buffer, BufferSize)) {
272 –      HydroProp* currProp = new HydroProp(buffer);
273 –      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
274 –    }
275 –
276 –    return props;
277 –  }
136    
137    void SMIPDForceManager::postCalculation(bool needStress){
138      SimInfo::MoleculeIterator i;
139      Molecule::IntegrableObjectIterator  j;
140      Molecule* mol;
141      StuntDouble* integrableObject;
284 –    RealType mass;
285 –    Vector3d pos;
286 –    Vector3d frc;
287 –    Mat3x3d A;
288 –    Mat3x3d Atrans;
289 –    Vector3d Tb;
290 –    Vector3d ji;
291 –    unsigned int index = 0;
292 –    int fdf;
293 –  
294 –    fdf = 0;
142    
143 <    /*Compute surface Mesh*/
143 >    // Compute surface Mesh
144      surfaceMesh_->computeHull(localSites_);
145  
146 <    /* Get area and number of surface stunt doubles and compute new variance */
147 <     RealType area = surfaceMesh_->getArea();
148 <     int nSurfaceSDs = surfaceMesh_->getNs();
302 <
303 <    
146 >    // Get total area and number of surface stunt doubles
147 >    RealType area = surfaceMesh_->getArea();
148 >    int nSurfaceSDs = surfaceMesh_->getNs();
149      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
150      int nTriangles = sMesh.size();
151  
152 +    // Generate all of the necessary random forces
153 +    std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
154  
155 +    RealType instaTemp = thermo->getTemperature();
156  
157 <     /* Compute variance for random forces */
158 <  
311 <    RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*area/nTriangles)
312 <       /OOPSEConstant::energyConvert;
313 <
314 <    gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * area * area * simParams->getDt()) /
315 <      (4.0 * nTriangles * nTriangles* OOPSEConstant::kB*simParams->getTargetTemp()) /OOPSEConstant::energyConvert;
316 <
317 <    std::vector<RealType>  randNums = genTriangleForces(nTriangles, sigma_t);
318 <
319 <    /* Loop over the mesh faces and apply random force to each of the faces*/
320 <    
321 <    
157 >    // Loop over the mesh faces and apply external pressure to each
158 >    // of the faces
159      std::vector<Triangle>::iterator face;
160      std::vector<StuntDouble*>::iterator vertex;
161 <    int thisNumber = 0;
161 >    int thisFacet = 0;
162      for (face = sMesh.begin(); face != sMesh.end(); ++face){
163      
164        Triangle thisTriangle = *face;
165        std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
166 <      
330 <      /* Get Random Force */
166 >      RealType thisArea = thisTriangle.getArea();
167        Vector3d unitNormal = thisTriangle.getNormal();
168        unitNormal.normalize();
333 –      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
169        Vector3d centroid = thisTriangle.getCentroid();
170 +      Vector3d facetVel = thisTriangle.getFacetVelocity();
171 +      RealType thisMass = thisTriangle.getFacetMass();
172  
173 <      Vector3d langevinForce = randomForce - gamma_t_*thisTriangle.getFacetVelocity();
173 >      // gamma is the drag coefficient normal to the face of the triangle      
174 >      RealType gamma = thermalConductivity_ * thisMass * thisArea  
175 >        / (2.0 * thermalLength_ * OOPSEConstant::kb);
176        
177 <      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
178 <        if ((*vertex) != NULL){
179 <          // mass = integrableObject->getMass();
180 <          Vector3d vertexForce = langevinForce/3.0;
342 <          (*vertex)->addFrc(vertexForce);
177 >      gamma *= fabs(1.0 - targetTemp_/instaTemp);      
178 >      
179 >      RealType extPressure = - (targetPressure_ * thisArea) /
180 >        OOPSEConstant::energyConvert;
181  
182 <          if ((*vertex)->isDirectional()){
182 >      RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
183 >      RealType dragForce = -gamma * dot(facetVel, unitNormal);
184 >
185 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) *
186 >        unitNormal;
187 >      
188 >      // Apply triangle force to stuntdouble vertices
189 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
190 >        if ((*vertex) != NULL){
191 >          Vector3d vertexForce = langevinForce / 3.0;
192 >          (*vertex)->addFrc(vertexForce);          
193 >          if ((*vertex)->isDirectional()){          
194              Vector3d vertexPos = (*vertex)->getPos();
195              Vector3d vertexCentroidVector = vertexPos - centroid;
196              (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
# Line 349 | Line 198 | namespace oopse {
198          }  
199        }
200      }
352 –
353 –    /* Now loop over all surface particles and apply the drag*/
354 –    /*
355 –    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
356 –    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
357 –      integrableObject = *vertex;
358 –      mass = integrableObject->getMass();
359 –      if (integrableObject->isDirectional()){
360 –        
361 –        // preliminaries for directional objects:
362 –        
363 –        A = integrableObject->getA();
364 –        Atrans = A.transpose();
365 –        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
366 –        //apply random force and torque at center of resistance
367 –        Mat3x3d I = integrableObject->getI();
368 –        Vector3d omegaBody;
369 –        
370 –        // What remains contains velocity explicitly, but the velocity required
371 –        // is at the full step: v(t + h), while we have initially the velocity
372 –        // at the half step: v(t + h/2).  We need to iterate to converge the
373 –        // friction force and friction torque vectors.
374 –        
375 –        // this is the velocity at the half-step:
376 –        
377 –        Vector3d vel =integrableObject->getVel();
378 –        Vector3d angMom = integrableObject->getJ();
379 –        
380 –        //estimate velocity at full-step using everything but friction forces:          
381 –        
382 –        frc = integrableObject->getFrc();
383 –        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
384 –        
385 –        Tb = integrableObject->lab2Body(integrableObject->getTrq());
386 –        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
387 –        
388 –        Vector3d omegaLab;
389 –        Vector3d vcdLab;
390 –        Vector3d vcdBody;
391 –        Vector3d frictionForceBody;
392 –        Vector3d frictionForceLab(0.0);
393 –        Vector3d oldFFL;  // used to test for convergence
394 –        Vector3d frictionTorqueBody(0.0);
395 –        Vector3d oldFTB;  // used to test for convergence
396 –        Vector3d frictionTorqueLab;
397 –        RealType fdot;
398 –        RealType tdot;
399 –
400 –        //iteration starts here:
401 –        
402 –        for (int k = 0; k < maxIterNum_; k++) {
403 –          
404 –          if (integrableObject->isLinear()) {
405 –            int linearAxis = integrableObject->linearAxis();
406 –            int l = (linearAxis +1 )%3;
407 –            int m = (linearAxis +2 )%3;
408 –            omegaBody[l] = angMomStep[l] /I(l, l);
409 –            omegaBody[m] = angMomStep[m] /I(m, m);
410 –            
411 –          } else {
412 –            omegaBody[0] = angMomStep[0] /I(0, 0);
413 –            omegaBody[1] = angMomStep[1] /I(1, 1);
414 –            omegaBody[2] = angMomStep[2] /I(2, 2);
415 –          }
416 –          
417 –          omegaLab = Atrans * omegaBody;
418 –          
419 –          // apply friction force and torque at center of resistance
420 –          
421 –          vcdLab = velStep + cross(omegaLab, rcrLab);      
422 –          vcdBody = A * vcdLab;
423 –          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
424 –          oldFFL = frictionForceLab;
425 –          frictionForceLab = Atrans * frictionForceBody;
426 –          oldFTB = frictionTorqueBody;
427 –          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
428 –          frictionTorqueLab = Atrans * frictionTorqueBody;
429 –          
430 –          // re-estimate velocities at full-step using friction forces:
431 –              
432 –          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
433 –          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
434 –
435 –          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
436 –              
437 –          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
438 –          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
439 –          
440 –          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
441 –            break; // iteration ends here
442 –        }
443 –        
444 –        integrableObject->addFrc(frictionForceLab);
445 –        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
446 –
447 –            
448 –      } else {
449 –        //spherical atom
450 –
451 –        // What remains contains velocity explicitly, but the velocity required
452 –        // is at the full step: v(t + h), while we have initially the velocity
453 –        // at the half step: v(t + h/2).  We need to iterate to converge the
454 –        // friction force vector.
455 –        
456 –        // this is the velocity at the half-step:
457 –        
458 –        Vector3d vel =integrableObject->getVel();
459 –        
460 –        //estimate velocity at full-step using everything but friction forces:          
461 –        
462 –        frc = integrableObject->getFrc();
463 –        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
464 –        
465 –        Vector3d frictionForce(0.0);
466 –        Vector3d oldFF;  // used to test for convergence
467 –        RealType fdot;
468 –        
469 –        //iteration starts here:
470 –        
471 –        for (int k = 0; k < maxIterNum_; k++) {
472 –          
473 –          oldFF = frictionForce;                            
474 –          frictionForce = -hydroProps_[index]->getXitt() * velStep;
475 –          //frictionForce = -gamma_t*velStep;
476 –          // re-estimate velocities at full-step using friction forces:
477 –          
478 –          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
479 –          
480 –          // check for convergence (if the vector has converged, fdot will be 1.0):
481 –          
482 –          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
483 –          
484 –          if (fabs(1.0 - fdot) <= forceTolerance_)
485 –            break; // iteration ends here
486 –        }
487 –        
488 –        integrableObject->addFrc(frictionForce);
489 –        
490 –        
491 –      }
492 –  
493 –      
494 –  }
495 –    */
496 –    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
497 –    currSnapshot->setVolume(surfaceMesh_->getVolume());
201      
202 +    veloMunge->removeComDrift();
203 +    veloMunge->removeAngularDrift();
204 +    
205 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
206 +    currSnapshot->setVolume(surfaceMesh_->getVolume());    
207      ForceManager::postCalculation(needStress);  
208    }
209 <
210 <  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
211 <
209 >  
210 >  
211 >  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
212 >                                                             RealType variance)
213 >  {
214      
505 –    Vector<RealType, 6> Z;
506 –    Vector<RealType, 6> generalForce;
507 –    
508 –
509 –    Z[0] = randNumGen_.randNorm(0, variance);
510 –    Z[1] = randNumGen_.randNorm(0, variance);
511 –    Z[2] = randNumGen_.randNorm(0, variance);
512 –    Z[3] = randNumGen_.randNorm(0, variance);
513 –    Z[4] = randNumGen_.randNorm(0, variance);
514 –    Z[5] = randNumGen_.randNorm(0, variance);
515 –    
516 –    generalForce = hydroProps_[index]->getS()*Z;
517 –    
518 –    force[0] = generalForce[0];
519 –    force[1] = generalForce[1];
520 –    force[2] = generalForce[2];
521 –    torque[0] = generalForce[3];
522 –    torque[1] = generalForce[4];
523 –    torque[2] = generalForce[5];
524 –    
525 – }
526 –  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
527 –
215      // zero fill the random vector before starting:
216      std::vector<RealType> gaussRand;
217      gaussRand.resize(nTriangles);
218      std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
219 <
533 <
219 >    
220   #ifdef IS_MPI
221      if (worldRank == 0) {
222   #endif
223        for (int i = 0; i < nTriangles; i++) {
224 <        gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
224 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
225        }
226   #ifdef IS_MPI
227      }
228   #endif
229 <
229 >    
230      // push these out to the other processors
231 <
231 >    
232   #ifdef IS_MPI
233      if (worldRank == 0) {
234 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
234 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
235      } else {
236 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
236 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
237      }
238   #endif
239 <
554 <    for (int i = 0; i < nTriangles; i++) {
555 <      gaussRand[i] = gaussRand[i] * variance;
556 <    }
557 <
239 >    
240      return gaussRand;
241    }
560 –
561 –
562 –
563 –
564 –
242   }

Diff Legend

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