diff --git a/Mesh/Levy3D.cpp b/Mesh/Levy3D.cpp
index 3d1fc515d5e608aaa2ed4ee371139b75e03763b3..480208c9c49ae8aa9f1f9728ad52cbf24f7c3315 100644
--- a/Mesh/Levy3D.cpp
+++ b/Mesh/Levy3D.cpp
@@ -1463,7 +1463,7 @@ void LpSmoother::improve_model(){
void LpSmoother::improve_region(GRegion* gr)
{
-#if defined(HAVE_BFGS)
+ #if defined(HAVE_BFGS)
unsigned int i;
int offset;
double epsg;
@@ -1490,7 +1490,7 @@ void LpSmoother::improve_region(GRegion* gr)
alglib::real_1d_array x;
alglib::real_1d_array alglib_scales;
- Frame_field::init_model();
+ Frame_field::init_region(gr);
octree = new MElementOctree(gr->model());
for(i=0;i<gr->getNumMeshElements();i++){
@@ -1621,7 +1621,7 @@ void LpSmoother::improve_region(GRegion* gr)
free(initial_conditions);
free(scales);
Frame_field::clear();
-#endif
+ #endif
}
int LpSmoother::get_nbr_interior_vertices(){
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index 3619d2907442127199537d2566f91b1227ff5328..ca7bcab8f1aa8d8a1c3020b29773dfe20941f78b 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -112,25 +112,27 @@ double Matrix::get_m33(){
Frame_field::Frame_field(){}
-void Frame_field::init_model(){
+void Frame_field::init_region(GRegion* gr){
#if defined(HAVE_ANN)
int index;
MVertex* vertex;
GFace* gf;
- GModel* model = GModel::current();
- GModel::fiter it;
+ std::list<GFace*> faces;
+ std::list<GFace*>::iterator it;
std::map<MVertex*,Matrix>::iterator it2;
Matrix m;
+ Nearest_point::init_region(gr);
+
+ faces = gr->faces();
+
field.clear();
random.clear();
- for(it=model->firstFace();it!=model->lastFace();it++)
+ for(it=faces.begin();it!=faces.end();it++)
{
gf = *it;
- if(gf->geomType()==GEntity::CompoundSurface){
- init_face(gf);
- }
+ init_face(gf);
}
duplicate = annAllocPts(field.size(),3);
@@ -260,7 +262,7 @@ bool Frame_field::translate(GFace* gf,MElementOctree* octree,MVertex* vertex,SPo
}
Matrix Frame_field::search(double x,double y,double z){
- int val;
+ int index;
double e;
#if defined(HAVE_ANN)
ANNpoint query;
@@ -277,14 +279,63 @@ Matrix Frame_field::search(double x,double y,double z){
e = 0.0;
kd_tree->annkSearch(query,1,indices,distances,e);
- val = indices[0];
+ index = indices[0];
annDeallocPt(query);
delete[] indices;
delete[] distances;
#endif
- return random[val].second;
+ return random[index].second;
+}
+
+Matrix Frame_field::combine(double x,double y,double z){
+ bool ok;
+ double val1,val2,val3;
+ SVector3 vec,other;
+ SVector3 vec1,vec2,vec3;
+ SVector3 final1,final2;
+ Matrix m,m2;
+
+ m = search(x,y,z);
+ m2 = m;
+ ok = Nearest_point::search(x,y,z,vec);
+
+ if(ok){
+ vec1 = SVector3(m.get_m11(),m.get_m21(),m.get_m31());
+ vec2 = SVector3(m.get_m12(),m.get_m22(),m.get_m32());
+ vec3 = SVector3(m.get_m13(),m.get_m23(),m.get_m33());
+
+ val1 = fabs(dot(vec,vec1));
+ val2 = fabs(dot(vec,vec2));
+ val3 = fabs(dot(vec,vec3));
+
+ if(val1<=val2 && val1<=val3){
+ other = vec1;
+ }
+ else if(val2<=val1 && val2<=val3){
+ other = vec2;
+ }
+ else{
+ other = vec3;
+ }
+
+ final1 = crossprod(vec,other);
+ final1.normalize();
+ final2 = crossprod(vec,final1);
+
+ m2.set_m11(vec.x());
+ m2.set_m21(vec.y());
+ m2.set_m31(vec.z());
+ m2.set_m12(final1.x());
+ m2.set_m22(final1.y());
+ m2.set_m32(final1.z());
+ m2.set_m13(final2.x());
+ m2.set_m23(final2.y());
+ m2.set_m33(final2.z());
+ }
+
+ return m2;
}
bool Frame_field::inside_domain(MElementOctree* octree,double x,double y){
@@ -314,8 +365,8 @@ void Frame_field::print_field1(){
std::map<MVertex*,Matrix>::iterator it;
Matrix m;
- k = 0.1;
- std::ofstream file("field1.pos");
+ k = 0.05;
+ std::ofstream file("frame1.pos");
file << "View \"test\" {\n";
for(it=field.begin();it!=field.end();it++){
@@ -341,7 +392,7 @@ void Frame_field::print_field1(){
file << "};\n";
}
-void Frame_field::print_field2(){
+void Frame_field::print_field2(GRegion* gr){
unsigned int i;
int j;
double k;
@@ -349,40 +400,33 @@ void Frame_field::print_field2(){
SPoint3 p1,p2,p3,p4,p5,p6;
MVertex* vertex;
MElement* element;
- GRegion* gr;
- GModel* model = GModel::current();
- GModel::riter it;
Matrix m;
k = 0.05;
- std::ofstream file("field2.pos");
+ std::ofstream file("frame2.pos");
file << "View \"test\" {\n";
- for(it=model->firstRegion();it!=model->lastRegion();it++)
- {
- gr = *it;
- for(i=0;i<gr->getNumMeshElements();i++){
- element = gr->getMeshElement(i);
- for(j=0;j<element->getNumVertices();j++){
- vertex = element->getVertex(j);
- if(vertex->onWhat()->dim()>2){
- m = search(vertex->x(),vertex->y(),vertex->z());
-
- p = SPoint3(vertex->x(),vertex->y(),vertex->z());
- p1 = SPoint3(vertex->x()+k*m.get_m11(),vertex->y()+k*m.get_m21(),vertex->z()+k*m.get_m31());
- p2 = SPoint3(vertex->x()-k*m.get_m11(),vertex->y()-k*m.get_m21(),vertex->z()-k*m.get_m31());
- p3 = SPoint3(vertex->x()+k*m.get_m12(),vertex->y()+k*m.get_m22(),vertex->z()+k*m.get_m32());
- p4 = SPoint3(vertex->x()-k*m.get_m12(),vertex->y()-k*m.get_m22(),vertex->z()-k*m.get_m32());
- p5 = SPoint3(vertex->x()+k*m.get_m13(),vertex->y()+k*m.get_m23(),vertex->z()+k*m.get_m33());
- p6 = SPoint3(vertex->x()-k*m.get_m13(),vertex->y()-k*m.get_m23(),vertex->z()-k*m.get_m33());
-
- print_segment(p,p1,file);
- print_segment(p,p2,file);
- print_segment(p,p3,file);
- print_segment(p,p4,file);
- print_segment(p,p5,file);
- print_segment(p,p6,file);
- }
+ for(i=0;i<gr->getNumMeshElements();i++){
+ element = gr->getMeshElement(i);
+ for(j=0;j<element->getNumVertices();j++){
+ vertex = element->getVertex(j);
+ if(vertex->onWhat()->dim()>2){
+ m = combine(vertex->x(),vertex->y(),vertex->z());
+
+ p = SPoint3(vertex->x(),vertex->y(),vertex->z());
+ p1 = SPoint3(vertex->x()+k*m.get_m11(),vertex->y()+k*m.get_m21(),vertex->z()+k*m.get_m31());
+ p2 = SPoint3(vertex->x()-k*m.get_m11(),vertex->y()-k*m.get_m21(),vertex->z()-k*m.get_m31());
+ p3 = SPoint3(vertex->x()+k*m.get_m12(),vertex->y()+k*m.get_m22(),vertex->z()+k*m.get_m32());
+ p4 = SPoint3(vertex->x()-k*m.get_m12(),vertex->y()-k*m.get_m22(),vertex->z()-k*m.get_m32());
+ p5 = SPoint3(vertex->x()+k*m.get_m13(),vertex->y()+k*m.get_m23(),vertex->z()+k*m.get_m33());
+ p6 = SPoint3(vertex->x()-k*m.get_m13(),vertex->y()-k*m.get_m23(),vertex->z()-k*m.get_m33());
+
+ print_segment(p,p1,file);
+ print_segment(p,p2,file);
+ print_segment(p,p3,file);
+ print_segment(p,p4,file);
+ print_segment(p,p5,file);
+ print_segment(p,p6,file);
}
}
}
@@ -397,7 +441,19 @@ void Frame_field::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){
<< "{10, 20};\n";
}
+GRegion* Frame_field::test(){
+ GRegion* gr;
+ GModel* model = GModel::current();
+
+ gr = *(model->firstRegion());
+
+ return gr;
+}
+
void Frame_field::clear(){
+ #if defined(HAVE_ANN)
+ Nearest_point::clear();
+ #endif
field.clear();
random.clear();
#if defined(HAVE_ANN)
@@ -411,7 +467,7 @@ void Frame_field::clear(){
Size_field::Size_field(){}
-void Size_field::init_model(){
+void Size_field::init_region(GRegion* gr){
size_t i;
int j;
double local_size;
@@ -421,36 +477,38 @@ void Size_field::init_model(){
MVertex* vertex;
GFace* gf;
GModel* model = GModel::current();
- GModel::fiter it;
+ std::list<GFace*> faces;
+ std::list<GFace*>::iterator it;
+ faces = gr->faces();
+
boundary.clear();
-
- for(it=model->firstFace();it!=model->lastFace();it++){
+
+ for(it=faces.begin();it!=faces.end();it++){
gf = *it;
- if(gf->geomType()==GEntity::CompoundSurface){
- backgroundMesh::set(gf);
- for(i=0;i<gf->getNumMeshElements();i++){
- element = gf->getMeshElement(i);
+ backgroundMesh::set(gf);
+
+ for(i=0;i<gf->getNumMeshElements();i++){
+ element = gf->getMeshElement(i);
- average_x = 0.0;
- average_y = 0.0;
+ average_x = 0.0;
+ average_y = 0.0;
- for(j=0;j<element->getNumVertices();j++){
- vertex = element->getVertex(j);
- reparamMeshVertexOnFace(vertex,gf,point);
- average_x = average_x + point.x();
- average_y = average_y + point.y();
- }
+ for(j=0;j<element->getNumVertices();j++){
+ vertex = element->getVertex(j);
+ reparamMeshVertexOnFace(vertex,gf,point);
+ average_x = average_x + point.x();
+ average_y = average_y + point.y();
+ }
- average_x = average_x/element->getNumVertices();
- average_y = average_y/element->getNumVertices();
+ average_x = average_x/element->getNumVertices();
+ average_y = average_y/element->getNumVertices();
- for(j=0;j<element->getNumVertices();j++){
- vertex = element->getVertex(j);
- reparamMeshVertexOnFace(vertex,gf,point);
- local_size = backgroundMesh::current()->operator()(point.x(),point.y(),0.0);
- boundary.insert(std::pair<MVertex*,double>(vertex,/*local_size*/0.2));
- }
+ for(j=0;j<element->getNumVertices();j++){
+ vertex = element->getVertex(j);
+ reparamMeshVertexOnFace(vertex,gf,point);
+ local_size = backgroundMesh::current()->operator()(point.x(),point.y(),0.0);
+ boundary.insert(std::pair<MVertex*,double>(vertex,/*local_size*/0.2));
}
}
}
@@ -458,25 +516,18 @@ void Size_field::init_model(){
octree = new MElementOctree(model);
}
-void Size_field::solve(){
-#if defined(HAVE_PETSC)
+void Size_field::solve(GRegion* gr){
+ #if defined(HAVE_PETSC)
linearSystem<double>* system = 0;
- //system = new linearSystemPETSc<double>;
- system = new linearSystemFull<double>;
+ system = new linearSystemPETSc<double>;
+ //system = new linearSystemFull<double>;
size_t i;
- int j;
double val;
- MElement* element;
- MVertex* vertex;
- MTetrahedron* temp;
- GRegion* gr;
- GModel* model = GModel::current();
- GModel::riter it;
- std::map<MVertex*,double>::iterator it2;
- std::set<MVertex*>::iterator it3;
std::set<MVertex*> interior;
-
+ std::set<MVertex*>::iterator it;
+ std::map<MVertex*,double>::iterator it2;
+
interior.clear();
dofManager<double> assembler(system);
@@ -485,45 +536,36 @@ void Size_field::solve(){
assembler.fixVertex(it2->first,0,1,it2->second);
}
- for(it=model->firstRegion();it!=model->lastRegion();it++){
- gr = *it;
- for(i=0;i<gr->getNumMeshElements();i++){
- element = gr->getMeshElement(i);
- for(j=0;j<element->getNumVertices();j++){
- vertex = element->getVertex(j);
- interior.insert(vertex);
- }
- }
+ for(i=0;i<gr->tetrahedra.size();i++){
+ interior.insert(gr->tetrahedra[i]->getVertex(0));
+ interior.insert(gr->tetrahedra[i]->getVertex(1));
+ interior.insert(gr->tetrahedra[i]->getVertex(2));
+ interior.insert(gr->tetrahedra[i]->getVertex(3));
}
-
- for(it3=interior.begin();it3!=interior.end();it3++){
- it2 = boundary.find(*it3);
+
+ for(it=interior.begin();it!=interior.end();it++){
+ it2 = boundary.find(*it);
if(it2==boundary.end()){
- assembler.numberVertex(*it3,0,1);
+ assembler.numberVertex(*it,0,1);
}
}
simpleFunction<double> ONE(1.0);
laplaceTerm term(0,1,&ONE);
- for(it=model->firstRegion();it!=model->lastRegion();it++){
- gr = *it;
- for(i=0;i<gr->getNumMeshElements();i++){
- element = gr->getMeshElement(i);
- temp = (MTetrahedron*)element;
- SElement se(temp);
- term.addToMatrix(assembler,&se);
- }
+ for(i=0;i<gr->tetrahedra.size();i++){
+ SElement se(gr->tetrahedra[i]);
+ term.addToMatrix(assembler,&se);
}
system->systemSolve();
- for(it3=interior.begin();it3!=interior.end();it3++){
- assembler.getDofValue(*it3,0,1,val);
- boundary.insert(std::pair<MVertex*,double>(*it3,val));
+ for(it=interior.begin();it!=interior.end();it++){
+ assembler.getDofValue(*it,0,1,val);
+ boundary.insert(std::pair<MVertex*,double>(*it,val));
}
delete system;
-#endif
+ #endif
}
double Size_field::search(double x,double y,double z){
@@ -591,11 +633,208 @@ void Size_field::print_field(){
printf("%zu\n",boundary.size());
}
+GRegion* Size_field::test(){
+ GRegion* gr;
+ GModel* model = GModel::current();
+
+ gr = *(model->firstRegion());
+
+ return gr;
+}
+
void Size_field::clear(){
delete octree;
boundary.clear();
}
+/****************class Nearest_point****************/
+
+Nearest_point::Nearest_point(){}
+
+void Nearest_point::init_region(GRegion* gr){
+ #if defined(HAVE_ANN)
+ unsigned int i;
+ int j;
+ int gauss_num;
+ double u,v;
+ double x,y,z;
+ double x1,y1,z1;
+ double x2,y2,z2;
+ double x3,y3,z3;
+ MElement* element;
+ GFace* gf;
+ std::list<GFace*> faces;
+ std::set<MVertex*> vertices;
+ std::list<GFace*>::iterator it;
+ std::set<MVertex*>::iterator it2;
+ fullMatrix<double> gauss_points;
+ fullVector<double> gauss_weights;
+
+ gaussIntegration::getTriangle(12,gauss_points,gauss_weights);
+ gauss_num = gauss_points.size1();
+
+ faces = gr->faces();
+
+ random.clear();
+ vertices.clear();
+
+ for(it=faces.begin();it!=faces.end();it++){
+ gf = *it;
+ for(i=0;i<gf->getNumMeshElements();i++){
+ element = gf->getMeshElement(i);
+
+ x1 = element->getVertex(0)->x();
+ y1 = element->getVertex(0)->y();
+ z1 = element->getVertex(0)->z();
+
+ x2 = element->getVertex(1)->x();
+ y2 = element->getVertex(1)->y();
+ z2 = element->getVertex(1)->z();
+
+ x3 = element->getVertex(2)->x();
+ y3 = element->getVertex(2)->y();
+ z3 = element->getVertex(2)->z();
+
+ for(j=0;j<gauss_num;j++){
+ u = gauss_points(j,0);
+ v = gauss_points(j,1);
+
+ x = T(u,v,x1,x2,x3);
+ y = T(u,v,y1,y2,y3);
+ z = T(u,v,z1,z2,z3);
+
+ random.push_back(SPoint3(x,y,z));
+ }
+
+ vertices.insert(element->getVertex(0));
+ vertices.insert(element->getVertex(1));
+ vertices.insert(element->getVertex(2));
+ }
+ }
+
+ for(it2=vertices.begin();it2!=vertices.end();it2++){
+ x = (*it2)->x();
+ y = (*it2)->y();
+ z = (*it2)->z();
+ random.push_back(SPoint3(x,y,z));
+ }
+
+ duplicate = annAllocPts(random.size(),3);
+
+ for(i=0;i<random.size();i++){
+ duplicate[i][0] = random[i].x();
+ duplicate[i][1] = random[i].y();
+ duplicate[i][2] = random[i].z();
+ }
+
+ kd_tree = new ANNkd_tree(duplicate,random.size(),3);
+ #endif
+}
+
+bool Nearest_point::search(double x,double y,double z,SVector3& vec){
+ int index;
+ bool val;
+ #if defined(HAVE_ANN)
+ double e;
+ double e2;
+ SPoint3 point;
+ ANNpoint query;
+ ANNidxArray indices;
+ ANNdistArray distances;
+
+ query = annAllocPt(3);
+ query[0] = x;
+ query[1] = y;
+ query[2] = z;
+
+ indices = new ANNidx[1];
+ distances = new ANNdist[1];
+
+ e = 0.0;
+ kd_tree->annkSearch(query,1,indices,distances,e);
+ index = indices[0];
+
+ annDeallocPt(query);
+ delete[] indices;
+ delete[] distances;
+
+ e2 = 0.000001;
+ point = random[index];
+
+ if(fabs(point.x()-x)>e2 || fabs(point.y()-y)>e2 || fabs(point.z()-z)>e2){
+ vec = SVector3(point.x()-x,point.y()-y,point.z()-z);
+ vec.normalize();
+ val = 1;
+ }
+ else{
+ vec = SVector3(1.0,0.0,0.0);
+ val = 0;
+ }
+ #endif
+
+ return val;
+}
+
+double Nearest_point::T(double u,double v,double val1,double val2,double val3){
+ return (1.0-u-v)*val1 + u*val2 + v*val3;
+}
+
+void Nearest_point::print_field(GRegion* gr){
+ unsigned int i;
+ int j;
+ bool val;
+ double k;
+ double x,y,z;
+ MElement* element;
+ MVertex* vertex;
+ SVector3 vec;
+
+ k = 0.05;
+ std::ofstream file("nearest.pos");
+ file << "View \"test\" {\n";
+
+ for(i=0;i<gr->getNumMeshElements();i++){
+ element = gr->getMeshElement(i);
+ for(j=0;j<element->getNumVertices();j++){
+ vertex = element->getVertex(j);
+ x = vertex->x();
+ y = vertex->y();
+ z = vertex->z();
+ val = search(x,y,z,vec);
+ if(val){
+ print_segment(SPoint3(x+k*vec.x(),y+k*vec.y(),z+k*vec.z()),SPoint3(x-k*vec.x(),y-k*vec.y(),z-k*vec.z()),file);
+ }
+ }
+ }
+
+ file << "};\n";
+}
+
+void Nearest_point::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){
+ file << "SL ("
+ << p1.x() << ", " << p1.y() << ", " << p1.z() << ", "
+ << p2.x() << ", " << p2.y() << ", " << p2.z() << ")"
+ << "{10, 20};\n";
+}
+
+GRegion* Nearest_point::test(){
+ GRegion* gr;
+ GModel* model = GModel::current();
+
+ gr = *(model->firstRegion());
+
+ return gr;
+}
+
+void Nearest_point::clear(){
+ random.clear();
+ #if defined(HAVE_ANN)
+ delete duplicate;
+ delete kd_tree;
+ annClose();
+ #endif
+}
+
/****************static declarations****************/
std::map<MVertex*,Matrix> Frame_field::field;
@@ -604,5 +843,12 @@ std::vector<std::pair<MVertex*,Matrix> > Frame_field::random;
ANNpointArray Frame_field::duplicate;
ANNkd_tree* Frame_field::kd_tree;
#endif
+
std::map<MVertex*,double> Size_field::boundary;
-MElementOctree* Size_field::octree;
\ No newline at end of file
+MElementOctree* Size_field::octree;
+
+std::vector<SPoint3> Nearest_point::random;
+#if defined(HAVE_ANN)
+ANNpointArray Nearest_point::duplicate;
+ANNkd_tree* Nearest_point::kd_tree;
+#endif
\ No newline at end of file
diff --git a/Mesh/directions3D.h b/Mesh/directions3D.h
index 59c4a65e3db2746f37dabdd1eb069cffcaeb12b7..2b583e166b8687151e487089c76eb143b8b24908 100644
--- a/Mesh/directions3D.h
+++ b/Mesh/directions3D.h
@@ -48,15 +48,17 @@ class Frame_field{
#endif
Frame_field();
public:
- static void init_model();
+ static void init_region(GRegion*);
static void init_face(GFace*);
static bool translate(GFace*,MElementOctree*,MVertex*,SPoint2,SVector3&,SVector3&);
static Matrix search(double,double,double);
+ static Matrix combine(double,double,double);
static bool inside_domain(MElementOctree*,double,double);
static double get_ratio(GFace*,SPoint2);
static void print_field1();
- static void print_field2();
+ static void print_field2(GRegion*);
static void print_segment(SPoint3,SPoint3,std::ofstream&);
+ static GRegion* test();
static void clear();
};
@@ -66,10 +68,29 @@ class Size_field{
static MElementOctree* octree;
Size_field();
public:
- static void init_model();
- static void solve();
+ static void init_region(GRegion*);
+ static void solve(GRegion*);
static double search(double,double,double);
static double get_ratio(GFace*,SPoint2);
static void print_field();
+ static GRegion* test();
static void clear();
};
+
+class Nearest_point{
+ private:
+ static std::vector<SPoint3> random;
+ #if defined(HAVE_ANN)
+ static ANNpointArray duplicate;
+ static ANNkd_tree* kd_tree;
+ #endif
+ Nearest_point();
+ public:
+ static void init_region(GRegion*);
+ static bool search(double,double,double,SVector3&);
+ static double T(double,double,double,double,double);
+ static void print_field(GRegion*);
+ static void print_segment(SPoint3,SPoint3,std::ofstream&);
+ static GRegion* test();
+ static void clear();
+};
\ No newline at end of file
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index cdd35e854639e22f91318799eeac4f730291ffad..5819da4a4ee0a0de34e1656e5ede0353d9f3c65a 100644
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -324,7 +324,7 @@ void Filler::treat_model(){
}
void Filler::treat_region(GRegion* gr){
-#if defined(HAVE_RTREE)
+ #if defined(HAVE_RTREE)
unsigned int i;
int j;
int count;
@@ -345,7 +345,7 @@ void Filler::treat_region(GRegion* gr){
std::set<MVertex*>::iterator it;
RTree<Node*,double,3,double> rtree;
- Frame_field::init_model();
+ Frame_field::init_region(gr);
octree = new MElementOctree(gr->model());
for(i=0;i<gr->getNumMeshElements();i++){
@@ -436,7 +436,7 @@ void Filler::treat_region(GRegion* gr){
for(i=0;i<new_vertices.size();i++) delete new_vertices[i];
new_vertices.clear();
Frame_field::clear();
-#endif
+ #endif
}
Metric Filler::get_metric(double x,double y,double z){