// Copyright (C) 2018-2019  OPEN CASCADE
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, 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_Slot.cxx
// Created   : Fri Nov 30 15:58:37 2018
// Author    : Edward AGAPOV (eap)

#include "SMESH_MeshAlgos.hxx"

#include "ObjectPool.hxx"
#include "SMDS_LinearEdge.hxx"
#include "SMDS_Mesh.hxx"
#include "SMDS_MeshGroup.hxx"

#include <IntAna_IntConicQuad.hxx>
#include <IntAna_Quadric.hxx>
#include <NCollection_DataMap.hxx>
#include <NCollection_Map.hxx>
#include <Precision.hxx>
#include <gp_Ax1.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Dir.hxx>
#include <gp_Lin.hxx>
#include <gp_Pln.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>

#include <Utils_SALOME_Exception.hxx>

#include <boost/container/flat_set.hpp>

namespace
{
  typedef SMESH_MeshAlgos::Edge TEdge;

  //================================================================================
  //! point of intersection of a face edge with the cylinder
  struct IntPoint
  {
    SMESH_NodeXYZ myNode;        // point and a node
    int           myEdgeIndex;   // face edge index
    bool          myIsOutPln[2]; // isOut of two planes
  };

  //================================================================================
  //! poly-line segment
  struct Segment
  {
    typedef boost::container::flat_set< const SMDS_MeshNode* > TNodeSet;
    //typedef std::list< TEdge > TEdgeList;

    const SMDS_MeshElement* myEdge;
    TNodeSet                myEndNodes; // ends of cut edges
    //TEdgeList               myCutEdges[2];
    

    // return its axis
    gp_Ax1 Ax1( bool reversed = false ) const
    {
      SMESH_NodeXYZ n1 = myEdge->GetNode(  reversed );
      SMESH_NodeXYZ n2 = myEdge->GetNode( !reversed );
      return gp_Ax1( n1, gp_Dir( n2 - n1 ));
    }
    // return a node
    const SMDS_MeshNode* Node(int i) const
    {
      return myEdge->GetNode( i % 2 );
    }
    // store an intersection edge forming the slot border
    void AddEdge( TEdge& e, double tol )
    {
      const SMDS_MeshNode** nodes = & e._node1;
      for ( int i = 0; i < 2; ++i )
      {
        std::pair< TNodeSet::iterator, bool > nItAdded = myEndNodes.insert( nodes[ i ]);
        if ( !nItAdded.second )
          myEndNodes.erase( nItAdded.first );
      }
    }
    // { -- PREV version
    //   int i = myCutEdges[0].empty() ? 0 : 1;
    //   std::insert_iterator< TEdgeList > where = inserter( myCutEdges[i], myCutEdges[i].begin() );

    //   //double minDist = 1e100;
    //   SMESH_NodeXYZ nNew[2] = { e._node1, e._node2 };
    //   int iNewMin = 0, iCurMin = 1;
    //   for ( i = 0; i < 2; ++i )
    //   {
    //     if ( myCutEdges[i].empty() )
    //       continue;
    //     SMESH_NodeXYZ nCur[2] = { myCutEdges[i].front()._node1,
    //                               myCutEdges[i].back()._node2 };
    //     for   ( int iN = 0; iN < 2; ++iN )
    //       for ( int iC = 0; iC < 2; ++iC )
    //       {
    //         if (( nCur[iC].Node() && nCur[iC] == nNew[iN] ) ||
    //             ( nCur[iC] - nNew[iN] ).SquareModulus() < tol * tol )
    //         {
    //           where   = inserter( myCutEdges[i], iC ? myCutEdges[i].end() : myCutEdges[i].begin() );
    //           iNewMin = iN;
    //           iCurMin = iC;
    //           //minDist = dist;
    //           iN = 2;
    //           break;
    //         }
    //       }
    //   }
    //   if ( iNewMin == iCurMin )
    //     std::swap( e._node1, e._node2 );

    //   where = e;
    // }
    Segment( const SMDS_MeshElement* e = 0 ): myEdge(e) { myEndNodes.reserve( 4 ); }
  };
  typedef ObjectPoolIterator<Segment> TSegmentIterator;

  
  //================================================================================
  /*!
   * \brief Intersect a face edge given by its nodes with a cylinder.
   */
  //================================================================================

  void intersectEdge( const gp_Cylinder&      cyl,
                      const SMESH_NodeXYZ&    n1,
                      const SMESH_NodeXYZ&    n2,
                      const double            tol,
                      std::vector< IntPoint >& intPoints )
  {
    gp_Lin line( gp_Ax1( n1, gp_Dir( n2 - n1 )));
    IntAna_IntConicQuad intersection( line, IntAna_Quadric( cyl ));

    if ( !intersection.IsDone()     ||
         intersection.IsParallel()  ||
         intersection.IsInQuadric() ||
         intersection.NbPoints() == 0 )
      return;

    gp_Vec edge( n1, n2 );

    size_t oldNbPnts = intPoints.size();
    for ( int iP = 1; iP <= intersection.NbPoints(); ++iP )
    {
      const gp_Pnt& p = intersection.Point( iP );

      gp_Vec n1p ( n1, p );
      const SMDS_MeshNode* n = 0;

      double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
      if ( u <= 0. ) {
        if ( p.SquareDistance( n1 ) < tol * tol )
          n = n1.Node();
        else
          continue;
      }
      else if ( u >= 1. ) {
        if ( p.SquareDistance( n2 ) < tol * tol )
          n = n2.Node();
        else
          continue;
      }
      else {
        if      ( p.SquareDistance( n1 ) < tol * tol )
          n = n1.Node();
        else if ( p.SquareDistance( n2 ) < tol * tol )
          n = n2.Node();
      }

      intPoints.push_back( IntPoint() );
      if ( n )
        intPoints.back().myNode.Set( n );
      else
        intPoints.back().myNode.SetCoord( p.X(),p.Y(),p.Z() );
    }

    // set points order along an edge
    if ( intPoints.size() - oldNbPnts == 2 &&
         intersection.ParamOnConic( 1 ) > intersection.ParamOnConic( 2 ))
    {
      int i = intPoints.size() - 1;
      std::swap( intPoints[ i ], intPoints[ i - 1 ]);
    }

    return;
  }

  //================================================================================
  /*!
   * \brief Return signed distance between a point and a plane
   */
  //================================================================================

  double signedDist( const gp_Pnt& p, const gp_Ax1& planeNormal )
  {
    const gp_Pnt& O = planeNormal.Location();
    gp_Vec Op( O, p );
    return Op * planeNormal.Direction();
  }

  //================================================================================
  /*!
   * \brief Check if a point is outside a segment domain bound by two planes
   */
  //================================================================================

  bool isOut( const gp_Pnt& p, const gp_Ax1* planeNormal, bool* isOutPtr )
  {
    isOutPtr[0] = isOutPtr[1] = false;

    for ( int i = 0; i < 2; ++i )
    {
      isOutPtr[i] = ( signedDist( p, planeNormal[i] ) <= 0. );
    }
    return ( isOutPtr[0] && isOutPtr[1] );
  }

  //================================================================================
  /*!
   * \brief Check if a segment between two points is outside a segment domain bound by two planes
   */
  //================================================================================

  bool isSegmentOut( bool* isOutPtr1, bool* isOutPtr2 )
  {
    return (( isOutPtr1[0] && isOutPtr2[0] ) ||
            ( isOutPtr1[1] && isOutPtr2[1] ));
  }

  //================================================================================
  /*!
   * \brief cut off ip1 from edge (ip1 - ip2) by a plane
   */
  //================================================================================

  void cutOff( IntPoint & ip1, const IntPoint & ip2, const gp_Ax1& planeNormal, double tol )
  {
    gp_Lin lin( ip1.myNode, ( ip2.myNode - ip1.myNode ));
    gp_Pln pln( planeNormal.Location(), planeNormal.Direction() );

    IntAna_IntConicQuad intersection( lin, pln, Precision::Angular/*Tolerance*/() );
    if ( intersection.IsDone()      &&
         !intersection.IsParallel()  &&
         !intersection.IsInQuadric() &&
         intersection.NbPoints() == 1 )
    {
      if ( intersection.Point( 1 ).SquareDistance( ip1.myNode ) > tol * tol )
      {
        static_cast< gp_XYZ& >( ip1.myNode ) = intersection.Point( 1 ).XYZ();
        ip1.myNode._node = 0;
        ip1.myEdgeIndex = -1;
      }
    }
  }

  //================================================================================
  /*!
   * \brief Assure that face normal is computed in faceNormals vector 
   */
  //================================================================================

  const gp_XYZ& computeNormal( const SMDS_MeshElement* face,
                               std::vector< gp_XYZ >&  faceNormals )
  {
    bool toCompute;
    if ((int) faceNormals.size() <= face->GetID() )
    {
      toCompute = true;
      faceNormals.resize( face->GetID() + 1 );
    }
    else
    {
      toCompute = faceNormals[ face->GetID() ].SquareModulus() == 0.;
    }
    if ( toCompute )
      SMESH_MeshAlgos::FaceNormal( face, faceNormals[ face->GetID() ], /*normalized=*/false );

    return faceNormals[ face->GetID() ];
  }

  typedef std::vector< SMDS_MeshGroup* > TGroupVec;

  //================================================================================
  /*!
   * \brief Fill theFaceID2Groups map for a given face
   *  \param [in] theFace - the face
   *  \param [in] theGroupsToUpdate - list of groups to treat
   *  \param [out] theFaceID2Groups - the map to fill in
   *  \param [out] theWorkGroups - a working buffer of groups
   */
  //================================================================================

  void findGroups( const SMDS_MeshElement *                theFace,
                   TGroupVec &                             theGroupsToUpdate,
                   NCollection_DataMap< int, TGroupVec > & theFaceID2Groups,
                   TGroupVec &                             theWorkGroups )
  {
    theWorkGroups.clear();
    for ( size_t i = 0; i < theGroupsToUpdate.size(); ++i )
      if ( theGroupsToUpdate[i]->Contains( theFace ))
        theWorkGroups.push_back( theGroupsToUpdate[i] );

    if ( !theWorkGroups.empty() )
      theFaceID2Groups.Bind( theFace->GetID(), theWorkGroups );
  }
}

//================================================================================
/*!
 * \brief Create a slot of given width around given 1D elements lying on a triangle mesh.
 * The slot is consrtucted by cutting faces by cylindrical surfaces made around each segment.
 * \return Edges located at the slot boundary
 */
//================================================================================

std::vector< SMESH_MeshAlgos::Edge >
SMESH_MeshAlgos::MakeSlot( SMDS_ElemIteratorPtr             theSegmentIt,
                           double                           theWidth,
                           SMDS_Mesh*                       theMesh,
                           std::vector< SMDS_MeshGroup* > & theGroupsToUpdate)
{
  std::vector< Edge > bndEdges;

  if ( !theSegmentIt || !theSegmentIt->more() || !theMesh || theWidth == 0.)
    return bndEdges;

  // put the input segments to a data map in order to be able finding neighboring ones

  typedef std::vector< Segment* >                                                TSegmentVec;
  typedef NCollection_DataMap< const SMDS_MeshNode*, TSegmentVec, SMESH_Hasher > TSegmentsOfNode;
  TSegmentsOfNode segmentsOfNode;
  ObjectPool< Segment > segmentPool;

  while( theSegmentIt->more() )
  {
    const SMDS_MeshElement* edge = theSegmentIt->next();
    if ( edge->GetType() != SMDSAbs_Edge )
      throw SALOME_Exception( "A segment is not a mesh edge");

    Segment* segment = segmentPool.getNew();
    segment->myEdge = edge;

    for ( SMDS_NodeIteratorPtr nIt = edge->nodeIterator(); nIt->more(); )
    {
      const SMDS_MeshNode* n = nIt->next();
      TSegmentVec* segVec = segmentsOfNode.ChangeSeek( n );
      if ( !segVec )
        segVec = segmentsOfNode.Bound( n, TSegmentVec() );
      segVec->reserve(2);
      segVec->push_back( segment );
    }
  }

  // Cut the mesh around the segments

  const double tol = Precision::Confusion();
  std::vector< gp_XYZ > faceNormals;
  SMESH_MeshAlgos::Intersector meshIntersector( theMesh, tol, faceNormals );
  std::unique_ptr< SMESH_ElementSearcher> faceSearcher;

  std::vector< NLink > startEdges;
  std::vector< const SMDS_MeshNode* > faceNodes(4), edgeNodes(2);
  std::vector<const SMDS_MeshElement *> faces(2);
  NCollection_Map<const SMDS_MeshElement*, SMESH_Hasher > checkedFaces;
  std::vector< IntPoint > intPoints, p(2);
  std::vector< SMESH_NodeXYZ > facePoints(4);
  std::vector< Intersector::TFace > cutFacePoints;

  NCollection_DataMap< int, TGroupVec > faceID2Groups;
  TGroupVec groupVec;

  std::vector< gp_Ax1 > planeNormalVec(2);
  gp_Ax1 * planeNormal = & planeNormalVec[0];
  
  for ( TSegmentIterator segIt( segmentPool ); segIt.more(); ) // loop on all segments
  {
    Segment* segment = const_cast< Segment* >( segIt.next() );

    gp_Lin      segLine( segment->Ax1() );
    gp_Ax3      cylAxis( segLine.Location(), segLine.Direction() );
    gp_Cylinder segCylinder( cylAxis, 0.5 * theWidth );
    double      radius2( segCylinder.Radius() * segCylinder.Radius() );

    // get normals of planes separating domains of neighboring segments
    for ( int i = 0; i < 2; ++i ) // loop on 2 segment ends
    {
      planeNormal[i] = segment->Ax1(i);

      const SMDS_MeshNode*    n = segment->Node( i );
      const TSegmentVec& segVec = segmentsOfNode( n );
      for ( size_t iS = 0; iS < segVec.size(); ++iS )
      {
        if ( segVec[iS] == segment )
          continue;

        gp_Ax1 axis2 = segVec[iS]->Ax1();
        if ( n != segVec[iS]->Node( 1 ))
          axis2.Reverse(); // along a wire

        planeNormal[i].SetDirection( planeNormal[i].Direction().XYZ() + axis2.Direction().XYZ() );
      }
    }

    // we explore faces around a segment starting from face edges;
    // initialize a list of starting edges
    startEdges.clear();
    {
      // get a face to start searching intersected faces from
      const SMDS_MeshNode*      n0 = segment->Node( 0 );
      SMDS_ElemIteratorPtr     fIt = n0->GetInverseElementIterator( SMDSAbs_Face );
      const SMDS_MeshElement* face = ( fIt->more() ) ? fIt->next() : 0;
      if ( !theMesh->Contains( face ))
      {
        if ( !faceSearcher )
          faceSearcher.reset( SMESH_MeshAlgos::GetElementSearcher( *theMesh ));
        face = faceSearcher->FindClosestTo( SMESH_NodeXYZ( n0 ), SMDSAbs_Face );
      }
      // collect face edges
      int nbNodes = face->NbCornerNodes();
      faceNodes.assign( face->begin_nodes(), face->end_nodes() );
      faceNodes.resize( nbNodes + 1 );
      faceNodes[ nbNodes ] = faceNodes[ 0 ];
      for ( int i = 0; i < nbNodes; ++i )
        startEdges.push_back( NLink( faceNodes[i], faceNodes[i+1] ));
    }

    // intersect faces located around a segment
    checkedFaces.Clear();
    while ( !startEdges.empty() )
    {
      edgeNodes[0] = startEdges[0].first;
      edgeNodes[1] = startEdges[0].second;

      theMesh->GetElementsByNodes( edgeNodes, faces, SMDSAbs_Face );
      for ( size_t iF = 0; iF < faces.size(); ++iF ) // loop on faces sharing a start edge
      {
        const SMDS_MeshElement* face = faces[iF];
        if ( !checkedFaces.Add( face ))
          continue;

        int nbNodes = face->NbCornerNodes();
        if ( nbNodes != 3 )
          throw SALOME_Exception( "MakeSlot() accepts triangles only" );
        facePoints.assign( face->begin_nodes(), face->end_nodes() );
        facePoints.resize( nbNodes + 1 );
        facePoints[ nbNodes ] = facePoints[ 0 ];

        // check if cylinder axis || face        
        const gp_XYZ& faceNorm = computeNormal( face, faceNormals );
        bool isCylinderOnFace  = ( Abs( faceNorm * cylAxis.Direction().XYZ() ) < tol );

        if ( !isCylinderOnFace )
        {
          if ( Intersector::CutByPlanes( face, planeNormalVec, tol, cutFacePoints ))
            continue; // whole face cut off
          facePoints.swap( cutFacePoints[0] );
          facePoints.push_back( facePoints[0] );
        }

        // find intersection points on face edges
        intPoints.clear();
        int nbPoints = facePoints.size()-1;
        int nbFarPoints = 0;
        for ( int i = 0; i < nbPoints; ++i )
        {
          const SMESH_NodeXYZ& n1 = facePoints[i];
          const SMESH_NodeXYZ& n2 = facePoints[i+1];

          size_t iP = intPoints.size();
          intersectEdge( segCylinder, n1, n2, tol, intPoints );

          // save edge index
          if ( isCylinderOnFace )
            for ( ; iP < intPoints.size(); ++iP )
              intPoints[ iP ].myEdgeIndex = i;
          else
            for ( ; iP < intPoints.size(); ++iP )
              if ( n1.Node() && n2.Node() )
                intPoints[ iP ].myEdgeIndex = face->GetNodeIndex( n1.Node() );
              else
                intPoints[ iP ].myEdgeIndex = -(i+1);

          nbFarPoints += ( segLine.SquareDistance( n1 ) > radius2 );
        }

        // feed startEdges
        if ( nbFarPoints < nbPoints || !intPoints.empty() )
          for ( int i = 0; i < nbPoints; ++i )
          {
            const SMESH_NodeXYZ& n1 = facePoints[i];
            const SMESH_NodeXYZ& n2 = facePoints[i+1];
            if ( n1.Node() && n2.Node() )
            {
              isOut( n1, planeNormal, p[0].myIsOutPln );
              isOut( n2, planeNormal, p[1].myIsOutPln );
              if ( !isSegmentOut( p[0].myIsOutPln, p[1].myIsOutPln ))
              {
                startEdges.push_back( NLink( n1.Node(), n2.Node() ));
              }
            }
          }

        if ( intPoints.size() < 2 )
          continue;

        // classify intPoints by planes
        for ( size_t i = 0; i < intPoints.size(); ++i )
          isOut( intPoints[i].myNode, planeNormal, intPoints[i].myIsOutPln );

        // cut the face

        if ( intPoints.size() > 2 )
          intPoints.push_back( intPoints[0] );

        for ( size_t iE = 1; iE < intPoints.size(); ++iE ) // 2 <= intPoints.size() <= 5
        {
          if (( intPoints[iE].myIsOutPln[0] && intPoints[iE].myIsOutPln[1]   ) ||
              ( isSegmentOut( intPoints[iE].myIsOutPln, intPoints[iE-1].myIsOutPln )))
            continue; // intPoint is out of domain

          // check if a cutting edge connecting two intPoints is on cylinder surface
          if ( intPoints[iE].myEdgeIndex == intPoints[iE-1].myEdgeIndex )
            continue; // on same edge
          if ( intPoints[iE].myNode.Node() &&
               intPoints[iE].myNode == intPoints[iE-1].myNode ) // coincide
            continue;

          gp_XYZ edegDir = intPoints[iE].myNode - intPoints[iE-1].myNode;

          bool toCut; // = edegDir.SquareModulus() > tol * tol;
          if ( intPoints.size() == 2 )
            toCut = true;
          else if ( isCylinderOnFace )
            toCut = cylAxis.Direction().IsParallel( edegDir, tol );
          else
          {
            SMESH_NodeXYZ nBetween;
            int eInd = intPoints[iE-1].myEdgeIndex;
            if ( eInd < 0 )
              nBetween = facePoints[( 1 - (eInd-1)) % nbPoints ];
            else
              nBetween = faceNodes[( 1 + eInd ) % nbNodes ];
            toCut = ( segLine.SquareDistance( nBetween ) > radius2 );
          }
          if ( !toCut )
            continue;

          // limit the edge by planes
          if ( intPoints[iE].myIsOutPln[0] ||
               intPoints[iE].myIsOutPln[1] )
            cutOff( intPoints[iE], intPoints[iE-1],
                    planeNormal[ intPoints[iE].myIsOutPln[1] ], tol );

          if ( intPoints[iE-1].myIsOutPln[0] ||
               intPoints[iE-1].myIsOutPln[1] )
            cutOff( intPoints[iE-1], intPoints[iE],
                    planeNormal[ intPoints[iE-1].myIsOutPln[1] ], tol );

          edegDir = intPoints[iE].myNode - intPoints[iE-1].myNode;
          if ( edegDir.SquareModulus() < tol * tol )
            continue; // fully cut off

          // face cut
          meshIntersector.Cut( face,
                               intPoints[iE-1].myNode, intPoints[iE-1].myEdgeIndex,
                               intPoints[iE  ].myNode, intPoints[iE  ].myEdgeIndex );

          Edge e = { intPoints[iE].myNode.Node(), intPoints[iE-1].myNode.Node(), 0 };
          segment->AddEdge( e, tol );
          bndEdges.push_back( e );

          findGroups( face, theGroupsToUpdate, faceID2Groups, groupVec );

        }
      }  // loop on faces sharing an edge

      startEdges[0] = startEdges.back();
      startEdges.pop_back();

    } // loop on startEdges
  } // loop on all input segments


  // Make cut at the end of group of segments

  std::vector<const SMDS_MeshElement*> polySegments;

  for ( TSegmentsOfNode::Iterator nSegsIt( segmentsOfNode ); nSegsIt.More(); nSegsIt.Next() )
  {
    const TSegmentVec& segVec = nSegsIt.Value();
    if ( segVec.size() != 1 )
      continue;

    const Segment*       segment = segVec[0];
    const SMDS_MeshNode* segNode = nSegsIt.Key();

    // find two end nodes of cut edges to make a cut between
    if ( segment->myEndNodes.size() != 4 )
      throw SALOME_Exception( "MakeSlot(): too short end edge?" );
    SMESH_MeshAlgos::PolySegment linkNodes;
    gp_Ax1 planeNorm = segment->Ax1( segNode != segment->Node(0) );
    double minDist[2] = { 1e100, 1e100 };
    Segment::TNodeSet::const_iterator nIt = segment->myEndNodes.begin();
    for ( ; nIt != segment->myEndNodes.end(); ++nIt )
    {
      SMESH_NodeXYZ n = *nIt;
      double d = Abs( signedDist( n, planeNorm ));
      double diff1 = minDist[0] - d, diff2 = minDist[1] - d;
      int i;
      if ( diff1 > 0 && diff2 > 0 )
      {
        i = ( diff1 < diff2 );
      }
      else if ( diff1 > 0 )
      {
        i = 0;
      }
      else if ( diff2 > 0 )
      {
        i = 1;
      }
      else
      {
        continue;
      }
      linkNodes.myXYZ[ i ] = n;
      minDist        [ i ] = d;
    }
    // for ( int iSide = 0; iSide < 2; ++iSide )
    // {
    //   if ( segment->myCutEdges[ iSide ].empty() )
    //     throw SALOME_Exception( "MakeSlot(): too short end edge?" );
    //   SMESH_NodeXYZ n1 = segment->myCutEdges[ iSide ].front()._node1;
    //   SMESH_NodeXYZ n2 = segment->myCutEdges[ iSide ].back ()._node2;
    //   double d1 = Abs( signedDist( n1, planeNorm ));
    //   double d2 = Abs( signedDist( n2, planeNorm ));
    //   linkNodes.myXYZ  [ iSide ] = ( d1 < d2 ) ? n1 : n2;
    //   linkNodes.myNode1[ iSide ] = linkNodes.myNode2[ iSide ] = 0;
    // }
    linkNodes.myVector = planeNorm.Direction() ^ (linkNodes.myXYZ[0] - linkNodes.myXYZ[1]);
    linkNodes.myNode1[ 0 ] = linkNodes.myNode2[ 0 ] = 0;
    linkNodes.myNode1[ 1 ] = linkNodes.myNode2[ 1 ] = 0;

    // create segments connecting linkNodes
    std::vector<const SMDS_MeshElement*> newSegments;
    std::vector<const SMDS_MeshNode*>    newNodes;
    SMESH_MeshAlgos::TListOfPolySegments polySegs(1, linkNodes);
    SMESH_MeshAlgos::MakePolyLine( theMesh, polySegs, newSegments, newNodes,
                                   /*group=*/0, faceSearcher.get() );
    // cut faces by newSegments
    intPoints.resize(2);
    for ( size_t i = 0; i < newSegments.size(); ++i )
    {
      intPoints[0].myNode = edgeNodes[0] = newSegments[i]->GetNode(0);
      intPoints[1].myNode = edgeNodes[1] = newSegments[i]->GetNode(1);

      // find an underlying face
      gp_XYZ                middle = 0.5 * ( intPoints[0].myNode + intPoints[1].myNode );
      const SMDS_MeshElement* face = faceSearcher->FindClosestTo( middle, SMDSAbs_Face );

      // find intersected edges of the face
      int nbNodes = face->NbCornerNodes();
      faceNodes.assign( face->begin_nodes(), face->end_nodes() );
      faceNodes.resize( nbNodes + 1 );
      faceNodes[ nbNodes ] = faceNodes[ 0 ];
      for ( int iP = 0; iP < 2; ++iP )
      {
        intPoints[iP].myEdgeIndex = -1;
        for ( int iN = 0; iN < nbNodes &&  intPoints[iP].myEdgeIndex < 0; ++iN )
        {
          SMDS_LinearEdge edge( faceNodes[iN], faceNodes[iN+1] );
          if ( SMESH_MeshAlgos::GetDistance( &edge, intPoints[iP].myNode) < tol )
            intPoints[iP].myEdgeIndex = iN;
        }
      }


      // face cut
      computeNormal( face, faceNormals );
      meshIntersector.Cut( face,
                           intPoints[0].myNode, intPoints[0].myEdgeIndex,
                           intPoints[1].myNode, intPoints[1].myEdgeIndex );

      Edge e = { intPoints[0].myNode.Node(), intPoints[1].myNode.Node(), 0 };
      bndEdges.push_back( e );

      findGroups( face, theGroupsToUpdate, faceID2Groups, groupVec );

      // add cut points to an adjacent face at ends of poly-line
      // if they fall onto face edges
      if (( i == 0                       && intPoints[0].myEdgeIndex >= 0 ) ||
          ( i == newSegments.size() - 1  && intPoints[1].myEdgeIndex >= 0 ))
      {
        for ( int iE = 0; iE < 2; ++iE ) // loop on ends of a new segment
        {
          if ( iE ? ( i != newSegments.size() - 1 ) : ( i != 0 ))
            continue;
          int iEdge = intPoints[ iE ].myEdgeIndex;
          edgeNodes[0] = faceNodes[ iEdge ];
          edgeNodes[1] = faceNodes[ iEdge+1 ];
          theMesh->GetElementsByNodes( edgeNodes, faces, SMDSAbs_Face );
          for ( size_t iF = 0; iF < faces.size(); ++iF )
            if ( faces[iF] != face )
            {
              int iN1 = faces[iF]->GetNodeIndex( edgeNodes[0] );
              int iN2 = faces[iF]->GetNodeIndex( edgeNodes[1] );
              intPoints[ iE ].myEdgeIndex = Abs( iN1 - iN2 ) == 1 ? Min( iN1, iN2 ) : 2;
              computeNormal( faces[iF], faceNormals );
              meshIntersector.Cut( faces[iF],
                                   intPoints[iE].myNode, intPoints[iE].myEdgeIndex,
                                   intPoints[iE].myNode, intPoints[iE].myEdgeIndex );

              findGroups( faces[iF], theGroupsToUpdate, faceID2Groups, groupVec );
            }
        }
      }

    } // loop on newSegments

    polySegments.insert( polySegments.end(), newSegments.begin(), newSegments.end() );

  } // loop on map of input segments

  // actual mesh splitting
  TElemIntPairVec new2OldFaces;
  TNodeIntPairVec new2OldNodes;
  meshIntersector.MakeNewFaces( new2OldFaces, new2OldNodes, /*sign=*/1, /*optimize=*/true );

  // add new faces to theGroupsToUpdate
  for ( size_t i = 0; i < new2OldFaces.size(); ++i )
  {
    const SMDS_MeshElement* newFace = new2OldFaces[i].first;
    const int             oldFaceID = new2OldFaces[i].second;
    if ( !newFace ) continue;

    if ( TGroupVec* groups = const_cast< TGroupVec* >( faceID2Groups.Seek( oldFaceID )))
      for ( size_t iG = 0; iG < groups->size(); ++iG )
        (*groups)[ iG ]->Add( newFace );
  }

  // remove poly-line edges
  for ( size_t i = 0; i < polySegments.size(); ++i )
  {
    edgeNodes[0] = polySegments[i]->GetNode(0);
    edgeNodes[1] = polySegments[i]->GetNode(1);

    theMesh->RemoveFreeElement( polySegments[i] );

    if ( edgeNodes[0]->NbInverseElements() == 0 )
      theMesh->RemoveNode( edgeNodes[0] );
    if ( edgeNodes[1]->NbInverseElements() == 0 )
      theMesh->RemoveNode( edgeNodes[1] );
  }

  return bndEdges;
}