GPUSPHGUI: add FillHole() operation

which is needed for NETGEN 2D remesher

idl/SMESH_MeshEditor.idl

@@ -754,6 +754,30 @@ module SMESH
     short GetPointState(in double x, in double y, in double z)
       raises (SALOME::SALOME_Exception);
 
+    /*!
+     * Check if a 2D mesh is manifold
+     */
+    boolean IsManifold()
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Check if orientation of 2D elements is coherent
+     */
+    boolean IsCoherentOrientation2D()
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Returns all or only closed FreeBorder's.
+     */
+    ListOfFreeBorders FindFreeBorders(in boolean closedOnly)
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Fill with 2D elements a hole defined by a FreeBorder.
+     */
+    void FillHole(in FreeBorder hole)
+      raises (SALOME::SALOME_Exception);
+
     /*!
      * Returns groups of FreeBorder's coincident within the given tolerance.
      * If the tolerance <= 0.0 then one tenth of an average size of elements adjacent

src/SMESHUtils/CMakeLists.txt

@@ -47,7 +47,7 @@ SET(_link_LIBRARIES
   ${CAS_TKMesh}
   ${Boost_LIBRARIES}
   SMDS
-)
+  )
 
 # --- headers ---
 
@@ -68,7 +68,7 @@ SET(SMESHUtils_HEADERS
   SMESH_MAT2d.hxx
   SMESH_ControlPnt.hxx
   SMESH_Delaunay.hxx
-)
+  )
 
 # --- sources ---
 
@@ -86,6 +86,7 @@ SET(SMESHUtils_SOURCES
   SMESH_ControlPnt.cxx
   SMESH_DeMerge.cxx
   SMESH_Delaunay.cxx
+  SMESH_FillHole.cxx
   )
 
 # --- rules ---

src/SMESHUtils/SMESH_FillHole.cxx

@@ -0,0 +1,500 @@
+// 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
+//
+// 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, 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
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : SMESH_FillHole.cxx
+// Created   : Tue Sep 26 15:11:17 2017
+// Author    : Edward AGAPOV (eap)
+//
+
+#include "SMESH_MeshAlgos.hxx"
+
+#include "SMESH_Comment.hxx"
+
+#include "SMESH_TypeDefs.hxx"
+#include "SMDS_Mesh.hxx"
+
+#include <Utils_SALOME_Exception.hxx>
+
+#include <boost/intrusive/circular_list_algorithms.hpp>
+#include <boost/container/flat_map.hpp>
+
+#include <Bnd_B3d.hxx>
+
+namespace
+{
+  bool isSmallAngle( double cos2 )
+  {
+    // cosine of min angle at which adjacent faces are considered overlapping
+    const double theMinCos2 = 0.996 * 0.996; // ~5 degrees
+    return ( cos2 > theMinCos2 );
+  }
+
+  struct BEdge;
+  typedef std::multimap< double, BEdge* >          TAngleMap;
+  typedef std::map< const SMDS_MeshElement*, int > TFaceIndMap;
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Edge of a free border
+   */
+  struct BEdge
+  {
+    const SMDS_MeshNode*    myNode1;
+    const SMDS_MeshNode*    myNode2;
+    const SMDS_MeshElement* myFace;   // face adjacent to the border
+
+    gp_XYZ                  myFaceNorm;
+    gp_XYZ                  myDir;     // myNode1 -> myNode2
+    double                  myDirCoef; // 1. or -1, to make myDir oriented as myNodes in myFace
+    double                  myLength;  // between nodes
+    double                  myAngleWithPrev; // between myDir and -myPrev->myDir
+    TAngleMap::iterator     myAngleMapPos;
+    double                  myOverlapAngle;  // angle delta due to overlapping
+    const SMDS_MeshNode*    myNode1Shift;    // nodes created to avoid overlapping of faces
+    const SMDS_MeshNode*    myNode2Shift;
+
+    BEdge*                  myPrev; // neighbors in the border
+    BEdge*                  myNext;
+
+    BEdge(): myNode1Shift(0), myNode2Shift(0) {}
+    void   Init( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2,
+                 const SMDS_MeshElement* f=0,
+                 const SMDS_MeshNode* nf1=0, const SMDS_MeshNode* nf2=0 );
+    void   ComputeAngle( bool reverseAngle = false );
+    void   ShiftOverlapped( const SMDS_MeshNode*                  oppNode,
+                            const TFaceIndMap&                    capFaceWithBordInd,
+                            SMDS_Mesh&                            mesh,
+                            std::vector<const SMDS_MeshElement*>& newFaces);
+    void   MakeShiftfFaces( SMDS_Mesh&                            mesh,
+                            std::vector<const SMDS_MeshElement*>& newFaces,
+                            const bool                            isReverse );
+    gp_XYZ GetInFaceDir() const { return myFaceNorm ^ myDir * myDirCoef; }
+    void   InsertSelf(TAngleMap& edgesByAngle, bool isReverseFaces, bool reBind, bool useOverlap )
+    {
+      if ( reBind ) edgesByAngle.erase( myAngleMapPos );
+      double key = (( isReverseFaces ? 2 * M_PI - myAngleWithPrev : myAngleWithPrev )
+                    + myOverlapAngle * useOverlap );
+      myAngleMapPos = edgesByAngle.insert( std::make_pair( key, this ));
+    }
+
+    // traits used by boost::intrusive::circular_list_algorithms
+    typedef BEdge         node;
+    typedef BEdge *       node_ptr;
+    typedef const BEdge * const_node_ptr;
+    static node_ptr get_next(const_node_ptr n)             {  return n->myNext;  }
+    static void     set_next(node_ptr n, node_ptr next)    {  n->myNext = next;  }
+    static node_ptr get_previous(const_node_ptr n)         {  return n->myPrev;  }
+    static void     set_previous(node_ptr n, node_ptr prev){  n->myPrev = prev;  }
+  };
+
+  //================================================================================
+  /*!
+   * \brief Initialize a border edge data
+   */
+  //================================================================================
+
+  void BEdge::Init( const SMDS_MeshNode*    n1,
+                    const SMDS_MeshNode*    n2,
+                    const SMDS_MeshElement* newFace, // new cap face
+                    const SMDS_MeshNode*    nf1,
+                    const SMDS_MeshNode*    nf2 )
+  {
+    myNode1  = n1;
+    myNode2  = n2;
+    myDir    = SMESH_NodeXYZ( n2 ) - SMESH_NodeXYZ( n1 );
+    myLength = myDir.Modulus();
+    if ( myLength > std::numeric_limits<double>::min() )
+      myDir /= myLength;
+
+    myFace = newFace;
+    if ( !myFace )
+    {
+      TIDSortedElemSet elemSet, avoidSet;
+      int ind1, ind2;
+      myFace = SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet, &ind1, &ind2 );
+      if ( !myFace )
+        throw SALOME_Exception( SMESH_Comment("No face sharing nodes #")
+                                << myNode1->GetID() << " and #" << myNode2->GetID());
+      avoidSet.insert( myFace );
+      if ( SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet ))
+        throw SALOME_Exception( SMESH_Comment("No free border between nodes #")
+                                << myNode1->GetID() << " and #" << myNode2->GetID());
+
+      myDirCoef = SMESH_MeshAlgos::IsRightOrder( myFace, myNode1, myNode2 ) ? 1. : -1.;
+    }
+
+    if (! SMESH_MeshAlgos::FaceNormal( myFace, myFaceNorm, /*normalized=*/false ))
+    {
+      SMDS_ElemIteratorPtr fIt = myNode1->GetInverseElementIterator( SMDSAbs_Face );
+      while ( fIt->more() )
+        if ( SMESH_MeshAlgos::FaceNormal( fIt->next(), myFaceNorm, /*normalized=*/false ))
+          break;
+    }
+
+    if ( newFace )
+    {
+      myFace    = 0;
+      myDirCoef = SMESH_MeshAlgos::IsRightOrder( newFace, nf1, nf2 ) ? 1. : -1.;
+      if ( myPrev->myNode2 == n1 )
+        myNode1Shift = myPrev->myNode2Shift;
+      if ( myNext->myNode1 == n2 )
+        myNode2Shift = myNext->myNode1Shift;
+    }
+    else if ( myDirCoef * myPrev->myDirCoef < 0 ) // different orientation of faces
+    {
+      myFaceNorm *= -1;
+      myDirCoef  *= -1;
+    }
+
+  }
+
+  //================================================================================
+  /*!
+   * \brief Compute myAngleWithPrev
+   */
+  //================================================================================
+
+  void BEdge::ComputeAngle( bool theReverseAngle )
+  {
+    myAngleWithPrev = ACos( myDir.Dot( myPrev->myDir.Reversed() ));
+
+    bool isObtuse;
+    gp_XYZ inNewFaceDir = myDir - myPrev->myDir;
+    double         dot1 = myDir.Dot( myPrev->myFaceNorm );
+    double         dot2 = myPrev->myDir.Dot( myFaceNorm );
+    bool     isOverlap1 = ( inNewFaceDir * myPrev->GetInFaceDir() > 0 );
+    bool     isOverlap2 = ( inNewFaceDir * GetInFaceDir()         > 0 );
+    if ( !myPrev->myFace )
+      isObtuse = isOverlap1;
+    else if  ( !myFace )
+      isObtuse = isOverlap2;
+    else
+    {
+      isObtuse = ( dot1 > 0 || dot2 < 0 ); // suppose face normals point outside the border
+      if ( theReverseAngle )
+        isObtuse = !isObtuse;
+    }
+    if ( isObtuse )
+    {
+      myAngleWithPrev = 2 * M_PI - myAngleWithPrev;
+    }
+
+    // if ( ! isObtuse )
+    //   isObtuse =
+    //     isSmallAngle( 1 - myDir.CrossSquareMagnitude( myPrev->myDir )); // edges co-directed
+
+    myOverlapAngle = 0;
+    //if ( !isObtuse )
+    {
+      // check if myFace and a triangle built on this and prev edges overlap
+      if ( isOverlap1 )
+      {
+        double cos2 = dot1 * dot1 / myPrev->myFaceNorm.SquareModulus();
+        myOverlapAngle += 0.5 * M_PI * ( 1 - cos2 );
+      }
+      if ( isOverlap2 )
+      {
+        double cos2 = dot2 * dot2 / myFaceNorm.SquareModulus();
+        myOverlapAngle += 0.5 * M_PI * ( 1 - cos2 );
+      }
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Check if myFace is overlapped by a triangle formed by myNode's and a
+   *        given node. If so, create shifted nodes to avoid overlapping
+   */
+  //================================================================================
+
+  void BEdge::ShiftOverlapped( const SMDS_MeshNode*                  theOppNode,
+                               const TFaceIndMap&                    theCapFaceWithBordInd,
+                               SMDS_Mesh&                            theMesh,
+                               std::vector<const SMDS_MeshElement*>& theNewFaces )
+  {
+    if ( myNode1Shift && myNode2Shift )
+      return;
+
+    gp_XYZ inNewFaceDir = SMESH_NodeXYZ( theOppNode ) - SMESH_NodeXYZ( myNode1 );
+    double          dot = inNewFaceDir.Dot( myFaceNorm );
+    double         cos2 = dot * dot / myFaceNorm.SquareModulus() / inNewFaceDir.SquareModulus();
+    bool      isOverlap = ( isSmallAngle( 1 - cos2 ) && GetInFaceDir() * inNewFaceDir > 0 );
+
+    if ( isOverlap )
+    {
+      gp_XYZ shift = myFaceNorm / myLength / 4;
+      if ( myFace )
+        shift.Reverse();
+      if ( !myNode1Shift )
+      {
+        gp_XYZ p = SMESH_NodeXYZ( myNode1 ) + shift;
+        myNode1Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() );
+        myPrev->myNode2Shift = myNode1Shift;
+      }
+      if ( !myNode2Shift )
+      {
+        gp_XYZ p = SMESH_NodeXYZ( myNode2 ) + shift;
+        myNode2Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() );
+        myNext->myNode1Shift = myNode2Shift;
+      }
+
+      // MakeShiftfFaces() for already created cap faces
+      for ( int is2nd = 0; is2nd < 2; ++is2nd )
+      {
+        const SMDS_MeshNode* ns = is2nd ? myNode2Shift : myNode1Shift;
+        const SMDS_MeshNode* n  = is2nd ? myNode2 : myNode1;
+        if ( !ns ) continue;
+
+        SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
+        while ( fIt->more() )
+        {
+          const SMDS_MeshElement* f = fIt->next();
+          if ( !f->isMarked() ) continue;
+
+          TFaceIndMap::const_iterator f2i = theCapFaceWithBordInd.find( f );
+          if ( f2i == theCapFaceWithBordInd.end() )
+            continue;
+          const SMDS_MeshNode* nf1 = f->GetNode(  f2i->second );
+          const SMDS_MeshNode* nf2 = f->GetNode(( f2i->second+1 ) % f->NbNodes() );
+          if ( nf1 == n || nf2 == n )
+          {
+            BEdge tmpE;
+            tmpE.myPrev = tmpE.myNext = this;
+            tmpE.Init( nf1, nf2, f, nf1, nf2 );
+            if ( !tmpE.myNode1Shift && !tmpE.myNode2Shift )
+              tmpE.Init( nf2, nf1, f, nf2, nf1 );
+            tmpE.myFace = f;
+            tmpE.MakeShiftfFaces( theMesh, theNewFaces, tmpE.myDirCoef < 0 );
+          }
+          std::vector< const SMDS_MeshNode* > nodes( f->begin_nodes(), f->end_nodes() );
+          nodes[ f->GetNodeIndex( n ) ] = ns;
+          theMesh.ChangeElementNodes( f, &nodes[0], nodes.size() );
+        }
+      }
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Create a triangle
+   */
+  //================================================================================
+
+  const SMDS_MeshElement* MakeTria( SMDS_Mesh&           mesh,
+                                    const SMDS_MeshNode* n1,
+                                    const SMDS_MeshNode* n2,
+                                    const SMDS_MeshNode* n3,
+                                    const bool           isReverse )
+  {
+    if ( isReverse )
+      return mesh.AddFace( n1, n3, n2 );
+    return mesh.AddFace( n1, n2, n3 );
+  }
+
+  //================================================================================
+  /*!
+   * \brief Create a quadrangle
+   */
+  //================================================================================
+
+  // const SMDS_MeshElement* MakeQuad( SMDS_Mesh&           mesh,
+  //                                   const SMDS_MeshNode* n1,
+  //                                   const SMDS_MeshNode* n2,
+  //                                   const SMDS_MeshNode* n3,
+  //                                   const SMDS_MeshNode* n4,
+  //                                   const bool           isReverse )
+  // {
+  //   if ( isReverse )
+  //     return mesh.AddFace( n4, n3, n2, n1 );
+  //   return mesh.AddFace( n1, n2, n3, n4 );
+  // }
+
+  //================================================================================
+  /*!
+   * \brief Create faces on myNode* and myNode*Shift
+   */
+  //================================================================================
+
+  void BEdge::MakeShiftfFaces(SMDS_Mesh&                            mesh,
+                              std::vector<const SMDS_MeshElement*>& newFaces,
+                              const bool                            isReverse )
+  {
+    if ( !myFace )
+      return;
+    if ( myNode1Shift && myNode2Shift )
+    {
+      newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse ));
+      newFaces.push_back( MakeTria( mesh, myNode1, myNode2Shift, myNode1Shift, isReverse ));
+    }
+    else if ( myNode1Shift )
+    {
+      newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode1Shift, isReverse ));
+    }
+    else if ( myNode2Shift )
+    {
+      newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse ));
+    }
+  }
+
+} // namespace
+
+//================================================================================
+/*!
+ * \brief Fill with 2D elements a hole defined by a TFreeBorder
+ */
+//================================================================================
+
+void SMESH_MeshAlgos::FillHole(const SMESH_MeshAlgos::TFreeBorder &  theFreeBorder,
+                               SMDS_Mesh&                            theMesh,
+                               std::vector<const SMDS_MeshElement*>& theNewFaces)
+{
+  if ( theFreeBorder.size() < 4 ||                // at least 3 nodes
+       theFreeBorder[0] != theFreeBorder.back() ) // the hole must be closed
+    return;
+
+  // prepare data of the border
+
+  ObjectPool< BEdge > edgeAllocator;
+  boost::intrusive::circular_list_algorithms< BEdge > circularList;
+  BEdge* edge;
+  BEdge* edge0 = edgeAllocator.getNew();
+  BEdge* edgePrev = edge0;
+  circularList.init_header( edge0 );
+  edge0->Init( theFreeBorder[0], theFreeBorder[1], 0 );
+  Bnd_B3d box;
+  box.Add( SMESH_NodeXYZ( edge0->myNode1 ));
+  for ( size_t i = 2; i < theFreeBorder.size(); ++i )
+  {
+    edge = edgeAllocator.getNew();
+    circularList.link_after( edgePrev, edge );
+    edge->Init( theFreeBorder[i-1], theFreeBorder[i] );
+    edge->ComputeAngle();
+    edgePrev = edge;
+    box.Add( SMESH_NodeXYZ( edge->myNode1 ));
+  }
+  edge0->ComputeAngle();
+
+  // check if face normals point outside the border
+
+  gp_XYZ hSize = 0.5 * ( box.CornerMax() - box.CornerMin() );
+  const double hDelta = 1e-6 * hSize.Modulus();
+  hSize -= gp_XYZ( hDelta, hDelta, hDelta );
+  if ( hSize.X() < 0 ) hSize.SetX(hDelta);
+  if ( hSize.Y() < 0 ) hSize.SetY(hDelta);
+  if ( hSize.Z() < 0 ) hSize.SetZ(hDelta);
+  box.SetHSize( hSize ); // decrease the box by hDelta
+
+  size_t nbEdges = theFreeBorder.size() - 1;
+  edge = edge0;
+  int nbRev = 0, nbFrw = 0;
+  double angTol = M_PI - ( nbEdges - 2 ) * M_PI / nbEdges, sumDirCoeff = 0;
+  for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
+  {
+    if ( box.IsOut( SMESH_NodeXYZ( edge->myNode1 )) &&
+         edge->myOverlapAngle < 0.1 * M_PI )
+    {
+      nbRev += edge->myAngleWithPrev > M_PI + angTol;
+      nbFrw += edge->myAngleWithPrev < M_PI - angTol;
+    }
+    sumDirCoeff += edge->myDirCoef;
+
+    // unmark all adjacent faces, new faces will be marked
+    SMDS_ElemIteratorPtr fIt = edge->myNode1->GetInverseElementIterator( SMDSAbs_Face );
+    while ( fIt->more() )
+      fIt->next()->setIsMarked( false );
+  }
+  bool isReverseAngle = ( nbRev > nbFrw ); // true == face normals point inside the border
+  //std::cout << "nbRev="<< nbRev << ", nbFrw="<< nbFrw<<std::endl;
+
+  // sort border edges by myAngleWithPrev
+
+  TAngleMap edgesByAngle;
+  bool useOverlap = true; // to add BEdge.myOverlapAngle when filling edgesByAngle
+  edge = edge0;
+  for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
+    edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap );
+
+  // create triangles to fill the hole
+
+  //compare order of nodes in the edges with their order in faces
+  bool isReverse = sumDirCoeff > 0.5 * nbEdges;
+
+  // faces filling the hole (cap faces) and indices of border edges in them
+  TFaceIndMap capFaceWithBordInd;
+
+  theNewFaces.reserve( nbEdges - 2 );
+  std::vector< const SMDS_MeshNode* > nodes(3);
+  while ( edgesByAngle.size() > 2 )
+  {
+    TAngleMap::iterator a2e = edgesByAngle.begin();
+    if ( useOverlap && a2e->first > M_PI - angTol ) // all new triangles need shift
+    {
+      // re-sort the edges
+      useOverlap = false;
+      edge    = a2e->second;
+      nbEdges = edgesByAngle.size();
+      edgesByAngle.clear();
+      for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
+        edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap );
+      a2e = edgesByAngle.begin();
+    }
+    edge     = a2e->second;
+    edgePrev = edge->myPrev;
+
+    // create shift nodes and faces
+    edgePrev->ShiftOverlapped( edge->myNode2, capFaceWithBordInd, theMesh, theNewFaces );
+    edge->ShiftOverlapped( edgePrev->myNode1, capFaceWithBordInd, theMesh, theNewFaces );
+    edge    ->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
+    edgePrev->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
+
+    // make a cap face
+    //nodes.resize( 3 );
+    nodes[0] = edgePrev->myNode1Shift ? edgePrev->myNode1Shift : edgePrev->myNode1;
+    nodes[1] = edgePrev->myNode2Shift ? edgePrev->myNode2Shift : edgePrev->myNode2;
+    nodes[2] = edge->myNode2Shift     ? edge->myNode2Shift     : edge->myNode2;
+    theNewFaces.push_back( MakeTria( theMesh, nodes[0], nodes[1], nodes[2], isReverse ));
+    // std::cout << nodes[1]->GetID() << " "  << nodes[0]->GetID() << " "  << nodes[2]->GetID()
+    //           << " " << edge->myAngleWithPrev << std::endl;
+
+    // remember a border edge within the new cap face
+    theNewFaces.back()->setIsMarked( true );
+    if ( edgePrev->myFace )
+      capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), isReverse ? 2 : 0 ));
+    if ( edge->myFace )
+      capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), 1 ));
+
+    // remove edgePrev from the list and update <edge>
+    edgesByAngle.erase( edgePrev->myAngleMapPos );
+    circularList.unlink( edgePrev ); // remove edgePrev from the border
+
+    edge->Init( edgePrev->myNode1, edge->myNode2, theNewFaces.back(), nodes[0], nodes[2] );
+    edge->ComputeAngle( isReverseAngle );
+    edge->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap );
+    edge->myNext->ComputeAngle( isReverseAngle );
+    edge->myNext->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap );
+    // std::cout << "A " << edge->myNode1->GetID() << " " << edge->myAngleWithPrev
+    //           << " " << edge->myNext->myNode1->GetID() << " " << edge->myNext->myAngleWithPrev
+    //           << std::endl;
+  }
+  edge = edgesByAngle.begin()->second;
+  edge->        MakeShiftfFaces( theMesh, theNewFaces, isReverse );
+  edge->myNext->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
+}

src/SMESHUtils/SMESH_FreeBorders.cxx

@@ -131,8 +131,11 @@ namespace
       if ( myID < 0 )
       {
         myID = id;
-        if ( myNext )
+
-          myNext->SetID( id + 1 );
+        for ( BEdge* be = myNext; be && be->myID < 0; be = be->myNext )
+        {
+          be->myID = ++id;
+        }
       }
     }
     //================================================================================
@@ -821,3 +824,136 @@ void SMESH_MeshAlgos::FindCoincidentFreeBorders(SMDS_Mesh&              mesh,
 
 } // SMESH_MeshAlgos::FindCoincidentFreeBorders()
 
+//================================================================================
+/*
+ * Returns all TFreeBorder's. Optionally check if the mesh is manifold
+ * and if faces are correctly oriented.
+ */
+//================================================================================
+
+void SMESH_MeshAlgos::FindFreeBorders(SMDS_Mesh&       theMesh,
+                                      TFreeBorderVec & theFoundFreeBordes,
+                                      const bool       theClosedOnly,
+                                      bool*            theIsManifold,
+                                      bool*            theIsGoodOri)
+{
+  bool isManifold = true;
+
+  // find free links
+  typedef NCollection_DataMap<SMESH_TLink, const SMDS_MeshElement*, SMESH_TLink > TLink2FaceMap;
+  TLink2FaceMap linkMap;
+  int nbSharedLinks = 0;
+  SMDS_FaceIteratorPtr faceIt = theMesh.facesIterator();
+  while ( faceIt->more() )
+  {
+    const SMDS_MeshElement* face = faceIt->next();
+    if ( !face ) continue;
+
+    const SMDS_MeshNode*     n0 = face->GetNode( face->NbNodes() - 1 );
+    SMDS_NodeIteratorPtr nodeIt = face->interlacedNodesIterator();
+    while ( nodeIt->more() )
+    {
+      const SMDS_MeshNode* n1 = nodeIt->next();
+      SMESH_TLink link( n0, n1 );
+      if ( const SMDS_MeshElement** faceInMap = linkMap.ChangeSeek( link ))
+      {
+        if ( *faceInMap )
+        {
+          if ( theIsGoodOri && *theIsGoodOri && !IsRightOrder( *faceInMap, n1, n0 ))
+            *theIsGoodOri = false;
+        }
+        else
+        {
+          isManifold = false;
+        }
+        nbSharedLinks += bool( *faceInMap );
+        *faceInMap = 0;
+      }
+      else
+      {
+        linkMap.Bind( link, face );
+      }
+      n0 = n1;
+    }
+  }
+  if ( theIsManifold )
+    *theIsManifold = isManifold;
+
+  if ( linkMap.Extent() == nbSharedLinks )
+    return;
+
+  // form free borders
+  std::set   < BNode > bNodes;
+  std::vector< BEdge > bEdges( linkMap.Extent() - nbSharedLinks );
+
+  TLink2FaceMap::Iterator linkIt( linkMap );
+  for ( int iEdge = 0; linkIt.More(); linkIt.Next() )
+  {
+    if ( !linkIt.Value() ) continue;
+    const SMESH_TLink & link = linkIt.Key();
+    std::set< BNode >::iterator n1 = bNodes.insert( BNode( link.node1() )).first;
+    std::set< BNode >::iterator n2 = bNodes.insert( BNode( link.node2() )).first;
+    bEdges[ iEdge ].Set( &*n1, &*n2, linkIt.Value(), iEdge+1 );
+    n1->AddLinked( & bEdges[ iEdge ] );
+    n2->AddLinked( & bEdges[ iEdge ] );
+    ++iEdge;
+  }
+  linkMap.Clear();
+
+  // assign IDs to borders
+  std::vector< BEdge* > borders; // 1st of connected (via myPrev and myNext) edges
+  std::set< BNode >::iterator bn = bNodes.begin();
+  for ( ; bn != bNodes.end(); ++bn )
+  {
+    for ( size_t i = 0; i < bn->myLinkedEdges.size(); ++i )
+    {
+      if ( bn->myLinkedEdges[i]->myBorderID < 0 )
+      {
+        BEdge* be = bn->myLinkedEdges[i];
+        int borderID = borders.size();
+        borders.push_back( be );
+        for ( ; be && be->myBorderID < 0; be = be->myNext )
+        {
+          be->myBorderID = borderID;
+          be->Orient();
+        }
+        bool isClosed = ( be == bn->myLinkedEdges[i] );
+        if ( !isClosed && theClosedOnly )
+        {
+          borders.pop_back();
+          continue;
+        }
+        be = bn->myLinkedEdges[i]->myPrev;
+        for ( ; be && be->myBorderID < 0; be = be->myPrev )
+        {
+          be->myBorderID = borderID;
+          be->Orient();
+        }
+        if ( !isClosed )
+          while ( borders.back()->myPrev )
+            borders.back() = borders.back()->myPrev;
+      }
+    }
+  }
+  theFoundFreeBordes.resize( borders.size() );
+  for ( size_t i = 0; i < borders.size(); ++i )
+  {
+    TFreeBorder & bordNodes = theFoundFreeBordes[ i ];
+    BEdge*               be = borders[i];
+
+    size_t cnt = 1;
+    for ( be = be->myNext; be && be != borders[i]; be = be->myNext )
+      ++cnt;
+    bordNodes.resize( cnt + 1 );
+
+    BEdge* beLast;
+    for ( be = borders[i], cnt = 0;
+          be && cnt < bordNodes.size()-1;
+          be = be->myNext, ++cnt )
+    {
+      bordNodes[ cnt ] = be->myBNode1->Node();
+      beLast = be;
+    }
+    bordNodes.back() = beLast->myBNode2->Node();
+  }
+}

src/SMESHUtils/SMESH_MeshAlgos.cxx

@@ -35,6 +35,8 @@
 #include "SMDS_VolumeTool.hxx"
 #include "SMESH_OctreeNode.hxx"
 
+#include <Utils_SALOME_Exception.hxx>
+
 #include <GC_MakeSegment.hxx>
 #include <GeomAPI_ExtremaCurveCurve.hxx>
 #include <Geom_Line.hxx>
@@ -228,12 +230,12 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
                       SMDSAbs_ElementType  elemType,
                       SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(),
                       double               tolerance = NodeRadius );
-    void prepare(); // !!!call it before calling the following methods!!!
     void getElementsNearPoint( const gp_Pnt& point, vector<const SMDS_MeshElement*>& foundElems );
     void getElementsNearLine ( const gp_Ax1& line,  vector<const SMDS_MeshElement*>& foundElems);
     void getElementsInSphere ( const gp_XYZ&                    center,
                                const double                     radius,
                                vector<const SMDS_MeshElement*>& foundElems);
+    ElementBndBoxTree* getLeafAtPoint( const gp_XYZ& point );
 
   protected:
     ElementBndBoxTree() {}
@@ -337,19 +339,6 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     }
   }
 
-  //================================================================================
-  /*!
-   * \brief Un-mark all elements
-   */
-  //================================================================================
-
-  void ElementBndBoxTree::prepare()
-  {
-    // TElementBoxPool& elBoPool = getElementBoxPool();
-    // for ( size_t i = 0; i < elBoPool.nbElements(); ++i )
-    //   const_cast< ElementBox* >( elBoPool[ i ])->_isMarked = false;
-  }
-
   //================================================================================
   /*!
    * \brief Return elements which can include the point
@@ -471,6 +460,30 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     }
   }
 
+  //================================================================================
+  /*!
+   * \brief Return a leaf including a point
+   */
+  //================================================================================
+
+  ElementBndBoxTree* ElementBndBoxTree::getLeafAtPoint( const gp_XYZ& point )
+  {
+    if ( getBox()->IsOut( point ))
+      return 0;
+
+    if ( isLeaf() )
+    {
+      return this;
+    }
+    else
+    {
+      for (int i = 0; i < 8; i++)
+        if ( ElementBndBoxTree* l = ((ElementBndBoxTree*) myChildren[i])->getLeafAtPoint( point ))
+          return l;
+    }
+    return 0;
+  }
+
   //================================================================================
   /*!
    * \brief Construct the element box
@@ -539,13 +552,16 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   virtual const SMDS_MeshElement* FindClosestTo( const gp_Pnt&       point,
                                                  SMDSAbs_ElementType type );
 
-  void GetElementsNearLine( const gp_Ax1&                      line,
+  virtual void GetElementsNearLine( const gp_Ax1&                      line,
-                            SMDSAbs_ElementType                type,
+                                    SMDSAbs_ElementType                type,
-                            vector< const SMDS_MeshElement* >& foundElems);
+                                    vector< const SMDS_MeshElement* >& foundElems);
-  void GetElementsInSphere( const gp_XYZ&                      center,
+  virtual void GetElementsInSphere( const gp_XYZ&                      center,
-                            const double                       radius,
+                                    const double                       radius,
-                            SMDSAbs_ElementType                type,
+                                    SMDSAbs_ElementType                type,
-                            vector< const SMDS_MeshElement* >& foundElems);
+                                    vector< const SMDS_MeshElement* >& foundElems);
+  virtual gp_XYZ Project(const gp_Pnt&            point,
+                         SMDSAbs_ElementType      type,
+                         const SMDS_MeshElement** closestElem);
   double getTolerance();
   bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
                             const double tolerance, double & param);
@@ -834,10 +850,6 @@ FindElementsByPoint(const gp_Pnt&                      point,
     {
       _ebbTree[_elementType] = new ElementBndBoxTree( *_mesh, type, _meshPartIt, tolerance );
     }
-    else
-    {
-      _ebbTree[ type ]->prepare();
-    }
     vector< const SMDS_MeshElement* > suspectElems;
     _ebbTree[ type ]->getElementsNearPoint( point, suspectElems );
     vector< const SMDS_MeshElement* >::iterator elem = suspectElems.begin();
@@ -863,13 +875,13 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
   const SMDS_MeshElement* closestElem = 0;
   _elementType = type;
 
-  if ( type == SMDSAbs_Face || type == SMDSAbs_Volume )
+  if ( type == SMDSAbs_Face ||
+       type == SMDSAbs_Volume ||
+       type == SMDSAbs_Edge )
   {
     ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
     if ( !ebbTree )
       ebbTree = new ElementBndBoxTree( *_mesh, type, _meshPartIt );
-    else
-      ebbTree->prepare();
 
     vector<const SMDS_MeshElement*> suspectElems;
     ebbTree->getElementsNearPoint( point, suspectElems );
@@ -885,7 +897,6 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
         radius = ebbTree->maxSize() / pow( 2., getTreeHeight()) / 2;
       while ( suspectElems.empty() )
       {
-        ebbTree->prepare();
         ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems );
         radius *= 1.1;
       }
@@ -956,8 +967,6 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
   ElementBndBoxTree*& ebbTree = _ebbTree[ SMDSAbs_Face ];
   if ( !ebbTree )
     ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
-  else
-    ebbTree->prepare();
 
   // Algo: analyse transition of a line starting at the point through mesh boundary;
   // try three lines parallel to axis of the coordinate system and perform rough
@@ -974,7 +983,6 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
     gp_Lin line    ( lineAxis );
 
     vector<const SMDS_MeshElement*> suspectFaces; // faces possibly intersecting the line
-    if ( axis > 0 ) ebbTree->prepare();
     ebbTree->getElementsNearLine( lineAxis, suspectFaces );
 
     // Intersect faces with the line
@@ -1187,8 +1195,6 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1&
   ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
   if ( !ebbTree )
     ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
-  else
-    ebbTree->prepare();
 
   ebbTree->getElementsNearLine( line, foundElems );
 }
@@ -1208,12 +1214,59 @@ void SMESH_ElementSearcherImpl::GetElementsInSphere( const gp_XYZ&
   ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
   if ( !ebbTree )
     ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
-  else
-    ebbTree->prepare();
 
   ebbTree->getElementsInSphere( center, radius, foundElems );
 }
 
+//=======================================================================
+/*
+ * \brief Return a projection of a given point to a mesh.
+ *        Optionally return the closest element
+ */
+//=======================================================================
+
+gp_XYZ SMESH_ElementSearcherImpl::Project(const gp_Pnt&            point,
+                                          SMDSAbs_ElementType      type,
+                                          const SMDS_MeshElement** closestElem)
+{
+  _elementType = type;
+  if ( _mesh->GetMeshInfo().NbElements( _elementType ) == 0 )
+    throw SALOME_Exception( LOCALIZED( "No elements of given type in the mesh" ));
+
+  ElementBndBoxTree*& ebbTree = _ebbTree[ _elementType ];
+  if ( !ebbTree )
+    ebbTree = new ElementBndBoxTree( *_mesh, _elementType );
+
+  gp_XYZ p = point.XYZ();
+  ElementBndBoxTree* ebbLeaf = ebbTree->getLeafAtPoint( p );
+  const Bnd_B3d* box = ebbLeaf->getBox();
+  double radius = ( box->CornerMax() - box->CornerMin() ).Modulus();
+
+  vector< const SMDS_MeshElement* > elems;
+  ebbTree->getElementsInSphere( p, radius, elems );
+  while ( elems.empty() )
+  {
+    radius *= 1.5;
+    ebbTree->getElementsInSphere( p, radius, elems );
+  }
+  gp_XYZ proj, bestProj;
+  const SMDS_MeshElement* elem = 0;
+  double minDist = 2 * radius;
+  for ( size_t i = 0; i < elems.size(); ++i )
+  {
+    double d = SMESH_MeshAlgos::GetDistance( elems[i], p, &proj );
+    if ( d < minDist )
+    {
+      bestProj = proj;
+      elem = elems[i];
+      minDist = d;
+    }
+  }
+  if ( closestElem ) *closestElem = elem;
+
+  return bestProj;
+}
+
 //=======================================================================
 /*!
  * \brief Return true if the point is IN or ON of the element
@@ -1461,17 +1514,19 @@ namespace
 //=======================================================================
 
 double SMESH_MeshAlgos::GetDistance( const SMDS_MeshElement* elem,
-                                     const gp_Pnt&           point )
+                                     const gp_Pnt&           point,
+                                     gp_XYZ*                 closestPnt )
 {
   switch ( elem->GetType() )
   {
   case SMDSAbs_Volume:
-    return GetDistance( dynamic_cast<const SMDS_MeshVolume*>( elem ), point);
+    return GetDistance( dynamic_cast<const SMDS_MeshVolume*>( elem ), point, closestPnt );
   case SMDSAbs_Face:
-    return GetDistance( dynamic_cast<const SMDS_MeshFace*>( elem ), point);
+    return GetDistance( dynamic_cast<const SMDS_MeshFace*>( elem ), point, closestPnt );
   case SMDSAbs_Edge:
-    return GetDistance( dynamic_cast<const SMDS_MeshEdge*>( elem ), point);
+    return GetDistance( dynamic_cast<const SMDS_MeshEdge*>( elem ), point, closestPnt );
   case SMDSAbs_Node:
+    if ( closestPnt ) *closestPnt = SMESH_TNodeXYZ( elem );
     return point.Distance( SMESH_TNodeXYZ( elem ));
   default:;
   }
@@ -1487,9 +1542,10 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshElement* elem,
 //=======================================================================
 
 double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face,
-                                     const gp_Pnt&        point )
+                                     const gp_Pnt&        point,
+                                     gp_XYZ*              closestPnt )
 {
-  double badDistance = -1;
+  const double badDistance = -1;
   if ( !face ) return badDistance;
 
   // coordinates of nodes (medium nodes, if any, ignored)
@@ -1533,7 +1589,7 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face,
   trsf.Transforms( tmpPnt );
   gp_XY point2D( tmpPnt.X(), tmpPnt.Z() );
 
-  // loop on segments of the face to analyze point position ralative to the face
+  // loop on edges of the face to analyze point position ralative to the face
   set< PointPos > pntPosSet;
   for ( size_t i = 1; i < xy.size(); ++i )
   {
@@ -1543,31 +1599,40 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face,
 
   // compute distance
   PointPos pos = *pntPosSet.begin();
-  // cout << "Face " << face->GetID() << " DIST: ";
   switch ( pos._name )
   {
-  case POS_LEFT: {
+  case POS_LEFT:
-    // point is most close to a segment
+  {
-    gp_Vec p0p1( point, xyz[ pos._index ] );
+    // point is most close to an edge
-    gp_Vec p1p2( xyz[ pos._index ], xyz[ pos._index+1 ]); // segment vector
+    gp_Vec edge( xyz[ pos._index ], xyz[ pos._index+1 ]);
-    p1p2.Normalize();
+    gp_Vec n1p ( xyz[ pos._index ], point  );
-    double projDist = p0p1 * p1p2; // distance projected to the segment
+    double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
-    gp_Vec projVec = p1p2 * projDist;
+    // projection of the point on the edge
-    gp_Vec distVec = p0p1 - projVec;
+    gp_XYZ proj = ( 1. - u ) * xyz[ pos._index ] + u * xyz[ pos._index+1 ];
-    // cout << distVec.Magnitude()  << ", SEG " << face->GetNode(pos._index)->GetID()
+    if ( closestPnt ) *closestPnt = proj;
-    //      << " - " << face->GetNodeWrap(pos._index+1)->GetID() << endl;
+    return point.Distance( proj );
-    return distVec.Magnitude();
   }
-  case POS_RIGHT: {
+  case POS_RIGHT:
+  {
     // point is inside the face
-    double distToFacePlane = tmpPnt.Y();
+    double distToFacePlane = Abs( tmpPnt.Y() );
-    // cout << distToFacePlane << ", INSIDE " << endl;
+    if ( closestPnt )
-    return Abs( distToFacePlane );
+    {
+      if ( distToFacePlane < std::numeric_limits<double>::min() ) {
+        *closestPnt = point.XYZ();
+      }
+      else {
+        tmpPnt.SetY( 0 );
+        trsf.Inverted().Transforms( tmpPnt );
+        *closestPnt = tmpPnt;
+      }
+    }
+    return distToFacePlane;
   }
-  case POS_VERTEX: {
+  case POS_VERTEX:
+  {
     // point is most close to a node
     gp_Vec distVec( point, xyz[ pos._index ]);
-    // cout << distVec.Magnitude()  << " VERTEX " << face->GetNode(pos._index)->GetID() << endl;
     return distVec.Magnitude();
   }
   default:;
@@ -1581,7 +1646,9 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face,
  */
 //=======================================================================
 
-double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg, const gp_Pnt& point )
+double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg,
+                                     const gp_Pnt&        point,
+                                     gp_XYZ*              closestPnt )
 {
   double dist = Precision::Infinite();
   if ( !seg ) return dist;
@@ -1597,16 +1664,25 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg, const gp_Pnt& poi
   {
     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
+    double d, u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
     if ( u <= 0. ) {
-      dist = Min( dist, n1p.SquareMagnitude() );
+      if (( d = n1p.SquareMagnitude() ) < dist ) {
+        dist = d;
+        if ( closestPnt ) *closestPnt = xyz[i-1];
+      }
     }
     else if ( u >= 1. ) {
-      dist = Min( dist, point.SquareDistance( xyz[i] ));
+      if (( d = point.SquareDistance( xyz[i] )) < dist ) {
+        dist = d;
+        if ( closestPnt ) *closestPnt = xyz[i];
+      }
     }
     else {
-      gp_XYZ proj = ( 1. - u ) * xyz[i-1] + u * xyz[i]; // projection of the point on the edge
+      gp_XYZ proj = xyz[i-1] + u * edge.XYZ(); // projection of the point on the edge
-      dist = Min( dist, point.SquareDistance( proj ));
+      if (( d = point.SquareDistance( proj )) < dist ) {
+        dist = d;
+        if ( closestPnt ) *closestPnt = proj;
+      }
     }
   }
   return Sqrt( dist );
@@ -1620,7 +1696,9 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg, const gp_Pnt& poi
  */
 //=======================================================================
 
-double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt& point )
+double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume,
+                                     const gp_Pnt&          point,
+                                     gp_XYZ*                closestPnt )
 {
   SMDS_VolumeTool vTool( volume );
   vTool.SetExternalNormal();
@@ -1628,6 +1706,8 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt
 
   double n[3], bc[3];
   double minDist = 1e100, dist;
+  gp_XYZ closeP = point.XYZ();
+  bool isOut = false;
   for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
   {
     // skip a facet with normal not "looking at" the point
@@ -1644,23 +1724,34 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt
     case 3:
     {
       SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ] );
-      dist = GetDistance( &tmpFace, point );
+      dist = GetDistance( &tmpFace, point, closestPnt );
       break;
     }
     case 4:
     {
       SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ], nodes[ 3*iQ ]);
-      dist = GetDistance( &tmpFace, point );
+      dist = GetDistance( &tmpFace, point, closestPnt );
       break;
     }
     default:
       vector<const SMDS_MeshNode *> nvec( nodes, nodes + vTool.NbFaceNodes( iF ));
       SMDS_PolygonalFaceOfNodes tmpFace( nvec );
-      dist = GetDistance( &tmpFace, point );
+      dist = GetDistance( &tmpFace, point, closestPnt );
+    }
+    if ( dist < minDist )
+    {
+      minDist = dist;
+      isOut = true;
+      if ( closestPnt ) closeP = *closestPnt;
     }
-    minDist = Min( minDist, dist );
   }
-  return minDist;
+  if ( isOut )
+  {
+    if ( closestPnt ) *closestPnt = closeP;
+    return minDist;
+  }
+
+  return 0; // point is inside the volume
 }
 
 //================================================================================
@@ -1804,6 +1895,34 @@ vector< const SMDS_MeshNode*> SMESH_MeshAlgos::GetCommonNodes(const SMDS_MeshEle
       common.push_back( e1->GetNode( i ));
   return common;
 }
+//================================================================================
+/*!
+ * \brief Return true if node1 encounters first in the face and node2, after
+ */
+//================================================================================
+
+bool SMESH_MeshAlgos::IsRightOrder( const SMDS_MeshElement* face,
+                                    const SMDS_MeshNode*    node0,
+                                    const SMDS_MeshNode*    node1 )
+{
+  int i0 = face->GetNodeIndex( node0 );
+  int i1 = face->GetNodeIndex( node1 );
+  if ( face->IsQuadratic() )
+  {
+    if ( face->IsMediumNode( node0 ))
+    {
+      i0 -= ( face->NbNodes()/2 - 1 );
+      i1 *= 2;
+    }
+    else
+    {
+      i1 -= ( face->NbNodes()/2 - 1 );
+      i0 *= 2;
+    }
+  }
+  int diff = i1 - i0;
+  return ( diff == 1 ) || ( diff == -face->NbNodes()+1 );
+}
 
 //=======================================================================
 /*!

src/SMESHUtils/SMESH_MeshAlgos.hxx

@@ -100,6 +100,15 @@ struct SMESHUtils_EXPORT SMESH_ElementSearcher
    * \brief Find out if the given point is out of closed 2D mesh.
    */
   virtual TopAbs_State GetPointState(const gp_Pnt& point) = 0;
+
+  /*!
+   * \brief Return a projection of a given point to a 2D mesh.
+   *        Optionally return the closest face
+   */
+  virtual gp_XYZ Project(const gp_Pnt&            point,
+                         SMDSAbs_ElementType      type,
+                         const SMDS_MeshElement** closestFace= 0) = 0;
+
   virtual ~SMESH_ElementSearcher();
 };
 
@@ -112,16 +121,16 @@ namespace SMESH_MeshAlgos
   bool IsOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol );
 
   SMESHUtils_EXPORT
-  double GetDistance( const SMDS_MeshElement* elem, const gp_Pnt& point );
+  double GetDistance( const SMDS_MeshElement* elem, const gp_Pnt& point, gp_XYZ* closestPnt = 0 );
 
   SMESHUtils_EXPORT
-  double GetDistance( const SMDS_MeshEdge* edge, const gp_Pnt& point );
+  double GetDistance( const SMDS_MeshEdge* edge, const gp_Pnt& point, gp_XYZ* closestPnt = 0 );
 
   SMESHUtils_EXPORT
-  double GetDistance( const SMDS_MeshFace* face, const gp_Pnt& point );
+  double GetDistance( const SMDS_MeshFace* face, const gp_Pnt& point, gp_XYZ* closestPnt = 0 );
 
   SMESHUtils_EXPORT
-  double GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt& point );
+  double GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt& point, gp_XYZ* closestPnt = 0 );
 
   SMESHUtils_EXPORT
   void GetBarycentricCoords( const gp_XY& point,
@@ -153,6 +162,14 @@ namespace SMESH_MeshAlgos
   SMESHUtils_EXPORT
   std::vector< const SMDS_MeshNode*> GetCommonNodes(const SMDS_MeshElement* e1,
                                                     const SMDS_MeshElement* e2);
+  /*!
+   * \brief Return true if node1 encounters first in the face and node2, after.
+   *        The nodes are supposed to be neighbor nodes in the face.
+   */
+  SMESHUtils_EXPORT
+  bool IsRightOrder( const SMDS_MeshElement* face,
+                     const SMDS_MeshNode*    node0,
+                     const SMDS_MeshNode*    node1 );
 
   /*!
    * \brief Return SMESH_NodeSearcher. The caller is responsible for deleteing it
@@ -204,6 +221,27 @@ namespace SMESH_MeshAlgos
   void FindCoincidentFreeBorders(SMDS_Mesh&              mesh,
                                  double                  tolerance,
                                  CoincidentFreeBorders & foundFreeBordes);
+  /*!
+   * Returns all or only closed TFreeBorder's.
+   * Optionally check if the mesh is manifold and if faces are correctly oriented.
+   *
+   * (Implemented in ./SMESH_FreeBorders.cxx)
+   */
+  SMESHUtils_EXPORT
+  void FindFreeBorders(SMDS_Mesh&       mesh,
+                       TFreeBorderVec & foundFreeBordes,
+                       const bool       closedOnly,
+                       bool*            isManifold = 0,
+                       bool*            isGoodOri = 0);
+  /*!
+   * Fill a hole defined by a TFreeBorder with 2D elements.
+   *
+   * (Implemented in ./SMESH_FillHole.cxx)
+   */
+  SMESHUtils_EXPORT
+  void FillHole(const TFreeBorder &                   freeBorder,
+                SMDS_Mesh&                            mesh,
+                std::vector<const SMDS_MeshElement*>& newFaces);
 
 
   /*!

src/SMESH_I/SMESH_2smeshpy.cxx

@@ -2430,7 +2430,7 @@ void _pyMeshEditor::Process( const Handle(_pyCommand)& theCommand)
       "ExtrusionAlongPathX","ExtrusionAlongPathObject1D","ExtrusionAlongPathObject2D",
       "ExtrusionSweepObjects","RotationSweepObjects","ExtrusionAlongPathObjects",
       "Mirror","MirrorObject","Translate","TranslateObject","Rotate","RotateObject",
-      "FindCoincidentNodes","MergeNodes","FindEqualElements",
+      "FindCoincidentNodes","MergeNodes","FindEqualElements","FillHole",
       "MergeElements","MergeEqualElements","SewFreeBorders","SewConformFreeBorders",
       "FindCoincidentFreeBorders", "SewCoincidentFreeBorders",
       "SewBorderToSide","SewSideElements","ChangeElemNodes","GetLastCreatedNodes",

src/SMESH_I/SMESH_MeshEditor_i.cxx

@@ -4006,7 +4006,7 @@ SMESH_MeshEditor_i::ScaleMakeMesh(SMESH::SMESH_IDSource_ptr  theObject,
     // and then "GetGroups" using SMESH_Mesh::GetGroups()
 
     TPythonDump pydump; // to prevent dump at mesh creation
-    mesh = makeMesh( theMeshName );
+    mesh   = makeMesh( theMeshName );
     mesh_i = SMESH::DownCast<SMESH_Mesh_i*>( mesh );
 
     if ( mesh_i )
@@ -4593,6 +4593,151 @@ CORBA::Short SMESH_MeshEditor_i::GetPointState(CORBA::Double x,
   return 0;
 }
 
+//=======================================================================
+//function : IsManifold
+//purpose  : Check if a 2D mesh is manifold
+//=======================================================================
+
+CORBA::Boolean SMESH_MeshEditor_i::IsManifold()
+  throw (SALOME::SALOME_Exception)
+{
+  bool isManifold = true;
+
+  SMESH_TRY;
+  SMESH_MeshAlgos::TFreeBorderVec foundFreeBordes;
+  SMESH_MeshAlgos::FindFreeBorders( *getMeshDS(),
+                                    foundFreeBordes,
+                                    /*closedOnly=*/true,
+                                    &isManifold );
+  SMESH_CATCH( SMESH::throwCorbaException );
+
+  return isManifold;
+}
+
+//=======================================================================
+//function : IsCoherentOrientation2D
+//purpose  : Check if orientation of 2D elements is coherent
+//=======================================================================
+
+CORBA::Boolean SMESH_MeshEditor_i::IsCoherentOrientation2D()
+  throw (SALOME::SALOME_Exception)
+{
+  bool isGoodOri = true;
+
+  SMESH_TRY;
+  SMESH_MeshAlgos::TFreeBorderVec foundFreeBordes;
+  SMESH_MeshAlgos::FindFreeBorders( *getMeshDS(),
+                                    foundFreeBordes,
+                                    /*closedOnly=*/true,
+                                    /*isManifold=*/0,
+                                    &isGoodOri);
+  SMESH_CATCH( SMESH::throwCorbaException );
+
+  return isGoodOri;
+}
+
+//=======================================================================
+//function : FindFreeBorders
+//purpose  : Returns all or only closed FreeBorder's.
+//=======================================================================
+
+SMESH::ListOfFreeBorders* SMESH_MeshEditor_i::FindFreeBorders(CORBA::Boolean closedOnly)
+  throw (SALOME::SALOME_Exception)
+{
+  SMESH::ListOfFreeBorders_var resBorders = new SMESH::ListOfFreeBorders;
+  SMESH_TRY;
+
+  SMESH_MeshAlgos::TFreeBorderVec foundFreeBordes;
+  SMESH_MeshAlgos::FindFreeBorders( *getMeshDS(), foundFreeBordes, closedOnly );
+
+  resBorders->length( foundFreeBordes.size() );
+  for ( size_t i = 0; i < foundFreeBordes.size(); ++i )
+  {
+    const SMESH_MeshAlgos::TFreeBorder& bordNodes = foundFreeBordes[i];
+    SMESH::FreeBorder&                    bordOut = resBorders[i];
+    bordOut.nodeIDs.length( bordNodes.size() );
+    for ( size_t iN = 0; iN < bordNodes.size(); ++iN )
+      bordOut.nodeIDs[ iN ] = bordNodes[ iN ]->GetID();
+  }
+
+  SMESH_CATCH( SMESH::throwCorbaException );
+
+  return resBorders._retn();
+}
+
+//=======================================================================
+//function : FillHole
+//purpose  : Fill with 2D elements a hole defined by a FreeBorder.
+//=======================================================================
+
+void SMESH_MeshEditor_i::FillHole(const SMESH::FreeBorder& theHole)
+  throw (SALOME::SALOME_Exception)
+{
+  initData();
+
+  if ( theHole.nodeIDs.length() < 4 )
+    THROW_SALOME_CORBA_EXCEPTION("A hole should be bound by at least 3 nodes", SALOME::BAD_PARAM);
+  if ( theHole.nodeIDs[0] != theHole.nodeIDs[ theHole.nodeIDs.length()-1 ] )
+    THROW_SALOME_CORBA_EXCEPTION("Not closed hole boundary. "
+                                 "First and last nodes must be same", SALOME::BAD_PARAM);
+
+  SMESH_MeshAlgos::TFreeBorder bordNodes;
+  bordNodes.resize( theHole.nodeIDs.length() );
+  for ( size_t iN = 0; iN < theHole.nodeIDs.length(); ++iN )
+  {
+    bordNodes[ iN ] = getMeshDS()->FindNode( theHole.nodeIDs[ iN ]);
+    if ( !bordNodes[ iN ] )
+      THROW_SALOME_CORBA_EXCEPTION(SMESH_Comment("Node #") << theHole.nodeIDs[ iN ]
+                                   << " does not exist", SALOME::BAD_PARAM);
+  }
+
+  SMESH_TRY;
+
+  MeshEditor_I::TPreviewMesh* previewMesh = 0;
+  SMDS_Mesh* meshDS = getMeshDS();
+  if ( myIsPreviewMode )
+  {
+    // copy faces sharing nodes of theHole
+    TIDSortedElemSet holeFaces;
+    previewMesh = getPreviewMesh( SMDSAbs_Face );
+    for ( size_t i = 0; i < bordNodes.size(); ++i )
+    {
+      SMDS_ElemIteratorPtr fIt = bordNodes[i]->GetInverseElementIterator( SMDSAbs_Face );
+      while ( fIt->more() )
+      {
+        const SMDS_MeshElement* face = fIt->next();
+        if ( holeFaces.insert( face ).second )
+          previewMesh->Copy( face );
+      }
+      bordNodes[i] = previewMesh->GetMeshDS()->FindNode( bordNodes[i]->GetID() );
+      ASSERT( bordNodes[i] );
+    }
+    meshDS = previewMesh->GetMeshDS();
+  }
+
+  std::vector<const SMDS_MeshElement*> newFaces;
+  SMESH_MeshAlgos::FillHole( bordNodes, *meshDS, newFaces );
+
+  if ( myIsPreviewMode )
+  {
+    previewMesh->Clear();
+    for ( size_t i = 0; i < newFaces.size(); ++i )
+      previewMesh->Copy( newFaces[i] );
+  }
+  else
+  {
+    getEditor().ClearLastCreated();
+    SMESH_SequenceOfElemPtr& aSeq =
+      const_cast<SMESH_SequenceOfElemPtr&>( getEditor().GetLastCreatedElems() );
+    for ( size_t i = 0; i < newFaces.size(); ++i )
+      aSeq.Append( newFaces[i] );
+
+    TPythonDump() << this << ".FillHole( SMESH.FreeBorder(" << theHole.nodeIDs << " ))";
+  }
+
+  SMESH_CATCH( SMESH::throwCorbaException );
+}
+
 //=======================================================================
 //function : convError
 //purpose  :

src/SMESH_I/SMESH_MeshEditor_i.hxx

@@ -545,6 +545,30 @@ public:
   CORBA::Short GetPointState(CORBA::Double x, CORBA::Double y, CORBA::Double z)
     throw (SALOME::SALOME_Exception);
 
+  /*!
+   * Check if a 2D mesh is manifold
+   */
+  CORBA::Boolean IsManifold()
+    throw (SALOME::SALOME_Exception);
+
+  /*!
+   * Check if orientation of 2D elements is coherent
+   */
+  CORBA::Boolean IsCoherentOrientation2D()
+    throw (SALOME::SALOME_Exception);
+
+  /*!
+   * Returns all or only closed FreeBorder's.
+   */
+  SMESH::ListOfFreeBorders* FindFreeBorders(CORBA::Boolean closedOnly)
+    throw (SALOME::SALOME_Exception);
+
+  /*!
+   * Fill with 2D elements a hole defined by a FreeBorder.
+   */
+  void FillHole(const SMESH::FreeBorder& hole)
+    throw (SALOME::SALOME_Exception);
+
   SMESH::CoincidentFreeBorders* FindCoincidentFreeBorders(CORBA::Double tolerance);
   CORBA::Short SewCoincidentFreeBorders(const SMESH::CoincidentFreeBorders& freeBorders,
                                         CORBA::Boolean                      createPolygons,
@@ -552,35 +576,35 @@ public:
     throw (SALOME::SALOME_Exception);
 
   SMESH::SMESH_MeshEditor::Sew_Error
-  SewFreeBorders(CORBA::Long FirstNodeID1,
+    SewFreeBorders(CORBA::Long FirstNodeID1,
-                 CORBA::Long SecondNodeID1,
+                   CORBA::Long SecondNodeID1,
-                 CORBA::Long LastNodeID1,
+                   CORBA::Long LastNodeID1,
-                 CORBA::Long FirstNodeID2,
+                   CORBA::Long FirstNodeID2,
-                 CORBA::Long SecondNodeID2,
+                   CORBA::Long SecondNodeID2,
-                 CORBA::Long LastNodeID2,
+                   CORBA::Long LastNodeID2,
-                 CORBA::Boolean CreatePolygons,
+                   CORBA::Boolean CreatePolygons,
-                 CORBA::Boolean CreatePolyedrs) throw (SALOME::SALOME_Exception);
+                   CORBA::Boolean CreatePolyedrs) throw (SALOME::SALOME_Exception);
   SMESH::SMESH_MeshEditor::Sew_Error
-  SewConformFreeBorders(CORBA::Long FirstNodeID1,
+    SewConformFreeBorders(CORBA::Long FirstNodeID1,
-                        CORBA::Long SecondNodeID1,
+                          CORBA::Long SecondNodeID1,
-                        CORBA::Long LastNodeID1,
+                          CORBA::Long LastNodeID1,
-                        CORBA::Long FirstNodeID2,
+                          CORBA::Long FirstNodeID2,
-                        CORBA::Long SecondNodeID2) throw (SALOME::SALOME_Exception);
+                          CORBA::Long SecondNodeID2) throw (SALOME::SALOME_Exception);
   SMESH::SMESH_MeshEditor::Sew_Error
-  SewBorderToSide(CORBA::Long FirstNodeIDOnFreeBorder,
+    SewBorderToSide(CORBA::Long FirstNodeIDOnFreeBorder,
-                  CORBA::Long SecondNodeIDOnFreeBorder,
+                    CORBA::Long SecondNodeIDOnFreeBorder,
-                  CORBA::Long LastNodeIDOnFreeBorder,
+                    CORBA::Long LastNodeIDOnFreeBorder,
-                  CORBA::Long FirstNodeIDOnSide,
+                    CORBA::Long FirstNodeIDOnSide,
-                  CORBA::Long LastNodeIDOnSide,
+                    CORBA::Long LastNodeIDOnSide,
-                  CORBA::Boolean CreatePolygons,
+                    CORBA::Boolean CreatePolygons,
-                  CORBA::Boolean CreatePolyedrs) throw (SALOME::SALOME_Exception);
+                    CORBA::Boolean CreatePolyedrs) throw (SALOME::SALOME_Exception);
   SMESH::SMESH_MeshEditor::Sew_Error
-  SewSideElements(const SMESH::long_array& IDsOfSide1Elements,
+    SewSideElements(const SMESH::long_array& IDsOfSide1Elements,
-                  const SMESH::long_array& IDsOfSide2Elements,
+                    const SMESH::long_array& IDsOfSide2Elements,
-                  CORBA::Long NodeID1OfSide1ToMerge,
+                    CORBA::Long NodeID1OfSide1ToMerge,
-                  CORBA::Long NodeID1OfSide2ToMerge,
+                    CORBA::Long NodeID1OfSide2ToMerge,
-                  CORBA::Long NodeID2OfSide1ToMerge,
+                    CORBA::Long NodeID2OfSide1ToMerge,
-                  CORBA::Long NodeID2OfSide2ToMerge) throw (SALOME::SALOME_Exception);
+                    CORBA::Long NodeID2OfSide2ToMerge) throw (SALOME::SALOME_Exception);
 
   /*!
    * Set new nodes for given element.