*** empty log message ***

Geo/GEdge.cpp

-// $Id: GEdge.cpp,v 1.16 2006-10-12 01:35:32 geuzaine Exp $
+// $Id: GEdge.cpp,v 1.17 2006-11-15 20:46:46 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -119,3 +119,30 @@ SVector3 GEdge::secondDer(double par) const
   SVector3 x2 = firstDer(par+eps);
   return 500*(x2-x1);
 }
+
+double GEdge::curvature(double par) const 
+{
+  double eps1 = 1.e-3;
+  double eps2 = 1.e-3;
+
+  Range<double> r = parBounds(0);
+  if (r.low() == par) eps2 = 0;
+  if (r.high() == par) eps1 = 0;
+
+
+  SVector3 n1 = firstDer(par-eps1);
+  SVector3 n2 = firstDer(par+eps2);
+
+  GPoint P1 = point(par-eps1);
+  GPoint P2 = point(par+eps2);
+
+  double D = sqrt ( (P1.x()-P2.x())*(P1.x()-P2.x())+
+		    (P1.y()-P2.y())*(P1.y()-P2.y())+
+		    (P1.z()-P2.z())*(P1.z()-P2.z()));  
+
+  n1.normalize();
+  n2.normalize();
+  const double one_over_D = 1./D;
+  SVector3 d = one_over_D*(n2-n1);
+  return norm(d);
+}

Geo/GEdge.h

@@ -73,6 +73,7 @@ class GEdge : public GEntity {
   // Get second derivative of edge at the given parameter.
   // Default implentation using central differences
   virtual SVector3 secondDer(double par) const ;
+  virtual double curvature (double par) const;  
 
   // Reparmaterize the point onto the given face.
   virtual SPoint2 reparamOnFace(GFace *face, double epar,int dir) const = 0;

Geo/GFace.cpp

-// $Id: GFace.cpp,v 1.18 2006-11-14 20:20:18 remacle Exp $
+// $Id: GFace.cpp,v 1.19 2006-11-15 20:46:46 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -481,3 +481,11 @@ SPoint2 GFace::parFromPoint(const SPoint3 &p) const
 
   return SPoint2(U,V);
 }
+
+void GFace::parFromPoint(const SPoint3 &p, std::list<double> &u, std::list<double> &v ) const
+{
+  double U,V;  
+  XYZtoUV(p.x(),p.y(),p.z(),U,V,1.0);
+  u.push_back(U);
+  v.push_back(V);
+}

Geo/GFace.h

@@ -28,6 +28,8 @@
 #include "Pair.h"
 #include "ExtrudeParams.h"
 
+#include <set>
+
 struct mean_plane
 {
   double plan[3][3];
@@ -41,6 +43,11 @@ class GRegion;
 class GFace : public GEntity 
 {
  protected: 
+  // in a manifold representation, edges can
+  // be taken twice in the topology in order 
+  // to represent the periodicity of the surface
+  std::set<GEdge *> edges_taken_twice;
+  // list of al the edges of the face
   std::list<GEdge *> l_edges;
   std::list<int> l_dirs;
   GRegion *r1, *r2;
@@ -86,7 +93,8 @@ class GFace : public GEntity
 
   // Return the parmater location on the face given a point in space
   // that is on the face.
-  virtual SPoint2 parFromPoint(const SPoint3 &) const;
+  virtual SPoint2            parFromPoint(const SPoint3 &) const;
+  virtual void               parFromPoint(const SPoint3 &, std::list<double> &u, std::list<double> &v ) const;
 
   // True if the parameter value is interior to the face.
   virtual int containsParam(const SPoint2 &pt) const = 0;

Geo/OCCEdge.cpp

-// $Id: OCCEdge.cpp,v 1.3 2006-11-15 15:06:45 geuzaine Exp $
+// $Id: OCCEdge.cpp,v 1.4 2006-11-15 20:46:46 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -75,11 +75,21 @@ GEntity::GeomType OCCEdge::geomType() const
 
 int OCCEdge::minimumMeshSegments () const
 {
-  return CTX.mesh.min_circ_points; ;
+  return 2 ;
 }
 
 int OCCEdge::minimumDrawSegments () const
 {
   return CTX.geom.circle_points;
 }
+
+double OCCEdge::curvature(double par) const 
+{
+  return GEdge::curvature(par);
+  BRepAdaptor_Curve brepc(c);
+  BRepLProp_CLProps prop(brepc, 1, 1e-5);
+  prop.SetParameter (par); 
+  printf("curvature = %12.5E\n",prop.Curvature()); 
+  return fabs(prop.Curvature());
+}
 #endif

Geo/OCCEdge.h

@@ -45,6 +45,7 @@ class OCCEdge : public GEdge {
   virtual int containsPoint(const SPoint3 &pt) const { throw; }
   virtual int containsParam(double pt) const;
   virtual SVector3 firstDer(double par) const;
+  virtual double curvature (double par) const;
   virtual SPoint2 reparamOnFace(GFace * face, double epar, int dir) const { throw; }
   ModelType getNativeType() const { return OpenCascadeModel; }
   void * getNativePtr() const { return (void*) &c; }

Geo/OCCFace.cpp

-// $Id: OCCFace.cpp,v 1.4 2006-11-15 15:06:45 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.5 2006-11-15 20:46:46 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -27,7 +27,7 @@
 #include "Message.h"
 
 OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape &emap)
-  : GFace(m, num), s(_s)
+  : GFace(m, num), s(_s),_periodic(false)
 {
   TopExp_Explorer exp0, exp01, exp1, exp2, exp3;
   for (exp2.Init (s, TopAbs_WIRE); exp2.More(); exp2.Next())
@@ -35,22 +35,28 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
       TopoDS_Shape wire = exp2.Current();
       Msg(INFO,"OCC Face %d - New Wire",num);
       std::list<GEdge*> l_wire;
+      std::set<GEdge*> testPeriodic;
       for (exp3.Init (wire, TopAbs_EDGE); exp3.More(); exp3.Next())
 	{	  
 	  TopoDS_Edge edge = TopoDS::Edge (exp3.Current());
 	  int index = emap.FindIndex(edge);
 	  GEdge *e = m->edgeByTag(index);
 	  if(!e) throw;
+	  if (testPeriodic.find(e) !=testPeriodic.end()){
+	    _periodic = true;
+	    edges_taken_twice.insert(e);
+	  }
+	  else testPeriodic.insert(e);
 	  l_edges.push_back(e);
 	  l_wire.push_back(e);
 	  e->addFace(this);
 	}      
       if (l_wire.size() == 1)l_dirs.push_back(1);
-      else
+      else if (!_periodic)
 	{
 	  GVertex *last;
 	  std::list<GEdge*>::iterator it = l_wire.begin();
-	  GEdge *e1 = *it;
+	  GEdge *e1 = *it;	  
 	  ++it;
 	  GEdge *e2 = *it;
 
@@ -83,7 +89,7 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
 		}
 	      else
 		{
-		  Msg(GERROR,"Incoherent surface %d",num);
+		  Msg(GERROR,"Incoherent surface %d Edge %d (%d,%d) ",num,e->tag(),e->getBeginVertex()->tag(),e->getEndVertex()->tag());
 		}
 	    }
 	}      
@@ -179,11 +185,36 @@ SPoint2 OCCFace::parFromPoint(const SPoint3 &qp) const
     {
       Msg(GERROR,"OCC Project Point on Surface FAIL");
       return GFace::parFromPoint(qp);
-    } 
-  pnt = proj.NearestPoint();
-  double pp[2];
-  proj.LowerDistanceParameters (pp[0], pp[1]);
-  return SPoint2(pp[0],pp[1]);
+    }   
+//   if (proj.NbPoints() != 1)
+//     {
+//       Msg(WARNING,"More than one points match the projection (%d) ", proj.NbPoints());
+//     }
+
+  double U,V;
+  proj.LowerDistanceParameters (U, V);
+  //  proj.Parameters (1,U,V);
+  return SPoint2(U,V);
+}
+
+void OCCFace::parFromPoint(const SPoint3 &qp, std::list<double> &u, std::list<double> &v ) const
+{
+  gp_Pnt pnt(qp.x(),qp.y(),qp.z());
+  GeomAPI_ProjectPointOnSurf proj(pnt, occface, umin, umax, vmin, vmax);
+  if (!proj.NbPoints())
+    {
+      Msg(GERROR,"OCC Project Point on Surface FAIL");
+      GFace::parFromPoint(qp,u,v);
+      return;
+    }   
+
+  for (int i = 1;i <= proj.NbPoints();i++)
+    {
+      double U,V;
+      proj.Parameters (i,U,V);
+      u.push_back(U);
+      v.push_back(V);
+    }
 }
 
 GEntity::GeomType OCCFace::geomType() const
@@ -191,6 +222,20 @@ GEntity::GeomType OCCFace::geomType() const
   return Unknown;
 }
 
+double OCCFace::curvature (const SPoint2 &param) const
+{
+  BRepAdaptor_Surface sf(s, Standard_True);
+  BRepLProp_SLProps prop(sf, 2, 1e-5);
+  prop.SetParameters (param.x(),param.y());
+
+  if (!prop.IsCurvatureDefined())
+    {
+      return GFace::curvature (param);
+    }
+  return std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
+}
+
+
 int OCCFace::containsPoint(const SPoint3 &pt) const
 { 
   Msg(GERROR,"Not Done Yet ...");

Geo/OCCFace.h

@@ -31,7 +31,7 @@ class OCCFace : public GFace {
   TopoDS_Face s;
   Handle(Geom_Surface) occface;
   double umin, umax, vmin, vmax;
-  void buildVisTriangulation ();
+  bool _periodic;
  public:
   OCCFace(GModel *m, TopoDS_Face s, int num, TopTools_IndexedMapOfShape &emap);
 
@@ -55,13 +55,14 @@ class OCCFace : public GFace {
   virtual int geomDirection() const { return 1; }
   
   virtual bool continuous(int dim) const { return true; }
-  virtual bool periodic(int dim) const { return false; }
   virtual bool degenerate(int dim) const { return false; }
   virtual double period(int dir) const {throw;}
   ModelType getNativeType() const { return OpenCascadeModel; }
   void * getNativePtr() const { return (void*)&s; }
-  virtual bool surfPeriodic(int dim) const {throw;}
+  virtual bool surfPeriodic(int dim) const {return _periodic;}
   virtual SPoint2 parFromPoint(const SPoint3 &) const;
+  virtual void parFromPoint(const SPoint3 &, std::list<double> &u, std::list<double> &v ) const;
+  virtual double curvature(const SPoint2 &param) const;
 };
 
 #endif

Geo/OCCVertex.h

@@ -30,10 +30,12 @@
 class OCCVertex : public GVertex {
  protected:
   TopoDS_Vertex v;
-
+  mutable double max_curvature;
+  double max_curvature_of_surfaces() const;
  public:
   OCCVertex(GModel *m, int num, TopoDS_Vertex _v) : GVertex(m, num), v(_v)
   {
+    max_curvature = -1;
     mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
   }
   virtual ~OCCVertex() {}
@@ -60,8 +62,11 @@ class OCCVertex : public GVertex {
   void * getNativePtr() const { return (void*) &v; }
   virtual double prescribedMeshSizeAtVertex() const { 
     SBoundingBox3d b = model()->bounds();
-    double lc = norm ( SVector3 ( b.max() , b.min() ) );
-    return lc*CTX.mesh.lc_factor; 
+    double lc = 0.1*norm ( SVector3 ( b.max() , b.min() ) ) * CTX.mesh.lc_factor;
+    double maxc = max_curvature_of_surfaces();
+    if (maxc !=0)       
+      lc = std::min (lc,6.28/(CTX.mesh.min_circ_points*maxc));
+    return lc;
   }
 };
 

Mesh/BDS.cpp

-// $Id: BDS.cpp,v 1.62 2006-09-14 15:23:29 remacle Exp $
+// $Id: BDS.cpp,v 1.63 2006-11-15 20:46:46 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -39,11 +39,9 @@ void outputScalarField(std::list < BDS_Face * >t, const char *iii)
     BDS_Point *pts[4];
     (*tit)->getNodes(pts);
     fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
-            pts[0]->X, pts[0]->Y, pts[0]->Z,
-            pts[1]->X, pts[1]->Y, pts[1]->Z,
-            pts[2]->X, pts[2]->Y, pts[2]->Z,
-            pts[0]->radius(), pts[1]->radius(),
-            pts[2]->radius());
+            pts[0]->u, pts[0]->v, 0.0,
+            pts[1]->u, pts[1]->v, 0.0,
+            pts[2]->u, pts[2]->v, 0.0,(double)pts[0]->iD,(double)pts[1]->iD,(double)pts[2]->iD);
     ++tit;
   }
   fprintf(f, "};\n");
@@ -225,17 +223,23 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2)
 	  e = (*it);
 	  //	  if (e->p1->iD >= 0 && e->p2->iD >= 0)Msg(INFO," %d %d %g %g - %g %g",e->p1->iD,e->p2->iD,e->p1->u,e->p1->v,e->p2->u,e->p2->v);
 	  if (!e->deleted && e->p1 != p1 && e->p1 != p2 && e->p2 != p1 && e->p2 != p2)
-	    if(Intersect_Edges_2d(e->p1->X, e->p1->Y,
-				  e->p2->X, e->p2->Y,
-				  p1->X, p1->Y,
-				  p2->X, p2->Y))
+	    if(Intersect_Edges_2d(e->p1->u, e->p1->v,
+				  e->p2->u, e->p2->v,
+				  p1->u, p1->v,
+				  p2->u, p2->v))
 	      intersected.push_back(e);	  
 	  ++it;
 	}
 
-      if (!intersected.size())
+      if (!intersected.size() || ix > 10000)
 	{
-	  return find_edge (num1, num2);
+	  BDS_Edge *eee = find_edge (num1, num2);
+	  if (!eee)
+	    {
+	      outputScalarField(triangles, "debug.pos");
+	      throw;
+	    }
+	  return eee;
 	}
       
       int ichoice = ix++ % intersected.size();

Mesh/meshGEdge.cpp

-// $Id: meshGEdge.cpp,v 1.14 2006-08-26 15:13:22 remacle Exp $
+// $Id: meshGEdge.cpp,v 1.15 2006-11-15 20:46:46 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -45,33 +45,15 @@ double F_LC_ANALY (double xx, double yy, double zz)
   return 0.05 + .1*fabs(xx*yy) ;
 }
 
-double max_surf_curvature ( GPoint & gp )
-{
-  std::list<GFace *> faces = _myGEdge->faces();
-  std::list<GFace *>::iterator it =  faces.begin();
-  double curv = 0;
-  while (it != faces.end())
-    {
-      SPoint2 par = (*it)->parFromPoint(SPoint3 (gp.x(),gp.y(),gp.z()));
-      curv = std::max(curv, (*it)->curvature ( par ) );					
-      ++it;
-    }  
-  return curv;
-}
-
 double F_Lc_bis(double t)
 {
   //  const double nb_points_per_radius_of_curv = 2;
   GPoint point = _myGEdge -> point(t) ;      
   const double fact = (t-t_begin)/(t_end-t_begin);
-  double lc_here = lc_begin + fact * (lc_end-lc_begin);
+  double lc_here   = lc_begin + fact * (lc_end-lc_begin);
   SVector3 der  = _myGEdge -> firstDer(t) ;
   const double d      = norm(der);
 
-  //  double curv = max_surf_curvature ( point );
-  //  if (curv != 0)
-  //    lc_here = std::min( 1./(curv * nb_points_per_radius_of_curv),lc_here);
-
   if(CTX.mesh.bgmesh_type == ONFILE) {
     const double Lc = BGMXYZ(point.x(), point.y(), point.z());
     if(CTX.mesh.constrained_bgmesh)
@@ -172,8 +154,13 @@ void meshGEdge :: operator() (GEdge *ge)
   // to pass an extra argument... 
   _myGEdge = ge;
     
-  // first compute the length of the curve by integrating one
+
+  // compute bounds
   _myGEdgeBounds = ge->parBounds(0) ;
+  t_begin = _myGEdgeBounds.low();
+  t_end   = _myGEdgeBounds.high();
+  
+  // first compute the length of the curve by integrating one
   _myGEdgeLength = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(), 
 			       F_One_bis, Points, 1.e-4);
   List_Reset(Points);
@@ -181,9 +168,6 @@ void meshGEdge :: operator() (GEdge *ge)
   lc_begin  =  _myGEdge->getBeginVertex()->prescribedMeshSizeAtVertex();
   lc_end    =  _myGEdge->getEndVertex()->prescribedMeshSizeAtVertex();
     
-  t_begin = _myGEdgeBounds.low();
-  t_end   = _myGEdgeBounds.high();
-
   // Integrate detJ/lc du 
   double a;
   int N;

Mesh/meshGFace.cpp

-// $Id: meshGFace.cpp,v 1.17 2006-11-14 17:11:33 remacle Exp $
+// $Id: meshGFace.cpp,v 1.18 2006-11-15 20:46:46 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -751,12 +751,15 @@ void gmsh2DMeshGenerator ( GFace *gf )
 
   // mesh the domain in the parametric space -> 
   // project all points in their parametric space
+
+  Msg(DEBUG1,"Calculation of local coordinates");
   fromCartesianToParametric c2p ( gf );
   std::for_each (all_vertices.begin(),all_vertices.end(),c2p);    
 
   // compute the bounding box in parametric space
   // I do not have SBoundingBox, so I use a 3D one...
   // At the same time, number the vertices locally
+  Msg(DEBUG1,"Calculation of the bounding box");
   SBoundingBox3d bbox;
   itv = all_vertices.begin();
   int NUM = 0;
@@ -819,6 +822,7 @@ void gmsh2DMeshGenerator ( GFace *gf )
   /// At this stage the triangulation is not what we need
   ///   -) It does not necessary recover the boundaries 
   ///   -) It contains triangles outside the domain (the first edge loop is the outer one)
+  Msg(DEBUG1,"Meshing of the convex hull (%d points)",all_vertices.size());
   Make_Mesh_With_Points(&doc,doc.points,all_vertices.size()+4);
 
   // Buid a BDS_Mesh structure that is convenient for doing the actual meshing procedure
@@ -1029,7 +1033,7 @@ void meshGFace :: operator() (GFace *gf)
   Msg(STATUS2, "Meshing surface %d", gf->tag());
 
   // TEST TEST 
-  // if (gf->tag() > 5) return;
+  if (gf->surfPeriodic(2)) return;
 
   // destroy the mesh if it exists
   deMeshGFace dem;
@@ -1046,11 +1050,16 @@ void meshGFace :: operator() (GFace *gf)
   // compute loops on the fly
   // indices indicate start and end points of a loop
   // loops are not yet oriented
+  Msg(DEBUG1, "Computing edge loops");
   computeEdgeLoops(gf, points, indices);
+  Msg(DEBUG1, "Computing mean plane");
+
+
 
   // compute the mean plane, this is sometimes useful 
   gf->computeMeanPlane(points);
 
+  Msg(DEBUG1, "Generating the mesh");
   // mesh the face
   gmsh2DMeshGenerator ( gf ) ;
 
@@ -1093,7 +1102,7 @@ void orientMeshGFace::operator()(GFace *gf)
   int sign = *ori.begin();
   MEdge ref(sign > 0 ? v1 : v2, sign > 0 ? v2 : v1);
   if(shouldRevert(ref, gf->triangles) || shouldRevert(ref, gf->quadrangles)){
-    Msg(DEBUG, "Reverting orientation of mesh in face %d", gf->tag());
+    Msg(DEBUG1, "Reverting orientation of mesh in face %d", gf->tag());
     for(unsigned int i = 0; i < gf->triangles.size(); i++)
       gf->triangles[i]->revert();
     for(unsigned int i = 0; i < gf->quadrangles.size(); i++)