--- branches/development/src/math/AlphaHull.cpp 2010/07/09 23:08:25 1465 +++ branches/development/src/math/AlphaHull.cpp 2013/04/03 21:32:13 1858 @@ -1,4 +1,4 @@ -/* Copyright (c) 2008, 2009, 2010 The University of Notre Dame. All Rights Reserved. +/* Copyright (c) 2010 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a * non-exclusive, royalty free, license to use, modify and @@ -32,20 +32,15 @@ * research, please cite the appropriate papers when you publish your * work. Good starting points are: * - * [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). + * [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, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). * - * * AlphaHull.cpp * - * Purpose: To calculate Alpha hull, hull volume libqhull. - * - * Created by Charles F. Vardeman II on 11 Dec 2006. - * @author Charles F. Vardeman II - * @version $Id$ - * + * Purpose: To calculate an alpha-shape hull. */ /* Standard includes independent of library */ @@ -55,7 +50,6 @@ #include #include #include -#include #include "math/AlphaHull.hpp" #include "utils/simError.h" @@ -63,42 +57,36 @@ #include #endif -using namespace OpenMD; +#include "math/qhull.hpp" #ifdef HAVE_QHULL -extern "C" -{ -#include -#include -#include -#include -#include -#include -#include -#include -} -double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim); +using namespace std; +using namespace OpenMD; -AlphaHull::AlphaHull(double alpha) : Hull(), dim_(3), alpha_(alpha), options_("qhull d QJ Tcv Pp") { +double calculate_circumradius(pointT* p0, pointT* p1, pointT* p2, int dim); + +AlphaHull::AlphaHull(double alpha) : Hull(), dim_(3), alpha_(alpha), + options_("qhull d QJ Tcv Pp") { } -void AlphaHull::computeHull(std::vector bodydoubles) { +void AlphaHull::computeHull(vector bodydoubles) { int numpoints = bodydoubles.size(); - bool alphashape=true; Triangles_.clear(); - vertexT *vertex, **vertexp; + vertexT *vertex; facetT *facet, *neighbor; - setT *vertices, *verticestop, *verticesbottom; - int curlong, totlong; pointT *interiorPoint; + int curlong, totlong; - std::vector ptArray(numpoints*dim_); - + Vector3d boxMax; + Vector3d boxMin; + + vector ptArray(numpoints*dim_); + // Copy the positon vector into a points vector for qhull. - std::vector::iterator SD; + vector::iterator SD; int i = 0; for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){ Vector3d pos = (*SD)->getPos(); @@ -107,21 +95,22 @@ void AlphaHull::computeHull(std::vector ptArray[dim_ * i + 2] = pos.z(); i++; } - - /* Clean up memory from previous convex hull calculations*/ + + /* Clean up memory from previous convex hull calculations*/ boolT ismalloc = False; - - int ridgesCount=0; + + /* compute the hull for our local points (or all the points for single + processor versions) */ if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc, const_cast(options_.c_str()), NULL, stderr)) { - + sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute convex hull"); painCave.isFatal = 1; simError(); } //qh_new_qhull - - + + #ifdef IS_MPI //If we are doing the mpi version, set up some vectors for data communication @@ -129,15 +118,15 @@ void AlphaHull::computeHull(std::vector int myrank = MPI::COMM_WORLD.Get_rank(); int localHullSites = 0; - std::vector hullSitesOnProc(nproc, 0); - std::vector coordsOnProc(nproc, 0); - std::vector displacements(nproc, 0); - std::vector vectorDisplacements(nproc, 0); + vector hullSitesOnProc(nproc, 0); + vector coordsOnProc(nproc, 0); + vector displacements(nproc, 0); + vector vectorDisplacements(nproc, 0); - std::vector coords; - std::vector vels; - std::vector indexMap; - std::vector masses; + vector coords; + vector vels; + vector indexMap; + vector masses; FORALLvertices{ localHullSites++; @@ -177,9 +166,9 @@ void AlphaHull::computeHull(std::vector vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1]; } - std::vector globalCoords(dim_ * globalHullSites); - std::vector globalVels(dim_ * globalHullSites); - std::vector globalMasses(globalHullSites); + vector globalCoords(dim_ * globalHullSites); + vector globalVels(dim_ * globalHullSites); + vector globalMasses(globalHullSites); int count = coordsOnProc[myrank]; @@ -198,20 +187,24 @@ void AlphaHull::computeHull(std::vector // Free previous hull qh_freeqhull(!qh_ALL); qh_memfreeshort(&curlong, &totlong); - if (curlong || totlong) - std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " - << totlong << curlong << std::endl; + if (curlong || totlong) { + sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n" + "\tdid not free %d bytes of long memory (%d pieces)", + totlong, curlong); + painCave.isFatal = 1; + simError(); + } if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc, const_cast(options_.c_str()), NULL, stderr)){ - sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute global convex hull"); + sprintf(painCave.errMsg, + "AlphaHull: Qhull failed to compute global convex hull"); painCave.isFatal = 1; simError(); - + } //qh_new_qhull - #endif //Set facet->center as the Voronoi center @@ -243,7 +236,7 @@ void AlphaHull::computeHull(std::vector qh visit_id++; int numFacets=0; - std::vector > facetlist; + vector > facetlist; interiorPoint = qh interior_point; FORALLfacet_(qh facet_list) { numFacets++; @@ -310,7 +303,7 @@ void AlphaHull::computeHull(std::vector //Filter the triangles (only the ones on the boundary of the alpha complex) and build the mesh - + int ridgesCount=0; ridgeT *ridge, **ridgep; FOREACHridge_(set) { @@ -325,7 +318,7 @@ void AlphaHull::computeHull(std::vector //coordT *center = qh_getcenter(ridge->vertices); - //std::cout << "Centers are " << center[0] << " " < RealType faceMass = 0.0; int ver = 0; - std::vector virtexlist; + vector virtexlist; FOREACHvertex_i_(ridge->vertices){ int id = qh_pointid(vertex->point); p[ver][0] = vertex->point[0]; @@ -345,7 +338,7 @@ void AlphaHull::computeHull(std::vector RealType mass; ver++; virtexlist.push_back(id); - // std::cout << "Ridge: " << ridgesCount << " Vertex " << id << std::endl; + // cout << "Ridge: " << ridgesCount << " Vertex " << id << endl; vel = bodydoubles[id]->getVel(); mass = bodydoubles[id]->getMass(); @@ -358,21 +351,21 @@ void AlphaHull::computeHull(std::vector facetlist.push_back(virtexlist); face.addVertices(p[0],p[1],p[2]); face.setFacetMass(faceMass); - face.setFacetVelocity(faceVel/3.0); + face.setFacetVelocity(faceVel / RealType(3.0)); RealType area = face.getArea(); area_ += area; Vector3d normal = face.getUnitNormal(); RealType offset = ((0.0-p[0][0])*normal[0] + (0.0-p[0][1])*normal[1] + (0.0-p[0][2])*normal[2]); RealType dist = normal[0] * interiorPoint[0] + normal[1]*interiorPoint[1] + normal[2]*interiorPoint[2]; - std::cout << "Dist and normal and area are: " << normal << std::endl; + cout << "Dist and normal and area are: " << normal << endl; volume_ += dist *area/qh hull_dim; Triangles_.push_back(face); } } - std::cout << "Volume is: " << volume_ << std::endl; + cout << "Volume is: " << volume_ << endl; //assert(pm.cm.fn == ridgesCount); /* @@ -466,27 +459,51 @@ void AlphaHull::computeHull(std::vector //volume_ = qh totvol; // area_ = qh totarea; + + int index = 0; + FORALLvertices { + Vector3d point(vertex->point[0], vertex->point[1], vertex->point[2]); + if (index == 0) { + boxMax = point; + boxMin = point; + } else { + for (int i = 0; i < 3; i++) { + boxMax[i] = max(boxMax[i], point[i]); + boxMin[i] = min(boxMin[i], point[i]); + } + } + index++; + } + boundingBox_ = Mat3x3d(0.0); + boundingBox_(0,0) = boxMax[0] - boxMin[0]; + boundingBox_(1,1) = boxMax[1] - boxMin[1]; + boundingBox_(2,2) = boxMax[2] - boxMin[2]; + qh_freeqhull(!qh_ALL); qh_memfreeshort(&curlong, &totlong); - if (curlong || totlong) - std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " - << totlong << curlong << std::endl; + if (curlong || totlong) { + sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n" + "\tdid not free %d bytes of long memory (%d pieces)", + totlong, curlong); + painCave.isFatal = 1; + simError(); + } } -void AlphaHull::printHull(const std::string& geomFileName) { - +void AlphaHull::printHull(const string& geomFileName) { + #ifdef IS_MPI if (worldRank == 0) { #endif - FILE *newGeomFile; - - //create new .md file based on old .md file - newGeomFile = fopen(geomFileName.c_str(), "w"); - qh_findgood_all(qh facet_list); - for (int i = 0; i < qh_PRINTEND; i++) - qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); - - fclose(newGeomFile); + FILE *newGeomFile; + + //create new .md file based on old .md file + newGeomFile = fopen(geomFileName.c_str(), "w"); + qh_findgood_all(qh facet_list); + for (int i = 0; i < qh_PRINTEND; i++) + qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); + + fclose(newGeomFile); #ifdef IS_MPI } #endif