diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index f882ec62c6d4adabe9528909a83723d6a735832e..cbfb2b4c076ab11cea53924386e14b4ff03fac51 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -1102,6 +1102,86 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
return;
}
+/**************************************************/
+/******************Temporary code******************/
+/**************************************************/
+void Centerline::operator()(double x,double y,double z,SVector3& v1,SVector3& v2,SVector3& v3,GEntity* ge){
+ if(update_needed){
+ std::ifstream input;
+ input.open(fileName.c_str());
+ if(StatFile(fileName)) Msg::Fatal("Centerline file '%s' does not exist", fileName.c_str());
+ importFile(fileName);
+ buildKdTree();
+ update_needed = false;
+ }
+
+ double xyz[3] = {x,y,z};
+
+ //take xyz = closest point on boundary in case we are on the planar IN/OUT FACES or in VOLUME
+ bool onTubularSurface = true;
+ double ds = 0.0;
+ bool isCompound = (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) ? true: false;
+ std::list<GFace*> cFaces;
+ if(isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
+ if(ge->dim()==3 || (ge->dim()==2&&ge->geomType()==GEntity::Plane) || (isCompound&&(*cFaces.begin())->geomType()==GEntity::Plane)){
+ onTubularSurface = false;
+ kdtreeR->annkSearch(xyz, 1, index, dist);
+ 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];
+ kdtree->annkSearch(xyz, 2, index2, dist2);
+ double radMax = sqrt(dist2[0]);
+ 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]);
+
+ //dir_a = direction along the centerline
+ //dir_n = normal direction of the disk
+ //dir_t = tangential direction of the disk
+ SVector3 dir_a = p1-p0;
+ SVector3 dir_n(x-p0.x(), y-p0.y(), z-p0.z());
+ SVector3 dir_t;
+ buildOrthoBasis2(dir_a, dir_n, dir_t);
+
+ double lc = 2*M_PI*radMax/nbPoints;
+ double lc_a = 3.*lc;
+ double lc_n, lc_t;
+ if(onTubularSurface){
+ lc_n = lc_t = lc;
+ }
+ else{
+ double e = radMax/3.;
+ double hn = e/10.;
+ double rm = std::max(radMax-ds, radMax-e);
+ lc_t = 2*M_PI*rm/nbPoints;
+ //lc_n = lc_t = lc;
+ //double ratio = 1.02; //1. + (lc_t-hn)/e;
+ //printf("ratio =%g \n", ratio);
+ //lc_n = ds*(ratio-1) + hn;
+ if(ds<e) lc_n = hn;
+ else lc_n = lc_t;
+ }
+ double lam_a = 1./(lc_a*lc_a);
+ double lam_n = 1./(lc_n*lc_n);
+ double lam_t = 1./(lc_t*lc_t);
+
+ v1 = dir_a;
+ v2 = dir_n;
+ v3 = dir_t;
+
+ delete[]index2;
+ delete[]dist2;
+
+ return;
+}
+/*************************************************/
+/*************************************************/
+/*************************************************/
+
void Centerline::printSplit() const{
FILE * f = fopen("mySPLIT.pos","w");
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index efba5de307a36dd0d07529ebcb816adf121bbe7a..6e82d6361fbd709ab12d7a23de85d4760add4fa3 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -115,7 +115,10 @@ class Centerline : public Field{
double operator() (double x, double y, double z, GEntity *ge=0);
//anisotropic operator
void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
-
+
+ //temporary operator where v1, v2 and v3 are three orthonormal directions
+ void operator()(double x,double y,double z,SVector3& v1,SVector3& v2,SVector3& v3,GEntity* ge=0);
+
//import the 1D mesh of the centerlines (in vtk format)
//and fill the vector of lines
void importFile(std::string fileName);
@@ -183,6 +186,9 @@ class Centerline : public Field{
//anisotropic operator
void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
+ //temporary operator where v1, v2 and v3 are three orthonormal directions
+ void operator()(double x,double y,double z,SVector3& v1,SVector3& v2,SVector3& v3,GEntity* ge=0);
+
};
#endif
diff --git a/Mesh/Field.h b/Mesh/Field.h
index 472e6337419a3a23d4efafc8c65300db3c7e0dc9..8b562efc5c4aa7c34f60d54dd81274575d94f09d 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -86,6 +86,10 @@ class Field {
virtual double operator() (double x, double y, double z, GEntity *ge=0) = 0;
// anisotropic
virtual void operator() (double x, double y, double z, SMetric3 &, GEntity *ge=0){}
+
+ //temporary
+ virtual void operator()(double x,double y,double z,SVector3& v1,SVector3& v2,SVector3& v3,GEntity* ge=0){}
+
bool update_needed;
//void update(){ printf("up f \n"); return;}
virtual const char *getName() = 0;
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index 1ea5e074ef275c2fd67b0a5e35e70240ac61f694..0a444e8d13e4335c8ee695fee3f837502651c330 100755
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -10,6 +10,7 @@
#include "meshGRegion.h"
#include <queue>
#include <fstream>
+#include "CenterlineField.h"
#if defined(HAVE_RTREE)
#include "rtree.h"
@@ -17,7 +18,7 @@
#define k1 0.7 //k1*h is the minimal distance between two nodes
#define k2 0.5 //k2*h is the minimal distance to the boundary
-#define sqrt2 1.4143
+#define sqrt3 1.73205081
/*********definitions*********/
@@ -371,7 +372,7 @@ void Filler::treat_region(GRegion* gr){
y = boundary_vertices[i]->y();
z = boundary_vertices[i]->z();
node = new Node(SPoint3(x,y,z));
- compute_parameters(node);
+ compute_parameters(node,gr);
rtree.Insert(node->min,node->max,node);
fifo.push(node);
print_node(node,file);
@@ -389,7 +390,7 @@ void Filler::treat_region(GRegion* gr){
n4 = new Node();
n5 = new Node();
n6 = new Node();
- offsprings(octree,parent,n1,n2,n3,n4,n5,n6);
+ offsprings(gr,octree,parent,n1,n2,n3,n4,n5,n6);
data.clear();
data.push_back(n1);
data.push_back(n2);
@@ -405,7 +406,7 @@ void Filler::treat_region(GRegion* gr){
y = point.y();
z = point.z();
if(inside_domain(octree,x,y,z)){
- compute_parameters(spawn);
+ compute_parameters(spawn,gr);
if(far_from_boundary(octree,spawn)){
wrapper.set_too_close(0);
wrapper.set_spawn(spawn);
@@ -466,6 +467,53 @@ double Filler::get_size(double x,double y,double z){
return 0.25;
}
+Metric Filler::get_clf_metric(double x,double y,double z,GEntity* ge){
+ SVector3 v1,v2,v3;
+ Metric m;
+ Field* field;
+ FieldManager* manager;
+
+ m = Metric();
+ manager = ge->model()->getFields();
+ if(manager->getBackgroundField()>0){
+ field = manager->get(manager->getBackgroundField());
+ if(field){
+ (*field)(x,y,z,v1,v2,v3,ge);
+
+ m.set_m11(v1.x());
+ m.set_m21(v1.y());
+ m.set_m31(v1.z());
+
+ m.set_m12(v2.x());
+ m.set_m22(v2.y());
+ m.set_m32(v2.z());
+
+ m.set_m13(v3.x());
+ m.set_m23(v3.y());
+ m.set_m33(v3.z());
+ }
+ }
+
+ return m;
+}
+
+double Filler::get_clf_size(double x,double y,double z,GEntity* ge){
+ double h;
+ Field* field;
+ FieldManager* manager;
+
+ h = 0.25;
+ manager = ge->model()->getFields();
+ if(manager->getBackgroundField()>0){
+ field = manager->get(manager->getBackgroundField());
+ if(field){
+ h = (*field)(x,y,z,ge);
+ }
+ }
+
+ return h;
+}
+
bool Filler::inside_domain(MElementOctree* octree,double x,double y,double z){
MElement* element;
element = (MElement*)octree->find(x,y,z,3,true);
@@ -496,7 +544,7 @@ bool Filler::far_from_boundary(MElementOctree* octree,Node* node){
else return 0;
}
-void Filler::compute_parameters(Node* node){
+void Filler::compute_parameters(Node* node,GEntity* ge){
double x,y,z;
double h;
Metric m;
@@ -510,15 +558,15 @@ void Filler::compute_parameters(Node* node){
h = get_size(x,y,z);
node->set_size(h);
node->set_metric(m);
- node->min[0] = x - sqrt2*k1*h;
- node->min[1] = y - sqrt2*k1*h;
- node->min[2] = z - sqrt2*k1*h;
- node->max[0] = x + sqrt2*k1*h;
- node->max[1] = y + sqrt2*k1*h;
- node->max[2] = z + sqrt2*k1*h;
+ node->min[0] = x - sqrt3*h;
+ node->min[1] = y - sqrt3*h;
+ node->min[2] = z - sqrt3*h;
+ node->max[0] = x + sqrt3*h;
+ node->max[1] = y + sqrt3*h;
+ node->max[2] = z + sqrt3*h;
}
-void Filler::offsprings(MElementOctree* octree,Node* node,Node* n1,Node* n2,Node* n3,Node* n4,Node* n5,Node* n6){
+void Filler::offsprings(GEntity* ge,MElementOctree* octree,Node* node,Node* n1,Node* n2,Node* n3,Node* n4,Node* n5,Node* n6){
double x,y,z;
double x1,y1,z1;
double x2,y2,z2;
@@ -538,32 +586,32 @@ void Filler::offsprings(MElementOctree* octree,Node* node,Node* n1,Node* n2,Node
h = node->get_size();
m = node->get_metric();
- h1 = improvement(octree,point,h,SVector3(m.get_m11(),m.get_m21(),m.get_m31()));
+ h1 = improvement(ge,octree,point,h,SVector3(m.get_m11(),m.get_m21(),m.get_m31()));
x1 = x + h1*m.get_m11();
y1 = y + h1*m.get_m21();
z1 = z + h1*m.get_m31();
- h2 = improvement(octree,point,h,SVector3(-m.get_m11(),-m.get_m21(),-m.get_m31()));
+ h2 = improvement(ge,octree,point,h,SVector3(-m.get_m11(),-m.get_m21(),-m.get_m31()));
x2 = x - h2*m.get_m11();
y2 = y - h2*m.get_m21();
z2 = z - h2*m.get_m31();
- h3 = improvement(octree,point,h,SVector3(m.get_m12(),m.get_m22(),m.get_m32()));
+ h3 = improvement(ge,octree,point,h,SVector3(m.get_m12(),m.get_m22(),m.get_m32()));
x3 = x + h3*m.get_m12();
y3 = y + h3*m.get_m22();
z3 = z + h3*m.get_m32();
- h4 = improvement(octree,point,h,SVector3(-m.get_m12(),-m.get_m22(),-m.get_m32()));
+ h4 = improvement(ge,octree,point,h,SVector3(-m.get_m12(),-m.get_m22(),-m.get_m32()));
x4 = x - h4*m.get_m12();
y4 = y - h4*m.get_m22();
z4 = z - h4*m.get_m32();
- h5 = improvement(octree,point,h,SVector3(m.get_m13(),m.get_m23(),m.get_m33()));
+ h5 = improvement(ge,octree,point,h,SVector3(m.get_m13(),m.get_m23(),m.get_m33()));
x5 = x + h5*m.get_m13();
y5 = y + h5*m.get_m23();
z5 = z + h5*m.get_m33();
- h6 = improvement(octree,point,h,SVector3(-m.get_m13(),-m.get_m23(),-m.get_m33()));
+ h6 = improvement(ge,octree,point,h,SVector3(-m.get_m13(),-m.get_m23(),-m.get_m33()));
x6 = x - h6*m.get_m13();
y6 = y - h6*m.get_m23();
z6 = z - h6*m.get_m33();
@@ -576,7 +624,7 @@ void Filler::offsprings(MElementOctree* octree,Node* node,Node* n1,Node* n2,Node
*n6 = Node(SPoint3(x6,y6,z6));
}
-double Filler::improvement(MElementOctree* octree,SPoint3 point,double h_nearer,SVector3 direction){
+double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,double h_nearer,SVector3 direction){
double x,y,z;
double average;
double h_farther;
diff --git a/Mesh/simple3D.h b/Mesh/simple3D.h
index 41e29a7af387772924862091092cf28ac33491e2..8d0538aae438450812584533db415666f7891f92 100755
--- a/Mesh/simple3D.h
+++ b/Mesh/simple3D.h
@@ -16,11 +16,13 @@ class Filler{
static std::vector<MVertex*> new_vertices;
Metric get_metric(double,double,double);
double get_size(double,double,double);
+ Metric get_clf_metric(double,double,double,GEntity*);
+ double get_clf_size(double,double,double,GEntity*);
bool inside_domain(MElementOctree*,double,double,double);
bool far_from_boundary(MElementOctree*,Node*);
- void compute_parameters(Node*);
- void offsprings(MElementOctree*,Node*,Node*,Node*,Node*,Node*,Node*,Node*);
- double improvement(MElementOctree*,SPoint3,double,SVector3);
+ void compute_parameters(Node*,GEntity*);
+ void offsprings(GEntity*,MElementOctree*,Node*,Node*,Node*,Node*,Node*,Node*,Node*);
+ double improvement(GEntity*,MElementOctree*,SPoint3,double,SVector3);
void print_segment(SPoint3,SPoint3,std::ofstream&);
void print_node(Node*,std::ofstream&);
public: