X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=1f6d46263460186259fe57f14cafee07b51facb4;hp=7c9588671f657f7ccd4e9a03b2d3e25387e720eb;hb=b1f58b701eb4ec6fb4a14b070041f2c89e44723c;hpb=3c53917e76b725efef5eb713a79d905140643774

diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx
index 7c9588671..1f6d46263 100644
--- 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
@@ -35,12 +35,14 @@
 #include "SMESH_ControlsDef.hxx"
 #include "SMESH_Gen.hxx"
 #include "SMESH_Group.hxx"
+#include "SMESH_HypoFilter.hxx"
 #include "SMESH_Mesh.hxx"
 #include "SMESH_MeshAlgos.hxx"
 #include "SMESH_MesherHelper.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
 #include "SMESH_subMeshEventListener.hxx"
+#include "StdMeshers_FaceSide.hxx"
 
 #include <BRepAdaptor_Curve2d.hxx>
 #include <BRep_Tool.hxx>
@@ -52,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>
@@ -62,6 +65,7 @@
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_IndexedMapOfShape.hxx>
+#include <TopTools_ListOfShape.hxx>
 #include <TopTools_MapOfShape.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Edge.hxx>
@@ -381,8 +385,10 @@ namespace VISCOUS_3D
   {
     TopoDS_Shape                    _solid;
     const StdMeshers_ViscousLayers* _hyp;
+    TopoDS_Shape                    _hypShape;
     _MeshOfSolid*                   _proxyMesh;
     set<TGeomID>                    _reversedFaceIds;
+    set<TGeomID>                    _ignoreFaceIds;
 
     double                          _stepSize, _stepSizeCoeff;
     const SMDS_MeshNode*            _stepSizeNodes[2];
@@ -392,10 +398,10 @@ namespace VISCOUS_3D
     // iteration over the map is 5 time longer than over the vector
     vector< _LayerEdge* >           _edges;
 
-    // key: an id of shape (EDGE or VERTEX) shared by a FACE with
-    // layers and a FACE w/o layers
+    // key:   an id of shape (EDGE or VERTEX) shared by a FACE with
+    //        layers and a FACE w/o layers
     // value: the shape (FACE or EDGE) to shrink mesh on.
-    // _LayerEdge's basing on nodes on key shape are inflated along the value shape
+    //       _LayerEdge's basing on nodes on key shape are inflated along the value shape
     map< TGeomID, TopoDS_Shape >     _shrinkShape2Shape;
 
     // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID
@@ -413,7 +419,9 @@ namespace VISCOUS_3D
 
     _SolidData(const TopoDS_Shape&             s=TopoDS_Shape(),
                const StdMeshers_ViscousLayers* h=0,
-               _MeshOfSolid*                   m=0) :_solid(s), _hyp(h), _proxyMesh(m) {}
+               const TopoDS_Shape&             hs=TopoDS_Shape(),
+               _MeshOfSolid*                   m=0)
+      :_solid(s), _hyp(h), _hypShape(hs), _proxyMesh(m) {}
     ~_SolidData();
 
     Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge&    E,
@@ -433,7 +441,7 @@ namespace VISCOUS_3D
     //vector<const SMDS_MeshNode*> _nodesAround;
     vector<_Simplex>             _simplices; // for quality check
 
-    enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR };
+    enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR, TFI };
 
     bool Smooth(int&                  badNb,
                 Handle(Geom_Surface)& surface,
@@ -499,7 +507,11 @@ namespace VISCOUS_3D
     bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
                               SMESH_MesherHelper& helper,
                               const SMESHDS_SubMesh* faceSubMesh );
-    void fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper);
+    void fixBadFaces(const TopoDS_Face&          F,
+                     SMESH_MesherHelper&         helper,
+                     const bool                  is2D,
+                     const int                   step,
+                     set<const SMDS_MeshNode*> * involvedNodes=NULL);
     bool addBoundaryElements();
 
     bool error( const string& text, int solidID=-1 );
@@ -512,7 +524,6 @@ namespace VISCOUS_3D
     SMESH_ComputeErrorPtr _error;
 
     vector< _SolidData >  _sdVec;
-    set<TGeomID>          _ignoreShapeIds;
     int                   _tmpFaceID;
   };
   //--------------------------------------------------------------------------------
@@ -548,7 +559,7 @@ namespace VISCOUS_3D
     virtual vtkIdType GetVtkType() const                      { return -1; }
     virtual SMDSAbs_EntityType   GetEntityType() const        { return SMDSEntity_Last; }
     virtual SMDSAbs_GeometryType GetGeomType() const          { return SMDSGeom_TRIANGLE; }
-virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
+    virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
     { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
   };
   //--------------------------------------------------------------------------------
@@ -567,6 +578,44 @@ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
       _nn[3]=_le2->_nodes[0];
     }
   };
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Retriever of node coordinates either directly of from a surface by node UV.
+   * \warning Location of a surface is ignored
+   */
+  struct NodeCoordHelper
+  {
+    SMESH_MesherHelper&        _helper;
+    const TopoDS_Face&         _face;
+    Handle(Geom_Surface)       _surface;
+    gp_XYZ (NodeCoordHelper::* _fun)(const SMDS_MeshNode* n) const;
+
+    NodeCoordHelper(const TopoDS_Face& F, SMESH_MesherHelper& helper, bool is2D)
+      : _helper( helper ), _face( F )
+    {
+      if ( is2D )
+      {
+        TopLoc_Location loc;
+        _surface = BRep_Tool::Surface( _face, loc );
+      }
+      if ( _surface.IsNull() )
+        _fun = & NodeCoordHelper::direct;
+      else
+        _fun = & NodeCoordHelper::byUV;
+    }
+    gp_XYZ operator()(const SMDS_MeshNode* n) const { return (this->*_fun)( n ); }
+
+  private:
+    gp_XYZ direct(const SMDS_MeshNode* n) const
+    {
+      return SMESH_TNodeXYZ( n );
+    }
+    gp_XYZ byUV  (const SMDS_MeshNode* n) const
+    {
+      gp_XY uv = _helper.GetNodeUV( _face, n );
+      return _surface->Value( uv.X(), uv.Y() ).XYZ();
+    }
+  };
 } // namespace VISCOUS_3D
 
 //================================================================================
@@ -574,7 +623,7 @@ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
 //
 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
   :SMESH_Hypothesis(hypId, studyId, gen),
-   _isToIgnoreShapes(18), _nbLayers(1), _thickness(1), _stretchFactor(1)
+   _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1)
 {
   _name = StdMeshers_ViscousLayers::GetHypType();
   _param_algo_dim = -3; // auxiliary hyp used by 3D algos
@@ -643,7 +692,7 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
        << " " << _thickness
        << " " << _stretchFactor
        << " " << _shapeIds.size();
-  for ( unsigned i = 0; i < _shapeIds.size(); ++i )
+  for ( size_t i = 0; i < _shapeIds.size(); ++i )
     save << " " << _shapeIds[i];
   save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies.
   return save;
@@ -745,7 +794,7 @@ namespace
     // get average dir of edges going fromV
     gp_XYZ edgeDir;
     //if ( edges.size() > 1 )
-      for ( unsigned i = 0; i < edges.size(); ++i )
+      for ( size_t i = 0; i < edges.size(); ++i )
       {
         edgeDir = getEdgeDir( edges[i], fromV );
         double size2 = edgeDir.SquareModulus();
@@ -780,6 +829,8 @@ namespace
 
   bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper )
   {
+    // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 )
+    //   return true;
     gp_Vec2d drv1, drv2;
     gp_Pnt2d p;
     TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
@@ -801,7 +852,7 @@ namespace
         double cross = drv2 ^ drv1;
         if ( E.Orientation() == TopAbs_REVERSED )
           cross = -cross;
-        isConvex = ( cross < 1e-9 );
+        isConvex = ( cross > 0.1 ); //-1e-9 );
       }
       // check if concavity is strong enough to care about it
       //const double maxAngle = 5 * Standard_PI180;
@@ -831,6 +882,26 @@ namespace
         // }
       }
     }
+    // check angles at VERTEXes
+    TError error;
+    TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error );
+    for ( size_t iW = 0; iW < wires.size(); ++iW )
+    {
+      const int nbEdges = wires[iW]->NbEdges();
+      if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0)))
+        continue;
+      for ( int iE1 = 0; iE1 < nbEdges; ++iE1 )
+      {
+        if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue;
+        int iE2 = ( iE1 + 1 ) % nbEdges;
+        while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
+          iE2 = ( iE2 + 1 ) % nbEdges;
+        double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
+                                        wires[iW]->Edge( iE2 ), F );
+        if ( angle < -5. * M_PI / 180. )
+          return true;
+      }
+    }
     return false;
   }
   //--------------------------------------------------------------------------------
@@ -842,13 +913,15 @@ namespace
       const char* fname = "/tmp/viscous.py";
       cout << "execfile('"<<fname<<"')"<<endl;
       py = new ofstream(fname);
-      *py << "from smesh import *" << endl
-          << "meshSO = GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
-          << "mesh = Mesh( meshSO.GetObject() )"<<endl;
+      *py << "import SMESH" << endl
+          << "from salome.smesh import smeshBuilder" << endl
+          << "smesh  = smeshBuilder.New(salome.myStudy)" << endl
+          << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
+          << "mesh   = smesh.Mesh( meshSO.GetObject() )"<<endl;
     }
     void Finish() {
       if (py)
-        *py << "mesh.MakeGroup('Viscous Prisms',VOLUME,FT_ElemGeomType,'=',Geom_PENTA)"<<endl;
+        *py << "mesh.MakeGroup('Viscous Prisms',SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA)"<<endl;
       delete py; py=0;
     }
     ~PyDump() { Finish(); }
@@ -1002,7 +1075,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
   if ( !findFacesWithLayers() )
     return _error;
 
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     if ( ! makeLayer(_sdVec[i]) )
       return _error;
@@ -1041,6 +1114,7 @@ bool _ViscousBuilder::findSolidsWithLayers()
   _sdVec.reserve( allSolids.Extent());
 
   SMESH_Gen* gen = _mesh->GetGen();
+  SMESH_HypoFilter filter;
   for ( int i = 1; i <= allSolids.Extent(); ++i )
   {
     // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
@@ -1055,10 +1129,14 @@ bool _ViscousBuilder::findSolidsWithLayers()
       viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
     if ( viscHyp )
     {
+      TopoDS_Shape hypShape;
+      filter.Init( filter.Is( viscHyp ));
+      _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
+
       _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
                                                                 allSolids(i),
                                                                 /*toCreate=*/true);
-      _sdVec.push_back( _SolidData( allSolids(i), viscHyp, proxyMesh ));
+      _sdVec.push_back( _SolidData( allSolids(i), viscHyp, hypShape, proxyMesh ));
       _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
     }
   }
@@ -1077,44 +1155,69 @@ bool _ViscousBuilder::findSolidsWithLayers()
 
 bool _ViscousBuilder::findFacesWithLayers()
 {
+  SMESH_MesherHelper helper( *_mesh );
+  TopExp_Explorer exp;
+  TopTools_IndexedMapOfShape solids;
+
   // collect all faces to ignore defined by hyp
-  vector<TopoDS_Shape> ignoreFaces;
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
+    solids.Add( _sdVec[i]._solid );
+
     vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
-    for ( unsigned i = 0; i < ids.size(); ++i )
+    if ( _sdVec[i]._hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
     {
-      const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[i] );
-      if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
+      for ( size_t ii = 0; ii < ids.size(); ++ii )
       {
-        _ignoreShapeIds.insert( ids[i] );
-        ignoreFaces.push_back( s );
+        const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
+        if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
+          _sdVec[i]._ignoreFaceIds.insert( ids[ii] );
+      }
+    }
+    else // FACEs with layers are given
+    {
+      exp.Init( _sdVec[i]._solid, TopAbs_FACE );
+      for ( ; exp.More(); exp.Next() )
+      {
+        TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
+        if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
+          _sdVec[i]._ignoreFaceIds.insert( faceInd );
       }
     }
-  }
 
-  // ignore internal faces
-  SMESH_MesherHelper helper( *_mesh );
-  TopExp_Explorer exp;
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
-  {
-    exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
-    for ( ; exp.More(); exp.Next() )
+    // ignore internal FACEs if inlets and outlets are specified
     {
-      TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
-      if ( helper.NbAncestors( exp.Current(), *_mesh, TopAbs_SOLID ) > 1 )
-      {     
-        _ignoreShapeIds.insert( faceInd );
-        ignoreFaces.push_back( exp.Current() );
-        if ( helper.IsReversedSubMesh( TopoDS::Face( exp.Current() )))
+      TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
+      if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
+        TopExp::MapShapesAndAncestors( _sdVec[i]._hypShape,
+                                       TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
+
+      exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
+      for ( ; exp.More(); exp.Next() )
+      {
+        const TopoDS_Face& face = TopoDS::Face( exp.Current() );
+        if ( helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
+          continue;
+
+        const TGeomID faceInd = getMeshDS()->ShapeToIndex( face );
+        if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
+        {
+          int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
+          if ( nbSolids > 1 )
+            _sdVec[i]._ignoreFaceIds.insert( faceInd );
+        }
+
+        if ( helper.IsReversedSubMesh( face ))
+        {
           _sdVec[i]._reversedFaceIds.insert( faceInd );
+        }
       }
     }
   }
 
   // Find faces to shrink mesh on (solution 2 in issue 0020832);
   TopTools_IndexedMapOfShape shapes;
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     shapes.Clear();
     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
@@ -1134,18 +1237,35 @@ bool _ViscousBuilder::findFacesWithLayers()
       // check presence of layers on them
       int ignore[2];
       for ( int j = 0; j < 2; ++j )
-        ignore[j] = _ignoreShapeIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
-      if ( ignore[0] == ignore[1] ) continue; // nothing interesting
+        ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
+      if ( ignore[0] == ignore[1] )
+        continue; // nothing interesting
       TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
+      // check presence of layers on fWOL within an adjacent SOLID
+      PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
+      while ( const TopoDS_Shape* solid = sIt->next() )
+        if ( !solid->IsSame( _sdVec[i]._solid ))
+        {
+          int iSolid = solids.FindIndex( *solid );
+          int  iFace = getMeshDS()->ShapeToIndex( fWOL );
+          if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
+          {
+            _sdVec[i]._noShrinkFaces.insert( iFace );
+            fWOL.Nullify();
+          }
+        }
       // add edge to maps
-      TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
-      _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
+      if ( !fWOL.IsNull())
+      {
+        TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
+        _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
+      }
     }
   }
   // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
   // the algo of the SOLID sharing the FACE does not support it
   set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     TopTools_MapOfShape noShrinkVertices;
     map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
@@ -1162,7 +1282,7 @@ bool _ViscousBuilder::findFacesWithLayers()
         SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
         if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
         notShrinkFace = true;
-        for ( unsigned j = 0; j < _sdVec.size(); ++j )
+        for ( size_t j = 0; j < _sdVec.size(); ++j )
         {
           if ( _sdVec[j]._solid.IsSame( *solid ) )
             if ( _sdVec[j]._shrinkShape2Shape.count( edgeID ))
@@ -1195,10 +1315,10 @@ bool _ViscousBuilder::findFacesWithLayers()
       }
     }
   }
-      
+
   // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
 
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     shapes.Clear();
     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
@@ -1212,11 +1332,12 @@ bool _ViscousBuilder::findFacesWithLayers()
       while ( fIt->more())
       {
         const TopoDS_Shape* f = fIt->next();
-        const int         fID = getMeshDS()->ShapeToIndex( *f );
         if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
         {
           totalNbFaces++;
-          if ( _ignoreShapeIds.count ( fID ) && ! _sdVec[i]._noShrinkFaces.count( fID ))
+          const int fID = getMeshDS()->ShapeToIndex( *f );
+          if ( _sdVec[i]._ignoreFaceIds.count ( fID ) &&
+               !_sdVec[i]._noShrinkFaces.count( fID ))
             facesWOL.push_back( *f );
         }
       }
@@ -1226,42 +1347,42 @@ bool _ViscousBuilder::findFacesWithLayers()
       switch ( facesWOL.size() )
       {
       case 1:
+      {
+        helper.SetSubShape( facesWOL[0] );
+        if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
         {
-          helper.SetSubShape( facesWOL[0] );
-          if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
+          TopoDS_Shape seamEdge;
+          PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
+          while ( eIt->more() && seamEdge.IsNull() )
           {
-            TopoDS_Shape seamEdge;
-            PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
-            while ( eIt->more() && seamEdge.IsNull() )
-            {
-              const TopoDS_Shape* e = eIt->next();
-              if ( helper.IsRealSeam( *e ) )
-                seamEdge = *e;
-            }
-            if ( !seamEdge.IsNull() )
-            {
-              _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
-              break;
-            }
+            const TopoDS_Shape* e = eIt->next();
+            if ( helper.IsRealSeam( *e ) )
+              seamEdge = *e;
+          }
+          if ( !seamEdge.IsNull() )
+          {
+            _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
+            break;
           }
-          _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
-          break;
         }
+        _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
+        break;
+      }
       case 2:
+      {
+        // find an edge shared by 2 faces
+        PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
+        while ( eIt->more())
         {
-          // find an edge shared by 2 faces
-          PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
-          while ( eIt->more())
+          const TopoDS_Shape* e = eIt->next();
+          if ( helper.IsSubShape( *e, facesWOL[0]) &&
+               helper.IsSubShape( *e, facesWOL[1]))
           {
-            const TopoDS_Shape* e = eIt->next();
-            if ( helper.IsSubShape( *e, facesWOL[0]) &&
-                 helper.IsSubShape( *e, facesWOL[1]))
-            {
-              _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
-            }
+            _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
           }
-          break;
         }
+        break;
+      }
       default:
         return error("Not yet supported case", _sdVec[i]._index);
       }
@@ -1284,10 +1405,10 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   subIds = data._noShrinkFaces;
   TopExp_Explorer exp( data._solid, TopAbs_FACE );
   for ( ; exp.More(); exp.Next() )
-    if ( ! _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
     {
       SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
-      faceIds.insert( fSubM->GetId() );
+      if ( ! data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
+        faceIds.insert( fSubM->GetId() );
       SMESH_subMeshIteratorPtr subIt =
         fSubM->getDependsOnIterator(/*includeSelf=*/true, /*complexShapeFirst=*/false);
       while ( subIt->more() )
@@ -1301,7 +1422,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
   {
     TGeomID shapeInd = s2s->first;
-    for ( unsigned i = 0; i < _sdVec.size(); ++i )
+    for ( size_t i = 0; i < _sdVec.size(); ++i )
     {
       if ( _sdVec[i]._index == data._index ) continue;
       map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
@@ -1408,7 +1529,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 
   // Set target nodes into _Simplex and _2NearEdges
   TNode2Edge::iterator n2e;
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     if ( data._edges[i]->IsOnEdge())
       for ( int j = 0; j < 2; ++j )
@@ -1422,7 +1543,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
         data._edges[i]->_2neibors->_edges[j] = n2e->second;
       }
     else
-      for ( unsigned j = 0; j < data._edges[i]->_simplices.size(); ++j )
+      for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
       {
         _Simplex& s = data._edges[i]->_simplices[j];
         s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
@@ -1506,7 +1627,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
   SMESH_MesherHelper helper( *_mesh );
   bool ok = true;
 
-  for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS )
+  for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
   {
     vector<_LayerEdge*>& eS = edgesByGeom[iS];
     if ( eS.empty() ) continue;
@@ -1550,7 +1671,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
         if ( eE.empty() ) continue;
         if ( eE[0]->_sWOL.IsNull() )
         {
-          for ( unsigned i = 0; i < eE.size() && !needSmooth; ++i )
+          for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
             needSmooth = ( eE[i]->_cosin > 0.1 );
         }
         else
@@ -1558,7 +1679,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
           const TopoDS_Face& F1 = TopoDS::Face( S );
           const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
           const TopoDS_Edge& E  = TopoDS::Edge( eExp.Current() );
-          for ( unsigned i = 0; i < eE.size() && !needSmooth; ++i )
+          for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
           {
             gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
             gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
@@ -1597,7 +1718,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
   }
 
   // then the rest _LayerEdge's
-  for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS )
+  for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
   {
     vector<_LayerEdge*>& eVec = edgesByGeom[iS];
     data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
@@ -1689,18 +1810,38 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
       if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
         continue;
       totalNbFaces++;
-      //nbLayerFaces += subIds.count( *id );
       F = TopoDS::Face( s );
 
       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 > numeric_limits<double>::min() )
-        geomNorm /= sqrt( size2 );
-      else
-        normOK = false;
+      {
+        gp_Dir normal;
+        if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 )
+        {
+          geomNorm = normal;
+          normOK = true;
+        }
+        else // hard singularity
+        {
+          SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
+          while ( fIt->more() )
+          {
+            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;
+        }
+      }
       if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
         geomNorm.Reverse();
       edge._normal += geomNorm.XYZ();
@@ -1774,9 +1915,9 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
 
     if ( posType == SMDS_TOP_FACE )
     {
-      getSimplices( node, edge._simplices, _ignoreShapeIds, &data );
+      getSimplices( node, edge._simplices, data._ignoreFaceIds, &data );
       double avgNormProj = 0, avgLen = 0;
-      for ( unsigned i = 0; i < edge._simplices.size(); ++i )
+      for ( size_t i = 0; i < edge._simplices.size(); ++i )
       {
         gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev );
         avgNormProj += edge._normal * vec;
@@ -1961,14 +2102,15 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
                                     const _SolidData*    dataToCheckOri,
                                     const bool           toSort)
 {
+  simplices.clear();
   SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
   while ( fIt->more() )
   {
     const SMDS_MeshElement* f = fIt->next();
-    const TGeomID shapeInd = f->getshapeId();
+    const TGeomID    shapeInd = f->getshapeId();
     if ( ingnoreShapes.count( shapeInd )) continue;
     const int nbNodes = f->NbCornerNodes();
-    int srcInd = f->GetNodeIndex( node );
+    const int  srcInd = f->GetNodeIndex( node );
     const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes ));
     const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes ));
     const SMDS_MeshNode* nOpp  = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+2, nbNodes ));
@@ -2006,7 +2148,7 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
 void _ViscousBuilder::makeGroupOfLE()
 {
 #ifdef _DEBUG_
-  for ( unsigned i = 0 ; i < _sdVec.size(); ++i )
+  for ( size_t i = 0 ; i < _sdVec.size(); ++i )
   {
     if ( _sdVec[i]._edges.empty() ) continue;
 //     string name = SMESH_Comment("_LayerEdge's_") << i;
@@ -2016,10 +2158,10 @@ void _ViscousBuilder::makeGroupOfLE()
 //     SMESHDS_Mesh* mDS = _mesh->GetMeshDS();
 
     dumpFunction( SMESH_Comment("make_LayerEdge_") << i );
-    for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j )
+    for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
     {
       _LayerEdge* le = _sdVec[i]._edges[j];
-      for ( unsigned iN = 1; iN < le->_nodes.size(); ++iN )
+      for ( size_t iN = 1; iN < le->_nodes.size(); ++iN )
         dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<le->_nodes[iN-1]->GetID()
                 << ", " << le->_nodes[iN]->GetID() <<"])");
       //gDS->SMDSGroup().Add( mDS->AddEdge( le->_nodes[iN-1], le->_nodes[iN]));
@@ -2027,7 +2169,7 @@ void _ViscousBuilder::makeGroupOfLE()
     dumpFunctionEnd();
 
     dumpFunction( SMESH_Comment("makeNormals") << i );
-    for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j )
+    for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
     {
       _LayerEdge& edge = *_sdVec[i]._edges[j];
       SMESH_TNodeXYZ nXYZ( edge._nodes[0] );
@@ -2081,7 +2223,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
   auto_ptr<SMESH_ElementSearcher> searcher
     ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
                                            data._proxyMesh->GetFaces( data._solid )) );
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     if ( data._edges[i]->IsOnEdge() ) continue;
     data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
@@ -2116,7 +2258,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
 
     // Elongate _LayerEdge's
     dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
-    for ( unsigned i = 0; i < data._edges.size(); ++i )
+    for ( size_t i = 0; i < data._edges.size(); ++i )
     {
       data._edges[i]->SetNewLength( curThick, helper );
     }
@@ -2132,7 +2274,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
       if ( nbSteps > 0 )
       {
         dumpFunction(SMESH_Comment("invalidate")<<data._index<<"_step"<<nbSteps); // debug
-        for ( unsigned i = 0; i < data._edges.size(); ++i )
+        for ( size_t i = 0; i < data._edges.size(); ++i )
         {
           data._edges[i]->InvalidateStep( nbSteps+1 );
         }
@@ -2144,7 +2286,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
 
     // Evaluate achieved thickness
     avgThick = 0;
-    for ( unsigned i = 0; i < data._edges.size(); ++i )
+    for ( size_t i = 0; i < data._edges.size(); ++i )
       avgThick += data._edges[i]->_len;
     avgThick /= data._edges.size();
 #ifdef __myDEBUG
@@ -2192,7 +2334,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
   TopoDS_Face F;
 
   int iBeg, iEnd = 0;
-  for ( unsigned iS = 0; iS < data._endEdgeToSmooth.size(); ++iS )
+  for ( size_t iS = 0; iS < data._endEdgeToSmooth.size(); ++iS )
   {
     iBeg = iEnd;
     iEnd = data._endEdgeToSmooth[ iS ];
@@ -2258,7 +2400,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         {
           _LayerEdge* edge = data._edges[i];
           SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
-          for ( unsigned j = 0; j < edge->_simplices.size(); ++j )
+          for ( size_t j = 0; j < edge->_simplices.size(); ++j )
             if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
             {
               cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
@@ -2287,7 +2429,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
   const SMDS_MeshElement* closestFace = 0;
   int iLE = 0;
 #endif
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
       return false;
@@ -2533,6 +2675,10 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData&           data,
 
     if ( F.IsNull() ) // 3D
     {
+      if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
+           data._edges[iTo-1]->_2neibors->_nodes[1] )
+        return true; // closed EDGE - nothing to do
+
       return false; // TODO ???
     }
     else // 2D
@@ -2594,7 +2740,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
     vector< const SMDS_MeshNode*> nodes(4); // of a tmp mesh face
 
     dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<<data._index);
-    for ( unsigned i = 0; i < data._edges.size(); ++i )
+    for ( size_t i = 0; i < data._edges.size(); ++i )
     {
       _LayerEdge* edge = data._edges[i];
       if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
@@ -2611,7 +2757,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         }
         // look for a _LayerEdge containg tgt2
 //         _LayerEdge* neiborEdge = 0;
-//         unsigned di = 0; // check _edges[i+di] and _edges[i-di]
+//         size_t di = 0; // check _edges[i+di] and _edges[i-di]
 //         while ( !neiborEdge && ++di <= data._edges.size() )
 //         {
 //           if ( i+di < data._edges.size() && data._edges[i+di]->_nodes.back() == tgt2 )
@@ -2650,7 +2796,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   TLEdge2LEdgeSet edge2CloseEdge;
 
   const double eps = data._epsilon * data._epsilon;
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     _LayerEdge* edge = data._edges[i];
     if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
@@ -2821,7 +2967,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   // 2) Check absence of intersections
   // TODO?
 
-  for ( unsigned i = 0 ; i < tmpFaces.size(); ++i )
+  for ( size_t i = 0 ; i < tmpFaces.size(); ++i )
     delete tmpFaces[i];
 
   return true;
@@ -2847,7 +2993,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
   bool segmentIntersected = false;
   distance = Precision::Infinite();
   int iFace = -1; // intersected face
-  for ( unsigned j = 0 ; j < suspectFaces.size() && !segmentIntersected; ++j )
+  for ( size_t j = 0 ; j < suspectFaces.size() && !segmentIntersected; ++j )
   {
     const SMDS_MeshElement* face = suspectFaces[j];
     if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
@@ -3135,7 +3281,7 @@ bool _LayerEdge::Smooth(int& badNb)
 
   // compute new position for the last _pos
   gp_XYZ newPos (0,0,0);
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
   newPos /= _simplices.size();
 
@@ -3158,11 +3304,11 @@ bool _LayerEdge::Smooth(int& badNb)
   // count quality metrics (orientation) of tetras around _tgtNode
   int nbOkBefore = 0;
   SMESH_TNodeXYZ tgtXYZ( _nodes.back() );
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ );
 
   int nbOkAfter = 0;
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos );
 
   if ( nbOkAfter < nbOkBefore )
@@ -3288,14 +3434,14 @@ bool _ViscousBuilder::refine(_SolidData& data)
   gp_XY uv;
   bool isOnEdge;
 
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     _LayerEdge& edge = *data._edges[i];
 
     // get accumulated length of segments
     vector< double > segLen( edge._pos.size() );
     segLen[0] = 0.0;
-    for ( unsigned j = 1; j < edge._pos.size(); ++j )
+    for ( size_t j = 1; j < edge._pos.size(); ++j )
       segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
 
     // allocate memory for new nodes if it is not yet refined
@@ -3343,8 +3489,8 @@ bool _ViscousBuilder::refine(_SolidData& data)
 
     // create intermediate nodes
     double hSum = 0, hi = h0/f;
-    unsigned iSeg = 1;
-    for ( unsigned iStep = 1; iStep < edge._nodes.size(); ++iStep )
+    size_t iSeg = 1;
+    for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep )
     {
       // compute an intermediate position
       hi *= f;
@@ -3411,7 +3557,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
 
   if ( !getMeshDS()->IsEmbeddedMode() )
     // Log node movement
-    for ( unsigned i = 0; i < data._edges.size(); ++i )
+    for ( size_t i = 0; i < data._edges.size(); ++i )
     {
       _LayerEdge& edge = *data._edges[i];
       SMESH_TNodeXYZ p ( edge._nodes.back() );
@@ -3425,7 +3571,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
   TopExp_Explorer exp( data._solid, TopAbs_FACE );
   for ( ; exp.More(); exp.Next() )
   {
-    if ( _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
+    if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
       continue;
     SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
     SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
@@ -3476,7 +3622,7 @@ bool _ViscousBuilder::shrink()
   // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
   // inflated along FACE or EDGE)
   map< TGeomID, _SolidData* > f2sdMap;
-  for ( unsigned i = 0 ; i < _sdVec.size(); ++i )
+  for ( size_t i = 0 ; i < _sdVec.size(); ++i )
   {
     _SolidData& data = _sdVec[i];
     TopTools_MapOfShape FFMap;
@@ -3575,9 +3721,15 @@ bool _ViscousBuilder::shrink()
       }
     }
 
+    dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
+    SMDS_ElemIteratorPtr fIt = smDS->GetElements();
+    while ( fIt->more() )
+      if ( const SMDS_MeshElement* f = fIt->next() )
+        dumpChangeNodes( f );
+
     // Replace source nodes by target nodes in mesh faces to shrink
     const SMDS_MeshNode* nodes[20];
-    for ( unsigned i = 0; i < lEdges.size(); ++i )
+    for ( size_t i = 0; i < lEdges.size(); ++i )
     {
       _LayerEdge& edge = *lEdges[i];
       const SMDS_MeshNode* srcNode = edge._nodes[0];
@@ -3588,10 +3740,10 @@ bool _ViscousBuilder::shrink()
         const SMDS_MeshElement* f = fIt->next();
         if ( !smDS->Contains( f ))
           continue;
-        SMDS_ElemIteratorPtr nIt = f->nodesIterator();
-        for ( int iN = 0; iN < f->NbNodes(); ++iN )
+        SMDS_NodeIteratorPtr nIt = f->nodeIterator();
+        for ( int iN = 0; nIt->more(); ++iN )
         {
-          const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
+          const SMDS_MeshNode* n = nIt->next();
           nodes[iN] = ( n == srcNode ? tgtNode : n );
         }
         helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
@@ -3604,9 +3756,8 @@ bool _ViscousBuilder::shrink()
     // Create _SmoothNode's on face F
     vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
     {
-      dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
       const bool sortSimplices = isConcaveFace;
-      for ( unsigned i = 0; i < smoothNodes.size(); ++i )
+      for ( size_t i = 0; i < smoothNodes.size(); ++i )
       {
         const SMDS_MeshNode* n = smoothNodes[i];
         nodesToSmooth[ i ]._node = n;
@@ -3616,14 +3767,13 @@ bool _ViscousBuilder::shrink()
         helper.GetNodeUV( F, n, 0, &isOkUV);
         dumpMove( n );
       }
-      dumpFunctionEnd();
     }
     //if ( nodesToSmooth.empty() ) continue;
 
-    // Find EDGE's to shrink
+    // Find EDGE's to shrink and set simpices to LayerEdge's
     set< _Shrinker1D* > eShri1D;
     {
-      for ( unsigned i = 0; i < lEdges.size(); ++i )
+      for ( size_t i = 0; i < lEdges.size(); ++i )
       {
         _LayerEdge* edge = lEdges[i];
         if ( edge->_sWOL.ShapeType() == TopAbs_EDGE )
@@ -3637,6 +3787,23 @@ bool _ViscousBuilder::shrink()
           // srinked while srinking another FACE
           srinker.RestoreParams();
         }
+        getSimplices( /*tgtNode=*/edge->_nodes.back(), edge->_simplices, ignoreShapes );
+      }
+    }
+
+    bool toFixTria = false; // to improve quality of trias by diagonal swap
+    if ( isConcaveFace )
+    {
+      const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
+      if ( hasTria != hasQuad ) {
+        toFixTria = hasTria;
+      }
+      else {
+        set<int> nbNodesSet;
+        SMDS_ElemIteratorPtr fIt = smDS->GetElements();
+        while ( fIt->more() && nbNodesSet.size() < 2 )
+          nbNodesSet.insert( fIt->next()->NbCornerNodes() );
+        toFixTria = ( *nbNodesSet.begin() == 3 );
       }
     }
 
@@ -3647,14 +3814,15 @@ bool _ViscousBuilder::shrink()
     bool shrinked = true;
     int badNb, shriStep=0, smooStep=0;
     _SmoothNode::SmoothType smoothType
-      = isConcaveFace ? _SmoothNode::CENTROIDAL : _SmoothNode::LAPLACIAN;
+      = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN;
     while ( shrinked )
     {
+      shriStep++;
       // Move boundary nodes (actually just set new UV)
       // -----------------------------------------------
-      dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep++ ); // debug
+      dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep ); // debug
       shrinked = false;
-      for ( unsigned i = 0; i < lEdges.size(); ++i )
+      for ( size_t i = 0; i < lEdges.size(); ++i )
       {
         shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
       }
@@ -3669,7 +3837,7 @@ bool _ViscousBuilder::shrink()
       // Smoothing in 2D
       // -----------------
       int nbNoImpSteps = 0;
-      bool moved = true;
+      bool       moved = true;
       badNb = 1;
       while (( nbNoImpSteps < 5 && badNb > 0) && moved)
       {
@@ -3678,7 +3846,7 @@ bool _ViscousBuilder::shrink()
         int oldBadNb = badNb;
         badNb = 0;
         moved = false;
-        for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
+        for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
         {
           moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
                                             smoothType, /*set3D=*/isConcaveFace);
@@ -3692,7 +3860,43 @@ bool _ViscousBuilder::shrink()
       }
       if ( badNb > 0 )
         return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first );
-    }
+      if ( shriStep > 200 )
+        return error(SMESH_Comment("Infinite loop at shrinking 2D mesh on face ") << f2sd->first );
+
+      // Fix narrow triangles by swapping diagonals
+      // ---------------------------------------
+      if ( toFixTria )
+      {
+        set<const SMDS_MeshNode*> usedNodes;
+        fixBadFaces( F, helper, /*is2D=*/true, shriStep, & usedNodes); // swap diagonals
+
+        // update working data
+        set<const SMDS_MeshNode*>::iterator n;
+        for ( size_t i = 0; i < nodesToSmooth.size() && !usedNodes.empty(); ++i )
+        {
+          n = usedNodes.find( nodesToSmooth[ i ]._node );
+          if ( n != usedNodes.end())
+          {
+            getSimplices( nodesToSmooth[ i ]._node,
+                          nodesToSmooth[ i ]._simplices,
+                          ignoreShapes, NULL,
+                          /*sortSimplices=*/ smoothType == _SmoothNode::ANGULAR );
+            usedNodes.erase( n );
+          }
+        }
+        for ( size_t i = 0; i < lEdges.size() && !usedNodes.empty(); ++i )
+        {
+          n = usedNodes.find( /*tgtNode=*/ lEdges[i]->_nodes.back() );
+          if ( n != usedNodes.end())
+          {
+            getSimplices( lEdges[i]->_nodes.back(),
+                          lEdges[i]->_simplices,
+                          ignoreShapes );
+            usedNodes.erase( n );
+          }
+        }
+      }
+    } // while ( shrinked )
 
     // No wrongly shaped faces remain; final smooth. Set node XYZ.
     bool isStructuredFixed = false;
@@ -3700,12 +3904,18 @@ bool _ViscousBuilder::shrink()
       isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F );
     if ( !isStructuredFixed )
     {
-      if ( isConcaveFace )
-        fixBadFaces( F, helper ); // fix narrow faces by swapping diagonals
-      for ( int st = /*highQuality ? 10 :*/ 3; st; --st )
+      if ( isConcaveFace ) // fix narrow faces by swapping diagonals
+        fixBadFaces( F, helper, /*is2D=*/false, ++shriStep );
+
+      for ( int st = 3; st; --st )
       {
+        switch( st ) {
+        case 1: smoothType = _SmoothNode::LAPLACIAN; break;
+        case 2: smoothType = _SmoothNode::LAPLACIAN; break;
+        case 3: smoothType = _SmoothNode::ANGULAR; break;
+        }
         dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
-        for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
+        for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
         {
           nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
                                    smoothType,/*set3D=*/st==1 );
@@ -3718,7 +3928,7 @@ bool _ViscousBuilder::shrink()
 
     if ( !getMeshDS()->IsEmbeddedMode() )
       // Log node movement
-      for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
+      for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
       {
         SMESH_TNodeXYZ p ( nodesToSmooth[i]._node );
         getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() );
@@ -3762,54 +3972,54 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0);
     edge._len = uvLen;
 
-    // IMPORTANT to have src nodes NOT yet REPLACED by tgt nodes in shrinked faces
-    vector<const SMDS_MeshElement*> faces;
-    multimap< double, const SMDS_MeshNode* > proj2node;
-    SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
-    while ( fIt->more() )
-    {
-      const SMDS_MeshElement* f = fIt->next();
-      if ( faceSubMesh->Contains( f ))
-        faces.push_back( f );
-    }
-    for ( unsigned i = 0; i < faces.size(); ++i )
-    {
-      const int nbNodes = faces[i]->NbCornerNodes();
-      for ( int j = 0; j < nbNodes; ++j )
-      {
-        const SMDS_MeshNode* n = faces[i]->GetNode(j);
-        if ( n == srcNode ) continue;
-        if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
-             ( faces.size() > 1 || nbNodes > 3 ))
-          continue;
-        gp_Pnt2d uv = helper.GetNodeUV( F, n );
-        gp_Vec2d uvDirN( srcUV, uv );
-        double proj = uvDirN * uvDir;
-        proj2node.insert( make_pair( proj, n ));
-      }
-    }
+    // // IMPORTANT to have src nodes NOT yet REPLACED by tgt nodes in shrinked faces
+    // vector<const SMDS_MeshElement*> faces;
+    // multimap< double, const SMDS_MeshNode* > proj2node;
+    // SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
+    // while ( fIt->more() )
+    // {
+    //   const SMDS_MeshElement* f = fIt->next();
+    //   if ( faceSubMesh->Contains( f ))
+    //     faces.push_back( f );
+    // }
+    // for ( size_t i = 0; i < faces.size(); ++i )
+    // {
+    //   const int nbNodes = faces[i]->NbCornerNodes();
+    //   for ( int j = 0; j < nbNodes; ++j )
+    //   {
+    //     const SMDS_MeshNode* n = faces[i]->GetNode(j);
+    //     if ( n == srcNode ) continue;
+    //     if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
+    //          ( faces.size() > 1 || nbNodes > 3 ))
+    //       continue;
+    //     gp_Pnt2d uv = helper.GetNodeUV( F, n );
+    //     gp_Vec2d uvDirN( srcUV, uv );
+    //     double proj = uvDirN * uvDir;
+    //     proj2node.insert( make_pair( proj, n ));
+    //   }
+    // }
 
-    multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd;
-    const double       minProj = p2n->first;
-    const double projThreshold = 1.1 * uvLen;
-    if ( minProj > projThreshold )
-    {
-      // tgtNode is located so that it does not make faces with wrong orientation
-      return true;
-    }
+    // multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd;
+    // const double       minProj = p2n->first;
+    // const double projThreshold = 1.1 * uvLen;
+    // if ( minProj > projThreshold )
+    // {
+    //   // tgtNode is located so that it does not make faces with wrong orientation
+    //   return true;
+    // }
     edge._pos.resize(1);
     edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
 
     // store most risky nodes in _simplices
-    p2nEnd = proj2node.lower_bound( projThreshold );
-    int nbSimpl = ( std::distance( p2n, p2nEnd ) + 1) / 2;
-    edge._simplices.resize( nbSimpl );
-    for ( int i = 0; i < nbSimpl; ++i )
-    {
-      edge._simplices[i]._nPrev = p2n->second;
-      if ( ++p2n != p2nEnd )
-        edge._simplices[i]._nNext = p2n->second;
-    }
+    // p2nEnd = proj2node.lower_bound( projThreshold );
+    // int nbSimpl = ( std::distance( p2n, p2nEnd ) + 1) / 2;
+    // edge._simplices.resize( nbSimpl );
+    // for ( int i = 0; i < nbSimpl; ++i )
+    // {
+    //   edge._simplices[i]._nPrev = p2n->second;
+    //   if ( ++p2n != p2nEnd )
+    //     edge._simplices[i]._nNext = p2n->second;
+    // }
     // set UV of source node to target node
     SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
     pos->SetUParameter( srcUV.X() );
@@ -3917,23 +4127,28 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
  */
 //================================================================================
 
-void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper)
+void _ViscousBuilder::fixBadFaces(const TopoDS_Face&          F,
+                                  SMESH_MesherHelper&         helper,
+                                  const bool                  is2D,
+                                  const int                   step,
+                                  set<const SMDS_MeshNode*> * involvedNodes)
 {
   SMESH::Controls::AspectRatio qualifier;
   SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3);
-  const double maxAspectRatio = 4.;
+  const double maxAspectRatio = is2D ? 4. : 2;
+  NodeCoordHelper xyz( F, helper, is2D );
 
   // find bad triangles
 
   vector< const SMDS_MeshElement* > badTrias;
   vector< double >                  badAspects;
-  SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F );
+  SMESHDS_SubMesh*      sm = helper.GetMeshDS()->MeshElements( F );
   SMDS_ElemIteratorPtr fIt = sm->GetElements();
   while ( fIt->more() )
   {
     const SMDS_MeshElement * f = fIt->next();
     if ( f->NbCornerNodes() != 3 ) continue;
-    for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = SMESH_TNodeXYZ( f->GetNode(iP));
+    for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = xyz( f->GetNode(iP));
     double aspect = qualifier.GetValue( points );
     if ( aspect > maxAspectRatio )
     {
@@ -3941,6 +4156,18 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
       badAspects.push_back( aspect );
     }
   }
+  if ( step == 1 )
+  {
+    dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
+    SMDS_ElemIteratorPtr fIt = sm->GetElements();
+    while ( fIt->more() )
+    {
+      const SMDS_MeshElement * f = fIt->next();
+      if ( f->NbCornerNodes() == 3 )
+        dumpChangeNodes( f );
+    }
+    dumpFunctionEnd();
+  }
   if ( badTrias.empty() )
     return;
 
@@ -3956,9 +4183,10 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
     double  aspRatio [3];
     int i1, i2, i3;
 
-    involvedFaces.insert( badTrias[iTia] );
+    if ( !involvedFaces.insert( badTrias[iTia] ).second )
+      continue;
     for ( int iP = 0; iP < 3; ++iP )
-      points(iP+1) = SMESH_TNodeXYZ( badTrias[iTia]->GetNode(iP));
+      points(iP+1) = xyz( badTrias[iTia]->GetNode(iP));
 
     // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping
     int bestCouple = -1;
@@ -3969,12 +4197,14 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
       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
       for ( int iP = 0; iP < 3; ++iP )
-        points2(iP+1) = SMESH_TNodeXYZ( trias[iSide].second->GetNode(iP));
+        points2(iP+1) = xyz( trias[iSide].second->GetNode(iP));
       double aspectInit = qualifier.GetValue( points2 );
 
       // arrange nodes as after diag-swaping
@@ -3991,6 +4221,12 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
       if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] )
         continue;
 
+      // prevent inversion of a triangle
+      gp_Vec norm1 = gp_Vec( points1(1), points1(3) ) ^ gp_Vec( points1(1), points1(2) );
+      gp_Vec norm2 = gp_Vec( points2(1), points2(3) ) ^ gp_Vec( points2(1), points2(2) );
+      if ( norm1 * norm2 < 0. && norm1.Angle( norm2 ) > 70./180.*M_PI )
+        continue;
+
       if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] )
         bestCouple = iSide;
     }
@@ -4011,17 +4247,25 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
   // swap diagonals
 
   SMESH_MeshEditor editor( helper.GetMesh() );
-  dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
+  dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
   for ( size_t i = 0; i < triaCouples.size(); ++i )
   {
     dumpChangeNodes( triaCouples[i].first );
     dumpChangeNodes( triaCouples[i].second );
     editor.InverseDiag( triaCouples[i].first, triaCouples[i].second );
   }
-  dumpFunctionEnd();
+
+  if ( involvedNodes )
+    for ( size_t i = 0; i < triaCouples.size(); ++i )
+    {
+      involvedNodes->insert( triaCouples[i].first->begin_nodes(),
+                             triaCouples[i].first->end_nodes() );
+      involvedNodes->insert( triaCouples[i].second->begin_nodes(),
+                             triaCouples[i].second->end_nodes() );
+    }
 
   // just for debug dump resulting triangles
-  dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID());
+  dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
   for ( size_t i = 0; i < triaCouples.size(); ++i )
   {
     dumpChangeNodes( triaCouples[i].first );
@@ -4047,41 +4291,40 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
   if ( _sWOL.ShapeType() == TopAbs_FACE )
   {
     gp_XY    curUV = helper.GetNodeUV( F, tgtNode );
-    gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y());
+    gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y() );
     gp_Vec2d uvDir( _normal.X(), _normal.Y() );
     const double uvLen = tgtUV.Distance( curUV );
+    const double kSafe = Max( 0.5, 1. - 0.1 * _simplices.size() );
 
     // Select shrinking step such that not to make faces with wrong orientation.
-    const double kSafe = 0.8;
-    const double minStepSize = uvLen / 10;
     double stepSize = uvLen;
-    for ( unsigned i = 0; i < _simplices.size(); ++i )
+    for ( size_t i = 0; i < _simplices.size(); ++i )
     {
-      const SMDS_MeshNode* nn[2] = { _simplices[i]._nPrev, _simplices[i]._nNext };
-      for ( int j = 0; j < 2; ++j )
-        if ( const SMDS_MeshNode* n = nn[j] )
-        {
-          gp_XY uv = helper.GetNodeUV( F, n );
-          gp_Vec2d uvDirN( curUV, uv );
-          double proj = uvDirN * uvDir * kSafe;
-          if ( proj < stepSize && proj > minStepSize )
-            stepSize = proj;
-          else if ( proj < minStepSize )
-            stepSize = minStepSize;
-        }
+      // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
+      gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev );
+      gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext );
+      gp_XY dirN = uvN2 - uvN1;
+      double det = uvDir.Crossed( dirN );
+      if ( Abs( det )  < std::numeric_limits<double>::min() ) continue;
+      gp_XY dirN2Cur = curUV - uvN1;
+      double step = dirN.Crossed( dirN2Cur ) / det;
+      if ( step > 0 )
+        stepSize = Min( step, stepSize );
     }
-
     gp_Pnt2d newUV;
-    if ( uvLen - stepSize < _len / 20. )
+    if ( uvLen - stepSize < _len / 200. )
     {
       newUV = tgtUV;
       _pos.clear();
     }
+    else if ( stepSize > 0 )
+    {
+      newUV = curUV + uvDir.XY() * stepSize * kSafe;
+    }
     else
     {
-      newUV = curUV + uvDir.XY() * stepSize;
+      return true;
     }
-
     SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
     pos->SetUParameter( newUV.X() );
     pos->SetVParameter( newUV.Y() );
@@ -4145,30 +4388,24 @@ bool _SmoothNode::Smooth(int&                  badNb,
 
   // compute new UV for the node
   gp_XY newPos (0,0);
-/*  if ( how == ANGULAR && _simplices.size() == 4 )
+  if ( how == TFI && _simplices.size() == 4 )
   {
-    vector<gp_XY> corners; corners.reserve(4);
+    gp_XY corners[4];
     for ( size_t i = 0; i < _simplices.size(); ++i )
       if ( _simplices[i]._nOpp )
-        corners.push_back( helper.GetNodeUV( face, _simplices[i]._nOpp, _node ));
-    if ( corners.size() == 4 )
-    {
-      newPos = helper.calcTFI
-        ( 0.5, 0.5,
-          corners[0], corners[1], corners[2], corners[3],
-          uv[1], uv[2], uv[3], uv[0] );
-    }
-    // vector<gp_XY> p( _simplices.size() * 2 + 1 );
-    // p.clear();
-    // for ( size_t i = 0; i < _simplices.size(); ++i )
-    // {
-    //   p.push_back( uv[i] );
-    //   if ( _simplices[i]._nOpp )
-    //     p.push_back( helper.GetNodeUV( face, _simplices[i]._nOpp, _node ));
-    // }
-    // newPos = computeAngularPos( p, helper.GetNodeUV( face, _node ), refSign );
+        corners[i] = helper.GetNodeUV( face, _simplices[i]._nOpp, _node );
+      else
+        throw SALOME_Exception(LOCALIZED("TFI smoothing: _Simplex::_nOpp not set!"));
+
+    newPos = helper.calcTFI ( 0.5, 0.5,
+                              corners[0], corners[1], corners[2], corners[3],
+                              uv[1], uv[2], uv[3], uv[0] );
+  }
+  else if ( how == ANGULAR )
+  {
+    newPos = computeAngularPos( uv, helper.GetNodeUV( face, _node ), refSign );
   }
-  else*/ if ( how == CENTROIDAL && _simplices.size() > 3 )
+  else if ( how == CENTROIDAL && _simplices.size() > 3 )
   {
     // average centers of diagonals wieghted with their reciprocal lengths
     if ( _simplices.size() == 4 )
@@ -4199,7 +4436,6 @@ bool _SmoothNode::Smooth(int&                  badNb,
   else
   {
     // Laplacian smooth
-    //isCentroidal = false;
     for ( size_t i = 0; i < _simplices.size(); ++i )
       newPos += uv[i];
     newPos /= _simplices.size();
@@ -4208,17 +4444,15 @@ bool _SmoothNode::Smooth(int&                  badNb,
   // count quality metrics (orientation) of triangles around the node
   int nbOkBefore = 0;
   gp_XY tgtUV = helper.GetNodeUV( face, _node );
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign );
 
   int nbOkAfter = 0;
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign );
 
   if ( nbOkAfter < nbOkBefore )
   {
-    // if ( isCentroidal )
-    //   return Smooth( badNb, surface, helper, refSign, !isCentroidal, set3D );
     badNb += _simplices.size() - nbOkBefore;
     return false;
   }
@@ -4253,22 +4487,22 @@ gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
 {
   uv.push_back( uv.front() );
 
-  vector< gp_XY > edgeDir( uv.size() );
+  vector< gp_XY >  edgeDir ( uv.size() );
   vector< double > edgeSize( uv.size() );
   for ( size_t i = 1; i < edgeDir.size(); ++i )
   {
-    edgeDir[i-1] = uv[i] - uv[i-1];
+    edgeDir [i-1] = uv[i] - uv[i-1];
     edgeSize[i-1] = edgeDir[i-1].Modulus();
     if ( edgeSize[i-1] < numeric_limits<double>::min() )
       edgeDir[i-1].SetX( 100 );
     else
       edgeDir[i-1] /= edgeSize[i-1] * refSign;
   }
-  edgeDir.back() = edgeDir.front();
+  edgeDir.back()  = edgeDir.front();
   edgeSize.back() = edgeSize.front();
 
-  gp_XY newPos(0,0);
-  int nbEdges = 0;
+  gp_XY  newPos(0,0);
+  int    nbEdges = 0;
   double sumSize = 0;
   for ( size_t i = 1; i < edgeDir.size(); ++i )
   {
@@ -4288,7 +4522,7 @@ gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
     }
     bisec /= bisecSize;
 
-    gp_XY  dirToN = uvToFix - p;
+    gp_XY  dirToN  = uvToFix - p;
     double distToN = dirToN.Modulus();
     if ( bisec * dirToN < 0 )
       distToN = -distToN;
@@ -4309,7 +4543,7 @@ gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
 
 _SolidData::~_SolidData()
 {
-  for ( unsigned i = 0; i < _edges.size(); ++i )
+  for ( size_t i = 0; i < _edges.size(); ++i )
   {
     if ( _edges[i] && _edges[i]->_2neibors )
       delete _edges[i]->_2neibors;
@@ -4382,7 +4616,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper )
   {
     // remove target node of the _LayerEdge from _nodes
     int nbFound = 0;
-    for ( unsigned i = 0; i < _nodes.size(); ++i )
+    for ( size_t i = 0; i < _nodes.size(); ++i )
       if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 )
         _nodes[i] = 0, nbFound++;
     if ( nbFound == _nodes.size() )
@@ -4420,7 +4654,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
       l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
     double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l );
 
-    for ( unsigned i = 0; i < _nodes.size(); ++i )
+    for ( size_t i = 0; i < _nodes.size(); ++i )
     {
       if ( !_nodes[i] ) continue;
       double len = totLen * _normPar[i];
@@ -4442,7 +4676,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
     if ( _edges[1] )
       l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
     
-    for ( unsigned i = 0; i < _nodes.size(); ++i )
+    for ( size_t i = 0; i < _nodes.size(); ++i )
     {
       if ( !_nodes[i] ) continue;
       double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
@@ -4461,7 +4695,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
 void _Shrinker1D::RestoreParams()
 {
   if ( _done )
-    for ( unsigned i = 0; i < _nodes.size(); ++i )
+    for ( size_t i = 0; i < _nodes.size(); ++i )
     {
       if ( !_nodes[i] ) continue;
       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
@@ -4514,7 +4748,7 @@ bool _ViscousBuilder::addBoundaryElements()
 {
   SMESH_MesherHelper helper( *_mesh );
 
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     _SolidData& data = _sdVec[i];
     TopTools_IndexedMapOfShape geomEdges;
@@ -4585,7 +4819,7 @@ bool _ViscousBuilder::addBoundaryElements()
         {
           const TopoDS_Shape* pF = fIt->next();
           if ( helper.IsSubShape( *pF, data._solid) &&
-               !_ignoreShapeIds.count( e2f->first ))
+               !data._ignoreFaceIds.count( e2f->first ))
             F = *pF;
         }
       }
@@ -4601,7 +4835,7 @@ bool _ViscousBuilder::addBoundaryElements()
       // Make faces
       const int dj1 = reverse ? 0 : 1;
       const int dj2 = reverse ? 1 : 0;
-      for ( unsigned j = 1; j < ledges.size(); ++j )
+      for ( size_t j = 1; j < ledges.size(); ++j )
       {
         vector< const SMDS_MeshNode*>&  nn1 = ledges[j-dj1]->_nodes;
         vector< const SMDS_MeshNode*>&  nn2 = ledges[j-dj2]->_nodes;