ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/math/ConvexHull.cpp
(Generate patch)

Comparing trunk/src/math/ConvexHull.cpp (file contents):
Revision 1137 by chuckv, Tue May 29 22:50:14 2007 UTC vs.
Revision 1242 by chuckv, Wed May 14 14:31:48 2008 UTC

# Line 44 | Line 44
44   *
45   *  Created by Charles F. Vardeman II on 11 Dec 2006.
46   *  @author  Charles F. Vardeman II
47 < *  @version $Id: ConvexHull.cpp,v 1.2 2007-05-29 22:50:14 chuckv Exp $
47 > *  @version $Id: ConvexHull.cpp,v 1.6 2008-05-14 14:31:48 chuckv Exp $
48   *
49   */
50  
51 < #include "math/ConvexHull.hpp"
51 > /* Standard includes independent of library */
52   #include <iostream>
53   #include <fstream>
54 + #include <list>
55 + #include <algorithm>
56 + #include <iterator>
57 + #include "math/ConvexHull.hpp"
58 + #include "utils/simError.h"
59 + using namespace oopse;
60  
61 < char options[] = "qhull Qt FA";
62 < int dim_ = 3;
61 > /* CGAL version of convex hull first then QHULL */
62 > #ifdef HAVE_CGAL
63  
64 < using namespace oopse;
64 > #include <CGAL/basic.h>
65 > #include <CGAL/Simple_cartesian.h>
66 > #include <CGAL/Convex_hull_traits_3.h>
67 > #include <CGAL/convex_hull_3.h>
68 > #include <CGAL/double.h>
69  
70 + #include <CGAL/Quotient.h>
71 + #include <CGAL/MP_Float.h>
72 + #include <CGAL/Lazy_exact_nt.h>
73 +
74 + typedef double RT;
75 + typedef CGAL::Simple_cartesian<RT>                K;
76 + typedef CGAL::Convex_hull_traits_3<K>             Traits;
77 + typedef Traits::Polyhedron_3                      Polyhedron_3;
78 + typedef K::Point_3                                Point;
79 +
80 +
81   ConvexHull::ConvexHull(){}
82  
83 + bool ConvexHull::genHull(std::vector<Vector3d> pos)
84 + {
85 +  
86 +  std::vector<Point> points;    
87 +  
88 +  
89 +  // Copy the positon vector into a points vector for cgal.
90 +  for (int i = 0; i < pos.size(); ++i)
91 +    {
92 +      Point pt(pos[i][0],pos[i][1],pos[i][2]);
93 +      points.push_back(pt);
94 +    }
95 +  
96 +  // define object to hold convex hull
97 +  Polyhedron_3 ch_object_;
98 +  // compute convex hull
99 +  CGAL::convex_hull_3(points.begin(), points.end(), ch_object_);
100 +  
101 +  for (Polyhedron_3::Vertex_iterator v = ch_object_.vertices_begin(); ch_object_.vertices_end(); ++v){
102 +    std::cout<< v.point()<<std::endl;
103 +  }
104 +  
105 +  
106 + }
107  
108 +
109 +
110 + #else
111 + #ifdef HAVE_QHULL
112 + /* Old options Qt Qu Qg QG0 FA */
113 + ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt  Qci Tcv Pp")
114 +                           //ConvexHull::ConvexHull() : dim_(3), options_("qhull d Qbb Qt i")
115 + {}
116 +
117   bool ConvexHull::genHull(std::vector<Vector3d> pos)
118   {
119 <        FILE *outfile = stdout;
120 <        FILE *errfile = stderr;
121 <        facetT *facet;
122 <        int exitcode;
123 <        boolT ismalloc = False;
124 <        int curlong,totlong;
119 >  
120 >  
121 >  int numpoints = pos.size();
122 >  coordT* pt_array;
123 >  coordT* surfpt_array;
124 >  std::list<int> surface_atoms;
125 >  FILE *outdummy = NULL;
126 >  FILE *errdummy = NULL;
127 >  
128 >  pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_));
129 >  
130 >  
131 >  for (int i = 0; i < numpoints; i++) {
132 >    pt_array[dim_ * i] = pos[i][0];
133 >    pt_array[dim_ * i + 1] = pos[i][1];
134 >    pt_array[dim_ * i + 2] = pos[i][2];
135 >  }
136 >  
137 >  
138 >  
139 >  
140 >  /*
141 >    qh_initflags(const_cast<char *>(options_.c_str()));
142 >    qh_init_B(pospoints, numpoints, dim_, ismalloc);
143 >    qh_qhull();
144 >    qh_check_output();
145 >
146 >    qh_produce_output();
147 >  */
148 >  boolT ismalloc = False;
149 >  
150 >  if (!qh_new_qhull(dim_, numpoints, pt_array, ismalloc,
151 >                    const_cast<char *>(options_.c_str()), NULL, stderr)) {
152 >    
153 >    vertexT *vertex, **vertexp;
154 >    facetT *facet;
155 >    setT *vertices;
156 >    unsigned int nf = qh num_facets;
157 >    
158 >    //Matrix idx(nf, dim);
159 >    /*
160 >      int j, i = 0, id = 0;
161 >      
162 >      int id2 = 0;
163 >      coordT *point,**pointp;
164 >      realT dist;
165 >      FORALLfacets {
166 >      j = 0;
167 >      
168 >      if (!facet->simplicial){
169 >      // should never happen with Qt
170 >      sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
171 >      painCave.isFatal = 0;
172 >      simError();
173 >      }
174 >                        
175 >      vertices = qh_facet3vertex(facet);
176 >      FOREACHvertex_(vertices){
177 >      id = qh_pointid(vertex->point);
178 >      surface_atoms.push_back(id);
179 >      //std::cout << "Ag  " << pos[id][0] << "    " << pos[id][1] << "    " << pos[id][2]<< std::endl;
180 >      }
181 >      qh_settempfree(&vertices);
182 >      
183 >      FOREACHpoint_(facet->coplanarset){
184 >      vertex= qh_nearvertex (facet, point, &dist);
185 >      //id= qh_pointid (vertex->point);
186 >      id2= qh_pointid (point);
187 >      surface_atoms.push_back(id2);
188 >      //std::cout << "Ag  " << pos[id2][0] << "    " << pos[id2][1] << "    " << pos[id2][2]<< std::endl;
189 >      //printf ("%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist);
190 >      //std::cout << "Neighbors are: %d $d %d\n" << id << id2 << facet->id;
191 >                                        
192 >      }
193 >      
194 >      }
195 >                
196 > */
197 >                
198 >                
199 >  }
200 >
201 >
202 >
203 >
204 >  qh_getarea(qh facet_list);
205 >  volume_ = qh totvol;
206 >  area_ = qh totarea;
207 >  //    FILE *newGeomFile;
208 >  
209 >  
210 >  /*
211 >    FORALLfacets {
212 >    for (int k=0; k < qh hull_dim; k++)
213 >    printf ("%6.2g ", facet->normal[k]);
214 >    printf ("\n");
215 >    }
216 >  */
217 >  
218 >  int curlong,totlong;
219 >  //    geomviewHull("junk.off");
220 >  qh_freeqhull(!qh_ALL);
221 >  qh_memfreeshort(&curlong, &totlong);
222 >  if (curlong || totlong)
223 >    std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
224 >              << totlong << curlong << std::endl;
225 >  free(pt_array);
226 >  /*
227 >    int j = 0;
228 >    surface_atoms.sort();
229 >    surface_atoms.unique();
230 >    surfpt_array = (coordT*) malloc(sizeof(coordT) * (surface_atoms.size() * dim_));
231 >    for(std::list<int>::iterator list_iter = surface_atoms.begin();
232 >    list_iter != surface_atoms.end(); list_iter++)
233 >    {
234 >    int i = *list_iter;
235 >    //surfpt_array[dim_ * j] = pos[i][0];
236 >    //surfpt_array[dim_ * j + 1] = pos[i][1];
237 >    //surfpt_array[dim_ * j + 2] = pos[i][2];
238 >    std::cout << "Ag  " << pos[i][0] << "  " << pos[i][1] << "  "<< pos[i][2] << std::endl;
239 >    j++;
240 >    }
241 >  */    
242 >  
243 >  /*    
244 >        std::string deloptions_ = "qhull d Qt";
245 >        facetT *facet, *neighbor;
246 >        ridgeT *ridge, **ridgep;
247          
248 <        int numpoints = pos.size();
248 >        if (!qh_new_qhull(dim_, surface_atoms.size(), surfpt_array, ismalloc,
249 >        const_cast<char *>(deloptions_.c_str()), NULL, stderr)) {
250          
251 <        coordT points[numpoints][dim_];
251 >        qh visit_id++;
252 >        FORALLfacets {
253 >        if (!facet->upperdelaunay) {
254 >        facet->visitid= qh visit_id;
255 >        qh_makeridges(facet);
256 >        FOREACHridge_(facet->ridges) {
257 >        neighbor= otherfacet_(ridge, facet);
258 >        if (neighbor->visitid != qh visit_id) {
259          
260 <        for (int i=0; i<numpoints; i++)
261 <        {
262 <     points[i][0] = pos[i][0];
79 <     points[i][1] = pos[i][1];
80 <     points[i][2] = pos[i][2];          
260 >        FOREACHvertex_(ridge->vertices)
261 >        int id2 = qh_pointid (vertex->point);
262 >        std::cout << "Ag  " << pos[id2][0] << "    " << pos[id2][1] << "    " << pos[id2][2]<< std::endl;
263          }
264 <  
264 >        }
265 >        }
266 >        
267  
268  
85 –        qh_initflags (options);
86 –    qh_init_B (points[0], numpoints, dim_, ismalloc);
87 –    qh_qhull();
88 –    qh_check_output();
89 –                                                        
90 –                                                        
91 –                                                        
92 –    qh_getarea(qh facet_list);
93 –    volume_ = qh totvol;                                                        
94 –        area_ = qh totarea;                                            
95 –                                                        
96 –                                                        
97 –                                                        
98 –    qh_freeqhull(!qh_ALL);
99 –        qh_memfreeshort (&curlong, &totlong);
100 –        if (curlong || totlong)
101 –                fprintf (errfile, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n",
102 –                                 totlong, curlong);
269          
270 +        }
271  
272 <        return true;
272 >        qh_freeqhull(!qh_ALL);
273 >        qh_memfreeshort(&curlong, &totlong);
274 >        if (curlong || totlong)
275 >        std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
276 >        << totlong << curlong << std::endl;
277 >        free(surfpt_array);
278 >  */            
279 >  return true;
280   }
281  
282  
283 +
284   RealType ConvexHull::getVolume()
285   {
286 <        return volume_;
286 >  return volume_;
287   }
288  
289   void ConvexHull::geomviewHull(const std::string& geomFileName)
290   {
291  
292 <  std::ofstream newGeomFile;
293 <
292 >  FILE *newGeomFile;
293 >  
294    //create new .md file based on old .md file
295 <  newGeomFile.open(geomFileName.c_str());
295 >  newGeomFile = fopen(geomFileName.c_str(), "w");
296 >  qh_findgood_all(qh facet_list);
297 >  for (int i = 0; i < qh_PRINTEND; i++)
298 >    qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
299 >  
300 >  fclose(newGeomFile);
301 > }
302 > #endif //QHULL
303 > #endif //CGAL
304  
305  
123 –  newGeomFile.close();
306  
125 –
126 – }

Diff Legend

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