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

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

Diff Legend

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