src/integrators/SMIPDForceManager.cpp

/*
 * Copyright (c) 2008 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
#include <fstream> 
#include <iostream>
#include "integrators/SMIPDForceManager.hpp"
#include "math/CholeskyDecomposition.hpp"
#include "utils/OOPSEConstant.hpp"
#include "hydrodynamics/Sphere.hpp"
#include "hydrodynamics/Ellipsoid.hpp"
#include "utils/ElementsTable.hpp"
#include "math/ConvexHull.hpp"
#include "math/Triangle.hpp"


namespace oopse {
  
  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
    simParams = info->getSimParams();
    veloMunge = new Velocitizer(info);
    
    // Create Hull, Convex Hull for now, other options later.
    surfaceMesh_ = new ConvexHull();

    
    /* Check that the simulation has target pressure and target
       temperature set*/
    
    if (!simParams->haveTargetTemp()) {
      sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
      painCave.isFatal = 1;
      painCave.severity = OOPSE_ERROR;
      simError();
    } else {
      targetTemp_ = simParams->getTargetTemp();
    }
    
    if (!simParams->haveTargetPressure()) {
      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "   without a targetPressure!\n");
      
      painCave.isFatal = 1;
      simError();
    } else {
      /* Convert pressure from atm -> amu/(fs^2*Ang)*/
      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
    }

   
    if (simParams->getUsePeriodicBoundaryConditions()) {
      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
              "   with periodic boundary conditions !\n");
      
      painCave.isFatal = 1;
      simError();
    } 

    if (!simParams->haveViscosity()) {
      sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
      painCave.isFatal = 1;
      painCave.severity = OOPSE_ERROR;
      simError();
    }else{
      viscosity_ = simParams->getViscosity();
    }

    dt_ = simParams->getDt();

    
    //Compute initial hull
    /*
    surfaceMesh_->computeHull(localSites_);
    Area0_ = surfaceMesh_->getArea();
    int nTriangles = surfaceMesh_->getNMeshElements();
    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
    //RealType eta0 = gamma_0/
    */


    Molecule* mol;
    StuntDouble* integrableObject;
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;              

    
    /* Compute hull first time through to get the area of t=0*/

    //Build a vector of integrable objects to determine if the are surface atoms 
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {   
        localSites_.push_back(integrableObject);
      }
    }


  }  

  std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
    std::map<std::string, HydroProp*> props;
    std::ifstream ifs(filename.c_str());
    if (ifs.is_open()) {
      
    }
    
    const unsigned int BufferSize = 65535;
    char buffer[BufferSize];   
    while (ifs.getline(buffer, BufferSize)) {
      HydroProp* currProp = new HydroProp(buffer);
      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
    }

    return props;
  }
   
  void SMIPDForceManager::postCalculation(bool needStress){
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    RealType mass;
    Vector3d pos;
    Vector3d frc;
    Mat3x3d A;
    Mat3x3d Atrans;
    Vector3d Tb;
    Vector3d ji;
    unsigned int index = 0;
    int fdf;
   
    fdf = 0;
  
    /*Compute surface Mesh*/
    surfaceMesh_->computeHull(localSites_);

    /* Get area and number of surface stunt doubles and compute new variance */
     RealType area = surfaceMesh_->getArea();
     int nSurfaceSDs = surfaceMesh_->getNs();

    
    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
    int nTriangles = sMesh.size();


     /* Compute variance for random forces */
    std::vector<RealType>  randNums = genTriangleForces(nTriangles, 1.0);

    /* Loop over the mesh faces and apply random force to each of the faces*/    
    std::vector<Triangle>::iterator face;
    std::vector<StuntDouble*>::iterator vertex;
    int thisNumber = 0;
    for (face = sMesh.begin(); face != sMesh.end(); ++face){
     
      Triangle thisTriangle = *face;
      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
      RealType thisArea = thisTriangle.getArea(); 
      // RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*thisArea) /OOPSEConstant::energyConvert;
      // gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());

      /* Get Random Force */
      Vector3d unitNormal = thisTriangle.getNormal();
      unitNormal.normalize();
      //Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
      Vector3d centroid = thisTriangle.getCentroid();
      Vector3d facetVel = thisTriangle.getFacetVelocity();
      RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/NumericConstant::PI;

      RealType f_normal = viscosity_*hydroLength*OOPSEConstant::viscoConvert;
      RealType extPressure = -(targetPressure_ * thisArea)/OOPSEConstant::energyConvert;
      RealType randomForce = randNums[thisNumber++] * sqrt(2.0 * f_normal * OOPSEConstant::kb*targetTemp_/dt_);

      RealType dragForce = -f_normal * dot(facetVel, unitNormal);


      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
      
      //      Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
      
      // std::cout << " " << randomForce << " " << f_normal <<   std::endl;

      /* Apply triangle force to stuntdouble vertices */
      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
        if ((*vertex) != NULL){
          // mass = integrableObject->getMass();
          Vector3d vertexForce = langevinForce/3.0;
          (*vertex)->addFrc(vertexForce);

          if ((*vertex)->isDirectional()){

            Vector3d vertexPos = (*vertex)->getPos();
            Vector3d vertexCentroidVector = vertexPos - centroid;
            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
          }
        }  
      }
    } 

    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    currSnapshot->setVolume(surfaceMesh_->getVolume());    
    ForceManager::postCalculation(needStress);   
  }

  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {

    
    Vector<RealType, 6> Z;
    Vector<RealType, 6> generalForce;
    

    Z[0] = randNumGen_.randNorm(0, variance);
    Z[1] = randNumGen_.randNorm(0, variance);
    Z[2] = randNumGen_.randNorm(0, variance);
    Z[3] = randNumGen_.randNorm(0, variance);
    Z[4] = randNumGen_.randNorm(0, variance);
    Z[5] = randNumGen_.randNorm(0, variance);
     
    generalForce = hydroProps_[index]->getS()*Z;
    
    force[0] = generalForce[0];
    force[1] = generalForce[1];
    force[2] = generalForce[2];
    torque[0] = generalForce[3];
    torque[1] = generalForce[4];
    torque[2] = generalForce[5];
    
} 
  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {

    // zero fill the random vector before starting:
    std::vector<RealType> gaussRand;
    gaussRand.resize(nTriangles);
    std::fill(gaussRand.begin(), gaussRand.end(), 0.0);

    
#ifdef IS_MPI
    if (worldRank == 0) {
#endif
      for (int i = 0; i < nTriangles; i++) {
        //gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));     
        gaussRand[i] = randNumGen_.randNorm(0.0, 1.0);     
      }
#ifdef IS_MPI
    }
#endif

    // push these out to the other processors

#ifdef IS_MPI
    if (worldRank == 0) {
      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
    } else {
      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
    }
#endif
   
    for (int i = 0; i < nTriangles; i++) {
      gaussRand[i] = gaussRand[i] * variance;
    }
   
    return gaussRand;
  }


}
Revision:	1323
Committed:	Tue Nov 25 16:32:47 2008 UTC (16 years, 5 months ago) by chuckv
Original Path:	*trunk/src/integrators/SMIPDForceManager.cpp*
File size:	10696 byte(s)
Log Message:	Code clean-up. Removed old langevin code.
#	User	Rev	Content
1	chuckv	1293	/*
2			* Copyright (c) 2008 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
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41			#include <fstream>
42			#include <iostream>
43			#include "integrators/SMIPDForceManager.hpp"
44			#include "math/CholeskyDecomposition.hpp"
45			#include "utils/OOPSEConstant.hpp"
46			#include "hydrodynamics/Sphere.hpp"
47			#include "hydrodynamics/Ellipsoid.hpp"
48			#include "utils/ElementsTable.hpp"
49			#include "math/ConvexHull.hpp"
50			#include "math/Triangle.hpp"
51
52
53			namespace oopse {
54	chuckv	1323
55	chuckv	1293	SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
56			simParams = info->getSimParams();
57			veloMunge = new Velocitizer(info);
58
59			// Create Hull, Convex Hull for now, other options later.
60			surfaceMesh_ = new ConvexHull();
61	chuckv	1307
62	chuckv	1293
63
64			/* Check that the simulation has target pressure and target
65			temperature set*/
66	chuckv	1323
67	chuckv	1293	if (!simParams->haveTargetTemp()) {
68			sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
69			painCave.isFatal = 1;
70			painCave.severity = OOPSE_ERROR;
71			simError();
72			} else {
73			targetTemp_ = simParams->getTargetTemp();
74			}
75	chuckv	1323
76	chuckv	1293	if (!simParams->haveTargetPressure()) {
77			sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
78			" without a targetPressure!\n");
79
80			painCave.isFatal = 1;
81			simError();
82			} else {
83	chuckv	1323	/* Convert pressure from atm -> amu/(fs^2Ang)/
84	chuckv	1307	targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
85	chuckv	1293	}
86
87
88			if (simParams->getUsePeriodicBoundaryConditions()) {
89			sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
90			" with periodic boundary conditions !\n");
91
92			painCave.isFatal = 1;
93			simError();
94			}
95
96	chuckv	1316	if (!simParams->haveViscosity()) {
97			sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
98			painCave.isFatal = 1;
99			painCave.severity = OOPSE_ERROR;
100			simError();
101	chuckv	1323	}else{
102			viscosity_ = simParams->getViscosity();
103	chuckv	1316	}
104	chuckv	1293
105	chuckv	1323	dt_ = simParams->getDt();
106	chuckv	1307
107
108
109			//Compute initial hull
110			/*
111			surfaceMesh_->computeHull(localSites_);
112			Area0_ = surfaceMesh_->getArea();
113			int nTriangles = surfaceMesh_->getNMeshElements();
114			// variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
115			gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
116			(4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
117			//RealType eta0 = gamma_0/
118			*/
119
120	chuckv	1293
121			Molecule* mol;
122			StuntDouble* integrableObject;
123			SimInfo::MoleculeIterator i;
124			Molecule::IntegrableObjectIterator j;
125	chuckv	1323
126	chuckv	1293
127
128			/* Compute hull first time through to get the area of t=0*/
129
130	chuckv	1307	//Build a vector of integrable objects to determine if the are surface atoms
131	chuckv	1293	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
132			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
133			integrableObject = mol->nextIntegrableObject(j)) {
134			localSites_.push_back(integrableObject);
135			}
136			}
137
138	chuckv	1307
139	chuckv	1293	}
140
141			std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
142			std::map<std::string, HydroProp*> props;
143			std::ifstream ifs(filename.c_str());
144			if (ifs.is_open()) {
145
146			}
147
148			const unsigned int BufferSize = 65535;
149			char buffer[BufferSize];
150			while (ifs.getline(buffer, BufferSize)) {
151			HydroProp* currProp = new HydroProp(buffer);
152			props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
153			}
154
155			return props;
156			}
157
158			void SMIPDForceManager::postCalculation(bool needStress){
159			SimInfo::MoleculeIterator i;
160			Molecule::IntegrableObjectIterator j;
161			Molecule* mol;
162			StuntDouble* integrableObject;
163			RealType mass;
164			Vector3d pos;
165			Vector3d frc;
166			Mat3x3d A;
167			Mat3x3d Atrans;
168			Vector3d Tb;
169			Vector3d ji;
170			unsigned int index = 0;
171			int fdf;
172
173			fdf = 0;
174	chuckv	1301
175	chuckv	1293	/Compute surface Mesh/
176			surfaceMesh_->computeHull(localSites_);
177
178			/* Get area and number of surface stunt doubles and compute new variance */
179	chuckv	1302	RealType area = surfaceMesh_->getArea();
180			int nSurfaceSDs = surfaceMesh_->getNs();
181	chuckv	1293
182	chuckv	1307
183	chuckv	1308	std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
184	chuckv	1307	int nTriangles = sMesh.size();
185
186
187
188	chuckv	1302	/* Compute variance for random forces */
189	chuckv	1316	std::vector<RealType> randNums = genTriangleForces(nTriangles, 1.0);
190
191	chuckv	1323	/* Loop over the mesh faces and apply random force to each of the faces*/
192	chuckv	1308	std::vector<Triangle>::iterator face;
193	chuckv	1302	std::vector<StuntDouble*>::iterator vertex;
194			int thisNumber = 0;
195	chuckv	1293	for (face = sMesh.begin(); face != sMesh.end(); ++face){
196
197	chuckv	1308	Triangle thisTriangle = *face;
198			std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
199	chuckv	1316	RealType thisArea = thisTriangle.getArea();
200			// RealType sigma_t = sqrt(NumericConstant::PI/2.0)((targetPressure_)thisArea) /OOPSEConstant::energyConvert;
201			// gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());
202
203	chuckv	1302	/* Get Random Force */
204	chuckv	1308	Vector3d unitNormal = thisTriangle.getNormal();
205	chuckv	1302	unitNormal.normalize();
206	chuckv	1316	//Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
207	chuckv	1308	Vector3d centroid = thisTriangle.getCentroid();
208	chuckv	1316	Vector3d facetVel = thisTriangle.getFacetVelocity();
209	chuckv	1323	RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/NumericConstant::PI;
210	chuckv	1304
211	chuckv	1323	RealType f_normal = viscosity_hydroLengthOOPSEConstant::viscoConvert;
212	chuckv	1318	RealType extPressure = -(targetPressure_ * thisArea)/OOPSEConstant::energyConvert;
213	chuckv	1323	RealType randomForce = randNums[thisNumber++] * sqrt(2.0 * f_normal * OOPSEConstant::kb*targetTemp_/dt_);
214	chuckv	1320
215			RealType dragForce = -f_normal * dot(facetVel, unitNormal);
216
217
218	chuckv	1316	Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
219	chuckv	1307
220	chuckv	1316	// Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
221
222	chuckv	1320	// std::cout << " " << randomForce << " " << f_normal << std::endl;
223	chuckv	1316
224	chuckv	1323	/* Apply triangle force to stuntdouble vertices */
225	chuckv	1302	for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
226	chuckv	1307	if ((*vertex) != NULL){
227			// mass = integrableObject->getMass();
228			Vector3d vertexForce = langevinForce/3.0;
229			(*vertex)->addFrc(vertexForce);
230	chuckv	1312
231			if ((*vertex)->isDirectional()){
232	chuckv	1316
233	chuckv	1307	Vector3d vertexPos = (*vertex)->getPos();
234			Vector3d vertexCentroidVector = vertexPos - centroid;
235			(*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
236			}
237			}
238	chuckv	1293	}
239	chuckv	1302	}
240	chuckv	1293
241	chuckv	1302	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
242	chuckv	1316	currSnapshot->setVolume(surfaceMesh_->getVolume());
243	chuckv	1302	ForceManager::postCalculation(needStress);
244	chuckv	1293	}
245
246	chuckv	1302	void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
247	chuckv	1293
248	chuckv	1302
249	chuckv	1293	Vector<RealType, 6> Z;
250			Vector<RealType, 6> generalForce;
251	chuckv	1302
252
253	chuckv	1293	Z[0] = randNumGen_.randNorm(0, variance);
254			Z[1] = randNumGen_.randNorm(0, variance);
255			Z[2] = randNumGen_.randNorm(0, variance);
256			Z[3] = randNumGen_.randNorm(0, variance);
257			Z[4] = randNumGen_.randNorm(0, variance);
258			Z[5] = randNumGen_.randNorm(0, variance);
259
260			generalForce = hydroProps_[index]->getS()*Z;
261
262			force[0] = generalForce[0];
263			force[1] = generalForce[1];
264			force[2] = generalForce[2];
265			torque[0] = generalForce[3];
266			torque[1] = generalForce[4];
267			torque[2] = generalForce[5];
268
269			}
270			std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
271
272			// zero fill the random vector before starting:
273			std::vector<RealType> gaussRand;
274			gaussRand.resize(nTriangles);
275			std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
276
277	chuckv	1320
278	chuckv	1293
279			#ifdef IS_MPI
280			if (worldRank == 0) {
281			#endif
282			for (int i = 0; i < nTriangles; i++) {
283	chuckv	1316	//gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));
284			gaussRand[i] = randNumGen_.randNorm(0.0, 1.0);
285	chuckv	1293	}
286			#ifdef IS_MPI
287			}
288			#endif
289
290			// push these out to the other processors
291
292			#ifdef IS_MPI
293			if (worldRank == 0) {
294			MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
295			} else {
296			MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
297			}
298			#endif
299	chuckv	1316
300	chuckv	1293	for (int i = 0; i < nTriangles; i++) {
301			gaussRand[i] = gaussRand[i] * variance;
302			}
303	chuckv	1316
304	chuckv	1293	return gaussRand;
305			}
306
307	chuckv	1306
308	chuckv	1293	}