--- trunk/src/math/AlphaHull.cpp 2010/01/08 17:15:27 1402 +++ trunk/src/math/AlphaHull.cpp 2014/02/26 14:14:50 1969 @@ -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,22 +32,21 @@ * 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: ConvexHull.cpp,v 1.21 2009-11-25 20:02:01 gezelter Exp $ - * + * Purpose: To calculate an alpha-shape hull. */ +#ifdef IS_MPI +#include +#endif + /* Standard includes independent of library */ #include @@ -55,49 +54,37 @@ #include #include #include -#include #include "math/AlphaHull.hpp" #include "utils/simError.h" -#ifdef IS_MPI -#include -#endif +#include "math/qhull.hpp" +#ifdef HAVE_QHULL +using namespace std; using namespace OpenMD; -#ifdef HAVE_QHULL -extern "C" -{ -#include -#include -#include -#include -#include -#include -#include -#include -} -double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim); +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 ") { +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; + pointT *interiorPoint; int curlong, totlong; - std::vector ptArray(numpoints*dim_); - + + 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(); @@ -106,37 +93,41 @@ 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 - int nproc = MPI::COMM_WORLD.Get_size(); - int myrank = MPI::COMM_WORLD.Get_rank(); + int nproc; + int myrank; + MPI_Comm_size( MPI_COMM_WORLD, &nproc); + MPI_Comm_rank( MPI_COMM_WORLD, &myrank); + 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++; @@ -159,8 +150,8 @@ void AlphaHull::computeHull(std::vector masses.push_back(sd->getMass()); } - MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], - 1, MPI::INT); + MPI_Allgather(&localHullSites, 1, MPI_INT, &hullSitesOnProc[0], + 1, MPI_INT, MPI_COMM_WORLD); int globalHullSites = 0; for (int iproc = 0; iproc < nproc; iproc++){ @@ -176,41 +167,45 @@ 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]; - MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0], - &coordsOnProc[0], &vectorDisplacements[0], - MPI::DOUBLE); - - MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0], - &coordsOnProc[0], &vectorDisplacements[0], - MPI::DOUBLE); - - MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE, - &globalMasses[0], &hullSitesOnProc[0], - &displacements[0], MPI::DOUBLE); - + MPI_Allgatherv(&coords[0], count, MPI_DOUBLE, &globalCoords[0], + &coordsOnProc[0], &vectorDisplacements[0], + MPI_DOUBLE, MPI_COMM_WORLD); + + MPI_Allgatherv(&vels[0], count, MPI_DOUBLE, &globalVels[0], + &coordsOnProc[0], &vectorDisplacements[0], + MPI_DOUBLE, MPI_COMM_WORLD); + + MPI_Allgatherv(&masses[0], localHullSites, MPI_DOUBLE, + &globalMasses[0], &hullSitesOnProc[0], + &displacements[0], MPI_DOUBLE, MPI_COMM_WORLD); + // 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 @@ -242,7 +237,8 @@ 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++; if (!facet->upperdelaunay) { @@ -307,6 +303,8 @@ void AlphaHull::computeHull(std::vector //assert(numFacets== qh num_facets); //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) { @@ -319,21 +317,19 @@ void AlphaHull::computeHull(std::vector // Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]); //face.setNormal(V3dNormal); - // RealType faceArea = qh_facetarea(facet); - //face.setArea(faceArea); - //vertices = qh_facet3vertex(facet); - - //coordT *center = qh_getcenter(vertices); + //coordT *center = qh_getcenter(ridge->vertices); + //cout << "Centers are " << center[0] << " " < virtexlist; + vector virtexlist; FOREACHvertex_i_(ridge->vertices){ int id = qh_pointid(vertex->point); p[ver][0] = vertex->point[0]; @@ -343,15 +339,37 @@ void AlphaHull::computeHull(std::vector RealType mass; ver++; virtexlist.push_back(id); - } - facetlist.push_back(virtexlist); + // cout << "Ridge: " << ridgesCount << " Vertex " << id << endl; + vel = bodydoubles[id]->getVel(); + mass = bodydoubles[id]->getMass(); + face.addVertexSD(bodydoubles[id]); + + + faceVel = faceVel + vel; + faceMass = faceMass + mass; + } //FOREACH Vertex + facetlist.push_back(virtexlist); + face.addVertices(p[0],p[1],p[2]); + face.setFacetMass(faceMass); + 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]; + cout << "Dist and normal and area are: " << normal << endl; + volume_ += dist *area/qh hull_dim; + + Triangles_.push_back(face); } } - -//assert(pm.cm.fn == ridgesCount); + cout << "Volume is: " << volume_ << endl; +//assert(pm.cm.fn == ridgesCount); +/* std::cout <<"OFF"< } std::cout << std::endl; } +*/ - /* FORALLfacets { Triangle face; @@ -438,31 +456,35 @@ void AlphaHull::computeHull(std::vector } //FORALLfacets */ - qh_getarea(qh facet_list); - volume_ = qh totvol; - area_ = qh totarea; + // qh_getarea(qh facet_list); + //volume_ = qh totvol; + // area_ = qh totarea; 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