--- trunk/src/math/ConvexHull.cpp 2008/10/15 18:26:01 1304 +++ trunk/src/math/ConvexHull.cpp 2009/10/19 17:44:18 1371 @@ -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,7 +44,7 @@ * * 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.15 2009-10-19 17:44:18 chuckv Exp $ * */ @@ -149,8 +149,8 @@ ConvexHull::ConvexHull() : Hull(){ nproc_ = MPI::COMM_WORLD.Get_size(); myrank_ = MPI::COMM_WORLD.Get_rank(); NstoProc_ = new int[nproc_]; - displs_ = new int[nproc_]; - + vecdispls_ = new int[nproc_]; + displs_ = new int[nproc_]; // Create a surface point type in MPI to send surfacePtType = MPI::DOUBLE.Create_contiguous(3); surfacePtType.Commit(); @@ -213,10 +213,10 @@ void ConvexHull::computeHull(std::vector /* Build a displacements array */ for (int i = 1; i < nproc_; i++){ - displs_[i] = displs_[i-1] + NstoProc_[i-1]; + vecdispls_[i] = vecdispls_[i-1] + NstoProc_[i-1]; } - int noffset = displs_[myrank_]; + int noffset = vecdispls_[myrank_]; /* gather the potential hull */ @@ -229,7 +229,7 @@ void ConvexHull::computeHull(std::vector surfacePtsLocal_.push_back(mpiSurfacePt); } - MPI::COMM_WORLD.Allgatherv(&surfacePtsLocal_[0],Ns_,surfacePtType,&surfacePtsGlobal_[0],NstoProc_,displs_,surfacePtType); + MPI::COMM_WORLD.Allgatherv(&surfacePtsLocal_[0],Ns_,surfacePtType,&surfacePtsGlobal_[0],NstoProc_,vecdispls_,surfacePtType); std::vector::iterator spt; std::vector gblpoints; @@ -297,7 +297,6 @@ void ConvexHull::computeHull(std::vector } - std::cout << "Number of surface atoms is: " << Ns_ << std::endl; @@ -370,7 +369,7 @@ void ConvexHull::printHull(const std::string& geomFile #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) { +ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp"), Ns_(200), nTriangles_(0) { //If we are doing the mpi version, set up some vectors for data communication #ifdef IS_MPI @@ -378,7 +377,9 @@ ConvexHull::ConvexHull() : Hull(), dim_(3), options_(" nproc_ = MPI::COMM_WORLD.Get_size(); myrank_ = MPI::COMM_WORLD.Get_rank(); NstoProc_ = new int[nproc_]; - displs_ = new int[nproc_]; + vecdispls_ = new int[nproc_]; + vecNstoProc_ = new int[nproc_]; + displs_ = new int[nproc_]; // Create a surface point type in MPI to send //surfacePtType = MPI::DOUBLE.Create_contiguous(3); @@ -438,6 +439,7 @@ void ConvexHull::computeHull(std::vector boolT ismalloc = False; /* Clean up memory from previous convex hull calculations*/ + Triangles_.clear(); surfaceSDs_.clear(); surfaceSDs_.reserve(Ns_); @@ -446,7 +448,7 @@ void ConvexHull::computeHull(std::vector const_cast(options_.c_str()), NULL, stderr)) { sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull"); - painCave.isFatal = 0; + painCave.isFatal = 1; simError(); } //qh_new_qhull @@ -454,6 +456,8 @@ void ConvexHull::computeHull(std::vector #ifdef IS_MPI std::vector localPts; + std::vector localVel; + std::vector localMass; int localPtArraySize; @@ -465,7 +469,7 @@ void ConvexHull::computeHull(std::vector if (!facet->simplicial){ // should never happen with Qt sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected"); - painCave.isFatal = 0; + painCave.isFatal = 1; simError(); } @@ -484,84 +488,83 @@ void ConvexHull::computeHull(std::vector - /* - 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]); + + Vector3d vel = bodydoubles[idx]->getVel(); + localVel.push_back(vel.x()); + localVel.push_back(vel.y()); + localVel.push_back(vel.z()); + + + RealType bdmass = bodydoubles[idx]->getMass(); + localMass.push_back(bdmass); + localPtsMap.push_back(idx); + } - localPtArraySize = localPts.size(); - + localPtArraySize = int(localPts.size()/3.0); 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]; + vecNstoProc_[i] = NstoProc_[i]*3; } - int nglobalPts = int(Nsglobal_/3); + int nglobalPts = Nsglobal_*3; - std::vector globalPts; - globalPts.resize(Nsglobal_); + std::vector globalPts(nglobalPts); + std::vector globalVel(nglobalPts); + std::vector globalMass(Nsglobal_); + + isSurfaceID.resize(nglobalPts); std::fill(globalPts.begin(),globalPts.end(),0.0); - displs_[0] = 0; + vecdispls_[0] = 0; /* Build a displacements array */ for (int i = 1; i < nproc_; i++){ - displs_[i] = displs_[i-1] + NstoProc_[i-1]; + vecdispls_[i] = vecdispls_[i-1] + vecNstoProc_[i-1]; } + displs_[0] = 0; + for (int i = 1; i < nproc_; i++){ + displs_[i] = displs_[i-1] + NstoProc_[i-1]; + } 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); + MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize*3,MPI::DOUBLE,&globalPts[0],&vecNstoProc_[0],&vecdispls_[0],MPI::DOUBLE); + MPI::COMM_WORLD.Allgatherv(&localVel[0],localPtArraySize*3,MPI::DOUBLE,&globalVel[0],&vecNstoProc_[0],&vecdispls_[0],MPI::DOUBLE); + MPI::COMM_WORLD.Allgatherv(&localMass[0],localPtArraySize,MPI::DOUBLE,&globalMass[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE); /* + int tmpidx = 0; + if (myrank_ == 0){ - for (i = 0; i < globalPts.size(); i++){ - std::cout << globalPts[i] << std::endl; + for (i = 0; i < nglobalPts-3; i++){ + std::cout << "Au " << globalPts[tmpidx] << " " << globalPts[tmpidx+1] << " " << globalPts[tmpidx +2] << std::endl; + tmpidx = tmpidx + 3; } } */ + // Free previous hull qh_freeqhull(!qh_ALL); qh_memfreeshort(&curlong, &totlong); @@ -569,7 +572,7 @@ void ConvexHull::computeHull(std::vector 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, + if (qh_new_qhull(dim_, Nsglobal_, &globalPts[0], ismalloc, const_cast(options_.c_str()), NULL, stderr)){ sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute global convex hull"); @@ -586,11 +589,11 @@ void ConvexHull::computeHull(std::vector unsigned int nf = qh num_facets; - + /* Build Surface SD list first */ std::fill(isSurfaceID.begin(),isSurfaceID.end(),false); - + int numvers = 0; FORALLfacets { if (!facet->simplicial){ @@ -600,31 +603,45 @@ void ConvexHull::computeHull(std::vector simError(); } //simplicical - Triangle* face = new Triangle(); + Triangle face; Vector3d V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]); - face->setNormal(V3dNormal); + face.setNormal(V3dNormal); - RealType faceArea = 0.5*V3dNormal.length(); - face->setArea(faceArea); + //RealType faceArea = 0.5*V3dNormal.length(); + 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); - + face.setCentroid(V3dCentroid); + Vector3d faceVel = V3Zero; + Vector3d p[3]; + RealType faceMass = 0.0; + int ver = 0; FOREACHvertex_(vertices){ id = qh_pointid(vertex->point); + p[ver][0] = vertex->point[0]; + p[ver][1] = vertex->point[1]; + p[ver][2] = vertex->point[2]; int localindex = id; #ifdef IS_MPI + Vector3d velVector(globalVel[dim_ * id],globalVel[dim_ * id + 1], globalVel[dim_ * id + 2]); - if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){ - localindex = localPtsMap[id-noffset/3]; + faceVel = faceVel + velVector; + faceMass = faceMass + globalMass[id]; + if (id >= noffset && id < (noffset + localPtArraySize) ){ + + localindex = localPtsMap[id-noffset]; +#else + faceVel = faceVel + bodydoubles[localindex]->getVel(); + faceMass = faceMass + bodydoubles[localindex]->getMass(); #endif - face->addVertex(bodydoubles[localindex]); + face.addVertexSD(bodydoubles[localindex]); if( !isSurfaceID[id] ){ isSurfaceID[id] = true; #ifdef IS_MPI @@ -632,15 +649,18 @@ void ConvexHull::computeHull(std::vector #endif surfaceSDs_.push_back(bodydoubles[localindex]); + // std::cout <<"This ID is: " << bodydoubles[localindex]->getGlobalIndex() << std::endl; } //IF isSurfaceID #ifdef IS_MPI }else{ - face->addVertex(NULL); + face.addVertexSD(NULL); } #endif + numvers++; + ver++; } //Foreachvertex /* if (!SETempty_(facet->coplanarset)){ @@ -649,13 +669,16 @@ void ConvexHull::computeHull(std::vector surfaceSDs_.push_back(bodydoubles[id]); } } - + */ + 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::cout << surfaceSDs_.size() << std::endl; for (SD = surfaceSDs_.begin(); SD != surfaceSDs_.end(); ++SD){ Vector3d thisatom = (*SD)->getPos(); @@ -664,10 +687,9 @@ void ConvexHull::computeHull(std::vector */ - Ns_ = surfaceSDs_.size(); + nTriangles_ = Triangles_.size(); - qh_getarea(qh facet_list); volume_ = qh totvol; area_ = qh totarea;