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]  Vardeman & Gezelter, in progress (2009).                        
 *
 *
 *  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: ConvexHull.cpp,v 1.21 2009-11-25 20:02:01 gezelter Exp $
 *
 */

/* 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/qhull.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 ") {
}

void AlphaHull::computeHull(std::vector<StuntDouble*> bodydoubles) { 
 
  int numpoints = bodydoubles.size();
  bool alphashape=true;
  
  Triangles_.clear();
  
  vertexT *vertex, **vertexp;
  facetT *facet, *neighbor;
  setT *vertices;
  int curlong, totlong;
  
  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;
  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);
      
      // 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;
      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);
      }
      facetlist.push_back(virtexlist);

    }
  }
                

//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:	1402
Committed:	Fri Jan 8 17:15:27 2010 UTC (15 years, 3 months ago) by chuckv
Original Path:	*trunk/src/math/AlphaHull.cpp*
File size:	15579 byte(s)
Log Message:	Added preliminary code for Alpha Hull calculation using qhull. Added preliminary support to SMIPD to support Alpha Hull. Alpha Hull does not yet add the correct things to triangle to be returned to SMPID. Preliminary changes for shadow hamiltonian integrator. Chages to md files so they will work in openMD.
#	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			* [4] Vardeman & Gezelter, in progress (2009).
39			*
40			*
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			* @version $Id: ConvexHull.cpp,v 1.21 2009-11-25 20:02:01 gezelter Exp $
48			*
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			#include <qhull/qhull.h>
72			#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			AlphaHull::AlphaHull(double alpha) : Hull(), dim_(3), alpha_(alpha), options_("qhull d QJ Tcv ") {
83			}
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			setT *vertices;
95			int curlong, totlong;
96
97			std::vector<double> ptArray(numpoints*dim_);
98
99			// Copy the positon vector into a points vector for qhull.
100			std::vector<StuntDouble*>::iterator SD;
101			int i = 0;
102			for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){
103			Vector3d pos = (*SD)->getPos();
104			ptArray[dim_ * i] = pos.x();
105			ptArray[dim_ * i + 1] = pos.y();
106			ptArray[dim_ * i + 2] = pos.z();
107			i++;
108			}
109
110			/* Clean up memory from previous convex hull calculations*/
111			boolT ismalloc = False;
112
113			int ridgesCount=0;
114			if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc,
115			const_cast<char *>(options_.c_str()), NULL, stderr)) {
116
117			sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute convex hull");
118			painCave.isFatal = 1;
119			simError();
120
121			} //qh_new_qhull
122
123
124			#ifdef IS_MPI
125			//If we are doing the mpi version, set up some vectors for data communication
126
127			int nproc = MPI::COMM_WORLD.Get_size();
128			int myrank = MPI::COMM_WORLD.Get_rank();
129			int localHullSites = 0;
130
131			std::vector<int> hullSitesOnProc(nproc, 0);
132			std::vector<int> coordsOnProc(nproc, 0);
133			std::vector<int> displacements(nproc, 0);
134			std::vector<int> vectorDisplacements(nproc, 0);
135
136			std::vector<double> coords;
137			std::vector<double> vels;
138			std::vector<int> indexMap;
139			std::vector<double> masses;
140
141			FORALLvertices{
142			localHullSites++;
143
144			int idx = qh_pointid(vertex->point);
145
146			indexMap.push_back(idx);
147
148			coords.push_back(ptArray[dim_ * idx]);
149			coords.push_back(ptArray[dim_ * idx + 1]);
150			coords.push_back(ptArray[dim_ * idx + 2]);
151
152			StuntDouble* sd = bodydoubles[idx];
153
154			Vector3d vel = sd->getVel();
155			vels.push_back(vel.x());
156			vels.push_back(vel.y());
157			vels.push_back(vel.z());
158
159			masses.push_back(sd->getMass());
160			}
161
162			MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0],
163			1, MPI::INT);
164
165			int globalHullSites = 0;
166			for (int iproc = 0; iproc < nproc; iproc++){
167			globalHullSites += hullSitesOnProc[iproc];
168			coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc];
169			}
170
171			displacements[0] = 0;
172			vectorDisplacements[0] = 0;
173
174			for (int iproc = 1; iproc < nproc; iproc++){
175			displacements[iproc] = displacements[iproc-1] + hullSitesOnProc[iproc-1];
176			vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1];
177			}
178
179			std::vector<double> globalCoords(dim_ * globalHullSites);
180			std::vector<double> globalVels(dim_ * globalHullSites);
181			std::vector<double> globalMasses(globalHullSites);
182
183			int count = coordsOnProc[myrank];
184
185			MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0],
186			&coordsOnProc[0], &vectorDisplacements[0],
187			MPI::DOUBLE);
188
189			MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0],
190			&coordsOnProc[0], &vectorDisplacements[0],
191			MPI::DOUBLE);
192
193			MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE,
194			&globalMasses[0], &hullSitesOnProc[0],
195			&displacements[0], MPI::DOUBLE);
196
197			// Free previous hull
198			qh_freeqhull(!qh_ALL);
199			qh_memfreeshort(&curlong, &totlong);
200			if (curlong \|\| totlong)
201			std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
202			<< totlong << curlong << std::endl;
203
204			if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc,
205			const_cast<char *>(options_.c_str()), NULL, stderr)){
206
207			sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute global convex hull");
208			painCave.isFatal = 1;
209			simError();
210
211			} //qh_new_qhull
212
213
214			#endif
215
216			//Set facet->center as the Voronoi center
217			qh_setvoronoi_all();
218
219
220			int convexNumVert = qh_setsize(qh_facetvertices (qh facet_list, NULL, false));
221			//Insert all the sample points, because, even with alpha=0, the alpha shape/alpha complex will
222			//contain them.
223
224			// tri::Allocator<CMeshO>::AddVertices(pm.cm,convexNumVert);
225
226			/*ivp length is 'qh num_vertices' because each vertex is accessed through its ID whose range is
227			0<=qh_pointid(vertex->point)<qh num_vertices*/
228			// vector<tri::Allocator<CMeshO>::VertexPointer> ivp(qh num_vertices);
229			/*i=0;
230			FORALLvertices{
231			if ((*vertex).point){
232			// pm.cm.vert[i].P()[0] = (*vertex).point[0];
233			// pm.cm.vert[i].P()[1] = (*vertex).point[1];
234			//pm.cm.vert[i].P()[2] = (*vertex).point[2];
235			// ivp[qh_pointid(vertex->point)] = &pm.cm.vert[i];
236			i++;
237			}
238			}
239			*/
240			//Set of alpha complex triangles for alphashape filtering
241			setT* set= qh_settemp(4* qh num_facets);
242
243			qh visit_id++;
244			int numFacets=0;
245			std::vector<std::vector <int> > facetlist;
246			FORALLfacet_(qh facet_list) {
247			numFacets++;
248			if (!facet->upperdelaunay) {
249			//For all facets (that are tetrahedrons)calculate the radius of the empty circumsphere considering
250			//the distance between the circumcenter and a vertex of the facet
251			vertexT* vertex = (vertexT *)(facet->vertices->e[0].p);
252			double* center = facet->center;
253			double radius = qh_pointdist(vertex->point,center,dim_);
254
255			if (radius>alpha_) // if the facet is not good consider the ridges
256			{
257			//if calculating the alphashape, unmark the facet ('good' is used as 'marked').
258			facet->good=false;
259
260			//Compute each ridge (triangle) once and test the cironference radius with alpha
261			facet->visitid= qh visit_id;
262			qh_makeridges(facet);
263			ridgeT ridge, *ridgep;
264			int goodTriangles=0;
265			FOREACHridge_(facet->ridges) {
266			neighbor= otherfacet_(ridge, facet);
267			if (( neighbor->visitid != qh visit_id)){
268			//Calculate the radius of the circumference
269			pointT* p0 = ((vertexT*) (ridge->vertices->e[0].p))->point;
270			pointT* p1 = ((vertexT*) (ridge->vertices->e[1].p))->point;
271			pointT* p2 = ((vertexT*) (ridge->vertices->e[2].p))->point;
272
273			radius = calculate_circumradius(p0,p1,p2, dim_);
274
275			if(radius <=alpha_){
276			goodTriangles++;
277			//save the triangle (ridge) for subsequent filtering
278			qh_setappend(&set, ridge);
279			}
280			}
281			}
282
283			//If calculating the alphashape, mark the facet('good' is used as 'marked').
284			//This facet will have some triangles hidden by the facet's neighbor.
285			if(goodTriangles==4)
286			facet->good=true;
287
288			}
289			else //the facet is good. Put all the triangles of the tetrahedron in the mesh
290			{
291			//Compute each ridge (triangle) once
292			facet->visitid= qh visit_id;
293			//If calculating the alphashape, mark the facet('good' is used as 'marked').
294			//This facet will have some triangles hidden by the facet's neighbor.
295			facet->good=true;
296			qh_makeridges(facet);
297			ridgeT ridge, *ridgep;
298			FOREACHridge_(facet->ridges) {
299			neighbor= otherfacet_(ridge, facet);
300			if ((neighbor->visitid != qh visit_id)){
301			qh_setappend(&set, ridge);
302			}
303			}
304			}
305			}
306			}
307			//assert(numFacets== qh num_facets);
308
309			//Filter the triangles (only the ones on the boundary of the alpha complex) and build the mesh
310
311			ridgeT ridge, *ridgep;
312			FOREACHridge_(set) {
313			if ((!ridge->top->good \|\| !ridge->bottom->good \|\| ridge->top->upperdelaunay \|\| ridge->bottom->upperdelaunay)){
314			// tri::Allocator<CMeshO>::FaceIterator fi=tri::Allocator<CMeshO>::AddFaces(pm.cm,1);
315			ridgesCount++;
316			int vertex_n, vertex_i;
317			Triangle face;
318
319			// Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]);
320			//face.setNormal(V3dNormal);
321
322			// RealType faceArea = qh_facetarea(facet);
323			//face.setArea(faceArea);
324
325			//vertices = qh_facet3vertex(facet);
326
327			//coordT *center = qh_getcenter(vertices);
328			//Vector3d V3dCentroid(center[0], center[1], center[2]);
329			//face.setCentroid(V3dCentroid);
330
331			//Vector3d faceVel = V3Zero;
332			Vector3d p[3];
333			//RealType faceMass = 0.0;
334
335			int ver = 0;
336			std::vector<int> virtexlist;
337			FOREACHvertex_i_(ridge->vertices){
338			int id = qh_pointid(vertex->point);
339			p[ver][0] = vertex->point[0];
340			p[ver][1] = vertex->point[1];
341			p[ver][2] = vertex->point[2];
342			Vector3d vel;
343			RealType mass;
344			ver++;
345			virtexlist.push_back(id);
346			}
347			facetlist.push_back(virtexlist);
348
349			}
350			}
351
352
353			//assert(pm.cm.fn == ridgesCount);
354
355			std::cout <<"OFF"<<std::endl;
356			std::cout << bodydoubles.size() << " " << facetlist.size() << " " << 3*facetlist.size() << std::endl;
357			for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){
358			Vector3d pos = (*SD)->getPos();
359			std::cout << pos.x() << " " << pos.y() << " " << pos.z() << std::endl;
360			}
361
362
363			std::vector<std::vector<int> >::iterator thisfacet;
364			std::vector<int>::iterator thisvertex;
365
366			for (thisfacet = facetlist.begin(); thisfacet != facetlist.end(); thisfacet++){
367			std::cout << (*thisfacet).size();
368			for (thisvertex = (thisfacet).begin(); thisvertex != (thisfacet).end(); thisvertex++){
369			std::cout << " " << *thisvertex;
370			}
371			std::cout << std::endl;
372			}
373
374
375
376
377			/*
378			FORALLfacets {
379			Triangle face;
380
381			Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]);
382			face.setNormal(V3dNormal);
383
384			RealType faceArea = qh_facetarea(facet);
385			face.setArea(faceArea);
386
387			vertices = qh_facet3vertex(facet);
388
389			coordT *center = qh_getcenter(vertices);
390			Vector3d V3dCentroid(center[0], center[1], center[2]);
391			face.setCentroid(V3dCentroid);
392
393			Vector3d faceVel = V3Zero;
394			Vector3d p[3];
395			RealType faceMass = 0.0;
396
397			int ver = 0;
398
399			FOREACHvertex_(vertices){
400			int id = qh_pointid(vertex->point);
401			p[ver][0] = vertex->point[0];
402			p[ver][1] = vertex->point[1];
403			p[ver][2] = vertex->point[2];
404
405			Vector3d vel;
406			RealType mass;
407
408			#ifdef IS_MPI
409			vel = Vector3d(globalVels[dim_ * id],
410			globalVels[dim_ * id + 1],
411			globalVels[dim_ * id + 2]);
412			mass = globalMasses[id];
413
414			// localID will be between 0 and hullSitesOnProc[myrank] if we
415			// own this guy.
416
417			int localID = id - displacements[myrank];
418
419			if (localID >= 0 && localID < hullSitesOnProc[myrank])
420			face.addVertexSD(bodydoubles[indexMap[localID]]);
421
422			#else
423			vel = bodydoubles[id]->getVel();
424			mass = bodydoubles[id]->getMass();
425			face.addVertexSD(bodydoubles[id]);
426			#endif
427
428			faceVel = faceVel + vel;
429			faceMass = faceMass + mass;
430			ver++;
431			} //Foreachvertex
432
433			face.addVertices(p[0], p[1], p[2]);
434			face.setFacetMass(faceMass);
435			face.setFacetVelocity(faceVel/3.0);
436			Triangles_.push_back(face);
437			qh_settempfree(&vertices);
438
439			} //FORALLfacets
440			*/
441			qh_getarea(qh facet_list);
442			volume_ = qh totvol;
443			area_ = qh totarea;
444
445			qh_freeqhull(!qh_ALL);
446			qh_memfreeshort(&curlong, &totlong);
447			if (curlong \|\| totlong)
448			std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
449			<< totlong << curlong << std::endl;
450			}
451
452			void AlphaHull::printHull(const std::string& geomFileName) {
453
454			#ifdef IS_MPI
455			if (worldRank == 0) {
456			#endif
457			FILE *newGeomFile;
458
459			//create new .md file based on old .md file
460			newGeomFile = fopen(geomFileName.c_str(), "w");
461			qh_findgood_all(qh facet_list);
462			for (int i = 0; i < qh_PRINTEND; i++)
463			qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
464
465			fclose(newGeomFile);
466			#ifdef IS_MPI
467			}
468			#endif
469			}
470
471			double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim){
472			coordT a = qh_pointdist(p0,p1,dim);
473			coordT b = qh_pointdist(p1,p2,dim);
474			coordT c = qh_pointdist(p2,p0,dim);
475
476			coordT sum =(a + b + c)*0.5;
477			coordT area = sum(a+b-sum)(a+c-sum)*(b+c-sum);
478			return (double) (abc)/(4*sqrt(area));
479			}
480
481			#endif //QHULL