diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp index 56c922f3a584ee9fd357bbeb657cb580a039aa57..97c9a858acfa6c4d6f7b586858b570099e5388c7 100644 --- a/Geo/GFaceCompound.cpp +++ b/Geo/GFaceCompound.cpp @@ -1217,10 +1217,6 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound, kdtree = NULL; uv_kdtree = NULL; - uv_nodes = NULL; - nodes = NULL; - index = new ANNidx[2]; - dist = new ANNdist[2]; } GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound, @@ -1270,10 +1266,6 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound, uv_kdtree = NULL; kdtree = NULL; - uv_nodes = NULL; - nodes = NULL; - index = new ANNidx[2]; - dist = new ANNdist[2]; } GFaceCompound::~GFaceCompound() @@ -1289,12 +1281,16 @@ GFaceCompound::~GFaceCompound() if (_lsys)delete _lsys; delete ONE; delete MONE; - if(uv_kdtree) delete uv_kdtree; - if(kdtree) delete kdtree; - if(uv_nodes) annDeallocPts(uv_nodes); - if(nodes) annDeallocPts(nodes); - delete[]index; - delete[]dist; + if(uv_kdtree) { + ANNpointArray uv_nodes = uv_kdtree->thePoints(); + if(uv_nodes) annDeallocPts(uv_nodes); + delete uv_kdtree; + } + if(kdtree){ + ANNpointArray nodes = kdtree->thePoints(); + if(nodes) annDeallocPts(nodes); + delete kdtree; + } } SPoint2 GFaceCompound::getCoordinates(MVertex *v) const @@ -1967,7 +1963,7 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const //build kdtree boundary nodes in parametric space printf("build bc kdtree \n"); int nbBCNodes = myBCNodes.size(); - uv_nodes = annAllocPts(nbBCNodes, 3); + ANNpointArray uv_nodes = annAllocPts(nbBCNodes, 3); std::set<SPoint2>::iterator itp = myBCNodes.begin(); int ind = 0; while (itp != myBCNodes.end()){ @@ -2004,7 +2000,10 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const printf("not found in new octree \n"); GPoint gp(50,50,50); double pt[3] = {par1,par2,0.0}; + ANNidx index[2]; + ANNdist dist[2]; uv_kdtree->annkSearch(pt, 2, index, dist); + ANNpointArray uv_nodes = uv_kdtree->thePoints(); SPoint3 p1(uv_nodes[index[0]][0], uv_nodes[index[0]][1], uv_nodes[index[0]][2]); SPoint3 p2(uv_nodes[index[1]][0], uv_nodes[index[1]][1], uv_nodes[index[1]][2]); SPoint3 pnew; double d; @@ -2064,6 +2063,8 @@ GPoint GFaceCompound::point(double par1, double par2) const else{ //printf("look in kdtree \n"); double pt[3] = {par1,par2, 0.0}; + ANNidx index[2]; + ANNdist dist[2]; kdtree->annkSearch(pt, 1, index, dist); lt = &(_gfct[index[0]]); @@ -2169,6 +2170,8 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const printf("FIRSTDER POINT NOT FOUND --> kdtree \n"); printf("uv=%g %g \n", param.x(), param.y()); double pt[3] = {param.x(), param.y(), 0.0}; + ANNidx index[2]; + ANNdist dist[2]; kdtree->annkSearch(pt, 1, index, dist); tri = (&_gfct[index[0]])->tri; } @@ -2399,7 +2402,7 @@ void GFaceCompound::buildOct() const //ANN octree std::set<MVertex*> allVS; - nodes = annAllocPts(count, 3); + ANNpointArray nodes = annAllocPts(count, 3); // make bounding box SPoint3 bbmin = bb.min(), bbmax = bb.max(); diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h index 0296fb48e4bee5663768a797feaecd5fd6784c2b..0022b889823c9ddd5d59729f343b72e7de51aa5d 100644 --- a/Geo/GFaceCompound.h +++ b/Geo/GFaceCompound.h @@ -21,7 +21,6 @@ template <class scalar> class simpleFunction; #include "linearSystem.h" #include "GRbf.h" #include "MElementOctree.h" -#include <ANN/ANN.h> class ANNkd_tree; #define AR_MAX 5 //maximal geometrical aspect ratio @@ -100,10 +99,6 @@ class GFaceCompound : public GFace { linearSystem <double> *_lsys; mutable ANNkd_tree *uv_kdtree; mutable ANNkd_tree *kdtree; - mutable ANNpointArray uv_nodes; - mutable ANNpointArray nodes; - ANNidxArray index; - ANNdistArray dist; void buildOct() const ; void buildAllNodes() const; diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp index 7d38dd5dfee086cdaa143eb71521dc1f0c33fcd0..fe161d2dc7e2c8728bd8c3d8fadc2fb19b4a6829 100644 --- a/Geo/gmshLevelset.cpp +++ b/Geo/gmshLevelset.cpp @@ -16,6 +16,10 @@ #include "MElement.h" #include "Numeric.h" #include "cartesian.h" +#include "GmshConfig.h" +#if defined(HAVE_ANN) +#include "ANN/ANN.h" +#endif void insertActiveCells(double x, double y, double z, double rmax, cartesianBox<double> &box) @@ -909,33 +913,35 @@ gLevelsetDistMesh::gLevelsetDistMesh(GModel *gm, std::string physical, int nbClo } } } - _nodes = annAllocPts(_all.size(), 3); + ANNpointArray nodes; + nodes = annAllocPts(_all.size(), 3); std::set<MVertex*>::iterator itp = _all.begin(); int ind = 0; while (itp != _all.end()){ MVertex* pt = *itp; - _nodes[ind][0] = pt->x(); - _nodes[ind][1] = pt->y(); - _nodes[ind][2] = pt->z(); + nodes[ind][0] = pt->x(); + nodes[ind][1] = pt->y(); + nodes[ind][2] = pt->z(); _vertices.push_back(pt); itp++; ind++; } - _kdtree = new ANNkd_tree(_nodes, _all.size(), 3); - _index = new ANNidx[_nbClose]; - _dist = new ANNdist[_nbClose]; + _kdtree = new ANNkd_tree(nodes, _all.size(), 3); } gLevelsetDistMesh::~gLevelsetDistMesh() { - delete [] _index; - delete [] _dist; - if (_kdtree) delete _kdtree; - if (_nodes) annDeallocPts (_nodes); + if (_kdtree) { + ANNpointArray nodes = _kdtree->thePoints(); + annDeallocPts(nodes); + delete _kdtree; + } } double gLevelsetDistMesh::operator () (double x, double y, double z) const { + ANNidx _index[_nbClose]; + ANNdist _dist[_nbClose]; double point[3] = {x,y,z}; _kdtree->annkSearch(point, _nbClose, _index, _dist); std::set<MElement*> elements; diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h index 4f11c0e061e947498eae358e89eb2a9140c6eca4..51efa0814e21f2b6e704fccc0f87433b8def8e62 100644 --- a/Geo/gmshLevelset.h +++ b/Geo/gmshLevelset.h @@ -23,7 +23,7 @@ #include "simpleFunction.h" #if defined(HAVE_ANN) -#include "ANN/ANN.h" +class ANNkd_tree; #endif #if defined(HAVE_POST) #include "PView.h" @@ -349,9 +349,6 @@ class gLevelsetDistMesh: public gLevelsetPrimitive std::vector<MVertex*> _vertices; std::multimap<MVertex*,MElement*> _v2e; ANNkd_tree *_kdtree; - ANNpointArray _nodes; - ANNidxArray _index; - ANNdistArray _dist; public : gLevelsetDistMesh(GModel *gm, std::string physical, int nbClose = 5, int tag=1); double operator () (double x, double y, double z) const; diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp index e0d6e32684fbf185237ef35532948d6ba7cde35e..fd4f9e5bc60b194286f78136712965feb650d01d 100644 --- a/Mesh/CenterlineField.cpp +++ b/Mesh/CenterlineField.cpp @@ -292,13 +292,10 @@ static void cutTriangle(MTriangle *tri, } } -Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0), nodes(0), nodesR(0) +Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0) { recombine = (CTX::instance()->mesh.recombineAll) || (CTX::instance()->mesh.recombine3DAll); - index = new ANNidx[1]; - dist = new ANNdist[1]; - printf("centerline filename =%s \n", fileName.c_str()); importFile(fileName); buildKdTree(); @@ -315,10 +312,8 @@ Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0), nodes(0), n } -Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0) +Centerline::Centerline(): kdtree(0), kdtreeR(0) { - index = new ANNidx[1]; - dist = new ANNdist[1]; recombine = (CTX::instance()->mesh.recombineAll) || (CTX::instance()->mesh.recombine3DAll); fileName = "centerlines.vtk";//default @@ -358,12 +353,16 @@ Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0) Centerline::~Centerline() { if (mod) delete mod; - if(kdtree) delete kdtree; - if(kdtreeR) delete kdtreeR; - if(nodes) annDeallocPts(nodes); - if(nodesR) annDeallocPts(nodesR); - delete[]index; - delete[]dist; + if(kdtree){ + ANNpointArray nodes = kdtree->thePoints(); + if(nodes) annDeallocPts(nodes); + delete kdtree; + } + if(kdtreeR){ + ANNpointArray nodesR = kdtreeR->thePoints(); + if(nodesR) annDeallocPts(nodesR); + delete kdtreeR; + } } void Centerline::importFile(std::string fileName) @@ -531,7 +530,7 @@ void Centerline::distanceToSurface() for(unsigned int j = 0; j < triangles.size(); j++) for(int k = 0; k<3; k++) allVS.insert(triangles[j]->getVertex(k)); int nbSNodes = allVS.size(); - nodesR = annAllocPts(nbSNodes, 3); + ANNpointArray nodesR = annAllocPts(nbSNodes, 3); vertices.resize(nbSNodes); int ind = 0; std::set<MVertex*>::iterator itp = allVS.begin(); @@ -550,6 +549,8 @@ void Centerline::distanceToSurface() MVertex *v1 = l->getVertex(0); MVertex *v2 = l->getVertex(1); double midp[3] = {0.5*(v1->x()+v2->x()), 0.5*(v1->y()+v1->y()),0.5*(v1->z()+v2->z())}; + ANNidx index[1]; + ANNdist dist[1]; kdtreeR->annkSearch(midp, 1, index, dist); double minRad = sqrt(dist[0]); radiusl.insert(std::make_pair(lines[i], minRad)); @@ -583,7 +584,7 @@ void Centerline::buildKdTree() //int nbNodes = (lines.size()+1) + (nbPL*lines.size()); int nbNodes = (colorp.size()) + (nbPL*lines.size()); - nodes = annAllocPts(nbNodes, 3); + ANNpointArray nodes = annAllocPts(nbNodes, 3); int ind = 0; std::map<MVertex*, int>::iterator itp = colorp.begin(); while (itp != colorp.end()){ @@ -771,7 +772,10 @@ void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundE SVector3 ne = e->getFace(0).normal(); SVector3 ps(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z()); double xyz[3] = {ps.x(), ps.y(), ps.z()}; + ANNidx index[1]; + ANNdist dist[1]; kdtree->annkSearch(xyz, 1, index, dist); + ANNpointArray nodes = kdtree->thePoints(); SVector3 pc(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]); SVector3 nc = ps-pc; if (dot(ne,nc) < 0) dir = 1; @@ -1077,7 +1081,10 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge) if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) || (isCompound && (*cFaces.begin())->geomType() == GEntity::Plane) ){ int num_neighbours = 1; + ANNidx index[num_neighbours]; + ANNdist dist[num_neighbours]; kdtreeR->annkSearch(xyz, num_neighbours, index, dist); + ANNpointArray nodesR = kdtreeR->thePoints(); xyz[0] = nodesR[index[0]][0]; xyz[1] = nodesR[index[0]][1]; xyz[2] = nodesR[index[0]][2]; @@ -1085,6 +1092,8 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge) } int num_neighbours = 1; + ANNidx index[num_neighbours]; + ANNdist dist[num_neighbours]; kdtree->annkSearch(xyz, num_neighbours, index, dist); double rad = sqrt(dist[0]); @@ -1130,24 +1139,25 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt onTubularSurface = false; } + ANNidx index[1]; + ANNdist dist[1]; kdtreeR->annkSearch(xyz, 1, index, dist); if (! onTubularSurface){ + ANNpointArray nodesR = kdtreeR->thePoints(); ds = sqrt(dist[0]); xyz[0] = nodesR[index[0]][0]; xyz[1] = nodesR[index[0]][1]; xyz[2] = nodesR[index[0]][2]; } - ANNidxArray index2 = new ANNidx[2]; - ANNdistArray dist2 = new ANNdist[2]; + ANNidx index2[2]; + ANNdist dist2[2]; kdtree->annkSearch(xyz, 2, index2, dist2); double radMax = sqrt(dist2[0]); + ANNpointArray nodes = kdtree->thePoints(); SVector3 p0(nodes[index2[0]][0], nodes[index2[0]][1], nodes[index2[0]][2]); SVector3 p1(nodes[index2[1]][0], nodes[index2[1]][1], nodes[index2[1]][2]); - delete[]index2; - delete[]dist2; - //dir_a = direction along the centerline //dir_n = normal direction of the disk //dir_t = tangential direction of the disk @@ -1185,7 +1195,7 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt double thickness = radMax/3.; double h_far = radMax/5.; double beta = (ds <= thickness) ? 1.2 : 2.1; //CTX::instance()->mesh.smoothRatio; - double dist = (ds <= thickness) ? ds: thickness; + double ddist = (ds <= thickness) ? ds: thickness; double h_n_0 = thickness/20.; double h_n = std::min( (h_n_0+ds*log(beta)), h_far); @@ -1198,10 +1208,10 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt (sqrt(1+ (4.*rhoMax*rhoMax*(betaMax*betaMax-1))/(h_n*h_n))-1.); double h_t1_0 = sqrt(1./oneOverD2_min); double h_t2_0 = sqrt(1./oneOverD2_max); - //double h_t1 = h_t1_0*(rhoMin+signMin*dist)/rhoMin ; - //double h_t2 = h_t2_0*(rhoMax+signMax*dist)/rhoMax ; - double h_t1 = std::min( (h_t1_0+(dist*log(beta))), radMax); - double h_t2 = std::min( (h_t2_0+(dist*log(beta))), h_far); + //double h_t1 = h_t1_0*(rhoMin+signMin*ddist)/rhoMin ; + //double h_t2 = h_t2_0*(rhoMax+signMax*ddist)/rhoMax ; + double h_t1 = std::min( (h_t1_0+(ddist*log(beta))), radMax); + double h_t2 = std::min( (h_t2_0+(ddist*log(beta))), h_far); double dCenter = radMax-ds; double h_a_0 = 0.5*radMax; diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h index abf5d57f9038ff58c3e2a61797a7612d76a6cbc7..525a1451c0495b69b893770844fa2498ce2e02a2 100644 --- a/Mesh/CenterlineField.h +++ b/Mesh/CenterlineField.h @@ -41,7 +41,6 @@ struct Branch{ }; #if defined(HAVE_ANN) -#include <ANN/ANN.h> class ANNkd_tree; // This class takes as input A 1D mesh which is the centerline @@ -58,9 +57,6 @@ class Centerline : public Field{ GModel *mod; //centerline GModel GModel *split; //split GModel ANNkd_tree *kdtree, *kdtreeR; - ANNpointArray nodes, nodesR; - ANNidxArray index; - ANNdistArray dist; std::string fileName; int nbPoints; double recombine; diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp index 9e6adac3882b1be986711f788f599b77f82ad3fc..7047e18772a75463dc363ae0e5b833fc29a98436 100644 --- a/Mesh/directions3D.cpp +++ b/Mesh/directions3D.cpp @@ -47,7 +47,7 @@ void Frame_field::init_region(GRegion* gr){ init_face(gf); } - duplicate = annAllocPts(temp.size(),3); + ANNpointArray duplicate = annAllocPts(temp.size(),3); index = 0; for(std::map<MVertex*,STensor3>::iterator it=temp.begin(); it != temp.end(); it++){ @@ -431,12 +431,12 @@ void Frame_field::clear(){ temp.clear(); field.clear(); #if defined(HAVE_ANN) - delete duplicate; + delete kd_tree->thePoints(); delete kd_tree; annClose(); #endif #if defined(HAVE_ANN) - if(annTreeData) delete annTreeData; + if(annTree && annTree->thePoints()) delete annTree->thePoints(); if(annTree) delete annTree; #endif } @@ -517,7 +517,7 @@ int Frame_field::buildAnnData(GEntity* ge, int dim){ build_listVertices(ge,dim); int n = listVertices.size(); #if defined(HAVE_ANN) - annTreeData = annAllocPts(n,3); + ANNpointArray annTreeData = annAllocPts(n,3); for(int i=0; i<n; i++){ MVertex* pVertex = listVertices[i]; annTreeData[i][0] = pVertex->x(); @@ -532,9 +532,8 @@ int Frame_field::buildAnnData(GEntity* ge, int dim){ void Frame_field::deleteAnnData(){ #if defined(HAVE_ANN) - if(annTreeData) delete annTreeData; + if(annTree && annTree->thePoints()) delete annTree->thePoints(); if(annTree) delete annTree; - annTreeData = NULL; annTree = NULL; #endif } @@ -1486,7 +1485,7 @@ void Nearest_point::init_region(GRegion* gr){ //vicinity.push_back(NULL); } - duplicate = annAllocPts(field.size(),3); + ANNpointArray duplicate = annAllocPts(field.size(),3); for(i=0;i<field.size();i++){ duplicate[i][0] = field[i].x(); @@ -1706,7 +1705,7 @@ void Nearest_point::clear(){ field.clear(); vicinity.clear(); #if defined(HAVE_ANN) - delete duplicate; + delete kd_tree->thePoints(); delete kd_tree; annClose(); #endif @@ -1722,9 +1721,7 @@ std::map<MVertex*,std::set<MVertex*> > Frame_field::vertex_to_vertices; std::map<MVertex*,std::set<MElement*> > Frame_field::vertex_to_elements; std::vector<MVertex*> Frame_field::listVertices; #if defined(HAVE_ANN) -ANNpointArray Frame_field::duplicate; ANNkd_tree* Frame_field::kd_tree; -ANNpointArray Frame_field::annTreeData; ANNkd_tree* Frame_field::annTree; #endif @@ -1734,6 +1731,5 @@ MElementOctree* Size_field::octree; std::vector<SPoint3> Nearest_point::field; std::vector<MElement*> Nearest_point::vicinity; #if defined(HAVE_ANN) -ANNpointArray Nearest_point::duplicate; ANNkd_tree* Nearest_point::kd_tree; #endif diff --git a/Mesh/directions3D.h b/Mesh/directions3D.h index a1f67511520df980b0fc3cec65bfd2016df3ab85..b52f7947b14dc49e6566e1513c2dd45357d9943e 100644 --- a/Mesh/directions3D.h +++ b/Mesh/directions3D.h @@ -13,7 +13,7 @@ #include "MEdge.h" #include "MElementOctree.h" #if defined(HAVE_ANN) -#include <ANN/ANN.h> +class ANNkd_tree; #endif #include "yamakawa.h" #include "STensor3.h" @@ -31,9 +31,7 @@ class Frame_field{ static std::map<MEdge, double, Less_Edge> crossDist; static std::vector<MVertex*> listVertices; #if defined(HAVE_ANN) - static ANNpointArray duplicate; static ANNkd_tree* kd_tree; - static ANNpointArray annTreeData; static ANNkd_tree* annTree; #endif Frame_field(); @@ -96,7 +94,6 @@ class Nearest_point{ static std::vector<SPoint3> field; static std::vector<MElement*> vicinity; #if defined(HAVE_ANN) - static ANNpointArray duplicate; static ANNkd_tree* kd_tree; #endif Nearest_point();