--- trunk/src/math/ConvexHull.cpp 2008/06/18 17:03:30 1261 +++ branches/development/src/math/ConvexHull.cpp 2013/01/09 19:27:52 1825 @@ -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,6 +28,15 @@ * 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 * @@ -44,11 +44,12 @@ * * Created by Charles F. Vardeman II on 11 Dec 2006. * @author Charles F. Vardeman II - * @version $Id: ConvexHull.cpp,v 1.7 2008-06-18 17:03:30 chuckv Exp $ + * @version $Id$ * */ /* Standard includes independent of library */ + #include #include #include @@ -56,250 +57,250 @@ #include #include "math/ConvexHull.hpp" #include "utils/simError.h" -using namespace oopse; -/* CGAL version of convex hull first then QHULL */ -#ifdef HAVE_CGAL +#ifdef IS_MPI +#include +#endif -#include -#include -#include -#include -#include +#include "math/qhull.hpp" -#include -#include -#include +#ifdef HAVE_QHULL +using namespace OpenMD; -typedef double RT; -typedef CGAL::Simple_cartesian K; -typedef CGAL::Convex_hull_traits_3 Traits; -typedef Traits::Polyhedron_3 Polyhedron_3; -typedef K::Point_3 Point_3; +ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp") { +} +void ConvexHull::computeHull(std::vector bodydoubles) { + + int numpoints = bodydoubles.size(); -ConvexHull::ConvexHull(){} - -bool ConvexHull::genHull(std::vector pos) -{ + Triangles_.clear(); - std::vector points; + vertexT *vertex, **vertexp; + facetT *facet; + setT *vertices; + int curlong, totlong; + // pointT *intPoint; + std::vector ptArray(numpoints*dim_); + + // 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++; + } - // 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); - } + /* Clean up memory from previous convex hull calculations */ + boolT ismalloc = False; - // define object to hold convex hull - CGAL::Object ch_object_; - Polyhedron_3 polyhedron; + /* 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 - std::cerr << "Creating hull" << std::endl; - CGAL::convex_hull_3(points.begin(), points.end(), ch_object_); - std::cerr << "Done Creating hull" << std::endl; - std::vector::const_iterator p_it; - for (p_it = points.begin(); p_it != points.end(); p_it++) - { - std::cerr << (*p_it).x() << std::endl; - } +#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; - /* - for (Polyhedron_3::Vertex_iterator v = ch_object_.vertices_begin(); - ch_object_.vertices_end(); ++v){ - std::cout<< v.point()< hullSitesOnProc(nproc, 0); + std::vector coordsOnProc(nproc, 0); + std::vector displacements(nproc, 0); + std::vector vectorDisplacements(nproc, 0); + std::vector coords; + std::vector vels; + std::vector indexMap; + std::vector masses; + FORALLvertices{ + localHullSites++; + + int idx = qh_pointid(vertex->point); -#else -#ifdef HAVE_QHULL -/* Old options Qt Qu Qg QG0 FA */ -ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt Qci Tcv Pp") - //ConvexHull::ConvexHull() : dim_(3), options_("qhull d Qbb Qt i") -{} + indexMap.push_back(idx); -bool ConvexHull::genHull(std::vector pos) -{ + coords.push_back(ptArray[dim_ * idx]); + coords.push_back(ptArray[dim_ * idx + 1]); + coords.push_back(ptArray[dim_ * idx + 2]); + + StuntDouble* sd = bodydoubles[idx]; + + Vector3d vel = sd->getVel(); + vels.push_back(vel.x()); + vels.push_back(vel.y()); + vels.push_back(vel.z()); + + masses.push_back(sd->getMass()); + } + + MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], + 1, MPI::INT); + + int globalHullSites = 0; + for (int iproc = 0; iproc < nproc; iproc++){ + globalHullSites += hullSitesOnProc[iproc]; + coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc]; + } + + displacements[0] = 0; + vectorDisplacements[0] = 0; - - int numpoints = pos.size(); - coordT* pt_array; - coordT* surfpt_array; - std::list surface_atoms; - FILE *outdummy = NULL; - FILE *errdummy = NULL; - - pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_)); - - - for (int i = 0; i < numpoints; i++) { - pt_array[dim_ * i] = pos[i][0]; - pt_array[dim_ * i + 1] = pos[i][1]; - pt_array[dim_ * i + 2] = pos[i][2]; + 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]; - - - - /* - qh_initflags(const_cast(options_.c_str())); - qh_init_B(pospoints, numpoints, dim_, ismalloc); - qh_qhull(); - qh_check_output(); + MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0], + &coordsOnProc[0], &vectorDisplacements[0], + MPI::DOUBLE); - qh_produce_output(); - */ - boolT ismalloc = False; + 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_, numpoints, pt_array, ismalloc, - const_cast(options_.c_str()), NULL, stderr)) { + if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc, + const_cast(options_.c_str()), NULL, stderr)){ - vertexT *vertex, **vertexp; - facetT *facet; - setT *vertices; - unsigned int nf = qh num_facets; + sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute global convex hull"); + painCave.isFatal = 1; + simError(); - //Matrix idx(nf, dim); - /* - int j, i = 0, id = 0; + } //qh_new_qhull + +#endif + // commented out below, so comment out here also. + // 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); - int id2 = 0; - coordT *point,**pointp; - realT dist; - FORALLfacets { - j = 0; - - if (!facet->simplicial){ - // should never happen with Qt - sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected"); - painCave.isFatal = 0; - simError(); - } - - vertices = qh_facet3vertex(facet); - FOREACHvertex_(vertices){ - id = qh_pointid(vertex->point); - surface_atoms.push_back(id); - //std::cout << "Ag " << pos[id][0] << " " << pos[id][1] << " " << pos[id][2]<< std::endl; - } - qh_settempfree(&vertices); - - FOREACHpoint_(facet->coplanarset){ - vertex= qh_nearvertex (facet, point, &dist); - //id= qh_pointid (vertex->point); - id2= qh_pointid (point); - surface_atoms.push_back(id2); - //std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl; - //printf ("%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist); - //std::cout << "Neighbors are: %d $d %d\n" << id << id2 << facet->id; - - } - - } - -*/ - - - } + 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; - qh_getarea(qh facet_list); - volume_ = qh totvol; - area_ = qh totarea; - // FILE *newGeomFile; - - - /* - FORALLfacets { - for (int k=0; k < qh hull_dim; k++) - printf ("%6.2g ", facet->normal[k]); - printf ("\n"); - } - */ - - int curlong,totlong; - // geomviewHull("junk.off"); - 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; - free(pt_array); - /* - int j = 0; - surface_atoms.sort(); - surface_atoms.unique(); - surfpt_array = (coordT*) malloc(sizeof(coordT) * (surface_atoms.size() * dim_)); - for(std::list::iterator list_iter = surface_atoms.begin(); - list_iter != surface_atoms.end(); list_iter++) - { - int i = *list_iter; - //surfpt_array[dim_ * j] = pos[i][0]; - //surfpt_array[dim_ * j + 1] = pos[i][1]; - //surfpt_array[dim_ * j + 2] = pos[i][2]; - std::cout << "Ag " << pos[i][0] << " " << pos[i][1] << " "<< pos[i][2] << std::endl; - j++; - } - */ - - /* - std::string deloptions_ = "qhull d Qt"; - facetT *facet, *neighbor; - ridgeT *ridge, **ridgep; - - if (!qh_new_qhull(dim_, surface_atoms.size(), surfpt_array, ismalloc, - const_cast(deloptions_.c_str()), NULL, stderr)) { - - qh visit_id++; - FORALLfacets { - if (!facet->upperdelaunay) { - facet->visitid= qh visit_id; - qh_makeridges(facet); - FOREACHridge_(facet->ridges) { - neighbor= otherfacet_(ridge, facet); - if (neighbor->visitid != qh visit_id) { - - FOREACHvertex_(ridge->vertices) - int id2 = qh_pointid (vertex->point); - std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl; - } - } - } - +#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]; - 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; - free(surfpt_array); - */ - return true; -} + 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); -RealType ConvexHull::getVolume() -{ - return volume_; + } //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::geomviewHull(const std::string& geomFileName) -{ +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 @@ -309,9 +310,8 @@ void ConvexHull::geomviewHull(const std::string& geomF qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); fclose(newGeomFile); +#ifdef IS_MPI + } +#endif } #endif //QHULL -#endif //CGAL - - -