From f46cd50624b7a9af2bab7228f2df6f5122cfe032 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 16 Oct 2007 20:00:07 +0000
Subject: [PATCH] OCCFace::containsPoint() did not correctly take orientation
 of edge loop into account (fixes bug reported by Jacques Kools about cross in
 OCC plane faces not drawn)

---
 Geo/GFace.cpp               |  6 +++---
 Geo/OCCFace.cpp             | 31 ++++++++++++++++++++-----------
 Mesh/meshGFace.cpp          |  4 ++--
 benchmarks/2d/naca12_2d.geo |  2 +-
 4 files changed, 26 insertions(+), 17 deletions(-)

diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 4fc415a5b2..635b2e8abf 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.37 2007-09-24 08:14:29 geuzaine Exp $
+// $Id: GFace.cpp,v 1.38 2007-10-16 20:00:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -92,9 +92,9 @@ std::list<GVertex*> GFace::vertices() const
   while (it != l_edges.end()){
     GVertex *v1 = (*it)->getBeginVertex();
     GVertex *v2 = (*it)->getEndVertex();
-    if(std::find(ret.begin(), ret.end(), v1) == ret.end())
+    if(v1 && std::find(ret.begin(), ret.end(), v1) == ret.end())
       ret.push_back(v1);
-    if(std::find(ret.begin(), ret.end(), v2) == ret.end())
+    if(v2 && std::find(ret.begin(), ret.end(), v2) == ret.end())
       ret.push_back(v2);
     ++it;
   }
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 06d41e59b5..6997c634d4 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.22 2007-03-16 10:03:40 remacle Exp $
+// $Id: OCCFace.cpp,v 1.23 2007-10-16 20:00:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -70,7 +70,8 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
   _periodic[1] = surface.IsVPeriodic();
 
   ShapeAnalysis::GetFaceUVBounds(_s, umin, umax, vmin, vmax);
-  Msg(DEBUG2, "OCC Face %d with %d edges bounds (%g,%g)(%g,%g)", num, l_edges.size(),umin,umax,vmin,vmax);
+  Msg(DEBUG2, "OCC Face %d with %d edges bounds (%g,%g)(%g,%g)", 
+      num, l_edges.size(),umin,umax,vmin,vmax);
   // we do that for the projections to converge on the 
   // borders of the surface
   const double du = umax-umin;
@@ -204,23 +205,31 @@ int OCCFace::containsPoint(const SPoint3 &pt) const
 { 
   if(geomType() == Plane){
     gp_Pln pl = Handle(Geom_Plane)::DownCast(occface)->Pln();
-    
-    // OK to use the normal from the mean plane here: we compensate
-    // for the (possibly wrong) orientation at the end
-    double n[3] , c;
-    pl.Coefficients (n[0], n[1], n[2], c);
+    double n[3], c;
+    pl.Coefficients(n[0], n[1], n[2], c);
     norme(n);
     double angle = 0.;
     double v[3] = {pt.x(), pt.y(), pt.z()};
-    for(std::list<GEdge*>::const_iterator  it = l_edges.begin(); it != l_edges.end(); it++) {
+
+    std::list<int>::const_iterator ito = l_dirs.begin();
+    for(std::list<GEdge*>::const_iterator it = l_edges.begin(); it != l_edges.end(); it++){
       GEdge *c = *it;
-      int N=10;
+      int ori = 1;
+      if(ito != l_dirs.end()){
+	ori = *ito;
+	++ito;
+      }
+      int N = 10;
       Range<double> range = c->parBounds(0);
       for(int j = 0; j < N ; j++) {
 	double u1 = (double)j / (double)N;
 	double u2 = (double)(j + 1) / (double)N;
-	GPoint pp1 = c->point(range.low() + u1 * (range.high() - range.low() ));
-	GPoint pp2 = c->point(range.low() + u2 * (range.high() - range.low() ));
+	if(ori < 0){
+	  u1 = 1. - u1;
+	  u2 = 1. - u2;
+	}
+	GPoint pp1 = c->point(range.low() + u1 * (range.high() - range.low()));
+	GPoint pp2 = c->point(range.low() + u2 * (range.high() - range.low()));
 	double v1[3] = {pp1.x(), pp1.y(), pp1.z()};
 	double v2[3] = {pp2.x(), pp2.y(), pp2.z()};
 	angle += angle_plan(v, v1, v2, n);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 296dce3b0e..60cd78c686 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.98 2007-10-14 19:54:16 remacle Exp $
+// $Id: meshGFace.cpp,v 1.99 2007-10-16 20:00:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -1820,7 +1820,7 @@ void meshGFace::operator() (GFace *gf)
   // temp fix until we create MEdgeLoops in gmshFace
   Msg(DEBUG1, "Generating the mesh");
   if(gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
-    gmsh2DMeshGenerator(gf,0, true);
+    gmsh2DMeshGenerator(gf,0, false);
   }
   else{
     if(!gmsh2DMeshGeneratorPeriodic(gf,false))
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index 34085f830a..dade16bbc5 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -225,7 +225,7 @@ Plane Surface(11) = {9,10};
 
 Point(9999) = {0.6,0,0,1};
 
-Attractor Point{9999} = {0.2,0.02,lc3*10,2};
+//Attractor Point{9999} = {0.2,0.02,lc3*10,2};
 
 //Attractor Line{1,2,3,4} = {1,.05,3};
 //Mesh.Algorithm = 2;
-- 
GitLab