--- branches/development/src/selection/HullFinder.cpp 2012/10/03 14:20:07 1803 +++ branches/development/src/selection/HullFinder.cpp 2013/02/20 15:39:39 1850 @@ -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). */ @@ -123,4 +123,39 @@ namespace OpenMD { return bsResult; } + OpenMDBitSet HullFinder::findHull(int frame) { + Snapshot* currSnapshot = info_->getSnapshotManager()->getSnapshot(frame); + OpenMDBitSet bsResult(nStuntDoubles_); +#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 sMesh = surfaceMesh_->getMesh(); + int nTriangles = sMesh.size(); + // Loop over the mesh faces + std::vector::iterator face; + std::vector::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 vertexSDs = thisTriangle.getVertices(); + for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) { + if ((*vertex) != NULL) { + bsResult.setBitOn((*vertex)->getGlobalIndex()); + } + } + } + return bsResult; + } + }