From abe0df8535d1c10e54f401f96d745975f02f394b Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Thu, 16 Nov 2006 21:14:10 +0000
Subject: [PATCH] *** empty log message ***

---
 Geo/GEdge.cpp        | 12 ++++++++++-
 Geo/GEdge.h          |  2 +-
 Geo/GModelIO_OCC.cpp | 17 +++++++++++++--
 Geo/GVertex.cpp      | 10 ++++++++-
 Geo/GVertex.h        |  2 ++
 Geo/OCCEdge.cpp      | 50 +++++++++++++++++++++++++++++++++++---------
 Geo/OCCEdge.h        |  9 +++++---
 Geo/OCCFace.cpp      |  7 ++++++-
 Geo/OCCVertex.cpp    | 34 ++++++++++++++++++++++++------
 Geo/OCCVertex.h      |  1 +
 Geo/gmshEdge.h       |  1 -
 Mesh/BDS.cpp         |  4 ++--
 Mesh/meshGFace.cpp   | 47 ++++++++++++++++++++++++++++++++++-------
 13 files changed, 161 insertions(+), 35 deletions(-)

diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 0c138c10d2..de27fb290b 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -1,4 +1,4 @@
-// $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;
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 905052aea7..8ab4a9d140 100644
--- a/Geo/GEdge.h
+++ b/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();
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 3895d67441..8308ab9e1d 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -1,4 +1,4 @@
-  // $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)
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index a8574f9097..cf847ffb91 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -1,4 +1,4 @@
-// $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];
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
index f93b341686..91b7a5c069 100644
--- a/Geo/GVertex.h
+++ b/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();
 };
 
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 4cd3930bbe..22724d5af3 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $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())
-    {
-      Msg(WARNING,"OCC Curve %d is not a 3D curve",tag());
-    }
 }
 
 Range<double> OCCEdge::parBounds(int i) const
 { 
-  double a,b;
-  BRep_Tool::Range (c,a,b); 
-  return(Range<double>(a,b));
+  //  double a,b;
+  //  BRep_Tool::Range (c,a,b); 
+  return(Range<double>(s0,s1));
+}
+
+void OCCEdge::setTrimmed (OCCFace *f)
+{
+  if (!trimmed)
+    {
+      trimmed = f;
+      const TopoDS_Face *s = (TopoDS_Face*) trimmed->getNativePtr();
+      curve2d = BRep_Tool::CurveOnSurface(c, *s, s0, s1);
+      if (curve2d.IsNull())  trimmed = 0;
+    }
 }
 
+SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar,int dir) const
+{
+  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,9 +73,16 @@ 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
     {
-      return GPoint(0,0,0);
+      Msg(WARNING,"OCC Curve %d is neither a 3D curve not a trimmed curve",tag());
+      return GPoint (0,0,0);
     }
 }
 
diff --git a/Geo/OCCEdge.h b/Geo/OCCEdge.h
index a9e0f17374..6a045ad6c1 100644
--- a/Geo/OCCEdge.h
+++ b/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
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index c339ea89f5..da730f051b 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $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);
diff --git a/Geo/OCCVertex.cpp b/Geo/OCCVertex.cpp
index 36fe163d84..8ac63a82de 100644
--- a/Geo/OCCVertex.cpp
+++ b/Geo/OCCVertex.cpp
@@ -1,4 +1,4 @@
-// $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,12 +47,34 @@ 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());
 
     }
   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
diff --git a/Geo/OCCVertex.h b/Geo/OCCVertex.h
index f9feb664b4..6fbd53d906 100644
--- a/Geo/OCCVertex.h
+++ b/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
diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h
index 10eb2ce711..bcfcb0d7e1 100644
--- a/Geo/gmshEdge.h
+++ b/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;
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 413f5d462e..323ad8afc0 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $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;
 	}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index c928b98bb8..0e04a234b8 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $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
 //
@@ -96,10 +96,39 @@ public :
     : gf(_gf){}
   void operator () (MVertex * v)
   {
-    SPoint2 param =  gf->parFromPoint (SPoint3(v->x(),v->y(),v->z()));
-    v->x() = param.x();  
-    v->y() = param.y();
-    v->z() = 0.0;
+
+    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;
+      }
   }
 };
 
@@ -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());
-- 
GitLab