--- trunk/src/math/ConvexHull.cpp	2008/10/15 18:26:01	1304
+++ trunk/src/math/ConvexHull.cpp	2008/10/21 16:44:00	1308
@@ -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.10 2008-10-15 18:26:01 chuckv Exp $
+ *  @version $Id: ConvexHull.cpp,v 1.12 2008-10-21 16:44:00 chuckv Exp $
  *
  */
 
@@ -370,7 +370,7 @@ void ConvexHull::printHull(const std::string& geomFile
 #ifdef HAVE_QHULL
 /* Old options Qt Qu Qg QG0 FA */
 /* More old opts Qc Qi Pp*/
-ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp"), Ns_(200) {
+ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp"), Ns_(200), nTriangles_(0) {
   //If we are doing the mpi version, set up some vectors for data communication
 #ifdef IS_MPI
 
@@ -438,6 +438,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
   
   boolT ismalloc = False;
   /* Clean up memory from previous convex hull calculations*/
+  
   Triangles_.clear();
   surfaceSDs_.clear();
   surfaceSDs_.reserve(Ns_);
@@ -446,7 +447,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
    		    const_cast<char *>(options_.c_str()), NULL, stderr)) {
 
       sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull");
-      painCave.isFatal = 0;
+      painCave.isFatal = 1;
       simError();
       
   } //qh_new_qhull
@@ -454,6 +455,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
 
 #ifdef IS_MPI
   std::vector<double> localPts;
+  std::vector<double> localVel;
   int localPtArraySize;
   
  
@@ -465,7 +467,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
     if (!facet->simplicial){
       // should never happen with Qt
       sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
-      painCave.isFatal = 0;
+      painCave.isFatal = 1;
       simError();
     }
     
@@ -517,6 +519,12 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
     localPts.push_back(ptArray[dim_ * idx]);     
     localPts.push_back(ptArray[dim_ * idx + 1]); 
     localPts.push_back(ptArray[dim_ * idx + 2]);
+
+    Vector3d vel = bodydoubles[idx]->getVel();
+    localVel.push_back(vel.x());
+    localVel.push_back(vel.y());
+    localVel.push_back(vel.z());
+
     localPtsMap.push_back(idx); 
   }
 
@@ -535,8 +543,8 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
   int nglobalPts = int(Nsglobal_/3);
  
 
-  std::vector<double> globalPts;
-  globalPts.resize(Nsglobal_);
+  std::vector<double> globalPts(Nsglobal_);
+  std::vector<double> globalVel(Nsglobal_);
 
   isSurfaceID.resize(nglobalPts);
 
@@ -554,6 +562,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
   /* gather the potential hull */
   
   MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
+  MPI::COMM_WORLD.Allgatherv(&localVel[0],localPtArraySize,MPI::DOUBLE,&globalVel[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
 
   /*
   if (myrank_ == 0){
@@ -600,31 +609,34 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
 	simError();
       } //simplicical
       
-      Triangle* face = new Triangle();
+      Triangle face;
       Vector3d  V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]);
-      face->setNormal(V3dNormal);
+      face.setNormal(V3dNormal);
  
       
 
       RealType faceArea = 0.5*V3dNormal.length();
-      face->setArea(faceArea);
+      face.setArea(faceArea);
 
 
       vertices = qh_facet3vertex(facet);
       
       coordT *center = qh_getcenter(vertices);
       Vector3d V3dCentroid(center[0], center[1], center[2]);
-      face->setCentroid(V3dCentroid);
-
+      face.setCentroid(V3dCentroid);
+      Vector3d faceVel = V3Zero;
       FOREACHvertex_(vertices){
 	id = qh_pointid(vertex->point);
 	int localindex = id;
 #ifdef IS_MPI
-	
+	Vector3d velVector(globalVel[dim_ * id],globalVel[dim_ * id + 1], globalVel[dim_ * id + 1]);
+	faceVel = faceVel + velVector;
 	if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){
 	  localindex = localPtsMap[id-noffset/3];
+#else
+	  faceVel = faceVel + bodydoubles[localindex]->getVel();
 #endif
-	  face->addVertex(bodydoubles[localindex]);
+	  face.addVertex(bodydoubles[localindex]);
 	  if( !isSurfaceID[id] ){
 	    isSurfaceID[id] = true;
 #ifdef IS_MPI	    
@@ -638,7 +650,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
 #ifdef IS_MPI
 	 
 	}else{
-	  face->addVertex(NULL);
+	  face.addVertex(NULL);
 	  }
 #endif
       } //Foreachvertex
@@ -649,10 +661,11 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
 	  surfaceSDs_.push_back(bodydoubles[id]);
 	}
       }
-
+      */
+      face.setFacetVelocity(faceVel/3.0);
       Triangles_.push_back(face);
       qh_settempfree(&vertices);      
-      */
+
     } //FORALLfacets
 
     /*
@@ -666,8 +679,8 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
 
 
     Ns_ = surfaceSDs_.size();
+    nTriangles_ = Triangles_.size();
     
-    
     qh_getarea(qh facet_list);
     volume_ = qh totvol;
     area_ = qh totarea;