22711: [CEA 1297] Regression of mesh with a projection on a spherical face

[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx

diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx
index 16f555477d17f08a72ef69ae07dfafb911bc7214..31f2d88b7b98a25c5d8339de5f334be10b009709 100644 (file)
--- a/src/SMESH/SMESH_MesherHelper.cxx
+++ b/src/SMESH/SMESH_MesherHelper.cxx
@@ -242,8 +242,6 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
   for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
   {
     const TopoDS_Face& face = TopoDS::Face( eF.Current() );
-    // TopLoc_Location loc;
-    // Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
 
     // if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
     //      surface->IsUClosed()   || surface->IsVClosed() )
@@ -382,10 +380,13 @@ void SMESH_MesherHelper::AddTLinkNode(const SMDS_MeshNode* n1,
  */
 //================================================================================
 
-void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
+bool SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
 {
-  if ( edge->IsQuadratic() )
+  if ( edge && edge->IsQuadratic() )
     AddTLinkNode(edge->GetNode(0), edge->GetNode(1), edge->GetNode(2));
+  else
+    return false;
+  return true;
 }
 
 //================================================================================
@@ -394,8 +395,9 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
  */
 //================================================================================
 
-void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
+bool SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
 {
+  bool isQuad = true;
   if ( !f->IsPoly() )
     switch ( f->NbNodes() ) {
     case 7:
@@ -419,7 +421,9 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
       AddTLinkNode(f->GetNode(2),f->GetNode(3),f->GetNode(6));
       AddTLinkNode(f->GetNode(3),f->GetNode(0),f->GetNode(7)); break;
     default:;
+      isQuad = false;
     }
+  return isQuad;
 }
 
 //================================================================================
@@ -428,7 +432,7 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
  */
 //================================================================================
 
-void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
+bool SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
 {
   if ( volume->IsQuadratic() )
   {
@@ -462,7 +466,9 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
                          nFCenter ));
       }
     }
+    return true;
   }
+  return false;
 }
 
 //================================================================================
@@ -2330,11 +2336,12 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
           while ( ++u_n != sortedBaseNN.end() && !isNodeInSubMesh( u_n->second, faceSubMesh ));
           sortedBaseNN.erase( sortedBaseNN.begin(), u_n );
         }
-        if ( u_n = --sortedBaseNN.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
-        {
-          while ( u_n != sortedBaseNN.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
-          sortedBaseNN.erase( ++u_n, sortedBaseNN.end() );
-        }
+        if ( !sortedBaseNN.empty() )
+          if ( u_n = --sortedBaseNN.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+          {
+            while ( u_n != sortedBaseNN.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
+            sortedBaseNN.erase( ++u_n, sortedBaseNN.end() );
+          }
         if ( sortedBaseNN.empty() ) continue;
       }
 
@@ -2351,7 +2358,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
         u2nn->second.push_back( u_n->second );
       }
     }
-    if ( theParam2ColumnMap.empty() )
+    if ( theParam2ColumnMap.size() < 2 )
       return false;
   }
 
@@ -2734,24 +2741,22 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
  */
 //================================================================================
 
-double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
-                                     const TopoDS_Edge & theE2,
-                                     const TopoDS_Face & theFace,
-                                     gp_Vec*             theFaceNormal)
+double SMESH_MesherHelper::GetAngle( const TopoDS_Edge &   theE1,
+                                     const TopoDS_Edge &   theE2,
+                                     const TopoDS_Face &   theFace,
+                                     const TopoDS_Vertex & theCommonV,
+                                     gp_Vec*               theFaceNormal)
 {
   double angle = 1e100;
   try
   {
-    TopoDS_Vertex vCommon;
-    if ( !TopExp::CommonVertex( theE1, theE2, vCommon ))
-      return angle;
     double f,l;
     Handle(Geom_Curve)     c1 = BRep_Tool::Curve( theE1, f,l );
     Handle(Geom_Curve)     c2 = BRep_Tool::Curve( theE2, f,l );
     Handle(Geom2d_Curve) c2d1 = BRep_Tool::CurveOnSurface( theE1, theFace, f,l );
     Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace );
-    double                 p1 = BRep_Tool::Parameter( vCommon, theE1 );
-    double                 p2 = BRep_Tool::Parameter( vCommon, theE2 );
+    double                 p1 = BRep_Tool::Parameter( theCommonV, theE1 );
+    double                 p2 = BRep_Tool::Parameter( theCommonV, theE2 );
     if ( c1.IsNull() || c2.IsNull() )
       return angle;
     gp_XY uv = c2d1->Value( p1 ).XY();
@@ -2760,10 +2765,10 @@ double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
     gp_Vec vec1, vec2, vecRef = du ^ dv;
     int  nbLoops = 0;
     double p1tmp = p1;
-    while ( vecRef.SquareMagnitude() < std::numeric_limits<double>::min() )
+    while ( vecRef.SquareMagnitude() < 1e-25 )
     {
       double dp = ( l - f ) / 1000.;
-      p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? +1. : -1.);
+      p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? -1. : +1.);
       uv = c2d1->Value( p1tmp ).XY();
       surf->D1( uv.X(), uv.Y(), p, du, dv );
       vecRef = du ^ dv;
@@ -4214,6 +4219,8 @@ namespace { // Structures used by FixQuadraticElements()
                 nInSolid = n;
               else if ( subIDs.count( n->getshapeId() ))
                 nOnFace.push_back( n );
+              else
+                nInSolid = n;
             }
             if ( !nInSolid || nOnFace.size() != nbN - 1 )
               continue;