src/math/AlphaHull.cpp

/* Copyright (c) 2008, 2009, 2010 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
 *
 *  AlphaHull.cpp
 *
 *  Purpose: To calculate Alpha hull, hull volume libqhull.
 *
 *  Created by Charles F. Vardeman II on 11 Dec 2006.
 *  @author  Charles F. Vardeman II
 *  @version $Id$
 *
 */

/* Standard includes independent of library */

#include <iostream>
#include <fstream>
#include <list>
#include <algorithm>
#include <iterator>
#include <utility>
#include "math/AlphaHull.hpp"
#include "utils/simError.h"

#ifdef IS_MPI
#include <mpi.h>
#endif

using namespace OpenMD;

#ifdef HAVE_QHULL
extern "C"
{
#include <qhull/libqhull.h>
#include <qhull/mem.h>
#include <qhull/qset.h>
#include <qhull/geom.h>
#include <qhull/merge.h>
#include <qhull/poly.h>
#include <qhull/io.h>
#include <qhull/stat.h>
}
double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim);

AlphaHull::AlphaHull(double alpha) : Hull(), dim_(3), alpha_(alpha), options_("qhull d QJ Tcv Pp") {
}

void AlphaHull::computeHull(std::vector<StuntDouble*> bodydoubles) { 
 
  int numpoints = bodydoubles.size();
  bool alphashape=true;
  
  Triangles_.clear();
  
  vertexT *vertex, **vertexp;
  facetT *facet, *neighbor;
  setT *vertices, *verticestop, *verticesbottom;
  int curlong, totlong;
  pointT *interiorPoint;
  
  std::vector<double> ptArray(numpoints*dim_);

  // Copy the positon vector into a points vector for qhull.
  std::vector<StuntDouble*>::iterator SD;
  int i = 0;
  for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){
    Vector3d pos = (*SD)->getPos();      
    ptArray[dim_ * i] = pos.x();
    ptArray[dim_ * i + 1] = pos.y();
    ptArray[dim_ * i + 2] = pos.z();
    i++;
  }

    /* Clean up memory from previous convex hull calculations*/
  boolT ismalloc = False;

  int ridgesCount=0;
  if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc,
                   const_cast<char *>(options_.c_str()), NULL, stderr)) {

    sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute convex hull");
    painCave.isFatal = 1;
    simError();
    
  } //qh_new_qhull


#ifdef IS_MPI
  //If we are doing the mpi version, set up some vectors for data communication
  
  int nproc = MPI::COMM_WORLD.Get_size();
  int myrank = MPI::COMM_WORLD.Get_rank();
  int localHullSites = 0;

  std::vector<int> hullSitesOnProc(nproc, 0);
  std::vector<int> coordsOnProc(nproc, 0);
  std::vector<int> displacements(nproc, 0);
  std::vector<int> vectorDisplacements(nproc, 0);

  std::vector<double> coords;
  std::vector<double> vels;
  std::vector<int> indexMap;
  std::vector<double> masses;

  FORALLvertices{
    localHullSites++;
    
    int idx = qh_pointid(vertex->point);

    indexMap.push_back(idx);

    coords.push_back(ptArray[dim_  * idx]);
    coords.push_back(ptArray[dim_  * idx + 1]);
    coords.push_back(ptArray[dim_  * idx + 2]);

    StuntDouble* sd = bodydoubles[idx];

    Vector3d vel = sd->getVel();
    vels.push_back(vel.x());
    vels.push_back(vel.y());
    vels.push_back(vel.z());

    masses.push_back(sd->getMass());
  }

  MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0],
                            1, MPI::INT);

  int globalHullSites = 0;
  for (int iproc = 0; iproc < nproc; iproc++){
    globalHullSites += hullSitesOnProc[iproc];
    coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc];
  }

  displacements[0] = 0;
  vectorDisplacements[0] = 0;
  
  for (int iproc = 1; iproc < nproc; iproc++){
    displacements[iproc] = displacements[iproc-1] + hullSitesOnProc[iproc-1];
    vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1]; 
  }

  std::vector<double> globalCoords(dim_ * globalHullSites);
  std::vector<double> globalVels(dim_ * globalHullSites);
  std::vector<double> globalMasses(globalHullSites);

  int count = coordsOnProc[myrank];
  
  MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0],
                             &coordsOnProc[0], &vectorDisplacements[0], 
                             MPI::DOUBLE);

  MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0], 
                             &coordsOnProc[0], &vectorDisplacements[0],
                             MPI::DOUBLE);

  MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE,
                             &globalMasses[0], &hullSitesOnProc[0], 
                             &displacements[0], MPI::DOUBLE);

  // Free previous hull
  qh_freeqhull(!qh_ALL);
  qh_memfreeshort(&curlong, &totlong);
  if (curlong || totlong)
    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
              << totlong << curlong << std::endl;
  
  if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc,
                   const_cast<char *>(options_.c_str()), NULL, stderr)){
    
    sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute global convex hull");
    painCave.isFatal = 1;
    simError();
    
  } //qh_new_qhull


#endif

  //Set facet->center as the Voronoi center
  qh_setvoronoi_all();
  
  
  int convexNumVert = qh_setsize(qh_facetvertices (qh facet_list, NULL, false));
  //Insert all the sample points, because, even with alpha=0, the alpha shape/alpha complex will
  //contain them.

  //  tri::Allocator<CMeshO>::AddVertices(pm.cm,convexNumVert);
  
  /*ivp length is 'qh num_vertices' because each vertex is accessed through its ID whose range is 
    0<=qh_pointid(vertex->point)<qh num_vertices*/
  //  vector<tri::Allocator<CMeshO>::VertexPointer> ivp(qh num_vertices);
  /*i=0;
  FORALLvertices{       
    if ((*vertex).point){
      //  pm.cm.vert[i].P()[0] = (*vertex).point[0];
      // pm.cm.vert[i].P()[1] = (*vertex).point[1];
      //pm.cm.vert[i].P()[2] = (*vertex).point[2];
      // ivp[qh_pointid(vertex->point)] = &pm.cm.vert[i];
      i++;
    }
  }
  */
  //Set of alpha complex triangles for alphashape filtering
  setT* set= qh_settemp(4* qh num_facets); 
  
  qh visit_id++;
  int numFacets=0;
  std::vector<std::vector <int> > facetlist;
  interiorPoint = qh interior_point;
  FORALLfacet_(qh facet_list) {
    numFacets++;
    if (!facet->upperdelaunay) {
      //For all facets (that are tetrahedrons)calculate the radius of the empty circumsphere considering 
      //the distance between the circumcenter and a vertex of the facet
      vertexT* vertex = (vertexT *)(facet->vertices->e[0].p);
      double* center = facet->center;
      double radius =  qh_pointdist(vertex->point,center,dim_);
      
      if (radius>alpha_) // if the facet is not good consider the ridges
        {
          //if calculating the alphashape, unmark the facet ('good' is used as 'marked'). 
          facet->good=false;
          
          //Compute each ridge (triangle) once and test the cironference radius with alpha
          facet->visitid= qh visit_id;
          qh_makeridges(facet);
          ridgeT *ridge, **ridgep;
          int goodTriangles=0;
          FOREACHridge_(facet->ridges) {
            neighbor= otherfacet_(ridge, facet);
            if (( neighbor->visitid != qh visit_id)){                   
              //Calculate the radius of the circumference 
              pointT* p0 = ((vertexT*) (ridge->vertices->e[0].p))->point;
              pointT* p1 = ((vertexT*) (ridge->vertices->e[1].p))->point;
              pointT* p2 = ((vertexT*) (ridge->vertices->e[2].p))->point;
              
              radius = calculate_circumradius(p0,p1,p2, dim_);
              
              if(radius <=alpha_){
                goodTriangles++;
                //save the triangle (ridge) for subsequent filtering
                qh_setappend(&set, ridge); 
              }
            }
          }

          //If calculating the alphashape, mark the facet('good' is used as 'marked'). 
          //This facet will have some triangles hidden by the facet's neighbor.
          if(goodTriangles==4)
            facet->good=true;
          
        }
      else //the facet is good. Put all the triangles of the tetrahedron in the mesh
        {
          //Compute each ridge (triangle) once
          facet->visitid= qh visit_id;
          //If calculating the alphashape, mark the facet('good' is used as 'marked').
          //This facet will have some triangles hidden by the facet's neighbor.
          facet->good=true;
          qh_makeridges(facet);
          ridgeT *ridge, **ridgep;
          FOREACHridge_(facet->ridges) {
            neighbor= otherfacet_(ridge, facet);
            if ((neighbor->visitid != qh visit_id)){
                 qh_setappend(&set, ridge);
            }   
          }
        }
    }
  }
  //assert(numFacets== qh num_facets);
  
  //Filter the triangles (only the ones on the boundary of the alpha complex) and build the mesh


  ridgeT *ridge, **ridgep;
  FOREACHridge_(set) {
    if ((!ridge->top->good || !ridge->bottom->good || ridge->top->upperdelaunay || ridge->bottom->upperdelaunay)){
      //        tri::Allocator<CMeshO>::FaceIterator fi=tri::Allocator<CMeshO>::AddFaces(pm.cm,1);
      ridgesCount++;
      int vertex_n, vertex_i;
      Triangle face;
      
      // Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]);
      //face.setNormal(V3dNormal);
      
      
      //coordT *center = qh_getcenter(ridge->vertices);
      //std::cout << "Centers are " << center[0] << "  " <<center[1] << "  " << center[2] << std::endl;
      //Vector3d V3dCentroid(center[0], center[1], center[2]);
      //face.setCentroid(V3dCentroid);

      
      Vector3d faceVel = V3Zero;
      Vector3d p[3];
      RealType faceMass = 0.0;
      
      int ver = 0;
      std::vector<int> virtexlist;
      FOREACHvertex_i_(ridge->vertices){
        int id = qh_pointid(vertex->point);
        p[ver][0] = vertex->point[0];
        p[ver][1] = vertex->point[1];
        p[ver][2] = vertex->point[2];
        Vector3d vel;
        RealType mass;
        ver++;
        virtexlist.push_back(id);
        // std::cout << "Ridge: " << ridgesCount << " Vertex " << id << std::endl; 

        vel = bodydoubles[id]->getVel();
        mass = bodydoubles[id]->getMass();
        face.addVertexSD(bodydoubles[id]);


        faceVel = faceVel + vel;
        faceMass = faceMass + mass;
      } //FOREACH Vertex 
      facetlist.push_back(virtexlist);
      face.addVertices(p[0],p[1],p[2]);
      face.setFacetMass(faceMass);
      face.setFacetVelocity(faceVel/3.0);
      
      RealType area = face.getArea();
      area_ += area;
      Vector3d normal = face.getUnitNormal();
      RealType offset =  ((0.0-p[0][0])*normal[0] + (0.0-p[0][1])*normal[1] + (0.0-p[0][2])*normal[2]);
      RealType dist =  normal[0] * interiorPoint[0] + normal[1]*interiorPoint[1] + normal[2]*interiorPoint[2];
      std::cout << "Dist and normal and area are: " << normal << std::endl;
      volume_ += dist *area/qh hull_dim;
      
      Triangles_.push_back(face);
    }
  }

  std::cout << "Volume is: " << volume_ << std::endl; 

//assert(pm.cm.fn == ridgesCount);
/*
  std::cout <<"OFF"<<std::endl; 
  std::cout << bodydoubles.size() << "  " << facetlist.size() << "  " << 3*facetlist.size() << std::endl; 
  for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){
    Vector3d pos = (*SD)->getPos();      
    std::cout << pos.x() << "  " << pos.y() << "  " << pos.z() << std::endl; 
  }

  
  std::vector<std::vector<int> >::iterator thisfacet;
  std::vector<int>::iterator thisvertex;

  for (thisfacet = facetlist.begin(); thisfacet != facetlist.end(); thisfacet++){
    std::cout << (*thisfacet).size(); 
    for (thisvertex = (*thisfacet).begin(); thisvertex != (*thisfacet).end(); thisvertex++){
      std::cout << "  " <<  *thisvertex;
    }
    std::cout << std::endl; 
  }
*/


/*  
  FORALLfacets {  
    Triangle face;

    Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]);
    face.setNormal(V3dNormal);
    
    RealType faceArea = qh_facetarea(facet);
    face.setArea(faceArea);
    
    vertices = qh_facet3vertex(facet);
      
    coordT *center = qh_getcenter(vertices);
    Vector3d V3dCentroid(center[0], center[1], center[2]);
    face.setCentroid(V3dCentroid);

    Vector3d faceVel = V3Zero;
    Vector3d p[3];
    RealType faceMass = 0.0;

    int ver = 0;

    FOREACHvertex_(vertices){
      int id = qh_pointid(vertex->point);
      p[ver][0] = vertex->point[0];
      p[ver][1] = vertex->point[1];
      p[ver][2] = vertex->point[2];
      
      Vector3d vel;
      RealType mass;

#ifdef IS_MPI
      vel = Vector3d(globalVels[dim_ * id],
                     globalVels[dim_ * id + 1], 
                     globalVels[dim_ * id + 2]);
      mass = globalMasses[id];

      // localID will be between 0 and hullSitesOnProc[myrank] if we
      // own this guy.

      int localID = id - displacements[myrank];

      if (localID >= 0 && localID < hullSitesOnProc[myrank])
        face.addVertexSD(bodydoubles[indexMap[localID]]);
      
#else
      vel = bodydoubles[id]->getVel();
      mass = bodydoubles[id]->getMass();
      face.addVertexSD(bodydoubles[id]);      
#endif
        
      faceVel = faceVel + vel;
      faceMass = faceMass + mass;
      ver++;      
    } //Foreachvertex

    face.addVertices(p[0], p[1], p[2]);
    face.setFacetMass(faceMass);
    face.setFacetVelocity(faceVel/3.0);
    Triangles_.push_back(face);
    qh_settempfree(&vertices);      

  } //FORALLfacets
*/
  // qh_getarea(qh facet_list);
  //volume_ = qh totvol;
  // area_ = qh totarea;

  qh_freeqhull(!qh_ALL);
  qh_memfreeshort(&curlong, &totlong);
  if (curlong || totlong)
    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
              << totlong << curlong << std::endl;    
}

void AlphaHull::printHull(const std::string& geomFileName) {

#ifdef IS_MPI
  if (worldRank == 0)  {
#endif
  FILE *newGeomFile;
  
  //create new .md file based on old .md file
  newGeomFile = fopen(geomFileName.c_str(), "w");
  qh_findgood_all(qh facet_list);
  for (int i = 0; i < qh_PRINTEND; i++)
    qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
  
  fclose(newGeomFile);
#ifdef IS_MPI
  }
#endif  
}

  double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim){
    coordT a = qh_pointdist(p0,p1,dim);
    coordT b = qh_pointdist(p1,p2,dim);
    coordT c = qh_pointdist(p2,p0,dim);
    
    coordT sum =(a + b + c)*0.5;
    coordT area = sum*(a+b-sum)*(a+c-sum)*(b+c-sum);
    return (double) (a*b*c)/(4*sqrt(area));
  }

#endif //QHULL
Revision:	1665
Committed:	Tue Nov 22 20:38:56 2011 UTC (13 years, 5 months ago) by gezelter
File size:	16618 byte(s)
Log Message:	updated copyright notices
#	User	Rev	Content
1	chuckv	1402	/* Copyright (c) 2008, 2009, 2010 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
5			* redistribute this software in source and binary code form, provided
6			* that the following conditions are met:
7			*
8			* 1. Redistributions of source code must retain the above copyright
9			* notice, this list of conditions and the following disclaimer.
10			*
11			* 2. Redistributions in binary form must reproduce the above copyright
12			* notice, this list of conditions and the following disclaimer in the
13			* documentation and/or other materials provided with the
14			* distribution.
15			*
16			* This software is provided "AS IS," without a warranty of any
17			* kind. All express or implied conditions, representations and
18			* warranties, including any implied warranty of merchantability,
19			* fitness for a particular purpose or non-infringement, are hereby
20			* excluded. The University of Notre Dame and its licensors shall not
21			* be liable for any damages suffered by licensee as a result of
22			* using, modifying or distributing the software or its
23			* derivatives. In no event will the University of Notre Dame or its
24			* licensors be liable for any lost revenue, profit or data, or for
25			* direct, indirect, special, consequential, incidental or punitive
26			* damages, however caused and regardless of the theory of liability,
27			* arising out of the use of or inability to use software, even if the
28			* University of Notre Dame has been advised of the possibility of
29			* such damages.
30			*
31			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32			* research, please cite the appropriate papers when you publish your
33			* work. Good starting points are:
34			*
35			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
38	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
39			* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
40	chuckv	1402	*
41			* AlphaHull.cpp
42			*
43			* Purpose: To calculate Alpha hull, hull volume libqhull.
44			*
45			* Created by Charles F. Vardeman II on 11 Dec 2006.
46			* @author Charles F. Vardeman II
47	gezelter	1442	* @version $Id$
48	chuckv	1402	*
49			*/
50
51			/* Standard includes independent of library */
52
53			#include <iostream>
54			#include <fstream>
55			#include <list>
56			#include <algorithm>
57			#include <iterator>
58			#include <utility>
59			#include "math/AlphaHull.hpp"
60			#include "utils/simError.h"
61
62			#ifdef IS_MPI
63			#include <mpi.h>
64			#endif
65
66			using namespace OpenMD;
67
68			#ifdef HAVE_QHULL
69			extern "C"
70			{
71	gezelter	1618	#include <qhull/libqhull.h>
72	chuckv	1402	#include <qhull/mem.h>
73			#include <qhull/qset.h>
74			#include <qhull/geom.h>
75			#include <qhull/merge.h>
76			#include <qhull/poly.h>
77			#include <qhull/io.h>
78			#include <qhull/stat.h>
79			}
80			double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim);
81
82	chuckv	1404	AlphaHull::AlphaHull(double alpha) : Hull(), dim_(3), alpha_(alpha), options_("qhull d QJ Tcv Pp") {
83	chuckv	1402	}
84
85			void AlphaHull::computeHull(std::vector<StuntDouble*> bodydoubles) {
86
87			int numpoints = bodydoubles.size();
88			bool alphashape=true;
89
90			Triangles_.clear();
91
92			vertexT vertex, *vertexp;
93			facetT facet, neighbor;
94	chuckv	1404	setT vertices, verticestop, *verticesbottom;
95	chuckv	1402	int curlong, totlong;
96	chuckv	1404	pointT *interiorPoint;
97	chuckv	1402
98			std::vector<double> ptArray(numpoints*dim_);
99
100			// Copy the positon vector into a points vector for qhull.
101			std::vector<StuntDouble*>::iterator SD;
102			int i = 0;
103			for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){
104			Vector3d pos = (*SD)->getPos();
105			ptArray[dim_ * i] = pos.x();
106			ptArray[dim_ * i + 1] = pos.y();
107			ptArray[dim_ * i + 2] = pos.z();
108			i++;
109			}
110
111			/* Clean up memory from previous convex hull calculations*/
112			boolT ismalloc = False;
113
114			int ridgesCount=0;
115			if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc,
116			const_cast<char *>(options_.c_str()), NULL, stderr)) {
117
118			sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute convex hull");
119			painCave.isFatal = 1;
120			simError();
121
122			} //qh_new_qhull
123
124
125			#ifdef IS_MPI
126			//If we are doing the mpi version, set up some vectors for data communication
127
128			int nproc = MPI::COMM_WORLD.Get_size();
129			int myrank = MPI::COMM_WORLD.Get_rank();
130			int localHullSites = 0;
131
132			std::vector<int> hullSitesOnProc(nproc, 0);
133			std::vector<int> coordsOnProc(nproc, 0);
134			std::vector<int> displacements(nproc, 0);
135			std::vector<int> vectorDisplacements(nproc, 0);
136
137			std::vector<double> coords;
138			std::vector<double> vels;
139			std::vector<int> indexMap;
140			std::vector<double> masses;
141
142			FORALLvertices{
143			localHullSites++;
144
145			int idx = qh_pointid(vertex->point);
146
147			indexMap.push_back(idx);
148
149			coords.push_back(ptArray[dim_ * idx]);
150			coords.push_back(ptArray[dim_ * idx + 1]);
151			coords.push_back(ptArray[dim_ * idx + 2]);
152
153			StuntDouble* sd = bodydoubles[idx];
154
155			Vector3d vel = sd->getVel();
156			vels.push_back(vel.x());
157			vels.push_back(vel.y());
158			vels.push_back(vel.z());
159
160			masses.push_back(sd->getMass());
161			}
162
163			MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0],
164			1, MPI::INT);
165
166			int globalHullSites = 0;
167			for (int iproc = 0; iproc < nproc; iproc++){
168			globalHullSites += hullSitesOnProc[iproc];
169			coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc];
170			}
171
172			displacements[0] = 0;
173			vectorDisplacements[0] = 0;
174
175			for (int iproc = 1; iproc < nproc; iproc++){
176			displacements[iproc] = displacements[iproc-1] + hullSitesOnProc[iproc-1];
177			vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1];
178			}
179
180			std::vector<double> globalCoords(dim_ * globalHullSites);
181			std::vector<double> globalVels(dim_ * globalHullSites);
182			std::vector<double> globalMasses(globalHullSites);
183
184			int count = coordsOnProc[myrank];
185
186			MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0],
187			&coordsOnProc[0], &vectorDisplacements[0],
188			MPI::DOUBLE);
189
190			MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0],
191			&coordsOnProc[0], &vectorDisplacements[0],
192			MPI::DOUBLE);
193
194			MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE,
195			&globalMasses[0], &hullSitesOnProc[0],
196			&displacements[0], MPI::DOUBLE);
197
198			// Free previous hull
199			qh_freeqhull(!qh_ALL);
200			qh_memfreeshort(&curlong, &totlong);
201			if (curlong \|\| totlong)
202			std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
203			<< totlong << curlong << std::endl;
204
205			if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc,
206			const_cast<char *>(options_.c_str()), NULL, stderr)){
207
208			sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute global convex hull");
209			painCave.isFatal = 1;
210			simError();
211
212			} //qh_new_qhull
213
214
215			#endif
216
217			//Set facet->center as the Voronoi center
218			qh_setvoronoi_all();
219
220
221			int convexNumVert = qh_setsize(qh_facetvertices (qh facet_list, NULL, false));
222			//Insert all the sample points, because, even with alpha=0, the alpha shape/alpha complex will
223			//contain them.
224
225			// tri::Allocator<CMeshO>::AddVertices(pm.cm,convexNumVert);
226
227			/*ivp length is 'qh num_vertices' because each vertex is accessed through its ID whose range is
228			0<=qh_pointid(vertex->point)<qh num_vertices*/
229			// vector<tri::Allocator<CMeshO>::VertexPointer> ivp(qh num_vertices);
230			/*i=0;
231			FORALLvertices{
232			if ((*vertex).point){
233			// pm.cm.vert[i].P()[0] = (*vertex).point[0];
234			// pm.cm.vert[i].P()[1] = (*vertex).point[1];
235			//pm.cm.vert[i].P()[2] = (*vertex).point[2];
236			// ivp[qh_pointid(vertex->point)] = &pm.cm.vert[i];
237			i++;
238			}
239			}
240			*/
241			//Set of alpha complex triangles for alphashape filtering
242			setT* set= qh_settemp(4* qh num_facets);
243
244			qh visit_id++;
245			int numFacets=0;
246			std::vector<std::vector <int> > facetlist;
247	chuckv	1404	interiorPoint = qh interior_point;
248	chuckv	1402	FORALLfacet_(qh facet_list) {
249			numFacets++;
250			if (!facet->upperdelaunay) {
251			//For all facets (that are tetrahedrons)calculate the radius of the empty circumsphere considering
252			//the distance between the circumcenter and a vertex of the facet
253			vertexT* vertex = (vertexT *)(facet->vertices->e[0].p);
254			double* center = facet->center;
255			double radius = qh_pointdist(vertex->point,center,dim_);
256
257			if (radius>alpha_) // if the facet is not good consider the ridges
258			{
259			//if calculating the alphashape, unmark the facet ('good' is used as 'marked').
260			facet->good=false;
261
262			//Compute each ridge (triangle) once and test the cironference radius with alpha
263			facet->visitid= qh visit_id;
264			qh_makeridges(facet);
265			ridgeT ridge, *ridgep;
266			int goodTriangles=0;
267			FOREACHridge_(facet->ridges) {
268			neighbor= otherfacet_(ridge, facet);
269			if (( neighbor->visitid != qh visit_id)){
270			//Calculate the radius of the circumference
271			pointT* p0 = ((vertexT*) (ridge->vertices->e[0].p))->point;
272			pointT* p1 = ((vertexT*) (ridge->vertices->e[1].p))->point;
273			pointT* p2 = ((vertexT*) (ridge->vertices->e[2].p))->point;
274
275			radius = calculate_circumradius(p0,p1,p2, dim_);
276
277			if(radius <=alpha_){
278			goodTriangles++;
279			//save the triangle (ridge) for subsequent filtering
280			qh_setappend(&set, ridge);
281			}
282			}
283			}
284
285			//If calculating the alphashape, mark the facet('good' is used as 'marked').
286			//This facet will have some triangles hidden by the facet's neighbor.
287			if(goodTriangles==4)
288			facet->good=true;
289
290			}
291			else //the facet is good. Put all the triangles of the tetrahedron in the mesh
292			{
293			//Compute each ridge (triangle) once
294			facet->visitid= qh visit_id;
295			//If calculating the alphashape, mark the facet('good' is used as 'marked').
296			//This facet will have some triangles hidden by the facet's neighbor.
297			facet->good=true;
298			qh_makeridges(facet);
299			ridgeT ridge, *ridgep;
300			FOREACHridge_(facet->ridges) {
301			neighbor= otherfacet_(ridge, facet);
302			if ((neighbor->visitid != qh visit_id)){
303			qh_setappend(&set, ridge);
304			}
305			}
306			}
307			}
308			}
309			//assert(numFacets== qh num_facets);
310
311			//Filter the triangles (only the ones on the boundary of the alpha complex) and build the mesh
312	chuckv	1404
313
314	chuckv	1402
315			ridgeT ridge, *ridgep;
316			FOREACHridge_(set) {
317			if ((!ridge->top->good \|\| !ridge->bottom->good \|\| ridge->top->upperdelaunay \|\| ridge->bottom->upperdelaunay)){
318			// tri::Allocator<CMeshO>::FaceIterator fi=tri::Allocator<CMeshO>::AddFaces(pm.cm,1);
319			ridgesCount++;
320			int vertex_n, vertex_i;
321			Triangle face;
322
323			// Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]);
324			//face.setNormal(V3dNormal);
325
326
327	chuckv	1404	//coordT *center = qh_getcenter(ridge->vertices);
328			//std::cout << "Centers are " << center[0] << " " <<center[1] << " " << center[2] << std::endl;
329	chuckv	1402	//Vector3d V3dCentroid(center[0], center[1], center[2]);
330			//face.setCentroid(V3dCentroid);
331	chuckv	1404
332	chuckv	1402
333	chuckv	1404	Vector3d faceVel = V3Zero;
334	chuckv	1402	Vector3d p[3];
335	chuckv	1404	RealType faceMass = 0.0;
336	chuckv	1402
337			int ver = 0;
338			std::vector<int> virtexlist;
339			FOREACHvertex_i_(ridge->vertices){
340			int id = qh_pointid(vertex->point);
341			p[ver][0] = vertex->point[0];
342			p[ver][1] = vertex->point[1];
343			p[ver][2] = vertex->point[2];
344			Vector3d vel;
345			RealType mass;
346			ver++;
347			virtexlist.push_back(id);
348	chuckv	1404	// std::cout << "Ridge: " << ridgesCount << " Vertex " << id << std::endl;
349
350			vel = bodydoubles[id]->getVel();
351			mass = bodydoubles[id]->getMass();
352			face.addVertexSD(bodydoubles[id]);
353
354
355			faceVel = faceVel + vel;
356			faceMass = faceMass + mass;
357			} //FOREACH Vertex
358	chuckv	1402	facetlist.push_back(virtexlist);
359	chuckv	1404	face.addVertices(p[0],p[1],p[2]);
360			face.setFacetMass(faceMass);
361			face.setFacetVelocity(faceVel/3.0);
362
363			RealType area = face.getArea();
364			area_ += area;
365			Vector3d normal = face.getUnitNormal();
366			RealType offset = ((0.0-p[0][0])normal[0] + (0.0-p[0][1])normal[1] + (0.0-p[0][2])*normal[2]);
367			RealType dist = normal[0] * interiorPoint[0] + normal[1]interiorPoint[1] + normal[2]interiorPoint[2];
368			std::cout << "Dist and normal and area are: " << normal << std::endl;
369			volume_ += dist *area/qh hull_dim;
370
371			Triangles_.push_back(face);
372	chuckv	1402	}
373			}
374
375	chuckv	1404	std::cout << "Volume is: " << volume_ << std::endl;
376
377	chuckv	1402	//assert(pm.cm.fn == ridgesCount);
378	chuckv	1404	/*
379	chuckv	1402	std::cout <<"OFF"<<std::endl;
380			std::cout << bodydoubles.size() << " " << facetlist.size() << " " << 3*facetlist.size() << std::endl;
381			for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){
382			Vector3d pos = (*SD)->getPos();
383			std::cout << pos.x() << " " << pos.y() << " " << pos.z() << std::endl;
384			}
385
386
387			std::vector<std::vector<int> >::iterator thisfacet;
388			std::vector<int>::iterator thisvertex;
389
390			for (thisfacet = facetlist.begin(); thisfacet != facetlist.end(); thisfacet++){
391			std::cout << (*thisfacet).size();
392			for (thisvertex = (thisfacet).begin(); thisvertex != (thisfacet).end(); thisvertex++){
393			std::cout << " " << *thisvertex;
394			}
395			std::cout << std::endl;
396			}
397	chuckv	1404	*/
398	chuckv	1402
399
400
401			/*
402			FORALLfacets {
403			Triangle face;
404
405			Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]);
406			face.setNormal(V3dNormal);
407
408			RealType faceArea = qh_facetarea(facet);
409			face.setArea(faceArea);
410
411			vertices = qh_facet3vertex(facet);
412
413			coordT *center = qh_getcenter(vertices);
414			Vector3d V3dCentroid(center[0], center[1], center[2]);
415			face.setCentroid(V3dCentroid);
416
417			Vector3d faceVel = V3Zero;
418			Vector3d p[3];
419			RealType faceMass = 0.0;
420
421			int ver = 0;
422
423			FOREACHvertex_(vertices){
424			int id = qh_pointid(vertex->point);
425			p[ver][0] = vertex->point[0];
426			p[ver][1] = vertex->point[1];
427			p[ver][2] = vertex->point[2];
428
429			Vector3d vel;
430			RealType mass;
431
432			#ifdef IS_MPI
433			vel = Vector3d(globalVels[dim_ * id],
434			globalVels[dim_ * id + 1],
435			globalVels[dim_ * id + 2]);
436			mass = globalMasses[id];
437
438			// localID will be between 0 and hullSitesOnProc[myrank] if we
439			// own this guy.
440
441			int localID = id - displacements[myrank];
442
443			if (localID >= 0 && localID < hullSitesOnProc[myrank])
444			face.addVertexSD(bodydoubles[indexMap[localID]]);
445
446			#else
447			vel = bodydoubles[id]->getVel();
448			mass = bodydoubles[id]->getMass();
449			face.addVertexSD(bodydoubles[id]);
450			#endif
451
452			faceVel = faceVel + vel;
453			faceMass = faceMass + mass;
454			ver++;
455			} //Foreachvertex
456
457			face.addVertices(p[0], p[1], p[2]);
458			face.setFacetMass(faceMass);
459			face.setFacetVelocity(faceVel/3.0);
460			Triangles_.push_back(face);
461			qh_settempfree(&vertices);
462
463			} //FORALLfacets
464			*/
465	chuckv	1404	// qh_getarea(qh facet_list);
466			//volume_ = qh totvol;
467			// area_ = qh totarea;
468	chuckv	1402
469			qh_freeqhull(!qh_ALL);
470			qh_memfreeshort(&curlong, &totlong);
471			if (curlong \|\| totlong)
472			std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
473			<< totlong << curlong << std::endl;
474			}
475
476			void AlphaHull::printHull(const std::string& geomFileName) {
477
478			#ifdef IS_MPI
479			if (worldRank == 0) {
480			#endif
481			FILE *newGeomFile;
482
483			//create new .md file based on old .md file
484			newGeomFile = fopen(geomFileName.c_str(), "w");
485			qh_findgood_all(qh facet_list);
486			for (int i = 0; i < qh_PRINTEND; i++)
487			qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
488
489			fclose(newGeomFile);
490			#ifdef IS_MPI
491			}
492			#endif
493			}
494
495			double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim){
496			coordT a = qh_pointdist(p0,p1,dim);
497			coordT b = qh_pointdist(p1,p2,dim);
498			coordT c = qh_pointdist(p2,p0,dim);
499
500			coordT sum =(a + b + c)*0.5;
501			coordT area = sum(a+b-sum)(a+c-sum)*(b+c-sum);
502			return (double) (abc)/(4*sqrt(area));
503			}
504
505			#endif //QHULL