# | Line 1 | Line 1 | |
---|---|---|
1 | < | /* Copyright (c) 2008, 2009, 2010 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 | |
# | Line 32 | Line 32 | |
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). |
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, 234107 (2008). |
38 | * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). | |
39 | < | * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). * |
39 | > | * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
40 | * | |
41 | * AlphaHull.cpp | |
42 | * | |
43 | < | * Purpose: To calculate Alpha hull, hull volume libqhull. |
44 | < | * |
45 | < | * Created by Charles F. Vardeman II on 11 Dec 2006. |
46 | < | * @author Charles F. Vardeman II |
47 | < | * @version $Id$ |
48 | < | * |
43 | > | * Purpose: To calculate an alpha-shape hull. |
44 | */ | |
45 | ||
46 | + | #ifdef IS_MPI |
47 | + | #include <mpi.h> |
48 | + | #endif |
49 | + | |
50 | /* Standard includes independent of library */ | |
51 | ||
52 | #include <iostream> | |
# | Line 55 | Line 54 | |
54 | #include <list> | |
55 | #include <algorithm> | |
56 | #include <iterator> | |
58 | – | #include <utility> |
57 | #include "math/AlphaHull.hpp" | |
58 | #include "utils/simError.h" | |
59 | < | #ifdef IS_MPI |
62 | < | #include <mpi.h> |
63 | < | #endif |
59 | > | |
60 | #include "math/qhull.hpp" | |
61 | ||
62 | + | #ifdef HAVE_QHULL |
63 | + | using namespace std; |
64 | using namespace OpenMD; | |
65 | ||
66 | < | #ifdef HAVE_QHULL |
69 | < | double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim); |
66 | > | double calculate_circumradius(pointT* p0, pointT* p1, pointT* p2, int dim); |
67 | ||
68 | < | AlphaHull::AlphaHull(double alpha) : Hull(), dim_(3), alpha_(alpha), options_("qhull d QJ Tcv Pp") { |
68 | > | AlphaHull::AlphaHull(double alpha) : Hull(), dim_(3), alpha_(alpha), |
69 | > | options_("qhull d QJ Tcv Pp") { |
70 | } | |
71 | ||
72 | < | void AlphaHull::computeHull(std::vector<StuntDouble*> bodydoubles) { |
72 | > | void AlphaHull::computeHull(vector<StuntDouble*> bodydoubles) { |
73 | ||
74 | int numpoints = bodydoubles.size(); | |
77 | – | bool alphashape=true; |
75 | ||
76 | Triangles_.clear(); | |
77 | ||
78 | < | vertexT *vertex, **vertexp; |
78 | > | vertexT *vertex; |
79 | facetT *facet, *neighbor; | |
83 | – | setT *vertices, *verticestop, *verticesbottom; |
84 | – | int curlong, totlong; |
80 | pointT *interiorPoint; | |
81 | + | int curlong, totlong; |
82 | ||
83 | < | std::vector<double> ptArray(numpoints*dim_); |
84 | < | |
83 | > | |
84 | > | vector<double> ptArray(numpoints*dim_); |
85 | > | |
86 | // Copy the positon vector into a points vector for qhull. | |
87 | < | std::vector<StuntDouble*>::iterator SD; |
87 | > | vector<StuntDouble*>::iterator SD; |
88 | int i = 0; | |
89 | for (SD =bodydoubles.begin(); SD != bodydoubles.end(); ++SD){ | |
90 | Vector3d pos = (*SD)->getPos(); | |
# | Line 96 | Line 93 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
93 | ptArray[dim_ * i + 2] = pos.z(); | |
94 | i++; | |
95 | } | |
96 | < | |
97 | < | /* Clean up memory from previous convex hull calculations*/ |
96 | > | |
97 | > | /* Clean up memory from previous convex hull calculations*/ |
98 | boolT ismalloc = False; | |
99 | < | |
100 | < | int ridgesCount=0; |
99 | > | |
100 | > | /* compute the hull for our local points (or all the points for single |
101 | > | processor versions) */ |
102 | if (qh_new_qhull(dim_, numpoints, &ptArray[0], ismalloc, | |
103 | const_cast<char *>(options_.c_str()), NULL, stderr)) { | |
104 | < | |
104 | > | |
105 | sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute convex hull"); | |
106 | painCave.isFatal = 1; | |
107 | simError(); | |
108 | ||
109 | } //qh_new_qhull | |
110 | < | |
111 | < | |
110 | > | |
111 | > | |
112 | #ifdef IS_MPI | |
113 | //If we are doing the mpi version, set up some vectors for data communication | |
114 | ||
115 | < | int nproc = MPI::COMM_WORLD.Get_size(); |
116 | < | int myrank = MPI::COMM_WORLD.Get_rank(); |
115 | > | int nproc; |
116 | > | int myrank; |
117 | > | MPI_Comm_size( MPI_COMM_WORLD, &nproc); |
118 | > | MPI_Comm_rank( MPI_COMM_WORLD, &myrank); |
119 | > | |
120 | int localHullSites = 0; | |
121 | ||
122 | < | std::vector<int> hullSitesOnProc(nproc, 0); |
123 | < | std::vector<int> coordsOnProc(nproc, 0); |
124 | < | std::vector<int> displacements(nproc, 0); |
125 | < | std::vector<int> vectorDisplacements(nproc, 0); |
122 | > | vector<int> hullSitesOnProc(nproc, 0); |
123 | > | vector<int> coordsOnProc(nproc, 0); |
124 | > | vector<int> displacements(nproc, 0); |
125 | > | vector<int> vectorDisplacements(nproc, 0); |
126 | ||
127 | < | std::vector<double> coords; |
128 | < | std::vector<double> vels; |
129 | < | std::vector<int> indexMap; |
130 | < | std::vector<double> masses; |
127 | > | vector<double> coords; |
128 | > | vector<double> vels; |
129 | > | vector<int> indexMap; |
130 | > | vector<double> masses; |
131 | ||
132 | FORALLvertices{ | |
133 | localHullSites++; | |
134 | ||
135 | int idx = qh_pointid(vertex->point); | |
136 | < | |
136 | > | |
137 | indexMap.push_back(idx); | |
138 | ||
139 | coords.push_back(ptArray[dim_ * idx]); | |
# | Line 149 | Line 150 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
150 | masses.push_back(sd->getMass()); | |
151 | } | |
152 | ||
153 | < | MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], |
154 | < | 1, MPI::INT); |
153 | > | MPI_Allgather(&localHullSites, 1, MPI_INT, &hullSitesOnProc[0], |
154 | > | 1, MPI_INT, MPI_COMM_WORLD); |
155 | ||
156 | int globalHullSites = 0; | |
157 | for (int iproc = 0; iproc < nproc; iproc++){ | |
# | Line 166 | Line 167 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
167 | vectorDisplacements[iproc] = vectorDisplacements[iproc-1] + coordsOnProc[iproc-1]; | |
168 | } | |
169 | ||
170 | < | std::vector<double> globalCoords(dim_ * globalHullSites); |
171 | < | std::vector<double> globalVels(dim_ * globalHullSites); |
172 | < | std::vector<double> globalMasses(globalHullSites); |
170 | > | vector<double> globalCoords(dim_ * globalHullSites); |
171 | > | vector<double> globalVels(dim_ * globalHullSites); |
172 | > | vector<double> globalMasses(globalHullSites); |
173 | ||
174 | int count = coordsOnProc[myrank]; | |
175 | ||
176 | < | MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0], |
177 | < | &coordsOnProc[0], &vectorDisplacements[0], |
178 | < | MPI::DOUBLE); |
179 | < | |
180 | < | MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0], |
181 | < | &coordsOnProc[0], &vectorDisplacements[0], |
182 | < | MPI::DOUBLE); |
183 | < | |
184 | < | MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE, |
185 | < | &globalMasses[0], &hullSitesOnProc[0], |
186 | < | &displacements[0], MPI::DOUBLE); |
187 | < | |
176 | > | MPI_Allgatherv(&coords[0], count, MPI_DOUBLE, &globalCoords[0], |
177 | > | &coordsOnProc[0], &vectorDisplacements[0], |
178 | > | MPI_DOUBLE, MPI_COMM_WORLD); |
179 | > | |
180 | > | MPI_Allgatherv(&vels[0], count, MPI_DOUBLE, &globalVels[0], |
181 | > | &coordsOnProc[0], &vectorDisplacements[0], |
182 | > | MPI_DOUBLE, MPI_COMM_WORLD); |
183 | > | |
184 | > | MPI_Allgatherv(&masses[0], localHullSites, MPI_DOUBLE, |
185 | > | &globalMasses[0], &hullSitesOnProc[0], |
186 | > | &displacements[0], MPI_DOUBLE, MPI_COMM_WORLD); |
187 | > | |
188 | // Free previous hull | |
189 | qh_freeqhull(!qh_ALL); | |
190 | qh_memfreeshort(&curlong, &totlong); | |
191 | < | if (curlong || totlong) |
192 | < | std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " |
193 | < | << totlong << curlong << std::endl; |
191 | > | if (curlong || totlong) { |
192 | > | sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n" |
193 | > | "\tdid not free %d bytes of long memory (%d pieces)", |
194 | > | totlong, curlong); |
195 | > | painCave.isFatal = 1; |
196 | > | simError(); |
197 | > | } |
198 | ||
199 | if (qh_new_qhull(dim_, globalHullSites, &globalCoords[0], ismalloc, | |
200 | const_cast<char *>(options_.c_str()), NULL, stderr)){ | |
201 | ||
202 | < | sprintf(painCave.errMsg, "AlphaHull: Qhull failed to compute global convex hull"); |
202 | > | sprintf(painCave.errMsg, |
203 | > | "AlphaHull: Qhull failed to compute global convex hull"); |
204 | painCave.isFatal = 1; | |
205 | simError(); | |
206 | < | |
206 | > | |
207 | } //qh_new_qhull | |
202 | – | |
208 | ||
209 | #endif | |
210 | ||
# | Line 232 | Line 237 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
237 | ||
238 | qh visit_id++; | |
239 | int numFacets=0; | |
240 | < | std::vector<std::vector <int> > facetlist; |
240 | > | vector<vector <int> > facetlist; |
241 | interiorPoint = qh interior_point; | |
242 | FORALLfacet_(qh facet_list) { | |
243 | numFacets++; | |
# | Line 299 | Line 304 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
304 | ||
305 | //Filter the triangles (only the ones on the boundary of the alpha complex) and build the mesh | |
306 | ||
307 | < | |
307 | > | int ridgesCount=0; |
308 | ||
309 | ridgeT *ridge, **ridgep; | |
310 | FOREACHridge_(set) { | |
# | Line 314 | Line 319 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
319 | ||
320 | ||
321 | //coordT *center = qh_getcenter(ridge->vertices); | |
322 | < | //std::cout << "Centers are " << center[0] << " " <<center[1] << " " << center[2] << std::endl; |
322 | > | //cout << "Centers are " << center[0] << " " <<center[1] << " " << center[2] << endl; |
323 | //Vector3d V3dCentroid(center[0], center[1], center[2]); | |
324 | //face.setCentroid(V3dCentroid); | |
325 | ||
# | Line 324 | Line 329 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
329 | RealType faceMass = 0.0; | |
330 | ||
331 | int ver = 0; | |
332 | < | std::vector<int> virtexlist; |
332 | > | vector<int> virtexlist; |
333 | FOREACHvertex_i_(ridge->vertices){ | |
334 | int id = qh_pointid(vertex->point); | |
335 | p[ver][0] = vertex->point[0]; | |
# | Line 334 | Line 339 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
339 | RealType mass; | |
340 | ver++; | |
341 | virtexlist.push_back(id); | |
342 | < | // std::cout << "Ridge: " << ridgesCount << " Vertex " << id << std::endl; |
342 | > | // cout << "Ridge: " << ridgesCount << " Vertex " << id << endl; |
343 | ||
344 | vel = bodydoubles[id]->getVel(); | |
345 | mass = bodydoubles[id]->getMass(); | |
# | Line 354 | Line 359 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
359 | Vector3d normal = face.getUnitNormal(); | |
360 | RealType offset = ((0.0-p[0][0])*normal[0] + (0.0-p[0][1])*normal[1] + (0.0-p[0][2])*normal[2]); | |
361 | RealType dist = normal[0] * interiorPoint[0] + normal[1]*interiorPoint[1] + normal[2]*interiorPoint[2]; | |
362 | < | std::cout << "Dist and normal and area are: " << normal << std::endl; |
362 | > | cout << "Dist and normal and area are: " << normal << endl; |
363 | volume_ += dist *area/qh hull_dim; | |
364 | ||
365 | Triangles_.push_back(face); | |
366 | } | |
367 | } | |
368 | ||
369 | < | std::cout << "Volume is: " << volume_ << std::endl; |
369 | > | cout << "Volume is: " << volume_ << endl; |
370 | ||
371 | //assert(pm.cm.fn == ridgesCount); | |
372 | /* | |
# | Line 457 | Line 462 | void AlphaHull::computeHull(std::vector<StuntDouble*> | |
462 | ||
463 | qh_freeqhull(!qh_ALL); | |
464 | qh_memfreeshort(&curlong, &totlong); | |
465 | < | if (curlong || totlong) |
466 | < | std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " |
467 | < | << totlong << curlong << std::endl; |
465 | > | if (curlong || totlong) { |
466 | > | sprintf(painCave.errMsg, "AlphaHull: qhull internal warning:\n" |
467 | > | "\tdid not free %d bytes of long memory (%d pieces)", |
468 | > | totlong, curlong); |
469 | > | painCave.isFatal = 1; |
470 | > | simError(); |
471 | > | } |
472 | } | |
473 | ||
474 | < | void AlphaHull::printHull(const std::string& geomFileName) { |
475 | < | |
474 | > | void AlphaHull::printHull(const string& geomFileName) { |
475 | > | |
476 | #ifdef IS_MPI | |
477 | if (worldRank == 0) { | |
478 | #endif | |
479 | < | FILE *newGeomFile; |
480 | < | |
481 | < | //create new .md file based on old .md file |
482 | < | newGeomFile = fopen(geomFileName.c_str(), "w"); |
483 | < | qh_findgood_all(qh facet_list); |
484 | < | for (int i = 0; i < qh_PRINTEND; i++) |
485 | < | qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); |
486 | < | |
487 | < | fclose(newGeomFile); |
479 | > | FILE *newGeomFile; |
480 | > | |
481 | > | //create new .md file based on old .md file |
482 | > | newGeomFile = fopen(geomFileName.c_str(), "w"); |
483 | > | qh_findgood_all(qh facet_list); |
484 | > | for (int i = 0; i < qh_PRINTEND; i++) |
485 | > | qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); |
486 | > | |
487 | > | fclose(newGeomFile); |
488 | #ifdef IS_MPI | |
489 | } | |
490 | #endif |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |