From a7fdb2ca93d2a1b69ff8089c14ca66c98e91292e Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@gmail.com>
Date: Tue, 14 Jun 2011 20:41:17 -0400
Subject: [PATCH] Use C++ QHull for finding tesselation/convex hull

---
 mytools/loadZobov.cpp | 51 +++++++++++++++++++++++++++++++++++++++++++
 mytools/loadZobov.hpp |  7 ++++++
 vozutil.c             |  4 ++--
 3 files changed, 60 insertions(+), 2 deletions(-)

diff --git a/mytools/loadZobov.cpp b/mytools/loadZobov.cpp
index f4ed8ad..7399874 100644
--- a/mytools/loadZobov.cpp
+++ b/mytools/loadZobov.cpp
@@ -6,8 +6,10 @@
 #include <sstream>
 #include <algorithm>
 #include "loadZobov.hpp"
+#include <CosmoTool/fortran.hpp>
 
 using namespace std;
+using namespace CosmoTool;
 
 bool loadZobov(const char *descName, const char *adjName, const char *voidsName, 
 	       const char *volName, ZobovRep& z)
@@ -18,6 +20,7 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName,
   int32_t numParticles, numZones, numPinZone;
   int32_t totalParticles;
   int32_t numVoids;
+  int32_t minParticlesInZone, maxParticlesInZone;
 
   adjFile.read((char *)&numParticles, sizeof(numParticles));
   adjFile.read((char *)&numZones, sizeof(numZones));
@@ -28,6 +31,10 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName,
   cout << "Number of zones = " << numZones << endl;
 
   totalParticles = 0;
+
+  minParticlesInZone = -1;
+  maxParticlesInZone = -1;
+
   z.allZones.resize(numZones);
   for (int zone = 0; zone < numZones; zone++)
     {
@@ -41,10 +48,19 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName,
       z.allZones[zone].pId.resize(numPinZone);
       adjFile.read((char *)&z.allZones[zone].pId[0], sizeof(int)*numPinZone);
 
+      if (maxParticlesInZone < 0 || numPinZone > maxParticlesInZone)
+	maxParticlesInZone = numPinZone;
+
+      if (minParticlesInZone < 0 || numPinZone < minParticlesInZone)
+	minParticlesInZone = numPinZone;
+
       totalParticles += numPinZone;      
     }
   cout << "Zoned " << totalParticles << endl;
 
+  cout << "Minimum number of particles in zone = " << minParticlesInZone << endl;
+  cout << "Maximum number of particles in zone = " << maxParticlesInZone << endl;
+  
   if (totalParticles != numParticles)
     {
       cerr << "The numbers of particles are inconsistent ! (" << totalParticles << " vs " << numParticles << ")"<< endl;
@@ -155,5 +171,40 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName,
 
   
 
+  return true;
+}
+
+bool loadZobovParticles(const char *fname, std::vector<ZobovParticle>& particles)
+{
+  UnformattedRead f(fname);
+  int N;
+
+  f.beginCheckpoint();
+  N = f.readInt32();
+  f.endCheckpoint();
+
+  particles.resize(N);
+
+  f.beginCheckpoint();
+  for (int i = 0; i < N; i++)
+    {
+      particles[i].x = f.readReal32();
+    }
+  f.endCheckpoint();
+
+  f.beginCheckpoint();
+  for (int i = 0; i < N; i++)
+    {
+      particles[i].y = f.readReal32();
+    }
+  f.endCheckpoint();
+
+  f.beginCheckpoint();
+  for (int i = 0; i < N; i++)
+    {
+      particles[i].z = f.readReal32();
+    }
+  f.endCheckpoint();
+
   return true;
 }
diff --git a/mytools/loadZobov.hpp b/mytools/loadZobov.hpp
index 8f6a9d6..fca01c8 100644
--- a/mytools/loadZobov.hpp
+++ b/mytools/loadZobov.hpp
@@ -23,8 +23,15 @@ struct ZobovRep
   std::vector<float> particleVolume;
 };
 
+struct ZobovParticle
+{
+  float x, y, z;
+};
+
 bool loadZobov(const char *descName,
 	       const char *adjName, const char *voidName,
 	       const char *volName, ZobovRep& z);
 
+bool loadZobovParticles(const char *fname, std::vector<ZobovParticle>& particles);
+
 #endif
diff --git a/vozutil.c b/vozutil.c
index 1e2aee0..3ab8fb8 100644
--- a/vozutil.c
+++ b/vozutil.c
@@ -53,10 +53,10 @@ int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs)
     /* 'qh facet_list' contains the convex hull */
     
     /* From qh_printvneighbors */
-    qh_countfacets(qh facet_list, NULL, NULL, &numfacets, &numsimplicial, 
+    qh_countfacets(qh facet_list, NULL, 0, &numfacets, &numsimplicial, 
 		   &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
     qh_vertexneighbors();
-    vertices= qh_facetvertices (qh facet_list, NULL, NULL);
+    vertices= qh_facetvertices (qh facet_list, NULL, 0);
     vertex_points= qh_settemp (nvpall);
     coplanar_points= qh_settemp (nvpall);
     qh_setzero (vertex_points, 0, nvpall);