--- trunk/src/math/ConvexHull.cpp 2006/12/14 19:32:32 1097 +++ branches/development/src/math/ConvexHull.cpp 2012/01/06 19:03:05 1668 @@ -1,23 +1,14 @@ -/* Copyright (c) 2006 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 * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,115 +28,298 @@ * University of Notre Dame has been advised of the possibility of * such damages. * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * 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] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). * * * ConvexHull.cpp * - * Purpose: To calculate convexhull, hull volume and radius - * using the CGAL library. + * Purpose: To calculate convexhull, hull volume libqhull. * * Created by Charles F. Vardeman II on 11 Dec 2006. * @author Charles F. Vardeman II - * @version $Id: ConvexHull.cpp,v 1.1 2006-12-14 19:32:32 chuckv Exp $ + * @version $Id$ * */ -#include "math/ConvexHull.hpp" +/* Standard includes independent of library */ + #include #include +#include +#include +#include +#include "math/ConvexHull.hpp" +#include "utils/simError.h" +#ifdef IS_MPI +#include +#endif -using namespace oopse; +using namespace OpenMD; +#ifdef HAVE_QHULL +extern "C" +{ +#include +#include +#include +#include +#include +#include +#include +#include +} +ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp") { +} +void ConvexHull::computeHull(std::vector bodydoubles) { + + int numpoints = bodydoubles.size(); -bool ConvexHull::genHull(std::vector pos) -{ + Triangles_.clear(); + + vertexT *vertex, **vertexp; + facetT *facet; + setT *vertices; + int curlong, totlong; + pointT *intPoint; + + std::vector ptArray(numpoints*dim_); - std::vector points; - points.reserve(pos.size()); - // Copy the positon vector into a points vector for cgal. - for (int i = 0; i < pos.size(); ++i) - { - Point_3 pt(pos[i][0],pos[i][1],pos[i][2]); - points.push_back(pt); - } + // Copy the positon vector into a points vector for qhull. + std::vector::iterator SD; + int i = 0; + for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){ + Vector3d pos = (*SD)->getPos(); + ptArray[dim_ * i] = pos.x(); + ptArray[dim_ * i + 1] = pos.y(); + ptArray[dim_ * i + 2] = pos.z(); + i++; + } + + /* Clean up memory from previous convex hull calculations */ + boolT ismalloc = False; + + /* 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, "ConvexHull: Qhull failed to compute convex hull"); + painCave.isFatal = 1; + simError(); + + } //qh_new_qhull - // compute convex hull - CGAL::convex_hull_3(points.begin(), points.end(), ch_object); - // Make sure hull is a polyhedron - if ( CGAL::assign(ch_polyhedron, ch_object) ) - { - return true; - } - else - { - return false; - } -} -RealType ConvexHull::getVolume() -{ - /* - std::list L; - L.push_front(Point(0,0,0)); - L.push_front(Point(1,0,0)); - L.push_front(Point(0,1,0)); +#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 localHullSites = 0; - Triangulation T(L.begin(), L.end()); + std::vector hullSitesOnProc(nproc, 0); + std::vector coordsOnProc(nproc, 0); + std::vector displacements(nproc, 0); + std::vector vectorDisplacements(nproc, 0); - int n = T.number_of_vertices(); + std::vector coords; + std::vector vels; + std::vector indexMap; + std::vector masses; - // insertion from a vector : - std::vector V(3); - V[0] = Point(0,0,1); - V[1] = Point(1,1,1); - V[2] = Point(2,2,2); + FORALLvertices{ + localHullSites++; + + int idx = qh_pointid(vertex->point); - n = n + T.insert(V.begin(), V.end()); + indexMap.push_back(idx); - assert( n == 6 ); // 6 points have been inserted - assert( T.is_valid() ); // checking validity of T + coords.push_back(ptArray[dim_ * idx]); + coords.push_back(ptArray[dim_ * idx + 1]); + coords.push_back(ptArray[dim_ * idx + 2]); - double sum_v = 0; - for(Triangulation::Cell_iterator c_it = T.cells_begin(); c_it != T.cells_end(); ++c_it) - { - sum_v += T.tetrahedron(c_it).volume(); - } - std::cout << "sum_v " << sum_v << std::endl; -*/ - return 0.0; -} + StuntDouble* sd = bodydoubles[idx]; -void ConvexHull::geomviewHull(const std::string& geomFileName) -{ + Vector3d vel = sd->getVel(); + vels.push_back(vel.x()); + vels.push_back(vel.y()); + vels.push_back(vel.z()); - std::ofstream newGeomFile; + masses.push_back(sd->getMass()); + } - //create new .md file based on old .md file - newGeomFile.open(geomFileName.c_str()); + MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], + 1, MPI::INT); - // Write polyhedron in Object File Format (OFF). - CGAL::set_ascii_mode( std::cout); - newGeomFile << "OFF" << std::endl << ch_polyhedron.size_of_vertices() << ' ' - << ch_polyhedron.size_of_facets() << " 0" << std::endl; - std::copy( ch_polyhedron.points_begin(), ch_polyhedron.points_end(), - std::ostream_iterator( newGeomFile, "\n")); - for ( Facet_iterator i = ch_polyhedron.facets_begin(); i != ch_polyhedron.facets_end(); ++i) - { - Halfedge_facet_circulator j = i->facet_begin(); - // Facets in polyhedral surfaces are at least triangles. - CGAL_assertion( CGAL::circulator_size(j) >= 3); - newGeomFile << CGAL::circulator_size(j) << ' '; - do - { - newGeomFile << ' ' << std::distance(ch_polyhedron.vertices_begin(), j->vertex()); - } - while ( ++j != i->facet_begin()); - newGeomFile << std::endl; - } + int globalHullSites = 0; + for (int iproc = 0; iproc < nproc; iproc++){ + globalHullSites += hullSitesOnProc[iproc]; + coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc]; + } - newGeomFile.close(); + displacements[0] = 0; + vectorDisplacements[0] = 0; + + for (int iproc = 1; iproc < nproc; iproc++){ + displacements[iproc] = displacements[iproc-1] + hullSitesOnProc[iproc-1]; + vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1]; + } + std::vector globalCoords(dim_ * globalHullSites); + std::vector globalVels(dim_ * globalHullSites); + std::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); + + // Free previous hull + qh_freeqhull(!qh_ALL); + qh_memfreeshort(&curlong, &totlong); + if (curlong || totlong) { + sprintf(painCave.errMsg, "ConvexHull: 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, "ConvexHull: Qhull failed to compute global convex hull"); + painCave.isFatal = 1; + simError(); + + } //qh_new_qhull + +#endif + intPoint = qh interior_point; + RealType calcvol = 0.0; + FORALLfacets { + Triangle face; + //Qhull sets the unit normal in facet->normal + Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]); + face.setUnitNormal(V3dNormal); + + RealType faceArea = qh_facetarea(facet); + face.setArea(faceArea); + + vertices = qh_facet3vertex(facet); + + coordT *center = qh_getcenter(vertices); + Vector3d V3dCentroid(center[0], center[1], center[2]); + face.setCentroid(V3dCentroid); + + Vector3d faceVel = V3Zero; + Vector3d p[3]; + RealType faceMass = 0.0; + + int ver = 0; + + FOREACHvertex_(vertices){ + int id = qh_pointid(vertex->point); + p[ver][0] = vertex->point[0]; + p[ver][1] = vertex->point[1]; + p[ver][2] = vertex->point[2]; + Vector3d vel; + RealType mass; + +#ifdef IS_MPI + vel = Vector3d(globalVels[dim_ * id], + globalVels[dim_ * id + 1], + globalVels[dim_ * id + 2]); + mass = globalMasses[id]; + + // localID will be between 0 and hullSitesOnProc[myrank] if we + // own this guy. + + int localID = id - displacements[myrank]; + + + if (localID >= 0 && localID < hullSitesOnProc[myrank]){ + face.addVertexSD(bodydoubles[indexMap[localID]]); + }else{ + face.addVertexSD(NULL); + } +#else + vel = bodydoubles[id]->getVel(); + mass = bodydoubles[id]->getMass(); + face.addVertexSD(bodydoubles[id]); +#endif + faceVel = faceVel + vel; + faceMass = faceMass + mass; + ver++; + } //Foreachvertex + + face.addVertices(p[0], p[1], p[2]); + face.setFacetMass(faceMass); + face.setFacetVelocity(faceVel / RealType(3.0)); + /* + RealType comparea = face.computeArea(); + realT calcarea = qh_facetarea (facet); + Vector3d V3dCompNorm = -face.computeUnitNormal(); + RealType thisOffset = ((0.0-p[0][0])*V3dCompNorm[0] + (0.0-p[0][1])*V3dCompNorm[1] + (0.0-p[0][2])*V3dCompNorm[2]); + RealType dist = facet->offset + intPoint[0]*V3dNormal[0] + intPoint[1]*V3dNormal[1] + intPoint[2]*V3dNormal[2]; + std::cout << "facet offset and computed offset: " << facet->offset << " " << thisOffset << std::endl; + calcvol += -dist*comparea/qh hull_dim; + */ + Triangles_.push_back(face); + qh_settempfree(&vertices); + + } //FORALLfacets + + qh_getarea(qh facet_list); + volume_ = qh totvol; + area_ = qh totarea; + // std::cout << "My volume is: " << calcvol << " qhull volume is:" << volume_ << std::endl; + qh_freeqhull(!qh_ALL); + qh_memfreeshort(&curlong, &totlong); + if (curlong || totlong) { + sprintf(painCave.errMsg, "ConvexHull: qhull internal warning:\n" + "\tdid not free %d bytes of long memory (%d pieces)", + totlong, curlong); + painCave.isFatal = 1; + simError(); + } } + +void ConvexHull::printHull(const std::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); +#ifdef IS_MPI + } +#endif +} +#endif //QHULL