# | Line 1 | Line 1 | |
---|---|---|
1 | < | /* Copyright (c) 2006 The University of Notre Dame. All Rights Reserved. |
1 | > | /* Copyright (c) 2008, 2009 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 | |
# | 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.4 2007-05-30 19:51:07 chuckv Exp $ |
47 | > | * @version $Id: ConvexHull.cpp,v 1.19 2009-10-21 15:48:12 gezelter Exp $ |
48 | * | |
49 | */ | |
50 | ||
51 | + | /* Standard includes independent of library */ |
52 | + | |
53 | #include <iostream> | |
54 | #include <fstream> | |
55 | + | #include <list> |
56 | + | #include <algorithm> |
57 | + | #include <iterator> |
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; | |
66 | ||
67 | < | ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt FA") { |
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 | < | bool ConvexHull::genHull(std::vector<Vector3d> pos) { |
81 | < | FILE *outfile = stdout; |
82 | < | FILE *errfile = stderr; |
80 | > | /* Old options Qt Qu Qg QG0 FA */ |
81 | > | /* More old opts Qc Qi Pp*/ |
82 | > | |
83 | > | ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp") { |
84 | > | } |
85 | > | |
86 | > | void ConvexHull::computeHull(std::vector<StuntDouble*> bodydoubles) { |
87 | > | |
88 | > | int numpoints = bodydoubles.size(); |
89 | > | |
90 | > | Triangles_.clear(); |
91 | > | |
92 | > | vertexT *vertex, **vertexp; |
93 | facetT *facet; | |
94 | < | int exitcode; |
95 | < | boolT ismalloc = False; |
66 | < | int curlong,totlong; |
94 | > | setT *vertices; |
95 | > | int curlong, totlong; |
96 | ||
97 | < | int numpoints = pos.size(); |
69 | < | |
70 | < | coordT* points; |
71 | < | points = (coordT*) malloc(sizeof(coordT) * (numpoints*dim_)); |
97 | > | std::vector<double> ptArray(numpoints*dim_); |
98 | ||
99 | < | for (int i=0; i<numpoints; i++) { |
100 | < | points[dim_ * i] = pos[i][0]; |
101 | < | points[dim_ * i + 1] = pos[i][1]; |
102 | < | points[dim_ * i + 2] = pos[i][2]; |
99 | > | // Copy the positon vector into a points vector for qhull. |
100 | > | std::vector<StuntDouble*>::iterator SD; |
101 | > | int i = 0; |
102 | > | for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){ |
103 | > | Vector3d pos = (*SD)->getPos(); |
104 | > | ptArray[dim_ * i] = pos.x(); |
105 | > | ptArray[dim_ * i + 1] = pos.y(); |
106 | > | ptArray[dim_ * i + 2] = pos.z(); |
107 | > | i++; |
108 | } | |
109 | ||
110 | + | boolT ismalloc = False; |
111 | + | /* Clean up memory from previous convex hull calculations*/ |
112 | ||
113 | < | qh_initflags (const_cast<char *>(options_.c_str())); |
114 | < | qh_init_B (points, numpoints, dim_, ismalloc); |
115 | < | qh_qhull(); |
116 | < | qh_check_output(); |
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 | > | |
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 | + | int* hullSitesOnProc = new int[nproc]; |
130 | + | int* coordsOnProc = new int[nproc]; |
131 | + | int* displacements = new int[nproc]; |
132 | + | int* vectorDisplacements = new int[nproc]; |
133 | + | |
134 | + | std::vector<double> coords; |
135 | + | std::vector<double> vels; |
136 | + | std::vector<int> indexMap; |
137 | + | std::vector<double> masses; |
138 | + | |
139 | + | FORALLvertices{ |
140 | + | localHullSites++; |
141 | + | |
142 | + | int idx = qh_pointid(vertex->point); |
143 | + | |
144 | + | indexMap.push_back(idx); |
145 | + | |
146 | + | coords.push_back(ptArray[dim_ * idx]); |
147 | + | coords.push_back(ptArray[dim_ * idx + 1]); |
148 | + | coords.push_back(ptArray[dim_ * idx + 2]); |
149 | + | |
150 | + | StuntDouble* sd = bodydoubles[idx]; |
151 | + | |
152 | + | Vector3d vel = sd->getVel(); |
153 | + | vels.push_back(vel.x()); |
154 | + | vels.push_back(vel.y()); |
155 | + | vels.push_back(vel.z()); |
156 | + | |
157 | + | masses.push_back(sd->getMass()); |
158 | + | } |
159 | + | |
160 | + | MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], |
161 | + | 1, MPI::INT); |
162 | + | |
163 | + | int globalHullSites = 0; |
164 | + | for (int iproc = 0; iproc < nproc; iproc++){ |
165 | + | globalHullSites += hullSitesOnProc[iproc]; |
166 | + | coordsOnProc[iproc] = dim_ * hullSitesOnProc[iproc]; |
167 | + | } |
168 | + | |
169 | + | displacements[0] = 0; |
170 | + | vectorDisplacements[0] = 0; |
171 | ||
172 | + | for (int iproc = 1; iproc < nproc; iproc++){ |
173 | + | displacements[iproc] = displacements[iproc-1] + hullSitesOnProc[iproc-1]; |
174 | + | vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1]; |
175 | + | } |
176 | + | |
177 | + | std::vector<double> globalCoords(dim_ * globalHullSites); |
178 | + | std::vector<double> globalVels(dim_ * globalHullSites); |
179 | + | std::vector<double> globalMasses(globalHullSites); |
180 | + | |
181 | + | int count = coordsOnProc[myrank]; |
182 | ||
183 | < | qh_getarea(qh facet_list); |
184 | < | volume_ = qh totvol; |
185 | < | area_ = qh totarea; |
186 | < | |
187 | < | |
188 | < | |
183 | > | MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0], |
184 | > | &coordsOnProc[0], &vectorDisplacements[0], |
185 | > | MPI::DOUBLE); |
186 | > | |
187 | > | MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0], |
188 | > | &coordsOnProc[0], &vectorDisplacements[0], |
189 | > | MPI::DOUBLE); |
190 | > | |
191 | > | MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE, |
192 | > | &globalMasses[0], &hullSitesOnProc[0], |
193 | > | &displacements[0], MPI::DOUBLE); |
194 | > | |
195 | > | // Free previous hull |
196 | qh_freeqhull(!qh_ALL); | |
197 | < | qh_memfreeshort (&curlong, &totlong); |
197 | > | qh_memfreeshort(&curlong, &totlong); |
198 | if (curlong || totlong) | |
199 | < | fprintf (errfile, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n", |
200 | < | totlong, curlong); |
199 | > | std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " |
200 | > | << totlong << curlong << std::endl; |
201 | ||
202 | < | |
203 | < | return true; |
204 | < | } |
202 | > | if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc, |
203 | > | const_cast<char *>(options_.c_str()), NULL, stderr)){ |
204 | > | |
205 | > | sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute global convex hull"); |
206 | > | painCave.isFatal = 1; |
207 | > | simError(); |
208 | > | |
209 | > | } //qh_new_qhull |
210 | ||
211 | + | #endif |
212 | ||
213 | < | RealType ConvexHull::getVolume() { |
214 | < | return volume_; |
106 | < | } |
213 | > | FORALLfacets { |
214 | > | Triangle face; |
215 | ||
216 | < | void ConvexHull::geomviewHull(const std::string& geomFileName) { |
216 | > | Vector3d V3dNormal(facet->normal[0], facet->normal[1], facet->normal[2]); |
217 | > | face.setNormal(V3dNormal); |
218 | > | |
219 | > | RealType faceArea = qh_facetarea(facet); |
220 | > | face.setArea(faceArea); |
221 | > | |
222 | > | vertices = qh_facet3vertex(facet); |
223 | > | |
224 | > | coordT *center = qh_getcenter(vertices); |
225 | > | Vector3d V3dCentroid(center[0], center[1], center[2]); |
226 | > | face.setCentroid(V3dCentroid); |
227 | > | |
228 | > | Vector3d faceVel = V3Zero; |
229 | > | Vector3d p[3]; |
230 | > | RealType faceMass = 0.0; |
231 | > | |
232 | > | int ver = 0; |
233 | > | |
234 | > | FOREACHvertex_(vertices){ |
235 | > | int id = qh_pointid(vertex->point); |
236 | > | p[ver][0] = vertex->point[0]; |
237 | > | p[ver][1] = vertex->point[1]; |
238 | > | p[ver][2] = vertex->point[2]; |
239 | > | |
240 | > | Vector3d vel; |
241 | > | RealType mass; |
242 | > | |
243 | > | #ifdef IS_MPI |
244 | > | vel = Vector3d(globalVels[dim_ * id], |
245 | > | globalVels[dim_ * id + 1], |
246 | > | globalVels[dim_ * id + 2]); |
247 | > | mass = globalMasses[id]; |
248 | > | |
249 | > | // localID will be between 0 and hullSitesOnProc[myrank] if we |
250 | > | // own this guy. |
251 | > | |
252 | > | int localID = id - displacements[myrank]; |
253 | > | |
254 | > | if (localID >= 0 && localID < hullSitesOnProc[myrank]) |
255 | > | face.addVertexSD(bodydoubles[indexMap[localID]]); |
256 | > | |
257 | > | #else |
258 | > | vel = bodydoubles[id]->getVel(); |
259 | > | mass = bodydoubles[id]->getMass(); |
260 | > | face.addVertexSD(bodydoubles[id]); |
261 | > | #endif |
262 | > | |
263 | > | faceVel = faceVel + vel; |
264 | > | faceMass = faceMass + mass; |
265 | > | ver++; |
266 | > | } //Foreachvertex |
267 | > | |
268 | > | face.addVertices(p[0], p[1], p[2]); |
269 | > | face.setFacetMass(faceMass); |
270 | > | face.setFacetVelocity(faceVel/3.0); |
271 | > | Triangles_.push_back(face); |
272 | > | qh_settempfree(&vertices); |
273 | > | |
274 | > | } //FORALLfacets |
275 | ||
276 | < | std::ofstream newGeomFile; |
276 | > | qh_getarea(qh facet_list); |
277 | > | volume_ = qh totvol; |
278 | > | area_ = qh totarea; |
279 | > | |
280 | > | #ifdef IS_MPI |
281 | > | delete [] hullSitesOnProc; |
282 | > | delete [] coordsOnProc; |
283 | > | delete [] displacements; |
284 | > | delete [] vectorDisplacements; |
285 | > | #endif |
286 | ||
287 | + | qh_freeqhull(!qh_ALL); |
288 | + | qh_memfreeshort(&curlong, &totlong); |
289 | + | if (curlong || totlong) |
290 | + | std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " |
291 | + | << totlong << curlong << std::endl; |
292 | + | } |
293 | + | |
294 | + | void ConvexHull::printHull(const std::string& geomFileName) { |
295 | + | FILE *newGeomFile; |
296 | + | |
297 | //create new .md file based on old .md file | |
298 | < | newGeomFile.open(geomFileName.c_str()); |
298 | > | newGeomFile = fopen(geomFileName.c_str(), "w"); |
299 | > | qh_findgood_all(qh facet_list); |
300 | > | for (int i = 0; i < qh_PRINTEND; i++) |
301 | > | qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); |
302 | ||
303 | < | |
116 | < | newGeomFile.close(); |
303 | > | fclose(newGeomFile); |
304 | } | |
305 | + | #endif //QHULL |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |