ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/math/ConvexHull.cpp
(Generate patch)

Comparing trunk/src/math/ConvexHull.cpp (file contents):
Revision 1308 by chuckv, Tue Oct 21 16:44:00 2008 UTC vs.
Revision 1316 by chuckv, Fri Nov 14 15:44:34 2008 UTC

# Line 44 | Line 44
44   *
45   *  Created by Charles F. Vardeman II on 11 Dec 2006.
46   *  @author  Charles F. Vardeman II
47 < *  @version $Id: ConvexHull.cpp,v 1.12 2008-10-21 16:44:00 chuckv Exp $
47 > *  @version $Id: ConvexHull.cpp,v 1.13 2008-11-14 15:44:34 chuckv Exp $
48   *
49   */
50  
# Line 149 | Line 149 | ConvexHull::ConvexHull() : Hull(){
149   nproc_ = MPI::COMM_WORLD.Get_size();
150   myrank_ = MPI::COMM_WORLD.Get_rank();
151   NstoProc_ = new int[nproc_];
152 < displs_   = new int[nproc_];
153 <
152 > vecdispls_   = new int[nproc_];
153 > displs_ = new int[nproc_];
154   // Create a surface point type in MPI to send
155   surfacePtType = MPI::DOUBLE.Create_contiguous(3);
156   surfacePtType.Commit();
# Line 213 | Line 213 | void ConvexHull::computeHull(std::vector<StuntDouble*>
213  
214    /* Build a displacements array */
215    for (int i = 1; i < nproc_; i++){
216 <    displs_[i] = displs_[i-1] + NstoProc_[i-1];
216 >    vecdispls_[i] = vecdispls_[i-1] + NstoProc_[i-1];
217    }
218    
219 <  int noffset = displs_[myrank_];
219 >  int noffset = vecdispls_[myrank_];
220    /* gather the potential hull */
221    
222    
# Line 229 | Line 229 | void ConvexHull::computeHull(std::vector<StuntDouble*>
229      surfacePtsLocal_.push_back(mpiSurfacePt);
230    }
231  
232 <  MPI::COMM_WORLD.Allgatherv(&surfacePtsLocal_[0],Ns_,surfacePtType,&surfacePtsGlobal_[0],NstoProc_,displs_,surfacePtType);
232 >  MPI::COMM_WORLD.Allgatherv(&surfacePtsLocal_[0],Ns_,surfacePtType,&surfacePtsGlobal_[0],NstoProc_,vecdispls_,surfacePtType);
233    std::vector<surfacePt_>::iterator spt;
234    std::vector<Enriched_Point_3> gblpoints;
235  
# Line 378 | Line 378 | ConvexHull::ConvexHull() : Hull(), dim_(3), options_("
378   nproc_ = MPI::COMM_WORLD.Get_size();
379   myrank_ = MPI::COMM_WORLD.Get_rank();
380   NstoProc_ = new int[nproc_];
381 < displs_   = new int[nproc_];
381 > vecdispls_   = new int[nproc_];
382 > vecNstoProc_ = new int[nproc_];
383 > displs_ = new int[nproc_];
384  
385   // Create a surface point type in MPI to send
386   //surfacePtType = MPI::DOUBLE.Create_contiguous(3);
# Line 456 | Line 458 | void ConvexHull::computeHull(std::vector<StuntDouble*>
458   #ifdef IS_MPI
459    std::vector<double> localPts;
460    std::vector<double> localVel;
461 +  std::vector<double> localMass;
462    int localPtArraySize;
463    
464  
# Line 485 | Line 488 | void ConvexHull::computeHull(std::vector<StuntDouble*>
488    } //FORALLfacets
489  
490  
488 –
489 –  /*
490 –  std::sort(surfaceIDs.begin(),surfaceIDs.end());
491 –  surfaceIDs.erase(std::unique(surfaceIDs.begin(), surfaceIDs.end()), surfaceIDs.end());
492 –  int localPtArraySize = surfaceIDs.size() * 3;
493 –  */
494 –
495 –  //localPts.resize(localPtArraySize);
496 –  //std::fill(localPts.begin(),localPts.end(),0.0);
491  
492  
493    int idx = 0;
494    int nIsIts = 0;
501 – /*
502 –  // Copy the surface points into an array.
503 –  for(std::vector<bool>::iterator list_iter = isSurfaceID.begin();
504 –      list_iter != isSurfaceID.end(); list_iter++)
505 –    {
506 –      bool isIt = *list_iter;
507 –      if (isIt){
508 –        localPts.push_back(ptArray[dim_ * idx]);    
509 –        localPts.push_back(ptArray[dim_ * idx + 1]);
510 –        localPts.push_back(ptArray[dim_ * idx + 2]);
511 –        localPtsMap.push_back(idx);
512 –        nIsIts++;
513 –      } //Isit
514 –      idx++;
515 –    } //isSurfaceID
516 –  */
495    FORALLvertices {
496      idx = qh_pointid(vertex->point);
497      localPts.push_back(ptArray[dim_ * idx]);    
# Line 525 | Line 503 | void ConvexHull::computeHull(std::vector<StuntDouble*>
503      localVel.push_back(vel.y());
504      localVel.push_back(vel.z());
505  
506 +    RealType bdmass = bodydoubles[idx]->getMass();
507 +    localMass.push_back(bdmass);
508 +
509      localPtsMap.push_back(idx);
510 +
511 +
512    }
513  
514  
532 –  localPtArraySize = localPts.size();
515  
516 +  localPtArraySize = int(localPts.size()/3.0);
517 +
518  
519    MPI::COMM_WORLD.Allgather(&localPtArraySize,1,MPI::INT,&NstoProc_[0],1,MPI::INT);
520 <
520 >  
521    Nsglobal_=0;
522    for (int i = 0; i < nproc_; i++){
523      Nsglobal_ += NstoProc_[i];
524 +    vecNstoProc_[i] = NstoProc_[i]*3;
525    }
526    
527  
528 <  int nglobalPts = int(Nsglobal_/3);
528 >  int nglobalPts = Nsglobal_*3;
529  
530  
531 <  std::vector<double> globalPts(Nsglobal_);
532 <  std::vector<double> globalVel(Nsglobal_);
531 >  std::vector<double> globalPts(nglobalPts);
532 >  std::vector<double> globalVel(nglobalPts);
533 >  std::vector<double> globalMass(Nsglobal_);
534  
535    isSurfaceID.resize(nglobalPts);
536  
537  
538    std::fill(globalPts.begin(),globalPts.end(),0.0);
539  
540 <  displs_[0] = 0;
540 >  vecdispls_[0] = 0;
541    /* Build a displacements array */
542    for (int i = 1; i < nproc_; i++){
543 <    displs_[i] = displs_[i-1] + NstoProc_[i-1];
543 >    vecdispls_[i] = vecdispls_[i-1] + vecNstoProc_[i-1];
544    }
545    
546 +  displs_[0] = 0;
547 +  for (int i = 1; i < nproc_; i++){
548 +    displs_[i] = displs_[i-1] + NstoProc_[i-1];
549 +  }
550    
551 <  int noffset = displs_[myrank_];
551 >  int noffset = vecdispls_[myrank_];
552    /* gather the potential hull */
553    
554 <  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
555 <  MPI::COMM_WORLD.Allgatherv(&localVel[0],localPtArraySize,MPI::DOUBLE,&globalVel[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
556 <
554 >  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&vecNstoProc_[0],&vecdispls_[0],MPI::DOUBLE);
555 >  MPI::COMM_WORLD.Allgatherv(&localVel[0],localPtArraySize,MPI::DOUBLE,&globalVel[0],&vecNstoProc_[0],&vecdispls_[0],MPI::DOUBLE);
556 >  MPI::COMM_WORLD.Allgatherv(&localMass[0],localPtArraySize,MPI::DOUBLE,&globalMass[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
557    /*
558    if (myrank_ == 0){
559      for (i = 0; i < globalPts.size(); i++){
# Line 615 | Line 605 | void ConvexHull::computeHull(std::vector<StuntDouble*>
605  
606        
607  
608 <      RealType faceArea = 0.5*V3dNormal.length();
608 >      //RealType faceArea = 0.5*V3dNormal.length();
609 >      RealType faceArea = qh_facetarea(facet);
610        face.setArea(faceArea);
611  
612  
# Line 625 | Line 616 | void ConvexHull::computeHull(std::vector<StuntDouble*>
616        Vector3d V3dCentroid(center[0], center[1], center[2]);
617        face.setCentroid(V3dCentroid);
618        Vector3d faceVel = V3Zero;
619 +      Vector3d p[3];
620 +      RealType faceMass = 0.0;
621 +      int ver = 0;
622        FOREACHvertex_(vertices){
623          id = qh_pointid(vertex->point);
624 +        p[ver][0] = vertex->point[0];
625 +        p[ver][1] = vertex->point[1];
626 +        p[ver][2] = vertex->point[2];
627          int localindex = id;
628   #ifdef IS_MPI
629          Vector3d velVector(globalVel[dim_ * id],globalVel[dim_ * id + 1], globalVel[dim_ * id + 1]);
630 +        
631          faceVel = faceVel + velVector;
632 +        faceMass = faceMass + globalMass[id];
633          if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){
634            localindex = localPtsMap[id-noffset/3];
635   #else
636            faceVel = faceVel + bodydoubles[localindex]->getVel();
637 +          faceMass = faceMass + bodydoubles[localindex]->getMass();
638   #endif
639 <          face.addVertex(bodydoubles[localindex]);
639 >          face.addVertexSD(bodydoubles[localindex]);
640            if( !isSurfaceID[id] ){
641              isSurfaceID[id] = true;
642   #ifdef IS_MPI      
# Line 650 | Line 650 | void ConvexHull::computeHull(std::vector<StuntDouble*>
650   #ifdef IS_MPI
651          
652          }else{
653 <          face.addVertex(NULL);
653 >          face.addVertexSD(NULL);
654            }
655   #endif
656 +        ver++;
657        } //Foreachvertex
658        /*
659        if (!SETempty_(facet->coplanarset)){
# Line 662 | Line 663 | void ConvexHull::computeHull(std::vector<StuntDouble*>
663          }
664        }
665        */
666 +      face.addVertices(p[0],p[1],p[2]);
667 +      face.setFacetMass(faceMass);
668        face.setFacetVelocity(faceVel/3.0);
669        Triangles_.push_back(face);
670        qh_settempfree(&vertices);      

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines