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 1097 by chuckv, Thu Dec 14 19:32:32 2006 UTC vs.
Revision 1188 by chuckv, Thu Nov 22 16:39:45 2007 UTC

# Line 40 | Line 40
40   *
41   *  ConvexHull.cpp
42   *
43 < *  Purpose: To calculate convexhull, hull volume and radius
44 < *  using the CGAL library.
43 > *  Purpose: To calculate convexhull, 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.1 2006-12-14 19:32:32 chuckv Exp $
47 > *  @version $Id: ConvexHull.cpp,v 1.5 2007-11-22 16:39:45 chuckv Exp $
48   *
49   */
50  
52 – #include "math/ConvexHull.hpp"
51   #include <iostream>
52   #include <fstream>
53 + #include <list>
54 + #include <algorithm>
55 + #include <iterator>
56 + #include "math/ConvexHull.hpp"
57 + #include "utils/simError.h"
58 + using namespace oopse;
59  
60 + /* Old options Qt Qu Qg QG0 FA */
61 + ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt  Qci Tcv Pp")
62 + //ConvexHull::ConvexHull() : dim_(3), options_("qhull d Qbb Qt i")
63 + {}
64  
65 < using namespace oopse;
65 > bool ConvexHull::genHull(std::vector<Vector3d> pos)
66 > {
67  
68  
69 +        int numpoints = pos.size();
70 +        coordT* pt_array;
71 +        coordT* surfpt_array;
72 +        std::list<int> surface_atoms;
73 +  FILE *outdummy = NULL;
74 +  FILE *errdummy = NULL;
75  
76 +        pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_));
77  
62 – bool ConvexHull::genHull(std::vector<Vector3d> pos)
63 – {
78  
79 <  std::vector<Point_3> points;
80 <  points.reserve(pos.size());
81 <  // Copy the positon vector into a points vector for cgal.
82 <  for (int i = 0; i < pos.size(); ++i)
83 <    {
70 <      Point_3 pt(pos[i][0],pos[i][1],pos[i][2]);
71 <      points.push_back(pt);
72 <    }
79 >        for (int i = 0; i < numpoints; i++) {
80 >                pt_array[dim_ * i] = pos[i][0];
81 >                pt_array[dim_ * i + 1] = pos[i][1];
82 >                pt_array[dim_ * i + 2] = pos[i][2];
83 >        }
84  
74 –  // compute convex hull
75 –  CGAL::convex_hull_3(points.begin(), points.end(), ch_object);
76 –  // Make sure hull is a polyhedron
77 –  if ( CGAL::assign(ch_polyhedron, ch_object) )
78 –    {
79 –      return true;
80 –    }
81 –  else
82 –    {
83 –      return false;
84 –    }
85 – }
85  
87 – RealType ConvexHull::getVolume()
88 – {
89 –        /*
90 –  std::list<Point> L;
91 –  L.push_front(Point(0,0,0));
92 –  L.push_front(Point(1,0,0));
93 –  L.push_front(Point(0,1,0));
86  
95 –  Triangulation T(L.begin(), L.end());
87  
88 <  int n = T.number_of_vertices();
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();
93  
94 <  // insertion from a vector :
95 <  std::vector<Point> V(3);
96 <  V[0] = Point(0,0,1);
102 <  V[1] = Point(1,1,1);
103 <  V[2] = Point(2,2,2);
94 >                qh_produce_output();
95 >        */
96 >        boolT ismalloc = False;
97  
98 <  n = n + T.insert(V.begin(), V.end());
98 >        if (!qh_new_qhull(dim_, numpoints, pt_array, ismalloc,
99 >                          const_cast<char *>(options_.c_str()), NULL, stderr)) {
100  
101 <  assert( n == 6 );       // 6 points have been inserted
102 <  assert( T.is_valid() ); // checking validity of T
101 >                vertexT *vertex, **vertexp;
102 >                facetT *facet;
103 >                setT *vertices;
104 >                unsigned int nf = qh num_facets;
105  
106 <  double sum_v = 0;
107 <  for(Triangulation::Cell_iterator c_it = T.cells_begin(); c_it != T.cells_end(); ++c_it)
108 <    {
109 <      sum_v += T.tetrahedron(c_it).volume();
110 <    }
111 <  std::cout << "sum_v " << sum_v << std::endl;
106 >                //Matrix idx(nf, dim);
107 > /*
108 >                int j, i = 0, id = 0;
109 >                        
110 >                        int id2 = 0;
111 >                        coordT *point,**pointp;
112 >                        realT dist;
113 >                        FORALLfacets {
114 >                            j = 0;
115 >                        
116 >                            if (!facet->simplicial){
117 >                            // should never happen with Qt
118 >                            sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
119 >                                    painCave.isFatal = 0;
120 >                                    simError();
121 >                            }
122 >                        
123 >                            vertices = qh_facet3vertex(facet);
124 >                            FOREACHvertex_(vertices){
125 >                                                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;
139 >                                        
140 >                                }
141 >                                
142 >                        }
143 >                
144   */
145 <        return 0.0;
145 >                
146 >                
147 >        }
148 >
149 >
150 >
151 >
152 >        qh_getarea(qh facet_list);
153 >        volume_ = qh totvol;
154 >        area_ = qh totarea;
155 > //      FILE *newGeomFile;
156 >
157 >
158 >        /*
159 >                FORALLfacets {
160 >                    for (int k=0; k < qh hull_dim; k++)
161 >                      printf ("%6.2g ", facet->normal[k]);
162 >                    printf ("\n");
163 >                  }
164 >        */
165 >
166 >        int curlong,totlong;
167 > //      geomviewHull("junk.off");
168 >        qh_freeqhull(!qh_ALL);
169 >        qh_memfreeshort(&curlong, &totlong);
170 >        if (curlong || totlong)
171 >                std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
172 >                << totlong << curlong << std::endl;
173 >        free(pt_array);
174 > /*
175 >        int j = 0;
176 >        surface_atoms.sort();
177 >        surface_atoms.unique();
178 >        surfpt_array = (coordT*) malloc(sizeof(coordT) * (surface_atoms.size() * dim_));
179 >        for(std::list<int>::iterator list_iter = surface_atoms.begin();
180 >            list_iter != surface_atoms.end(); list_iter++)
181 >        {
182 >                int i = *list_iter;
183 >                        //surfpt_array[dim_ * j] = pos[i][0];
184 >                        //surfpt_array[dim_ * j + 1] = pos[i][1];
185 >                        //surfpt_array[dim_ * j + 2] = pos[i][2];
186 >                        std::cout << "Ag  " << pos[i][0] << "  " << pos[i][1] << "  "<< pos[i][2] << std::endl;
187 >                        j++;
188 >        }
189 > */      
190 >        
191 > /*      
192 >        std::string deloptions_ = "qhull d Qt";
193 >        facetT *facet, *neighbor;
194 >        ridgeT *ridge, **ridgep;
195 >
196 >        if (!qh_new_qhull(dim_, surface_atoms.size(), surfpt_array, ismalloc,
197 >        const_cast<char *>(deloptions_.c_str()), NULL, stderr)) {
198 >
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 >                }
214 >
215 >
216 >
217 >
218   }
219  
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 + }
229 +
230 +
231 +
232 + RealType ConvexHull::getVolume()
233 + {
234 +        return volume_;
235 + }
236 +
237   void ConvexHull::geomviewHull(const std::string& geomFileName)
238   {
239  
240 <  std::ofstream newGeomFile;
240 >        FILE *newGeomFile;
241  
242 <  //create new .md file based on old .md file
243 <  newGeomFile.open(geomFileName.c_str());
242 >        //create new .md file based on old .md file
243 >        newGeomFile = fopen(geomFileName.c_str(), "w");
244 >        qh_findgood_all(qh facet_list);
245 >        for (int i = 0; i < qh_PRINTEND; i++)
246 >                qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
247  
248 <  // Write polyhedron in Object File Format (OFF).
249 <  CGAL::set_ascii_mode( std::cout);
130 <  newGeomFile << "OFF" << std::endl << ch_polyhedron.size_of_vertices() << ' '
131 <  << ch_polyhedron.size_of_facets() << " 0" << std::endl;
132 <  std::copy( ch_polyhedron.points_begin(), ch_polyhedron.points_end(),
133 <             std::ostream_iterator<Point_3>( newGeomFile, "\n"));
134 <  for (  Facet_iterator i = ch_polyhedron.facets_begin(); i != ch_polyhedron.facets_end(); ++i)
135 <    {
136 <      Halfedge_facet_circulator j = i->facet_begin();
137 <      // Facets in polyhedral surfaces are at least triangles.
138 <      CGAL_assertion( CGAL::circulator_size(j) >= 3);
139 <      newGeomFile << CGAL::circulator_size(j) << ' ';
140 <      do
141 <        {
142 <          newGeomFile << ' ' << std::distance(ch_polyhedron.vertices_begin(), j->vertex());
143 <        }
144 <      while ( ++j != i->facet_begin());
145 <      newGeomFile << std::endl;
146 <    }
248 >        fclose(newGeomFile);
249 > }
250  
148 –  newGeomFile.close();
251  
252  
253 < }
253 >

Diff Legend

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