--- trunk/src/math/ConvexHull.cpp 2008/06/18 17:03:30 1261 +++ trunk/src/math/ConvexHull.cpp 2009/10/20 20:36:56 1376 @@ -1,4 +1,4 @@ -/* Copyright (c) 2006 The University of Notre Dame. All Rights Reserved. +/* Copyright (c) 2008, 2009 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 @@ -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: ConvexHull.cpp,v 1.17 2009-10-20 20:36:56 gezelter Exp $ * */ /* Standard includes independent of library */ + #include #include #include @@ -56,250 +57,240 @@ #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 +using namespace oopse; -#include -#include -#include +#ifdef HAVE_QHULL +extern "C" +{ +#include +#include +#include +#include +#include +#include +#include +#include +} -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; +/* Old options Qt Qu Qg QG0 FA */ +/* More old opts Qc Qi Pp*/ +ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp") { +} -ConvexHull::ConvexHull(){} +void ConvexHull::computeHull(std::vector bodydoubles) { + + int numpoints = bodydoubles.size(); -bool ConvexHull::genHull(std::vector pos) -{ + Triangles_.clear(); - std::vector points; + vertexT *vertex, **vertexp; + facetT *facet; + setT *vertices; + int curlong, totlong; + std::vector ptArray(numpoints*3); + std::vector isSurfaceID(numpoints); + + // 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); - } + boolT ismalloc = False; + /* Clean up memory from previous convex hull calculations*/ - // define object to hold convex hull - CGAL::Object ch_object_; - Polyhedron_3 polyhedron; + if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc, + const_cast(options_.c_str()), NULL, stderr)) { - // 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; + sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull"); + painCave.isFatal = 1; + simError(); + + } //qh_new_qhull - for (p_it = points.begin(); p_it != points.end(); p_it++) - { - std::cerr << (*p_it).x() << std::endl; - } - /* - for (Polyhedron_3::Vertex_iterator v = ch_object_.vertices_begin(); - ch_object_.vertices_end(); ++v){ - std::cout<< v.point()< coords; + std::vector vels; + std::vector objectIDs; + std::vector masses; + + FORALLvertices{ + localHullSites++; + + int idx = qh_pointid(vertex->point); + 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()); } - */ -} -#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") -{} + MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], + 1, MPI::INT); -bool ConvexHull::genHull(std::vector pos) -{ - - - 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]; + int globalHullSites = 0; + for (int iproc = 0; iproc < nproc; iproc++){ + std::cerr << "iproc = " << iproc << " sites = " << hullSitesOnProc[iproc] << "\n"; + globalHullSites += hullSitesOnProc[iproc]; + coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc]; } - - - - - /* - qh_initflags(const_cast(options_.c_str())); - qh_init_B(pospoints, numpoints, dim_, ismalloc); - qh_qhull(); - qh_check_output(); - qh_produce_output(); - */ - boolT ismalloc = False; + displacements[0] = 0; + vectorDisplacements[0] = 0; - if (!qh_new_qhull(dim_, numpoints, pt_array, ismalloc, - const_cast(options_.c_str()), NULL, stderr)) { - - vertexT *vertex, **vertexp; - facetT *facet; - setT *vertices; - unsigned int nf = qh num_facets; - - //Matrix idx(nf, dim); - /* - int j, i = 0, id = 0; - - 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; - - } - - } - -*/ - - + 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); - 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"); + // 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; - 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; - } - } - } - + 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 - - } + FORALLfacets { + Triangle face; - 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; -} + 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); + 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; -RealType ConvexHull::getVolume() -{ - return volume_; +#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 (id >= 0 && id < hullSitesOnProc[myrank]) + face.addVertexSD(bodydoubles[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/3.0); + Triangles_.push_back(face); + qh_settempfree(&vertices); + + } //FORALLfacets + + qh_getarea(qh facet_list); + volume_ = qh totvol; + area_ = qh totarea; + +#ifdef IS_MPI + delete [] hullSitesOnProc; + delete [] coordsOnProc; + delete [] displacements; + delete [] vectorDisplacements; +#endif + + 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; } -void ConvexHull::geomviewHull(const std::string& geomFileName) -{ + +void ConvexHull::printHull(const std::string& geomFileName) { FILE *newGeomFile; //create new .md file based on old .md file @@ -311,7 +302,3 @@ void ConvexHull::geomviewHull(const std::string& geomF fclose(newGeomFile); } #endif //QHULL -#endif //CGAL - - -