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

Comparing trunk/src/math/ConvexHull.cpp (file contents):
Revision 1302 by chuckv, Tue Oct 7 17:12:48 2008 UTC vs.
Revision 1371 by chuckv, Mon Oct 19 17:44:18 2009 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 2008 The University of Notre Dame. All Rights Reserved.
1 > /* Copyright (c) 2008, 2009 The University of Notre Dame. All Rights Reserved.
2   *
3   * The University of Notre Dame grants you ("Licensee") a
4   * non-exclusive, royalty free, license to use, modify and
# 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.9 2008-10-07 17:12:48 chuckv Exp $
47 > *  @version $Id: ConvexHull.cpp,v 1.15 2009-10-19 17:44:18 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 297 | Line 297 | void ConvexHull::computeHull(std::vector<StuntDouble*>
297  
298    }
299    
300 –  std::cout << "Number of surface atoms is: " << Ns_ << std::endl;
300    
301  
302  
# Line 369 | Line 368 | void ConvexHull::printHull(const std::string& geomFile
368   #else
369   #ifdef HAVE_QHULL
370   /* Old options Qt Qu Qg QG0 FA */
371 < ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt  Qci Tcv Pp"), Ns_(200) {
371 > /* More old opts Qc Qi Pp*/
372 > ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp"), Ns_(200), nTriangles_(0) {
373    //If we are doing the mpi version, set up some vectors for data communication
374   #ifdef IS_MPI
375  
# Line 377 | Line 377 | ConvexHull::ConvexHull() : Hull(), dim_(3), options_("
377   nproc_ = MPI::COMM_WORLD.Get_size();
378   myrank_ = MPI::COMM_WORLD.Get_rank();
379   NstoProc_ = new int[nproc_];
380 < displs_   = new int[nproc_];
380 > vecdispls_   = new int[nproc_];
381 > vecNstoProc_ = new int[nproc_];
382 > displs_ = new int[nproc_];
383  
384   // Create a surface point type in MPI to send
385   //surfacePtType = MPI::DOUBLE.Create_contiguous(3);
# Line 437 | Line 439 | void ConvexHull::computeHull(std::vector<StuntDouble*>
439    
440    boolT ismalloc = False;
441    /* Clean up memory from previous convex hull calculations*/
442 +  
443    Triangles_.clear();
444    surfaceSDs_.clear();
445    surfaceSDs_.reserve(Ns_);
# Line 445 | Line 448 | void ConvexHull::computeHull(std::vector<StuntDouble*>
448                      const_cast<char *>(options_.c_str()), NULL, stderr)) {
449  
450        sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull");
451 <      painCave.isFatal = 0;
451 >      painCave.isFatal = 1;
452        simError();
453        
454    } //qh_new_qhull
# Line 453 | Line 456 | void ConvexHull::computeHull(std::vector<StuntDouble*>
456  
457   #ifdef IS_MPI
458    std::vector<double> localPts;
459 +  std::vector<double> localVel;
460 +  std::vector<double> localMass;
461    int localPtArraySize;
462    
463  
# Line 464 | Line 469 | void ConvexHull::computeHull(std::vector<StuntDouble*>
469      if (!facet->simplicial){
470        // should never happen with Qt
471        sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
472 <      painCave.isFatal = 0;
472 >      painCave.isFatal = 1;
473        simError();
474      }
475      
# Line 483 | Line 488 | void ConvexHull::computeHull(std::vector<StuntDouble*>
488  
489  
490  
486 –  /*
487 –  std::sort(surfaceIDs.begin(),surfaceIDs.end());
488 –  surfaceIDs.erase(std::unique(surfaceIDs.begin(), surfaceIDs.end()), surfaceIDs.end());
489 –  int localPtArraySize = surfaceIDs.size() * 3;
490 –  */
491  
492 –  //localPts.resize(localPtArraySize);
493 –  //std::fill(localPts.begin(),localPts.end(),0.0);
494 –
495 –
492    int idx = 0;
493    int nIsIts = 0;
498 – /*
499 –  // Copy the surface points into an array.
500 –  for(std::vector<bool>::iterator list_iter = isSurfaceID.begin();
501 –      list_iter != isSurfaceID.end(); list_iter++)
502 –    {
503 –      bool isIt = *list_iter;
504 –      if (isIt){
505 –        localPts.push_back(ptArray[dim_ * idx]);    
506 –        localPts.push_back(ptArray[dim_ * idx + 1]);
507 –        localPts.push_back(ptArray[dim_ * idx + 2]);
508 –        localPtsMap.push_back(idx);
509 –        nIsIts++;
510 –      } //Isit
511 –      idx++;
512 –    } //isSurfaceID
513 –  */
494    FORALLvertices {
495      idx = qh_pointid(vertex->point);
496      localPts.push_back(ptArray[dim_ * idx]);    
497      localPts.push_back(ptArray[dim_ * idx + 1]);
498      localPts.push_back(ptArray[dim_ * idx + 2]);
499 +
500 +    Vector3d vel = bodydoubles[idx]->getVel();
501 +    localVel.push_back(vel.x());
502 +    localVel.push_back(vel.y());
503 +    localVel.push_back(vel.z());
504 +
505 +
506 +    RealType bdmass = bodydoubles[idx]->getMass();
507 +    localMass.push_back(bdmass);
508 +
509      localPtsMap.push_back(idx);
510 +
511    }
512  
513  
514 <  localPtArraySize = localPts.size();
524 <
514 >  localPtArraySize = int(localPts.size()/3.0);
515  
516    MPI::COMM_WORLD.Allgather(&localPtArraySize,1,MPI::INT,&NstoProc_[0],1,MPI::INT);
517 <
517 >  
518    Nsglobal_=0;
519    for (int i = 0; i < nproc_; i++){
520      Nsglobal_ += NstoProc_[i];
521 +    vecNstoProc_[i] = NstoProc_[i]*3;
522    }
523    
524  
525 <  int nglobalPts = int(Nsglobal_/3);
525 >  int nglobalPts = Nsglobal_*3;
526  
527  
528 <  std::vector<double> globalPts;
529 <  globalPts.resize(Nsglobal_);
528 >  std::vector<double> globalPts(nglobalPts);
529 >  std::vector<double> globalVel(nglobalPts);
530 >  std::vector<double> globalMass(Nsglobal_);
531  
532 +
533 +  
534    isSurfaceID.resize(nglobalPts);
535  
536  
537    std::fill(globalPts.begin(),globalPts.end(),0.0);
538  
539 <  displs_[0] = 0;
539 >  vecdispls_[0] = 0;
540    /* Build a displacements array */
541    for (int i = 1; i < nproc_; i++){
542 <    displs_[i] = displs_[i-1] + NstoProc_[i-1];
542 >    vecdispls_[i] = vecdispls_[i-1] + vecNstoProc_[i-1];
543    }
544    
545 +  displs_[0] = 0;
546 +  for (int i = 1; i < nproc_; i++){
547 +    displs_[i] = displs_[i-1] + NstoProc_[i-1];
548 +  }
549    
550    int noffset = displs_[myrank_];
551    /* gather the potential hull */
552    
553 <  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
553 >  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize*3,MPI::DOUBLE,&globalPts[0],&vecNstoProc_[0],&vecdispls_[0],MPI::DOUBLE);
554 >  MPI::COMM_WORLD.Allgatherv(&localVel[0],localPtArraySize*3,MPI::DOUBLE,&globalVel[0],&vecNstoProc_[0],&vecdispls_[0],MPI::DOUBLE);
555 >  MPI::COMM_WORLD.Allgatherv(&localMass[0],localPtArraySize,MPI::DOUBLE,&globalMass[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
556  
557    /*
558 +  int tmpidx = 0;
559 +  
560    if (myrank_ == 0){
561 <    for (i = 0; i < globalPts.size(); i++){
562 <      std::cout << globalPts[i] << std::endl;
561 >    for (i = 0; i < nglobalPts-3; i++){      
562 >      std::cout << "Au   " << globalPts[tmpidx] << "  " << globalPts[tmpidx+1] << "  " << globalPts[tmpidx +2] << std::endl;
563 >      tmpidx = tmpidx + 3;
564      }
565    }
566    */
567 +  
568    // Free previous hull
569    qh_freeqhull(!qh_ALL);
570    qh_memfreeshort(&curlong, &totlong);
# Line 568 | Line 572 | void ConvexHull::computeHull(std::vector<StuntDouble*>
572      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
573                << totlong << curlong << std::endl;
574  
575 <  if (qh_new_qhull(dim_, nglobalPts, &globalPts[0], ismalloc,
575 >  if (qh_new_qhull(dim_, Nsglobal_, &globalPts[0], ismalloc,
576                      const_cast<char *>(options_.c_str()), NULL, stderr)){
577  
578        sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute global convex hull");
# Line 585 | Line 589 | void ConvexHull::computeHull(std::vector<StuntDouble*>
589  
590  
591      unsigned int nf = qh num_facets;
592 <    
592 >    
593      /* Build Surface SD list first */
594  
595      std::fill(isSurfaceID.begin(),isSurfaceID.end(),false);
596 <
596 >    int numvers = 0;
597      FORALLfacets {
598        
599        if (!facet->simplicial){
# Line 599 | Line 603 | void ConvexHull::computeHull(std::vector<StuntDouble*>
603          simError();
604        } //simplicical
605        
606 <      Triangle* face = new Triangle();
606 >      Triangle face;
607        Vector3d  V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]);
608 <      face->setNormal(V3dNormal);
609 <      //face->setCentroid(V3dCentroid);
610 <      RealType faceArea = 0.5*V3dNormal.length();
607 <      //face->setArea(faceArea);
608 >      face.setNormal(V3dNormal);
609 >
610 >      
611  
612 +      //RealType faceArea = 0.5*V3dNormal.length();
613 +      RealType faceArea = qh_facetarea(facet);
614 +      face.setArea(faceArea);
615  
616 +
617        vertices = qh_facet3vertex(facet);
618 +      
619 +      coordT *center = qh_getcenter(vertices);
620 +      Vector3d V3dCentroid(center[0], center[1], center[2]);
621 +      face.setCentroid(V3dCentroid);
622 +      Vector3d faceVel = V3Zero;
623 +      Vector3d p[3];
624 +      RealType faceMass = 0.0;
625 +      int ver = 0;
626        FOREACHvertex_(vertices){
627          id = qh_pointid(vertex->point);
628 +        p[ver][0] = vertex->point[0];
629 +        p[ver][1] = vertex->point[1];
630 +        p[ver][2] = vertex->point[2];
631          int localindex = id;
632   #ifdef IS_MPI
633 +        Vector3d velVector(globalVel[dim_ * id],globalVel[dim_ * id + 1], globalVel[dim_ * id + 2]);
634          
635 <        if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){
636 <          localindex = localPtsMap[id-noffset/3];
635 >        faceVel = faceVel + velVector;
636 >        faceMass = faceMass + globalMass[id];
637 >        if (id >= noffset && id < (noffset + localPtArraySize) ){
638 >          
639 >          localindex = localPtsMap[id-noffset];
640 > #else
641 >          faceVel = faceVel + bodydoubles[localindex]->getVel();
642 >          faceMass = faceMass + bodydoubles[localindex]->getMass();
643   #endif
644 <          face->addVertex(bodydoubles[localindex]);
644 >          face.addVertexSD(bodydoubles[localindex]);
645            if( !isSurfaceID[id] ){
646              isSurfaceID[id] = true;
647   #ifdef IS_MPI      
# Line 624 | Line 649 | void ConvexHull::computeHull(std::vector<StuntDouble*>
649   #endif
650              
651              surfaceSDs_.push_back(bodydoubles[localindex]);
652 +            //   std::cout <<"This ID is: " << bodydoubles[localindex]->getGlobalIndex() << std::endl;
653              
654            } //IF isSurfaceID
655  
656   #ifdef IS_MPI
657          
658          }else{
659 <          face->addVertex(NULL);
659 >          face.addVertexSD(NULL);
660            }
661   #endif
662 +        numvers++;
663 +        ver++;
664        } //Foreachvertex
665 <
665 >      /*
666 >      if (!SETempty_(facet->coplanarset)){
667 >        FOREACHpoint_(facet->coplanarset){
668 >          id = qh_pointid(point);
669 >          surfaceSDs_.push_back(bodydoubles[id]);
670 >        }
671 >      }
672 >      */
673 >      face.addVertices(p[0],p[1],p[2]);
674 >      face.setFacetMass(faceMass);
675 >      face.setFacetVelocity(faceVel/3.0);
676        Triangles_.push_back(face);
677        qh_settempfree(&vertices);      
678 <      
678 >
679      } //FORALLfacets
680  
681 <                        
681 >    /*    
682 >    std::cout << surfaceSDs_.size() << std::endl;
683 >    for (SD = surfaceSDs_.begin(); SD != surfaceSDs_.end(); ++SD){
684 >      Vector3d thisatom = (*SD)->getPos();
685 >      std::cout << "Au " << thisatom.x() << "  " << thisatom.y() << " " << thisatom.z() << std::endl;
686 >    }
687 >    */
688  
645 –  Ns_ = surfaceSDs_.size();
689  
690 <
691 <  qh_getarea(qh facet_list);
692 <  volume_ = qh totvol;
693 <  area_ = qh totarea;
694 <
695 <  
696 <
697 <  qh_freeqhull(!qh_ALL);
698 <  qh_memfreeshort(&curlong, &totlong);
699 <  if (curlong || totlong)
700 <    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
701 <              << totlong << curlong << std::endl;
702 <  
703 <
704 <  
690 >    Ns_ = surfaceSDs_.size();
691 >    nTriangles_ = Triangles_.size();
692 >    
693 >    qh_getarea(qh facet_list);
694 >    volume_ = qh totvol;
695 >    area_ = qh totarea;
696 >    
697 >    
698 >    
699 >    qh_freeqhull(!qh_ALL);
700 >    qh_memfreeshort(&curlong, &totlong);
701 >    if (curlong || totlong)
702 >      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
703 >                << totlong << curlong << std::endl;
704 >    
705 >    
706 >    
707   }
708  
709  

Diff Legend

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