--- trunk/src/math/AlphaHull.cpp 2013/10/31 15:32:17 1938 +++ trunk/src/math/AlphaHull.cpp 2015/03/05 16:30:23 2069 @@ -112,8 +112,11 @@ void AlphaHull::computeHull(vector bodyd #ifdef IS_MPI //If we are doing the mpi version, set up some vectors for data communication - int nproc = MPI::COMM_WORLD.Get_size(); - int myrank = MPI::COMM_WORLD.Get_rank(); + int nproc; + int myrank; + MPI_Comm_size( MPI_COMM_WORLD, &nproc); + MPI_Comm_rank( MPI_COMM_WORLD, &myrank); + int localHullSites = 0; vector hullSitesOnProc(nproc, 0); @@ -130,7 +133,7 @@ void AlphaHull::computeHull(vector bodyd localHullSites++; int idx = qh_pointid(vertex->point); - + indexMap.push_back(idx); coords.push_back(ptArray[dim_ * idx]); @@ -147,8 +150,8 @@ void AlphaHull::computeHull(vector bodyd masses.push_back(sd->getMass()); } - MPI::COMM_WORLD.Allgather(&localHullSites, 1, MPI::INT, &hullSitesOnProc[0], - 1, MPI::INT); + MPI_Allgather(&localHullSites, 1, MPI_INT, &hullSitesOnProc[0], + 1, MPI_INT, MPI_COMM_WORLD); int globalHullSites = 0; for (int iproc = 0; iproc < nproc; iproc++){ @@ -170,18 +173,18 @@ void AlphaHull::computeHull(vector bodyd int count = coordsOnProc[myrank]; - MPI::COMM_WORLD.Allgatherv(&coords[0], count, MPI::DOUBLE, &globalCoords[0], - &coordsOnProc[0], &vectorDisplacements[0], - MPI::DOUBLE); - - MPI::COMM_WORLD.Allgatherv(&vels[0], count, MPI::DOUBLE, &globalVels[0], - &coordsOnProc[0], &vectorDisplacements[0], - MPI::DOUBLE); - - MPI::COMM_WORLD.Allgatherv(&masses[0], localHullSites, MPI::DOUBLE, - &globalMasses[0], &hullSitesOnProc[0], - &displacements[0], MPI::DOUBLE); - + MPI_Allgatherv(&coords[0], count, MPI_DOUBLE, &globalCoords[0], + &coordsOnProc[0], &vectorDisplacements[0], + MPI_DOUBLE, MPI_COMM_WORLD); + + MPI_Allgatherv(&vels[0], count, MPI_DOUBLE, &globalVels[0], + &coordsOnProc[0], &vectorDisplacements[0], + MPI_DOUBLE, MPI_COMM_WORLD); + + MPI_Allgatherv(&masses[0], localHullSites, MPI_DOUBLE, + &globalMasses[0], &hullSitesOnProc[0], + &displacements[0], MPI_DOUBLE, MPI_COMM_WORLD); + // Free previous hull qh_freeqhull(!qh_ALL); qh_memfreeshort(&curlong, &totlong); @@ -209,9 +212,9 @@ void AlphaHull::computeHull(vector bodyd qh_setvoronoi_all(); - int convexNumVert = qh_setsize(qh_facetvertices (qh facet_list, NULL, false)); - //Insert all the sample points, because, even with alpha=0, the alpha shape/alpha complex will - //contain them. + // int convexNumVert = qh_setsize(qh_facetvertices (qh facet_list, NULL, false)); + //Insert all the sample points, because, even with alpha=0, the + //alpha shape/alpha complex will contain them. // tri::Allocator::AddVertices(pm.cm,convexNumVert); @@ -239,18 +242,21 @@ void AlphaHull::computeHull(vector bodyd FORALLfacet_(qh facet_list) { numFacets++; if (!facet->upperdelaunay) { - //For all facets (that are tetrahedrons)calculate the radius of the empty circumsphere considering - //the distance between the circumcenter and a vertex of the facet + //For all facets (that are tetrahedrons)calculate the radius of + //the empty circumsphere considering the distance between the + //circumcenter and a vertex of the facet vertexT* vertex = (vertexT *)(facet->vertices->e[0].p); double* center = facet->center; double radius = qh_pointdist(vertex->point,center,dim_); if (radius>alpha_) // if the facet is not good consider the ridges { - //if calculating the alphashape, unmark the facet ('good' is used as 'marked'). + //if calculating the alphashape, unmark the facet ('good' is + //used as 'marked'). facet->good=false; - //Compute each ridge (triangle) once and test the cironference radius with alpha + //Compute each ridge (triangle) once and test the + //cironference radius with alpha facet->visitid= qh visit_id; qh_makeridges(facet); ridgeT *ridge, **ridgep; @@ -273,18 +279,21 @@ void AlphaHull::computeHull(vector bodyd } } - //If calculating the alphashape, mark the facet('good' is used as 'marked'). - //This facet will have some triangles hidden by the facet's neighbor. + //If calculating the alphashape, mark the facet('good' is + //used as 'marked'). This facet will have some triangles + //hidden by the facet's neighbor. if(goodTriangles==4) facet->good=true; } - else //the facet is good. Put all the triangles of the tetrahedron in the mesh + else //the facet is good. Put all the triangles of the + //tetrahedron in the mesh { //Compute each ridge (triangle) once facet->visitid= qh visit_id; - //If calculating the alphashape, mark the facet('good' is used as 'marked'). - //This facet will have some triangles hidden by the facet's neighbor. + //If calculating the alphashape, mark the facet('good' is + //used as 'marked'). This facet will have some triangles + //hidden by the facet's neighbor. facet->good=true; qh_makeridges(facet); ridgeT *ridge, **ridgep; @@ -299,7 +308,8 @@ void AlphaHull::computeHull(vector bodyd } //assert(numFacets== qh num_facets); - //Filter the triangles (only the ones on the boundary of the alpha complex) and build the mesh + //Filter the triangles (only the ones on the boundary of the alpha + //complex) and build the mesh int ridgesCount=0; @@ -354,16 +364,14 @@ void AlphaHull::computeHull(vector bodyd RealType area = face.getArea(); area_ += area; Vector3d normal = face.getUnitNormal(); - RealType offset = ((0.0-p[0][0])*normal[0] + (0.0-p[0][1])*normal[1] + (0.0-p[0][2])*normal[2]); + // RealType offset = ((0.0-p[0][0])*normal[0] + (0.0-p[0][1])*normal[1] + (0.0-p[0][2])*normal[2]); RealType dist = normal[0] * interiorPoint[0] + normal[1]*interiorPoint[1] + normal[2]*interiorPoint[2]; - cout << "Dist and normal and area are: " << normal << endl; volume_ += dist *area/qh hull_dim; Triangles_.push_back(face); } } - cout << "Volume is: " << volume_ << endl; //assert(pm.cm.fn == ridgesCount); /*