diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index a9631630a9419b27588a58ac802d4616c16e9215..6c19b7f7bd7fc4594672398d2c83bcb314b657fe 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -2059,10 +2059,10 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const {
//now find point
double uvw[3]={par1,par2, 0.0};
double UV[3];
- double initialTol = MElement::getTolerance();
- MElement::setTolerance(1.e-2);
- MElement *e = octNew->find(par1,par2, 0.0,-1,false);
- MElement::setTolerance(initialTol);
+ //double initialTol = MElement::getTolerance();
+ //MElement::setTolerance(1.e-2);
+ MElement *e = octNew->find(par1,par2, 0.0,-1,true);
+ //MElement::setTolerance(initialTol);
if (e){
e->xyz2uvw(uvw,UV);
double valX[8], valY[8], valZ[8];
@@ -2087,7 +2087,7 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const {
#if defined(HAVE_ANN)
double pt[3] = {par1,par2,0.0};
- uv_kdtree->annkSearch(pt, 1, index, dist);
+ uv_kdtree->annkSearch(pt, 2, index, dist);
SPoint3 p1(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
SPoint3 p2(nodes[index[1]][0], nodes[index[1]][1], nodes[index[1]][2]);
SPoint3 pnew; double d;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index e11a62f94b9c72095d6f4c9245154f49f6fff91f..d5516a3600f07fa6eb6e1e452d10d3ad1b18ec5e 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -460,19 +460,22 @@ backgroundMesh::backgroundMesh(GFace *_gf)
// those triangles are local to the backgroundMesh so that
// they do not depend on the actual mesh that can be deleted
+ std::set<SPoint2> myBCNodes;
for (unsigned int i = 0; i < _gf->triangles.size(); i++){
MTriangle *e = _gf->triangles[i];
MVertex *news[3];
for (int j=0;j<3;j++){
- std::map<MVertex*,MVertex*>::iterator it = _3Dto2D.find(e->getVertex(j));
+ MVertex *v = e->getVertex(j);
+ std::map<MVertex*,MVertex*>::iterator it = _3Dto2D.find(v);
MVertex *newv =0;
if (it == _3Dto2D.end()){
SPoint2 p;
- bool success = reparamMeshVertexOnFace(e->getVertex(j), _gf, p);
+ bool success = reparamMeshVertexOnFace(v, _gf, p);
newv = new MVertex (p.x(), p.y(), 0.0);
_vertices.push_back(newv);
- _3Dto2D[e->getVertex(j)] = newv;
- _2Dto3D[newv] = e->getVertex(j);
+ _3Dto2D[v] = newv;
+ _2Dto3D[newv] = v;
+ if(v->onWhat()->dim()<2) myBCNodes.insert(p);
}
else newv = it->second;
news[j] = newv;
@@ -480,6 +483,24 @@ backgroundMesh::backgroundMesh(GFace *_gf)
_triangles.push_back(new MTriangle(news[0],news[1],news[2]));
}
+#if defined(HAVE_ANN)
+ //printf("creating uv kdtree %d \n", myBCNodes.size());
+ index = new ANNidx[2];
+ dist = new ANNdist[2];
+ nodes = annAllocPts(myBCNodes.size(), 3);
+ std::set<SPoint2>::iterator itp = myBCNodes.begin();
+ int ind = 0;
+ while (itp != myBCNodes.end()){
+ SPoint2 pt = *itp;
+ //fprintf(of, "SP(%g,%g,%g){%g};\n", pt.x(), pt.y(), 0.0, 10000);
+ nodes[ind][0] = pt.x();
+ nodes[ind][1] = pt.y();
+ nodes[ind][2] = 0.0;
+ itp++; ind++;
+ }
+ uv_kdtree = new ANNkd_tree(nodes, myBCNodes.size(), 3);
+#endif
+
// build a search structure
_octree = new MElementOctree(_triangles);
@@ -507,6 +528,12 @@ backgroundMesh::~backgroundMesh()
for (unsigned int i = 0; i < _vertices.size(); i++) delete _vertices[i];
for (unsigned int i = 0; i < _triangles.size(); i++) delete _triangles[i];
delete _octree;
+#if defined(HAVE_ANN)
+ if(uv_kdtree) delete uv_kdtree;
+ if(nodes) annDeallocPts(nodes);
+ delete[]index;
+ delete[]dist;
+#endif
}
static void propagateValuesOnFace(GFace *_gf,
@@ -772,8 +799,22 @@ double backgroundMesh::operator() (double u, double v, double w) const
double uv2[3];
MElement *e = _octree->find(u, v, w, 2, true);
if (!e) {
- Msg::Error("BGM octree: cannot find %g %g %g", u, v, w);
- return -1000.0;//0.4;
+#if defined(HAVE_ANN)
+ //printf("BGM octree not found --> find in kdtree \n");
+ double pt[3] = {u,v,0.0};
+ uv_kdtree->annkSearch(pt, 2, index, dist);
+ SPoint3 p1(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
+ SPoint3 p2(nodes[index[1]][0], nodes[index[1]][1], nodes[index[1]][2]);
+ SPoint3 pnew; double d;
+ signedDistancePointLine(p1,p2, SPoint3(u,v,0.0), d, pnew);
+ double uvw[3]={pnew.x(),pnew.y(), 0.0};
+ double UV[3];
+ e = _octree->find(pnew.x(), pnew.y(), 0.0, 2, true);
+#endif
+ if(!e){
+ Msg::Error("BGM octree: cannot find UVW=%g %g %g", u, v, w);
+ return -1000.0;//0.4;
+ }
}
e->xyz2uvw(uv, uv2);
std::map<MVertex*,double>::const_iterator itv1 = _sizes.find(e->getVertex(0));
@@ -799,8 +840,22 @@ double backgroundMesh::getAngle(double u, double v, double w) const
double uv2[3];
MElement *e = _octree->find(u, v, w, 2, true);
if (!e) {
- Msg::Error("BGM octree : cannot find %g %g %g", u, v, w);
- return 0.0;
+#if defined(HAVE_ANN)
+ //printf("BGM octree not found --> find in kdtree \n");
+ double pt[3] = {u,v,0.0};
+ uv_kdtree->annkSearch(pt, 2, index, dist);
+ SPoint3 p1(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
+ SPoint3 p2(nodes[index[1]][0], nodes[index[1]][1], nodes[index[1]][2]);
+ SPoint3 pnew; double d;
+ signedDistancePointLine(p1,p2, SPoint3(u,v,0.0), d, pnew);
+ double uvw[3]={pnew.x(),pnew.y(), 0.0};
+ double UV[3];
+ e = _octree->find(pnew.x(), pnew.y(), 0.0, 2,true);
+#endif
+ if(!e){
+ Msg::Error("BGM octree angle: cannot find UVW=%g %g %g", u, v, w);
+ return -1000.0;
+ }
}
e->xyz2uvw(uv, uv2);
std::map<MVertex*,double>::const_iterator itv1 = _angles.find(e->getVertex(0));
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index d7d0ae7c73d595d5e93966e364ec61f1dece77d2..8d5030c6856ea392b668e927e8364d37748384c1 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -11,6 +11,11 @@
#include <list>
#include "simpleFunction.h"
+#if defined(HAVE_ANN)
+#include <ANN/ANN.h>
+class ANNkd_tree;
+#endif
+
class MElementOctree;
class GFace;
class GVertex;
@@ -47,6 +52,12 @@ class backgroundMesh : public simpleFunction<double>
static backgroundMesh * _current;
backgroundMesh(GFace *);
~backgroundMesh();
+#if defined(HAVE_ANN)
+ mutable ANNkd_tree *uv_kdtree;
+ mutable ANNpointArray nodes;
+ ANNidxArray index;
+ ANNdistArray dist;
+#endif
public:
static void set(GFace *);
static void unset();
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 3f7838de4ab7a478cd0c2ae6cc647c51cd442794..8b12aef059497b165272c6ff96f145f9d7aecf31 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1876,6 +1876,7 @@ void meshGFace::operator() (GFace *gf)
twoPassesMesh--;
if (backgroundMesh::current()){
backgroundMesh::unset();
+ //backgroundMesh::set(gf);
}
if (CTX::instance()->mesh.saveAll){
backgroundMesh::set(gf);