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 3461 by chuckv, Wed Oct 15 18:26:01 2008 UTC vs.
Revision 3516 by gezelter, Wed Jul 22 15:00:21 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();
59      
62 –    
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 (K)!\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 (atm)!\n");      
78        painCave.isFatal = 1;
79        simError();
80      } else {
81 <      targetPressure_ = simParams->getTargetPressure();
81 >      // Convert pressure from atm -> amu/(fs^2*Ang)
82 >      targetPressure_ = simParams->getTargetPressure() /
83 >        OOPSEConstant::pressureConvert;
84      }
84 –
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 (W m^-1 K^-1)!\n");
98 +      painCave.isFatal = 1;
99 +      painCave.severity = OOPSE_ERROR;
100 +      simError();
101 +    }else{
102 +      thermalConductivity_ = simParams->getThermalConductivity() *
103 +        OOPSEConstant::thermalConductivityConvert;
104 +    }
105  
106 +    if (!simParams->haveThermalLength()) {
107 +      sprintf(painCave.errMsg,
108 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
109 +              "   without a thermalLength (Angstroms)!\n");
110 +      painCave.isFatal = 1;
111 +      painCave.severity = OOPSE_ERROR;
112 +      simError();
113 +    }else{
114 +      thermalLength_ = simParams->getThermalLength();
115 +    }
116 +    
117 +    dt_ = simParams->getDt();
118 +  
119 +    variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;
120  
121 <    // Build the hydroProp map:
122 <    std::map<std::string, HydroProp*> hydroPropMap;
97 <
121 >    // Build a vector of integrable objects to determine if the are
122 >    // surface atoms
123      Molecule* mol;
124      StuntDouble* integrableObject;
125      SimInfo::MoleculeIterator i;
126 <    Molecule::IntegrableObjectIterator  j;              
127 <    bool needHydroPropFile = false;
128 <    
129 <    for (mol = info->beginMolecule(i); mol != NULL;
105 <         mol = info->nextMolecule(i)) {
126 >    Molecule::IntegrableObjectIterator  j;
127 >
128 >    for (mol = info_->beginMolecule(i); mol != NULL;
129 >         mol = info_->nextMolecule(i)) {          
130        for (integrableObject = mol->beginIntegrableObject(j);
131             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 SMIPDynamics\n"
125 –                 "\tis specified for rigidBodies which contain more than one atom\n"
126 –                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
127 –                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
128 –                 "\tset to 1.0 because the friction and random forces will be\n"
129 –                 "\tdynamically re-set assuming this is true.\n");
130 –        painCave.severity = OOPSE_ERROR;
131 –        painCave.isFatal = 1;
132 –        simError();  
133 –      }      
134 –
135 –      for (mol = info->beginMolecule(i); mol != NULL;
136 –           mol = info->nextMolecule(i)) {
137 –        for (integrableObject = mol->beginIntegrableObject(j);
138 –             integrableObject != NULL;
139 –             integrableObject = mol->nextIntegrableObject(j)) {
140 –
141 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
142 –          if (iter != hydroPropMap.end()) {
143 –            hydroProps_.push_back(iter->second);
144 –          } else {
145 –            sprintf( painCave.errMsg,
146 –                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
147 –            painCave.severity = OOPSE_ERROR;
148 –            painCave.isFatal = 1;
149 –            simError();  
150 –          }        
151 –        }
152 –      }
153 –    } else {
154 –      
155 –      std::map<std::string, HydroProp*> hydroPropMap;
156 –      for (mol = info->beginMolecule(i); mol != NULL;
157 –           mol = info->nextMolecule(i)) {
158 –        for (integrableObject = mol->beginIntegrableObject(j);
159 –             integrableObject != NULL;
160 –             integrableObject = mol->nextIntegrableObject(j)) {
161 –          Shape* currShape = NULL;
162 –
163 –          if (integrableObject->isAtom()){
164 –            Atom* atom = static_cast<Atom*>(integrableObject);
165 –            AtomType* atomType = atom->getAtomType();
166 –            if (atomType->isGayBerne()) {
167 –              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
168 –              GenericData* data = dAtomType->getPropertyByName("GayBerne");
169 –              if (data != NULL) {
170 –                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
171 –                
172 –                if (gayBerneData != NULL) {  
173 –                  GayBerneParam gayBerneParam = gayBerneData->getData();
174 –                  currShape = new Ellipsoid(V3Zero,
175 –                                            gayBerneParam.GB_l / 2.0,
176 –                                            gayBerneParam.GB_d / 2.0,
177 –                                            Mat3x3d::identity());
178 –                } else {
179 –                  sprintf( painCave.errMsg,
180 –                           "Can not cast GenericData to GayBerneParam\n");
181 –                  painCave.severity = OOPSE_ERROR;
182 –                  painCave.isFatal = 1;
183 –                  simError();  
184 –                }
185 –              } else {
186 –                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
187 –                painCave.severity = OOPSE_ERROR;
188 –                painCave.isFatal = 1;
189 –                simError();    
190 –              }
191 –            } else {
192 –              if (atomType->isLennardJones()){
193 –                GenericData* data = atomType->getPropertyByName("LennardJones");
194 –                if (data != NULL) {
195 –                  LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
196 –                  if (ljData != NULL) {
197 –                    LJParam ljParam = ljData->getData();
198 –                    currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
199 –                  } else {
200 –                    sprintf( painCave.errMsg,
201 –                             "Can not cast GenericData to LJParam\n");
202 –                    painCave.severity = OOPSE_ERROR;
203 –                    painCave.isFatal = 1;
204 –                    simError();          
205 –                  }      
206 –                }
207 –              } else {
208 –                int aNum = etab.GetAtomicNum((atom->getType()).c_str());
209 –                if (aNum != 0) {
210 –                  currShape = new Sphere(atom->getPos(), etab.GetVdwRad(aNum));
211 –                } else {
212 –                  sprintf( painCave.errMsg,
213 –                           "Could not find atom type in default element.txt\n");
214 –                  painCave.severity = OOPSE_ERROR;
215 –                  painCave.isFatal = 1;
216 –                  simError();          
217 –                }
218 –              }
219 –            }
220 –          }
221 –          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
222 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
223 –          if (iter != hydroPropMap.end())
224 –            hydroProps_.push_back(iter->second);
225 –          else {
226 –            currHydroProp->complete();
227 –            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
228 –            hydroProps_.push_back(currHydroProp);
229 –          }
230 –        }
231 –      }
232 –    }
233 –
234 –    /* Compute hull first time through to get the area of t=0*/
235 –
236 –    /* Build a vector of integrable objects to determine if the are surface atoms */
237 –    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
238 –      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
132             integrableObject = mol->nextIntegrableObject(j)) {  
133          localSites_.push_back(integrableObject);
134        }
135 <    }
243 <
244 <    surfaceMesh_->computeHull(localSites_);
245 <    Area0_ = surfaceMesh_->getArea();
246 <    //variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
247 <    
135 >    }  
136    }  
249 –
250 –  std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
251 –    std::map<std::string, HydroProp*> props;
252 –    std::ifstream ifs(filename.c_str());
253 –    if (ifs.is_open()) {
254 –      
255 –    }
256 –    
257 –    const unsigned int BufferSize = 65535;
258 –    char buffer[BufferSize];  
259 –    while (ifs.getline(buffer, BufferSize)) {
260 –      HydroProp* currProp = new HydroProp(buffer);
261 –      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
262 –    }
263 –
264 –    return props;
265 –  }
137    
138    void SMIPDForceManager::postCalculation(bool needStress){
139      SimInfo::MoleculeIterator i;
140      Molecule::IntegrableObjectIterator  j;
141      Molecule* mol;
142      StuntDouble* integrableObject;
272 –    RealType mass;
273 –    Vector3d pos;
274 –    Vector3d frc;
275 –    Mat3x3d A;
276 –    Mat3x3d Atrans;
277 –    Vector3d Tb;
278 –    Vector3d ji;
279 –    unsigned int index = 0;
280 –    int fdf;
281 –  
282 –    fdf = 0;
143    
144 <    /*Compute surface Mesh*/
144 >    // Compute surface Mesh
145      surfaceMesh_->computeHull(localSites_);
146  
147 <    /* Get area and number of surface stunt doubles and compute new variance */
148 <     RealType area = surfaceMesh_->getArea();
149 <     int nSurfaceSDs = surfaceMesh_->getNs();
147 >    // Get total area and number of surface stunt doubles
148 >    RealType area = surfaceMesh_->getArea();
149 >    int nSurfaceSDs = surfaceMesh_->getNs();
150 >    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
151 >    int nTriangles = sMesh.size();
152  
153 <     /* Compute variance for random forces */
154 <    
155 <     variance_ = sqrt(2.0*NumericConstant::PI)*((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs)
156 <       /OOPSEConstant::energyConvert;
157 <    
158 <    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
159 <    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
160 <    
299 <    /* Loop over the mesh faces and apply random force to each of the faces*/
300 <    
301 <    
302 <    std::vector<Triangle*>::iterator face;
153 >    // Generate all of the necessary random forces
154 >    std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
155 >
156 >    RealType instaTemp = thermo->getTemperature();
157 >
158 >    // Loop over the mesh faces and apply external pressure to each
159 >    // of the faces
160 >    std::vector<Triangle>::iterator face;
161      std::vector<StuntDouble*>::iterator vertex;
162 <    int thisNumber = 0;
162 >    int thisFacet = 0;
163      for (face = sMesh.begin(); face != sMesh.end(); ++face){
164      
165 <      Triangle* thisTriangle = *face;
166 <      std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
167 <      
168 <      /* Get Random Force */
311 <      Vector3d unitNormal = thisTriangle->getNormal();
165 >      Triangle thisTriangle = *face;
166 >      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
167 >      RealType thisArea = thisTriangle.getArea();
168 >      Vector3d unitNormal = thisTriangle.getNormal();
169        unitNormal.normalize();
170 <      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
171 <      Vector3d centroid = thisTriangle->getCentroid();
170 >      Vector3d centroid = thisTriangle.getCentroid();
171 >      Vector3d facetVel = thisTriangle.getFacetVelocity();
172 >      RealType thisMass = thisTriangle.getFacetMass();
173  
174 <      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
174 >      // gamma is the drag coefficient normal to the face of the triangle      
175 >      RealType gamma = thermalConductivity_ * thisMass * thisArea  
176 >        / (2.0 * thermalLength_ * OOPSEConstant::kB);
177 >      
178 >      gamma *= fabs(1.0 - targetTemp_/instaTemp);      
179 >      
180 >      RealType extPressure = - (targetPressure_ * thisArea) /
181 >        OOPSEConstant::energyConvert;
182  
183 <         // mass = integrableObject->getMass();
184 <        Vector3d vertexForce = randomForce/3.0;
320 <        (*vertex)->addFrc(vertexForce);
321 <        if (integrableObject->isDirectional()){
322 <          Vector3d vertexPos = (*vertex)->getPos();
323 <          Vector3d vertexCentroidVector = vertexPos - centroid;
324 <          (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
325 <        }
326 <          
327 <      }
328 <    }
183 >      RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
184 >      RealType dragForce = -gamma * dot(facetVel, unitNormal);
185  
186 <    /* Now loop over all surface particles and apply the drag*/
187 <
188 <    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
189 <    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
190 <      integrableObject = *vertex;
191 <      mass = integrableObject->getMass();
192 <      if (integrableObject->isDirectional()){
193 <        
194 <        // preliminaries for directional objects:
195 <        
196 <        A = integrableObject->getA();
197 <        Atrans = A.transpose();
342 <        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
343 <
344 <        
345 <        Mat3x3d I = integrableObject->getI();
346 <        Vector3d omegaBody;
347 <        
348 <        // What remains contains velocity explicitly, but the velocity required
349 <        // is at the full step: v(t + h), while we have initially the velocity
350 <        // at the half step: v(t + h/2).  We need to iterate to converge the
351 <        // friction force and friction torque vectors.
352 <        
353 <        // this is the velocity at the half-step:
354 <        
355 <        Vector3d vel =integrableObject->getVel();
356 <        Vector3d angMom = integrableObject->getJ();
357 <        
358 <        //estimate velocity at full-step using everything but friction forces:          
359 <        
360 <        frc = integrableObject->getFrc();
361 <        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
362 <        
363 <        Tb = integrableObject->lab2Body(integrableObject->getTrq());
364 <        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
365 <        
366 <        Vector3d omegaLab;
367 <        Vector3d vcdLab;
368 <        Vector3d vcdBody;
369 <        Vector3d frictionForceBody;
370 <        Vector3d frictionForceLab(0.0);
371 <        Vector3d oldFFL;  // used to test for convergence
372 <        Vector3d frictionTorqueBody(0.0);
373 <        Vector3d oldFTB;  // used to test for convergence
374 <        Vector3d frictionTorqueLab;
375 <        RealType fdot;
376 <        RealType tdot;
377 <
378 <        //iteration starts here:
379 <        
380 <        for (int k = 0; k < maxIterNum_; k++) {
381 <          
382 <          if (integrableObject->isLinear()) {
383 <            int linearAxis = integrableObject->linearAxis();
384 <            int l = (linearAxis +1 )%3;
385 <            int m = (linearAxis +2 )%3;
386 <            omegaBody[l] = angMomStep[l] /I(l, l);
387 <            omegaBody[m] = angMomStep[m] /I(m, m);
388 <            
389 <          } else {
390 <            omegaBody[0] = angMomStep[0] /I(0, 0);
391 <            omegaBody[1] = angMomStep[1] /I(1, 1);
392 <            omegaBody[2] = angMomStep[2] /I(2, 2);
186 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) *
187 >        unitNormal;
188 >      
189 >      // Apply triangle force to stuntdouble vertices
190 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
191 >        if ((*vertex) != NULL){
192 >          Vector3d vertexForce = langevinForce / 3.0;
193 >          (*vertex)->addFrc(vertexForce);          
194 >          if ((*vertex)->isDirectional()){          
195 >            Vector3d vertexPos = (*vertex)->getPos();
196 >            Vector3d vertexCentroidVector = vertexPos - centroid;
197 >            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
198            }
199 <          
395 <          omegaLab = Atrans * omegaBody;
396 <          
397 <          // apply friction force and torque at center of resistance
398 <          
399 <          vcdLab = velStep + cross(omegaLab, rcrLab);      
400 <          vcdBody = A * vcdLab;
401 <          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
402 <          oldFFL = frictionForceLab;
403 <          frictionForceLab = Atrans * frictionForceBody;
404 <          oldFTB = frictionTorqueBody;
405 <          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
406 <          frictionTorqueLab = Atrans * frictionTorqueBody;
407 <          
408 <          // re-estimate velocities at full-step using friction forces:
409 <              
410 <          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
411 <          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
412 <
413 <          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
414 <              
415 <          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
416 <          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
417 <          
418 <          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
419 <            break; // iteration ends here
420 <        }
421 <        
422 <        integrableObject->addFrc(frictionForceLab);
423 <        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
424 <
425 <            
426 <      } else {
427 <        //spherical atom
428 <
429 <        
430 <        // What remains contains velocity explicitly, but the velocity required
431 <        // is at the full step: v(t + h), while we have initially the velocity
432 <        // at the half step: v(t + h/2).  We need to iterate to converge the
433 <        // friction force vector.
434 <        
435 <        // this is the velocity at the half-step:
436 <        
437 <        Vector3d vel =integrableObject->getVel();
438 <        
439 <        //estimate velocity at full-step using everything but friction forces:          
440 <        
441 <        frc = integrableObject->getFrc();
442 <        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
443 <        
444 <        Vector3d frictionForce(0.0);
445 <        Vector3d oldFF;  // used to test for convergence
446 <        RealType fdot;
447 <        
448 <        //iteration starts here:
449 <        
450 <        for (int k = 0; k < maxIterNum_; k++) {
451 <          
452 <          oldFF = frictionForce;                            
453 <          frictionForce = -hydroProps_[index]->getXitt() * velStep;
454 <          
455 <          // re-estimate velocities at full-step using friction forces:
456 <          
457 <          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
458 <          
459 <          // check for convergence (if the vector has converged, fdot will be 1.0):
460 <          
461 <          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
462 <          
463 <          if (fabs(1.0 - fdot) <= forceTolerance_)
464 <            break; // iteration ends here
465 <        }
466 <        
467 <        integrableObject->addFrc(frictionForce);
468 <        
469 <        
199 >        }  
200        }
201 <      
472 <      
473 <    }
474 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
475 <    currSnapshot->setVolume(surfaceMesh_->getVolume());
201 >    }
202      
203 +    veloMunge->removeComDrift();
204 +    veloMunge->removeAngularDrift();
205 +    
206 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
207 +    currSnapshot->setVolume(surfaceMesh_->getVolume());    
208      ForceManager::postCalculation(needStress);  
209    }
210 <
211 <  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
212 <
210 >  
211 >  
212 >  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
213 >                                                             RealType variance)
214 >  {
215      
483 –    Vector<RealType, 6> Z;
484 –    Vector<RealType, 6> generalForce;
485 –    
486 –
487 –    Z[0] = randNumGen_.randNorm(0, variance);
488 –    Z[1] = randNumGen_.randNorm(0, variance);
489 –    Z[2] = randNumGen_.randNorm(0, variance);
490 –    Z[3] = randNumGen_.randNorm(0, variance);
491 –    Z[4] = randNumGen_.randNorm(0, variance);
492 –    Z[5] = randNumGen_.randNorm(0, variance);
493 –    
494 –    generalForce = hydroProps_[index]->getS()*Z;
495 –    
496 –    force[0] = generalForce[0];
497 –    force[1] = generalForce[1];
498 –    force[2] = generalForce[2];
499 –    torque[0] = generalForce[3];
500 –    torque[1] = generalForce[4];
501 –    torque[2] = generalForce[5];
502 –    
503 – }
504 –  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
505 –
216      // zero fill the random vector before starting:
217      std::vector<RealType> gaussRand;
218      gaussRand.resize(nTriangles);
219      std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
220 <
511 <
220 >    
221   #ifdef IS_MPI
222      if (worldRank == 0) {
223   #endif
224        for (int i = 0; i < nTriangles; i++) {
225 <        gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
225 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
226        }
227   #ifdef IS_MPI
228      }
229   #endif
230 <
230 >    
231      // push these out to the other processors
232 <
232 >    
233   #ifdef IS_MPI
234      if (worldRank == 0) {
235 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
235 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
236      } else {
237 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
237 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
238      }
239   #endif
240 <
532 <    for (int i = 0; i < nTriangles; i++) {
533 <      gaussRand[i] = gaussRand[i] * variance;
534 <    }
535 <
240 >    
241      return gaussRand;
242    }
538 –
243   }

Diff Legend

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