--- trunk/src/math/ConvexHull.cpp	2008/10/07 17:12:48	1302
+++ 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.9 2008-10-07 17:12:48 chuckv Exp $
+ *  @version $Id: ConvexHull.cpp,v 1.12 2008-10-21 16:44:00 chuckv Exp $
  *
  */
 
@@ -369,7 +369,8 @@ void ConvexHull::printHull(const std::string& geomFile
 #else
 #ifdef HAVE_QHULL
 /* Old options Qt Qu Qg QG0 FA */
-ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt  Qci Tcv Pp"), Ns_(200) {
+/* More old opts Qc Qi Pp*/
+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
 
@@ -437,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_);
@@ -445,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
@@ -453,6 +455,7 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
 
 #ifdef IS_MPI
   std::vector<double> localPts;
+  std::vector<double> localVel;
   int localPtArraySize;
   
  
@@ -464,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();
     }
     
@@ -516,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); 
   }
 
@@ -534,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);
 
@@ -553,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){
@@ -599,24 +609,34 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
 	simError();
       } //simplicical
       
-      Triangle* face = new Triangle();
-      Vector3d  V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]);
-      face->setNormal(V3dNormal);
-      //face->setCentroid(V3dCentroid);
+      Triangle face;
+      Vector3d  V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]);
+      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);
+      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	    
@@ -630,35 +650,51 @@ void ConvexHull::computeHull(std::vector<StuntDouble*>
 #ifdef IS_MPI
 	 
 	}else{
-	  face->addVertex(NULL);
+	  face.addVertex(NULL);
 	  }
 #endif
       } //Foreachvertex
-
+      /*
+      if (!SETempty_(facet->coplanarset)){
+	FOREACHpoint_(facet->coplanarset){
+	  id = qh_pointid(point);
+	  surfaceSDs_.push_back(bodydoubles[id]);
+	}
+      }
+      */
+      face.setFacetVelocity(faceVel/3.0);
       Triangles_.push_back(face);
       qh_settempfree(&vertices);      
-      
+
     } //FORALLfacets
 
-        		
+    /*
+    std::cout << surfaceSDs_.size() << std::endl;
+    for (SD = surfaceSDs_.begin(); SD != surfaceSDs_.end(); ++SD){
+      Vector3d thisatom = (*SD)->getPos();
+      std::cout << "Au " << thisatom.x() << "  " << thisatom.y() << " " << thisatom.z() << std::endl;
+    }
+    */
 
-  Ns_ = surfaceSDs_.size();
 
 
-  qh_getarea(qh facet_list);
-  volume_ = qh totvol;
-  area_ = qh totarea;
-
-  
-
-  qh_freeqhull(!qh_ALL);
-  qh_memfreeshort(&curlong, &totlong);
-  if (curlong || totlong)
-    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
-	      << totlong << curlong << std::endl;
-  
-
-  
+    Ns_ = surfaceSDs_.size();
+    nTriangles_ = Triangles_.size();
+    
+    qh_getarea(qh facet_list);
+    volume_ = qh totvol;
+    area_ = qh totarea;
+    
+    
+    
+    qh_freeqhull(!qh_ALL);
+    qh_memfreeshort(&curlong, &totlong);
+    if (curlong || totlong)
+      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
+		<< totlong << curlong << std::endl;
+    
+    
+    
 }