--- trunk/src/math/ConvexHull.cpp 2008/10/20 19:36:32 1307 +++ 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.11 2008-10-20 19:36:32 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; @@ -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_); @@ -455,6 +457,7 @@ void ConvexHull::computeHull(std::vector #ifdef IS_MPI std::vector localPts; std::vector localVel; + std::vector localMass; int localPtArraySize; @@ -485,34 +488,9 @@ 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]); @@ -524,52 +502,69 @@ void ConvexHull::computeHull(std::vector 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(Nsglobal_); - std::vector globalVel(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(&localVel[0],localPtArraySize,MPI::DOUBLE,&globalVel[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); @@ -577,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"); @@ -594,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){ @@ -608,34 +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 + 1]); + Vector3d velVector(globalVel[dim_ * id],globalVel[dim_ * id + 1], globalVel[dim_ * id + 2]); + faceVel = faceVel + velVector; - if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){ - localindex = localPtsMap[id-noffset/3]; + 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 @@ -643,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)){ @@ -661,13 +670,15 @@ void ConvexHull::computeHull(std::vector } } */ - face->setFacetVelocity(faceVel/3.0); + 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(); @@ -676,7 +687,6 @@ void ConvexHull::computeHull(std::vector */ - Ns_ = surfaceSDs_.size(); nTriangles_ = Triangles_.size();