Salome HOME
Minor doc changes
[modules/smesh.git] / src / SMESHUtils / SMESH_MeshAlgos.cxx
index d0d6008d373f426e51dbfc49d9891fef1ebfb276..e65e20114c42723fd8ec3ff435c115a73205332c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -6,7 +6,7 @@
 // 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
 
 #include "SMESH_MeshAlgos.hxx"
 
+#include "SMDS_FaceOfNodes.hxx"
 #include "SMDS_LinearEdge.hxx"
-#include "SMDS_VolumeTool.hxx"
 #include "SMDS_Mesh.hxx"
+#include "SMDS_PolygonalFaceOfNodes.hxx"
+#include "SMDS_VolumeTool.hxx"
 #include "SMESH_OctreeNode.hxx"
 
 #include <GC_MakeSegment.hxx>
@@ -58,7 +60,8 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
   /*!
    * \brief Constructor
    */
-  SMESH_NodeSearcherImpl( const SMDS_Mesh* theMesh )
+  SMESH_NodeSearcherImpl( const SMDS_Mesh*     theMesh   = 0,
+                          SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr() )
   {
     myMesh = ( SMDS_Mesh* ) theMesh;
 
@@ -68,6 +71,14 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
       while ( nIt->more() )
         nodes.insert( nodes.end(), nIt->next() );
     }
+    else if ( theElemIt )
+    {
+      while ( theElemIt->more() )
+      {
+        const SMDS_MeshElement* e = theElemIt->next();
+        nodes.insert( e->begin_nodes(), e->end_nodes() );
+      }
+    }
     myOctreeNode = new SMESH_OctreeNode(nodes) ;
 
     // get max size of a leaf box
@@ -166,6 +177,18 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
     return closestNode;
   }
 
+  //---------------------------------------------------------------------
+  /*!
+   * \brief Finds nodes located within a tolerance near a point 
+   */
+  int FindNearPoint(const gp_Pnt&                        point,
+                    const double                         tolerance,
+                    std::vector< const SMDS_MeshNode* >& foundNodes)
+  {
+    myOctreeNode->NodesAround( point.Coord(), foundNodes, tolerance );
+    return foundNodes.size();
+  }
+
   //---------------------------------------------------------------------
   /*!
    * \brief Destructor
@@ -210,7 +233,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     void getElementsInSphere ( const gp_XYZ& center,
                                const double  radius, TIDSortedElemSet& foundElems);
     size_t getSize() { return std::max( _size, _elements.size() ); }
-    ~ElementBndBoxTree();
+    virtual ~ElementBndBoxTree();
 
   protected:
     ElementBndBoxTree():_size(0) {}
@@ -256,7 +279,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 
   ElementBndBoxTree::~ElementBndBoxTree()
   {
-    for ( int i = 0; i < _elements.size(); ++i )
+    for ( size_t i = 0; i < _elements.size(); ++i )
       if ( --_elements[i]->_refCount <= 0 )
         delete _elements[i];
   }
@@ -270,7 +293,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
   Bnd_B3d* ElementBndBoxTree::buildRootBox()
   {
     Bnd_B3d* box = new Bnd_B3d;
-    for ( int i = 0; i < _elements.size(); ++i )
+    for ( size_t i = 0; i < _elements.size(); ++i )
       box->Add( *_elements[i] );
     return box;
   }
@@ -283,7 +306,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 
   void ElementBndBoxTree::buildChildrenData()
   {
-    for ( int i = 0; i < _elements.size(); ++i )
+    for ( size_t i = 0; i < _elements.size(); ++i )
     {
       for (int j = 0; j < 8; j++)
       {
@@ -301,7 +324,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     for (int j = 0; j < 8; j++)
     {
       ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j]);
-      if ( child->_elements.size() <= MaxNbElemsInLeaf )
+      if ((int) child->_elements.size() <= MaxNbElemsInLeaf )
         child->myIsLeaf = true;
 
       if ( child->_elements.capacity() - child->_elements.size() > 1000 )
@@ -323,7 +346,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 
     if ( isLeaf() )
     {
-      for ( int i = 0; i < _elements.size(); ++i )
+      for ( size_t i = 0; i < _elements.size(); ++i )
         if ( !_elements[i]->IsOut( point.XYZ() ))
           foundElems.insert( _elements[i]->_element );
     }
@@ -348,7 +371,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 
     if ( isLeaf() )
     {
-      for ( int i = 0; i < _elements.size(); ++i )
+      for ( size_t i = 0; i < _elements.size(); ++i )
         if ( !_elements[i]->IsOut( line ))
           foundElems.insert( _elements[i]->_element );
     }
@@ -374,7 +397,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 
     if ( isLeaf() )
     {
-      for ( int i = 0; i < _elements.size(); ++i )
+      for ( size_t i = 0; i < _elements.size(); ++i )
         if ( !_elements[i]->IsOut( center, radius ))
           foundElems.insert( _elements[i]->_element );
     }
@@ -410,6 +433,10 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
  */
 //=======================================================================
 
+SMESH_ElementSearcher::~SMESH_ElementSearcher()
+{
+}
+
 struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
 {
   SMDS_Mesh*                   _mesh;
@@ -421,9 +448,11 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   bool                         _outerFacesFound;
   set<const SMDS_MeshElement*> _outerFaces; // empty means "no internal faces at all"
 
-  SMESH_ElementSearcherImpl( SMDS_Mesh& mesh, SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr())
-    : _mesh(&mesh),_meshPartIt(elemIt),_ebbTree(0),_nodeSearcher(0),_tolerance(-1),_outerFacesFound(false) {}
-  ~SMESH_ElementSearcherImpl()
+  SMESH_ElementSearcherImpl( SMDS_Mesh&           mesh,
+                             double               tol=-1,
+                             SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr())
+    : _mesh(&mesh),_meshPartIt(elemIt),_ebbTree(0),_nodeSearcher(0),_tolerance(tol),_outerFacesFound(false) {}
+  virtual ~SMESH_ElementSearcherImpl()
   {
     if ( _ebbTree )      delete _ebbTree;      _ebbTree      = 0;
     if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0;
@@ -438,6 +467,10 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   void GetElementsNearLine( const gp_Ax1&                      line,
                             SMDSAbs_ElementType                type,
                             vector< const SMDS_MeshElement* >& foundElems);
+  void GetElementsInSphere( const gp_XYZ&                      center,
+                            const double                       radius,
+                            SMDSAbs_ElementType                type,
+                            vector< const SMDS_MeshElement* >& foundElems);
   double getTolerance();
   bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
                             const double tolerance, double & param);
@@ -446,6 +479,7 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   {
     return _outerFaces.empty() || _outerFaces.count(face);
   }
+
   struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
   {
     const SMDS_MeshElement* _face;
@@ -550,7 +584,7 @@ bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin&           lin
   {
     GC_MakeSegment edge( SMESH_TNodeXYZ( face->GetNode( i )),
                          SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) ));
-    anExtCC.Init( lineCurve, edge);
+    anExtCC.Init( lineCurve, edge.Value() );
     if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol)
     {
       Quantity_Parameter pl, pe;
@@ -626,6 +660,7 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF
         set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin();
         for ( ; face != faces.end(); ++face )
         {
+          if ( *face == outerFace ) continue;
           if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false ))
             continue;
           gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2;
@@ -640,8 +675,8 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF
     // store the found outer face and add its links to continue seaching from
     if ( outerFace2 )
     {
-      _outerFaces.insert( outerFace );
-      int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 );
+      _outerFaces.insert( outerFace2 );
+      int nbNodes = outerFace2->NbCornerNodes();
       for ( int i = 0; i < nbNodes; ++i )
       {
         SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes));
@@ -688,23 +723,27 @@ FindElementsByPoint(const gp_Pnt&                      point,
   if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement || type == SMDSAbs_Ball)
   {
     if ( !_nodeSearcher )
-      _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
-
-    const SMDS_MeshNode* closeNode = _nodeSearcher->FindClosestTo( point );
-    if ( !closeNode ) return foundElements.size();
-
-    if ( point.Distance( SMESH_TNodeXYZ( closeNode )) > tolerance )
-      return foundElements.size(); // to far from any node
+    {
+      if ( _meshPartIt )
+        _nodeSearcher = new SMESH_NodeSearcherImpl( 0, _meshPartIt );
+      else
+        _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
+    }
+    std::vector< const SMDS_MeshNode* > foundNodes;
+    _nodeSearcher->FindNearPoint( point, tolerance, foundNodes );
 
     if ( type == SMDSAbs_Node )
     {
-      foundElements.push_back( closeNode );
+      foundElements.assign( foundNodes.begin(), foundNodes.end() );
     }
     else
     {
-      SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( type );
-      while ( elemIt->more() )
-        foundElements.push_back( elemIt->next() );
+      for ( size_t i = 0; i < foundNodes.size(); ++i )
+      {
+        SMDS_ElemIteratorPtr elemIt = foundNodes[i]->GetInverseElementIterator( type );
+        while ( elemIt->more() )
+          foundElements.push_back( elemIt->next() );
+      }
     }
   }
   // =================================================================================
@@ -739,7 +778,7 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
 {
   const SMDS_MeshElement* closestElem = 0;
 
-  if ( type == SMDSAbs_Face )
+  if ( type == SMDSAbs_Face || type == SMDSAbs_Volume )
   {
     if ( !_ebbTree || _elementType != type )
     {
@@ -753,10 +792,10 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
     {
       gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox()->CornerMin() +
                                  _ebbTree->getBox()->CornerMax() );
-      double radius;
+      double radius = -1;
       if ( _ebbTree->getBox()->IsOut( point.XYZ() ))
         radius = point.Distance( boxCenter ) - 0.5 * _ebbTree->maxSize();
-      else
+      if ( radius < 0 )
         radius = _ebbTree->maxSize() / pow( 2., _ebbTree->getHeight()) / 2;
       while ( suspectElems.empty() )
       {
@@ -769,8 +808,7 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
     TIDSortedElemSet::iterator elem = suspectElems.begin();
     for ( ; elem != suspectElems.end(); ++elem )
     {
-      double dist = SMESH_MeshAlgos::GetDistance( dynamic_cast<const SMDS_MeshFace*>(*elem),
-                                                   point );
+      double dist = SMESH_MeshAlgos::GetDistance( *elem, point );
       if ( dist < minDist + 1e-10)
       {
         minDist = dist;
@@ -868,8 +906,9 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
       }
       else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 )
       {
+        double tol = 1e-4 * Sqrt( fNorm.Modulus() );
         gp_Pnt intersectionPoint = intersection.Point(1);
-        if ( !SMESH_MeshAlgos::IsOut( *face, intersectionPoint, tolerance ))
+        if ( !SMESH_MeshAlgos::IsOut( *face, intersectionPoint, tol ))
           u2inters.insert(make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm )));
       }
     }
@@ -1059,6 +1098,27 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1&
   foundElems.assign( suspectFaces.begin(), suspectFaces.end());
 }
 
+//=======================================================================
+/*
+ * Return elements whose bounding box intersects a sphere
+ */
+//=======================================================================
+
+void SMESH_ElementSearcherImpl::GetElementsInSphere( const gp_XYZ&                      center,
+                                                     const double                       radius,
+                                                     SMDSAbs_ElementType                type,
+                                                     vector< const SMDS_MeshElement* >& foundElems)
+{
+  if ( !_ebbTree || _elementType != type )
+  {
+    if ( _ebbTree ) delete _ebbTree;
+    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
+  }
+  TIDSortedElemSet suspectFaces; // elements possibly intersecting the line
+  _ebbTree->getElementsInSphere( center, radius, suspectFaces );
+  foundElems.assign( suspectFaces.begin(), suspectFaces.end() );
+}
+
 //=======================================================================
 /*!
  * \brief Return true if the point is IN or ON of the element
@@ -1074,32 +1134,19 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
 
   // get ordered nodes
 
-  vector< gp_XYZ > xyz;
-  vector<const SMDS_MeshNode*> nodeList;
+  vector< SMESH_TNodeXYZ > xyz; xyz.reserve( element->NbNodes()+1 );
 
-  SMDS_ElemIteratorPtr nodeIt = element->nodesIterator();
-  if ( element->IsQuadratic() ) {
-    nodeIt = element->interlacedNodesElemIterator();
-    // if (const SMDS_VtkFace* f=dynamic_cast<const SMDS_VtkFace*>(element))
-    //   nodeIt = f->interlacedNodesElemIterator();
-    // else if (const SMDS_VtkEdge*  e =dynamic_cast<const SMDS_VtkEdge*>(element))
-    //   nodeIt = e->interlacedNodesElemIterator();
-  }
-  while ( nodeIt->more() )
-  {
-    SMESH_TNodeXYZ node = nodeIt->next();
-    xyz.push_back( node );
-    nodeList.push_back(node._node);
-  }
+  SMDS_ElemIteratorPtr nodeIt = element->interlacedNodesElemIterator();
+  for ( int i = 0; nodeIt->more(); ++i )
+    xyz.push_back( SMESH_TNodeXYZ( nodeIt->next() ));
 
-  int i, nbNodes = (int) nodeList.size(); // central node of biquadratic is missing
+  int i, nbNodes = (int) xyz.size(); // central node of biquadratic is missing
 
   if ( element->GetType() == SMDSAbs_Face ) // --------------------------------------------------
   {
     // compute face normal
     gp_Vec faceNorm(0,0,0);
     xyz.push_back( xyz.front() );
-    nodeList.push_back( nodeList.front() );
     for ( i = 0; i < nbNodes; ++i )
     {
       gp_Vec edge1( xyz[i+1], xyz[i]);
@@ -1112,7 +1159,7 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
       // degenerated face: point is out if it is out of all face edges
       for ( i = 0; i < nbNodes; ++i )
       {
-        SMDS_LinearEdge edge( nodeList[i], nodeList[i+1] );
+        SMDS_LinearEdge edge( xyz[i]._node, xyz[i+1]._node );
         if ( !IsOut( &edge, point, tol ))
           return false;
       }
@@ -1122,8 +1169,23 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
 
     // check if the point lays on face plane
     gp_Vec n2p( xyz[0], point );
-    if ( fabs( n2p * faceNorm ) > tol )
-      return true; // not on face plane
+    double dot = n2p * faceNorm;
+    if ( Abs( dot ) > tol ) // not on face plane
+    {
+      bool isOut = true;
+      if ( nbNodes > 3 ) // maybe the face is not planar
+      {
+        double elemThick = 0;
+        for ( i = 1; i < nbNodes; ++i )
+        {
+          gp_Vec n2n( xyz[0], xyz[i] );
+          elemThick = Max( elemThick, Abs( n2n * faceNorm ));
+        }
+        isOut = Abs( dot ) > elemThick + tol;
+      }
+      if ( isOut )
+        return isOut;
+    }
 
     // check if point is out of face boundary:
     // define it by closest transition of a ray point->infinity through face boundary
@@ -1135,7 +1197,7 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     normSize = plnNorm.Magnitude();
     if ( normSize <= tol ) return false; // point coincides with the first node
     plnNorm /= normSize;
-    // for each node of the face, compute its signed distance to the plane
+    // for each node of the face, compute its signed distance to the cutting plane
     vector<double> dist( nbNodes + 1);
     for ( i = 0; i < nbNodes; ++i )
     {
@@ -1145,12 +1207,12 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     dist.back() = dist.front();
     // find the closest intersection
     int    iClosest = -1;
-    double rClosest, distClosest = 1e100;;
+    double rClosest = 0, distClosest = 1e100;
     gp_Pnt pClosest;
     for ( i = 0; i < nbNodes; ++i )
     {
       double r;
-      if ( fabs( dist[i]) < tol )
+      if ( fabs( dist[i] ) < tol )
         r = 0.;
       else if ( fabs( dist[i+1]) < tol )
         r = 1.;
@@ -1159,17 +1221,14 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
       else
         continue; // no intersection
       gp_Pnt pInt = xyz[i] * (1.-r) + xyz[i+1] * r;
-      gp_Vec p2int ( point, pInt);
-      if ( p2int * ray > -tol ) // right half-space
+      gp_Vec p2int( point, pInt);
+      double intDist = p2int.SquareMagnitude();
+      if ( intDist < distClosest )
       {
-        double intDist = p2int.SquareMagnitude();
-        if ( intDist < distClosest )
-        {
-          iClosest = i;
-          rClosest = r;
-          pClosest = pInt;
-          distClosest = intDist;
-        }
+        iClosest = i;
+        rClosest = r;
+        pClosest = pInt;
+        distClosest = intDist;
       }
     }
     if ( iClosest < 0 )
@@ -1193,28 +1252,39 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     bool covexCorner = ( edgeNorm * edgeAdjacent * (rClosest==1. ? 1. : -1.)) < 0;
     return covexCorner ? (out || out2) : (out && out2);
   }
+
   if ( element->GetType() == SMDSAbs_Edge ) // --------------------------------------------------
   {
     // point is out of edge if it is NOT ON any straight part of edge
     // (we consider quadratic edge as being composed of two straight parts)
     for ( i = 1; i < nbNodes; ++i )
     {
-      gp_Vec edge( xyz[i-1], xyz[i]);
-      gp_Vec n1p ( xyz[i-1], point);
-      double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude();
-      if ( dist > tol )
+      gp_Vec edge( xyz[i-1], xyz[i] );
+      gp_Vec n1p ( xyz[i-1], point  );
+      double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
+      if ( u <= 0. ) {
+        if ( n1p.SquareMagnitude() < tol * tol )
+          return false;
         continue;
-      gp_Vec n2p( xyz[i], point );
-      if ( fabs( edge.Magnitude() - n1p.Magnitude() - n2p.Magnitude()) > tol )
+      }
+      if ( u >= 1. ) {
+        if ( point.SquareDistance( xyz[i] ) < tol * tol )
+          return false;
+        continue;
+      }
+      gp_XYZ proj = ( 1. - u ) * xyz[i-1] + u * xyz[i]; // projection of the point on the edge
+      double dist2 = point.SquareDistance( proj );
+      if ( dist2 > tol * tol )
         continue;
       return false; // point is ON this part
     }
     return true;
   }
+
   // Node or 0D element -------------------------------------------------------------------------
   {
     gp_Vec n2p ( xyz[0], point );
-    return n2p.Magnitude() <= tol;
+    return n2p.SquareMagnitude() > tol * tol;
   }
   return true;
 }
@@ -1286,6 +1356,32 @@ namespace
   }
 }
 
+//=======================================================================
+/*!
+ * \brief Return minimal distance from a point to an element
+ *
+ * Currently we ignore non-planarity and 2nd order of face
+ */
+//=======================================================================
+
+double SMESH_MeshAlgos::GetDistance( const SMDS_MeshElement* elem,
+                                     const gp_Pnt&           point )
+{
+  switch ( elem->GetType() )
+  {
+  case SMDSAbs_Volume:
+    return GetDistance( dynamic_cast<const SMDS_MeshVolume*>( elem ), point);
+  case SMDSAbs_Face:
+    return GetDistance( dynamic_cast<const SMDS_MeshFace*>( elem ), point);
+  case SMDSAbs_Edge:
+    return GetDistance( dynamic_cast<const SMDS_MeshEdge*>( elem ), point);
+  case SMDSAbs_Node:
+    return point.Distance( SMESH_TNodeXYZ( elem ));
+  default:;
+  }
+  return -1;
+}
+
 //=======================================================================
 /*!
  * \brief Return minimal distance from a point to a face
@@ -1378,10 +1474,132 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face,
     // cout << distVec.Magnitude()  << " VERTEX " << face->GetNode(pos._index)->GetID() << endl;
     return distVec.Magnitude();
   }
+  default:;
   }
   return badDistance;
 }
 
+//=======================================================================
+/*!
+ * \brief Return minimal distance from a point to an edge
+ */
+//=======================================================================
+
+double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg, const gp_Pnt& point )
+{
+  double dist = Precision::Infinite();
+  if ( !seg ) return dist;
+
+  int i = 0, nbNodes = seg->NbNodes();
+
+  vector< SMESH_TNodeXYZ > xyz( nbNodes );
+  SMDS_ElemIteratorPtr nodeIt = seg->interlacedNodesElemIterator();
+  while ( nodeIt->more() )
+    xyz[ i++ ].Set( nodeIt->next() );
+
+  for ( i = 1; i < nbNodes; ++i )
+  {
+    gp_Vec edge( xyz[i-1], xyz[i] );
+    gp_Vec n1p ( xyz[i-1], point  );
+    double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
+    if ( u <= 0. ) {
+      dist = Min( dist, n1p.SquareMagnitude() );
+    }
+    else if ( u >= 1. ) {
+      dist = Min( dist, point.SquareDistance( xyz[i] ));
+    }
+    else {
+      gp_XYZ proj = ( 1. - u ) * xyz[i-1] + u * xyz[i]; // projection of the point on the edge
+      dist = Min( dist, point.SquareDistance( proj ));
+    }
+  }
+  return Sqrt( dist );
+}
+
+//=======================================================================
+/*!
+ * \brief Return minimal distance from a point to a volume
+ *
+ * Currently we ignore non-planarity and 2nd order
+ */
+//=======================================================================
+
+double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt& point )
+{
+  SMDS_VolumeTool vTool( volume );
+  vTool.SetExternalNormal();
+  const int iQ = volume->IsQuadratic() ? 2 : 1;
+
+  double n[3], bc[3];
+  double minDist = 1e100, dist;
+  for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+  {
+    // skip a facet with normal not "looking at" the point
+    if ( !vTool.GetFaceNormal( iF, n[0], n[1], n[2] ) ||
+         !vTool.GetFaceBaryCenter( iF, bc[0], bc[1], bc[2] ))
+      continue;
+    gp_XYZ bcp = point.XYZ() - gp_XYZ( bc[0], bc[1], bc[2] );
+    if ( gp_XYZ( n[0], n[1], n[2] ) * bcp < 1e-6 )
+      continue;
+
+    // find distance to a facet
+    const SMDS_MeshNode** nodes = vTool.GetFaceNodes( iF );
+    switch ( vTool.NbFaceNodes( iF ) / iQ ) {
+    case 3:
+    {
+      SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ] );
+      dist = GetDistance( &tmpFace, point );
+      break;
+    }
+    case 4:
+    {
+      SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ], nodes[ 3*iQ ]);
+      dist = GetDistance( &tmpFace, point );
+      break;
+    }
+    default:
+      vector<const SMDS_MeshNode *> nvec( nodes, nodes + vTool.NbFaceNodes( iF ));
+      SMDS_PolygonalFaceOfNodes tmpFace( nvec );
+      dist = GetDistance( &tmpFace, point );
+    }
+    minDist = Min( minDist, dist );
+  }
+  return minDist;
+}
+
+//================================================================================
+/*!
+ * \brief Returns barycentric coordinates of a point within a triangle.
+ *        A not returned bc2 = 1. - bc0 - bc1.
+ *        The point lies within the triangle if ( bc0 >= 0 && bc1 >= 0 && bc0+bc1 <= 1 )
+ */
+//================================================================================
+
+void SMESH_MeshAlgos::GetBarycentricCoords( const gp_XY& p,
+                                            const gp_XY& t0,
+                                            const gp_XY& t1,
+                                            const gp_XY& t2,
+                                            double &     bc0,
+                                            double &     bc1)
+{
+  const double // matrix 2x2
+    T11 = t0.X()-t2.X(), T12 = t1.X()-t2.X(),
+    T21 = t0.Y()-t2.Y(), T22 = t1.Y()-t2.Y();
+  const double Tdet = T11*T22 - T12*T21; // matrix determinant
+  if ( Abs( Tdet ) < std::numeric_limits<double>::min() )
+  {
+    bc0 = bc1 = 2.;
+    return;
+  }
+  // matrix inverse
+  const double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
+  // vector
+  const double r11 = p.X()-t2.X(), r12 = p.Y()-t2.Y();
+  // barycentric coordinates: mutiply matrix by vector
+  bc0 = (t11 * r11 + t12 * r12)/Tdet;
+  bc1 = (t21 * r11 + t22 * r12)/Tdet;
+}
+
 //=======================================================================
 //function : FindFaceInSet
 //purpose  : Return a face having linked nodes n1 and n2 and which is
@@ -1399,14 +1617,12 @@ SMESH_MeshAlgos::FindFaceInSet(const SMDS_MeshNode*    n1,
                                int*                    n2ind)
 
 {
-  int i1, i2;
+  int i1 = 0, i2 = 0;
   const SMDS_MeshElement* face = 0;
 
   SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
-  //MESSAGE("n1->GetInverseElementIterator(SMDSAbs_Face) " << invElemIt);
   while ( invElemIt->more() && !face ) // loop on inverse faces of n1
   {
-    //MESSAGE("in while ( invElemIt->more() && !face )");
     const SMDS_MeshElement* elem = invElemIt->next();
     if (avoidSet.count( elem ))
       continue;
@@ -1425,9 +1641,6 @@ SMESH_MeshAlgos::FindFaceInSet(const SMDS_MeshNode*    n1,
     if ( !face && elem->IsQuadratic())
     {
       // analysis for quadratic elements using all nodes
-      // const SMDS_VtkFace* F = dynamic_cast<const SMDS_VtkFace*>(elem);
-      // if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
-      // use special nodes iterator
       SMDS_ElemIteratorPtr anIter = elem->interlacedNodesElemIterator();
       const SMDS_MeshNode* prevN = static_cast<const SMDS_MeshNode*>( anIter->next() );
       for ( i1 = -1, i2 = 0; anIter->more() && !face; i1++, i2++ )
@@ -1462,7 +1675,7 @@ bool SMESH_MeshAlgos::FaceNormal(const SMDS_MeshElement* F, gp_XYZ& normal, bool
     return false;
 
   normal.SetCoord(0,0,0);
-  int nbNodes = F->IsQuadratic() ? F->NbNodes()/2 : F->NbNodes();
+  int nbNodes = F->NbCornerNodes();
   for ( int i = 0; i < nbNodes-2; ++i )
   {
     gp_XYZ p[3];
@@ -1507,15 +1720,27 @@ SMESH_NodeSearcher* SMESH_MeshAlgos::GetNodeSearcher(SMDS_Mesh& mesh)
   return new SMESH_NodeSearcherImpl( &mesh );
 }
 
+//=======================================================================
+/*!
+ * \brief Return SMESH_NodeSearcher
+ */
+//=======================================================================
+
+SMESH_NodeSearcher* SMESH_MeshAlgos::GetNodeSearcher(SMDS_ElemIteratorPtr elemIt)
+{
+  return new SMESH_NodeSearcherImpl( 0, elemIt );
+}
+
 //=======================================================================
 /*!
  * \brief Return SMESH_ElementSearcher
  */
 //=======================================================================
 
-SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh& mesh)
+SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh& mesh,
+                                                           double     tolerance)
 {
-  return new SMESH_ElementSearcherImpl( mesh );
+  return new SMESH_ElementSearcherImpl( mesh, tolerance );
 }
 
 //=======================================================================
@@ -1525,7 +1750,8 @@ SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh& mesh)
 //=======================================================================
 
 SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh&           mesh,
-                                                           SMDS_ElemIteratorPtr elemIt)
+                                                           SMDS_ElemIteratorPtr elemIt,
+                                                           double               tolerance)
 {
-  return new SMESH_ElementSearcherImpl( mesh, elemIt );
+  return new SMESH_ElementSearcherImpl( mesh, tolerance, elemIt );
 }