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 1293 by chuckv, Sun Sep 14 01:32:26 2008 UTC vs.
Revision 1307 by chuckv, Mon Oct 20 19:36:32 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.11 2008-10-20 19:36:32 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), nTriangles_(0) {
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");
449 <      painCave.isFatal = 0;
449 >      painCave.isFatal = 1;
450        simError();
451        
452    } //qh_new_qhull
# Line 450 | Line 454 | void ConvexHull::computeHull(std::vector<StuntDouble*>
454  
455   #ifdef IS_MPI
456    std::vector<double> localPts;
457 +  std::vector<double> localVel;
458    int localPtArraySize;
459    
460 +
461    std::fill(isSurfaceID.begin(),isSurfaceID.end(),false);
462 +
463  
464    FORALLfacets {
465      
466      if (!facet->simplicial){
467        // should never happen with Qt
468        sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
469 <      painCave.isFatal = 0;
469 >      painCave.isFatal = 1;
470        simError();
471      }
472      
# Line 484 | Line 491 | void ConvexHull::computeHull(std::vector<StuntDouble*>
491    int localPtArraySize = surfaceIDs.size() * 3;
492    */
493  
494 <  localPts.resize(localPtArraySize);
495 <  //  std::fill(localPts.begin(),globalPts.end(),0.0);
494 >  //localPts.resize(localPtArraySize);
495 >  //std::fill(localPts.begin(),localPts.end(),0.0);
496  
497  
498    int idx = 0;
499 +  int nIsIts = 0;
500 + /*
501    // Copy the surface points into an array.
502    for(std::vector<bool>::iterator list_iter = isSurfaceID.begin();
503        list_iter != isSurfaceID.end(); list_iter++)
# Line 499 | Line 508 | void ConvexHull::computeHull(std::vector<StuntDouble*>
508          localPts.push_back(ptArray[dim_ * idx + 1]);
509          localPts.push_back(ptArray[dim_ * idx + 2]);
510          localPtsMap.push_back(idx);
511 +        nIsIts++;
512        } //Isit
513        idx++;
514      } //isSurfaceID
515 <  
515 >  */
516 >  FORALLvertices {
517 >    idx = qh_pointid(vertex->point);
518 >    localPts.push_back(ptArray[dim_ * idx]);    
519 >    localPts.push_back(ptArray[dim_ * idx + 1]);
520 >    localPts.push_back(ptArray[dim_ * idx + 2]);
521  
522 +    Vector3d vel = bodydoubles[idx]->getVel();
523 +    localVel.push_back(vel.x());
524 +    localVel.push_back(vel.y());
525 +    localVel.push_back(vel.z());
526 +
527 +    localPtsMap.push_back(idx);
528 +  }
529 +
530 +
531    localPtArraySize = localPts.size();
532 +
533 +
534    MPI::COMM_WORLD.Allgather(&localPtArraySize,1,MPI::INT,&NstoProc_[0],1,MPI::INT);
535  
536 +  Nsglobal_=0;
537    for (int i = 0; i < nproc_; i++){
538      Nsglobal_ += NstoProc_[i];
539    }
540 +  
541 +
542 +  int nglobalPts = int(Nsglobal_/3);
543 +
544  
545 <  std::vector<double> globalPts;
546 <  globalPts.resize(Nsglobal_);
547 <  isSurfaceID.resize(int(Nsglobal_/3));
545 >  std::vector<double> globalPts(Nsglobal_);
546 >  std::vector<double> globalVel(Nsglobal_);
547 >
548 >  isSurfaceID.resize(nglobalPts);
549 >
550 >
551    std::fill(globalPts.begin(),globalPts.end(),0.0);
552 +
553 +  displs_[0] = 0;
554    /* Build a displacements array */
555    for (int i = 1; i < nproc_; i++){
556      displs_[i] = displs_[i-1] + NstoProc_[i-1];
557    }
558    
559 +  
560    int noffset = displs_[myrank_];
561    /* gather the potential hull */
562    
563 <  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],NstoProc_,displs_,MPI::DOUBLE);
563 >  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
564 >  MPI::COMM_WORLD.Allgatherv(&localVel[0],localPtArraySize,MPI::DOUBLE,&globalVel[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
565  
566 <
567 <
568 <
566 >  /*
567 >  if (myrank_ == 0){
568 >    for (i = 0; i < globalPts.size(); i++){
569 >      std::cout << globalPts[i] << std::endl;
570 >    }
571 >  }
572 >  */
573    // Free previous hull
574    qh_freeqhull(!qh_ALL);
575    qh_memfreeshort(&curlong, &totlong);
# Line 535 | Line 577 | void ConvexHull::computeHull(std::vector<StuntDouble*>
577      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
578                << totlong << curlong << std::endl;
579  
580 <  if (qh_new_qhull(dim_, numpoints, &globalPts[0], ismalloc,
580 >  if (qh_new_qhull(dim_, nglobalPts, &globalPts[0], ismalloc,
581                      const_cast<char *>(options_.c_str()), NULL, stderr)){
582  
583        sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute global convex hull");
584 <      painCave.isFatal = 0;
584 >      painCave.isFatal = 1;
585        simError();
586        
587    } //qh_new_qhull
# Line 552 | Line 594 | void ConvexHull::computeHull(std::vector<StuntDouble*>
594  
595  
596      unsigned int nf = qh num_facets;
597 <    
597 >    
598      /* Build Surface SD list first */
599  
600      std::fill(isSurfaceID.begin(),isSurfaceID.end(),false);
# Line 562 | Line 604 | void ConvexHull::computeHull(std::vector<StuntDouble*>
604        if (!facet->simplicial){
605        // should never happen with Qt
606          sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
607 <        painCave.isFatal = 0;
607 >        painCave.isFatal = 1;
608          simError();
609        } //simplicical
610        
611        Triangle* face = new Triangle();
612        Vector3d  V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]);
613        face->setNormal(V3dNormal);
614 <      //face->setCentroid(V3dCentroid);
614 >
615 >      
616 >
617        RealType faceArea = 0.5*V3dNormal.length();
618        face->setArea(faceArea);
619  
620  
621        vertices = qh_facet3vertex(facet);
622 +      
623 +      coordT *center = qh_getcenter(vertices);
624 +      Vector3d V3dCentroid(center[0], center[1], center[2]);
625 +      face->setCentroid(V3dCentroid);
626 +      Vector3d faceVel = V3Zero;
627        FOREACHvertex_(vertices){
628          id = qh_pointid(vertex->point);
629 <        face->addVertex(bodydoubles[id]);
630 <        if( !isSurfaceID[id] ){
631 <          isSurfaceID[id] = true;
632 <          surfaceSDs_.push_back(bodydoubles[id]);
633 <        } //IF isSurfaceID
634 <      } //Foreachvertex
629 >        int localindex = id;
630 > #ifdef IS_MPI
631 >        Vector3d velVector(globalVel[dim_ * id],globalVel[dim_ * id + 1], globalVel[dim_ * id + 1]);
632 >        faceVel = faceVel + velVector;
633 >        if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){
634 >          localindex = localPtsMap[id-noffset/3];
635 > #else
636 >          faceVel = faceVel + bodydoubles[localindex]->getVel();
637 > #endif
638 >          face->addVertex(bodydoubles[localindex]);
639 >          if( !isSurfaceID[id] ){
640 >            isSurfaceID[id] = true;
641 > #ifdef IS_MPI      
642 >            
643 > #endif
644 >            
645 >            surfaceSDs_.push_back(bodydoubles[localindex]);
646 >            
647 >          } //IF isSurfaceID
648  
649 <
649 > #ifdef IS_MPI
650 >        
651 >        }else{
652 >          face->addVertex(NULL);
653 >          }
654 > #endif
655 >      } //Foreachvertex
656 >      /*
657 >      if (!SETempty_(facet->coplanarset)){
658 >        FOREACHpoint_(facet->coplanarset){
659 >          id = qh_pointid(point);
660 >          surfaceSDs_.push_back(bodydoubles[id]);
661 >        }
662 >      }
663 >      */
664 >      face->setFacetVelocity(faceVel/3.0);
665        Triangles_.push_back(face);
666        qh_settempfree(&vertices);      
667 <      
667 >
668      } //FORALLfacets
592 –                
669  
594 –                
595 –                
670      /*
671 <  std::sort(surfaceIDs.begin(),surfaceIDs.end());
672 <  surfaceIDs.erase(std::unique(surfaceIDs.begin(), surfaceIDs.end()), surfaceIDs.end());
673 <  for(std::vector<int>::iterator list_iter = surfaceIDs.begin();
674 <      list_iter != surfaceIDs.end(); list_iter++)
601 <    {
602 <    int i = *list_iter;
603 <    surfaceSDs_.push_back(bodydoubles[i]);
671 >    std::cout << surfaceSDs_.size() << std::endl;
672 >    for (SD = surfaceSDs_.begin(); SD != surfaceSDs_.end(); ++SD){
673 >      Vector3d thisatom = (*SD)->getPos();
674 >      std::cout << "Au " << thisatom.x() << "  " << thisatom.y() << " " << thisatom.z() << std::endl;
675      }
676      */
677  
678  
679  
680 <
681 <
682 <
683 <  Ns_ = surfaceSDs_.size();
684 <
685 <
686 <  qh_getarea(qh facet_list);
687 <  volume_ = qh totvol;
688 <  area_ = qh totarea;
689 <
690 <  
691 <
692 <  qh_freeqhull(!qh_ALL);
693 <  qh_memfreeshort(&curlong, &totlong);
694 <  if (curlong || totlong)
695 <    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
696 <              << totlong << curlong << std::endl;
626 <  
627 <  
628 <  // free(pt_array);
629 <  
680 >    Ns_ = surfaceSDs_.size();
681 >    nTriangles_ = Triangles_.size();
682 >    
683 >    qh_getarea(qh facet_list);
684 >    volume_ = qh totvol;
685 >    area_ = qh totarea;
686 >    
687 >    
688 >    
689 >    qh_freeqhull(!qh_ALL);
690 >    qh_memfreeshort(&curlong, &totlong);
691 >    if (curlong || totlong)
692 >      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
693 >                << totlong << curlong << std::endl;
694 >    
695 >    
696 >    
697   }
698  
699  

Diff Legend

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