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 1293 by chuckv, Sun Sep 14 01:32:26 2008 UTC vs.
Revision 1304 by chuckv, Wed Oct 15 18:26:01 2008 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 2006 The University of Notre Dame. All Rights Reserved.
1 > /* Copyright (c) 2008 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.8 2008-09-14 01:32:25 chuckv Exp $
47 > *  @version $Id: ConvexHull.cpp,v 1.10 2008-10-15 18:26:01 chuckv Exp $
48   *
49   */
50  
# Line 369 | Line 369 | void ConvexHull::printHull(const std::string& geomFile
369   #else
370   #ifdef HAVE_QHULL
371   /* Old options Qt Qu Qg QG0 FA */
372 < ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt  Qci Tcv Pp") {
372 > /* More old opts Qc Qi Pp*/
373 > ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp"), Ns_(200) {
374    //If we are doing the mpi version, set up some vectors for data communication
375   #ifdef IS_MPI
376  
# Line 380 | Line 381 | ConvexHull::ConvexHull() : Hull(), dim_(3), options_("
381   displs_   = new int[nproc_];
382  
383   // Create a surface point type in MPI to send
384 < surfacePtType = MPI::DOUBLE.Create_contiguous(3);
385 < surfacePtType.Commit();
384 > //surfacePtType = MPI::DOUBLE.Create_contiguous(3);
385 > // surfacePtType.Commit();
386  
387  
388   #endif
# Line 413 | Line 414 | void ConvexHull::computeHull(std::vector<StuntDouble*>
414    
415    //pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_));
416  
417 <  double* ptArray = new double[numpoints * 3];
417 > //  double* ptArray = new double[numpoints * 3];
418 >  std::vector<double> ptArray(numpoints*3);
419    std::vector<bool> isSurfaceID(numpoints);
420  
421    // Copy the positon vector into a points vector for qhull.
# Line 435 | Line 437 | void ConvexHull::computeHull(std::vector<StuntDouble*>
437    
438    
439    boolT ismalloc = False;
440 +  /* Clean up memory from previous convex hull calculations*/
441    Triangles_.clear();
442    surfaceSDs_.clear();
443 <  
444 <  if (qh_new_qhull(dim_, numpoints, ptArray, ismalloc,
443 >  surfaceSDs_.reserve(Ns_);
444 >
445 >  if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc,
446                      const_cast<char *>(options_.c_str()), NULL, stderr)) {
447  
448        sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull");
# Line 452 | Line 456 | void ConvexHull::computeHull(std::vector<StuntDouble*>
456    std::vector<double> localPts;
457    int localPtArraySize;
458    
459 +
460    std::fill(isSurfaceID.begin(),isSurfaceID.end(),false);
461 +
462  
463    FORALLfacets {
464      
# Line 484 | Line 490 | void ConvexHull::computeHull(std::vector<StuntDouble*>
490    int localPtArraySize = surfaceIDs.size() * 3;
491    */
492  
493 <  localPts.resize(localPtArraySize);
494 <  //  std::fill(localPts.begin(),globalPts.end(),0.0);
493 >  //localPts.resize(localPtArraySize);
494 >  //std::fill(localPts.begin(),localPts.end(),0.0);
495  
496  
497    int idx = 0;
498 +  int nIsIts = 0;
499 + /*
500    // Copy the surface points into an array.
501    for(std::vector<bool>::iterator list_iter = isSurfaceID.begin();
502        list_iter != isSurfaceID.end(); list_iter++)
# Line 499 | Line 507 | void ConvexHull::computeHull(std::vector<StuntDouble*>
507          localPts.push_back(ptArray[dim_ * idx + 1]);
508          localPts.push_back(ptArray[dim_ * idx + 2]);
509          localPtsMap.push_back(idx);
510 +        nIsIts++;
511        } //Isit
512        idx++;
513      } //isSurfaceID
514 <  
514 >  */
515 >  FORALLvertices {
516 >    idx = qh_pointid(vertex->point);
517 >    localPts.push_back(ptArray[dim_ * idx]);    
518 >    localPts.push_back(ptArray[dim_ * idx + 1]);
519 >    localPts.push_back(ptArray[dim_ * idx + 2]);
520 >    localPtsMap.push_back(idx);
521 >  }
522  
523 +
524    localPtArraySize = localPts.size();
525 +
526 +
527    MPI::COMM_WORLD.Allgather(&localPtArraySize,1,MPI::INT,&NstoProc_[0],1,MPI::INT);
528  
529 +  Nsglobal_=0;
530    for (int i = 0; i < nproc_; i++){
531      Nsglobal_ += NstoProc_[i];
532    }
533 +  
534 +
535 +  int nglobalPts = int(Nsglobal_/3);
536 +
537  
538    std::vector<double> globalPts;
539    globalPts.resize(Nsglobal_);
540 <  isSurfaceID.resize(int(Nsglobal_/3));
540 >
541 >  isSurfaceID.resize(nglobalPts);
542 >
543 >
544    std::fill(globalPts.begin(),globalPts.end(),0.0);
545 +
546 +  displs_[0] = 0;
547    /* Build a displacements array */
548    for (int i = 1; i < nproc_; i++){
549      displs_[i] = displs_[i-1] + NstoProc_[i-1];
550    }
551    
552 +  
553    int noffset = displs_[myrank_];
554    /* gather the potential hull */
555    
556 <  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],NstoProc_,displs_,MPI::DOUBLE);
556 >  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
557  
558 <
559 <
560 <
558 >  /*
559 >  if (myrank_ == 0){
560 >    for (i = 0; i < globalPts.size(); i++){
561 >      std::cout << globalPts[i] << std::endl;
562 >    }
563 >  }
564 >  */
565    // Free previous hull
566    qh_freeqhull(!qh_ALL);
567    qh_memfreeshort(&curlong, &totlong);
# Line 535 | Line 569 | void ConvexHull::computeHull(std::vector<StuntDouble*>
569      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
570                << totlong << curlong << std::endl;
571  
572 <  if (qh_new_qhull(dim_, numpoints, &globalPts[0], ismalloc,
572 >  if (qh_new_qhull(dim_, nglobalPts, &globalPts[0], ismalloc,
573                      const_cast<char *>(options_.c_str()), NULL, stderr)){
574  
575        sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute global convex hull");
576 <      painCave.isFatal = 0;
576 >      painCave.isFatal = 1;
577        simError();
578        
579    } //qh_new_qhull
# Line 552 | Line 586 | void ConvexHull::computeHull(std::vector<StuntDouble*>
586  
587  
588      unsigned int nf = qh num_facets;
589 <    
589 >    
590      /* Build Surface SD list first */
591  
592      std::fill(isSurfaceID.begin(),isSurfaceID.end(),false);
# Line 562 | Line 596 | void ConvexHull::computeHull(std::vector<StuntDouble*>
596        if (!facet->simplicial){
597        // should never happen with Qt
598          sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
599 <        painCave.isFatal = 0;
599 >        painCave.isFatal = 1;
600          simError();
601        } //simplicical
602        
603        Triangle* face = new Triangle();
604        Vector3d  V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]);
605        face->setNormal(V3dNormal);
606 <      //face->setCentroid(V3dCentroid);
606 >
607 >      
608 >
609        RealType faceArea = 0.5*V3dNormal.length();
610        face->setArea(faceArea);
611  
612  
613        vertices = qh_facet3vertex(facet);
614 +      
615 +      coordT *center = qh_getcenter(vertices);
616 +      Vector3d V3dCentroid(center[0], center[1], center[2]);
617 +      face->setCentroid(V3dCentroid);
618 +
619        FOREACHvertex_(vertices){
620          id = qh_pointid(vertex->point);
621 <        face->addVertex(bodydoubles[id]);
622 <        if( !isSurfaceID[id] ){
623 <          isSurfaceID[id] = true;
624 <          surfaceSDs_.push_back(bodydoubles[id]);
625 <        } //IF isSurfaceID
626 <      } //Foreachvertex
621 >        int localindex = id;
622 > #ifdef IS_MPI
623 >        
624 >        if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){
625 >          localindex = localPtsMap[id-noffset/3];
626 > #endif
627 >          face->addVertex(bodydoubles[localindex]);
628 >          if( !isSurfaceID[id] ){
629 >            isSurfaceID[id] = true;
630 > #ifdef IS_MPI      
631 >            
632 > #endif
633 >            
634 >            surfaceSDs_.push_back(bodydoubles[localindex]);
635 >            
636 >          } //IF isSurfaceID
637  
638 + #ifdef IS_MPI
639 +        
640 +        }else{
641 +          face->addVertex(NULL);
642 +          }
643 + #endif
644 +      } //Foreachvertex
645 +      /*
646 +      if (!SETempty_(facet->coplanarset)){
647 +        FOREACHpoint_(facet->coplanarset){
648 +          id = qh_pointid(point);
649 +          surfaceSDs_.push_back(bodydoubles[id]);
650 +        }
651 +      }
652  
653        Triangles_.push_back(face);
654        qh_settempfree(&vertices);      
655 <      
655 >      */
656      } //FORALLfacets
592 –                
657  
594 –                
595 –                
658      /*
659 <  std::sort(surfaceIDs.begin(),surfaceIDs.end());
660 <  surfaceIDs.erase(std::unique(surfaceIDs.begin(), surfaceIDs.end()), surfaceIDs.end());
661 <  for(std::vector<int>::iterator list_iter = surfaceIDs.begin();
662 <      list_iter != surfaceIDs.end(); list_iter++)
601 <    {
602 <    int i = *list_iter;
603 <    surfaceSDs_.push_back(bodydoubles[i]);
659 >    std::cout << surfaceSDs_.size() << std::endl;
660 >    for (SD = surfaceSDs_.begin(); SD != surfaceSDs_.end(); ++SD){
661 >      Vector3d thisatom = (*SD)->getPos();
662 >      std::cout << "Au " << thisatom.x() << "  " << thisatom.y() << " " << thisatom.z() << std::endl;
663      }
664      */
665  
666  
667  
668 <
669 <
670 <
671 <  Ns_ = surfaceSDs_.size();
672 <
673 <
674 <  qh_getarea(qh facet_list);
675 <  volume_ = qh totvol;
676 <  area_ = qh totarea;
677 <
678 <  
679 <
680 <  qh_freeqhull(!qh_ALL);
681 <  qh_memfreeshort(&curlong, &totlong);
682 <  if (curlong || totlong)
683 <    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
684 <              << totlong << curlong << std::endl;
626 <  
627 <  
628 <  // free(pt_array);
629 <  
668 >    Ns_ = surfaceSDs_.size();
669 >    
670 >    
671 >    qh_getarea(qh facet_list);
672 >    volume_ = qh totvol;
673 >    area_ = qh totarea;
674 >    
675 >    
676 >    
677 >    qh_freeqhull(!qh_ALL);
678 >    qh_memfreeshort(&curlong, &totlong);
679 >    if (curlong || totlong)
680 >      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
681 >                << totlong << curlong << std::endl;
682 >    
683 >    
684 >    
685   }
686  
687  

Diff Legend

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