--- trunk/src/math/ConvexHull.cpp 2006/12/14 19:32:32 1097 +++ trunk/src/math/ConvexHull.cpp 2009/10/20 20:05:28 1374 @@ -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 @@ -40,112 +40,262 @@ * * 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: ConvexHull.cpp,v 1.16 2009-10-20 20:05:28 chuckv Exp $ * */ -#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; +#ifdef HAVE_QHULL +extern "C" +{ +#include +#include +#include +#include +#include +#include +#include +#include +} +/* Old options Qt Qu Qg QG0 FA */ +/* More old opts Qc Qi Pp*/ +ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp") { +} -bool ConvexHull::genHull(std::vector pos) -{ +void ConvexHull::computeHull(std::vector bodydoubles) { + + int numpoints = bodydoubles.size(); - 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); - } + Triangles_.clear(); + + vertexT *vertex, **vertexp; + facetT *facet; + setT *vertices; + int curlong, totlong; + + std::vector ptArray(numpoints*3); + std::vector isSurfaceID(numpoints); - // 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; - } -} + // 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++; + } + + boolT ismalloc = False; + /* Clean up memory from previous convex hull calculations*/ + + if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc, + const_cast(options_.c_str()), NULL, stderr)) { -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)); + sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull"); + painCave.isFatal = 1; + simError(); + + } //qh_new_qhull - Triangulation T(L.begin(), L.end()); - int n = T.number_of_vertices(); +#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; + int* hullSitesOnProc = new int[nproc]; + int* coordsOnProc = new int[nproc]; + int* displacements = new int[nproc]; + int* vectorDisplacements = new int[nproc]; - // 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); + std::vector coords; + std::vector vels; + std::vector objectIDs; + std::vector masses; - n = n + T.insert(V.begin(), V.end()); + 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]); - assert( n == 6 ); // 6 points have been inserted - assert( T.is_valid() ); // checking validity of T + StuntDouble* sd = bodydoubles[idx]; - 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; -} + Vector3d vel = sd->getVel(); + vels.push_back(vel.x()); + vels.push_back(vel.y()); + vels.push_back(vel.z()); -void ConvexHull::geomviewHull(const std::string& geomFileName) -{ + masses.push_back(sd->getMass()); + } - std::ofstream newGeomFile; + MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], + 1, MPI::INT); - //create new .md file based on old .md file - newGeomFile.open(geomFileName.c_str()); + int globalHullSites = 0; + for (int iproc = 0; i < nproc; iproc++){ + globalHullSites += hullSitesOnProc[iproc]; + coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc]; + } - // 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; - } + displacements[0] = 0; + vectorDisplacements[0] = 0; + + for (int iproc = 1; i < nproc; iproc++){ + displacements[iproc] = displacements[iproc-1] + hullSitesOnProc[iproc-1]; + vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1]; + } - newGeomFile.close(); + 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) + std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " + << totlong << curlong << 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; + + 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; + +#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::printHull(const std::string& geomFileName) { + 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); +} +#endif //QHULL