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 1983 by gezelter, Tue Apr 15 20:36:19 2014 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    
115 <  int nproc = MPI::COMM_WORLD.Get_size();
116 <  int myrank = MPI::COMM_WORLD.Get_rank();
115 >  int nproc;
116 >  int myrank;
117 >  MPI_Comm_size( MPI_COMM_WORLD, &nproc);
118 >  MPI_Comm_rank( MPI_COMM_WORLD, &myrank);
119 >
120    int localHullSites = 0;
121  
122 <  std::vector<int> hullSitesOnProc(nproc, 0);
123 <  std::vector<int> coordsOnProc(nproc, 0);
124 <  std::vector<int> displacements(nproc, 0);
125 <  std::vector<int> vectorDisplacements(nproc, 0);
122 >  vector<int> hullSitesOnProc(nproc, 0);
123 >  vector<int> coordsOnProc(nproc, 0);
124 >  vector<int> displacements(nproc, 0);
125 >  vector<int> vectorDisplacements(nproc, 0);
126  
127 <  std::vector<double> coords;
128 <  std::vector<double> vels;
129 <  std::vector<int> indexMap;
130 <  std::vector<double> masses;
127 >  vector<double> coords;
128 >  vector<double> vels;
129 >  vector<int> indexMap;
130 >  vector<double> masses;
131  
132    FORALLvertices{
133      localHullSites++;
134      
135      int idx = qh_pointid(vertex->point);
136 <
136 >    
137      indexMap.push_back(idx);
138  
139      coords.push_back(ptArray[dim_  * idx]);
# Line 149 | Line 150 | void AlphaHull::computeHull(std::vector<StuntDouble*>
150      masses.push_back(sd->getMass());
151    }
152  
153 <  MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0],
154 <                            1, MPI::INT);
153 >  MPI_Allgather(&localHullSites, 1, MPI_INT, &hullSitesOnProc[0],
154 >                1, MPI_INT, MPI_COMM_WORLD);
155  
156    int globalHullSites = 0;
157    for (int iproc = 0; iproc < nproc; iproc++){
# Line 166 | Line 167 | void AlphaHull::computeHull(std::vector<StuntDouble*>
167      vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1];
168    }
169  
170 <  std::vector<double> globalCoords(dim_ * globalHullSites);
171 <  std::vector<double> globalVels(dim_ * globalHullSites);
172 <  std::vector<double> globalMasses(globalHullSites);
170 >  vector<double> globalCoords(dim_ * globalHullSites);
171 >  vector<double> globalVels(dim_ * globalHullSites);
172 >  vector<double> globalMasses(globalHullSites);
173  
174    int count = coordsOnProc[myrank];
175    
176 <  MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0],
177 <                             &coordsOnProc[0], &vectorDisplacements[0],
178 <                             MPI::DOUBLE);
179 <
180 <  MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0],
181 <                             &coordsOnProc[0], &vectorDisplacements[0],
182 <                             MPI::DOUBLE);
183 <
184 <  MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE,
185 <                             &globalMasses[0], &hullSitesOnProc[0],
186 <                             &displacements[0], MPI::DOUBLE);
187 <
176 >  MPI_Allgatherv(&coords[0], count, MPI_DOUBLE, &globalCoords[0],
177 >                 &coordsOnProc[0], &vectorDisplacements[0],
178 >                 MPI_DOUBLE, MPI_COMM_WORLD);
179 >  
180 >  MPI_Allgatherv(&vels[0], count, MPI_DOUBLE, &globalVels[0],
181 >                 &coordsOnProc[0], &vectorDisplacements[0],
182 >                 MPI_DOUBLE, MPI_COMM_WORLD);
183 >  
184 >  MPI_Allgatherv(&masses[0], localHullSites, MPI_DOUBLE,
185 >                 &globalMasses[0], &hullSitesOnProc[0],
186 >                 &displacements[0], MPI_DOUBLE, MPI_COMM_WORLD);
187 >  
188    // Free previous hull
189    qh_freeqhull(!qh_ALL);
190    qh_memfreeshort(&curlong, &totlong);
191 <  if (curlong || totlong)
192 <    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
193 <              << totlong << curlong << std::endl;
191 >  if (curlong || totlong) {
192 >    sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n"
193 >            "\tdid not free %d bytes of long memory (%d pieces)",
194 >            totlong, curlong);
195 >    painCave.isFatal = 1;
196 >    simError();
197 >  }
198    
199    if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc,
200                     const_cast<char *>(options_.c_str()), NULL, stderr)){
201      
202 <    sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute global convex hull");
202 >    sprintf(painCave.errMsg,
203 >            "AlphaHull: Qhull failed to compute global convex hull");
204      painCave.isFatal = 1;
205      simError();
206 <    
206 >        
207    } //qh_new_qhull
202 –
208  
209   #endif
210  
# Line 232 | Line 237 | void AlphaHull::computeHull(std::vector<StuntDouble*>
237    
238    qh visit_id++;
239    int numFacets=0;
240 <  std::vector<std::vector <int> > facetlist;
240 >  vector<vector <int> > facetlist;
241    interiorPoint = qh interior_point;
242    FORALLfacet_(qh facet_list) {
243      numFacets++;
# Line 299 | Line 304 | void AlphaHull::computeHull(std::vector<StuntDouble*>
304    
305    //Filter the triangles (only the ones on the boundary of the alpha complex) and build the mesh
306  
307 <
307 >  int ridgesCount=0;
308    
309    ridgeT *ridge, **ridgep;
310    FOREACHridge_(set) {
# Line 314 | Line 319 | void AlphaHull::computeHull(std::vector<StuntDouble*>
319        
320        
321        //coordT *center = qh_getcenter(ridge->vertices);
322 <      //std::cout << "Centers are " << center[0] << "  " <<center[1] << "  " << center[2] << std::endl;
322 >      //cout << "Centers are " << center[0] << "  " <<center[1] << "  " << center[2] << endl;
323        //Vector3d V3dCentroid(center[0], center[1], center[2]);
324        //face.setCentroid(V3dCentroid);
325  
# Line 324 | Line 329 | void AlphaHull::computeHull(std::vector<StuntDouble*>
329        RealType faceMass = 0.0;
330        
331        int ver = 0;
332 <      std::vector<int> virtexlist;
332 >      vector<int> virtexlist;
333        FOREACHvertex_i_(ridge->vertices){
334          int id = qh_pointid(vertex->point);
335          p[ver][0] = vertex->point[0];
# Line 334 | Line 339 | void AlphaHull::computeHull(std::vector<StuntDouble*>
339          RealType mass;
340          ver++;
341          virtexlist.push_back(id);
342 <        // std::cout << "Ridge: " << ridgesCount << " Vertex " << id << std::endl;
342 >        // cout << "Ridge: " << ridgesCount << " Vertex " << id << endl;
343  
344          vel = bodydoubles[id]->getVel();
345          mass = bodydoubles[id]->getMass();
# Line 354 | Line 359 | void AlphaHull::computeHull(std::vector<StuntDouble*>
359        Vector3d normal = face.getUnitNormal();
360        RealType offset =  ((0.0-p[0][0])*normal[0] + (0.0-p[0][1])*normal[1] + (0.0-p[0][2])*normal[2]);
361        RealType dist =  normal[0] * interiorPoint[0] + normal[1]*interiorPoint[1] + normal[2]*interiorPoint[2];
362 <      std::cout << "Dist and normal and area are: " << normal << std::endl;
362 >     cout << "Dist and normal and area are: " << normal << endl;
363        volume_ += dist *area/qh hull_dim;
364        
365        Triangles_.push_back(face);
366      }
367    }
368  
369 <  std::cout << "Volume is: " << volume_ << std::endl;
369 >  cout << "Volume is: " << volume_ << endl;
370  
371   //assert(pm.cm.fn == ridgesCount);
372   /*
# Line 457 | Line 462 | void AlphaHull::computeHull(std::vector<StuntDouble*>
462  
463    qh_freeqhull(!qh_ALL);
464    qh_memfreeshort(&curlong, &totlong);
465 <  if (curlong || totlong)
466 <    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
467 <              << totlong << curlong << std::endl;    
465 >  if (curlong || totlong) {
466 >    sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n"
467 >            "\tdid not free %d bytes of long memory (%d pieces)",
468 >            totlong, curlong);
469 >    painCave.isFatal = 1;
470 >    simError();
471 >  }
472   }
473  
474 < void AlphaHull::printHull(const std::string& geomFileName) {
475 <
474 > void AlphaHull::printHull(const string& geomFileName) {
475 >  
476   #ifdef IS_MPI
477    if (worldRank == 0)  {
478   #endif
479 <  FILE *newGeomFile;
480 <  
481 <  //create new .md file based on old .md file
482 <  newGeomFile = fopen(geomFileName.c_str(), "w");
483 <  qh_findgood_all(qh facet_list);
484 <  for (int i = 0; i < qh_PRINTEND; i++)
485 <    qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
486 <  
487 <  fclose(newGeomFile);
479 >    FILE *newGeomFile;
480 >    
481 >    //create new .md file based on old .md file
482 >    newGeomFile = fopen(geomFileName.c_str(), "w");
483 >    qh_findgood_all(qh facet_list);
484 >    for (int i = 0; i < qh_PRINTEND; i++)
485 >      qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
486 >    
487 >    fclose(newGeomFile);
488   #ifdef IS_MPI
489    }
490   #endif  

Diff Legend

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