Skip to content
Snippets Groups Projects
Commit 702a3b8f authored by Tristan Carrier Baudouin's avatar Tristan Carrier Baudouin
Browse files

centerline

parent 66784a07
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Mesh/CenterlineField.cpp
+ 80
− 0
View file @ 702a3b8f
...@@ -1102,6 +1102,86 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt ...@@ -1102,6 +1102,86 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
return; 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{ void Centerline::printSplit() const{
FILE * f = fopen("mySPLIT.pos","w"); FILE * f = fopen("mySPLIT.pos","w");
......
Mesh/CenterlineField.h
+ 7
− 1
View file @ 702a3b8f
...@@ -116,6 +116,9 @@ class Centerline : public Field{ ...@@ -116,6 +116,9 @@ class Centerline : public Field{
//anisotropic operator //anisotropic operator
void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0); 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) //import the 1D mesh of the centerlines (in vtk format)
//and fill the vector of lines //and fill the vector of lines
void importFile(std::string fileName); void importFile(std::string fileName);
...@@ -183,6 +186,9 @@ class Centerline : public Field{ ...@@ -183,6 +186,9 @@ class Centerline : public Field{
//anisotropic operator //anisotropic operator
void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0); 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 #endif
......
Mesh/Field.h
+ 4
− 0
View file @ 702a3b8f
...@@ -86,6 +86,10 @@ class Field { ...@@ -86,6 +86,10 @@ class Field {
virtual double operator() (double x, double y, double z, GEntity *ge=0) = 0; virtual double operator() (double x, double y, double z, GEntity *ge=0) = 0;
// anisotropic // anisotropic
virtual void operator() (double x, double y, double z, SMetric3 &, GEntity *ge=0){} 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; bool update_needed;
//void update(){ printf("up f \n"); return;} //void update(){ printf("up f \n"); return;}
virtual const char *getName() = 0; virtual const char *getName() = 0;
......
Mesh/simple3D.cpp
+ 67
− 19
View file @ 702a3b8f
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include "meshGRegion.h" #include "meshGRegion.h"
#include <queue> #include <queue>
#include <fstream> #include <fstream>
#include "CenterlineField.h"
#if defined(HAVE_RTREE) #if defined(HAVE_RTREE)
#include "rtree.h" #include "rtree.h"
...@@ -17,7 +18,7 @@ ...@@ -17,7 +18,7 @@
#define k1 0.7 //k1*h is the minimal distance between two nodes #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 k2 0.5 //k2*h is the minimal distance to the boundary
#define sqrt2 1.4143 #define sqrt3 1.73205081
/*********definitions*********/ /*********definitions*********/
...@@ -371,7 +372,7 @@ void Filler::treat_region(GRegion* gr){ ...@@ -371,7 +372,7 @@ void Filler::treat_region(GRegion* gr){
y = boundary_vertices[i]->y(); y = boundary_vertices[i]->y();
z = boundary_vertices[i]->z(); z = boundary_vertices[i]->z();
node = new Node(SPoint3(x,y,z)); node = new Node(SPoint3(x,y,z));
compute_parameters(node); compute_parameters(node,gr);
rtree.Insert(node->min,node->max,node); rtree.Insert(node->min,node->max,node);
fifo.push(node); fifo.push(node);
print_node(node,file); print_node(node,file);
...@@ -389,7 +390,7 @@ void Filler::treat_region(GRegion* gr){ ...@@ -389,7 +390,7 @@ void Filler::treat_region(GRegion* gr){
n4 = new Node(); n4 = new Node();
n5 = new Node(); n5 = new Node();
n6 = 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.clear();
data.push_back(n1); data.push_back(n1);
data.push_back(n2); data.push_back(n2);
...@@ -405,7 +406,7 @@ void Filler::treat_region(GRegion* gr){ ...@@ -405,7 +406,7 @@ void Filler::treat_region(GRegion* gr){
y = point.y(); y = point.y();
z = point.z(); z = point.z();
if(inside_domain(octree,x,y,z)){ if(inside_domain(octree,x,y,z)){
compute_parameters(spawn); compute_parameters(spawn,gr);
if(far_from_boundary(octree,spawn)){ if(far_from_boundary(octree,spawn)){
wrapper.set_too_close(0); wrapper.set_too_close(0);
wrapper.set_spawn(spawn); wrapper.set_spawn(spawn);
...@@ -466,6 +467,53 @@ double Filler::get_size(double x,double y,double z){ ...@@ -466,6 +467,53 @@ double Filler::get_size(double x,double y,double z){
return 0.25; 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){ bool Filler::inside_domain(MElementOctree* octree,double x,double y,double z){
MElement* element; MElement* element;
element = (MElement*)octree->find(x,y,z,3,true); element = (MElement*)octree->find(x,y,z,3,true);
...@@ -496,7 +544,7 @@ bool Filler::far_from_boundary(MElementOctree* octree,Node* node){ ...@@ -496,7 +544,7 @@ bool Filler::far_from_boundary(MElementOctree* octree,Node* node){
else return 0; else return 0;
} }
void Filler::compute_parameters(Node* node){ void Filler::compute_parameters(Node* node,GEntity* ge){
double x,y,z; double x,y,z;
double h; double h;
Metric m; Metric m;
...@@ -510,15 +558,15 @@ void Filler::compute_parameters(Node* node){ ...@@ -510,15 +558,15 @@ void Filler::compute_parameters(Node* node){
h = get_size(x,y,z); h = get_size(x,y,z);
node->set_size(h); node->set_size(h);
node->set_metric(m); node->set_metric(m);
node->min[0] = x - sqrt2*k1*h; node->min[0] = x - sqrt3*h;
node->min[1] = y - sqrt2*k1*h; node->min[1] = y - sqrt3*h;
node->min[2] = z - sqrt2*k1*h; node->min[2] = z - sqrt3*h;
node->max[0] = x + sqrt2*k1*h; node->max[0] = x + sqrt3*h;
node->max[1] = y + sqrt2*k1*h; node->max[1] = y + sqrt3*h;
node->max[2] = z + sqrt2*k1*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 x,y,z;
double x1,y1,z1; double x1,y1,z1;
double x2,y2,z2; double x2,y2,z2;
...@@ -538,32 +586,32 @@ void Filler::offsprings(MElementOctree* octree,Node* node,Node* n1,Node* n2,Node ...@@ -538,32 +586,32 @@ void Filler::offsprings(MElementOctree* octree,Node* node,Node* n1,Node* n2,Node
h = node->get_size(); h = node->get_size();
m = node->get_metric(); 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(); x1 = x + h1*m.get_m11();
y1 = y + h1*m.get_m21(); y1 = y + h1*m.get_m21();
z1 = z + h1*m.get_m31(); 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(); x2 = x - h2*m.get_m11();
y2 = y - h2*m.get_m21(); y2 = y - h2*m.get_m21();
z2 = z - h2*m.get_m31(); 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(); x3 = x + h3*m.get_m12();
y3 = y + h3*m.get_m22(); y3 = y + h3*m.get_m22();
z3 = z + h3*m.get_m32(); 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(); x4 = x - h4*m.get_m12();
y4 = y - h4*m.get_m22(); y4 = y - h4*m.get_m22();
z4 = z - h4*m.get_m32(); 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(); x5 = x + h5*m.get_m13();
y5 = y + h5*m.get_m23(); y5 = y + h5*m.get_m23();
z5 = z + h5*m.get_m33(); 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(); x6 = x - h6*m.get_m13();
y6 = y - h6*m.get_m23(); y6 = y - h6*m.get_m23();
z6 = z - h6*m.get_m33(); z6 = z - h6*m.get_m33();
...@@ -576,7 +624,7 @@ void Filler::offsprings(MElementOctree* octree,Node* node,Node* n1,Node* n2,Node ...@@ -576,7 +624,7 @@ void Filler::offsprings(MElementOctree* octree,Node* node,Node* n1,Node* n2,Node
*n6 = Node(SPoint3(x6,y6,z6)); *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 x,y,z;
double average; double average;
double h_farther; double h_farther;
......
Mesh/simple3D.h
+ 5
− 3
View file @ 702a3b8f
...@@ -16,11 +16,13 @@ class Filler{ ...@@ -16,11 +16,13 @@ class Filler{
static std::vector<MVertex*> new_vertices; static std::vector<MVertex*> new_vertices;
Metric get_metric(double,double,double); Metric get_metric(double,double,double);
double get_size(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 inside_domain(MElementOctree*,double,double,double);
bool far_from_boundary(MElementOctree*,Node*); bool far_from_boundary(MElementOctree*,Node*);
void compute_parameters(Node*); void compute_parameters(Node*,GEntity*);
void offsprings(MElementOctree*,Node*,Node*,Node*,Node*,Node*,Node*,Node*); void offsprings(GEntity*,MElementOctree*,Node*,Node*,Node*,Node*,Node*,Node*,Node*);
double improvement(MElementOctree*,SPoint3,double,SVector3); double improvement(GEntity*,MElementOctree*,SPoint3,double,SVector3);
void print_segment(SPoint3,SPoint3,std::ofstream&); void print_segment(SPoint3,SPoint3,std::ofstream&);
void print_node(Node*,std::ofstream&); void print_node(Node*,std::ofstream&);
public: public:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment