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 3464 by chuckv, Mon Oct 20 19:36:32 2008 UTC vs.
Revision 3481 by gezelter, Wed Nov 26 14:26:17 2008 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) {
51  
55 –  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
52      simParams = info->getSimParams();
57 –    veloMunge = new Velocitizer(info);
53      
54      // Create Hull, Convex Hull for now, other options later.
55 +    
56      surfaceMesh_ = new ConvexHull();
61 –
57      
63 –    
58      /* Check that the simulation has target pressure and target
59         temperature set*/
60 <
60 >    
61      if (!simParams->haveTargetTemp()) {
62 <      sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
62 >      sprintf(painCave.errMsg,
63 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
64 >              "   without a targetTemp!\n");      
65        painCave.isFatal = 1;
66        painCave.severity = OOPSE_ERROR;
67        simError();
68      } else {
69        targetTemp_ = simParams->getTargetTemp();
70      }
71 <
71 >    
72      if (!simParams->haveTargetPressure()) {
73 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
74 <              "   without a targetPressure!\n");
75 <      
73 >      sprintf(painCave.errMsg,
74 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
75 >              "   without a targetPressure!\n");      
76        painCave.isFatal = 1;
77        simError();
78      } else {
79 <      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
79 >      // Convert pressure from atm -> amu/(fs^2*Ang)
80 >      targetPressure_ = simParams->getTargetPressure() /
81 >        OOPSEConstant::pressureConvert;
82      }
85 –
83    
84      if (simParams->getUsePeriodicBoundaryConditions()) {
85 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
86 <              "   with periodic boundary conditions !\n");
87 <      
85 >      sprintf(painCave.errMsg,
86 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
87 >              "   with periodic boundary conditions!\n");    
88        painCave.isFatal = 1;
89        simError();
90      }
94 –
95 –
96 –
91      
92 +    if (!simParams->haveViscosity()) {
93 +      sprintf(painCave.errMsg,
94 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
95 +              "   without a viscosity!\n");
96 +      painCave.isFatal = 1;
97 +      painCave.severity = OOPSE_ERROR;
98 +      simError();
99 +    }else{
100 +      viscosity_ = simParams->getViscosity();
101 +    }
102 +    
103 +    dt_ = simParams->getDt();
104 +  
105 +    variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;
106  
99 –    //Compute initial hull
100 –    /*
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 –
107      Molecule* mol;
108      StuntDouble* integrableObject;
109      SimInfo::MoleculeIterator i;
110 <    Molecule::IntegrableObjectIterator  j;              
117 <    bool needHydroPropFile = false;
118 <    
119 <    for (mol = info->beginMolecule(i); mol != NULL;
120 <         mol = info->nextMolecule(i)) {
121 <      for (integrableObject = mol->beginIntegrableObject(j);
122 <           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 <        
110 >    Molecule::IntegrableObjectIterator  j;              
111  
112 <    if (needHydroPropFile) {              
113 <      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 <      }      
112 >    // Build a vector of integrable objects to determine if the are
113 >    // surface atoms
114  
115 <      for (mol = info->beginMolecule(i); mol != NULL;
116 <           mol = info->nextMolecule(i)) {
117 <        for (integrableObject = mol->beginIntegrableObject(j);
118 <             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;
115 >    for (mol = info_->beginMolecule(i); mol != NULL;
116 >         mol = info_->nextMolecule(i)) {          
117 >      for (integrableObject = mol->beginIntegrableObject(j);
118 >           integrableObject != NULL;
119             integrableObject = mol->nextIntegrableObject(j)) {  
120          localSites_.push_back(integrableObject);
121        }
122 <    }
258 <
259 <
122 >    }  
123    }  
124 <
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 <  }
278 <  
124 >  
125    void SMIPDForceManager::postCalculation(bool needStress){
126      SimInfo::MoleculeIterator i;
127      Molecule::IntegrableObjectIterator  j;
128      Molecule* mol;
129      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;
130    
131 <    fdf = 0;
295 <  
296 <    /*Compute surface Mesh*/
131 >    // Compute surface Mesh
132      surfaceMesh_->computeHull(localSites_);
298 –
299 –    /* Get area and number of surface stunt doubles and compute new variance */
300 –     RealType area = surfaceMesh_->getArea();
301 –     int nSurfaceSDs = surfaceMesh_->getNs();
302 –
133      
134 <    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
134 >    // Get total area and number of surface stunt doubles
135 >    RealType area = surfaceMesh_->getArea();
136 >    int nSurfaceSDs = surfaceMesh_->getNs();        
137 >    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
138      int nTriangles = sMesh.size();
139 +        
140 +    // Generate all of the necessary random forces
141 +    std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
142  
143 <
144 <
309 <     /* Compute variance for random forces */
310 <  
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 <    
322 <    std::vector<Triangle*>::iterator face;
143 >    // Loop over the mesh faces and apply random force to each of the faces
144 >    std::vector<Triangle>::iterator face;
145      std::vector<StuntDouble*>::iterator vertex;
146      int thisNumber = 0;
147      for (face = sMesh.begin(); face != sMesh.end(); ++face){
148      
149 <      Triangle* thisTriangle = *face;
150 <      std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
151 <      
152 <      /* Get Random Force */
331 <      Vector3d unitNormal = thisTriangle->getNormal();
149 >      Triangle thisTriangle = *face;
150 >      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
151 >      RealType thisArea = thisTriangle.getArea();
152 >      Vector3d unitNormal = thisTriangle.getNormal();
153        unitNormal.normalize();
333 –      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
334 –      Vector3d centroid = thisTriangle->getCentroid();
154  
155 <      Vector3d langevinForce = randomForce - gamma_t_*thisTriangle->getFacetVelocity();
155 >      Vector3d centroid = thisTriangle.getCentroid();
156 >      Vector3d facetVel = thisTriangle.getFacetVelocity();
157 >      RealType hydroLength = thisTriangle.getIncircleRadius() * 2.0 /
158 >        NumericConstant::PI;
159        
160 <      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
161 <        if ((*vertex) != NULL){
162 <          // mass = integrableObject->getMass();
163 <          Vector3d vertexForce = langevinForce/3.0;
164 <          (*vertex)->addFrc(vertexForce);
165 <          if (integrableObject->isDirectional()){
160 >      // gamma is the drag coefficient normal to the face of the triangle      
161 >      RealType gamma = viscosity_ * hydroLength *
162 >        OOPSEConstant::viscoConvert;
163 >
164 >      RealType extPressure = - (targetPressure_ * thisArea) /
165 >        OOPSEConstant::energyConvert;
166 >
167 >      RealType randomForce = randNums[thisNumber++] * sqrt(gamma);
168 >      RealType dragForce = -gamma * dot(facetVel, unitNormal);
169 >
170 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
171 >      
172 >      // Apply triangle force to stuntdouble vertices
173 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) {
174 >        if ((*vertex) != NULL) {
175 >          Vector3d vertexForce = langevinForce / 3.0;
176 >          (*vertex)->addFrc(vertexForce);          
177 >          if ((*vertex)->isDirectional()) {
178              Vector3d vertexPos = (*vertex)->getPos();
179              Vector3d vertexCentroidVector = vertexPos - centroid;
180              (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
# Line 348 | Line 182 | namespace oopse {
182          }  
183        }
184      }
351 –
352 –    /* Now loop over all surface particles and apply the drag*/
353 –    /*
354 –    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
355 –    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
356 –      integrableObject = *vertex;
357 –      mass = integrableObject->getMass();
358 –      if (integrableObject->isDirectional()){
359 –        
360 –        // preliminaries for directional objects:
361 –        
362 –        A = integrableObject->getA();
363 –        Atrans = A.transpose();
364 –        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
365 –        //apply random force and torque at center of resistance
366 –        Mat3x3d I = integrableObject->getI();
367 –        Vector3d omegaBody;
368 –        
369 –        // What remains contains velocity explicitly, but the velocity required
370 –        // is at the full step: v(t + h), while we have initially the velocity
371 –        // at the half step: v(t + h/2).  We need to iterate to converge the
372 –        // friction force and friction torque vectors.
373 –        
374 –        // this is the velocity at the half-step:
375 –        
376 –        Vector3d vel =integrableObject->getVel();
377 –        Vector3d angMom = integrableObject->getJ();
378 –        
379 –        //estimate velocity at full-step using everything but friction forces:          
380 –        
381 –        frc = integrableObject->getFrc();
382 –        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
383 –        
384 –        Tb = integrableObject->lab2Body(integrableObject->getTrq());
385 –        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
386 –        
387 –        Vector3d omegaLab;
388 –        Vector3d vcdLab;
389 –        Vector3d vcdBody;
390 –        Vector3d frictionForceBody;
391 –        Vector3d frictionForceLab(0.0);
392 –        Vector3d oldFFL;  // used to test for convergence
393 –        Vector3d frictionTorqueBody(0.0);
394 –        Vector3d oldFTB;  // used to test for convergence
395 –        Vector3d frictionTorqueLab;
396 –        RealType fdot;
397 –        RealType tdot;
398 –
399 –        //iteration starts here:
400 –        
401 –        for (int k = 0; k < maxIterNum_; k++) {
402 –          
403 –          if (integrableObject->isLinear()) {
404 –            int linearAxis = integrableObject->linearAxis();
405 –            int l = (linearAxis +1 )%3;
406 –            int m = (linearAxis +2 )%3;
407 –            omegaBody[l] = angMomStep[l] /I(l, l);
408 –            omegaBody[m] = angMomStep[m] /I(m, m);
409 –            
410 –          } else {
411 –            omegaBody[0] = angMomStep[0] /I(0, 0);
412 –            omegaBody[1] = angMomStep[1] /I(1, 1);
413 –            omegaBody[2] = angMomStep[2] /I(2, 2);
414 –          }
415 –          
416 –          omegaLab = Atrans * omegaBody;
417 –          
418 –          // apply friction force and torque at center of resistance
419 –          
420 –          vcdLab = velStep + cross(omegaLab, rcrLab);      
421 –          vcdBody = A * vcdLab;
422 –          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
423 –          oldFFL = frictionForceLab;
424 –          frictionForceLab = Atrans * frictionForceBody;
425 –          oldFTB = frictionTorqueBody;
426 –          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
427 –          frictionTorqueLab = Atrans * frictionTorqueBody;
428 –          
429 –          // re-estimate velocities at full-step using friction forces:
430 –              
431 –          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
432 –          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
433 –
434 –          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
435 –              
436 –          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
437 –          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
438 –          
439 –          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
440 –            break; // iteration ends here
441 –        }
442 –        
443 –        integrableObject->addFrc(frictionForceLab);
444 –        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
445 –
446 –            
447 –      } else {
448 –        //spherical atom
449 –
450 –        // What remains contains velocity explicitly, but the velocity required
451 –        // is at the full step: v(t + h), while we have initially the velocity
452 –        // at the half step: v(t + h/2).  We need to iterate to converge the
453 –        // friction force vector.
454 –        
455 –        // this is the velocity at the half-step:
456 –        
457 –        Vector3d vel =integrableObject->getVel();
458 –        
459 –        //estimate velocity at full-step using everything but friction forces:          
460 –        
461 –        frc = integrableObject->getFrc();
462 –        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
463 –        
464 –        Vector3d frictionForce(0.0);
465 –        Vector3d oldFF;  // used to test for convergence
466 –        RealType fdot;
467 –        
468 –        //iteration starts here:
469 –        
470 –        for (int k = 0; k < maxIterNum_; k++) {
471 –          
472 –          oldFF = frictionForce;                            
473 –          frictionForce = -hydroProps_[index]->getXitt() * velStep;
474 –          //frictionForce = -gamma_t*velStep;
475 –          // re-estimate velocities at full-step using friction forces:
476 –          
477 –          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
478 –          
479 –          // check for convergence (if the vector has converged, fdot will be 1.0):
480 –          
481 –          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
482 –          
483 –          if (fabs(1.0 - fdot) <= forceTolerance_)
484 –            break; // iteration ends here
485 –        }
486 –        
487 –        integrableObject->addFrc(frictionForce);
488 –        
489 –        
490 –      }
491 –  
492 –      
493 –  }
494 –    */
495 –    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
496 –    currSnapshot->setVolume(surfaceMesh_->getVolume());
185      
186 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
187 +    currSnapshot->setVolume(surfaceMesh_->getVolume());    
188      ForceManager::postCalculation(needStress);  
189    }
190  
191 <  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
191 >  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
192 >                                                             RealType variance) {
193  
503 –    
504 –    Vector<RealType, 6> Z;
505 –    Vector<RealType, 6> generalForce;
506 –    
507 –
508 –    Z[0] = randNumGen_.randNorm(0, variance);
509 –    Z[1] = randNumGen_.randNorm(0, variance);
510 –    Z[2] = randNumGen_.randNorm(0, variance);
511 –    Z[3] = randNumGen_.randNorm(0, variance);
512 –    Z[4] = randNumGen_.randNorm(0, variance);
513 –    Z[5] = randNumGen_.randNorm(0, variance);
514 –    
515 –    generalForce = hydroProps_[index]->getS()*Z;
516 –    
517 –    force[0] = generalForce[0];
518 –    force[1] = generalForce[1];
519 –    force[2] = generalForce[2];
520 –    torque[0] = generalForce[3];
521 –    torque[1] = generalForce[4];
522 –    torque[2] = generalForce[5];
523 –    
524 – }
525 –  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
526 –
194      // zero fill the random vector before starting:
195      std::vector<RealType> gaussRand;
196      gaussRand.resize(nTriangles);
197      std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
198 <
532 <
198 >  
199   #ifdef IS_MPI
200      if (worldRank == 0) {
201   #endif
202        for (int i = 0; i < nTriangles; i++) {
203 <        gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
203 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
204        }
205   #ifdef IS_MPI
206      }
# Line 544 | Line 210 | namespace oopse {
210  
211   #ifdef IS_MPI
212      if (worldRank == 0) {
213 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
213 >      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
214      } else {
215 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
215 >      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
216      }
217   #endif
218 <
553 <    for (int i = 0; i < nTriangles; i++) {
554 <      gaussRand[i] = gaussRand[i] * variance;
555 <    }
556 <
218 >  
219      return gaussRand;
220    }
559 –
560 –
561 –
562 –
563 –
221   }

Diff Legend

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