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 3473 by chuckv, Fri Nov 14 15:44:34 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();
61 –
58      
63 –    
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()/OOPSEConstant::pressureConvert;
80 >      // Convert pressure from atm -> amu/(fs^2*Ang)
81 >      targetPressure_ = simParams->getTargetPressure() /
82 >        OOPSEConstant::pressureConvert;
83      }
85 –
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->haveViscosity()) {
94 <      sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
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 <
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 <    //Compute initial hull
121 <    /*
107 <    surfaceMesh_->computeHull(localSites_);
108 <    Area0_ = surfaceMesh_->getArea();
109 <    int nTriangles = surfaceMesh_->getNMeshElements();
110 <    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
111 <    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
112 <      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
113 <    //RealType eta0 = gamma_0/
114 <    */
115 <
116 <    // Build the hydroProp map:
117 <    std::map<std::string, HydroProp*> hydroPropMap;
118 <
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;
126 <         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;
129 –           integrableObject = mol->nextIntegrableObject(j)) {
130 –        
131 –        if (integrableObject->isRigidBody()) {
132 –          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
133 –          if (rb->getNumAtoms() > 1) needHydroPropFile = true;
134 –        }
135 –        
136 –      }
137 –    }
138 –        
139 –
140 –    if (needHydroPropFile) {              
141 –      if (simParams->haveHydroPropFile()) {
142 –        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
143 –      } else {              
144 –        sprintf( painCave.errMsg,
145 –                 "HydroPropFile must be set to a file name if SMIPDynamics\n"
146 –                 "\tis specified for rigidBodies which contain more than one atom\n"
147 –                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
148 –                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
149 –                 "\tset to 1.0 because the friction and random forces will be\n"
150 –                 "\tdynamically re-set assuming this is true.\n");
151 –        painCave.severity = OOPSE_ERROR;
152 –        painCave.isFatal = 1;
153 –        simError();  
154 –      }      
155 –
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 –
162 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
163 –          if (iter != hydroPropMap.end()) {
164 –            hydroProps_.push_back(iter->second);
165 –          } else {
166 –            sprintf( painCave.errMsg,
167 –                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
168 –            painCave.severity = OOPSE_ERROR;
169 –            painCave.isFatal = 1;
170 –            simError();  
171 –          }        
172 –        }
173 –      }
174 –    } else {
175 –      
176 –      std::map<std::string, HydroProp*> hydroPropMap;
177 –      for (mol = info->beginMolecule(i); mol != NULL;
178 –           mol = info->nextMolecule(i)) {
179 –        for (integrableObject = mol->beginIntegrableObject(j);
180 –             integrableObject != NULL;
181 –             integrableObject = mol->nextIntegrableObject(j)) {
182 –          Shape* currShape = NULL;
183 –
184 –          if (integrableObject->isAtom()){
185 –            Atom* atom = static_cast<Atom*>(integrableObject);
186 –            AtomType* atomType = atom->getAtomType();
187 –            if (atomType->isGayBerne()) {
188 –              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
189 –              GenericData* data = dAtomType->getPropertyByName("GayBerne");
190 –              if (data != NULL) {
191 –                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
192 –                
193 –                if (gayBerneData != NULL) {  
194 –                  GayBerneParam gayBerneParam = gayBerneData->getData();
195 –                  currShape = new Ellipsoid(V3Zero,
196 –                                            gayBerneParam.GB_l / 2.0,
197 –                                            gayBerneParam.GB_d / 2.0,
198 –                                            Mat3x3d::identity());
199 –                } else {
200 –                  sprintf( painCave.errMsg,
201 –                           "Can not cast GenericData to GayBerneParam\n");
202 –                  painCave.severity = OOPSE_ERROR;
203 –                  painCave.isFatal = 1;
204 –                  simError();  
205 –                }
206 –              } else {
207 –                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
208 –                painCave.severity = OOPSE_ERROR;
209 –                painCave.isFatal = 1;
210 –                simError();    
211 –              }
212 –            } else {
213 –              if (atomType->isLennardJones()){
214 –                GenericData* data = atomType->getPropertyByName("LennardJones");
215 –                if (data != NULL) {
216 –                  LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
217 –                  if (ljData != NULL) {
218 –                    LJParam ljParam = ljData->getData();
219 –                    currShape = new Sphere(atom->getPos(), 2.0);
220 –                  } else {
221 –                    sprintf( painCave.errMsg,
222 –                             "Can not cast GenericData to LJParam\n");
223 –                    painCave.severity = OOPSE_ERROR;
224 –                    painCave.isFatal = 1;
225 –                    simError();          
226 –                  }      
227 –                }
228 –              } else {
229 –                int aNum = etab.GetAtomicNum((atom->getType()).c_str());
230 –                if (aNum != 0) {
231 –                  currShape = new Sphere(atom->getPos(), 2.0);
232 –                } else {
233 –                  sprintf( painCave.errMsg,
234 –                           "Could not find atom type in default element.txt\n");
235 –                  painCave.severity = OOPSE_ERROR;
236 –                  painCave.isFatal = 1;
237 –                  simError();          
238 –                }
239 –              }
240 –            }
241 –          }
242 –          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
243 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
244 –          if (iter != hydroPropMap.end())
245 –            hydroProps_.push_back(iter->second);
246 –          else {
247 –            currHydroProp->complete();
248 –            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
249 –            hydroProps_.push_back(currHydroProp);
250 –          }
251 –        }
252 –      }
253 –    }
254 –
255 –    /* Compute hull first time through to get the area of t=0*/
256 –
257 –    //Build a vector of integrable objects to determine if the are surface atoms
258 –    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
259 –      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
131             integrableObject = mol->nextIntegrableObject(j)) {  
132          localSites_.push_back(integrableObject);
133        }
134 <    }
264 <
265 <
134 >    }  
135    }  
267 –
268 –  std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
269 –    std::map<std::string, HydroProp*> props;
270 –    std::ifstream ifs(filename.c_str());
271 –    if (ifs.is_open()) {
272 –      
273 –    }
274 –    
275 –    const unsigned int BufferSize = 65535;
276 –    char buffer[BufferSize];  
277 –    while (ifs.getline(buffer, BufferSize)) {
278 –      HydroProp* currProp = new HydroProp(buffer);
279 –      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
280 –    }
281 –
282 –    return props;
283 –  }
136    
137    void SMIPDForceManager::postCalculation(bool needStress){
138      SimInfo::MoleculeIterator i;
139      Molecule::IntegrableObjectIterator  j;
140      Molecule* mol;
141      StuntDouble* integrableObject;
290 –    RealType mass;
291 –    Vector3d pos;
292 –    Vector3d frc;
293 –    Mat3x3d A;
294 –    Mat3x3d Atrans;
295 –    Vector3d Tb;
296 –    Vector3d ji;
297 –    unsigned int index = 0;
298 –    int fdf;
299 –  
300 –    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();
307 <     int nSurfaceSDs = surfaceMesh_->getNs();
308 <
309 <    
146 >    // Get total area and number of surface stunt doubles
147 >    RealType area = surfaceMesh_->getArea();
148      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
149      int nTriangles = sMesh.size();
150  
151 <
152 <
315 <     /* Compute variance for random forces */
316 <  
317 <
151 >    // Generate all of the necessary random forces
152 >    std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
153  
154 <
155 <    std::vector<RealType>  randNums = genTriangleForces(nTriangles, 1.0);
321 <
322 <    /* Loop over the mesh faces and apply random force to each of the faces*/
323 <    
324 <    
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 thisNumber = 0;
158 >    int thisFacet = 0;
159      for (face = sMesh.begin(); face != sMesh.end(); ++face){
329 –    
160        Triangle thisTriangle = *face;
161        std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
162        RealType thisArea = thisTriangle.getArea();
333 –      // RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*thisArea) /OOPSEConstant::energyConvert;
334 –      // gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());
335 –
336 –      /* Get Random Force */
163        Vector3d unitNormal = thisTriangle.getNormal();
164        unitNormal.normalize();
339 –      //Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
165        Vector3d centroid = thisTriangle.getCentroid();
166        Vector3d facetVel = thisTriangle.getFacetVelocity();
167 <      RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/3.14;
167 >      RealType thisMass = thisTriangle.getFacetMass();
168  
169 <      RealType f_normal = simParams->getViscosity()*hydroLength*1.439326479e4*OOPSEConstant::energyConvert;
170 <      RealType extPressure = -(targetPressure_ * thisArea);
171 <      RealType randomForce = randNums[thisNumber] * f_normal * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
347 <      RealType dragForce = -f_normal * dot(facetVel, unitNormal);      
348 <      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
349 <      
350 <      //      Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
351 <      
352 <      //std::cout << "randomForce " << randomForce << " dragForce " << dragForce <<  " hydro  " << hydroLength << std::endl;
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 +      RealType extPressure = - (targetPressure_ * thisArea) /
174 +        OOPSEConstant::energyConvert;
175 +      RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
176 +      RealType dragForce = -gamma * dot(facetVel, unitNormal);
177  
178 +      Vector3d langevinForce = (extPressure + randomForce + dragForce) *
179 +        unitNormal;
180 +      
181 +      // Apply triangle force to stuntdouble vertices
182        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
183          if ((*vertex) != NULL){
184 <          // mass = integrableObject->getMass();
185 <          Vector3d vertexForce = langevinForce/3.0;
359 <          (*vertex)->addFrc(vertexForce);
360 <
361 <          if ((*vertex)->isDirectional()){
362 <
363 <            Vector3d vertexPos = (*vertex)->getPos();
364 <            Vector3d vertexCentroidVector = vertexPos - centroid;
365 <            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
366 <          }
184 >          Vector3d vertexForce = langevinForce / 3.0;
185 >          (*vertex)->addFrc(vertexForce);          
186          }  
187        }
188      }
189 <
190 <    /* Now loop over all surface particles and apply the drag*/
191 <    /*
192 <    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
374 <    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
375 <      integrableObject = *vertex;
376 <      mass = integrableObject->getMass();
377 <      if (integrableObject->isDirectional()){
378 <        
379 <        // preliminaries for directional objects:
380 <        
381 <        A = integrableObject->getA();
382 <        Atrans = A.transpose();
383 <        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
384 <        //apply random force and torque at center of resistance
385 <        Mat3x3d I = integrableObject->getI();
386 <        Vector3d omegaBody;
387 <        
388 <        // What remains contains velocity explicitly, but the velocity required
389 <        // is at the full step: v(t + h), while we have initially the velocity
390 <        // at the half step: v(t + h/2).  We need to iterate to converge the
391 <        // friction force and friction torque vectors.
392 <        
393 <        // this is the velocity at the half-step:
394 <        
395 <        Vector3d vel =integrableObject->getVel();
396 <        Vector3d angMom = integrableObject->getJ();
397 <        
398 <        //estimate velocity at full-step using everything but friction forces:          
399 <        
400 <        frc = integrableObject->getFrc();
401 <        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
402 <        
403 <        Tb = integrableObject->lab2Body(integrableObject->getTrq());
404 <        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
405 <        
406 <        Vector3d omegaLab;
407 <        Vector3d vcdLab;
408 <        Vector3d vcdBody;
409 <        Vector3d frictionForceBody;
410 <        Vector3d frictionForceLab(0.0);
411 <        Vector3d oldFFL;  // used to test for convergence
412 <        Vector3d frictionTorqueBody(0.0);
413 <        Vector3d oldFTB;  // used to test for convergence
414 <        Vector3d frictionTorqueLab;
415 <        RealType fdot;
416 <        RealType tdot;
417 <
418 <        //iteration starts here:
419 <        
420 <        for (int k = 0; k < maxIterNum_; k++) {
421 <          
422 <          if (integrableObject->isLinear()) {
423 <            int linearAxis = integrableObject->linearAxis();
424 <            int l = (linearAxis +1 )%3;
425 <            int m = (linearAxis +2 )%3;
426 <            omegaBody[l] = angMomStep[l] /I(l, l);
427 <            omegaBody[m] = angMomStep[m] /I(m, m);
428 <            
429 <          } else {
430 <            omegaBody[0] = angMomStep[0] /I(0, 0);
431 <            omegaBody[1] = angMomStep[1] /I(1, 1);
432 <            omegaBody[2] = angMomStep[2] /I(2, 2);
433 <          }
434 <          
435 <          omegaLab = Atrans * omegaBody;
436 <          
437 <          // apply friction force and torque at center of resistance
438 <          
439 <          vcdLab = velStep + cross(omegaLab, rcrLab);      
440 <          vcdBody = A * vcdLab;
441 <          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
442 <          oldFFL = frictionForceLab;
443 <          frictionForceLab = Atrans * frictionForceBody;
444 <          oldFTB = frictionTorqueBody;
445 <          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
446 <          frictionTorqueLab = Atrans * frictionTorqueBody;
447 <          
448 <          // re-estimate velocities at full-step using friction forces:
449 <              
450 <          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
451 <          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
452 <
453 <          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
454 <              
455 <          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
456 <          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
457 <          
458 <          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
459 <            break; // iteration ends here
460 <        }
461 <        
462 <        integrableObject->addFrc(frictionForceLab);
463 <        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
464 <
465 <            
466 <      } else {
467 <        //spherical atom
468 <
469 <        // What remains contains velocity explicitly, but the velocity required
470 <        // is at the full step: v(t + h), while we have initially the velocity
471 <        // at the half step: v(t + h/2).  We need to iterate to converge the
472 <        // friction force vector.
473 <        
474 <        // this is the velocity at the half-step:
475 <        
476 <        Vector3d vel =integrableObject->getVel();
477 <        
478 <        //estimate velocity at full-step using everything but friction forces:          
479 <        
480 <        frc = integrableObject->getFrc();
481 <        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
482 <        
483 <        Vector3d frictionForce(0.0);
484 <        Vector3d oldFF;  // used to test for convergence
485 <        RealType fdot;
486 <        
487 <        //iteration starts here:
488 <        
489 <        for (int k = 0; k < maxIterNum_; k++) {
490 <          
491 <          oldFF = frictionForce;                            
492 <          frictionForce = -hydroProps_[index]->getXitt() * velStep;
493 <          //frictionForce = -gamma_t*velStep;
494 <          // re-estimate velocities at full-step using friction forces:
495 <          
496 <          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
497 <          
498 <          // check for convergence (if the vector has converged, fdot will be 1.0):
499 <          
500 <          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
501 <          
502 <          if (fabs(1.0 - fdot) <= forceTolerance_)
503 <            break; // iteration ends here
504 <        }
505 <        
506 <        integrableObject->addFrc(frictionForce);
507 <        
508 <        
509 <      }
510 <  
511 <      
512 <  }
513 <    */
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 <
197 >  
198 >  
199 >  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
200 >                                                             RealType variance)
201 >  {
202      
522 –    Vector<RealType, 6> Z;
523 –    Vector<RealType, 6> generalForce;
524 –    
525 –
526 –    Z[0] = randNumGen_.randNorm(0, variance);
527 –    Z[1] = randNumGen_.randNorm(0, variance);
528 –    Z[2] = randNumGen_.randNorm(0, variance);
529 –    Z[3] = randNumGen_.randNorm(0, variance);
530 –    Z[4] = randNumGen_.randNorm(0, variance);
531 –    Z[5] = randNumGen_.randNorm(0, variance);
532 –    
533 –    generalForce = hydroProps_[index]->getS()*Z;
534 –    
535 –    force[0] = generalForce[0];
536 –    force[1] = generalForce[1];
537 –    force[2] = generalForce[2];
538 –    torque[0] = generalForce[3];
539 –    torque[1] = generalForce[4];
540 –    torque[2] = generalForce[5];
541 –    
542 – }
543 –  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
544 –
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 <
550 <
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));    
556 <        gaussRand[i] = 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 <  
572 <    for (int i = 0; i < nTriangles; i++) {
573 <      gaussRand[i] = gaussRand[i] * variance;
574 <    }
575 <  
227 >    
228      return gaussRand;
229    }
578 –
579 –
580 –
581 –
582 –
230   }

Diff Legend

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