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 3450 by chuckv, Sun Sep 14 01:32:26 2008 UTC vs.
Revision 3550 by gezelter, Wed Oct 21 02:49:43 2009 UTC

# Line 1 | Line 1
1   /*
2 < * Copyright (c) 2008 The University of Notre Dame. All Rights Reserved.
2 > * Copyright (c) 2008, 2009 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
5   * non-exclusive, royalty free, license to use, modify and
# 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      veloMunge = new Velocitizer(info);
54      
55      // Create Hull, Convex Hull for now, other options later.
56 +    
57      surfaceMesh_ = new ConvexHull();
58      
62 –    
59      /* Check that the simulation has target pressure and target
60 <       temperature set*/
61 <
60 >       temperature set */
61 >    
62      if (!simParams->haveTargetTemp()) {
63 <      sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
63 >      sprintf(painCave.errMsg,
64 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
65 >              "\twithout a targetTemp (K)!\n");      
66        painCave.isFatal = 1;
67        painCave.severity = OOPSE_ERROR;
68        simError();
69      } else {
70        targetTemp_ = simParams->getTargetTemp();
71      }
72 <
72 >    
73      if (!simParams->haveTargetPressure()) {
74 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
75 <              "   without a targetPressure!\n");
76 <      
74 >      sprintf(painCave.errMsg,
75 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
76 >              "\twithout a targetPressure (atm)!\n");      
77        painCave.isFatal = 1;
78        simError();
79      } else {
80 <      targetPressure_ = simParams->getTargetPressure();
80 >      // Convert pressure from atm -> amu/(fs^2*Ang)
81 >      targetPressure_ = simParams->getTargetPressure() /
82 >        OOPSEConstant::pressureConvert;
83      }
84 –
84    
85      if (simParams->getUsePeriodicBoundaryConditions()) {
86 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
87 <              "   with periodic boundary conditions !\n");
88 <      
86 >      sprintf(painCave.errMsg,
87 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
88 >              "\twith periodic boundary conditions!\n");    
89        painCave.isFatal = 1;
90        simError();
91      }
92 +    
93 +    if (!simParams->haveThermalConductivity()) {
94 +      sprintf(painCave.errMsg,
95 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
96 +              "\twithout a thermalConductivity (W m^-1 K^-1)!\n");
97 +      painCave.isFatal = 1;
98 +      painCave.severity = OOPSE_ERROR;
99 +      simError();
100 +    }else{
101 +      thermalConductivity_ = simParams->getThermalConductivity() *
102 +        OOPSEConstant::thermalConductivityConvert;
103 +    }
104  
105 +    if (!simParams->haveThermalLength()) {
106 +      sprintf(painCave.errMsg,
107 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
108 +              "\twithout a thermalLength (Angstroms)!\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 <    // Build the hydroProp map:
121 <    std::map<std::string, HydroProp*> hydroPropMap;
97 <
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;
105 <         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;
108 –           integrableObject = mol->nextIntegrableObject(j)) {
109 –        
110 –        if (integrableObject->isRigidBody()) {
111 –          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
112 –          if (rb->getNumAtoms() > 1) needHydroPropFile = true;
113 –        }
114 –        
115 –      }
116 –    }
117 –        
118 –
119 –    if (needHydroPropFile) {              
120 –      if (simParams->haveHydroPropFile()) {
121 –        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
122 –      } else {              
123 –        sprintf( painCave.errMsg,
124 –                 "HydroPropFile must be set to a file name if Langevin Dynamics\n"
125 –                 "\tis specified for rigidBodies which contain more than one atom\n"
126 –                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n");
127 –        painCave.severity = OOPSE_ERROR;
128 –        painCave.isFatal = 1;
129 –        simError();  
130 –      }      
131 –
132 –      for (mol = info->beginMolecule(i); mol != NULL;
133 –           mol = info->nextMolecule(i)) {
134 –        for (integrableObject = mol->beginIntegrableObject(j);
135 –             integrableObject != NULL;
136 –             integrableObject = mol->nextIntegrableObject(j)) {
137 –
138 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
139 –          if (iter != hydroPropMap.end()) {
140 –            hydroProps_.push_back(iter->second);
141 –          } else {
142 –            sprintf( painCave.errMsg,
143 –                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
144 –            painCave.severity = OOPSE_ERROR;
145 –            painCave.isFatal = 1;
146 –            simError();  
147 –          }        
148 –        }
149 –      }
150 –    } else {
151 –      
152 –      std::map<std::string, HydroProp*> hydroPropMap;
153 –      for (mol = info->beginMolecule(i); mol != NULL;
154 –           mol = info->nextMolecule(i)) {
155 –        for (integrableObject = mol->beginIntegrableObject(j);
156 –             integrableObject != NULL;
157 –             integrableObject = mol->nextIntegrableObject(j)) {
158 –          Shape* currShape = NULL;
159 –
160 –          if (integrableObject->isAtom()){
161 –            Atom* atom = static_cast<Atom*>(integrableObject);
162 –            AtomType* atomType = atom->getAtomType();
163 –            if (atomType->isGayBerne()) {
164 –              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
165 –              GenericData* data = dAtomType->getPropertyByName("GayBerne");
166 –              if (data != NULL) {
167 –                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
168 –                
169 –                if (gayBerneData != NULL) {  
170 –                  GayBerneParam gayBerneParam = gayBerneData->getData();
171 –                  currShape = new Ellipsoid(V3Zero,
172 –                                            gayBerneParam.GB_l / 2.0,
173 –                                            gayBerneParam.GB_d / 2.0,
174 –                                            Mat3x3d::identity());
175 –                } else {
176 –                  sprintf( painCave.errMsg,
177 –                           "Can not cast GenericData to GayBerneParam\n");
178 –                  painCave.severity = OOPSE_ERROR;
179 –                  painCave.isFatal = 1;
180 –                  simError();  
181 –                }
182 –              } else {
183 –                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
184 –                painCave.severity = OOPSE_ERROR;
185 –                painCave.isFatal = 1;
186 –                simError();    
187 –              }
188 –            } else {
189 –              if (atomType->isLennardJones()){
190 –                GenericData* data = atomType->getPropertyByName("LennardJones");
191 –                if (data != NULL) {
192 –                  LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
193 –                  if (ljData != NULL) {
194 –                    LJParam ljParam = ljData->getData();
195 –                    currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
196 –                  } else {
197 –                    sprintf( painCave.errMsg,
198 –                             "Can not cast GenericData to LJParam\n");
199 –                    painCave.severity = OOPSE_ERROR;
200 –                    painCave.isFatal = 1;
201 –                    simError();          
202 –                  }      
203 –                }
204 –              } else {
205 –                int aNum = etab.GetAtomicNum((atom->getType()).c_str());
206 –                if (aNum != 0) {
207 –                  currShape = new Sphere(atom->getPos(), etab.GetVdwRad(aNum));
208 –                } else {
209 –                  sprintf( painCave.errMsg,
210 –                           "Could not find atom type in default element.txt\n");
211 –                  painCave.severity = OOPSE_ERROR;
212 –                  painCave.isFatal = 1;
213 –                  simError();          
214 –                }
215 –              }
216 –            }
217 –          }
218 –          HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
219 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
220 –          if (iter != hydroPropMap.end())
221 –            hydroProps_.push_back(iter->second);
222 –          else {
223 –            currHydroProp->complete();
224 –            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
225 –            hydroProps_.push_back(currHydroProp);
226 –          }
227 –        }
228 –      }
229 –    }
230 –
231 –    /* Compute hull first time through to get the area of t=0*/
232 –
233 –    /* Build a vector of integrable objects to determine if the are surface atoms */
234 –    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
235 –      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
131             integrableObject = mol->nextIntegrableObject(j)) {  
132          localSites_.push_back(integrableObject);
133        }
134 <    }
240 <
241 <    surfaceMesh_->computeHull(localSites_);
242 <    Area0_ = surfaceMesh_->getArea();
243 <
244 <
134 >    }  
135    }  
246 –
247 –  std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
248 –    std::map<std::string, HydroProp*> props;
249 –    std::ifstream ifs(filename.c_str());
250 –    if (ifs.is_open()) {
251 –      
252 –    }
253 –    
254 –    const unsigned int BufferSize = 65535;
255 –    char buffer[BufferSize];  
256 –    while (ifs.getline(buffer, BufferSize)) {
257 –      HydroProp* currProp = new HydroProp(buffer);
258 –      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
259 –    }
260 –
261 –    return props;
262 –  }
136    
137    void SMIPDForceManager::postCalculation(bool needStress){
138      SimInfo::MoleculeIterator i;
139      Molecule::IntegrableObjectIterator  j;
140      Molecule* mol;
141      StuntDouble* integrableObject;
142 <    RealType mass;
143 <    Vector3d pos;
271 <    Vector3d frc;
272 <    Mat3x3d A;
273 <    Mat3x3d Atrans;
274 <    Vector3d Tb;
275 <    Vector3d ji;
276 <    unsigned int index = 0;
277 <    int fdf;
278 <  
279 <    fdf = 0;
280 <
281 <  
282 <    /*Compute surface Mesh*/
142 >  
143 >    // Compute surface Mesh
144      surfaceMesh_->computeHull(localSites_);
145  
146 <    /* Get area and number of surface stunt doubles and compute new variance */
146 >    // Get total area and number of surface stunt doubles
147      RealType area = surfaceMesh_->getArea();
148 <    RealType nSurfaceSDs = surfaceMesh_->getNs();
149 <    
289 <    std::cerr << "Surface Area is: " << area << " nSurfaceSDs is: " << nSurfaceSDs << std::endl;
148 >    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
149 >    int nTriangles = sMesh.size();
150  
151 <    /* Compute variance for random forces */
152 <    
293 <    variance_ = sqrt(2.0*NumericConstant::PI)*(targetPressure_*area/nSurfaceSDs);
151 >    // Generate all of the necessary random forces
152 >    std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
153  
154 <    /* Loop over the mesh faces and apply random force to each of the faces*/
155 <
156 <    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
298 <    std::vector<Triangle*>::iterator face;
154 >    // Loop over the mesh faces and apply external pressure to each
155 >    // of the faces
156 >    std::vector<Triangle>::iterator face;
157      std::vector<StuntDouble*>::iterator vertex;
158 +    int thisFacet = 0;
159      for (face = sMesh.begin(); face != sMesh.end(); ++face){
160 <    
161 <      Triangle* thisTriangle = *face;
162 <      std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
160 >      Triangle thisTriangle = *face;
161 >      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
162 >      RealType thisArea = thisTriangle.getArea();
163 >      Vector3d unitNormal = thisTriangle.getNormal();
164 >      unitNormal.normalize();
165 >      Vector3d centroid = thisTriangle.getCentroid();
166 >      Vector3d facetVel = thisTriangle.getFacetVelocity();
167 >      RealType thisMass = thisTriangle.getFacetMass();
168  
169 <      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
170 <         Vector3d randomForce;
171 <         Vector3d randomTorque;
308 <         genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
309 <         mass = integrableObject->getMass();
169 >      // gamma is the drag coefficient normal to the face of the triangle      
170 >      RealType gamma = thermalConductivity_ * thisMass * thisArea
171 >        / (2.0 * thermalLength_ * OOPSEConstant::kB);
172  
173 <         integrableObject->addFrc(randomForce);          
174 <      }
173 >      RealType extPressure = - (targetPressure_ * thisArea) /
174 >        OOPSEConstant::energyConvert;
175 >      RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
176 >      RealType dragForce = -gamma * dot(facetVel, unitNormal);
177  
178 <
179 <    }
316 <
317 <
318 <
319 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
320 <
178 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) *
179 >        unitNormal;
180        
181 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
182 <           integrableObject = mol->nextIntegrableObject(j)) {
183 <          
184 <          mass = integrableObject->getMass();
185 <          if (integrableObject->isDirectional()){
186 <
328 <            // preliminaries for directional objects:
329 <
330 <            A = integrableObject->getA();
331 <            Atrans = A.transpose();
332 <            Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
333 <
334 <            //apply random force and torque at center of resistance
335 <
336 <            Vector3d randomForceBody;
337 <            Vector3d randomTorqueBody;
338 <            genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
339 <            Vector3d randomForceLab = Atrans * randomForceBody;
340 <            Vector3d randomTorqueLab = Atrans * randomTorqueBody;
341 <            integrableObject->addFrc(randomForceLab);            
342 <            integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
343 <
344 <            Mat3x3d I = integrableObject->getI();
345 <            Vector3d omegaBody;
346 <
347 <            // What remains contains velocity explicitly, but the velocity required
348 <            // is at the full step: v(t + h), while we have initially the velocity
349 <            // at the half step: v(t + h/2).  We need to iterate to converge the
350 <            // friction force and friction torque vectors.
351 <
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 <              }
393 <              
394 <              omegaLab = Atrans * omegaBody;
395 <              
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);
411 <
412 <              // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
413 <              
414 <              fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
415 <              tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
416 <              
417 <              if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
418 <                break; // iteration ends here
419 <            }
420 <
421 <            integrableObject->addFrc(frictionForceLab);
422 <            integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
423 <
424 <            
425 <          } else {
426 <            //spherical atom
427 <
428 <            Vector3d randomForce;
429 <            Vector3d randomTorque;
430 <            genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
431 <            integrableObject->addFrc(randomForce);            
432 <
433 <            // What remains contains velocity explicitly, but the velocity required
434 <            // is at the full step: v(t + h), while we have initially the velocity
435 <            // at the half step: v(t + h/2).  We need to iterate to converge the
436 <            // friction force vector.
437 <
438 <            // this is the velocity at the half-step:
439 <            
440 <            Vector3d vel =integrableObject->getVel();
441 <
442 <            //estimate velocity at full-step using everything but friction forces:          
443 <
444 <            frc = integrableObject->getFrc();
445 <            Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
446 <
447 <            Vector3d frictionForce(0.0);
448 <            Vector3d oldFF;  // used to test for convergence
449 <            RealType fdot;
450 <
451 <            //iteration starts here:
452 <
453 <            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 <
472 <          
473 <          }
474 <          
475 <        ++index;
476 <    
181 >      // Apply triangle force to stuntdouble vertices
182 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
183 >        if ((*vertex) != NULL){
184 >          Vector3d vertexForce = langevinForce / 3.0;
185 >          (*vertex)->addFrc(vertexForce);          
186 >        }  
187        }
188 <    }    
189 <
190 <    info_->setFdf(fdf);
191 <    // veloMunge->removeComDrift();
192 <    // Remove angular drift if we are not using periodic boundary conditions.
193 <    //if(!simParams->getUsePeriodicBoundaryConditions())
194 <    //  veloMunge->removeAngularDrift();
485 <
188 >    }
189 >    
190 >    veloMunge->removeComDrift();
191 >    veloMunge->removeAngularDrift();
192 >    
193 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
194 >    currSnapshot->setVolume(surfaceMesh_->getVolume());    
195      ForceManager::postCalculation(needStress);  
196    }
197 <
198 < void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
199 <
200 <
201 <    Vector<RealType, 6> Z;
493 <    Vector<RealType, 6> generalForce;
494 <        
495 <    Z[0] = randNumGen_.randNorm(0, variance);
496 <    Z[1] = randNumGen_.randNorm(0, variance);
497 <    Z[2] = randNumGen_.randNorm(0, variance);
498 <    Z[3] = randNumGen_.randNorm(0, variance);
499 <    Z[4] = randNumGen_.randNorm(0, variance);
500 <    Z[5] = randNumGen_.randNorm(0, variance);
501 <    
502 <    generalForce = hydroProps_[index]->getS()*Z;
197 >  
198 >  
199 >  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
200 >                                                             RealType variance)
201 >  {
202      
504 –    force[0] = generalForce[0];
505 –    force[1] = generalForce[1];
506 –    force[2] = generalForce[2];
507 –    torque[0] = generalForce[3];
508 –    torque[1] = generalForce[4];
509 –    torque[2] = generalForce[5];
510 –    
511 – }
512 –  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
513 –
203      // zero fill the random vector before starting:
204      std::vector<RealType> gaussRand;
205      gaussRand.resize(nTriangles);
206      std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
207 <
519 <
207 >    
208   #ifdef IS_MPI
209      if (worldRank == 0) {
210   #endif
211        for (int i = 0; i < nTriangles; i++) {
212 <        gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
212 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
213        }
214   #ifdef IS_MPI
215      }
216   #endif
217 <
217 >    
218      // push these out to the other processors
219 <
219 >    
220   #ifdef IS_MPI
221      if (worldRank == 0) {
222 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
222 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
223      } else {
224 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
224 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
225      }
226   #endif
227 <
540 <    for (int i = 0; i < nTriangles; i++) {
541 <      gaussRand[i] = gaussRand[i] * variance;
542 <    }
543 <
227 >    
228      return gaussRand;
229    }
546 –
230   }

Diff Legend

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