+//================================================================================
+/*!
+ * \brief Try to get normal at a singularity of a surface basing on it's nature
+ */
+//================================================================================
+
+bool _ViscousBuilder::getFaceNormalAtSingularity( const gp_XY& uv,
+ const TopoDS_Face& face,
+ SMESH_MesherHelper& helper,
+ gp_Dir& normal )
+{
+ BRepAdaptor_Surface surface( face );
+ gp_Dir axis;
+ if ( !getRovolutionAxis( surface, axis ))
+ return false;
+
+ double f,l, d, du, dv;
+ f = surface.FirstUParameter();
+ l = surface.LastUParameter();
+ d = ( uv.X() - f ) / ( l - f );
+ du = ( d < 0.5 ? +1. : -1 ) * 1e-5 * ( l - f );
+ f = surface.FirstVParameter();
+ l = surface.LastVParameter();
+ d = ( uv.Y() - f ) / ( l - f );
+ dv = ( d < 0.5 ? +1. : -1 ) * 1e-5 * ( l - f );
+
+ gp_Dir refDir;
+ gp_Pnt2d testUV = uv;
+ enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
+ double tol = 1e-5;
+ Handle(Geom_Surface) geomsurf = surface.Surface().Surface();
+ for ( int iLoop = 0; true ; ++iLoop )
+ {
+ testUV.SetCoord( testUV.X() + du, testUV.Y() + dv );
+ if ( GeomLib::NormEstim( geomsurf, testUV, tol, refDir ) == REGULAR )
+ break;
+ if ( iLoop > 20 )
+ return false;
+ tol /= 10.;
+ }
+
+ if ( axis * refDir < 0. )
+ axis.Reverse();
+
+ normal = axis;
+
+ return true;
+}
+