From 9f1bd050932c66e50e0444d1a9423061cd7cc963 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Mon, 11 Dec 2006 20:12:33 +0000
Subject: [PATCH] *** empty log message ***

---
 Geo/OCCEdge.cpp    | 31 ++++++++++++++++----------
 Mesh/BDS.cpp       | 12 +++++++---
 Mesh/meshGFace.cpp | 55 +++++++++++++++++++++++++++-------------------
 3 files changed, 61 insertions(+), 37 deletions(-)

diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index f9fcddb10a..c5136f5e37 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.14 2006-11-29 16:57:01 remacle Exp $
+// $Id: OCCEdge.cpp,v 1.15 2006-12-11 20:12:33 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -23,6 +23,8 @@
 #include "Message.h"
 #include "OCCEdge.h"
 #include "OCCFace.h"
+#include "Context.h"
+extern Context_T CTX;
 
 #if defined(HAVE_OCC)
 #include "Geom2dLProp_CLProps2d.hxx"
@@ -62,26 +64,31 @@ SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar,int dir) const
   double t0,t1;
   Handle(Geom2d_Curve) c2d = BRep_Tool::CurveOnSurface(c, *s, t0, t1);
   if (c2d.IsNull())
-    return GEdge::reparamOnFace(face, epar,dir);
+    {
+      Msg(GERROR,"Reparam on face failed : curve %d is not on surface %d\n",tag(),face->tag());
+      return GEdge::reparamOnFace(face, epar,dir);
+    }
   double u,v;
   c2d->Value(epar).Coord(u,v);
+
   if (! isSeam ( face ) )
     {
+      // sometimes OCC miserably fails ...
+      GPoint p1 = point(epar);
+      GPoint p2 = face->point(u,v);
+      const double dx = p1.x()-p2.x();
+      const double dy = p1.y()-p2.y();
+      const double dz = p1.z()-p2.z();
+      if (sqrt(dx*dx+dy*dy+dz*dz) > 1.e-7 * CTX.lc)
+	{
+	  GPoint ppp = face->closestPoint(SPoint3(p1.x(),p1.y(),p1.z()));
+	  return SPoint2(ppp.u(),ppp.v());
+	}
       return SPoint2 (u,v);
     }
   else
     {
       BRepAdaptor_Surface surface( *s );
-
-//       printf ("surface %d (%d %d) firstu %g lastu %g firstv %g lastv %g\n",
-// 	      surface.IsUPeriodic() ,
-// 	      surface.IsVPeriodic() ,
-// 	      face->tag(),
-// 	      surface.FirstUParameter(),
-// 	      surface.LastUParameter(),
-// 	      surface.FirstVParameter(),
-// 	      surface.LastVParameter());
-
       if ( surface.IsUPeriodic() )
 	{
 	  if (dir == -1) 
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 2bdac31691..588bdaf7b2 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.70 2006-12-05 14:22:05 remacle Exp $
+// $Id: BDS.cpp,v 1.71 2006-12-11 20:12:33 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -991,6 +991,9 @@ bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS
   double b[2] = {pc[0]-pa[0],pc[1]-pa[1]};
 
   double area_init = fabs(a[1] * b[2] - a[2] * b[1]);
+
+  if (area_init == 0.0) return true;
+
   double ori_init = gmsh::orient2d(pa, pb, pc);
 
   if (p == pts[0])
@@ -1006,7 +1009,7 @@ bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS
   if (area_final < 0.1 * area_init)return false;
   double ori_final = gmsh::orient2d(pa, pb, pc);
   // allow to move a point when a triangle was flat
-  return ori_init*ori_final >= 0;
+  return ori_init*ori_final > 0;
 }
 
 bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
@@ -1025,15 +1028,18 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
   while(eit != p->edges.end()) {
     if ((*eit)->numfaces() == 1) return false;
     BDS_Point *op = ((*eit)->p1 == p) ? (*eit)->p2 : (*eit)->p1;
+    //    const double l_e = (*eit)->length();     
     U += op->u; 
     V += op->v;
-    //    tot_length += (*eit)->length(); 
+    //    tot_length += l_e;
     LC += op->lc();
     ++eit;
   }
   
   U /= (p->edges.size());
   V /= (p->edges.size());
+  //  U /= tot_length;
+  //  V /= tot_length;
   LC /= p->edges.size();
 
   std::list < BDS_Face * >ts;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 5b3e4dad40..5c57a060ac 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.38 2006-12-05 14:22:05 remacle Exp $
+// $Id: meshGFace.cpp,v 1.39 2006-12-11 20:12:33 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -427,7 +427,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
 	}
 	// recompute mesh sizes takin into account curvature , BGMESH & co
 	m.cleanup();  
-	if (IT % 5 == 0 && (CTX.mesh.lc_from_curvature || BGMExists()))
+	if (IT == 5 && (CTX.mesh.lc_from_curvature || BGMExists()))
 	{
 	  std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
 	  while (itp != m.points.end())
@@ -441,7 +441,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
 		}
 	      ++itp;
 	    }
-	  for (int ITERA = 0;ITERA< 6; ITERA++);
+	  for (int ITERA = 0;ITERA< 0; ITERA++);
 	  {
 	    it = m.edges.begin();
 	    while (it != m.edges.end())
@@ -556,7 +556,7 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge)
 // domain, including embedded points 
 // and surfaces
 
-void gmsh2DMeshGenerator ( GFace *gf )
+bool gmsh2DMeshGenerator ( GFace *gf )
 {
 
   typedef std::set<MVertex*> v_container ;
@@ -570,6 +570,8 @@ void gmsh2DMeshGenerator ( GFace *gf )
   it = edges.begin();
   while(it != edges.end())
     {
+      if ((*it)->isSeam(gf))return false;
+
       all_vertices.insert ( (*it)->mesh_vertices.begin() , (*it)->mesh_vertices.end() );
       all_vertices.insert ( (*it)->getBeginVertex()->mesh_vertices.begin() , (*it)->getBeginVertex()->mesh_vertices.end() );
       all_vertices.insert ( (*it)->getEndVertex()->mesh_vertices.begin() , (*it)->getEndVertex()->mesh_vertices.end() );
@@ -720,7 +722,7 @@ void gmsh2DMeshGenerator ( GFace *gf )
       if (!recover_medge ( m, *it))
 	{
 	  Msg(WARNING,"Face not meshed");
-	  return;
+	  return false;
 	}
       ++it;
     }
@@ -796,8 +798,8 @@ void gmsh2DMeshGenerator ( GFace *gf )
    char name[245];
    //sprintf(name,"param%d.pos",gf->tag());
    //outputScalarField(m->triangles, name,1);
-//   sprintf(name,"real%d.pos",gf->tag());
-//   outputScalarField(m->triangles, name,0);
+   //   sprintf(name,"real%d.pos",gf->tag());
+   //   outputScalarField(m->triangles, name,0);
 
 
   m->cleanup();
@@ -863,14 +865,13 @@ void gmsh2DMeshGenerator ( GFace *gf )
 
   // delete the mesh
 
-  //  char name[245];
-  //  sprintf(name,"s%d.pos",gf->tag());
-  //  outputScalarField(m->triangles, name);
   delete m;
 
   delete [] U_;
   delete [] V_;
 
+  return true;
+
 }
 
 // this function buils a list of vertices (BDS) that 
@@ -1082,10 +1083,11 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
      result.insert(result.end(),edgeLoop_BDS.begin(),edgeLoop_BDS.end());	         
    }
 
-//    for (int i=0;i<result.size();i++)
-//      {
-//        printf("point %3d (%8.5f %8.5f) (%2d,%2d)\n",i,result[i]->u,result[i]->v,result[i]->g->classif_tag,result[i]->g->classif_degree);
-//      }
+  if (gf->tag() == 280)
+    for (int i=0;i<result.size();i++)
+      {
+        printf("point %3d (%8.5f %8.5f) (%2d,%2d)\n",i,result[i]->u,result[i]->v,result[i]->g->classif_tag,result[i]->g->classif_degree);
+      }
 
   return true;
 }
@@ -1093,6 +1095,8 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
 
 bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
 {
+  //  if(gf->tag() != 6)return true;
+
 
   std::map<BDS_Point*,MVertex*> recover_map;
 
@@ -1108,8 +1112,8 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
 
   // Buid a BDS_Mesh structure that is convenient for doing the actual meshing procedure    
   BDS_Mesh *m = new BDS_Mesh;
-  m->scalingU = fabs(du);
-  m->scalingV = fabs(dv);
+  //  m->scalingU = fabs(du);
+  //  m->scalingV = fabs(dv);
   std::vector< std::vector<BDS_Point* > > edgeLoops_BDS;
   SBoundingBox3d bbox;
   int nbPointsTotal = 0;
@@ -1186,6 +1190,8 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
   // Free stuff
   free (doc.points);
   free (doc.delaunay);
+
+
  
   // Recover the boundary edges
   // and compute characteristic lenghts using mesh edge spacing
@@ -1281,6 +1287,8 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
   m->del_point(m->find_point(-3));
   m->del_point(m->find_point(-4));
     
+
+
   // goto hhh;
   // start mesh generation
   
@@ -1337,12 +1345,12 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
   
   
   // delete the mesh
+    char name[245];
+   sprintf(name,"param%d.pos",gf->tag());
+    outputScalarField(m->triangles, name,1);
+    sprintf(name,"real%d.pos",gf->tag());
+    outputScalarField(m->triangles, name,0);
   
-//   char name[245];
-//   sprintf(name,"param%d.pos",gf->tag());
-//   outputScalarField(m->triangles, name,1);
-//   sprintf(name,"real%d.pos",gf->tag());
-//   outputScalarField(m->triangles, name,0);
   delete m; 
   return true;
 }
@@ -1385,7 +1393,10 @@ void meshGFace :: operator() (GFace *gf)
   if(gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty())
     gmsh2DMeshGenerator ( gf ) ;
   else
-    gmsh2DMeshGeneratorPeriodic ( gf ) ;
+    {
+      if (!gmsh2DMeshGeneratorPeriodic ( gf ))
+	Msg(GERROR, "Impossible to mesh face %d",gf->tag());
+    }
 
   Msg(DEBUG1, "type %d %d triangles generated, %d internal vertices",
       gf->geomType(),gf->triangles.size(),gf->mesh_vertices.size());
-- 
GitLab