--- trunk/src/selection/HullFinder.cpp	2012/10/01 18:21:15	1801
+++ trunk/src/selection/HullFinder.cpp	2015/02/20 15:12:07	2056
@@ -35,7 +35,7 @@
  *                                                                      
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
- * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
+ * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
  * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
  * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
@@ -48,21 +48,43 @@ namespace OpenMD {
 
   HullFinder::HullFinder(SimInfo* info) : info_(info) {
 
-    nStuntDoubles_ = info_->getNGlobalAtoms() + info_->getNGlobalRigidBodies();
-    stuntdoubles_.resize(nStuntDoubles_);
+    nObjects_.push_back(info_->getNGlobalAtoms()+info_->getNGlobalRigidBodies());
+    nObjects_.push_back(info_->getNGlobalBonds());
+    nObjects_.push_back(info_->getNGlobalBends());
+    nObjects_.push_back(info_->getNGlobalTorsions());
+    nObjects_.push_back(info_->getNGlobalInversions());
+    nObjects_.push_back(info_->getNGlobalMolecules());
+
+    stuntdoubles_.resize(nObjects_[STUNTDOUBLE]);
+    bonds_.resize(nObjects_[BOND]);
+    bends_.resize(nObjects_[BEND]);
+    torsions_.resize(nObjects_[TORSION]);
+    inversions_.resize(nObjects_[INVERSION]);
+    molecules_.resize(nObjects_[MOLECULE]);
     
     SimInfo::MoleculeIterator mi;
+    Molecule::IntegrableObjectIterator ioi;
+    Molecule::AtomIterator ai;
+    Molecule::RigidBodyIterator rbIter;
+    Molecule::BondIterator bondIter;
+    Molecule::BendIterator bendIter;
+    Molecule::TorsionIterator torsionIter;
+    Molecule::InversionIterator inversionIter;
+    
     Molecule* mol;
     StuntDouble* sd;
-    Molecule::IntegrableObjectIterator  ioi;
-    Molecule::AtomIterator ai;
     Atom* atom;
-    Molecule::RigidBodyIterator rbIter;
     RigidBody* rb;
-    
+    Bond* bond;
+    Bend* bend;
+    Torsion* torsion;
+    Inversion* inversion;    
+
     for (mol = info_->beginMolecule(mi); mol != NULL; 
          mol = info_->nextMolecule(mi)) {
-        
+
+      molecules_[mol->getGlobalIndex()] = mol;
+      
       // Hull is constructed from all known integrable objects.
       for (sd = mol->beginIntegrableObject(ioi);
            sd != NULL;
@@ -81,29 +103,83 @@ namespace OpenMD {
            rb = mol->nextRigidBody(rbIter)) {
 	stuntdoubles_[rb->getGlobalIndex()] = rb;
       }
-        
+
+      // These others are going to be inferred from the objects on the hull:
+      for (bond = mol->beginBond(bondIter); bond != NULL; 
+           bond = mol->nextBond(bondIter)) {
+        bonds_[bond->getGlobalIndex()] = bond;
+      }   
+      for (bend = mol->beginBend(bendIter); bend != NULL; 
+           bend = mol->nextBend(bendIter)) {
+        bends_[bend->getGlobalIndex()] = bend;
+      }   
+      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; 
+           torsion = mol->nextTorsion(torsionIter)) {
+        torsions_[torsion->getGlobalIndex()] = torsion;
+      }   
+      for (inversion = mol->beginInversion(inversionIter); inversion != NULL; 
+           inversion = mol->nextInversion(inversionIter)) {
+        inversions_[inversion->getGlobalIndex()] = inversion;
+      }   
+
     }
 #ifdef HAVE_QHULL
     surfaceMesh_ = new ConvexHull();
 #endif
   }
 
-  OpenMDBitSet HullFinder::findHull() {
-    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
-    OpenMDBitSet bsResult(nStuntDoubles_);
+  HullFinder::~HullFinder() {
+    delete surfaceMesh_;
+  }
+
+  SelectionSet HullFinder::findHull() {
+    SelectionSet ssResult(nObjects_);
 #ifdef HAVE_QHULL
     surfaceMesh_->computeHull(localSites_);
 #else
     sprintf( painCave.errMsg,
              "HullFinder : Hull calculation is not possible without libqhull.\n"
              "\tPlease rebuild OpenMD with qhull enabled.");
+    painCave.severity = OPENMD_ERROR;
+    painCave.isFatal = 1;
+    simError();
+#endif
+    
+    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
+    // Loop over the mesh faces
+    std::vector<Triangle>::iterator face;
+    std::vector<StuntDouble*>::iterator vertex;
+
+    // This will work in parallel because the triangles returned by the mesh
+    // have a NULL stuntDouble if this processor doesn't own the 
+    
+    for (face = sMesh.begin(); face != sMesh.end(); ++face) {
+      Triangle thisTriangle = *face;
+      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
+      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) {
+        if ((*vertex) != NULL) {
+          ssResult.bitsets_[STUNTDOUBLE].setBitOn((*vertex)->getGlobalIndex());
+        }
+      }
+    }
+    surfaceArea_ = surfaceMesh_->getArea();
+    return ssResult;
+  }
+
+  SelectionSet HullFinder::findHull(int frame) {
+    SelectionSet ssResult(nObjects_);
+#ifdef HAVE_QHULL
+    surfaceMesh_->computeHull(localSites_);
+#else
+    sprintf( painCave.errMsg,
+             "HullFinder : Hull calculation is not possible without libqhull.\n"
+             "\tPlease rebuild OpenMD with qhull enabled.");
       painCave.severity = OPENMD_ERROR;
       painCave.isFatal = 1;
       simError();
 #endif
     
     std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
-    int nTriangles = sMesh.size();
     // Loop over the mesh faces
     std::vector<Triangle>::iterator face;
     std::vector<StuntDouble*>::iterator vertex;
@@ -116,11 +192,12 @@ namespace OpenMD {
       std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
       for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) {
         if ((*vertex) != NULL) {
-          bsResult.setBitOn((*vertex)->getGlobalIndex());          
+          ssResult.bitsets_[STUNTDOUBLE].setBitOn((*vertex)->getGlobalIndex());
         }
       }
     }
-    return bsResult;
+    surfaceArea_ = surfaceMesh_->getArea();
+    return ssResult;
   }
 
 }