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 1316 by chuckv, Fri Nov 14 15:44:34 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.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 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 377 | 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);
387 < surfacePtType.Commit();
386 > //surfacePtType = MPI::DOUBLE.Create_contiguous(3);
387 > // surfacePtType.Commit();
388  
389  
390   #endif
# Line 413 | Line 416 | void ConvexHull::computeHull(std::vector<StuntDouble*>
416    
417    //pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_));
418  
419 <  double* ptArray = new double[numpoints * 3];
419 > //  double* ptArray = new double[numpoints * 3];
420 >  std::vector<double> ptArray(numpoints*3);
421    std::vector<bool> isSurfaceID(numpoints);
422  
423    // Copy the positon vector into a points vector for qhull.
# Line 435 | Line 439 | void ConvexHull::computeHull(std::vector<StuntDouble*>
439    
440    
441    boolT ismalloc = False;
442 +  /* Clean up memory from previous convex hull calculations*/
443 +  
444    Triangles_.clear();
445    surfaceSDs_.clear();
446 <  
447 <  if (qh_new_qhull(dim_, numpoints, ptArray, ismalloc,
446 >  surfaceSDs_.reserve(Ns_);
447 >
448 >  if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc,
449                      const_cast<char *>(options_.c_str()), NULL, stderr)) {
450  
451        sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull");
452 <      painCave.isFatal = 0;
452 >      painCave.isFatal = 1;
453        simError();
454        
455    } //qh_new_qhull
# Line 450 | Line 457 | void ConvexHull::computeHull(std::vector<StuntDouble*>
457  
458   #ifdef IS_MPI
459    std::vector<double> localPts;
460 +  std::vector<double> localVel;
461 +  std::vector<double> localMass;
462    int localPtArraySize;
463    
464 +
465    std::fill(isSurfaceID.begin(),isSurfaceID.end(),false);
466 +
467  
468    FORALLfacets {
469      
470      if (!facet->simplicial){
471        // should never happen with Qt
472        sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
473 <      painCave.isFatal = 0;
473 >      painCave.isFatal = 1;
474        simError();
475      }
476      
# Line 478 | Line 489 | void ConvexHull::computeHull(std::vector<StuntDouble*>
489  
490  
491  
481 –  /*
482 –  std::sort(surfaceIDs.begin(),surfaceIDs.end());
483 –  surfaceIDs.erase(std::unique(surfaceIDs.begin(), surfaceIDs.end()), surfaceIDs.end());
484 –  int localPtArraySize = surfaceIDs.size() * 3;
485 –  */
492  
487 –  localPts.resize(localPtArraySize);
488 –  //  std::fill(localPts.begin(),globalPts.end(),0.0);
489 –
490 –
493    int idx = 0;
494 <  // Copy the surface points into an array.
495 <  for(std::vector<bool>::iterator list_iter = isSurfaceID.begin();
496 <      list_iter != isSurfaceID.end(); list_iter++)
497 <    {
498 <      bool isIt = *list_iter;
499 <      if (isIt){
498 <        localPts.push_back(ptArray[dim_ * idx]);    
499 <        localPts.push_back(ptArray[dim_ * idx + 1]);
500 <        localPts.push_back(ptArray[dim_ * idx + 2]);
501 <        localPtsMap.push_back(idx);
502 <      } //Isit
503 <      idx++;
504 <    } //isSurfaceID
505 <  
506 <
507 <  localPtArraySize = localPts.size();
508 <  MPI::COMM_WORLD.Allgather(&localPtArraySize,1,MPI::INT,&NstoProc_[0],1,MPI::INT);
494 >  int nIsIts = 0;
495 >  FORALLvertices {
496 >    idx = qh_pointid(vertex->point);
497 >    localPts.push_back(ptArray[dim_ * idx]);    
498 >    localPts.push_back(ptArray[dim_ * idx + 1]);
499 >    localPts.push_back(ptArray[dim_ * idx + 2]);
500  
501 +    Vector3d vel = bodydoubles[idx]->getVel();
502 +    localVel.push_back(vel.x());
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 +
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 +  
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 = Nsglobal_*3;
529 +
530  
531 <  std::vector<double> globalPts;
532 <  globalPts.resize(Nsglobal_);
533 <  isSurfaceID.resize(int(Nsglobal_/3));
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 +  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 <  int noffset = displs_[myrank_];
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 = vecdispls_[myrank_];
552    /* gather the potential hull */
553    
554 <  MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],NstoProc_,displs_,MPI::DOUBLE);
555 <
556 <
557 <
558 <
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++){
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        
602 <      Triangle* face = new Triangle();
602 >      Triangle face;
603        Vector3d  V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]);
604 <      face->setNormal(V3dNormal);
605 <      //face->setCentroid(V3dCentroid);
606 <      RealType faceArea = 0.5*V3dNormal.length();
574 <      face->setArea(faceArea);
604 >      face.setNormal(V3dNormal);
605 >
606 >      
607  
608 +      //RealType faceArea = 0.5*V3dNormal.length();
609 +      RealType faceArea = qh_facetarea(facet);
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 +      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 <        face->addVertex(bodydoubles[id]);
625 <        if( !isSurfaceID[id] ){
626 <          isSurfaceID[id] = true;
627 <          surfaceSDs_.push_back(bodydoubles[id]);
628 <        } //IF isSurfaceID
629 <      } //Foreachvertex
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.addVertexSD(bodydoubles[localindex]);
640 >          if( !isSurfaceID[id] ){
641 >            isSurfaceID[id] = true;
642 > #ifdef IS_MPI      
643 >            
644 > #endif
645 >            
646 >            surfaceSDs_.push_back(bodydoubles[localindex]);
647 >            
648 >          } //IF isSurfaceID
649  
650 <
650 > #ifdef IS_MPI
651 >        
652 >        }else{
653 >          face.addVertexSD(NULL);
654 >          }
655 > #endif
656 >        ver++;
657 >      } //Foreachvertex
658 >      /*
659 >      if (!SETempty_(facet->coplanarset)){
660 >        FOREACHpoint_(facet->coplanarset){
661 >          id = qh_pointid(point);
662 >          surfaceSDs_.push_back(bodydoubles[id]);
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);      
671 <      
671 >
672      } //FORALLfacets
592 –                
673  
594 –                
595 –                
674      /*
675 <  std::sort(surfaceIDs.begin(),surfaceIDs.end());
676 <  surfaceIDs.erase(std::unique(surfaceIDs.begin(), surfaceIDs.end()), surfaceIDs.end());
677 <  for(std::vector<int>::iterator list_iter = surfaceIDs.begin();
678 <      list_iter != surfaceIDs.end(); list_iter++)
601 <    {
602 <    int i = *list_iter;
603 <    surfaceSDs_.push_back(bodydoubles[i]);
675 >    std::cout << surfaceSDs_.size() << std::endl;
676 >    for (SD = surfaceSDs_.begin(); SD != surfaceSDs_.end(); ++SD){
677 >      Vector3d thisatom = (*SD)->getPos();
678 >      std::cout << "Au " << thisatom.x() << "  " << thisatom.y() << " " << thisatom.z() << std::endl;
679      }
680      */
681  
682  
683  
684 <
685 <
686 <
687 <  Ns_ = surfaceSDs_.size();
688 <
689 <
690 <  qh_getarea(qh facet_list);
691 <  volume_ = qh totvol;
692 <  area_ = qh totarea;
693 <
694 <  
695 <
696 <  qh_freeqhull(!qh_ALL);
697 <  qh_memfreeshort(&curlong, &totlong);
698 <  if (curlong || totlong)
699 <    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
700 <              << totlong << curlong << std::endl;
626 <  
627 <  
628 <  // free(pt_array);
629 <  
684 >    Ns_ = surfaceSDs_.size();
685 >    nTriangles_ = Triangles_.size();
686 >    
687 >    qh_getarea(qh facet_list);
688 >    volume_ = qh totvol;
689 >    area_ = qh totarea;
690 >    
691 >    
692 >    
693 >    qh_freeqhull(!qh_ALL);
694 >    qh_memfreeshort(&curlong, &totlong);
695 >    if (curlong || totlong)
696 >      std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
697 >                << totlong << curlong << std::endl;
698 >    
699 >    
700 >    
701   }
702  
703  

Diff Legend

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