SALOME_TESTS/Grids/smesh/2D_mesh_QuadranglePreference_00/A1

[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx

diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx
index a012d0292b9cca955567585276c1b7a8f097fe90..11dec7d1e436f1d4992eb11871e7ca2efa1e33d3 100644 (file)
--- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx
+++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx
@@ -54,6 +54,7 @@
 #include <Geom2d_Line.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
 #include <GeomAdaptor_Curve.hxx>
+#include <GeomLib.hxx>
 #include <Geom_Circle.hxx>
 #include <Geom_Curve.hxx>
 #include <Geom_Line.hxx>
@@ -1811,45 +1812,36 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
       totalNbFaces++;
       F = TopoDS::Face( s );
 
-      // IDEA: if there is a problem with finding a normal,
-      // we can compute an area-weighted sum of normals of all faces sharing the node
       gp_XY uv = helper.GetNodeUV( F, node, 0, &normOK );
       Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
-      surface->D1( uv.X(), uv.Y(), p, du,dv );
-      geomNorm = du ^ dv;
-      double size2 = geomNorm.SquareMagnitude();
-      if ( size2 < 1e-10 ) // singularity
       {
-        SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
-        while ( fIt->more() )
+        gp_Dir normal;
+        if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 )
+        {
+          geomNorm = normal;
+          normOK = true;
+        }
+        else // hard singularity
         {
-          const SMDS_MeshElement* f = fIt->next();
-          if ( editor.FindShape( f ) == *id )
+          SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
+          while ( fIt->more() )
           {
-            SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) geomNorm.XYZ(), /*normalized=*/false );
-            size2 = geomNorm.SquareMagnitude();
-            break;
+            const SMDS_MeshElement* f = fIt->next();
+            if ( editor.FindShape( f ) == *id )
+            {
+              SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) geomNorm.XYZ(), /*normalized=*/false );
+              if ( helper.IsReversedSubMesh( F ))
+                geomNorm.Reverse();
+              break;
+            }
           }
+          double size2 = geomNorm.SquareMagnitude();
+          if ( size2 > numeric_limits<double>::min() )
+            geomNorm /= sqrt( size2 );
+          else
+            normOK = false;
         }
-        // double ddu = 0, ddv = 0;
-        // if ( du.SquareMagnitude() > dv.SquareMagnitude() )
-        //   ddu = 1e-3;
-        // else
-        //   ddv = 1e-3;
-        // surface->D1( uv.X()+ddu, uv.Y()+ddv, p, du,dv );
-        // geomNorm = du ^ dv;
-        // size2 = geomNorm.SquareMagnitude();
-        // if ( size2 < 1e-10 )
-        // {
-        //   surface->D1( uv.X()-ddu, uv.Y()-ddv, p, du,dv );
-        //   geomNorm = du ^ dv;
-        //   size2 = geomNorm.SquareMagnitude();
-        // }
       }
-      if ( size2 > numeric_limits<double>::min() )
-        geomNorm /= sqrt( size2 );
-      else
-        normOK = false;
       if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
         geomNorm.Reverse();
       edge._normal += geomNorm.XYZ();
@@ -2387,8 +2379,8 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
     else
     {
       // smooth on FACE's
-      int step = 0, badNb = 0; moved = true;
-      while (( ++step <= 5 && moved ) || improved )
+      int step = 0, stepLimit = 5, badNb = 0; moved = true;
+      while (( ++step <= stepLimit && moved ) || improved )
       {
         dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
                      <<"_InfStep"<<nbSteps<<"_"<<step); // debug
@@ -2399,6 +2391,10 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
           moved |= data._edges[i]->Smooth(badNb);
         improved = ( badNb < oldBadNb );
 
+        // issue 22576. no bad faces but still there are intersections to fix
+        if ( improved && badNb == 0 )
+          stepLimit = step + 3;
+
         dumpFunctionEnd();
       }
       if ( badNb > 0 )
@@ -4205,7 +4201,9 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face&          F,
       trias [iSide].first  = badTrias[iTia];
       trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces,
                                                              & i1, & i2 );
-      if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 )
+      if (( ! trias[iSide].second ) ||
+          ( trias[iSide].second->NbCornerNodes() != 3 ) ||
+          ( ! sm->Contains( trias[iSide].second )))
         continue;
 
       // aspect ratio of an adjacent tria