*** empty log message ***

Mesh/BDS.cpp

@@ -12,6 +12,27 @@
   
 
 */
+double BDS_Quadric:: normalCurv ( double x, double y, double z ) const
+{
+    // K = div n = div ( Grad(f) / ||Grad(f)||)
+    //   = Laplacian (f) /  ||Grad(f)|| - Grad^t(f) Hessian Grad(f) / ||Grad(f)||^3
+    BDS_Vector g = Gradient ( x , y , z );
+    double normG = sqrt ( g.x*g.x +g.y*g.y +g.z*g.z); 
+
+    double c1 = 2*(a + b + c)/normG;
+
+
+    double c2 = 
+	2*a * g.x * g.x +  
+	2*b * g.y * g.y +  
+	2*c * g.z * g.z +
+	4*d *  g.x * g.y +    
+	4*e *  g.x * g.z +  
+	4*f *  g.y * g.z ;
+
+    //printf("%g %g %g %g %g %g\n",a,b,c,normG,c1, c2 / (normG*normG*normG));
+    return fabs(c1 - c2 / (normG*normG*normG));
+}
 
 void BDS_Quadric::projection ( double xa, double ya, double za,
 			       double &x, double &y, double &z) const 
@@ -47,7 +68,7 @@ void BDS_Quadric::projection ( double xa, double ya, double za,
 	    grad.x * g +grad.y * h +grad.z * i;
 	double delta = coef2*coef2 - 4.*coef1*residual;
 	if (delta < 0){
-	    printf ("aaargh error projection");
+//	    printf ("aaargh error projection");
 	}
 	else
 	{
@@ -432,9 +453,9 @@ void BDS_Mesh :: reverseEngineerCAD ( )
 		    int k=0;
 		    while (pit!=pite)
 		    {
-			PLANE  ( k , 0 ) = (*pit).x;
-			PLANE  ( k , 1 ) = (*pit).y;
-			PLANE  ( k , 2 ) = (*pit).z;
+			PLANE  ( k , 0 ) = (*pit).x + (rand() % 1000) * LC / 1.e12;
+			PLANE  ( k , 1 ) = (*pit).y+ (rand() % 1000) * LC / 1.e12;
+			PLANE  ( k , 2 ) = (*pit).z+ (rand() % 1000) * LC / 1.e12;
 			ONES  ( k )     = -1;
 			k++;
 			++pit;
@@ -447,7 +468,7 @@ void BDS_Mesh :: reverseEngineerCAD ( )
 		    for (int i=0;i<pts.size();i++)
 		    {
 			const double dist = PLANE  ( i , 0 ) * RSLT(0)+PLANE  ( i , 1 ) * RSLT(1)+PLANE  ( i , 2 ) * RSLT(2)+1;
-			if (fabs(dist) > 1.e-5 * LC)plane = false;
+			if (fabs(dist) > 1.e-5 * LC * sqrt (RSLT(0)*RSLT(0)+RSLT(1)*RSLT(1)+RSLT(2)*RSLT(2)) )plane = false;
 		    }
 		    
 		    if ( plane ) 
@@ -456,7 +477,7 @@ void BDS_Mesh :: reverseEngineerCAD ( )
 			(*it)->surf = new BDS_Plane (RSLT(0),RSLT(1),RSLT(2));
 		    }
 		}
-		if (!(*it)->surf && pts.size() > 10)
+		if (!(*it)->surf && pts.size() > 20)
 		{
 		    printf("coucou quadrique\n");
 		    Double_Matrix QUADRIC  ( pts.size() , 9 );
@@ -482,20 +503,20 @@ void BDS_Mesh :: reverseEngineerCAD ( )
 		    }
 		    QUADRIC.least_squares (ONES, RSLT);
 		    bool quad = true;
-		    for (int i=0;i<pts.size();i++)
+		    BDS_Quadric temp ( RSLT(0),RSLT(1),RSLT(2),RSLT(3),RSLT(4),RSLT(5),RSLT(6),RSLT(7),RSLT(8));
+		    pit  = pts.begin();
+		    while (pit!=pite)
 		    {
-			const double dist = 
-			    QUADRIC  ( i , 0 ) * RSLT(0) +
-			    QUADRIC  ( i , 1 ) * RSLT(1) +
-			    QUADRIC  ( i , 2 ) * RSLT(2) +
-			    QUADRIC  ( i , 3 ) * RSLT(3) +
-			    QUADRIC  ( i , 4 ) * RSLT(4) +
-			    QUADRIC  ( i , 5 ) * RSLT(5) +
-			    QUADRIC  ( i , 6 ) * RSLT(6) +
-			    QUADRIC  ( i , 7 ) * RSLT(7) +
-			    QUADRIC  ( i , 8 ) * RSLT(8) - 1;
-//			printf("dist = %g LC %g\n",dist,LC);
-			if (fabs(dist) > 1.e-3 * LC )quad = false;
+			double x,y,z;
+			temp.projection ((*pit).x , (*pit).y , (*pit).z ,
+					 x,y,z);
+			const double dist = sqrt ( ((*pit).x-x)*((*pit).x-x) +
+						   ((*pit).y-y)*((*pit).y-y) +
+						   ((*pit).z-z)*((*pit).z-z));
+//			printf("%g %g %g ... %g %g %g\n",(*pit).x , (*pit).y , (*pit).z ,x,y,z);
+			       
+			if (dist > 1.e-2 * LC )quad = false;
+			++pit;
 		    }		    
 		    if ( quad ) 
 		    {
@@ -1860,6 +1881,7 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
 {
     int nb_modif = 0;
 
+    double CURV_FACTOR = 5;
 
     // add initial set of edges in a list
     std::set<BDS_Edge*, EdgeLessThan> small_to_long;
@@ -1887,12 +1909,26 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
 	      double ll = sqrt((op[0]->X-op[1]->X)*(op[0]->X-op[1]->X)+
 			       (op[0]->Y-op[1]->Y)*(op[0]->Y-op[1]->Y)+
 			       (op[0]->Z-op[1]->Z)*(op[0]->Z-op[1]->Z));
-	      double curv = 0.5 * (*it)->p1->radius_of_curvature + 
-		  0.5 * (*it)->p2->radius_of_curvature;
 	      double LLL = l;
-///	      if (l > .50 * curv ) LLL = .50*curv;
+
+	      if ((*it)->g && (*it)->g->classif_degree == 2 &&(*it)->g->surf)
+	      {
+		  double curvature = (*it)->g->surf->normalCurv (0.5*((*it)->p1->X+(*it)->p2->X),
+								 0.5*((*it)->p1->Y+(*it)->p2->Y),
+								 0.5*((*it)->p1->Z+(*it)->p2->Z));
+		  if (curvature > 0.0)
+		  {
+
+		      double radius_of_curvature = 1./curvature;
+//		      printf("surface %d radius %g\n",(*it)->g->classif_tag,radius_of_curvature);
+		      double maxL = radius_of_curvature / CURV_FACTOR;
+
+		      if ( LLL > maxL ) LLL = maxL;
+		  }
+	      }
 
 	      if (!(*it)->deleted && (*it)->length() > LLL/0.7 && ll > 0.25 * LLL){
+//		  printf("we schpilt\n");
 		  split_edge (*it, 0.5 );
 		  nb_modif++;
 	      }
@@ -1919,10 +1955,20 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
 	std::set<BDS_Edge*, EdgeLessThan>::iterator ite  = small_to_long.end();
 	while (it != ite)
 	{
-	    double curv = 0.5 * (*it)->p1->radius_of_curvature + 
-		0.5 * (*it)->p2->radius_of_curvature;
+
 	    double LLL = l;
-//	    if (l > .50 * curv ) LLL = .50*curv;
+	    if ((*it)->g && (*it)->g->classif_degree == 2 &&(*it)->g->surf)
+	    {
+		double curvature = (*it)->g->surf->normalCurv (0.5*((*it)->p1->X+(*it)->p2->X),
+							       0.5*((*it)->p1->Y+(*it)->p2->Y),
+							       0.5*((*it)->p1->Z+(*it)->p2->Z));
+		if (curvature > 0.0)
+		{
+		    double radius_of_curvature = 1./curvature;
+		    double maxL = radius_of_curvature / CURV_FACTOR;
+		    if ( LLL > maxL ) LLL = maxL;
+		}
+	    }
 
 	    if (!(*it)->deleted && (*it)->length() < 0.7 * LLL ){
 		if (nb_modif %2 )
@@ -2033,6 +2079,8 @@ BDS_Mesh::BDS_Mesh (const BDS_Mesh &other)
 	BDS_Point *p =  add_point((*it)->iD,(*it)->X,(*it)->Y,(*it)->Z);
 	p->g = ((*it)->g)? get_geom  ((*it)->g->classif_tag,(*it)->g->classif_degree) : 0;
 	p->radius_of_curvature = (*it)->radius_of_curvature ;
+	if (p->g->classif_degree == 0)
+	    p->g->p = p;
     }
     for ( std::set<BDS_Edge*, EdgeLessThan>::const_iterator it  = other.edges.begin();
 	  it != other.edges.end();
@@ -2040,6 +2088,8 @@ BDS_Mesh::BDS_Mesh (const BDS_Mesh &other)
     {
 	BDS_Edge  * e = add_edge ((*it)->p1->iD,(*it)->p2->iD);
 	e->g = ((*it)->g)? get_geom  ((*it)->g->classif_tag,(*it)->g->classif_degree) : 0;
+	if (e->g->classif_degree == 1)
+	    e->g->e.push_back(e);
     }
     for (std::list<BDS_Triangle*>::const_iterator it  = other.triangles.begin();
 	 it != other.triangles.end();
@@ -2049,5 +2099,6 @@ BDS_Mesh::BDS_Mesh (const BDS_Mesh &other)
         (*it)->getNodes (n);
 	BDS_Triangle *t = add_triangle(n[0]->iD,n[1]->iD,n[2]->iD);
 	t->g = get_geom  ((*it)->g->classif_tag,(*it)->g->classif_degree);
+	t->g->t.push_back(t);
     }
 }

Mesh/BDS.h

@@ -23,6 +23,8 @@ public :
     virtual void projection ( double xa, double ya, double za,
 			      double &x, double &y, double &z) const =0;
     virtual std::string nameOf () const = 0;
+//    virtual BDS_Vector Gradient ( double x, double y, double z ) const = 0;
+    virtual double normalCurv ( double x, double y, double z ) const = 0;
 };
 
 
@@ -295,7 +297,7 @@ public:
 
 class BDS_Plane : public  BDS_Surface
 {
-    double a,b,c;
+    double a,b,c,d;
 public :
     BDS_Plane (const double &A, const double &B, const double &C)
 	: a(A),b(B),c(C)
@@ -311,7 +313,14 @@ public :
 	    z = za + k * c;
 	}
     virtual std::string nameOf () const {return std::string("Plane");}
-
+    virtual BDS_Vector Gradient ( double x, double y, double z ) const
+	{
+	    return BDS_Vector ( a , b , c );
+	} 
+    virtual double normalCurv ( double x, double y, double z ) const
+	{
+	    return 0.0;
+	}
 };
 
 class BDS_Quadric : public  BDS_Surface
@@ -330,6 +339,8 @@ public :
 				2* ( e * x + f * y + c * z ) + i );
 	}
 
+    virtual double normalCurv ( double x, double y, double z ) const;
+
     virtual double signedDistanceTo (  double x, double y, double z ) const {
 	const double q = 
 	    a * x * x +  

Mesh/DiscreteSurface.cpp

-// $Id: DiscreteSurface.cpp,v 1.17 2005-07-04 15:07:41 remacle Exp $
+// $Id: DiscreteSurface.cpp,v 1.18 2005-07-07 07:19:27 remacle Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -63,6 +63,8 @@ void Mesh_To_BDS(Surface *s, BDS_Mesh *m)
 
 void BDS_To_Mesh_2(Mesh *m)
 {
+    Msg(STATUS2, "Moving the surface mesh in the old gmsh structure\n");
+    
     Tree_Action(m->Vertices, Free_Vertex);  
     Tree_Delete(m->Vertices);
     m->Vertices = Tree_Create(sizeof(Vertex *), compareVertex);
@@ -77,6 +79,7 @@ void BDS_To_Mesh_2(Mesh *m)
 	    ++it;
 	}
     }
+
     {
 	std::set<BDS_Edge*, EdgeLessThan>::iterator it  = m->bds_mesh->edges.begin();
 	std::set<BDS_Edge*, EdgeLessThan>::iterator ite = m->bds_mesh->edges.end();
@@ -89,6 +92,7 @@ void BDS_To_Mesh_2(Mesh *m)
 		Vertex *v2 = FindVertex((*it)->p2->iD, m);
 		SimplexBase *simp = Create_SimplexBase(v1,v2,NULL, NULL);
 		Curve *c = FindCurve (g->classif_tag,m);
+		if (c)
 		    simp->iEnt = g->classif_tag;
 		Tree_Insert(c->SimplexesBase, &simp);
 	    }
@@ -107,13 +111,15 @@ void BDS_To_Mesh_2(Mesh *m)
 	    SimplexBase *simp = Create_SimplexBase(v1,v2,v3, NULL);
 	    BDS_GeomEntity *g = (*it)->g;
 	    Surface *s = FindSurface (g->classif_tag,m);
+	    if(s)
 		simp->iEnt = g->classif_tag;
-	    Tree_Add(s->Simplexes, &simp);
+	    else
+		printf("argh\n");
+	    Tree_Add(s->SimplexesBase, &simp);
 	    ++it;
 	}
     }
-
-
+    Msg(STATUS2, "Ready");
 }
 void BDS_To_Mesh(Mesh *m)
 {
@@ -166,14 +172,13 @@ void BDS_To_Mesh(Mesh *m)
 
     CTX.mesh.changed = 1;
 
-    printf ("%d surfaces %d curves\n",Tree_Nbr(m->Surfaces), Tree_Nbr(m->Curves) );
-
 }
 
 // Public interface for discrete surface/curve mesh algo
 
 int MeshDiscreteSurface(Surface *s)
 { 
+    Msg(STATUS2, "Discrete Surface Mesh Generator...");
     if(s->bds){
 	// s->bds is the discrete surface that defines the geometry
 	if(!THEM->bds_mesh){
@@ -181,12 +186,12 @@ int MeshDiscreteSurface(Surface *s)
 	    int iter = 0;
 	    while(iter < 20 && THEM->bds_mesh->adapt_mesh(CTX.mesh.lc_factor * THEM->bds->LC, 
 						    true, THEM->bds)){
-	printf("iter %d done\n",iter);
+		Msg(STATUS2, "Iteration %2d/20 done (%d triangles)\n",iter, THEM->bds_mesh->triangles.size());
 		iter ++;
 	    }
 	    BDS_To_Mesh_2(THEM);
-      printf("mesh has %d vertices (%d)\n",Tree_Nbr(THEM->Vertices),THEM->bds->points.size());
-      THEM->bds_mesh->save_gmsh_format ( "3.msh" );
+	    Msg(STATUS2, "Mesh has %d vertices (%d)\n",Tree_Nbr(THEM->Vertices),THEM->bds->points.size());
+//	    THEM->bds_mesh->save_gmsh_format ( "3.msh" );
 	}
 	return 1;
     }