Salome HOME
- Modifing Geometry and Mesh Python scripts from SALOME 6 and before
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index a855c24a35a5a3182d9c9bcdc3f96b054ef9a581..6a6d871c43d65d25128c00a3f466aa903662e998 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2013  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
 #include "SMESH_Gen.hxx"
 #include "SMESH_Group.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 "utilities.h"
+#include "StdMeshers_FaceSide.hxx"
 
 #include <BRepAdaptor_Curve2d.hxx>
 #include <BRep_Tool.hxx>
@@ -224,8 +224,11 @@ namespace VISCOUS_3D
   struct _Simplex
   {
     const SMDS_MeshNode *_nPrev, *_nNext; // nodes on a smoothed mesh surface
-    _Simplex(const SMDS_MeshNode* nPrev=0, const SMDS_MeshNode* nNext=0)
-      : _nPrev(nPrev), _nNext(nNext) {}
+    const SMDS_MeshNode *_nOpp; // in 2D case, a node opposite to a smoothed node in QUAD
+    _Simplex(const SMDS_MeshNode* nPrev=0,
+             const SMDS_MeshNode* nNext=0,
+             const SMDS_MeshNode* nOpp=0)
+      : _nPrev(nPrev), _nNext(nNext), _nOpp(nOpp) {}
     bool IsForward(const SMDS_MeshNode* nSrc, const gp_XYZ* pntTgt) const
     {
       const double M[3][3] =
@@ -431,12 +434,18 @@ namespace VISCOUS_3D
     //vector<const SMDS_MeshNode*> _nodesAround;
     vector<_Simplex>             _simplices; // for quality check
 
+    enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR };
+
     bool Smooth(int&                  badNb,
                 Handle(Geom_Surface)& surface,
                 SMESH_MesherHelper&   helper,
                 const double          refSign,
-                bool                  isCentroidal,
+                SmoothType            how,
                 bool                  set3D);
+
+    gp_XY computeAngularPos(vector<gp_XY>& uv,
+                            const gp_XY&   uvToFix,
+                            const double   refSign );
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -566,20 +575,17 @@ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
 //
 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
   :SMESH_Hypothesis(hypId, studyId, gen),
-   _nbLayers(1), _thickness(1), _stretchFactor(1)
+   _isToIgnoreShapes(18), _nbLayers(1), _thickness(1), _stretchFactor(1)
 {
   _name = StdMeshers_ViscousLayers::GetHypType();
   _param_algo_dim = -3; // auxiliary hyp used by 3D algos
 } // --------------------------------------------------------------------------------
-void StdMeshers_ViscousLayers::SetBndShapesToIgnore(const std::vector<int>& faceIds)
+void StdMeshers_ViscousLayers::SetBndShapes(const std::vector<int>& faceIds, bool toIgnore)
 {
-  if ( faceIds != _ignoreBndShapeIds )
-    _ignoreBndShapeIds = faceIds, NotifySubMeshesHypothesisModification();
-} // --------------------------------------------------------------------------------
-bool StdMeshers_ViscousLayers::IsIgnoredShape(const int shapeID) const
-{
-  return ( find( _ignoreBndShapeIds.begin(), _ignoreBndShapeIds.end(), shapeID )
-           != _ignoreBndShapeIds.end() );
+  if ( faceIds != _shapeIds )
+    _shapeIds = faceIds, NotifySubMeshesHypothesisModification();
+  if ( _isToIgnoreShapes != toIgnore )
+    _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification();
 } // --------------------------------------------------------------------------------
 void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
 {
@@ -637,17 +643,22 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
   save << " " << _nbLayers
        << " " << _thickness
        << " " << _stretchFactor
-       << " " << _ignoreBndShapeIds.size();
-  for ( unsigned i = 0; i < _ignoreBndShapeIds.size(); ++i )
-    save << " " << _ignoreBndShapeIds[i];
+       << " " << _shapeIds.size();
+  for ( unsigned i = 0; i < _shapeIds.size(); ++i )
+    save << " " << _shapeIds[i];
+  save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies.
   return save;
 } // --------------------------------------------------------------------------------
 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
 {
-  int nbFaces, faceID;
+  int nbFaces, faceID, shapeToTreat;
   load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
-  while ( _ignoreBndShapeIds.size() < nbFaces && load >> faceID )
-    _ignoreBndShapeIds.push_back( faceID );
+  while ( _shapeIds.size() < nbFaces && load >> faceID )
+    _shapeIds.push_back( faceID );
+  if ( load >> shapeToTreat )
+    _isToIgnoreShapes = !shapeToTreat;
+  else
+    _isToIgnoreShapes = true; // old behavior
   return load;
 } // --------------------------------------------------------------------------------
 bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh*   theMesh,
@@ -770,13 +781,15 @@ 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 );
     for ( ; eExp.More(); eExp.Next() )
     {
       const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
-      if ( BRep_Tool::Degenerated( E )) continue;
+      if ( SMESH_Algo::isDegenerated( E )) continue;
       // check if 2D curve is concave
       BRepAdaptor_Curve2d curve( E, F );
       const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
@@ -788,10 +801,10 @@ namespace
         double u1 = intervals( i );
         double u2 = intervals( i+1 );
         curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 );
-        double cross = drv2 ^ drv1;
+        double cross = drv2 * drv1; //drv2 ^ drv1;
         if ( E.Orientation() == TopAbs_REVERSED )
           cross = -cross;
-        isConvex = ( cross 1e-9 );
+        isConvex = ( cross > -1e-9 );
       }
       // check if concavity is strong enough to care about it
       //const double maxAngle = 5 * Standard_PI180;
@@ -821,6 +834,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;
   }
   //--------------------------------------------------------------------------------
@@ -832,9 +865,11 @@ 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)
@@ -1071,7 +1106,7 @@ bool _ViscousBuilder::findFacesWithLayers()
   vector<TopoDS_Shape> ignoreFaces;
   for ( unsigned i = 0; i < _sdVec.size(); ++i )
   {
-    vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapesToIgnore();
+    vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
     for ( unsigned i = 0; i < ids.size(); ++i )
     {
       const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[i] );
@@ -1096,7 +1131,7 @@ bool _ViscousBuilder::findFacesWithLayers()
       {     
         _ignoreShapeIds.insert( faceInd );
         ignoreFaces.push_back( exp.Current() );
-        if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face( exp.Current() ), getMeshDS()))
+        if ( helper.IsReversedSubMesh( TopoDS::Face( exp.Current() )))
           _sdVec[i]._reversedFaceIds.insert( faceInd );
       }
     }
@@ -1961,9 +1996,10 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
     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 ));
     if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd ))
       std::swap( nPrev, nNext );
-    simplices.push_back( _Simplex( nPrev, nNext ));
+    simplices.push_back( _Simplex( nPrev, nNext, nOpp ));
   }
 
   if ( toSort )
@@ -2067,14 +2103,14 @@ bool _ViscousBuilder::inflate(_SolidData& data)
   // Limit inflation step size by geometry size found by itersecting
   // normals of _LayerEdge's with mesh faces
   double geomSize = Precision::Infinite(), intersecDist;
-  SMESH_MeshEditor editor( _mesh );
   auto_ptr<SMESH_ElementSearcher> searcher
-    ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) );
+    ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
+                                           data._proxyMesh->GetFaces( data._solid )) );
   for ( unsigned i = 0; i < data._edges.size(); ++i )
   {
     if ( data._edges[i]->IsOnEdge() ) continue;
     data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
-    if ( geomSize > intersecDist )
+    if ( geomSize > intersecDist && intersecDist > 0 )
       geomSize = intersecDist;
   }
   if ( data._stepSize > 0.3 * geomSize )
@@ -2265,9 +2301,9 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
   // Check if the last segments of _LayerEdge intersects 2D elements;
   // checked elements are either temporary faces or faces on surfaces w/o the layers
 
-  SMESH_MeshEditor editor( _mesh );
   auto_ptr<SMESH_ElementSearcher> searcher
-    ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) );
+    ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
+                                           data._proxyMesh->GetFaces( data._solid )) );
 
   distToIntersection = Precision::Infinite();
   double dist;
@@ -2522,6 +2558,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
@@ -2627,10 +2667,10 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   // 1) to find and fix intersection
   // 2) to check that no new intersection appears as result of 1)
 
-  SMESH_MeshEditor editor( _mesh );
   SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
                                                             tmpFaces.end()));
-  auto_ptr<SMESH_ElementSearcher> searcher ( editor.GetElementSearcher( fIt ));
+  auto_ptr<SMESH_ElementSearcher> searcher
+    ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt ));
 
   // 1) Find intersections
   double dist;
@@ -3550,7 +3590,7 @@ bool _ViscousBuilder::shrink()
         sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false);
       while ( subIt->more() )
       {
-        SMESH_subMesh* sub = subIt->next();
+        SMESH_subMesh*     sub = subIt->next();
         SMESHDS_SubMesh* subDS = sub->GetSubMeshDS();
         if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() ))
           continue;
@@ -3594,12 +3634,13 @@ bool _ViscousBuilder::shrink()
     vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
     {
       dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
+      const bool sortSimplices = isConcaveFace;
       for ( unsigned i = 0; i < smoothNodes.size(); ++i )
       {
         const SMDS_MeshNode* n = smoothNodes[i];
         nodesToSmooth[ i ]._node = n;
         // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
-        getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, isConcaveFace );
+        getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, sortSimplices );
         // fix up incorrect uv of nodes on the FACE
         helper.GetNodeUV( F, n, 0, &isOkUV);
         dumpMove( n );
@@ -3634,11 +3675,14 @@ bool _ViscousBuilder::shrink()
 
     bool shrinked = true;
     int badNb, shriStep=0, smooStep=0;
+    _SmoothNode::SmoothType smoothType
+      = isConcaveFace ? _SmoothNode::CENTROIDAL : _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 )
       {
@@ -3647,6 +3691,7 @@ bool _ViscousBuilder::shrink()
       dumpFunctionEnd();
 
       // Move nodes on EDGE's
+      // (XYZ is set as soon as a needed length reached in SetNewLength2d())
       set< _Shrinker1D* >::iterator shr = eShri1D.begin();
       for ( ; shr != eShri1D.end(); ++shr )
         (*shr)->Compute( /*set3D=*/false, helper );
@@ -3666,8 +3711,7 @@ bool _ViscousBuilder::shrink()
         for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
         {
           moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
-                                            /*isCentroidal=*/isConcaveFace,
-                                            /*set3D=*/isConcaveFace);
+                                            smoothType, /*set3D=*/isConcaveFace);
         }
         if ( badNb < oldBadNb )
           nbNoImpSteps = 0;
@@ -3678,33 +3722,29 @@ 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 );
     }
+
     // No wrongly shaped faces remain; final smooth. Set node XYZ.
-    // First, find out a needed quality of smoothing (high for quadrangles only)
-    bool highQuality;
+    bool isStructuredFixed = false;
+    if ( SMESH_2D_Algo* algo = dynamic_cast<SMESH_2D_Algo*>( sm->GetAlgo() ))
+      isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F );
+    if ( !isStructuredFixed )
     {
-      const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
-      if ( hasTria != hasQuad ) {
-        highQuality = hasQuad;
-      }
-      else {
-        set<int> nbNodesSet;
-        SMDS_ElemIteratorPtr fIt = smDS->GetElements();
-        while ( fIt->more() && nbNodesSet.size() < 2 )
-          nbNodesSet.insert( fIt->next()->NbCornerNodes() );
-        highQuality = ( *nbNodesSet.begin() == 4 );
+      if ( isConcaveFace )
+        fixBadFaces( F, helper ); // fix narrow faces by swapping diagonals
+      for ( int st = /*highQuality ? 10 :*/ 3; st; --st )
+      {
+        dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
+        for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
+        {
+          nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
+                                   smoothType,/*set3D=*/st==1 );
+        }
+        dumpFunctionEnd();
       }
     }
-    if ( !highQuality && isConcaveFace )
-      fixBadFaces( F, helper ); // fix narrow faces by swaping diagonals
-    for ( int st = highQuality ? 10 : 3; st; --st )
-    {
-      dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
-      for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
-        nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
-                                 /*isCentroidal=*/isConcaveFace,/*set3D=*/st==1 );
-      dumpFunctionEnd();
-    }
     // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
     VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
 
@@ -3752,6 +3792,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     double uvLen = uvDir.Magnitude();
     uvDir /= uvLen;
     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;
@@ -3781,7 +3822,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     }
 
     multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd;
-    const double minProj = p2n->first;
+    const double       minProj = p2n->first;
     const double projThreshold = 1.1 * uvLen;
     if ( minProj > projThreshold )
     {
@@ -3958,8 +3999,8 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
       const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide );
       const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 );
       trias [iSide].first  = badTrias[iTia];
-      trias [iSide].second = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, involvedFaces,
-                                                              & i1, & i2 );
+      trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces,
+                                                             & i1, & i2 );
       if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 )
         continue;
 
@@ -4057,11 +4098,13 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
           double proj = uvDirN * uvDir * kSafe;
           if ( proj < stepSize && proj > minStepSize )
             stepSize = proj;
+          else if ( proj < minStepSize )
+            stepSize = minStepSize;
         }
     }
 
     gp_Pnt2d newUV;
-    if ( stepSize == uvLen )
+    if ( uvLen - stepSize < _len / 20. )
     {
       newUV = tgtUV;
       _pos.clear();
@@ -4085,22 +4128,22 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
   {
     TopoDS_Edge E = TopoDS::Edge( _sWOL );
     const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
+    SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
 
     const double u2 = helper.GetNodeU( E, n2, tgtNode );
     const double uSrc   = _pos[0].Coord( U_SRC );
     const double lenTgt = _pos[0].Coord( LEN_TGT );
 
     double newU = _pos[0].Coord( U_TGT );
-    if ( lenTgt < 0.99 * fabs( uSrc-u2 ))
+    if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range
     {
       _pos.clear();
     }
     else
     {
-      newU = 0.1 * uSrc + 0.9 * u2;
+      newU = 0.1 * tgtPos->GetUParameter() + 0.9 * u2;
     }
-    SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
-    pos->SetUParameter( newU );
+    tgtPos->SetUParameter( newU );
 #ifdef __myDEBUG
     gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]);
     gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
@@ -4122,7 +4165,7 @@ bool _SmoothNode::Smooth(int&                  badNb,
                          Handle(Geom_Surface)& surface,
                          SMESH_MesherHelper&   helper,
                          const double          refSign,
-                         bool                  isCentroidal,
+                         SmoothType            how,
                          bool                  set3D)
 {
   const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
@@ -4134,7 +4177,30 @@ bool _SmoothNode::Smooth(int&                  badNb,
 
   // compute new UV for the node
   gp_XY newPos (0,0);
-  if ( isCentroidal && _simplices.size() > 3 )
+/*  if ( how == ANGULAR && _simplices.size() == 4 )
+  {
+    vector<gp_XY> corners; corners.reserve(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 );
+  }
+  else*/ if ( how == CENTROIDAL && _simplices.size() > 3 )
   {
     // average centers of diagonals wieghted with their reciprocal lengths
     if ( _simplices.size() == 4 )
@@ -4159,13 +4225,13 @@ bool _SmoothNode::Smooth(int&                  badNb,
           newPos += w * ( uv[i]+uv[i2] );
         }
       }
-      newPos /= 2 * sumWeight;
+      newPos /= 2 * sumWeight; // 2 is to get a middle between uv's
     }
   }
   else
   {
     // Laplacian smooth
-    isCentroidal = false;
+    //isCentroidal = false;
     for ( size_t i = 0; i < _simplices.size(); ++i )
       newPos += uv[i];
     newPos /= _simplices.size();
@@ -4207,6 +4273,66 @@ bool _SmoothNode::Smooth(int&                  badNb,
   return ( (tgtUV-newPos).SquareModulus() > 1e-10 );
 }
 
+//================================================================================
+/*!
+ * \brief Computes new UV using angle based smoothing technic
+ */
+//================================================================================
+
+gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
+                                     const gp_XY&   uvToFix,
+                                     const double   refSign)
+{
+  uv.push_back( uv.front() );
+
+  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];
+    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();
+  edgeSize.back() = edgeSize.front();
+
+  gp_XY newPos(0,0);
+  int nbEdges = 0;
+  double sumSize = 0;
+  for ( size_t i = 1; i < edgeDir.size(); ++i )
+  {
+    if ( edgeDir[i-1].X() > 1. ) continue;
+    int i1 = i-1;
+    while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() );
+    if ( i == edgeDir.size() ) break;
+    gp_XY p = uv[i];
+    gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() );
+    gp_XY norm2( -edgeDir[i].Y(),  edgeDir[i].X() );
+    gp_XY bisec = norm1 + norm2;
+    double bisecSize = bisec.Modulus();
+    if ( bisecSize < numeric_limits<double>::min() )
+    {
+      bisec = -edgeDir[i1] + edgeDir[i];
+      bisecSize = bisec.Modulus();
+    }
+    bisec /= bisecSize;
+
+    gp_XY  dirToN = uvToFix - p;
+    double distToN = dirToN.Modulus();
+    if ( bisec * dirToN < 0 )
+      distToN = -distToN;
+
+    newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] );
+    ++nbEdges;
+    sumSize += edgeSize[i1] + edgeSize[i];
+  }
+  newPos /= /*nbEdges * */sumSize;
+  return newPos;
+}
+
 //================================================================================
 /*!
  * \brief Delete _SolidData
@@ -4480,7 +4606,7 @@ bool _ViscousBuilder::addBoundaryElements()
         reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
         if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
           reverse = !reverse, F.Reverse();
-        if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face(F), getMeshDS() ))
+        if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
           reverse = !reverse;
       }
       else