Salome HOME
22525: EDF 2863 SMESH : problem between BLSurf with viscous layer and projection
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index 65a19e286c96bad7e4e20d3433964b09c56187f9..6a97d47246cc5928909c9a7037536ebf849d245b 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -6,7 +6,7 @@
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -245,8 +245,8 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
     TopLoc_Location loc;
     Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
 
-    if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
-         surface->IsUClosed()   || surface->IsVClosed() )
+    // if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
+    //      surface->IsUClosed()   || surface->IsVClosed() )
     {
       //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
       //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface();
@@ -281,7 +281,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
         }
 
         // look for a degenerated edge
-        if ( BRep_Tool::Degenerated( edge )) {
+        if ( SMESH_Algo::isDegenerated( edge )) {
           myDegenShapeIds.insert( meshDS->ShapeToIndex( edge ));
           for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
             myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
@@ -478,7 +478,10 @@ bool SMESH_MesherHelper::toCheckPosOnShape(int shapeID ) const
 
 void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const
 {
-  ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok));
+  std::map< int,bool >::iterator sh_ok = 
+    ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok)).first;
+  if ( !ok )
+    sh_ok->second = ok;
 }
 
 //=======================================================================
@@ -600,10 +603,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
         for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() )
           uvOK = ( V == vert.Current() );
         if ( !uvOK ) {
-#ifdef _DEBUG_
           MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
-               << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
-#endif
+                    << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
           // get UV of a vertex closest to the node
           double dist = 1e100;
           gp_Pnt pn = XYZ( n );
@@ -663,9 +664,10 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
                                      const bool           force,
                                      double               distXYZ[4]) const
 {
-  int shapeID = n->getshapeId();
+  int  shapeID = n->getshapeId();
   bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() ));
-  if ( force || toCheckPosOnShape( shapeID ) || infinit )
+  bool zero    = ( uv.X() == 0. && uv.Y() == 0. );
+  if ( force || toCheckPosOnShape( shapeID ) || infinit || zero )
   {
     // check that uv is correct
     TopLoc_Location loc;
@@ -898,8 +900,10 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
                                     const bool           force,
                                     double               distXYZ[4]) const
 {
-  int shapeID = n->getshapeId();
-  if ( force || toCheckPosOnShape( shapeID ))
+  int  shapeID = n->getshapeId();
+  bool infinit = Precision::IsInfinite( u );
+  bool zero    = ( u == 0. );
+  if ( force || toCheckPosOnShape( shapeID ) || infinit || zero )
   {
     TopLoc_Location loc; double f,l;
     Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l );
@@ -915,12 +919,17 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
     {
       gp_Pnt nodePnt = SMESH_TNodeXYZ( n );
       if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
-      gp_Pnt curvPnt = curve->Value( u );
-      double dist = nodePnt.Distance( curvPnt );
-      if ( distXYZ ) {
-        curvPnt.Transform( loc );
-        distXYZ[0] = dist;
-        distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
+      gp_Pnt curvPnt;
+      double dist = u;
+      if ( !infinit )
+      {
+        curvPnt = curve->Value( u );
+        dist    = nodePnt.Distance( curvPnt );
+        if ( distXYZ ) {
+          curvPnt.Transform( loc );
+          distXYZ[0] = dist;
+          distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
+        }
       }
       if ( dist > tol )
       {
@@ -944,6 +953,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
         }
         Quantity_Parameter U = projector->LowerDistanceParameter();
         u = double( U );
+        MESSAGE(" f " << f << " l " << l << " u " << u);
         curvPnt = curve->Value( u );
         dist = nodePnt.Distance( curvPnt );
         if ( distXYZ ) {
@@ -968,6 +978,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
       }
       if (( u < f-tol || u > l+tol ) && force )
       {
+        MESSAGE("u < f-tol || u > l+tol  ; u " << u << " f " << f << " l " << l);
         // node is on vertex but is set on periodic but trimmed edge (issue 0020890)
         try
         {
@@ -1172,11 +1183,12 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
 
   gp_XY  uvAvg;
   gp_Pnt P;
+  bool toCheck = true;
   if ( !F.IsNull() && !force3d )
   {
     uvAvg = calcTFI (0.5, 0.5,
-                     GetNodeUV(F,n1,n3),  GetNodeUV(F,n2,n4),
-                     GetNodeUV(F,n3,n1),  GetNodeUV(F,n4,n2), 
+                     GetNodeUV(F,n1,n3,&toCheck), GetNodeUV(F,n2,n4,&toCheck),
+                     GetNodeUV(F,n3,n1,&toCheck), GetNodeUV(F,n4,n2,&toCheck), 
                      GetNodeUV(F,n12,n3), GetNodeUV(F,n23,n4),
                      GetNodeUV(F,n34,n2), GetNodeUV(F,n41,n2));
     TopLoc_Location loc;
@@ -1196,10 +1208,10 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
 
     if ( !F.IsNull() ) // force3d
     {
-      uvAvg = (GetNodeUV(F,n1,n3) +
-               GetNodeUV(F,n2,n4) +
-               GetNodeUV(F,n3,n1) +
-               GetNodeUV(F,n4,n2)) / 4;
+      uvAvg = (GetNodeUV(F,n1,n3,&toCheck) +
+               GetNodeUV(F,n2,n4,&toCheck) +
+               GetNodeUV(F,n3,n1,&toCheck) +
+               GetNodeUV(F,n4,n2,&toCheck)) / 4;
       //CheckNodeUV( F, centralNode, uvAvg, 2*BRep_Tool::Tolerance( F ), /*force=*/true);
       meshDS->SetNodeOnFace( centralNode, faceID, uvAvg.X(), uvAvg.Y() );
     }
@@ -1416,8 +1428,15 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
       return getMediumNodeOnComposedWire(n1,n2,force3d);
     }
     E = TopoDS::Edge(meshDS->IndexToShape( edgeID = pos.first ));
-    u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]);
-    u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]);
+    try {
+      u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]);
+      u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]);
+    }
+    catch ( Standard_Failure& f )
+    {
+      // issue 22502 / a node is on VERTEX not belonging to E
+      return getMediumNodeOnComposedWire(n1,n2,force3d);
+    }
   }
 
   if ( !force3d & uvOK[0] && uvOK[1] )
@@ -2489,9 +2508,16 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace)
     {
       SMESH_TNodeXYZ nPnt[3];
       SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator();
+      int iNodeOnFace = 0, iPosDim = SMDS_TOP_VERTEX;
       for ( int iN = 0; nodesIt->more() && iN < 3; ++iN) // loop on nodes
+      {
         nPnt[ iN ] = nodesIt->next();
-
+        if ( nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition() > iPosDim )
+        {
+          iNodeOnFace = iN;
+          iPosDim = nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition();
+        }
+      }
       // compute normal
       gp_Vec v01( nPnt[0], nPnt[1] ), v02( nPnt[0], nPnt[2] );
       if ( v01.SquareMagnitude() > RealSmall() &&
@@ -2499,7 +2525,7 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace)
       {
         Ne = v01 ^ v02;
         if (( normalOK = ( Ne.SquareMagnitude() > RealSmall() )))
-          uv = GetNodeUV( theFace, nPnt[0]._node, nPnt[2]._node, &normalOK );
+          uv = GetNodeUV( theFace, nPnt[iNodeOnFace]._node, 0, &normalOK );
       }
     }
   }
@@ -2646,6 +2672,76 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
   return tol;
 }
 
+//================================================================================
+/*!
+ * \brief Return an angle between two EDGEs sharing a common VERTEX with reference
+ *        of the FACE normal
+ *  \return double - the angle (between -Pi and Pi), negative if the angle is concave,
+ *                   1e100 in case of failure
+ *  \waring Care about order of the EDGEs and their orientation to be as they are
+ *          within the FACE! Don't pass degenerated EDGEs neither!
+ */
+//================================================================================
+
+double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
+                                     const TopoDS_Edge & theE2,
+                                     const TopoDS_Face & theFace)
+{
+  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 );
+    if ( c1.IsNull() || c2.IsNull() )
+      return angle;
+    gp_XY uv = c2d1->Value( p1 ).XY();
+    gp_Vec du, dv; gp_Pnt p;
+    surf->D1( uv.X(), uv.Y(), p, du, dv );
+    gp_Vec vec1, vec2, vecRef = du ^ dv;
+    int  nbLoops = 0;
+    double p1tmp = p1;
+    while ( vecRef.SquareMagnitude() < std::numeric_limits<double>::min() )
+    {
+      double dp = ( l - f ) / 1000.;
+      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;
+      if ( ++nbLoops > 10 )
+      {
+#ifdef _DEBUG_
+        cout << "SMESH_MesherHelper::GetAngle(): Captured in a sigularity" << endl;
+#endif
+        return angle;
+      }
+    }
+    if ( theFace.Orientation() == TopAbs_REVERSED )
+      vecRef.Reverse();
+    c1->D1( p1, p, vec1 );
+    c2->D1( p2, p, vec2 );
+    TopoDS_Face F = theFace;
+    if ( F.Orientation() == TopAbs_INTERNAL )
+      F.Orientation( TopAbs_FORWARD );
+    if ( theE1.Orientation() /*GetSubShapeOri( F, theE1 )*/ == TopAbs_REVERSED )
+      vec1.Reverse();
+    if ( theE2.Orientation() /*GetSubShapeOri( F, theE2 )*/ == TopAbs_REVERSED )
+      vec2.Reverse();
+    angle = vec1.AngleWithRef( vec2, vecRef );
+  }
+  catch (...)
+  {
+  }
+  return angle;
+}
+
 //================================================================================
 /*!
  * \brief Check if the first and last vertices of an edge are the same
@@ -3818,7 +3914,7 @@ namespace { // Structures used by FixQuadraticElements()
       {
         // check if the EDGE needs checking
         const TopoDS_Edge& edge = TopoDS::Edge( edgeIt.Current() );
-        if ( BRep_Tool::Degenerated( edge ) )
+        if ( SMESH_Algo::isDegenerated( edge ) )
           continue;
         if ( theHelper.IsRealSeam( edge ) &&
              edge.Orientation() == TopAbs_REVERSED )
@@ -4666,4 +4762,3 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
     }
   }
 }
-