src/math/AlphaHull.cpp

/* Copyright (c) 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, 234107 (2008).          
 * [4] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 *
 *  AlphaHull.cpp
 *
 *  Purpose: To calculate an alpha-shape hull.
 */

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

/* Standard includes independent of library */

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

#include "math/qhull.hpp"

#ifdef HAVE_QHULL
using namespace std;
using namespace OpenMD;

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(vector<StuntDouble*> bodydoubles) { 
 
  int numpoints = bodydoubles.size();
  
  Triangles_.clear();
  
  vertexT *vertex;
  facetT *facet, *neighbor;
  pointT *interiorPoint;
  int curlong, totlong;
  
  
  vector<double> ptArray(numpoints*dim_);
  
  // Copy the positon vector into a points vector for qhull.
  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;
  
  /* compute the hull for our local points (or all the points for single
     processor versions) */
  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;
  int myrank;
  MPI_Comm_size( MPI_COMM_WORLD, &nproc);
  MPI_Comm_rank( MPI_COMM_WORLD, &myrank);

  int localHullSites = 0;

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

  vector<double> coords;
  vector<double> vels;
  vector<int> indexMap;
  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_Allgather(&localHullSites, 1, MPI_INT, &hullSitesOnProc[0],
                1, MPI_INT, MPI_COMM_WORLD);

  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]; 
  }

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

  int count = coordsOnProc[myrank];
  
  MPI_Allgatherv(&coords[0], count, MPI_DOUBLE, &globalCoords[0],
                 &coordsOnProc[0], &vectorDisplacements[0], 
                 MPI_DOUBLE, MPI_COMM_WORLD);
  
  MPI_Allgatherv(&vels[0], count, MPI_DOUBLE, &globalVels[0], 
                 &coordsOnProc[0], &vectorDisplacements[0],
                 MPI_DOUBLE, MPI_COMM_WORLD);
  
  MPI_Allgatherv(&masses[0], localHullSites, MPI_DOUBLE,
                 &globalMasses[0], &hullSitesOnProc[0], 
                 &displacements[0], MPI_DOUBLE, MPI_COMM_WORLD);
  
  // Free previous hull
  qh_freeqhull(!qh_ALL);
  qh_memfreeshort(&curlong, &totlong);
  if (curlong || totlong) {
    sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n"
            "\tdid not free %d bytes of long memory (%d pieces)", 
            totlong, curlong);
    painCave.isFatal = 1;
    simError();
  }
  
  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;
  vector<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

  int ridgesCount=0;
  
  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);
      //cout << "Centers are " << center[0] << "  " <<center[1] << "  " << center[2] << 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;
      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);
        // cout << "Ridge: " << ridgesCount << " Vertex " << id << 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 / RealType(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];
     cout << "Dist and normal and area are: " << normal << endl;
      volume_ += dist *area/qh hull_dim;
      
      Triangles_.push_back(face);
    }
  }

  cout << "Volume is: " << volume_ << 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) {
    sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n"
            "\tdid not free %d bytes of long memory (%d pieces)", 
            totlong, curlong);
    painCave.isFatal = 1;
    simError();
  }
}

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