--- trunk/src/math/ConvexHull.cpp 2008/10/15 18:26:01 1304 +++ trunk/src/math/ConvexHull.cpp 2009/10/20 20:36:56 1376 @@ -1,4 +1,4 @@ -/* Copyright (c) 2008 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.10 2008-10-15 18:26:01 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 @@ -57,638 +58,239 @@ #include "math/ConvexHull.hpp" #include "utils/simError.h" +#ifdef IS_MPI +#include +#endif using namespace oopse; -/* CGAL version of convex hull first then QHULL */ -#ifdef HAVE_CGAL -//#include -#include -//#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#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*/ -//#include -#include -//#include +ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp") { +} - - -typedef CGAL::MP_Float RT; -//typedef double RT; -//typedef CGAL::Homogeneous K; -typedef CGAL::Exact_predicates_exact_constructions_kernel K; -typedef K::Vector_3 Vector_3; -//typedef CGAL::Convex_hull_traits_3 Traits; -typedef CGAL::Polyhedron_traits_with_normals_3 Traits; -//typedef Traits::Polyhedron_3 Polyhedron_3; -typedef CGAL::Polyhedron_3 Polyhedron_3; -typedef K::Point_3 Point_3; - - -typedef Polyhedron_3::HalfedgeDS HalfedgeDS; -typedef Polyhedron_3::Facet_iterator Facet_iterator; -typedef Polyhedron_3::Halfedge_around_facet_circulator Halfedge_facet_circulator; -typedef Polyhedron_3::Halfedge_handle Halfedge_handle; -typedef Polyhedron_3::Facet_iterator Facet_iterator; -typedef Polyhedron_3::Plane_iterator Plane_iterator; -typedef Polyhedron_3::Vertex_iterator Vertex_iterator; -typedef Polyhedron_3::Vertex_handle Vertex_handle; -typedef Polyhedron_3::Point_iterator Point_iterator; +void ConvexHull::computeHull(std::vector bodydoubles) { + int numpoints = bodydoubles.size(); + Triangles_.clear(); + + vertexT *vertex, **vertexp; + facetT *facet; + setT *vertices; + int curlong, totlong; + + std::vector ptArray(numpoints*3); + std::vector isSurfaceID(numpoints); -class Enriched_Point_3 : public K::Point_3{ -public: - Enriched_Point_3(double x,double y,double z) : K::Point_3(x,y,z), yupMyPoint(false), mySD(NULL) {} - - bool isMyPoint() const{ return yupMyPoint; } - void myPoint(){ yupMyPoint = true; } - void setSD(StuntDouble* SD){mySD = SD;} - StuntDouble* getStuntDouble(){return mySD;} -private: - bool yupMyPoint; - StuntDouble* mySD; - -}; - - - - - - // compare Point_3's... used in setting up the STL map from points to indices -template -struct Point_3_comp { - bool operator() (const Pt3 & p, const Pt3 & q) const { - return CGAL::lexicographically_xyz_smaller(p,q); // this is defined inline & hence we had to create fn object & not ptrfun + // 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)) { -// coordinate-based hashing inefficient but can we do better if pts are copied? -typedef std::map > ptMapType; + sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull"); + painCave.isFatal = 1; + simError(); + + } //qh_new_qhull -#ifdef IS_MPI -struct { - double x,y,z; -} surfacePt; -#endif -ConvexHull::ConvexHull() : Hull(){ - //If we are doing the mpi version, set up some vectors for data communication #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]; + std::vector coords; + std::vector vels; + std::vector objectIDs; + std::vector masses; - nproc_ = MPI::COMM_WORLD.Get_size(); - myrank_ = MPI::COMM_WORLD.Get_rank(); - NstoProc_ = new int[nproc_]; - displs_ = new int[nproc_]; + 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]); - // Create a surface point type in MPI to send - surfacePtType = MPI::DOUBLE.Create_contiguous(3); - surfacePtType.Commit(); - + StuntDouble* sd = bodydoubles[idx]; -#endif -} + Vector3d vel = sd->getVel(); + vels.push_back(vel.x()); + vels.push_back(vel.y()); + vels.push_back(vel.z()); -void ConvexHull::computeHull(std::vector bodydoubles) -{ - - std::vector points; - ptMapType myMap; - Point_iterator hc; - - // Copy the positon vector into a points vector for cgal. - std::vector::iterator SD; + masses.push_back(sd->getMass()); + } - for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD) - { - Vector3d pos = (*SD)->getPos(); - Enriched_Point_3* pt = new Enriched_Point_3(pos.x(),pos.y(),pos.z()); - pt->setSD(*SD); - points.push_back(*pt); - // myMap[pt]=(*SD); - } - - // define object to hold convex hull - CGAL::Object ch_object_; - Polyhedron_3 polyhedron; - // compute convex hull - - std::vector::iterator testpt; - - - - CGAL::convex_hull_3(points.begin(), points.end(), polyhedron); - - - Ns_ = polyhedron.size_of_vertices(); + MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], + 1, MPI::INT); -#ifdef IS_MPI - /* Gather an array of the number of verticies on each processor */ - - - surfacePtsGlobal_.clear(); - surfacePtsLocal_.clear(); - - MPI::COMM_WORLD.Allgather(&Ns_,1,MPI::INT,&NstoProc_[0],1,MPI::INT); - - for (int i = 0; i < nproc_; i++){ - Nsglobal_ += NstoProc_[i]; + 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]; } - /*Reminder ideally, we would like to reserve size for the vectors here*/ - surfacePtsLocal_.reserve(Ns_); - surfacePtsGlobal_.resize(Nsglobal_); - // std::fill(surfacePtsGlobal_.begin(),surfacePtsGlobal_.end(),0); - /* Build a displacements array */ - for (int i = 1; i < nproc_; i++){ - displs_[i] = displs_[i-1] + NstoProc_[i-1]; - } + displacements[0] = 0; + vectorDisplacements[0] = 0; - int noffset = displs_[myrank_]; - /* gather the potential hull */ - - - for (hc =polyhedron.points_begin();hc != polyhedron.points_end(); ++hc){ - Point_3 mypoint = *hc; - surfacePt_ mpiSurfacePt; - mpiSurfacePt.x = CGAL::to_double(mypoint.x()); - mpiSurfacePt.y = CGAL::to_double(mypoint.y()); - mpiSurfacePt.z = CGAL::to_double(mypoint.z()); - surfacePtsLocal_.push_back(mpiSurfacePt); + for (int iproc = 1; iproc < nproc; iproc++){ + displacements[iproc] = displacements[iproc-1] + hullSitesOnProc[iproc-1]; + vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1]; } - MPI::COMM_WORLD.Allgatherv(&surfacePtsLocal_[0],Ns_,surfacePtType,&surfacePtsGlobal_[0],NstoProc_,displs_,surfacePtType); - std::vector::iterator spt; - std::vector gblpoints; - - int mine = 0; - int pointidx = 0; - for (spt = surfacePtsGlobal_.begin(); spt != surfacePtsGlobal_.end(); ++spt) - { - surfacePt_ thispos = *spt; - Enriched_Point_3 ept(thispos.x,thispos.y,thispos.z); - if (mine >= noffset && mine < noffset + Ns_){ - ept.myPoint(); - ept.setSD(points[pointidx].getStuntDouble()); - pointidx++; - } - gblpoints.push_back(ept); - - mine++; - } - - /* Compute the global hull */ - polyhedron.clear(); - CGAL::convex_hull_3(gblpoints.begin(), gblpoints.end(), polyhedron); - - -#endif - - + std::vector globalCoords(dim_*globalHullSites); + std::vector globalVels(dim_*globalHullSites); + std::vector globalMasses(globalHullSites); + int count = coordsOnProc[myrank]; - /* Loop over all of the surface triangles and build data structures for atoms and normals*/ - Facet_iterator j; - area_ = 0; - for ( j = polyhedron.facets_begin(); j !=polyhedron.facets_end(); ++j) { - Halfedge_handle h = j->halfedge(); + MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, + &globalCoords[0], &coordsOnProc[0], &vectorDisplacements[0], + MPI::DOUBLE); - Point_3 r0=h->vertex()->point(); - Point_3 r1=h->next()->vertex()->point(); - Point_3 r2=h->next()->next()->vertex()->point(); + MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, + &globalVels[0], &coordsOnProc[0], &vectorDisplacements[0], + MPI::DOUBLE); - Point_3* pr0 = &r0; - Point_3* pr1 = &r1; - Point_3* pr2 = &r2; + MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE, + &globalMasses[0], &hullSitesOnProc[0], &displacements[0], + MPI::DOUBLE); - Enriched_Point_3* er0 = static_cast(pr0); - Enriched_Point_3* er1 = static_cast(pr1); - Enriched_Point_3* er2 = static_cast(pr2); - - // StuntDouble* sd = er0->getStuntDouble(); - std::cerr << "sd globalIndex = " << to_double(er0->x()) << "\n"; - - Point_3 thisCentroid = CGAL::centroid(r0,r1,r2); - - Vector_3 normal = CGAL::cross_product(r1-r0,r2-r0); - - Triangle* face = new Triangle(); - Vector3d V3dNormal(CGAL::to_double(normal.x()),CGAL::to_double(normal.y()),CGAL::to_double(normal.z())); - Vector3d V3dCentroid(CGAL::to_double(thisCentroid.x()),CGAL::to_double(thisCentroid.y()),CGAL::to_double(thisCentroid.z())); - face->setNormal(V3dNormal); - face->setCentroid(V3dCentroid); - RealType faceArea = 0.5*V3dNormal.length(); - face->setArea(faceArea); - area_ += faceArea; - Triangles_.push_back(face); - // ptMapType::const_iterator locn=myMap.find(mypoint); - // int myIndex = locn->second; - - } + // 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; - std::cout << "Number of surface atoms is: " << Ns_ << std::endl; - - - -} -void ConvexHull::printHull(const std::string& geomFileName) -{ - /* - std::ofstream newGeomFile; - - //create new .md file based on old .md file - newGeomFile.open("testhull.off"); - - // Write polyhedron in Object File Format (OFF). - CGAL::set_ascii_mode( std::cout); - newGeomFile << "OFF" << std::endl << polyhedron.size_of_vertices() << ' ' - << polyhedron.size_of_facets() << " 0" << std::endl; - std::copy( polyhedron.points_begin(), polyhedron.points_end(), - std::ostream_iterator( newGeomFile, "\n")); - for ( Facet_iterator i = polyhedron.facets_begin(); i != 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(polyhedron.vertices_begin(), j->vertex()); - } while ( ++j != i->facet_begin()); - newGeomFile << std::endl; - } - - newGeomFile.close(); - */ -/* - std::ofstream newGeomFile; - - //create new .md file based on old .md file - newGeomFile.open(geomFileName.c_str()); - - // 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; - } - - newGeomFile.close(); -*/ - -} - + 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 - - - - - -#else -#ifdef HAVE_QHULL -/* Old options Qt Qu Qg QG0 FA */ -/* More old opts Qc Qi Pp*/ -ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp"), Ns_(200) { - //If we are doing the mpi version, set up some vectors for data communication -#ifdef IS_MPI - - - nproc_ = MPI::COMM_WORLD.Get_size(); - myrank_ = MPI::COMM_WORLD.Get_rank(); - NstoProc_ = new int[nproc_]; - displs_ = new int[nproc_]; - - // Create a surface point type in MPI to send - //surfacePtType = MPI::DOUBLE.Create_contiguous(3); - // surfacePtType.Commit(); - - #endif -} + FORALLfacets { + Triangle face; - -void ConvexHull::computeHull(std::vector bodydoubles) -{ - - std::vector surfaceIDs; - std::vector surfaceIDsGlobal; - std::vector localPtsMap; - int numpoints = bodydoubles.size(); - - //coordT* pt_array; - coordT* surfpt_array; - vertexT *vertex, **vertexp; - facetT *facet; - setT *vertices; - int curlong,totlong; - int id; - - coordT *point,**pointp; - - - FILE *outdummy = NULL; - FILE *errdummy = NULL; - - //pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_)); - -// double* ptArray = new double[numpoints * 3]; - 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(); + 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); - ptArray[dim_ * i] = pos.x(); - ptArray[dim_ * i + 1] = pos.y(); - ptArray[dim_ * i + 2] = pos.z(); - i++; - } - + coordT *center = qh_getcenter(vertices); + Vector3d V3dCentroid(center[0], center[1], center[2]); + face.setCentroid(V3dCentroid); - - - - - boolT ismalloc = False; - /* Clean up memory from previous convex hull calculations*/ - Triangles_.clear(); - surfaceSDs_.clear(); - surfaceSDs_.reserve(Ns_); + Vector3d faceVel = V3Zero; + Vector3d p[3]; + RealType faceMass = 0.0; + int ver = 0; - 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 = 0; - simError(); + 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]; - } //qh_new_qhull + Vector3d vel; + RealType mass; - #ifdef IS_MPI - std::vector localPts; - int localPtArraySize; - - - std::fill(isSurfaceID.begin(),isSurfaceID.end(),false); - + vel = Vector3d(globalVels[dim_ * id], + globalVels[dim_ * id + 1], + globalVels[dim_ * id + 2]); + mass = globalMasses[id]; - FORALLfacets { - - 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); + // 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 - if( !isSurfaceID[id] ){ - isSurfaceID[id] = true; - } - } + face.addVertices(p[0], p[1], p[2]); + face.setFacetMass(faceMass); + face.setFacetVelocity(faceVel/3.0); + Triangles_.push_back(face); qh_settempfree(&vertices); - - } //FORALLfacets - - - /* - std::sort(surfaceIDs.begin(),surfaceIDs.end()); - surfaceIDs.erase(std::unique(surfaceIDs.begin(), surfaceIDs.end()), surfaceIDs.end()); - int localPtArraySize = surfaceIDs.size() * 3; - */ - - //localPts.resize(localPtArraySize); - //std::fill(localPts.begin(),localPts.end(),0.0); - - - int idx = 0; - int nIsIts = 0; -/* - // Copy the surface points into an array. - for(std::vector::iterator list_iter = isSurfaceID.begin(); - list_iter != isSurfaceID.end(); list_iter++) - { - bool isIt = *list_iter; - if (isIt){ - localPts.push_back(ptArray[dim_ * idx]); - localPts.push_back(ptArray[dim_ * idx + 1]); - localPts.push_back(ptArray[dim_ * idx + 2]); - localPtsMap.push_back(idx); - nIsIts++; - } //Isit - idx++; - } //isSurfaceID - */ - FORALLvertices { - idx = qh_pointid(vertex->point); - localPts.push_back(ptArray[dim_ * idx]); - localPts.push_back(ptArray[dim_ * idx + 1]); - localPts.push_back(ptArray[dim_ * idx + 2]); - localPtsMap.push_back(idx); - } - - - localPtArraySize = localPts.size(); - - - MPI::COMM_WORLD.Allgather(&localPtArraySize,1,MPI::INT,&NstoProc_[0],1,MPI::INT); - - Nsglobal_=0; - for (int i = 0; i < nproc_; i++){ - Nsglobal_ += NstoProc_[i]; - } + } //FORALLfacets - - int nglobalPts = int(Nsglobal_/3); - + qh_getarea(qh facet_list); + volume_ = qh totvol; + area_ = qh totarea; - std::vector globalPts; - globalPts.resize(Nsglobal_); - - isSurfaceID.resize(nglobalPts); - - - std::fill(globalPts.begin(),globalPts.end(),0.0); - - displs_[0] = 0; - /* Build a displacements array */ - for (int i = 1; i < nproc_; i++){ - displs_[i] = displs_[i-1] + NstoProc_[i-1]; - } +#ifdef IS_MPI + delete [] hullSitesOnProc; + delete [] coordsOnProc; + delete [] displacements; + delete [] vectorDisplacements; +#endif - - int noffset = displs_[myrank_]; - /* gather the potential hull */ - - MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE); - - /* - if (myrank_ == 0){ - for (i = 0; i < globalPts.size(); i++){ - std::cout << globalPts[i] << std::endl; - } - } - */ - // 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_, nglobalPts, &globalPts[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 - - - - - - - unsigned int nf = qh num_facets; - - /* Build Surface SD list first */ - - std::fill(isSurfaceID.begin(),isSurfaceID.end(),false); - - FORALLfacets { - - if (!facet->simplicial){ - // should never happen with Qt - sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected"); - painCave.isFatal = 1; - simError(); - } //simplicical - - Triangle* face = new Triangle(); - Vector3d V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]); - face->setNormal(V3dNormal); - - - - RealType faceArea = 0.5*V3dNormal.length(); - face->setArea(faceArea); - - - vertices = qh_facet3vertex(facet); - - coordT *center = qh_getcenter(vertices); - Vector3d V3dCentroid(center[0], center[1], center[2]); - face->setCentroid(V3dCentroid); - - FOREACHvertex_(vertices){ - id = qh_pointid(vertex->point); - int localindex = id; -#ifdef IS_MPI - - if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){ - localindex = localPtsMap[id-noffset/3]; -#endif - face->addVertex(bodydoubles[localindex]); - if( !isSurfaceID[id] ){ - isSurfaceID[id] = true; -#ifdef IS_MPI - -#endif - - surfaceSDs_.push_back(bodydoubles[localindex]); - - } //IF isSurfaceID - -#ifdef IS_MPI - - }else{ - face->addVertex(NULL); - } -#endif - } //Foreachvertex - /* - if (!SETempty_(facet->coplanarset)){ - FOREACHpoint_(facet->coplanarset){ - id = qh_pointid(point); - surfaceSDs_.push_back(bodydoubles[id]); - } - } - - Triangles_.push_back(face); - qh_settempfree(&vertices); - */ - } //FORALLfacets - - /* - std::cout << surfaceSDs_.size() << std::endl; - for (SD = surfaceSDs_.begin(); SD != surfaceSDs_.end(); ++SD){ - Vector3d thisatom = (*SD)->getPos(); - std::cout << "Au " << thisatom.x() << " " << thisatom.y() << " " << thisatom.z() << std::endl; - } - */ - - - - Ns_ = surfaceSDs_.size(); - - - 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; - - - + << totlong << curlong << std::endl; } -void ConvexHull::printHull(const std::string& geomFileName) -{ - +void ConvexHull::printHull(const std::string& geomFileName) { FILE *newGeomFile; //create new .md file based on old .md file @@ -700,7 +302,3 @@ void ConvexHull::printHull(const std::string& geomFile fclose(newGeomFile); } #endif //QHULL -#endif //CGAL - - -