# | Line 1 | Line 1 | |
---|---|---|
1 | < | /* Copyright (c) 2008 The University of Notre Dame. All Rights Reserved. |
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. Acknowledgement of the program authors must be made in any |
9 | < | * publication of scientific results based in part on use of the |
10 | < | * program. An acceptable form of acknowledgement is citation of |
11 | < | * the article in which the program was described (Matthew |
12 | < | * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
13 | < | * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
14 | < | * Parallel Simulation Engine for Molecular Dynamics," |
15 | < | * J. Comput. Chem. 26, pp. 252-271 (2005)) |
16 | < | * |
17 | < | * 2. Redistributions of source code must retain the above copyright |
8 | > | * 1. Redistributions of source code must retain the above copyright |
9 | * notice, this list of conditions and the following disclaimer. | |
10 | * | |
11 | < | * 3. Redistributions in binary form must reproduce the above copyright |
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. | |
# | Line 37 | Line 28 | |
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, 24107 (2008). |
38 | + | * [4] Vardeman & Gezelter, in progress (2009). |
39 | * | |
40 | + | * |
41 | * ConvexHull.cpp | |
42 | * | |
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.9 2008-10-07 17:12:48 chuckv Exp $ |
47 | > | * @version $Id: ConvexHull.cpp,v 1.21 2009-11-25 20:02:01 gezelter Exp $ |
48 | * | |
49 | */ | |
50 | ||
51 | /* Standard includes independent of library */ | |
52 | + | |
53 | #include <iostream> | |
54 | #include <fstream> | |
55 | #include <list> | |
# | Line 57 | Line 58 | |
58 | #include "math/ConvexHull.hpp" | |
59 | #include "utils/simError.h" | |
60 | ||
61 | + | #ifdef IS_MPI |
62 | + | #include <mpi.h> |
63 | + | #endif |
64 | ||
65 | < | using namespace oopse; |
65 | > | using namespace OpenMD; |
66 | ||
67 | < | /* CGAL version of convex hull first then QHULL */ |
68 | < | #ifdef HAVE_CGAL |
69 | < | //#include <CGAL/Homogeneous.h> |
70 | < | #include <CGAL/basic.h> |
71 | < | //#include <CGAL/Simple_cartesian.h> |
72 | < | #include <CGAL/Cartesian.h> |
73 | < | #include <CGAL/Origin.h> |
74 | < | #include <CGAL/Exact_predicates_exact_constructions_kernel.h> |
75 | < | #include <CGAL/Convex_hull_traits_3.h> |
76 | < | #include <CGAL/convex_hull_3.h> |
77 | < | #include <CGAL/Polyhedron_traits_with_normals_3.h> |
78 | < | #include <CGAL/Polyhedron_3.h> |
75 | < | #include <CGAL/double.h> |
76 | < | #include <CGAL/number_utils.h> |
67 | > | #ifdef HAVE_QHULL |
68 | > | extern "C" |
69 | > | { |
70 | > | #include <qhull/qhull.h> |
71 | > | #include <qhull/mem.h> |
72 | > | #include <qhull/qset.h> |
73 | > | #include <qhull/geom.h> |
74 | > | #include <qhull/merge.h> |
75 | > | #include <qhull/poly.h> |
76 | > | #include <qhull/io.h> |
77 | > | #include <qhull/stat.h> |
78 | > | } |
79 | ||
80 | + | ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp") { |
81 | + | } |
82 | ||
83 | < | //#include <CGAL/Quotient.h> |
84 | < | #include <CGAL/MP_Float.h> |
85 | < | //#include <CGAL/Lazy_exact_nt.h> |
83 | > | void ConvexHull::computeHull(std::vector<StuntDouble*> bodydoubles) { |
84 | > | |
85 | > | int numpoints = bodydoubles.size(); |
86 | ||
87 | + | Triangles_.clear(); |
88 | + | |
89 | + | vertexT *vertex, **vertexp; |
90 | + | facetT *facet; |
91 | + | setT *vertices; |
92 | + | int curlong, totlong; |
93 | + | pointT *intPoint; |
94 | + | |
95 | + | std::vector<double> ptArray(numpoints*dim_); |
96 | ||
97 | + | // Copy the positon vector into a points vector for qhull. |
98 | + | std::vector<StuntDouble*>::iterator SD; |
99 | + | int i = 0; |
100 | + | for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){ |
101 | + | Vector3d pos = (*SD)->getPos(); |
102 | + | ptArray[dim_ * i] = pos.x(); |
103 | + | ptArray[dim_ * i + 1] = pos.y(); |
104 | + | ptArray[dim_ * i + 2] = pos.z(); |
105 | + | i++; |
106 | + | } |
107 | + | |
108 | + | /* Clean up memory from previous convex hull calculations */ |
109 | + | boolT ismalloc = False; |
110 | + | |
111 | + | /* compute the hull for our local points (or all the points for single |
112 | + | processor versions) */ |
113 | + | if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc, |
114 | + | const_cast<char *>(options_.c_str()), NULL, stderr)) { |
115 | + | |
116 | + | sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull"); |
117 | + | painCave.isFatal = 1; |
118 | + | simError(); |
119 | + | |
120 | + | } //qh_new_qhull |
121 | ||
85 | – | typedef CGAL::MP_Float RT; |
86 | – | //typedef double RT; |
87 | – | //typedef CGAL::Homogeneous<RT> K; |
88 | – | typedef CGAL::Exact_predicates_exact_constructions_kernel K; |
89 | – | typedef K::Vector_3 Vector_3; |
90 | – | //typedef CGAL::Convex_hull_traits_3<K> Traits; |
91 | – | typedef CGAL::Polyhedron_traits_with_normals_3<K> Traits; |
92 | – | //typedef Traits::Polyhedron_3 Polyhedron_3; |
93 | – | typedef CGAL::Polyhedron_3<Traits> Polyhedron_3; |
94 | – | typedef K::Point_3 Point_3; |
122 | ||
123 | + | #ifdef IS_MPI |
124 | + | //If we are doing the mpi version, set up some vectors for data communication |
125 | + | |
126 | + | int nproc = MPI::COMM_WORLD.Get_size(); |
127 | + | int myrank = MPI::COMM_WORLD.Get_rank(); |
128 | + | int localHullSites = 0; |
129 | ||
130 | < | typedef Polyhedron_3::HalfedgeDS HalfedgeDS; |
131 | < | typedef Polyhedron_3::Facet_iterator Facet_iterator; |
132 | < | typedef Polyhedron_3::Halfedge_around_facet_circulator Halfedge_facet_circulator; |
133 | < | typedef Polyhedron_3::Halfedge_handle Halfedge_handle; |
101 | < | typedef Polyhedron_3::Facet_iterator Facet_iterator; |
102 | < | typedef Polyhedron_3::Plane_iterator Plane_iterator; |
103 | < | typedef Polyhedron_3::Vertex_iterator Vertex_iterator; |
104 | < | typedef Polyhedron_3::Vertex_handle Vertex_handle; |
105 | < | typedef Polyhedron_3::Point_iterator Point_iterator; |
106 | < | |
130 | > | std::vector<int> hullSitesOnProc(nproc, 0); |
131 | > | std::vector<int> coordsOnProc(nproc, 0); |
132 | > | std::vector<int> displacements(nproc, 0); |
133 | > | std::vector<int> vectorDisplacements(nproc, 0); |
134 | ||
135 | + | std::vector<double> coords; |
136 | + | std::vector<double> vels; |
137 | + | std::vector<int> indexMap; |
138 | + | std::vector<double> masses; |
139 | ||
140 | < | class Enriched_Point_3 : public K::Point_3{ |
141 | < | public: |
142 | < | Enriched_Point_3(double x,double y,double z) : K::Point_3(x,y,z), yupMyPoint(false), mySD(NULL) {} |
140 | > | FORALLvertices{ |
141 | > | localHullSites++; |
142 | > | |
143 | > | int idx = qh_pointid(vertex->point); |
144 | ||
145 | < | bool isMyPoint() const{ return yupMyPoint; } |
114 | < | void myPoint(){ yupMyPoint = true; } |
115 | < | void setSD(StuntDouble* SD){mySD = SD;} |
116 | < | StuntDouble* getStuntDouble(){return mySD;} |
117 | < | private: |
118 | < | bool yupMyPoint; |
119 | < | StuntDouble* mySD; |
145 | > | indexMap.push_back(idx); |
146 | ||
147 | < | }; |
147 | > | coords.push_back(ptArray[dim_ * idx]); |
148 | > | coords.push_back(ptArray[dim_ * idx + 1]); |
149 | > | coords.push_back(ptArray[dim_ * idx + 2]); |
150 | ||
151 | + | StuntDouble* sd = bodydoubles[idx]; |
152 | ||
153 | + | Vector3d vel = sd->getVel(); |
154 | + | vels.push_back(vel.x()); |
155 | + | vels.push_back(vel.y()); |
156 | + | vels.push_back(vel.z()); |
157 | ||
158 | + | masses.push_back(sd->getMass()); |
159 | + | } |
160 | ||
161 | + | MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], |
162 | + | 1, MPI::INT); |
163 | ||
164 | < | // compare Point_3's... used in setting up the STL map from points to indices |
165 | < | template <typename Pt3> |
166 | < | struct Point_3_comp { |
167 | < | bool operator() (const Pt3 & p, const Pt3 & q) const { |
131 | < | return CGAL::lexicographically_xyz_smaller(p,q); // this is defined inline & hence we had to create fn object & not ptrfun |
164 | > | int globalHullSites = 0; |
165 | > | for (int iproc = 0; iproc < nproc; iproc++){ |
166 | > | globalHullSites += hullSitesOnProc[iproc]; |
167 | > | coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc]; |
168 | } | |
133 | – | }; |
169 | ||
170 | < | // coordinate-based hashing inefficient but can we do better if pts are copied? |
171 | < | typedef std::map<Point_3, StuntDouble* ,Point_3_comp<Point_3> > ptMapType; |
170 | > | displacements[0] = 0; |
171 | > | vectorDisplacements[0] = 0; |
172 | > | |
173 | > | for (int iproc = 1; iproc < nproc; iproc++){ |
174 | > | displacements[iproc] = displacements[iproc-1] + hullSitesOnProc[iproc-1]; |
175 | > | vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1]; |
176 | > | } |
177 | ||
178 | < | #ifdef IS_MPI |
179 | < | struct { |
180 | < | double x,y,z; |
141 | < | } surfacePt; |
142 | < | #endif |
178 | > | std::vector<double> globalCoords(dim_ * globalHullSites); |
179 | > | std::vector<double> globalVels(dim_ * globalHullSites); |
180 | > | std::vector<double> globalMasses(globalHullSites); |
181 | ||
182 | < | ConvexHull::ConvexHull() : Hull(){ |
145 | < | //If we are doing the mpi version, set up some vectors for data communication |
146 | < | #ifdef IS_MPI |
147 | < | |
148 | < | |
149 | < | nproc_ = MPI::COMM_WORLD.Get_size(); |
150 | < | myrank_ = MPI::COMM_WORLD.Get_rank(); |
151 | < | NstoProc_ = new int[nproc_]; |
152 | < | displs_ = new int[nproc_]; |
153 | < | |
154 | < | // Create a surface point type in MPI to send |
155 | < | surfacePtType = MPI::DOUBLE.Create_contiguous(3); |
156 | < | surfacePtType.Commit(); |
157 | < | |
158 | < | |
159 | < | #endif |
160 | < | } |
161 | < | |
162 | < | void ConvexHull::computeHull(std::vector<StuntDouble*> bodydoubles) |
163 | < | { |
164 | < | |
165 | < | std::vector<Enriched_Point_3> points; |
166 | < | ptMapType myMap; |
167 | < | Point_iterator hc; |
182 | > | int count = coordsOnProc[myrank]; |
183 | ||
184 | < | // Copy the positon vector into a points vector for cgal. |
185 | < | std::vector<StuntDouble*>::iterator SD; |
184 | > | MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0], |
185 | > | &coordsOnProc[0], &vectorDisplacements[0], |
186 | > | MPI::DOUBLE); |
187 | ||
188 | < | for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD) |
189 | < | { |
190 | < | Vector3d pos = (*SD)->getPos(); |
175 | < | Enriched_Point_3* pt = new Enriched_Point_3(pos.x(),pos.y(),pos.z()); |
176 | < | pt->setSD(*SD); |
177 | < | points.push_back(*pt); |
178 | < | // myMap[pt]=(*SD); |
179 | < | } |
180 | < | |
181 | < | // define object to hold convex hull |
182 | < | CGAL::Object ch_object_; |
183 | < | Polyhedron_3 polyhedron; |
184 | < | |
185 | < | // compute convex hull |
186 | < | |
187 | < | std::vector<Enriched_Point_3>::iterator testpt; |
188 | < | |
189 | < | |
190 | < | |
191 | < | CGAL::convex_hull_3(points.begin(), points.end(), polyhedron); |
192 | < | |
193 | < | |
194 | < | |
195 | < | Ns_ = polyhedron.size_of_vertices(); |
196 | < | |
197 | < | #ifdef IS_MPI |
198 | < | /* Gather an array of the number of verticies on each processor */ |
199 | < | |
200 | < | |
201 | < | surfacePtsGlobal_.clear(); |
202 | < | surfacePtsLocal_.clear(); |
203 | < | |
204 | < | MPI::COMM_WORLD.Allgather(&Ns_,1,MPI::INT,&NstoProc_[0],1,MPI::INT); |
205 | < | |
206 | < | for (int i = 0; i < nproc_; i++){ |
207 | < | Nsglobal_ += NstoProc_[i]; |
208 | < | } |
209 | < | /*Reminder ideally, we would like to reserve size for the vectors here*/ |
210 | < | surfacePtsLocal_.reserve(Ns_); |
211 | < | surfacePtsGlobal_.resize(Nsglobal_); |
212 | < | // std::fill(surfacePtsGlobal_.begin(),surfacePtsGlobal_.end(),0); |
213 | < | |
214 | < | /* Build a displacements array */ |
215 | < | for (int i = 1; i < nproc_; i++){ |
216 | < | displs_[i] = displs_[i-1] + NstoProc_[i-1]; |
217 | < | } |
218 | < | |
219 | < | int noffset = displs_[myrank_]; |
220 | < | /* gather the potential hull */ |
221 | < | |
222 | < | |
223 | < | for (hc =polyhedron.points_begin();hc != polyhedron.points_end(); ++hc){ |
224 | < | Point_3 mypoint = *hc; |
225 | < | surfacePt_ mpiSurfacePt; |
226 | < | mpiSurfacePt.x = CGAL::to_double(mypoint.x()); |
227 | < | mpiSurfacePt.y = CGAL::to_double(mypoint.y()); |
228 | < | mpiSurfacePt.z = CGAL::to_double(mypoint.z()); |
229 | < | surfacePtsLocal_.push_back(mpiSurfacePt); |
230 | < | } |
231 | < | |
232 | < | MPI::COMM_WORLD.Allgatherv(&surfacePtsLocal_[0],Ns_,surfacePtType,&surfacePtsGlobal_[0],NstoProc_,displs_,surfacePtType); |
233 | < | std::vector<surfacePt_>::iterator spt; |
234 | < | std::vector<Enriched_Point_3> gblpoints; |
235 | < | |
236 | < | int mine = 0; |
237 | < | int pointidx = 0; |
238 | < | for (spt = surfacePtsGlobal_.begin(); spt != surfacePtsGlobal_.end(); ++spt) |
239 | < | { |
240 | < | surfacePt_ thispos = *spt; |
241 | < | Enriched_Point_3 ept(thispos.x,thispos.y,thispos.z); |
242 | < | if (mine >= noffset && mine < noffset + Ns_){ |
243 | < | ept.myPoint(); |
244 | < | ept.setSD(points[pointidx].getStuntDouble()); |
245 | < | pointidx++; |
246 | < | } |
247 | < | gblpoints.push_back(ept); |
248 | < | |
249 | < | mine++; |
250 | < | } |
251 | < | |
252 | < | /* Compute the global hull */ |
253 | < | polyhedron.clear(); |
254 | < | CGAL::convex_hull_3(gblpoints.begin(), gblpoints.end(), polyhedron); |
255 | < | |
256 | < | |
257 | < | #endif |
258 | < | |
259 | < | |
260 | < | |
261 | < | /* Loop over all of the surface triangles and build data structures for atoms and normals*/ |
262 | < | Facet_iterator j; |
263 | < | area_ = 0; |
264 | < | for ( j = polyhedron.facets_begin(); j !=polyhedron.facets_end(); ++j) { |
265 | < | Halfedge_handle h = j->halfedge(); |
266 | < | |
267 | < | Point_3 r0=h->vertex()->point(); |
268 | < | Point_3 r1=h->next()->vertex()->point(); |
269 | < | Point_3 r2=h->next()->next()->vertex()->point(); |
270 | < | |
271 | < | Point_3* pr0 = &r0; |
272 | < | Point_3* pr1 = &r1; |
273 | < | Point_3* pr2 = &r2; |
274 | < | |
275 | < | Enriched_Point_3* er0 = static_cast<Enriched_Point_3*>(pr0); |
276 | < | Enriched_Point_3* er1 = static_cast<Enriched_Point_3*>(pr1); |
277 | < | Enriched_Point_3* er2 = static_cast<Enriched_Point_3*>(pr2); |
278 | < | |
279 | < | // StuntDouble* sd = er0->getStuntDouble(); |
280 | < | std::cerr << "sd globalIndex = " << to_double(er0->x()) << "\n"; |
281 | < | |
282 | < | Point_3 thisCentroid = CGAL::centroid(r0,r1,r2); |
283 | < | |
284 | < | Vector_3 normal = CGAL::cross_product(r1-r0,r2-r0); |
285 | < | |
286 | < | Triangle* face = new Triangle(); |
287 | < | Vector3d V3dNormal(CGAL::to_double(normal.x()),CGAL::to_double(normal.y()),CGAL::to_double(normal.z())); |
288 | < | Vector3d V3dCentroid(CGAL::to_double(thisCentroid.x()),CGAL::to_double(thisCentroid.y()),CGAL::to_double(thisCentroid.z())); |
289 | < | face->setNormal(V3dNormal); |
290 | < | face->setCentroid(V3dCentroid); |
291 | < | RealType faceArea = 0.5*V3dNormal.length(); |
292 | < | face->setArea(faceArea); |
293 | < | area_ += faceArea; |
294 | < | Triangles_.push_back(face); |
295 | < | // ptMapType::const_iterator locn=myMap.find(mypoint); |
296 | < | // int myIndex = locn->second; |
297 | < | |
298 | < | } |
299 | < | |
300 | < | std::cout << "Number of surface atoms is: " << Ns_ << std::endl; |
301 | < | |
302 | < | |
303 | < | |
304 | < | } |
305 | < | void ConvexHull::printHull(const std::string& geomFileName) |
306 | < | { |
307 | < | /* |
308 | < | std::ofstream newGeomFile; |
309 | < | |
310 | < | //create new .md file based on old .md file |
311 | < | newGeomFile.open("testhull.off"); |
312 | < | |
313 | < | // Write polyhedron in Object File Format (OFF). |
314 | < | CGAL::set_ascii_mode( std::cout); |
315 | < | newGeomFile << "OFF" << std::endl << polyhedron.size_of_vertices() << ' ' |
316 | < | << polyhedron.size_of_facets() << " 0" << std::endl; |
317 | < | std::copy( polyhedron.points_begin(), polyhedron.points_end(), |
318 | < | std::ostream_iterator<Point_3>( newGeomFile, "\n")); |
319 | < | for ( Facet_iterator i = polyhedron.facets_begin(); i != polyhedron.facets_end(); ++i) { |
320 | < | Halfedge_facet_circulator j = i->facet_begin(); |
321 | < | // Facets in polyhedral surfaces are at least triangles. |
322 | < | CGAL_assertion( CGAL::circulator_size(j) >= 3); |
323 | < | newGeomFile << CGAL::circulator_size(j) << ' '; |
324 | < | do { |
325 | < | newGeomFile << ' ' << std::distance(polyhedron.vertices_begin(), j->vertex()); |
326 | < | } while ( ++j != i->facet_begin()); |
327 | < | newGeomFile << std::endl; |
328 | < | } |
329 | < | |
330 | < | newGeomFile.close(); |
331 | < | */ |
332 | < | /* |
333 | < | std::ofstream newGeomFile; |
334 | < | |
335 | < | //create new .md file based on old .md file |
336 | < | newGeomFile.open(geomFileName.c_str()); |
337 | < | |
338 | < | // Write polyhedron in Object File Format (OFF). |
339 | < | CGAL::set_ascii_mode( std::cout); |
340 | < | newGeomFile << "OFF" << std::endl << ch_polyhedron.size_of_vertices() << ' ' |
341 | < | << ch_polyhedron.size_of_facets() << " 0" << std::endl; |
342 | < | std::copy( ch_polyhedron.points_begin(), ch_polyhedron.points_end(), |
343 | < | std::ostream_iterator<Point_3>( newGeomFile, "\n")); |
344 | < | for ( Facet_iterator i = ch_polyhedron.facets_begin(); i != ch_polyhedron.facets_end(); ++i) |
345 | < | { |
346 | < | Halfedge_facet_circulator j = i->facet_begin(); |
347 | < | // Facets in polyhedral surfaces are at least triangles. |
348 | < | CGAL_assertion( CGAL::circulator_size(j) >= 3); |
349 | < | newGeomFile << CGAL::circulator_size(j) << ' '; |
350 | < | do |
351 | < | { |
352 | < | newGeomFile << ' ' << std::distance(ch_polyhedron.vertices_begin(), j->vertex()); |
353 | < | } |
354 | < | while ( ++j != i->facet_begin()); |
355 | < | newGeomFile << std::endl; |
356 | < | } |
357 | < | |
358 | < | newGeomFile.close(); |
359 | < | */ |
360 | < | |
361 | < | } |
188 | > | MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0], |
189 | > | &coordsOnProc[0], &vectorDisplacements[0], |
190 | > | MPI::DOUBLE); |
191 | ||
192 | + | MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE, |
193 | + | &globalMasses[0], &hullSitesOnProc[0], |
194 | + | &displacements[0], MPI::DOUBLE); |
195 | ||
196 | < | |
197 | < | |
198 | < | |
199 | < | |
200 | < | |
201 | < | #else |
202 | < | #ifdef HAVE_QHULL |
203 | < | /* Old options Qt Qu Qg QG0 FA */ |
204 | < | ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Qci Tcv Pp"), Ns_(200) { |
205 | < | //If we are doing the mpi version, set up some vectors for data communication |
374 | < | #ifdef IS_MPI |
375 | < | |
376 | < | |
377 | < | nproc_ = MPI::COMM_WORLD.Get_size(); |
378 | < | myrank_ = MPI::COMM_WORLD.Get_rank(); |
379 | < | NstoProc_ = new int[nproc_]; |
380 | < | displs_ = new int[nproc_]; |
381 | < | |
382 | < | // Create a surface point type in MPI to send |
383 | < | //surfacePtType = MPI::DOUBLE.Create_contiguous(3); |
384 | < | // surfacePtType.Commit(); |
385 | < | |
386 | < | |
387 | < | #endif |
388 | < | } |
389 | < | |
390 | < | |
391 | < | |
392 | < | void ConvexHull::computeHull(std::vector<StuntDouble*> bodydoubles) |
393 | < | { |
196 | > | // Free previous hull |
197 | > | qh_freeqhull(!qh_ALL); |
198 | > | qh_memfreeshort(&curlong, &totlong); |
199 | > | if (curlong || totlong) { |
200 | > | sprintf(painCave.errMsg, "ConvexHull: qhull internal warning:\n" |
201 | > | "\tdid not free %d bytes of long memory (%d pieces)", |
202 | > | totlong, curlong); |
203 | > | painCave.isFatal = 1; |
204 | > | simError(); |
205 | > | } |
206 | ||
207 | < | std::vector<int> surfaceIDs; |
208 | < | std::vector<int> surfaceIDsGlobal; |
209 | < | std::vector<int> localPtsMap; |
210 | < | int numpoints = bodydoubles.size(); |
211 | < | |
212 | < | //coordT* pt_array; |
213 | < | coordT* surfpt_array; |
402 | < | vertexT *vertex, **vertexp; |
403 | < | facetT *facet; |
404 | < | setT *vertices; |
405 | < | int curlong,totlong; |
406 | < | int id; |
407 | < | |
408 | < | coordT *point,**pointp; |
409 | < | |
410 | < | |
411 | < | FILE *outdummy = NULL; |
412 | < | FILE *errdummy = NULL; |
413 | < | |
414 | < | //pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_)); |
415 | < | |
416 | < | // double* ptArray = new double[numpoints * 3]; |
417 | < | std::vector<double> ptArray(numpoints*3); |
418 | < | std::vector<bool> isSurfaceID(numpoints); |
419 | < | |
420 | < | // Copy the positon vector into a points vector for qhull. |
421 | < | std::vector<StuntDouble*>::iterator SD; |
422 | < | int i = 0; |
423 | < | for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD) |
424 | < | { |
425 | < | Vector3d pos = (*SD)->getPos(); |
426 | < | |
427 | < | ptArray[dim_ * i] = pos.x(); |
428 | < | ptArray[dim_ * i + 1] = pos.y(); |
429 | < | ptArray[dim_ * i + 2] = pos.z(); |
430 | < | i++; |
431 | < | } |
432 | < | |
433 | < | |
434 | < | |
435 | < | |
436 | < | |
437 | < | |
438 | < | boolT ismalloc = False; |
439 | < | /* Clean up memory from previous convex hull calculations*/ |
440 | < | Triangles_.clear(); |
441 | < | surfaceSDs_.clear(); |
442 | < | surfaceSDs_.reserve(Ns_); |
443 | < | |
444 | < | if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc, |
445 | < | const_cast<char *>(options_.c_str()), NULL, stderr)) { |
446 | < | |
447 | < | sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull"); |
448 | < | painCave.isFatal = 0; |
449 | < | simError(); |
450 | < | |
207 | > | if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc, |
208 | > | const_cast<char *>(options_.c_str()), NULL, stderr)){ |
209 | > | |
210 | > | sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute global convex hull"); |
211 | > | painCave.isFatal = 1; |
212 | > | simError(); |
213 | > | |
214 | } //qh_new_qhull | |
215 | ||
216 | < | |
217 | < | #ifdef IS_MPI |
218 | < | std::vector<double> localPts; |
219 | < | int localPtArraySize; |
220 | < | |
221 | < | |
222 | < | std::fill(isSurfaceID.begin(),isSurfaceID.end(),false); |
223 | < | |
461 | < | |
462 | < | FORALLfacets { |
216 | > | #endif |
217 | > | intPoint = qh interior_point; |
218 | > | RealType calcvol = 0.0; |
219 | > | FORALLfacets { |
220 | > | Triangle face; |
221 | > | //Qhull sets the unit normal in facet->normal |
222 | > | Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]); |
223 | > | face.setUnitNormal(V3dNormal); |
224 | ||
225 | < | if (!facet->simplicial){ |
226 | < | // should never happen with Qt |
466 | < | sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected"); |
467 | < | painCave.isFatal = 0; |
468 | < | simError(); |
469 | < | } |
225 | > | RealType faceArea = qh_facetarea(facet); |
226 | > | face.setArea(faceArea); |
227 | ||
471 | – | |
228 | vertices = qh_facet3vertex(facet); | |
473 | – | FOREACHvertex_(vertices){ |
474 | – | id = qh_pointid(vertex->point); |
475 | – | |
476 | – | if( !isSurfaceID[id] ){ |
477 | – | isSurfaceID[id] = true; |
478 | – | } |
479 | – | } |
480 | – | qh_settempfree(&vertices); |
229 | ||
230 | < | } //FORALLfacets |
230 | > | coordT *center = qh_getcenter(vertices); |
231 | > | Vector3d V3dCentroid(center[0], center[1], center[2]); |
232 | > | face.setCentroid(V3dCentroid); |
233 | ||
234 | < | |
234 | > | Vector3d faceVel = V3Zero; |
235 | > | Vector3d p[3]; |
236 | > | RealType faceMass = 0.0; |
237 | ||
238 | < | /* |
487 | < | std::sort(surfaceIDs.begin(),surfaceIDs.end()); |
488 | < | surfaceIDs.erase(std::unique(surfaceIDs.begin(), surfaceIDs.end()), surfaceIDs.end()); |
489 | < | int localPtArraySize = surfaceIDs.size() * 3; |
490 | < | */ |
238 | > | int ver = 0; |
239 | ||
240 | < | //localPts.resize(localPtArraySize); |
241 | < | //std::fill(localPts.begin(),localPts.end(),0.0); |
240 | > | FOREACHvertex_(vertices){ |
241 | > | int id = qh_pointid(vertex->point); |
242 | > | p[ver][0] = vertex->point[0]; |
243 | > | p[ver][1] = vertex->point[1]; |
244 | > | p[ver][2] = vertex->point[2]; |
245 | > | Vector3d vel; |
246 | > | RealType mass; |
247 | ||
495 | – | |
496 | – | int idx = 0; |
497 | – | int nIsIts = 0; |
498 | – | /* |
499 | – | // Copy the surface points into an array. |
500 | – | for(std::vector<bool>::iterator list_iter = isSurfaceID.begin(); |
501 | – | list_iter != isSurfaceID.end(); list_iter++) |
502 | – | { |
503 | – | bool isIt = *list_iter; |
504 | – | if (isIt){ |
505 | – | localPts.push_back(ptArray[dim_ * idx]); |
506 | – | localPts.push_back(ptArray[dim_ * idx + 1]); |
507 | – | localPts.push_back(ptArray[dim_ * idx + 2]); |
508 | – | localPtsMap.push_back(idx); |
509 | – | nIsIts++; |
510 | – | } //Isit |
511 | – | idx++; |
512 | – | } //isSurfaceID |
513 | – | */ |
514 | – | FORALLvertices { |
515 | – | idx = qh_pointid(vertex->point); |
516 | – | localPts.push_back(ptArray[dim_ * idx]); |
517 | – | localPts.push_back(ptArray[dim_ * idx + 1]); |
518 | – | localPts.push_back(ptArray[dim_ * idx + 2]); |
519 | – | localPtsMap.push_back(idx); |
520 | – | } |
521 | – | |
522 | – | |
523 | – | localPtArraySize = localPts.size(); |
524 | – | |
525 | – | |
526 | – | MPI::COMM_WORLD.Allgather(&localPtArraySize,1,MPI::INT,&NstoProc_[0],1,MPI::INT); |
527 | – | |
528 | – | Nsglobal_=0; |
529 | – | for (int i = 0; i < nproc_; i++){ |
530 | – | Nsglobal_ += NstoProc_[i]; |
531 | – | } |
532 | – | |
533 | – | |
534 | – | int nglobalPts = int(Nsglobal_/3); |
535 | – | |
536 | – | |
537 | – | std::vector<double> globalPts; |
538 | – | globalPts.resize(Nsglobal_); |
539 | – | |
540 | – | isSurfaceID.resize(nglobalPts); |
541 | – | |
542 | – | |
543 | – | std::fill(globalPts.begin(),globalPts.end(),0.0); |
544 | – | |
545 | – | displs_[0] = 0; |
546 | – | /* Build a displacements array */ |
547 | – | for (int i = 1; i < nproc_; i++){ |
548 | – | displs_[i] = displs_[i-1] + NstoProc_[i-1]; |
549 | – | } |
550 | – | |
551 | – | |
552 | – | int noffset = displs_[myrank_]; |
553 | – | /* gather the potential hull */ |
554 | – | |
555 | – | MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE); |
556 | – | |
557 | – | /* |
558 | – | if (myrank_ == 0){ |
559 | – | for (i = 0; i < globalPts.size(); i++){ |
560 | – | std::cout << globalPts[i] << std::endl; |
561 | – | } |
562 | – | } |
563 | – | */ |
564 | – | // Free previous hull |
565 | – | qh_freeqhull(!qh_ALL); |
566 | – | qh_memfreeshort(&curlong, &totlong); |
567 | – | if (curlong || totlong) |
568 | – | std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " |
569 | – | << totlong << curlong << std::endl; |
570 | – | |
571 | – | if (qh_new_qhull(dim_, nglobalPts, &globalPts[0], ismalloc, |
572 | – | const_cast<char *>(options_.c_str()), NULL, stderr)){ |
573 | – | |
574 | – | sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute global convex hull"); |
575 | – | painCave.isFatal = 1; |
576 | – | simError(); |
577 | – | |
578 | – | } //qh_new_qhull |
579 | – | |
580 | – | #endif |
581 | – | |
582 | – | |
583 | – | |
584 | – | |
585 | – | |
586 | – | |
587 | – | unsigned int nf = qh num_facets; |
588 | – | |
589 | – | /* Build Surface SD list first */ |
590 | – | |
591 | – | std::fill(isSurfaceID.begin(),isSurfaceID.end(),false); |
592 | – | |
593 | – | FORALLfacets { |
594 | – | |
595 | – | if (!facet->simplicial){ |
596 | – | // should never happen with Qt |
597 | – | sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected"); |
598 | – | painCave.isFatal = 1; |
599 | – | simError(); |
600 | – | } //simplicical |
601 | – | |
602 | – | Triangle* face = new Triangle(); |
603 | – | Vector3d V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]); |
604 | – | face->setNormal(V3dNormal); |
605 | – | //face->setCentroid(V3dCentroid); |
606 | – | RealType faceArea = 0.5*V3dNormal.length(); |
607 | – | //face->setArea(faceArea); |
608 | – | |
609 | – | |
610 | – | vertices = qh_facet3vertex(facet); |
611 | – | FOREACHvertex_(vertices){ |
612 | – | id = qh_pointid(vertex->point); |
613 | – | int localindex = id; |
248 | #ifdef IS_MPI | |
249 | < | |
250 | < | if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){ |
251 | < | localindex = localPtsMap[id-noffset/3]; |
252 | < | #endif |
619 | < | face->addVertex(bodydoubles[localindex]); |
620 | < | if( !isSurfaceID[id] ){ |
621 | < | isSurfaceID[id] = true; |
622 | < | #ifdef IS_MPI |
623 | < | |
624 | < | #endif |
625 | < | |
626 | < | surfaceSDs_.push_back(bodydoubles[localindex]); |
627 | < | |
628 | < | } //IF isSurfaceID |
249 | > | vel = Vector3d(globalVels[dim_ * id], |
250 | > | globalVels[dim_ * id + 1], |
251 | > | globalVels[dim_ * id + 2]); |
252 | > | mass = globalMasses[id]; |
253 | ||
254 | < | #ifdef IS_MPI |
255 | < | |
632 | < | }else{ |
633 | < | face->addVertex(NULL); |
634 | < | } |
635 | < | #endif |
636 | < | } //Foreachvertex |
254 | > | // localID will be between 0 and hullSitesOnProc[myrank] if we |
255 | > | // own this guy. |
256 | ||
257 | < | Triangles_.push_back(face); |
639 | < | qh_settempfree(&vertices); |
640 | < | |
641 | < | } //FORALLfacets |
257 | > | int localID = id - displacements[myrank]; |
258 | ||
643 | – | |
259 | ||
260 | < | Ns_ = surfaceSDs_.size(); |
260 | > | if (localID >= 0 && localID < hullSitesOnProc[myrank]){ |
261 | > | face.addVertexSD(bodydoubles[indexMap[localID]]); |
262 | > | }else{ |
263 | > | face.addVertexSD(NULL); |
264 | > | } |
265 | > | #else |
266 | > | vel = bodydoubles[id]->getVel(); |
267 | > | mass = bodydoubles[id]->getMass(); |
268 | > | face.addVertexSD(bodydoubles[id]); |
269 | > | #endif |
270 | > | faceVel = faceVel + vel; |
271 | > | faceMass = faceMass + mass; |
272 | > | ver++; |
273 | > | } //Foreachvertex |
274 | ||
275 | + | face.addVertices(p[0], p[1], p[2]); |
276 | + | face.setFacetMass(faceMass); |
277 | + | face.setFacetVelocity(faceVel/3.0); |
278 | + | /* |
279 | + | RealType comparea = face.computeArea(); |
280 | + | realT calcarea = qh_facetarea (facet); |
281 | + | Vector3d V3dCompNorm = -face.computeUnitNormal(); |
282 | + | RealType thisOffset = ((0.0-p[0][0])*V3dCompNorm[0] + (0.0-p[0][1])*V3dCompNorm[1] + (0.0-p[0][2])*V3dCompNorm[2]); |
283 | + | RealType dist = facet->offset + intPoint[0]*V3dNormal[0] + intPoint[1]*V3dNormal[1] + intPoint[2]*V3dNormal[2]; |
284 | + | std::cout << "facet offset and computed offset: " << facet->offset << " " << thisOffset << std::endl; |
285 | + | calcvol += -dist*comparea/qh hull_dim; |
286 | + | */ |
287 | + | Triangles_.push_back(face); |
288 | + | qh_settempfree(&vertices); |
289 | ||
290 | + | } //FORALLfacets |
291 | + | |
292 | qh_getarea(qh facet_list); | |
293 | volume_ = qh totvol; | |
294 | area_ = qh totarea; | |
295 | < | |
652 | < | |
653 | < | |
295 | > | // std::cout << "My volume is: " << calcvol << " qhull volume is:" << volume_ << std::endl; |
296 | qh_freeqhull(!qh_ALL); | |
297 | qh_memfreeshort(&curlong, &totlong); | |
298 | < | if (curlong || totlong) |
299 | < | std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " |
300 | < | << totlong << curlong << std::endl; |
301 | < | |
302 | < | |
303 | < | |
298 | > | if (curlong || totlong) { |
299 | > | sprintf(painCave.errMsg, "ConvexHull: qhull internal warning:\n" |
300 | > | "\tdid not free %d bytes of long memory (%d pieces)", |
301 | > | totlong, curlong); |
302 | > | painCave.isFatal = 1; |
303 | > | simError(); |
304 | > | } |
305 | } | |
306 | ||
307 | + | void ConvexHull::printHull(const std::string& geomFileName) { |
308 | ||
309 | < | |
310 | < | void ConvexHull::printHull(const std::string& geomFileName) |
311 | < | { |
668 | < | |
309 | > | #ifdef IS_MPI |
310 | > | if (worldRank == 0) { |
311 | > | #endif |
312 | FILE *newGeomFile; | |
313 | ||
314 | //create new .md file based on old .md file | |
# | Line 675 | Line 318 | void ConvexHull::printHull(const std::string& geomFile | |
318 | qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); | |
319 | ||
320 | fclose(newGeomFile); | |
321 | + | #ifdef IS_MPI |
322 | + | } |
323 | + | #endif |
324 | } | |
325 | #endif //QHULL | |
680 | – | #endif //CGAL |
681 | – | |
682 | – | |
683 | – |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |