--- trunk/src/math/ConvexHull.cpp	2009/10/20 20:05:28	1374
+++ trunk/src/math/ConvexHull.cpp	2009/10/22 19:43:10	1384
@@ -44,7 +44,7 @@
  *
  *  Created by Charles F. Vardeman II on 11 Dec 2006.
  *  @author  Charles F. Vardeman II
- *  @version $Id: ConvexHull.cpp,v 1.16 2009-10-20 20:05:28 chuckv Exp $
+ *  @version $Id: ConvexHull.cpp,v 1.20 2009-10-22 19:43:10 gezelter Exp $
  *
  */
 
@@ -77,9 +77,6 @@ extern "C"
 #include <qhull/stat.h>
 }
 
-/* Old options Qt Qu Qg QG0 FA */
-/* More old opts Qc Qi Pp*/
-
 ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp") {
 }
 
@@ -94,8 +91,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
   setT *vertices;
   int curlong, totlong;
   
-  std::vector<double> ptArray(numpoints*3);
-  std::vector<bool> isSurfaceID(numpoints); 
+  std::vector<double> ptArray(numpoints*dim_);
 
   // Copy the positon vector into a points vector for qhull.
   std::vector<StuntDouble*>::iterator SD;
@@ -127,20 +123,24 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
   int nproc = MPI::COMM_WORLD.Get_size();
   int myrank = MPI::COMM_WORLD.Get_rank();
   int localHullSites = 0;
-  int* hullSitesOnProc = new int[nproc];
-  int* coordsOnProc = new int[nproc];
-  int* displacements = new int[nproc];
-  int* vectorDisplacements = new int[nproc];
 
+  std::vector<int> hullSitesOnProc(nproc, 0);
+  std::vector<int> coordsOnProc(nproc, 0);
+  std::vector<int> displacements(nproc, 0);
+  std::vector<int> vectorDisplacements(nproc, 0);
+
   std::vector<double> coords;
   std::vector<double> vels;
-  std::vector<int> objectIDs;
+  std::vector<int> indexMap;
   std::vector<double> masses;
 
   FORALLvertices{
     localHullSites++;
     
     int idx = qh_pointid(vertex->point);
+
+    indexMap.push_back(idx);
+
     coords.push_back(ptArray[dim_  * idx]);
     coords.push_back(ptArray[dim_  * idx + 1]);
     coords.push_back(ptArray[dim_  * idx + 2]);
@@ -159,7 +159,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
                             1, MPI::INT);
 
   int globalHullSites = 0;
-  for (int iproc = 0; i < nproc; iproc++){
+  for (int iproc = 0; iproc < nproc; iproc++){
     globalHullSites += hullSitesOnProc[iproc];
     coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc];
   }
@@ -167,27 +167,28 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
   displacements[0] = 0;
   vectorDisplacements[0] = 0;
   
-  for (int iproc = 1; i < nproc; iproc++){
+  for (int iproc = 1; iproc < nproc; iproc++){
     displacements[iproc] = displacements[iproc-1] + hullSitesOnProc[iproc-1];
     vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1]; 
   }
 
-  std::vector<double> globalCoords(dim_*globalHullSites);
-  std::vector<double> globalVels(dim_*globalHullSites);
+  std::vector<double> globalCoords(dim_ * globalHullSites);
+  std::vector<double> globalVels(dim_ * globalHullSites);
   std::vector<double> globalMasses(globalHullSites);
+
   int count = coordsOnProc[myrank];
   
-  MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE,
-                             &globalCoords[0], &coordsOnProc[0], &vectorDisplacements[0],
+  MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0],
+                             &coordsOnProc[0], &vectorDisplacements[0], 
                              MPI::DOUBLE);
 
-  MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE,
-                             &globalVels[0], &coordsOnProc[0], &vectorDisplacements[0],
+  MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0], 
+                             &coordsOnProc[0], &vectorDisplacements[0],
                              MPI::DOUBLE);
 
   MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE,
-                             &globalMasses[0], &hullSitesOnProc[0], &displacements[0],
-                             MPI::DOUBLE);
+                             &globalMasses[0], &hullSitesOnProc[0], 
+                             &displacements[0], MPI::DOUBLE);
 
   // Free previous hull
   qh_freeqhull(!qh_ALL);
@@ -225,6 +226,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
     Vector3d faceVel = V3Zero;
     Vector3d p[3];
     RealType faceMass = 0.0;
+
     int ver = 0;
 
     FOREACHvertex_(vertices){
@@ -242,12 +244,14 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
                      globalVels[dim_ * id + 2]);
       mass = globalMasses[id];
 
-      // localID will be between 0 and hullSitesOnProc[myrank] if we own this guy.
+      // localID will be between 0 and hullSitesOnProc[myrank] if we
+      // own this guy.
+
       int localID = id - displacements[myrank];
-      if (id >= 0 && id < hullSitesOnProc[myrank])
-        face.addVertexSD(bodydoubles[localID]);
-      else
-        face.addVertexSD(NULL);
+
+      if (localID >= 0 && localID < hullSitesOnProc[myrank])
+        face.addVertexSD(bodydoubles[indexMap[localID]]);
+      
 #else
       vel = bodydoubles[id]->getVel();
       mass = bodydoubles[id]->getMass();
@@ -271,13 +275,6 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
   volume_ = qh totvol;
   area_ = qh totarea;
 
-#ifdef IS_MPI
-  delete [] hullSitesOnProc;
-  delete [] coordsOnProc;
-  delete [] displacements;
-  delete [] vectorDisplacements;
-#endif
-  
   qh_freeqhull(!qh_ALL);
   qh_memfreeshort(&curlong, &totlong);
   if (curlong || totlong)
@@ -285,9 +282,11 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
               << totlong << curlong << std::endl;    
 }
 
-
-
 void ConvexHull::printHull(const std::string& geomFileName) {
+
+#ifdef IS_MPI
+  if (worldRank == 0)  {
+#endif
   FILE *newGeomFile;
   
   //create new .md file based on old .md file
@@ -297,5 +296,8 @@ void ConvexHull::printHull(const std::string& geomFile
     qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
   
   fclose(newGeomFile);
+#ifdef IS_MPI
+  }
+#endif
 }
 #endif //QHULL