From b4bf46dc41a33307078eb66420e4343d0214c69d Mon Sep 17 00:00:00 2001
From: Paul-Emile Bernard <paul-emile.bernard@uclouvain.be>
Date: Sun, 29 Jun 2014 01:25:45 +0000
Subject: [PATCH]
---
Mesh/BackgroundMesh.cpp | 119 +++++++++++++++++++++-------------------
Mesh/BackgroundMesh.h | 10 ++++
2 files changed, 74 insertions(+), 55 deletions(-)
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 59c31b377c..70509ca5f0 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -504,7 +504,7 @@ backgroundMesh::backgroundMesh(GFace *_gf, bool cfd)
_vertices.push_back(newv);
_3Dto2D[v] = newv;
_2Dto3D[newv] = v;
- if(v->onWhat()->dim()<2) myBCNodes.insert(p);
+ if(v->onWhat()->dim()<2) myBCNodes.insert(p);
}
else newv = it->second;
news[j] = newv;
@@ -514,20 +514,20 @@ backgroundMesh::backgroundMesh(GFace *_gf, bool cfd)
#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);
+ 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
@@ -569,8 +569,8 @@ backgroundMesh::~backgroundMesh()
}
static void propagateValuesOnFace(GFace *_gf,
- std::map<MVertex*,double> &dirichlet,
- bool in_parametric_plane = false)
+ std::map<MVertex*,double> &dirichlet,
+ bool in_parametric_plane = false)
{
#if defined(HAVE_SOLVER)
linearSystem<double> *_lsys = 0;
@@ -656,19 +656,19 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
MVertex *v1 = (*it)->lines[i]->getVertex(0);
MVertex *v2 = (*it)->lines[i]->getVertex(1);
- if (v1 != v2){
- double d = sqrt((v1->x() - v2->x()) * (v1->x() - v2->x()) +
- (v1->y() - v2->y()) * (v1->y() - v2->y()) +
- (v1->z() - v2->z()) * (v1->z() -v2->z()));
- for (int k=0;k<2;k++){
- MVertex *v = (*it)->lines[i]->getVertex(k);
- std::map<MVertex*, double>::iterator itv = sizes.find(v);
- if (itv == sizes.end())
- sizes[v] = log(d);
- else
- itv->second = 0.5 * (itv->second + log(d));
- }
- }
+ if (v1 != v2){
+ double d = sqrt((v1->x() - v2->x()) * (v1->x() - v2->x()) +
+ (v1->y() - v2->y()) * (v1->y() - v2->y()) +
+ (v1->z() - v2->z()) * (v1->z() -v2->z()));
+ for (int k=0;k<2;k++){
+ MVertex *v = (*it)->lines[i]->getVertex(k);
+ std::map<MVertex*, double>::iterator itv = sizes.find(v);
+ if (itv == sizes.end())
+ sizes[v] = log(d);
+ else
+ itv->second = 0.5 * (itv->second + log(d));
+ }
+ }
}
}
}
@@ -715,14 +715,14 @@ void backgroundMesh::propagateCrossFieldByDistance(GFace *_gf)
v[1] = (*it)->lines[i]->getVertex(1);
SPoint2 p1,p2;
reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
- /* a correct way of computing angles */
- Pair<SVector3, SVector3> der = _gf->firstDer((p1+p2)*.5);
- SVector3 t1 = der.first();
- SVector3 t2 (v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
- t1.normalize();
- t2.normalize();
- double _angle = angle (t1,t2);
- // double angle = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
+ /* a correct way of computing angles */
+ Pair<SVector3, SVector3> der = _gf->firstDer((p1+p2)*.5);
+ SVector3 t1 = der.first();
+ SVector3 t2 (v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
+ t1.normalize();
+ t2.normalize();
+ double _angle = angle (t1,t2);
+ // double angle = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
crossField2d::normalizeAngle (_angle);
for (int i=0;i<2;i++){
std::map<MVertex*,double>::iterator itc = _cosines4.find(v[i]);
@@ -732,7 +732,7 @@ void backgroundMesh::propagateCrossFieldByDistance(GFace *_gf)
its->second = 0.5*(its->second + sin(4*_angle));
}
else {
- _param[v[i]] = (i==0) ? p1 : p2;
+ _param[v[i]] = (i==0) ? p1 : p2;
_cosines4[v[i]] = cos(4*_angle);
_sines4[v[i]] = sin(4*_angle);
}
@@ -790,15 +790,15 @@ void backgroundMesh::propagateCrossField(GFace *_gf)
v[1] = (*it)->lines[i]->getVertex(1);
SPoint2 p1,p2;
reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
- Pair<SVector3, SVector3> der = _gf->firstDer((p1+p2)*.5);
- SVector3 t1 = der.first();
- SVector3 t2 = der.second();
- SVector3 n = crossprod(t1,t2);
- n.normalize();
- SVector3 d1(v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
- t1.normalize();
- d1.normalize();
- double _angle = myAngle (t1,d1,n);
+ Pair<SVector3, SVector3> der = _gf->firstDer((p1+p2)*.5);
+ SVector3 t1 = der.first();
+ SVector3 t2 = der.second();
+ SVector3 n = crossprod(t1,t2);
+ n.normalize();
+ SVector3 d1(v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
+ t1.normalize();
+ d1.normalize();
+ double _angle = myAngle (t1,d1,n);
crossField2d::normalizeAngle (_angle);
for (int i=0;i<2;i++){
std::map<MVertex*,double>::iterator itc = _cosines4.find(v[i]);
@@ -980,7 +980,7 @@ double backgroundMesh::getAngle(double u, double v, double w) const
}
void backgroundMesh::print(const std::string &filename, GFace *gf,
- const std::map<MVertex*,double> &_whatToPrint) const
+ const std::map<MVertex*,double> &_whatToPrint) const
{
FILE *f = Fopen (filename.c_str(),"w");
fprintf(f,"View \"Background Mesh\"{\n");
@@ -993,9 +993,9 @@ void backgroundMesh::print(const std::string &filename, GFace *gf,
std::map<MVertex*,double>::const_iterator itv3 = _whatToPrint.find(v3);
if (!gf){
fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
- v1->x(),v1->y(),v1->z(),
- v2->x(),v2->y(),v2->z(),
- v3->x(),v3->y(),v3->z(),itv1->second,itv2->second,itv3->second);
+ v1->x(),v1->y(),v1->z(),
+ v2->x(),v2->y(),v2->z(),
+ v3->x(),v3->y(),v3->z(),itv1->second,itv2->second,itv3->second);
}
else {
@@ -1003,9 +1003,9 @@ void backgroundMesh::print(const std::string &filename, GFace *gf,
GPoint p2 = gf->point(SPoint2(v2->x(),v2->y()));
GPoint p3 = gf->point(SPoint2(v3->x(),v3->y()));
fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
- p1.x(),p1.y(),p1.z(),
- p2.x(),p2.y(),p2.z(),
- p3.x(),p3.y(),p3.z(),itv1->second,itv2->second,itv3->second);
+ p1.x(),p1.y(),p1.z(),
+ p2.x(),p2.y(),p2.z(),
+ p3.x(),p3.y(),p3.z(),itv1->second,itv2->second,itv3->second);
}
}
fprintf(f,"};\n");
@@ -1016,4 +1016,13 @@ MElementOctree* backgroundMesh::get_octree(){
return _octree;
}
+MElement *backgroundMesh::getMeshElementByCoord(double u, double v, double w, bool strict)
+{
+ if(!_octree){
+ Msg::Debug("Rebuilding BackgroundMesh element octree");
+ _octree = new MElementOctree(_triangles);
+ }
+ return _octree->find(u,v,w, 2, strict);
+}
+
backgroundMesh* backgroundMesh::_current = 0;
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index a46f5b1ebc..69526e37ca 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -85,6 +85,16 @@ class backgroundMesh : public simpleFunction<double>
}
}
MElementOctree* get_octree();
+ MElement *getMeshElementByCoord(double u, double v, double w, bool strict=true);
+ int getNumMeshElements()const{return _triangles.size();}
+ std::vector<MVertex*>::iterator begin_vertices(){return _vertices.begin();}
+ std::vector<MVertex*>::iterator end_vertices(){return _vertices.end();}
+ std::vector<MVertex*>::const_iterator begin_vertices()const{return _vertices.begin();}
+ std::vector<MVertex*>::const_iterator end_vertices()const{return _vertices.end();}
+ std::vector<MElement*>::iterator begin_triangles(){return _triangles.begin();}
+ std::vector<MElement*>::iterator end_triangles(){return _triangles.end();}
+ std::vector<MElement*>::const_iterator begin_triangles()const{return _triangles.begin();}
+ std::vector<MElement*>::const_iterator end_triangles()const{return _triangles.end();}
};
SMetric3 buildMetricTangentToCurve (SVector3 &t, double l_t, double l_n);
--
GitLab