reverting to the isotropic version of the RTree algorithm

Mesh/simple3D.cpp

@@ -17,7 +17,9 @@
 #include "rtree.h"
 #endif 
 
-#define baseline 0.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 sqrt3 1.73205081
 
 /*********definitions*********/
 
@@ -49,46 +51,43 @@ class Metric{
 
 class Node{
  private:
-  double h1;
-  double h2;
-  double h3;
+  double h;
   Metric m;
   SPoint3 point;
  public:
+  double min[3];
+  double max[3];
   Node();
   Node(SPoint3);
   ~Node();
-  void set_h1(double);
-  void set_h2(double);
-  void set_h3(double);
+  void set_size(double);
   void set_metric(Metric);
   void set_point(SPoint3);
-  double get_h1();
-  double get_h2();
-  double get_h3();
+  double get_size();
   Metric get_metric();
   SPoint3 get_point();
-  double min[3];
-  double max[3];
 };
 
 class Wrapper{
  private:
   bool too_close;
   Node* spawn;
+  Node* parent;
  public:
   Wrapper();
-  Wrapper(Node*);
+  Wrapper(Node*,Node*);
   ~Wrapper();
   void set_too_close(bool);
   void set_spawn(Node*);
+  void set_parent(Node*);
   bool get_too_close();
   Node* get_spawn();
+  Node* get_parent();
 };
 
 /*********functions*********/
 
-double infinity_distance(SPoint3 p1,SPoint3 p2,Metric m,double h1,double h2,double h3){
+double infinity_distance(SPoint3 p1,SPoint3 p2,Metric m){
   double distance;
   double x1,y1,z1;
   double x2,y2,z2;
@@ -103,39 +102,38 @@ double infinity_distance(SPoint3 p1,SPoint3 p2,Metric m,double h1,double h2,doub
   g = m.get_m13();
   h = m.get_m23();
   i = m.get_m33();
+		
+  x1 = a*p1.x() + b*p1.y() + c*p1.z();
+  y1 = d*p1.x() + e*p1.y() + f*p1.z();
+  z1 = g*p1.x() + h*p1.y() + i*p1.z();
 	
-  x1 = (1.0/h1)*(a*p1.x() + b*p1.y() + c*p1.z());
-  y1 = (1.0/h2)*(d*p1.x() + e*p1.y() + f*p1.z());
-  z1 = (1.0/h3)*(g*p1.x() + h*p1.y() + i*p1.z());
-	
-  x2 = (1.0/h1)*(a*p2.x() + b*p2.y() + c*p2.z());
-  y2 = (1.0/h2)*(d*p2.x() + e*p2.y() + f*p2.z());
-  z2 = (1.0/h3)*(g*p2.x() + h*p2.y() + i*p2.z());
+  x2 = a*p2.x() + b*p2.y() + c*p2.z();
+  y2 = d*p2.x() + e*p2.y() + f*p2.z();
+  z2 = g*p2.x() + h*p2.y() + i*p2.z();
 	
   distance = std::max(std::max(fabs(x2-x1),fabs(y2-y1)),fabs(z2-z1)); 		
   return distance;
 }
 
 bool rtree_callback(Node* neighbour,void* w){
-  double h1;
-  double h2;
-  double h3;
+  double h;
   double distance;
   Metric m;
-  Node *spawn;
+  Node *spawn,*parent;
   Wrapper* wrapper;
 	
   wrapper = static_cast<Wrapper*>(w);
   spawn = wrapper->get_spawn();
-  h1 = spawn->get_h1();
-  h2 = spawn->get_h2();
-  h3 = spawn->get_h3();
+  parent = wrapper->get_parent();
+  h = spawn->get_size();
   m = spawn->get_metric();
-
-  distance = infinity_distance(spawn->get_point(),neighbour->get_point(),m,h1,h2,h3);
-  if(distance<baseline){
-    wrapper->set_too_close(1);
-	return false;
+	
+  if(neighbour!=parent){
+    distance = infinity_distance(spawn->get_point(),neighbour->get_point(),m);
+	if(distance<k1*h){
+	  wrapper->set_too_close(1);
+	  return false;
+	}
   }
 	
   return true;
@@ -239,16 +237,8 @@ Node::Node(SPoint3 new_point){
 
 Node::~Node(){}
 
-void Node::set_h1(double new_h1){
-  h1 = new_h1;
-}
-
-void Node::set_h2(double new_h2){
-  h2 = new_h2;
-}
-
-void Node::set_h3(double new_h3){
-  h3 = new_h3;
+void Node::set_size(double new_h){
+  h = new_h;
 }
 
 void Node::set_metric(Metric new_m){
@@ -259,16 +249,8 @@ void Node::set_point(SPoint3 new_point){
   point = new_point;
 }
 
-double Node::get_h1(){
-  return h1;
-}
-
-double Node::get_h2(){
-  return h2;
-}
-
-double Node::get_h3(){
-  return h3;
+double Node::get_size(){
+  return h;
 }
 
 Metric Node::get_metric(){
@@ -285,9 +267,10 @@ Wrapper::Wrapper(){
   too_close = 0;
 }
 
-Wrapper::Wrapper(Node* new_spawn){
+Wrapper::Wrapper(Node* new_spawn,Node* new_parent){
   too_close = 0;
   spawn = new_spawn;
+  parent = new_parent;
 }
 
 Wrapper::~Wrapper(){}
@@ -300,6 +283,10 @@ void Wrapper::set_spawn(Node* new_spawn){
   spawn = new_spawn;
 }
 
+void Wrapper::set_parent(Node* new_parent){
+  parent = new_parent;
+}
+
 bool Wrapper::get_too_close(){
   return too_close;
 }
@@ -308,6 +295,10 @@ Node* Wrapper::get_spawn(){
   return spawn;
 }
 
+Node* Wrapper::get_parent(){
+  return parent;
+}
+
 /*********class Filler*********/
 
 Filler::Filler(){}
@@ -335,8 +326,7 @@ void Filler::treat_region(GRegion* gr){
   bool ok;
   double x,y,z;
   SPoint3 point;
-  Node *node,*spawn,*parent;
-  Node *n1,*n2,*n3,*n4,*n5,*n6;
+  Node *node,*spawn,*parent,*n1,*n2,*n3,*n4,*n5,*n6;
   MVertex* vertex;
   MElement* element;
   MElementOctree* octree;
@@ -391,7 +381,7 @@ void Filler::treat_region(GRegion* gr){
 	n4 = new Node();
 	n5 = new Node();
 	n6 = new Node();
-	offsprings(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);
@@ -408,16 +398,19 @@ void Filler::treat_region(GRegion* gr){
 	  z = point.z();
 	  if(inside_domain(octree,x,y,z)){
 		compute_parameters(spawn,gr);
-		wrapper.set_too_close(0);
-		wrapper.set_spawn(spawn);
-		rtree.Search(spawn->min,spawn->max,rtree_callback,&wrapper);
-		if(!wrapper.get_too_close()){
-	      fifo.push(spawn);
-		  rtree.Insert(spawn->min,spawn->max,spawn);
-		  vertex = new MVertex(x,y,z,gr,0);
-		  new_vertices.push_back(vertex);
-		  ok = 1;
-		}
+	    if(far_from_boundary(octree,spawn)){
+		  wrapper.set_too_close(0);
+		  wrapper.set_spawn(spawn);
+		  wrapper.set_parent(parent);
+	      rtree.Search(spawn->min,spawn->max,rtree_callback,&wrapper);
+		  if(!wrapper.get_too_close()){
+		    fifo.push(spawn);
+		    rtree.Insert(spawn->min,spawn->max,spawn);
+			vertex = new MVertex(x,y,z,gr,0);
+			new_vertices.push_back(vertex);
+			ok = 1;
+		  }
+	    }
 	  }
 	  if(!ok) delete spawn;
 	}
@@ -461,15 +454,7 @@ Metric Filler::get_metric(double x,double y,double z){
   return m;
 }
 
-double Filler::get_h1(double x,double y,double z){
-  return 0.25;
-}
-
-double Filler::get_h2(double x,double y,double z){
-  return 0.25;
-}
-
-double Filler::get_h3(double x,double y,double z){
+double Filler::get_size(double x,double y,double z){
   return 0.25;
 }
 
@@ -503,7 +488,7 @@ Metric Filler::get_clf_metric(double x,double y,double z,GEntity* ge){
   return m;
 }
 
-double Filler::get_clf_h(double x,double y,double z,GEntity* ge){
+double Filler::get_clf_size(double x,double y,double z,GEntity* ge){
   double h;
   Field* field;
   FieldManager* manager;
@@ -527,11 +512,32 @@ bool Filler::inside_domain(MElementOctree* octree,double x,double y,double z){
   else return 0;
 }
 
+bool Filler::far_from_boundary(MElementOctree* octree,Node* node){
+  double x,y,z;
+  double h;
+  SPoint3 point;
+  MElement *e1,*e2,*e3,*e4,*e5,*e6;
+	
+  point = node->get_point();
+  x = point.x();
+  y = point.y();
+  z = point.z();
+  h = node->get_size();
+	
+  e1 = (MElement*)octree->find(x+k2*h,y,z,3,true);
+  e2 = (MElement*)octree->find(x-k2*h,y,z,3,true);
+  e3 = (MElement*)octree->find(x,y+k2*h,z,3,true);
+  e4 = (MElement*)octree->find(x,y-k2*h,z,3,true);
+  e5 = (MElement*)octree->find(x,y,z+k2*h,3,true);
+  e6 = (MElement*)octree->find(x,y,z-k2*h,3,true);
+	
+  if(e1!=NULL && e2!=NULL && e3!=NULL && e4!=NULL && e5!=NULL && e6!=NULL) return 1;
+  else return 0;
+}
+
 void Filler::compute_parameters(Node* node,GEntity* ge){
   double x,y,z;
-  double h1;
-  double h2;
-  double h3;
+  double h;
   Metric m;
   SPoint3 point;
 	
@@ -540,23 +546,19 @@ void Filler::compute_parameters(Node* node,GEntity* ge){
   y = point.y();
   z = point.z();
   m = get_metric(x,y,z);
-  h1 = get_h1(x,y,z);
-  h2 = get_h2(x,y,z);
-  h3 = get_h3(x,y,z);
+  h = get_size(x,y,z);
   
+  node->set_size(h);
   node->set_metric(m);
-  node->set_h1(h1);
-  node->set_h2(h2);
-  node->set_h3(h3);
-  node->min[0] = x - sqrt(h1*h1 + h2*h2 + h3*h3);
-  node->min[1] = y - sqrt(h1*h1 + h2*h2 + h3*h3);
-  node->min[2] = z - sqrt(h1*h1 + h2*h2 + h3*h3);
-  node->max[0] = x + sqrt(h1*h1 + h2*h2 + h3*h3);
-  node->max[1] = y + sqrt(h1*h1 + h2*h2 + h3*h3);
-  node->max[2] = z + sqrt(h1*h1 + h2*h2 + h3*h3);
-}
-
-void Filler::offsprings(Node* node,Node* n1,Node* n2,Node* n3,Node* n4,Node* n5,Node* n6){
+  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(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;
@@ -564,9 +566,8 @@ void Filler::offsprings(Node* node,Node* n1,Node* n2,Node* n3,Node* n4,Node* n5,
   double x4,y4,z4;
   double x5,y5,z5;
   double x6,y6,z6;
-  double h1;
-  double h2;
-  double h3;
+  double h;
+  double h1,h2,h3,h4,h5,h6;
   Metric m;
   SPoint3 point;
 	
@@ -574,34 +575,38 @@ void Filler::offsprings(Node* node,Node* n1,Node* n2,Node* n3,Node* n4,Node* n5,
   x = point.x();
   y = point.y();
   z = point.z();
-  h1 = node->get_h1();
-  h2 = node->get_h2();
-  h3 = node->get_h3();
+  h = node->get_size();
   m = node->get_metric();
 	
+  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();
   
-  x2 = x - h1*m.get_m11();
-  y2 = y - h1*m.get_m21();
-  z2 = z - h1*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();
   
-  x3 = x + h2*m.get_m12();
-  y3 = y + h2*m.get_m22();
-  z3 = z + h2*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();
   
-  x4 = x - h2*m.get_m12();
-  y4 = y - h2*m.get_m22();
-  z4 = z - h2*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();
   
-  x5 = x + h3*m.get_m13();
-  y5 = y + h3*m.get_m23();
-  z5 = z + h3*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();
   
-  x6 = x - h3*m.get_m13();
-  y6 = y - h3*m.get_m23();
-  z6 = z - h3*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();
 	
   *n1 = Node(SPoint3(x1,y1,z1));
   *n2 = Node(SPoint3(x2,y2,z2));
@@ -611,6 +616,31 @@ void Filler::offsprings(Node* node,Node* n1,Node* n2,Node* n3,Node* n4,Node* n5,
   *n6 = Node(SPoint3(x6,y6,z6));
 }
 
+double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,double h_nearer,SVector3 direction){
+  double x,y,z;
+  double average;
+  double h_farther;
+  double coeffA,coeffB;
+	
+  x = point.x() + h_nearer*direction.x();
+  y = point.y() + h_nearer*direction.y();
+  z = point.z() + h_nearer*direction.z();
+  if(inside_domain(octree,x,y,z)){
+    h_farther = get_size(x,y,z);
+  }
+  else h_farther = h_nearer;
+  
+  coeffA = 1.0;
+  coeffB = 0.2;
+  if(h_farther>h_nearer){
+    average = coeffA*h_nearer + (1.0-coeffA)*h_farther;
+  }
+  else{
+    average = coeffB*h_nearer + (1.0-coeffB)*h_farther;
+  }
+  return average;
+}
+
 int Filler::get_nbr_new_vertices(){
   return new_vertices.size();
 }
@@ -634,9 +664,7 @@ void Filler::print_node(Node* node,std::ofstream& file){
   double x4,y4,z4;
   double x5,y5,z5;
   double x6,y6,z6;
-  double h1;
-  double h2;
-  double h3;
+  double h;
   Metric m;
   SPoint3 point;
 	
@@ -644,34 +672,32 @@ void Filler::print_node(Node* node,std::ofstream& file){
   x = point.x();
   y = point.y();
   z = point.z();
-  h1 = node->get_h1();
-  h2 = node->get_h2();
-  h3 = node->get_h3();
+  h = node->get_size();
   m = node->get_metric();
 	
-  x1 = x + h1*m.get_m11();
-  y1 = y + h1*m.get_m21();
-  z1 = z + h1*m.get_m31();
+  x1 = x + h*m.get_m11();
+  y1 = y + h*m.get_m21();
+  z1 = z + h*m.get_m31();
   
-  x2 = x - h1*m.get_m11();
-  y2 = y - h1*m.get_m21();
-  z2 = z - h1*m.get_m31();
+  x2 = x - h*m.get_m11();
+  y2 = y - h*m.get_m21();
+  z2 = z - h*m.get_m31();
   
-  x3 = x + h2*m.get_m12();
-  y3 = y + h2*m.get_m22();
-  z3 = z + h2*m.get_m32();
+  x3 = x + h*m.get_m12();
+  y3 = y + h*m.get_m22();
+  z3 = z + h*m.get_m32();
   
-  x4 = x - h2*m.get_m12();
-  y4 = y - h2*m.get_m22();
-  z4 = z - h2*m.get_m32();
+  x4 = x - h*m.get_m12();
+  y4 = y - h*m.get_m22();
+  z4 = z - h*m.get_m32();
   
-  x5 = x + h3*m.get_m13();
-  y5 = y + h3*m.get_m23();
-  z5 = z + h3*m.get_m33();
+  x5 = x + h*m.get_m13();
+  y5 = y + h*m.get_m23();
+  z5 = z + h*m.get_m33();
   
-  x6 = x - h3*m.get_m13();
-  y6 = y - h3*m.get_m23();
-  z6 = z - h3*m.get_m33();
+  x6 = x - h*m.get_m13();
+  y6 = y - h*m.get_m23();
+  z6 = z - h*m.get_m33();
 
   print_segment(SPoint3(x,y,z),SPoint3(x1,y1,z1),file);
   print_segment(SPoint3(x,y,z),SPoint3(x2,y2,z2),file);

Mesh/simple3D.h

@@ -15,14 +15,14 @@ class Filler{
  private:
   static std::vector<MVertex*> new_vertices;
   Metric get_metric(double,double,double);
-  double get_h1(double,double,double);
-  double get_h2(double,double,double);
-  double get_h3(double,double,double);
+  double get_size(double,double,double);
   Metric get_clf_metric(double,double,double,GEntity*);
-  double get_clf_h(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*,GEntity*);
-  void offsprings(Node*,Node*,Node*,Node*,Node*,Node*,Node*);
+  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: