*** empty log message ***

Geo/GEdge.cpp

-// $Id: GEdge.cpp,v 1.17 2006-11-15 20:46:46 remacle Exp $
+// $Id: GEdge.cpp,v 1.18 2006-11-16 21:14:10 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -22,6 +22,7 @@
 #include <algorithm>
 #include "GModel.h"
 #include "GEdge.h"
+#include "GFace.h"
 #include "GmshDefines.h"
 
 GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
@@ -120,6 +121,15 @@ SVector3 GEdge::secondDer(double par) const
   return 500*(x2-x1);
 }
 
+// Reparmaterize the point onto the given face.
+SPoint2 GEdge::reparamOnFace(GFace *face, double epar,int dir) const
+{
+  const GPoint p3 = point (epar);
+  SPoint3 sp3(p3.x(),p3.y(),p3.z());
+  return face->parFromPoint(sp3);
+}
+
+
 double GEdge::curvature(double par) const 
 {
   double eps1 = 1.e-3;

Geo/GEdge.h

@@ -76,7 +76,7 @@ class GEdge : public GEntity {
   virtual double curvature (double par) const;  
 
   // Reparmaterize the point onto the given face.
-  virtual SPoint2 reparamOnFace(GFace *face, double epar,int dir) const = 0;
+  virtual SPoint2 reparamOnFace(GFace *face, double epar,int dir) const ;
 
   // Recompute the mesh partitions defined on this edge.
   void recomputeMeshPartitions();

Geo/GModelIO_OCC.cpp

-  // $Id: GModelIO_OCC.cpp,v 1.7 2006-11-16 18:48:00 geuzaine Exp $
+  // $Id: GModelIO_OCC.cpp,v 1.8 2006-11-16 21:14:10 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -67,6 +67,12 @@ public:
 void OCC_Internals :: buildLists ()
 {
   TopExp_Explorer exp0, exp1, exp2, exp3, exp4, exp5;
+  somap.Clear();
+  shmap.Clear();
+  fmap.Clear();
+  wmap.Clear();
+  emap.Clear();
+  vmap.Clear();
   
   for (exp0.Init(shape, TopAbs_SOLID);
        exp0.More(); exp0.Next())
@@ -485,7 +491,7 @@ void OCC_Internals :: loadBREP (const char *fn)
   Standard_Boolean result = BRepTools::Read( shape, (char*)fn, aBuilder );
   BRepTools::Clean (shape);
   buildLists();
-  //  HealGeometry();
+  HealGeometry();
   BRepTools::Clean (shape);
 }
 
@@ -497,8 +503,12 @@ void OCC_Internals :: loadSTEP (const char *fn)
   reader.TransferRoots (); 
   shape = reader.OneShape();  
   BRepTools::Clean (shape);
+  buildLists();
+  HealGeometry();
+  BRepTools::Clean (shape);
 }
 
+
 void OCC_Internals :: loadIGES (const char *fn)
 {
   IGESControl_Reader reader;
@@ -507,6 +517,9 @@ void OCC_Internals :: loadIGES (const char *fn)
   reader.TransferRoots (); 
   shape = reader.OneShape();  
   BRepTools::Clean (shape);
+  buildLists();
+  HealGeometry();
+  BRepTools::Clean (shape);
 }
 
 void OCC_Internals :: buildGModel (GModel *model)

Geo/GVertex.cpp

-// $Id: GVertex.cpp,v 1.6 2006-10-12 01:35:32 geuzaine Exp $
+// $Id: GVertex.cpp,v 1.7 2006-11-16 21:14:10 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -20,6 +20,7 @@
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "GVertex.h"
+#include "GFace.h"
 #include <algorithm>
 
 GVertex::GVertex(GModel *m, int tag) : GEntity (m, tag) 
@@ -43,6 +44,13 @@ void GVertex::delEdge(GEdge *e)
   l_edges.erase(std::find(l_edges.begin(), l_edges.end(), e));
 }
 
+
+SPoint2 GVertex::reparamOnFace ( GFace *gf ) const
+{
+  return gf->parFromPoint ( SPoint3(x(),y(),z() ));
+}
+
+
 std::string GVertex::getAdditionalInfoString()
 {
   char str[256];

Geo/GVertex.h

@@ -23,6 +23,7 @@
 #include "GEntity.h"
 #include "MVertex.h"
 #include "GPoint.h"
+#include "SPoint2.h"
 
 // A model vertex
 class GVertex : public GEntity 
@@ -43,6 +44,7 @@ class GVertex : public GEntity
   virtual GeomType geomType() const {return Point;}
   virtual double prescribedMeshSizeAtVertex() const {return 0;}
   virtual SBoundingBox3d bounds(){ return SBoundingBox3d(SPoint3(x(), y(), z())); }
+  virtual SPoint2 reparamOnFace ( GFace *gf ) const;
   virtual std::string getAdditionalInfoString();
 };
 

Geo/OCCEdge.cpp

-// $Id: OCCEdge.cpp,v 1.6 2006-11-16 18:48:00 geuzaine Exp $
+// $Id: OCCEdge.cpp,v 1.7 2006-11-16 21:14:10 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -22,26 +22,49 @@
 #include "GModel.h"
 #include "Message.h"
 #include "OCCEdge.h"
+#include "OCCFace.h"
 
 #if defined(HAVE_OCC)
 
 OCCEdge::OCCEdge(GModel *model, TopoDS_Edge edge, int num, GVertex *v1, GVertex *v2)
-  : GEdge(model, num, v1, v2), c(edge)
+  : GEdge(model, num, v1, v2), trimmed(0),c(edge)
 {
   curve = BRep_Tool::Curve(c, s0, s1);
-  if (curve.IsNull())
+}
+
+Range<double> OCCEdge::parBounds(int i) const
+{ 
+  //  double a,b;
+  //  BRep_Tool::Range (c,a,b); 
+  return(Range<double>(s0,s1));
+}
+
+void OCCEdge::setTrimmed (OCCFace *f)
+{
+  if (!trimmed)
     {
-      Msg(WARNING,"OCC Curve %d is not a 3D curve",tag());
+      trimmed = f;
+      const TopoDS_Face *s = (TopoDS_Face*) trimmed->getNativePtr();
+      curve2d = BRep_Tool::CurveOnSurface(c, *s, s0, s1);
+      if (curve2d.IsNull())  trimmed = 0;
     }
 }
 
-Range<double> OCCEdge::parBounds(int i) const
+SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar,int dir) const
 {
-  double a,b;
-  BRep_Tool::Range (c,a,b); 
-  return(Range<double>(a,b));
+  double t0,t1;
+  const TopoDS_Face *s = (TopoDS_Face*) face->getNativePtr();
+  Handle(Geom2d_Curve) c2d = BRep_Tool::CurveOnSurface(c, *s, t0, t1);
+  if (c2d.IsNull())
+    return GEdge::reparamOnFace(face, epar,dir);
+
+  double u,v;
+  c2d->Value(epar).Coord(u,v);
+  return SPoint2 (u,v);
 }
 
+
+
 GPoint OCCEdge::point(double par) const
 {
   double s0,s1;  
@@ -50,8 +73,15 @@ GPoint OCCEdge::point(double par) const
       gp_Pnt pnt = curve->Value (par);
       return GPoint(pnt.X(),pnt.Y(),pnt.Z());
     }
+  else if (trimmed)
+    {
+      double u,v;
+      curve2d->Value(par).Coord(u,v);
+      return trimmed->point(u,v);
+    }
   else
     {
+      Msg(WARNING,"OCC Curve %d is neither a 3D curve not a trimmed curve",tag());
       return GPoint (0,0,0);
     }
 }

Geo/OCCEdge.h

@@ -26,6 +26,8 @@
 #include "Mesh.h"
 #include "Range.h"
 
+class OCCFace;
+
 #if defined(HAVE_OCC)
 
 class OCCEdge : public GEdge {
@@ -33,8 +35,8 @@ class OCCEdge : public GEdge {
   TopoDS_Edge c;
   double s0,s1;
   Handle(Geom_Curve) curve;
-  Handle(Geom2d_Curve) curve2d;
-
+  mutable Handle(Geom2d_Curve) curve2d;
+  mutable GFace *trimmed;
  public:
   OCCEdge(GModel *model, TopoDS_Edge _e, int num, GVertex *v1, GVertex *v2);
   virtual ~OCCEdge() {}
@@ -50,13 +52,14 @@ class OCCEdge : public GEdge {
   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; }
+  virtual SPoint2 reparamOnFace(GFace * face, double epar, int dir) const ;
   ModelType getNativeType() const { return OpenCascadeModel; }
   void * getNativePtr() const { return (void*) &c; }
   virtual double parFromPoint(const SPoint3 &pt) const;
   virtual int minimumMeshSegments () const;
   virtual int minimumDrawSegments () const;
   bool is3D() const {return !curve.IsNull();}
+  void setTrimmed (OCCFace *);
 };
 
 #endif

Geo/OCCFace.cpp

-// $Id: OCCFace.cpp,v 1.7 2006-11-16 18:48:00 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.8 2006-11-16 21:14:10 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -45,6 +45,11 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
 	  if(!e) throw;
 	  l_wire.push_back(e);
 	  e->addFace(this);
+	  if (!e->is3D())
+	    {
+	      OCCEdge *occe = (OCCEdge*)e;
+	      occe->setTrimmed(this);
+	    }
 	}      
 
       GEdgeLoop el (l_wire);

Geo/OCCVertex.cpp

-// $Id: OCCVertex.cpp,v 1.2 2006-11-16 18:48:00 geuzaine Exp $
+// $Id: OCCVertex.cpp,v 1.3 2006-11-16 21:14:10 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -26,16 +26,16 @@
 
 #if defined(HAVE_OCC)
 
-double max_surf_curvature ( double x, double y, double z , const GEdge *_myGEdge)
+double max_surf_curvature ( const GVertex *gv, double x, double y, double z , const GEdge *_myGEdge)
 {
   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 (x,y,z));
+      SPoint2 par = gv->reparamOnFace((*it));
       double cc = (*it)->curvature ( par );
-      if (cc < 1.e4) curv = std::max(curv, cc );					
+      if (cc < 1.e2) curv = std::max(curv, cc );					
       ++it;
     }  
   return curv;
@@ -47,7 +47,7 @@ double OCCVertex::max_curvature_of_surfaces() const
     {
       for (std::list<GEdge*> :: const_iterator it = l_edges.begin() ; it != l_edges.end() ; ++it )
 	{
-	  max_curvature = std::max ( max_surf_curvature (x(), y(), z(), *it) , max_curvature);  
+	  max_curvature = std::max ( max_surf_curvature (this, x(), y(), z(), *it) , max_curvature);  
 	}
       //      printf("max curvature (%d) = %12.5E lc = %12.5E\n",tag(),max_curvature,prescribedMeshSizeAtVertex());
 
@@ -55,4 +55,26 @@ double OCCVertex::max_curvature_of_surfaces() const
   return max_curvature;
 }
 
+SPoint2 OCCVertex::reparamOnFace ( GFace *gf ) const
+{
+  std::list<GEdge*>::const_iterator it = l_edges.begin();
+  while (it != l_edges.end())
+    {
+      std::list<GEdge*> l_edges = gf->edges();
+      
+      if (std::find(l_edges.begin(),l_edges.end(),*it) != l_edges.end())
+	{
+	  const TopoDS_Face *s = (TopoDS_Face*) gf->getNativePtr();
+	  const TopoDS_Edge *c = (TopoDS_Edge*) (*it)->getNativePtr();
+	  double s1,s0;
+	  Handle(Geom2d_Curve) curve2d = BRep_Tool::CurveOnSurface(*c, *s, s0, s1);
+	  if ((*it)->getBeginVertex() == this)
+	    return (*it)->reparamOnFace(gf,s0,1);
+	  else if ((*it)->getEndVertex() == this)
+	    return (*it)->reparamOnFace(gf,s1,1);
+	}
+      ++it;
+    }
+}
+
 #endif

Geo/OCCVertex.h

@@ -69,6 +69,7 @@ class OCCVertex : public GVertex {
       lc = std::min (lc,6.28/(CTX.mesh.min_circ_points*maxc));
     return lc;
   }
+  virtual SPoint2 reparamOnFace ( GFace *gf ) const;
 };
 
 #endif

Geo/gmshEdge.h

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

Mesh/BDS.cpp

-// $Id: BDS.cpp,v 1.63 2006-11-15 20:46:46 remacle Exp $
+// $Id: BDS.cpp,v 1.64 2006-11-16 21:14:10 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -237,7 +237,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2)
 	  if (!eee)
 	    {
 	      outputScalarField(triangles, "debug.pos");
-	      throw;
+	      return 0;
 	    }
 	  return eee;
 	}

Mesh/meshGFace.cpp

-// $Id: meshGFace.cpp,v 1.18 2006-11-15 20:46:46 remacle Exp $
+// $Id: meshGFace.cpp,v 1.19 2006-11-16 21:14:10 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -95,12 +95,41 @@ public :
   fromCartesianToParametric ( GFace *_gf )  
     : gf(_gf){}
   void operator () (MVertex * v)
+  {
+
+    GEntity *ge = v->onWhat();
+
+    // here, the point is classified on a model edge. So
+    // it is possible that the CAD can easily compute
+    // parametric coordinates of the point on the face using the
+    // parametric coordinate of the point on the edge. By default,
+    // the model will use parFromPoint
+    if (ge->dim() == 0)
+      {
+	GVertex *gve = (GVertex*)ge;
+	SPoint2 param = gve->reparamOnFace(gf);
+	v->x() = param.x();  
+	v->y() = param.y();
+	v->z() = 0.0;
+      }
+    else if (ge->dim() == 1)
+      {
+	GEdge *ged = (GEdge*)ge;
+	double u;
+	v->getParameter(0,u);
+	SPoint2 param = ged->reparamOnFace(gf,u,1);
+	v->x() = param.x();  
+	v->y() = param.y();
+	v->z() = 0.0;
+      }
+    else
       {
 	SPoint2 param =  gf->parFromPoint (SPoint3(v->x(),v->y(),v->z()));
 	v->x() = param.x();  
 	v->y() = param.y();
 	v->z() = 0.0;
       }
+  }
 };
 
 fromParametricToCartesian::fromParametricToCartesian ( GFace *_gf )  
@@ -656,7 +685,7 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge)
 
   BDS_Edge * e = m->recover_edge ( vstart->getNum(), vend->getNum());
   if (e)e->g = g;
-  else throw;
+  else return false;
 
   for (unsigned int i=1;i<ge->mesh_vertices.size();i++)
     {
@@ -854,7 +883,11 @@ void gmsh2DMeshGenerator ( GFace *gf )
   it = edges.begin();
   while(it != edges.end())
     {
-      recover_medge ( m, *it);
+      if (!recover_medge ( m, *it))
+	{
+	  Msg(WARNING,"Face not meshed");
+	  return;
+	}
       ++it;
     }
   //  Msg(INFO,"Boundary Edges recovered for surface %d",gf->tag());