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

Comparing trunk/src/math/AlphaHull.cpp (file contents):
Revision 1938 by gezelter, Thu Oct 31 15:32:17 2013 UTC vs.
Revision 2056 by gezelter, Fri Feb 20 15:12:07 2015 UTC

# Line 112 | Line 112 | void AlphaHull::computeHull(vector<StuntDouble*> bodyd
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    vector<int> hullSitesOnProc(nproc, 0);
# Line 130 | Line 133 | void AlphaHull::computeHull(vector<StuntDouble*> bodyd
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 147 | Line 150 | void AlphaHull::computeHull(vector<StuntDouble*> bodyd
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 170 | Line 173 | void AlphaHull::computeHull(vector<StuntDouble*> bodyd
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);
# Line 210 | Line 213 | void AlphaHull::computeHull(vector<StuntDouble*> bodyd
213    
214    
215    int convexNumVert = qh_setsize(qh_facetvertices (qh facet_list, NULL, false));
216 <  //Insert all the sample points, because, even with alpha=0, the alpha shape/alpha complex will
217 <  //contain them.
216 >  //Insert all the sample points, because, even with alpha=0, the
217 >  //alpha shape/alpha complex will contain them.
218  
219    //  tri::Allocator<CMeshO>::AddVertices(pm.cm,convexNumVert);
220    
# Line 239 | Line 242 | void AlphaHull::computeHull(vector<StuntDouble*> bodyd
242    FORALLfacet_(qh facet_list) {
243      numFacets++;
244      if (!facet->upperdelaunay) {
245 <      //For all facets (that are tetrahedrons)calculate the radius of the empty circumsphere considering
246 <      //the distance between the circumcenter and a vertex of the facet
245 >      //For all facets (that are tetrahedrons)calculate the radius of
246 >      //the empty circumsphere considering the distance between the
247 >      //circumcenter and a vertex of the facet
248        vertexT* vertex = (vertexT *)(facet->vertices->e[0].p);
249        double* center = facet->center;
250        double radius =  qh_pointdist(vertex->point,center,dim_);
251        
252        if (radius>alpha_) // if the facet is not good consider the ridges
253          {
254 <          //if calculating the alphashape, unmark the facet ('good' is used as 'marked').
254 >          //if calculating the alphashape, unmark the facet ('good' is
255 >          //used as 'marked').
256            facet->good=false;
257            
258 <          //Compute each ridge (triangle) once and test the cironference radius with alpha
258 >          //Compute each ridge (triangle) once and test the
259 >          //cironference radius with alpha
260            facet->visitid= qh visit_id;
261            qh_makeridges(facet);
262            ridgeT *ridge, **ridgep;
# Line 273 | Line 279 | void AlphaHull::computeHull(vector<StuntDouble*> bodyd
279              }
280            }
281  
282 <          //If calculating the alphashape, mark the facet('good' is used as 'marked').
283 <          //This facet will have some triangles hidden by the facet's neighbor.
282 >          //If calculating the alphashape, mark the facet('good' is
283 >          //used as 'marked').  This facet will have some triangles
284 >          //hidden by the facet's neighbor.
285            if(goodTriangles==4)
286              facet->good=true;
287            
288          }
289 <      else //the facet is good. Put all the triangles of the tetrahedron in the mesh
289 >      else //the facet is good. Put all the triangles of the
290 >           //tetrahedron in the mesh
291          {
292            //Compute each ridge (triangle) once
293            facet->visitid= qh visit_id;
294 <          //If calculating the alphashape, mark the facet('good' is used as 'marked').
295 <          //This facet will have some triangles hidden by the facet's neighbor.
294 >          //If calculating the alphashape, mark the facet('good' is
295 >          //used as 'marked').  This facet will have some triangles
296 >          //hidden by the facet's neighbor.
297            facet->good=true;
298            qh_makeridges(facet);
299            ridgeT *ridge, **ridgep;
# Line 299 | Line 308 | void AlphaHull::computeHull(vector<StuntDouble*> bodyd
308    }
309    //assert(numFacets== qh num_facets);
310    
311 <  //Filter the triangles (only the ones on the boundary of the alpha complex) and build the mesh
311 >  //Filter the triangles (only the ones on the boundary of the alpha
312 >  //complex) and build the mesh
313  
314    int ridgesCount=0;
315    
# Line 356 | Line 366 | void AlphaHull::computeHull(vector<StuntDouble*> bodyd
366        Vector3d normal = face.getUnitNormal();
367        RealType offset =  ((0.0-p[0][0])*normal[0] + (0.0-p[0][1])*normal[1] + (0.0-p[0][2])*normal[2]);
368        RealType dist =  normal[0] * interiorPoint[0] + normal[1]*interiorPoint[1] + normal[2]*interiorPoint[2];
359 –     cout << "Dist and normal and area are: " << normal << endl;
369        volume_ += dist *area/qh hull_dim;
370        
371        Triangles_.push_back(face);
372      }
373    }
374  
366 –  cout << "Volume is: " << volume_ << endl;
375  
376   //assert(pm.cm.fn == ridgesCount);
377   /*

Diff Legend

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