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 1301 by chuckv, Sun Sep 14 01:32:26 2008 UTC vs.
Revision 1302 by chuckv, Tue Oct 7 17:12:48 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.9 2008-10-07 17:12:48 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 > ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt  Qci Tcv Pp"), Ns_(200) {
373    //If we are doing the mpi version, set up some vectors for data communication
374   #ifdef IS_MPI
375  
# Line 380 | Line 380 | ConvexHull::ConvexHull() : Hull(), dim_(3), options_("
380   displs_   = new int[nproc_];
381  
382   // Create a surface point type in MPI to send
383 < surfacePtType = MPI::DOUBLE.Create_contiguous(3);
384 < surfacePtType.Commit();
383 > //surfacePtType = MPI::DOUBLE.Create_contiguous(3);
384 > // surfacePtType.Commit();
385  
386  
387   #endif
# Line 413 | Line 413 | void ConvexHull::computeHull(std::vector<StuntDouble*>
413    
414    //pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_));
415  
416 <  double* ptArray = new double[numpoints * 3];
416 > //  double* ptArray = new double[numpoints * 3];
417 >  std::vector<double> ptArray(numpoints*3);
418    std::vector<bool> isSurfaceID(numpoints);
419  
420    // Copy the positon vector into a points vector for qhull.
# Line 435 | Line 436 | void ConvexHull::computeHull(std::vector<StuntDouble*>
436    
437    
438    boolT ismalloc = False;
439 +  /* Clean up memory from previous convex hull calculations*/
440    Triangles_.clear();
441    surfaceSDs_.clear();
442 <  
443 <  if (qh_new_qhull(dim_, numpoints, ptArray, ismalloc,
442 >  surfaceSDs_.reserve(Ns_);
443 >
444 >  if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc,
445                      const_cast<char *>(options_.c_str()), NULL, stderr)) {
446  
447        sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull");
# Line 452 | Line 455 | void ConvexHull::computeHull(std::vector<StuntDouble*>
455    std::vector<double> localPts;
456    int localPtArraySize;
457    
458 +
459    std::fill(isSurfaceID.begin(),isSurfaceID.end(),false);
460 +
461  
462    FORALLfacets {
463      
# Line 484 | Line 489 | void ConvexHull::computeHull(std::vector<StuntDouble*>
489    int localPtArraySize = surfaceIDs.size() * 3;
490    */
491  
492 <  localPts.resize(localPtArraySize);
493 <  //  std::fill(localPts.begin(),globalPts.end(),0.0);
492 >  //localPts.resize(localPtArraySize);
493 >  //std::fill(localPts.begin(),localPts.end(),0.0);
494  
495  
496    int idx = 0;
497 +  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++)
# Line 499 | Line 506 | void ConvexHull::computeHull(std::vector<StuntDouble*>
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 <  
513 >  */
514 >  FORALLvertices {
515 >    idx = qh_pointid(vertex->point);
516 >    localPts.push_back(ptArray[dim_ * idx]);    
517 >    localPts.push_back(ptArray[dim_ * idx + 1]);
518 >    localPts.push_back(ptArray[dim_ * idx + 2]);
519 >    localPtsMap.push_back(idx);
520 >  }
521  
522 +
523    localPtArraySize = localPts.size();
524 +
525 +
526    MPI::COMM_WORLD.Allgather(&localPtArraySize,1,MPI::INT,&NstoProc_[0],1,MPI::INT);
527  
528 +  Nsglobal_=0;
529    for (int i = 0; i < nproc_; i++){
530      Nsglobal_ += NstoProc_[i];
531    }
532 +  
533 +
534 +  int nglobalPts = int(Nsglobal_/3);
535 +
536  
537    std::vector<double> globalPts;
538    globalPts.resize(Nsglobal_);
539 <  isSurfaceID.resize(int(Nsglobal_/3));
539 >
540 >  isSurfaceID.resize(nglobalPts);
541 >
542 >
543    std::fill(globalPts.begin(),globalPts.end(),0.0);
544 +
545 +  displs_[0] = 0;
546    /* Build a displacements array */
547    for (int i = 1; i < nproc_; i++){
548      displs_[i] = displs_[i-1] + NstoProc_[i-1];
549    }
550    
551 +  
552    int noffset = displs_[myrank_];
553    /* gather the potential hull */
554    
555 <  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],NstoProc_,displs_,MPI::DOUBLE);
555 >  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE);
556  
557 <
558 <
559 <
557 >  /*
558 >  if (myrank_ == 0){
559 >    for (i = 0; i < globalPts.size(); i++){
560 >      std::cout << globalPts[i] << std::endl;
561 >    }
562 >  }
563 >  */
564    // Free previous hull
565    qh_freeqhull(!qh_ALL);
566    qh_memfreeshort(&curlong, &totlong);
# Line 535 | Line 568 | void ConvexHull::computeHull(std::vector<StuntDouble*>
568      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
569                << totlong << curlong << std::endl;
570  
571 <  if (qh_new_qhull(dim_, numpoints, &globalPts[0], ismalloc,
571 >  if (qh_new_qhull(dim_, nglobalPts, &globalPts[0], ismalloc,
572                      const_cast<char *>(options_.c_str()), NULL, stderr)){
573  
574        sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute global convex hull");
575 <      painCave.isFatal = 0;
575 >      painCave.isFatal = 1;
576        simError();
577        
578    } //qh_new_qhull
# Line 552 | Line 585 | void ConvexHull::computeHull(std::vector<StuntDouble*>
585  
586  
587      unsigned int nf = qh num_facets;
588 <    
588 >    
589      /* Build Surface SD list first */
590  
591      std::fill(isSurfaceID.begin(),isSurfaceID.end(),false);
# Line 562 | Line 595 | void ConvexHull::computeHull(std::vector<StuntDouble*>
595        if (!facet->simplicial){
596        // should never happen with Qt
597          sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
598 <        painCave.isFatal = 0;
598 >        painCave.isFatal = 1;
599          simError();
600        } //simplicical
601        
# Line 571 | Line 604 | void ConvexHull::computeHull(std::vector<StuntDouble*>
604        face->setNormal(V3dNormal);
605        //face->setCentroid(V3dCentroid);
606        RealType faceArea = 0.5*V3dNormal.length();
607 <      face->setArea(faceArea);
607 >      //face->setArea(faceArea);
608  
609  
610        vertices = qh_facet3vertex(facet);
611        FOREACHvertex_(vertices){
612          id = qh_pointid(vertex->point);
613 <        face->addVertex(bodydoubles[id]);
614 <        if( !isSurfaceID[id] ){
615 <          isSurfaceID[id] = true;
616 <          surfaceSDs_.push_back(bodydoubles[id]);
617 <        } //IF isSurfaceID
613 >        int localindex = id;
614 > #ifdef IS_MPI
615 >        
616 >        if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){
617 >          localindex = localPtsMap[id-noffset/3];
618 > #endif
619 >          face->addVertex(bodydoubles[localindex]);
620 >          if( !isSurfaceID[id] ){
621 >            isSurfaceID[id] = true;
622 > #ifdef IS_MPI      
623 >            
624 > #endif
625 >            
626 >            surfaceSDs_.push_back(bodydoubles[localindex]);
627 >            
628 >          } //IF isSurfaceID
629 >
630 > #ifdef IS_MPI
631 >        
632 >        }else{
633 >          face->addVertex(NULL);
634 >          }
635 > #endif
636        } //Foreachvertex
637  
587 –
638        Triangles_.push_back(face);
639        qh_settempfree(&vertices);      
640        
641      } //FORALLfacets
592 –                
642  
643 <                
595 <                
596 <    /*
597 <  std::sort(surfaceIDs.begin(),surfaceIDs.end());
598 <  surfaceIDs.erase(std::unique(surfaceIDs.begin(), surfaceIDs.end()), surfaceIDs.end());
599 <  for(std::vector<int>::iterator list_iter = surfaceIDs.begin();
600 <      list_iter != surfaceIDs.end(); list_iter++)
601 <    {
602 <    int i = *list_iter;
603 <    surfaceSDs_.push_back(bodydoubles[i]);
604 <    }
605 <    */
643 >                        
644  
607 –
608 –
609 –
610 –
611 –
645    Ns_ = surfaceSDs_.size();
646  
647  
# Line 624 | Line 657 | void ConvexHull::computeHull(std::vector<StuntDouble*>
657      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
658                << totlong << curlong << std::endl;
659    
660 +
661    
628 –  // free(pt_array);
629 –  
662   }
663  
664  

Diff Legend

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