Merge branch 'rnc/zcracks_doc_integration'

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

diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx
index c2a3c954b19c39410eb4105ca4160548b2cdd588..1f6d46263460186259fe57f14cafee07b51facb4 100644 (file)
--- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx
+++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx
@@ -1,9 +1,9 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // 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
@@ -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();
@@ -4205,7 +4197,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