ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/selection/HullFinder.cpp
(Generate patch)

Comparing trunk/src/selection/HullFinder.cpp (file contents):
Revision 1801 by gezelter, Mon Oct 1 18:21:15 2012 UTC vs.
Revision 2052 by gezelter, Fri Jan 9 19:06:35 2015 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 48 | Line 48 | namespace OpenMD {
48  
49    HullFinder::HullFinder(SimInfo* info) : info_(info) {
50  
51 <    nStuntDoubles_ = info_->getNGlobalAtoms() + info_->getNGlobalRigidBodies();
52 <    stuntdoubles_.resize(nStuntDoubles_);
51 >    nObjects_.push_back(info_->getNGlobalAtoms()+info_->getNGlobalRigidBodies());
52 >    nObjects_.push_back(info_->getNGlobalBonds());
53 >    nObjects_.push_back(info_->getNGlobalBends());
54 >    nObjects_.push_back(info_->getNGlobalTorsions());
55 >    nObjects_.push_back(info_->getNGlobalInversions());
56 >    nObjects_.push_back(info_->getNGlobalMolecules());
57 >
58 >    stuntdoubles_.resize(nObjects_[STUNTDOUBLE]);
59 >    bonds_.resize(nObjects_[BOND]);
60 >    bends_.resize(nObjects_[BEND]);
61 >    torsions_.resize(nObjects_[TORSION]);
62 >    inversions_.resize(nObjects_[INVERSION]);
63 >    molecules_.resize(nObjects_[MOLECULE]);
64      
65      SimInfo::MoleculeIterator mi;
66 +    Molecule::IntegrableObjectIterator ioi;
67 +    Molecule::AtomIterator ai;
68 +    Molecule::RigidBodyIterator rbIter;
69 +    Molecule::BondIterator bondIter;
70 +    Molecule::BendIterator bendIter;
71 +    Molecule::TorsionIterator torsionIter;
72 +    Molecule::InversionIterator inversionIter;
73 +    
74      Molecule* mol;
75      StuntDouble* sd;
57 –    Molecule::IntegrableObjectIterator  ioi;
58 –    Molecule::AtomIterator ai;
76      Atom* atom;
60 –    Molecule::RigidBodyIterator rbIter;
77      RigidBody* rb;
78 <    
78 >    Bond* bond;
79 >    Bend* bend;
80 >    Torsion* torsion;
81 >    Inversion* inversion;    
82 >
83      for (mol = info_->beginMolecule(mi); mol != NULL;
84           mol = info_->nextMolecule(mi)) {
85 <        
85 >
86 >      molecules_[mol->getGlobalIndex()] = mol;
87 >      
88        // Hull is constructed from all known integrable objects.
89        for (sd = mol->beginIntegrableObject(ioi);
90             sd != NULL;
# Line 81 | Line 103 | namespace OpenMD {
103             rb = mol->nextRigidBody(rbIter)) {
104          stuntdoubles_[rb->getGlobalIndex()] = rb;
105        }
106 <        
106 >
107 >      // These others are going to be inferred from the objects on the hull:
108 >      for (bond = mol->beginBond(bondIter); bond != NULL;
109 >           bond = mol->nextBond(bondIter)) {
110 >        bonds_[bond->getGlobalIndex()] = bond;
111 >      }  
112 >      for (bend = mol->beginBend(bendIter); bend != NULL;
113 >           bend = mol->nextBend(bendIter)) {
114 >        bends_[bend->getGlobalIndex()] = bend;
115 >      }  
116 >      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
117 >           torsion = mol->nextTorsion(torsionIter)) {
118 >        torsions_[torsion->getGlobalIndex()] = torsion;
119 >      }  
120 >      for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
121 >           inversion = mol->nextInversion(inversionIter)) {
122 >        inversions_[inversion->getGlobalIndex()] = inversion;
123 >      }  
124 >
125      }
126   #ifdef HAVE_QHULL
127      surfaceMesh_ = new ConvexHull();
128   #endif
129    }
130  
131 <  OpenMDBitSet HullFinder::findHull() {
132 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
133 <    OpenMDBitSet bsResult(nStuntDoubles_);
131 >  HullFinder::~HullFinder() {
132 >    delete surfaceMesh_;
133 >  }
134 >
135 >  SelectionSet HullFinder::findHull() {
136 >    SelectionSet ssResult(nObjects_);
137   #ifdef HAVE_QHULL
138      surfaceMesh_->computeHull(localSites_);
139   #else
140      sprintf( painCave.errMsg,
141               "HullFinder : Hull calculation is not possible without libqhull.\n"
142               "\tPlease rebuild OpenMD with qhull enabled.");
143 +    painCave.severity = OPENMD_ERROR;
144 +    painCave.isFatal = 1;
145 +    simError();
146 + #endif
147 +    
148 +    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
149 +    // Loop over the mesh faces
150 +    std::vector<Triangle>::iterator face;
151 +    std::vector<StuntDouble*>::iterator vertex;
152 +
153 +    // This will work in parallel because the triangles returned by the mesh
154 +    // have a NULL stuntDouble if this processor doesn't own the
155 +    
156 +    for (face = sMesh.begin(); face != sMesh.end(); ++face) {
157 +      Triangle thisTriangle = *face;
158 +      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
159 +      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) {
160 +        if ((*vertex) != NULL) {
161 +          ssResult.bitsets_[STUNTDOUBLE].setBitOn((*vertex)->getGlobalIndex());
162 +        }
163 +      }
164 +    }
165 +    return ssResult;
166 +  }
167 +
168 +  SelectionSet HullFinder::findHull(int frame) {
169 +    SelectionSet ssResult(nObjects_);
170 + #ifdef HAVE_QHULL
171 +    surfaceMesh_->computeHull(localSites_);
172 + #else
173 +    sprintf( painCave.errMsg,
174 +             "HullFinder : Hull calculation is not possible without libqhull.\n"
175 +             "\tPlease rebuild OpenMD with qhull enabled.");
176        painCave.severity = OPENMD_ERROR;
177        painCave.isFatal = 1;
178        simError();
179   #endif
180      
181      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
106 –    int nTriangles = sMesh.size();
182      // Loop over the mesh faces
183      std::vector<Triangle>::iterator face;
184      std::vector<StuntDouble*>::iterator vertex;
# Line 116 | Line 191 | namespace OpenMD {
191        std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
192        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) {
193          if ((*vertex) != NULL) {
194 <          bsResult.setBitOn((*vertex)->getGlobalIndex());          
194 >          ssResult.bitsets_[STUNTDOUBLE].setBitOn((*vertex)->getGlobalIndex());
195          }
196        }
197      }
198 <    return bsResult;
198 >    surfaceArea_ = surfaceMesh_->getArea();
199 >    return ssResult;
200    }
201  
202   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines