3D fields

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(){

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;
+MElementOctree* Size_field::octree;
\ No newline at end of file
+
+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

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

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){