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

Comparing trunk/src/math/AlphaHull.cpp (file contents):
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC vs.
Revision 1938 by gezelter, Thu Oct 31 15:32:17 2013 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 2008, 2009, 2010 The University of Notre Dame. All Rights Reserved.
1 > /* Copyright (c) 2010 The University of Notre Dame. All Rights Reserved.
2   *
3   * The University of Notre Dame grants you ("Licensee") a
4   * non-exclusive, royalty free, license to use, modify and
# Line 32 | Line 32
32   * research, please cite the appropriate papers when you publish your
33   * work.  Good starting points are:
34   *                                                                      
35 < * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
36 < * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
37 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
35 > * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
36 > * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
37 > * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
38   * [4] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
39 < * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
39 > * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
40   *
41   *  AlphaHull.cpp
42   *
43 < *  Purpose: To calculate Alpha hull, hull volume libqhull.
44 < *
45 < *  Created by Charles F. Vardeman II on 11 Dec 2006.
46 < *  @author  Charles F. Vardeman II
47 < *  @version $Id$
48 < *
43 > *  Purpose: To calculate an alpha-shape hull.
44   */
45  
46 + #ifdef IS_MPI
47 + #include <mpi.h>
48 + #endif
49 +
50   /* Standard includes independent of library */
51  
52   #include <iostream>
# Line 55 | Line 54
54   #include <list>
55   #include <algorithm>
56   #include <iterator>
58 – #include <utility>
57   #include "math/AlphaHull.hpp"
58   #include "utils/simError.h"
59 < #ifdef IS_MPI
62 < #include <mpi.h>
63 < #endif
59 >
60   #include "math/qhull.hpp"
61  
62 + #ifdef HAVE_QHULL
63 + using namespace std;
64   using namespace OpenMD;
65  
66 < #ifdef HAVE_QHULL
69 < double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim);
66 > double calculate_circumradius(pointT* p0, pointT* p1, pointT* p2, int dim);
67  
68 < AlphaHull::AlphaHull(double alpha) : Hull(), dim_(3), alpha_(alpha), options_("qhull d QJ Tcv Pp") {
68 > AlphaHull::AlphaHull(double alpha) : Hull(), dim_(3), alpha_(alpha),
69 >                                     options_("qhull d QJ Tcv Pp") {
70   }
71  
72 < void AlphaHull::computeHull(std::vector<StuntDouble*> bodydoubles) {
72 > void AlphaHull::computeHull(vector<StuntDouble*> bodydoubles) {
73  
74    int numpoints = bodydoubles.size();
77 –  bool alphashape=true;
75    
76    Triangles_.clear();
77    
78 <  vertexT *vertex, **vertexp;
78 >  vertexT *vertex;
79    facetT *facet, *neighbor;
83 –  setT *vertices, *verticestop, *verticesbottom;
84 –  int curlong, totlong;
80    pointT *interiorPoint;
81 +  int curlong, totlong;
82    
83 <  std::vector<double> ptArray(numpoints*dim_);
84 <
83 >  
84 >  vector<double> ptArray(numpoints*dim_);
85 >  
86    // Copy the positon vector into a points vector for qhull.
87 <  std::vector<StuntDouble*>::iterator SD;
87 >  vector<StuntDouble*>::iterator SD;
88    int i = 0;
89    for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){
90      Vector3d pos = (*SD)->getPos();      
# Line 96 | Line 93 | void AlphaHull::computeHull(std::vector<StuntDouble*>
93      ptArray[dim_ * i + 2] = pos.z();
94      i++;
95    }
96 <
97 <    /* Clean up memory from previous convex hull calculations*/
96 >  
97 >  /* Clean up memory from previous convex hull calculations*/
98    boolT ismalloc = False;
99 <
100 <  int ridgesCount=0;
99 >  
100 >  /* compute the hull for our local points (or all the points for single
101 >     processor versions) */
102    if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc,
103                     const_cast<char *>(options_.c_str()), NULL, stderr)) {
104 <
104 >    
105      sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute convex hull");
106      painCave.isFatal = 1;
107      simError();
108      
109    } //qh_new_qhull
110 <
111 <
110 >  
111 >  
112   #ifdef IS_MPI
113    //If we are doing the mpi version, set up some vectors for data communication
114    
# Line 118 | Line 116 | void AlphaHull::computeHull(std::vector<StuntDouble*>
116    int myrank = MPI::COMM_WORLD.Get_rank();
117    int localHullSites = 0;
118  
119 <  std::vector<int> hullSitesOnProc(nproc, 0);
120 <  std::vector<int> coordsOnProc(nproc, 0);
121 <  std::vector<int> displacements(nproc, 0);
122 <  std::vector<int> vectorDisplacements(nproc, 0);
119 >  vector<int> hullSitesOnProc(nproc, 0);
120 >  vector<int> coordsOnProc(nproc, 0);
121 >  vector<int> displacements(nproc, 0);
122 >  vector<int> vectorDisplacements(nproc, 0);
123  
124 <  std::vector<double> coords;
125 <  std::vector<double> vels;
126 <  std::vector<int> indexMap;
127 <  std::vector<double> masses;
124 >  vector<double> coords;
125 >  vector<double> vels;
126 >  vector<int> indexMap;
127 >  vector<double> masses;
128  
129    FORALLvertices{
130      localHullSites++;
# Line 166 | Line 164 | void AlphaHull::computeHull(std::vector<StuntDouble*>
164      vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1];
165    }
166  
167 <  std::vector<double> globalCoords(dim_ * globalHullSites);
168 <  std::vector<double> globalVels(dim_ * globalHullSites);
169 <  std::vector<double> globalMasses(globalHullSites);
167 >  vector<double> globalCoords(dim_ * globalHullSites);
168 >  vector<double> globalVels(dim_ * globalHullSites);
169 >  vector<double> globalMasses(globalHullSites);
170  
171    int count = coordsOnProc[myrank];
172    
# Line 187 | Line 185 | void AlphaHull::computeHull(std::vector<StuntDouble*>
185    // Free previous hull
186    qh_freeqhull(!qh_ALL);
187    qh_memfreeshort(&curlong, &totlong);
188 <  if (curlong || totlong)
189 <    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
190 <              << totlong << curlong << std::endl;
188 >  if (curlong || totlong) {
189 >    sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n"
190 >            "\tdid not free %d bytes of long memory (%d pieces)",
191 >            totlong, curlong);
192 >    painCave.isFatal = 1;
193 >    simError();
194 >  }
195    
196    if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc,
197                     const_cast<char *>(options_.c_str()), NULL, stderr)){
198      
199 <    sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute global convex hull");
199 >    sprintf(painCave.errMsg,
200 >            "AlphaHull: Qhull failed to compute global convex hull");
201      painCave.isFatal = 1;
202      simError();
203 <    
203 >        
204    } //qh_new_qhull
205  
203 –
206   #endif
207  
208    //Set facet->center as the Voronoi center
# Line 232 | Line 234 | void AlphaHull::computeHull(std::vector<StuntDouble*>
234    
235    qh visit_id++;
236    int numFacets=0;
237 <  std::vector<std::vector <int> > facetlist;
237 >  vector<vector <int> > facetlist;
238    interiorPoint = qh interior_point;
239    FORALLfacet_(qh facet_list) {
240      numFacets++;
# Line 299 | Line 301 | void AlphaHull::computeHull(std::vector<StuntDouble*>
301    
302    //Filter the triangles (only the ones on the boundary of the alpha complex) and build the mesh
303  
304 <
304 >  int ridgesCount=0;
305    
306    ridgeT *ridge, **ridgep;
307    FOREACHridge_(set) {
# Line 314 | Line 316 | void AlphaHull::computeHull(std::vector<StuntDouble*>
316        
317        
318        //coordT *center = qh_getcenter(ridge->vertices);
319 <      //std::cout << "Centers are " << center[0] << "  " <<center[1] << "  " << center[2] << std::endl;
319 >      //cout << "Centers are " << center[0] << "  " <<center[1] << "  " << center[2] << endl;
320        //Vector3d V3dCentroid(center[0], center[1], center[2]);
321        //face.setCentroid(V3dCentroid);
322  
# Line 324 | Line 326 | void AlphaHull::computeHull(std::vector<StuntDouble*>
326        RealType faceMass = 0.0;
327        
328        int ver = 0;
329 <      std::vector<int> virtexlist;
329 >      vector<int> virtexlist;
330        FOREACHvertex_i_(ridge->vertices){
331          int id = qh_pointid(vertex->point);
332          p[ver][0] = vertex->point[0];
# Line 334 | Line 336 | void AlphaHull::computeHull(std::vector<StuntDouble*>
336          RealType mass;
337          ver++;
338          virtexlist.push_back(id);
339 <        // std::cout << "Ridge: " << ridgesCount << " Vertex " << id << std::endl;
339 >        // cout << "Ridge: " << ridgesCount << " Vertex " << id << endl;
340  
341          vel = bodydoubles[id]->getVel();
342          mass = bodydoubles[id]->getMass();
# Line 354 | Line 356 | void AlphaHull::computeHull(std::vector<StuntDouble*>
356        Vector3d normal = face.getUnitNormal();
357        RealType offset =  ((0.0-p[0][0])*normal[0] + (0.0-p[0][1])*normal[1] + (0.0-p[0][2])*normal[2]);
358        RealType dist =  normal[0] * interiorPoint[0] + normal[1]*interiorPoint[1] + normal[2]*interiorPoint[2];
359 <      std::cout << "Dist and normal and area are: " << normal << std::endl;
359 >     cout << "Dist and normal and area are: " << normal << endl;
360        volume_ += dist *area/qh hull_dim;
361        
362        Triangles_.push_back(face);
363      }
364    }
365  
366 <  std::cout << "Volume is: " << volume_ << std::endl;
366 >  cout << "Volume is: " << volume_ << endl;
367  
368   //assert(pm.cm.fn == ridgesCount);
369   /*
# Line 457 | Line 459 | void AlphaHull::computeHull(std::vector<StuntDouble*>
459  
460    qh_freeqhull(!qh_ALL);
461    qh_memfreeshort(&curlong, &totlong);
462 <  if (curlong || totlong)
463 <    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
464 <              << totlong << curlong << std::endl;    
462 >  if (curlong || totlong) {
463 >    sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n"
464 >            "\tdid not free %d bytes of long memory (%d pieces)",
465 >            totlong, curlong);
466 >    painCave.isFatal = 1;
467 >    simError();
468 >  }
469   }
470  
471 < void AlphaHull::printHull(const std::string& geomFileName) {
472 <
471 > void AlphaHull::printHull(const string& geomFileName) {
472 >  
473   #ifdef IS_MPI
474    if (worldRank == 0)  {
475   #endif
476 <  FILE *newGeomFile;
477 <  
478 <  //create new .md file based on old .md file
479 <  newGeomFile = fopen(geomFileName.c_str(), "w");
480 <  qh_findgood_all(qh facet_list);
481 <  for (int i = 0; i < qh_PRINTEND; i++)
482 <    qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
483 <  
484 <  fclose(newGeomFile);
476 >    FILE *newGeomFile;
477 >    
478 >    //create new .md file based on old .md file
479 >    newGeomFile = fopen(geomFileName.c_str(), "w");
480 >    qh_findgood_all(qh facet_list);
481 >    for (int i = 0; i < qh_PRINTEND; i++)
482 >      qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
483 >    
484 >    fclose(newGeomFile);
485   #ifdef IS_MPI
486    }
487   #endif  

Diff Legend

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