--- trunk/src/selection/HullFinder.cpp 2010/03/22 18:45:39 1412 +++ trunk/src/selection/HullFinder.cpp 2012/10/01 18:21:15 1801 @@ -36,7 +36,8 @@ * [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). - * [4] Vardeman & Gezelter, in progress (2009). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include "selection/HullFinder.hpp" @@ -52,7 +53,7 @@ namespace OpenMD { SimInfo::MoleculeIterator mi; Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; Molecule::IntegrableObjectIterator ioi; Molecule::AtomIterator ai; Atom* atom; @@ -63,10 +64,10 @@ namespace OpenMD { mol = info_->nextMolecule(mi)) { // Hull is constructed from all known integrable objects. - for (integrableObject = mol->beginIntegrableObject(ioi); - integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ioi)) { - localSites_.push_back(integrableObject); + for (sd = mol->beginIntegrableObject(ioi); + sd != NULL; + sd = mol->nextIntegrableObject(ioi)) { + localSites_.push_back(sd); } // selection can include atoms (which may be a subset of the IOs) @@ -81,22 +82,35 @@ namespace OpenMD { stuntdoubles_[rb->getGlobalIndex()] = rb; } - } - surfaceMesh_ = new ConvexHull(); + } +#ifdef HAVE_QHULL + surfaceMesh_ = new ConvexHull(); +#endif } OpenMDBitSet HullFinder::findHull() { - StuntDouble* sd; Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); 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();