Salome HOME
GPUSPHGUI: Offset transformation
authoreap <eap@opencascade.com>
Tue, 13 Feb 2018 12:16:32 +0000 (15:16 +0300)
committereap <eap@opencascade.com>
Tue, 13 Feb 2018 12:16:32 +0000 (15:16 +0300)
add missing files

src/SMESHGUI/SMESHGUI_OffsetDlg.h [new file with mode: 0644]
src/SMESHUtils/SMESH_MeshAlgos.cxx
src/SMESHUtils/SMESH_Offset.cxx [new file with mode: 0644]
src/SMESHUtils/SMESH_Triangulate.cxx [new file with mode: 0644]
src/SMESHUtils/SMESH_Triangulate.hxx [new file with mode: 0644]

diff --git a/src/SMESHGUI/SMESHGUI_OffsetDlg.h b/src/SMESHGUI/SMESHGUI_OffsetDlg.h
new file mode 100644 (file)
index 0000000..e8f4b9c
--- /dev/null
@@ -0,0 +1,134 @@
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// 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
+//
+
+// SMESH SMESHGUI : GUI for SMESH component
+// File   : SMESHGUI_OffsetDlg.h
+//
+#ifndef SMESHGUI_OffsetDLG_H
+#define SMESHGUI_OffsetDLG_H
+
+// SMESH includes
+#include "SMESH_SMESHGUI.hxx"
+#include "SMESHGUI_PreviewDlg.h"
+
+// IDL includes
+#include <SALOMEconfig.h>
+#include CORBA_SERVER_HEADER(SMESH_Mesh)
+
+class QButtonGroup;
+class QGroupBox;
+class QLabel;
+class QLineEdit;
+class QPushButton;
+class QRadioButton;
+class QCheckBox;
+class SMESHGUI;
+class SMESHGUI_IdValidator;
+class SMESHGUI_SpinBox;
+class SMESHGUI_FilterDlg;
+class SMESH_Actor;
+class SVTK_Selector;
+class LightApp_SelectionMgr;
+class SMESH_LogicalFilter;
+
+//=================================================================================
+// class    : SMESHGUI_OffsetDlg
+// purpose  :
+//=================================================================================
+class SMESHGUI_EXPORT SMESHGUI_OffsetDlg : public SMESHGUI_MultiPreviewDlg
+{
+  Q_OBJECT
+
+ public:
+  SMESHGUI_OffsetDlg( SMESHGUI* );
+  ~SMESHGUI_OffsetDlg();
+
+ private:
+  void                   Init( bool = true );
+  void                   enterEvent( QEvent* );           /* mouse enter the QWidget */
+  void                   keyPressEvent( QKeyEvent* );
+  int                    GetConstructorId();
+  void                   setNewMeshName();
+
+  bool                   isValid();
+  void                   getOffset( SMESH::PointStruct& thePoint,
+                                    SMESH::double_array_var& theOffsetFact);
+
+  SMESHGUI_IdValidator*  myIdValidator;
+  LightApp_SelectionMgr* mySelectionMgr;          /* User shape selection */
+  QString                myElementsId;
+  int                    myNbOkElements;          /* to check when elements are defined */
+
+  SVTK_Selector*         mySelector;
+
+  QWidget*               myEditCurrentArgument;
+
+  bool                   myBusy;
+  SMESH_Actor*           myActor;
+  SMESH_LogicalFilter*   myMeshOrSubMeshOrGroupFilter;
+
+  QList<SMESH::SMESH_IDSource_var> myObjects;
+  QList<QString>                   myObjectsNames;
+  QList<SMESH::SMESH_Mesh_var>     myMeshes;
+
+
+  QGroupBox*             ConstructorsBox;
+  QGroupBox*             GroupButtons;
+  QPushButton*           buttonOk;
+  QPushButton*           buttonCancel;
+  QPushButton*           buttonApply;
+  QPushButton*           buttonHelp;
+  QGroupBox*             GroupArguments;
+  QLabel*                TextLabelElements;
+  //QPushButton*           SelectElementsButton;
+  QLineEdit*             LineEditElements;
+  QCheckBox*             CheckBoxMesh;
+  SMESHGUI_SpinBox*      SpinBox;
+  QGroupBox*             ActionBox;
+  QButtonGroup*          ActionGroup;
+  QCheckBox*             MakeGroupsCheck;
+  QLineEdit*             LineEditNewMesh;
+
+  QString                myHelpFileName;
+
+  QPushButton*           myFilterBtn;
+  SMESHGUI_FilterDlg*    myFilterDlg;
+
+
+protected slots:
+  virtual void           onDisplaySimulation( bool );
+  virtual void           reject();
+  void                   onFilterAccepted();
+   
+private slots:
+  void                   ClickOnOk();
+  bool                   ClickOnApply();
+  void                   ClickOnHelp();
+  void                   SelectionIntoArgument();
+  void                   DeactivateActiveDialog();
+  void                   ActivateThisDialog();
+  void                   onTextChange( const QString& );
+  void                   onSelectMesh( bool );
+  void                   onActionClicked( int );
+  void                   onOpenView();
+  void                   onCloseView();
+  void                   setFilters();
+};
+
+#endif // SMESHGUI_OffsetDLG_H
index 1de6ba8..fa65474 100644 (file)
@@ -258,7 +258,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     //!< allocator of ElementBox's and SMESH_TreeLimit
     struct LimitAndPool : public SMESH_TreeLimit
     {
-      TElementBoxPool            _elBoPool;
+      TElementBoxPool _elBoPool;
       LimitAndPool():SMESH_TreeLimit( MaxLevel, /*minSize=*/0. ) {}
     };
     LimitAndPool* getLimitAndPool() const
diff --git a/src/SMESHUtils/SMESH_Offset.cxx b/src/SMESHUtils/SMESH_Offset.cxx
new file mode 100644 (file)
index 0000000..f571535
--- /dev/null
@@ -0,0 +1,2912 @@
+// 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_Offset.cxx
+// Created   : Mon Dec 25 15:52:38 2017
+// Author    : Edward AGAPOV (eap)
+
+#include "SMESH_MeshAlgos.hxx"
+
+#include <SMDS_PolygonalFaceOfNodes.hxx>
+#include "SMDS_Mesh.hxx"
+
+#include <Utils_SALOME_Exception.hxx>
+
+#include <Bnd_B3d.hxx>
+#include <NCollection_Map.hxx>
+#include <gp_Lin.hxx>
+#include <gp_Pln.hxx>
+
+#include <boost/container/flat_set.hpp>
+#include <boost/dynamic_bitset.hpp>
+
+namespace
+{
+  const size_t theMaxNbFaces = 256; // max number of faces sharing a node
+
+  typedef NCollection_DataMap< Standard_Address, const SMDS_MeshNode* > TNNMap;
+  typedef NCollection_Map< SMESH_Link, SMESH_Link >                     TLinkMap;
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Intersected face side storing a node created at this intersection
+   *        and a intersected face
+   */
+  struct CutLink
+  {
+    bool                     myReverse;
+    const SMDS_MeshNode*     myNode[2]; // side nodes
+    mutable SMESH_NodeXYZ    myIntNode; // intersection node
+    const SMDS_MeshElement*  myFace;    // cutter face
+    int                      myIndex;   // index of a node on the same link
+
+    CutLink(const SMDS_MeshNode*    node1=0,
+            const SMDS_MeshNode*    node2=0,
+            const SMDS_MeshElement* face=0,
+            const int               index=0) { Set ( node1, node2, face, index ); }
+
+    void Set( const SMDS_MeshNode*    node1,
+              const SMDS_MeshNode*    node2,
+              const SMDS_MeshElement* face,
+              const int               index=0)
+    {
+      myNode[0] = node1; myNode[1] = node2; myFace = face; myIndex = index; myReverse = false;
+      if ( myNode[0] && ( myReverse = ( myNode[0]->GetID() > myNode[1]->GetID() )))
+        std::swap( myNode[0], myNode[1] );
+    }
+    const SMDS_MeshNode* IntNode() const { return myIntNode.Node(); }
+    const SMDS_MeshNode* Node1() const { return myNode[ myReverse ]; }
+    const SMDS_MeshNode* Node2() const { return myNode[ !myReverse ]; }
+
+    static Standard_Integer HashCode(const CutLink&         link,
+                                     const Standard_Integer upper)
+    {
+      Standard_Integer n = ( link.myNode[0]->GetID() +
+                             link.myNode[1]->GetID() +
+                             link.myIndex );
+      return ::HashCode( n, upper );
+    }
+    static Standard_Boolean IsEqual(const CutLink& link1, const CutLink& link2 )
+    {
+      return ( link1.myNode[0] == link2.myNode[0] &&
+               link1.myNode[1] == link2.myNode[1] &&
+               link1.myIndex == link2.myIndex );
+    }
+  };
+
+  typedef NCollection_Map< CutLink, CutLink > TCutLinkMap;
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Part of a divided face edge
+   */
+  struct EdgePart
+  {
+    const SMDS_MeshNode*    myNode1;
+    const SMDS_MeshNode*    myNode2;
+    int                     myIndex; // positive -> side index, negative -> State
+    const SMDS_MeshElement* myFace;
+
+    enum State { _INTERNAL = -1, _COPLANAR = -2 };
+
+    void Set( const SMDS_MeshNode*    Node1,
+              const SMDS_MeshNode*    Node2,
+              const SMDS_MeshElement* Face = 0,
+              int                     EdgeIndex = _INTERNAL )
+    { myNode1 = Node1; myNode2 = Node2; myIndex = EdgeIndex; myFace = Face; }
+
+    // bool HasSameNode( const EdgePart& other ) { return ( myNode1 == other.myNode1 ||
+    //                                                      myNode1 == other.myNode2 ||
+    //                                                      myNode2 == other.myNode1 ||
+    //                                                      myNode2 == other.myNode2 );
+    // }
+    bool IsInternal() const { return myIndex < 0; }
+    bool IsTwin( const EdgePart& e ) const { return myNode1 == e.myNode2 && myNode2 == e.myNode1; }
+    bool IsSame( const EdgePart& e ) const {
+      return (( myNode1 == e.myNode2 && myNode2 == e.myNode1 ) ||
+              ( myNode1 == e.myNode1 && myNode2 == e.myNode2 )); }
+    bool ReplaceCoplanar( const EdgePart& e );
+    operator SMESH_Link() const { return SMESH_Link( myNode1, myNode2 ); }
+    operator gp_Vec() const { return SMESH_NodeXYZ( myNode2 ) - SMESH_NodeXYZ( myNode1 ); }
+  };
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Loop of EdgePart's forming a new face which is a part of CutFace
+   */
+  struct EdgeLoop : public SMDS_PolygonalFaceOfNodes
+  {
+    std::vector< const EdgePart* > myLinks;
+    bool                           myIsBndConnected; //!< is there a path to CutFace side edges
+    bool                           myHasPending;     //!< an edge encounters twice
+
+    EdgeLoop() : SMDS_PolygonalFaceOfNodes( std::vector<const SMDS_MeshNode *>() ) {}
+    void Clear() { myLinks.clear(); myIsBndConnected = false; myHasPending = false; }
+    bool SetConnected() { bool was = myIsBndConnected; myIsBndConnected = true; return !was; }
+    bool Contains( const SMDS_MeshNode* n ) const
+    {
+      for ( size_t i = 0; i < myLinks.size(); ++i )
+        if ( myLinks[i]->myNode1 == n ) return true;
+      return false;
+    }
+    virtual int NbNodes() const { return myLinks.size(); }
+    virtual SMDS_ElemIteratorPtr nodesIterator() const
+    {
+      return setNodes(), SMDS_PolygonalFaceOfNodes::nodesIterator();
+    }
+    virtual SMDS_NodeIteratorPtr nodeIterator() const
+    {
+      return setNodes(), SMDS_PolygonalFaceOfNodes::nodeIterator();
+    }
+    void setNodes() const //!< set nodes to SMDS_PolygonalFaceOfNodes
+    {
+      EdgeLoop* me = const_cast<EdgeLoop*>( this );
+      me->myNodes.resize( NbNodes() );
+      size_t iMin = 0;
+      for ( size_t i = 1; i < myNodes.size(); ++i ) {
+        if ( myLinks[ i ]->myNode1->GetID() < myLinks[ iMin ]->myNode1->GetID() )
+          iMin = i;
+      }
+      for ( size_t i = 0; i < myNodes.size(); ++i )
+        me->myNodes[ i ] = myLinks[ ( iMin + i ) % myNodes.size() ]->myNode1;
+    }
+  };
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Set of EdgeLoop's constructed from a CutFace
+   */
+  struct EdgeLoopSet
+  {
+    std::vector< EdgeLoop >  myLoops;       //!< buffer of EdgeLoop's
+    size_t                   myNbLoops;     //!< number of constructed loops
+
+    const EdgePart*          myEdge0;       //!< & CutFace.myLinks[0]
+    size_t                   myNbUsedEdges; //!< nb of EdgePart's added to myLoops
+    boost::dynamic_bitset<>  myIsUsedEdge;  //!< is i-th EdgePart of CutFace is in any EdgeLoop
+    std::vector< EdgeLoop* > myLoopOfEdge;  //!< EdgeLoop of CutFace.myLinks[i]
+    std::vector< EdgePart* > myCandidates;  //!< EdgePart's starting at the same node
+
+    EdgeLoopSet(): myLoops(100) {}
+
+    void Init( const std::vector< EdgePart >& edges )
+    {
+      size_t nb = edges.size();
+      myEdge0 = & edges[0];
+      myNbLoops = 0;
+      myNbUsedEdges = 0;
+      myIsUsedEdge.reset();
+      myIsUsedEdge.resize( nb, false );
+      myLoopOfEdge.clear();
+      myLoopOfEdge.resize( nb, (EdgeLoop*) 0 );
+    }
+    EdgeLoop& AddNewLoop()
+    {
+      if ( ++myNbLoops >= myLoops.size() )
+        myLoops.resize( myNbLoops + 10 );
+      myLoops[ myNbLoops-1 ].Clear();
+      return myLoops[ myNbLoops-1 ];
+    }
+    bool AllEdgesUsed() const { return myNbUsedEdges == myLoopOfEdge.size(); }
+
+    bool AddEdge( EdgePart& edge )
+    {
+      size_t i = Index( edge );
+      if ( myIsUsedEdge[ i ])
+        return false;
+      myLoops[ myNbLoops-1 ].myLinks.push_back( &edge );
+      myLoopOfEdge[ i ] = & myLoops[ myNbLoops-1 ];
+      myIsUsedEdge[ i ] = true;
+      ++myNbUsedEdges;
+      return true;
+    }
+    void Erase( EdgeLoop* loop )
+    {
+      for ( size_t iE = 0; iE < loop->myLinks.size(); ++iE )
+        myLoopOfEdge[ Index( *loop->myLinks[ iE ] )] = 0;
+      loop->Clear();
+    }
+    size_t    Index( const EdgePart& edge ) const { return &edge - myEdge0; }
+    EdgeLoop* GetLoopOf( const EdgePart* edge ) { return myLoopOfEdge[ Index( *edge )]; }
+  };
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Intersections of a face
+   */
+  struct CutFace
+  {
+    mutable std::vector< EdgePart > myLinks;
+    const SMDS_MeshElement*         myInitFace;
+
+    CutFace( const SMDS_MeshElement* face ): myInitFace( face ) {}
+    void AddEdge( const CutLink&          p1,
+                  const CutLink&          p2,
+                  const SMDS_MeshElement* cutter,
+                  const int               nbOnPlane,
+                  const int               iNotOnPlane = -1) const;
+    void AddPoint( const CutLink& p1, const CutLink& p2, double tol ) const;
+    bool ReplaceNodes( const TNNMap& theRm2KeepMap ) const;
+    bool IsCut() const;
+    int  NbInternalEdges() const;
+    void MakeLoops( EdgeLoopSet& loops, const gp_XYZ& theFaceNorm ) const;
+    bool RemoveInternalLoops( EdgeLoopSet& theLoops ) const;
+    void CutOffLoops( EdgeLoopSet&                 theLoops,
+                      const double                 theSign,
+                      const std::vector< gp_XYZ >& theNormals,
+                      std::vector< EdgePart >&     theCutOffLinks,
+                      TLinkMap&                    theCutOffCoplanarLinks) const;
+    void InitLinks() const;
+    bool IsCoplanar( const EdgePart* edge ) const;
+
+    static Standard_Integer HashCode(const CutFace& f, const Standard_Integer upper)
+    {
+      return ::HashCode( f.myInitFace->GetID(), upper );
+    }
+    static Standard_Boolean IsEqual(const CutFace& f1, const CutFace& f2 )
+    {
+      return f1.myInitFace == f2.myInitFace;
+    }
+    void Dump() const;
+
+  private:
+
+    EdgePart* getTwin( const EdgePart* edge ) const;
+  };
+
+  typedef NCollection_Map< CutFace, CutFace > TCutFaceMap;
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Intersection point of two edges of co-planar triangles
+   */
+  struct IntPoint2D
+  {
+    size_t        myEdgeInd[2]; //!< edge indices of triangles
+    double        myU      [2]; //!< parameter [0,1] on edges of triangles
+    SMESH_NodeXYZ myNode;       //!< intersection node
+    bool          myIsCollinear;//!< edges are collinear
+
+    IntPoint2D() : myIsCollinear( false ) {}
+
+    void InitLink( CutLink& link, int iFace, const std::vector< SMESH_NodeXYZ >& nodes ) const
+    {
+      link.Set( nodes[  myEdgeInd[ iFace ]                      ].Node(),
+                nodes[( myEdgeInd[ iFace ] + 1 ) % nodes.size() ].Node(),
+                link.myFace );
+      link.myIntNode = myNode;
+    }
+    const SMDS_MeshNode* Node() const { return myNode.Node(); }
+  };
+  struct IntPoint2DCompare
+  {
+    int myI;
+    IntPoint2DCompare( int iFace=0 ): myI( iFace ) {}
+    bool operator() ( const IntPoint2D* ip1, const IntPoint2D* ip2 ) const
+    {
+      return ip1->myU[ myI ] < ip2->myU[ myI ];
+    }
+    bool operator() ( const IntPoint2D& ip1, const IntPoint2D& ip2 ) const
+    {
+      return ip1.myU[ myI ] < ip2.myU[ myI ];
+    }
+  };
+  typedef boost::container::flat_set< IntPoint2D, IntPoint2DCompare >  TIntPointSet;
+  typedef boost::container::flat_set< IntPoint2D*, IntPoint2DCompare > TIntPointPtrSet;
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Face used to find translated position of the node
+   */
+  struct Face
+  {
+    const SMDS_MeshElement* myFace;
+    SMESH_TNodeXYZ          myNode1; //!< nodes neighboring another node of myFace
+    SMESH_TNodeXYZ          myNode2;
+    const gp_XYZ*           myNorm;
+    bool                    myNodeRightOrder;
+    void operator=(const SMDS_MeshElement* f) { myFace = f; }
+    const SMDS_MeshElement* operator->() { return myFace; }
+    void SetNodes( int i0, int i1 ) //!< set myNode's
+    {
+      myNode1.Set( myFace->GetNode( i1 ));
+      int i2 = ( i0 - 1 + myFace->NbCornerNodes() ) % myFace->NbCornerNodes();
+      if ( i2 == i1 )
+        i2 = ( i0 + 1 ) % myFace->NbCornerNodes();
+      myNode2.Set( myFace->GetNode( i2 ));
+      myNodeRightOrder = ( Abs( i2-i1 ) == 1 ) ?  i2 > i1  :  i2 < i1;
+    }
+    void SetOldNodes( const SMDS_Mesh& theSrcMesh )
+    {
+      myNode1.Set( theSrcMesh.FindNode( myNode1->GetID() ));
+      myNode2.Set( theSrcMesh.FindNode( myNode2->GetID() ));
+    }
+    bool SetNormal( const std::vector< gp_XYZ >& faceNormals )
+    {
+      myNorm = & faceNormals[ myFace->GetID() ];
+      return ( myNorm->SquareModulus() > gp::Resolution() * gp::Resolution() );
+    }
+    const gp_XYZ& Norm() const { return *myNorm; }
+  };
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Offset plane used to find translated position of the node
+   */
+  struct OffsetPlane
+  {
+    gp_XYZ myNode;
+    Face*  myFace;
+    gp_Pln myPln;
+    gp_Lin myLines[2]; //!< line of intersection with neighbor OffsetPlane's
+    bool   myIsLineOk[2];
+    double myWeight[2];
+
+    void   Init( const gp_XYZ& node, Face& tria, double offset )
+    {
+      myNode = node;
+      myFace = & tria;
+      myPln  = gp_Pln( node + tria.Norm() * offset, tria.Norm() );
+      myIsLineOk[0] = myIsLineOk[1] = false;
+      myWeight[0] = myWeight[1] = 0;
+    }
+    bool   ComputeIntersectionLine( OffsetPlane& pln );
+    void   SetSkewLine( const gp_Lin& line );
+    gp_XYZ GetCommonPoint( int & nbOkPoints, double& sumWeight );
+    gp_XYZ ProjectNodeOnLine( int & nbOkPoints );
+    double Weight() const { return myWeight[0] + myWeight[1]; }
+  };
+
+  //================================================================================
+  /*!
+   * \brief Set the second line
+   */
+  //================================================================================
+
+  void OffsetPlane::SetSkewLine( const gp_Lin& line )
+  {
+    myLines[1] = line;
+    gp_XYZ n = myLines[0].Direction().XYZ() ^ myLines[1].Direction().XYZ();
+    if (( myIsLineOk[1] = n.SquareModulus() > gp::Resolution() ))
+      myPln = gp_Pln( myPln.Location(), n );
+  }
+
+  //================================================================================
+  /*!
+   * \brief Project myNode on myLine[0]
+   */
+  //================================================================================
+
+  gp_XYZ OffsetPlane::ProjectNodeOnLine( int & nbOkPoints )
+  {
+    gp_XYZ p = gp::Origin().XYZ();
+    if ( myIsLineOk[0] )
+    {
+      gp_Vec l2n( myLines[0].Location(), myNode );
+      double u = l2n * myLines[0].Direction();
+      p = myLines[0].Location().XYZ() + u * myLines[0].Direction().XYZ();
+      ++nbOkPoints;
+    }
+    return p;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Computes intersection point of myLines
+   */
+  //================================================================================
+
+  gp_XYZ OffsetPlane::GetCommonPoint( int & nbOkPoints, double& sumWeight )
+  {
+    if ( !myIsLineOk[0] || !myIsLineOk[1] )
+    {
+      // sumWeight += myWeight[0];
+      // return ProjectNodeOnLine( nbOkPoints ) * myWeight[0];
+      return gp::Origin().XYZ();
+    }
+
+    gp_XYZ p;
+
+    gp_Vec lPerp0 = myLines[0].Direction().XYZ() ^ myPln.Axis().Direction().XYZ();
+    double  dot01 = lPerp0 * myLines[1].Direction().XYZ();
+    if ( Abs( dot01 ) > 0.05 )
+    {
+      gp_Vec l0l1 = myLines[1].Location().XYZ() - myLines[0].Location().XYZ();
+      double   u1 = - ( lPerp0 * l0l1 ) / dot01;
+      p = ( myLines[1].Location().XYZ() + myLines[1].Direction().XYZ() * u1 );
+    }
+    else
+    {
+      gp_Vec  lv0( myLines[0].Location(), myNode),  lv1(myLines[1].Location(), myNode );
+      double dot0( lv0 * myLines[0].Direction() ), dot1( lv1 * myLines[1].Direction() );
+      p  = 0.5 * ( myLines[0].Location().XYZ() + myLines[0].Direction().XYZ() * dot0 );
+      p += 0.5 * ( myLines[1].Location().XYZ() + myLines[1].Direction().XYZ() * dot1 );
+    }
+
+    sumWeight += Weight();
+    ++nbOkPoints;
+
+    return p * Weight();
+  }
+
+  //================================================================================
+  /*!
+   * \brief Compute line of intersection of 2 planes
+   */
+  //================================================================================
+
+  bool OffsetPlane::ComputeIntersectionLine( OffsetPlane& theNextPln )
+  {
+    const gp_XYZ& n1 = myFace->Norm();
+    const gp_XYZ& n2 = theNextPln.myFace->Norm();
+
+    gp_XYZ lineDir = n1 ^ n2;
+    gp_Pnt linePos;
+
+    double x = Abs( lineDir.X() );
+    double y = Abs( lineDir.Y() );
+    double z = Abs( lineDir.Z() );
+
+    int cooMax; // max coordinate
+    if (x > y) {
+      if (x > z) cooMax = 1;
+      else       cooMax = 3;
+    }
+    else {
+      if (y > z) cooMax = 2;
+      else       cooMax = 3;
+    }
+
+    bool ok = true;
+    if ( Abs( lineDir.Coord( cooMax )) < 0.05 )
+    {
+      // parallel planes - intersection is an offset of the common edge
+      linePos  = 0.5 * ( myPln.Location().XYZ() + theNextPln.myPln.Location().XYZ() );
+      lineDir  = myNode - myFace->myNode2;
+      ok       = false;
+      myWeight[0] = 0;
+    }
+    else
+    {
+      // the constants in the 2 plane equations
+      double d1 = - ( n1 * myPln.Location().XYZ() );
+      double d2 = - ( n2 * theNextPln.myPln.Location().XYZ() );
+
+      switch ( cooMax ) {
+      case 1:
+        linePos.SetX(  0 );
+        linePos.SetY(( d2*n1.Z() - d1*n2.Z()) / lineDir.X() );
+        linePos.SetZ(( d1*n2.Y() - d2*n1.Y()) / lineDir.X() );
+        break;
+      case 2:
+        linePos.SetX(( d1*n2.Z() - d2*n1.Z()) / lineDir.Y() );
+        linePos.SetY(  0 );
+        linePos.SetZ(( d2*n1.X() - d1*n2.X()) / lineDir.Y() );
+        break;
+      case 3:
+        linePos.SetX(( d2*n1.Y() - d1*n2.Y()) / lineDir.Z() );
+        linePos.SetY(( d1*n2.X() - d2*n1.X()) / lineDir.Z() );
+        linePos.SetZ(  0 );
+      }
+      myWeight[0] = lineDir.SquareModulus();
+      if ( n1 * n2 < 0 )
+        myWeight[0] = 2. - myWeight[0];
+    }
+    myLines   [ 0 ].SetDirection( lineDir );
+    myLines   [ 0 ].SetLocation ( linePos );
+    myIsLineOk[ 0 ] = ok;
+
+    theNextPln.myLines   [ 1 ] = myLines[ 0 ];
+    theNextPln.myIsLineOk[ 1 ] = ok;
+    theNextPln.myWeight  [ 1 ] = myWeight[ 0 ];
+
+    return ok;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Return a translated position of a node
+   *  \param [in] new2OldNodes - new and old nodes
+   *  \param [in] faceNormals - normals to input faces
+   *  \param [in] theSrcMesh - initial mesh
+   *  \param [in] theNewPos - a computed normal
+   *  \return bool - true if theNewPos is computed
+   */
+  //================================================================================
+
+  bool getTranslatedPosition( const SMDS_MeshNode*         theNewNode,
+                              const double                 theOffset,
+                              const double                 theTol,
+                              const double                 theSign,
+                              const std::vector< gp_XYZ >& theFaceNormals,
+                              SMDS_Mesh&                   theSrcMesh,
+                              gp_XYZ&                      theNewPos)
+  {
+    bool useOneNormal = true;
+
+    // check if theNewNode needs an average position, i.e. theNewNode is convex
+    // SMDS_ElemIteratorPtr faceIt = theNewNode->GetInverseElementIterator();
+    // const SMDS_MeshElement*  f0 = faceIt->next();
+    // const gp_XYZ&         norm0 = theFaceNormals[ f0->GetID() ];
+    // const SMESH_NodeXYZ nodePos = theNewNode;
+    // while ( faceIt->more() )
+    // {
+    //   const SMDS_MeshElement* f = faceIt->next();
+    //   const int         nodeInd = f->GetNodeIndex( theNewNode );
+    //   SMESH_NodeXYZ    nodePos2 = f->GetWrapNode( nodeInd + 1 );
+    //   try {
+    //     const gp_XYZ      nnDir = ( nodePos2 - nodePos ).Normalized();
+    //   }
+    //   catch {
+    //     continue;
+    //   }
+    //   const double dot = norm0 * nnDir;
+    //   bool isConvex = 
+
+
+
+    // get faces surrounding theNewNode and sort them
+    Face faces[ theMaxNbFaces ];
+    SMDS_ElemIteratorPtr faceIt = theNewNode->GetInverseElementIterator();
+    faces[0] = faceIt->next();
+    while ( !faces[0].SetNormal( theFaceNormals ) && faceIt->more() )
+      faces[0] = faceIt->next();
+    int i0 = faces[0]->GetNodeIndex( theNewNode );
+    int i1 = ( i0 + 1 ) % faces[0]->NbCornerNodes();
+    faces[0].SetNodes( i0, i1 );
+    TIDSortedElemSet elemSet, avoidSet;
+    int iFace = 0;
+    const SMDS_MeshElement* f;
+    for ( ; faceIt->more(); faceIt->next() )
+    {
+      avoidSet.insert( faces[ iFace ].myFace );
+      f = SMESH_MeshAlgos::FindFaceInSet( theNewNode, faces[ iFace ].myNode2.Node(),
+                                          elemSet, avoidSet, &i0, &i1 );
+      if ( !f )
+      {
+        std::reverse( &faces[0], &faces[0] + iFace + 1 );
+        for ( int i = 0; i <= iFace; ++i )
+        {
+          std::swap( faces[i].myNode1, faces[i].myNode2 );
+          faces[i].myNodeRightOrder = !faces[i].myNodeRightOrder;
+        }
+        f = SMESH_MeshAlgos::FindFaceInSet( theNewNode, faces[ iFace ].myNode2.Node(),
+                                            elemSet, avoidSet, &i0, &i1 );
+        if ( !f )
+          break;
+      }
+      faces[ ++iFace ] = f;
+      faces[ iFace ].SetNodes( i0, i1 );
+      faces[ iFace ].SetNormal( theFaceNormals );
+    }
+    int nbFaces = Min( iFace + 1, (int)theMaxNbFaces );
+
+    theNewPos.SetCoord( 0, 0, 0 );
+    gp_XYZ oldXYZ = SMESH_NodeXYZ( theNewNode );
+
+    // check if all faces are co-planar
+    bool isPlanar = true;
+    const double tol = 1e-2;
+    for ( int i = 1; i < nbFaces &&  isPlanar;  ++i )
+      isPlanar = ( faces[i].Norm() - faces[i-1].Norm() ).SquareModulus() < tol*tol;
+
+    if ( isPlanar )
+    {
+      theNewPos = oldXYZ + faces[0].Norm() * theOffset;
+      return useOneNormal;
+    }
+
+    // prepare OffsetPlane's
+    OffsetPlane pln[ theMaxNbFaces ];
+    for ( int i = 0; i < nbFaces; ++i )
+    {
+      faces[i].SetOldNodes( theSrcMesh );
+      pln[i].Init( oldXYZ, faces[i], theOffset );
+    }
+    // intersect neighboring OffsetPlane's
+    int nbOkPoints = 0;
+    for ( int i = 1; i < nbFaces; ++i )
+      nbOkPoints += pln[ i-1 ].ComputeIntersectionLine( pln[ i ]);
+    nbOkPoints += pln[ nbFaces-1 ].ComputeIntersectionLine( pln[ 0 ]);
+
+    // move intersection lines to over parallel planes
+    if ( nbOkPoints > 1 )
+      for ( int i = 0; i < nbFaces; ++i )
+        if ( pln[i].myIsLineOk[0] && !pln[i].myIsLineOk[1] )
+          for ( int j = 1; j < nbFaces &&  !pln[i].myIsLineOk[1]; ++j )
+          {
+            int i2 = ( i + j ) % nbFaces;
+            if ( pln[i2].myIsLineOk[0] )
+              pln[i].SetSkewLine( pln[i2].myLines[0] );
+          }
+
+    // get the translated position
+    nbOkPoints = 0;
+    double sumWegith = 0;
+    const double minWeight = Sin( 30 * M_PI / 180. ) * Sin( 30 * M_PI / 180. );
+    for ( int i = 0; i < nbFaces; ++i )
+      if ( pln[ i ].Weight() > minWeight )
+        theNewPos += pln[ i ].GetCommonPoint( nbOkPoints, sumWegith );
+
+    if ( nbOkPoints == 0 )
+    {
+      // there is only one feature edge;
+      // find the theNewPos by projecting oldXYZ to any intersection line
+      for ( int i = 0; i < nbFaces; ++i )
+        theNewPos += pln[ i ].ProjectNodeOnLine( nbOkPoints );
+
+      if ( nbOkPoints == 0 )
+      {
+        theNewPos = oldXYZ + faces[0].Norm() * theOffset;
+        return useOneNormal;
+      }
+      sumWegith = nbOkPoints;
+    }
+    theNewPos /= sumWegith;
+
+
+    // mark theNewNode if it is concave
+    useOneNormal = false;
+    gp_Vec moveVec( oldXYZ, theNewPos );
+    for ( int i = 0, iPrev = nbFaces-1; i < nbFaces; iPrev = i++ )
+    {
+      gp_Vec nodeVec( oldXYZ, faces[ i ].myNode1 );
+      double u = ( moveVec * nodeVec ) / nodeVec.SquareMagnitude();
+      if ( u > 0.5 ) // param [0,1] on nodeVec
+      {
+        theNewNode->setIsMarked( true );
+      }
+      if ( !useOneNormal )
+      {
+        gp_XYZ inFaceVec = faces[ i ].Norm() ^ nodeVec.XYZ();
+        double       dot = inFaceVec * faces[ iPrev ].Norm();
+        if ( !faces[ i ].myNodeRightOrder )
+          dot *= -1;
+        if ( dot * theSign < 0 )
+        {
+          gp_XYZ p1 = oldXYZ + faces[ i ].Norm()     * theOffset;
+          gp_XYZ p2 = oldXYZ + faces[ iPrev ].Norm() * theOffset;
+          useOneNormal = ( p1 - p2 ).SquareModulus() > theTol * theTol;
+        }
+      }
+      if ( useOneNormal && theNewNode->isMarked() )
+        break;
+    }
+
+    return useOneNormal;
+  }
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Intersect faces of a mesh
+   */
+  struct Intersector
+  {
+    SMDS_Mesh*                   myMesh;
+    double                       myTol, myEps;
+    const std::vector< gp_XYZ >& myNormals;
+    TCutLinkMap                  myCutLinks; //!< assure sharing of new nodes
+    TCutFaceMap                  myCutFaces;
+    TNNMap                       myRemove2KeepNodes; //!< node merge map
+
+    // data to intersect 2 faces
+    const SMDS_MeshElement*      myFace1;
+    const SMDS_MeshElement*      myFace2;
+    std::vector< SMESH_NodeXYZ > myNodes1, myNodes2;
+    std::vector< double >        myDist1,  myDist2;
+    int                          myInd1, myInd2; // coordinate indices on an axis-aligned plane
+    int                          myNbOnPlane1, myNbOnPlane2;
+    TIntPointSet                 myIntPointSet;
+
+    Intersector( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals )
+      : myMesh( mesh ),
+        myTol( tol ),
+        myEps( 1e-100 ),
+        //myEps( Sqrt( std::numeric_limits<double>::min() )),
+        //myEps( gp::Resolution() ),
+        myNormals( normals )
+    {}
+    void Cut( const SMDS_MeshElement* face1,
+              const SMDS_MeshElement* face2,
+              const int               nbCommonNodes );
+    void MakeNewFaces( SMESH_MeshAlgos::TEPairVec& theNew2OldFaces,
+                       SMESH_MeshAlgos::TNPairVec& theNew2OldNodes,
+                       const double                theSign );
+
+  private:
+
+    bool isPlaneIntersected( const gp_XYZ&                       n2,
+                             const double                        d2,
+                             const std::vector< SMESH_NodeXYZ >& nodes1,
+                             std::vector< double > &             dist1,
+                             int &                               nbOnPlane1,
+                             int &                               iNotOnPlane1);
+    void computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes,
+                           const std::vector< double >&        dist,
+                           const int                           nbOnPln,
+                           const int                           iMaxCoo,
+                           double *                            u,
+                           int*                                iE);
+    void cutCoplanar();
+    void addLink ( CutLink& link );
+    bool findLink( CutLink& link );
+    bool coincide( const gp_XYZ& p1, const gp_XYZ& p2, const double tol ) const
+    {
+      return ( p1 - p2 ).SquareModulus() < tol * tol;
+    }
+    gp_XY p2D( const gp_XYZ& p ) const { return gp_XY( p.Coord( myInd1 ), p.Coord( myInd2 )); }
+
+    void intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1,
+                        const std::vector< double > &       dist1,
+                        const int                           iEdge1,
+                        const SMDS_MeshElement*             face2,
+                        CutLink&                            link1);
+    void findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes,
+                              const std::vector< double > &       dist,
+                              CutLink&                            link );
+    void replaceIntNode( const SMDS_MeshNode* nToKeep, const SMDS_MeshNode* nToRemove );
+    void computeIntPoint( const double           u1,
+                          const double           u2,
+                          const int              iE1,
+                          const int              iE2,
+                          CutLink &              link,
+                          const SMDS_MeshNode* & node1,
+                          const SMDS_MeshNode* & node2);
+    void cutCollinearLink( const int                           iNotOnPlane1,
+                           const std::vector< SMESH_NodeXYZ >& nodes1,
+                           const SMDS_MeshElement*             face2,
+                           const CutLink&                      link1,
+                           const CutLink&                      link2);
+    void setPlaneIndices( const gp_XYZ& planeNorm );
+    bool intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1,
+                            const gp_XY s2p0, const gp_XY s2p1,
+                            double &    t1,   double &    t2,
+                            bool &      isCollinear  );
+    bool intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint );
+    bool isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes );
+    void intersectNewEdges( const CutFace& theCFace );
+    const SMDS_MeshNode* createNode( const gp_XYZ& p );
+  };
+
+  //================================================================================
+  /*!
+   * \brief Return coordinate index with maximal abs value
+   */
+  //================================================================================
+
+  int MaxIndex( const gp_XYZ& x )
+  {
+    int iMaxCoo = ( Abs( x.X()) < Abs( x.Y() )) + 1;
+    if ( Abs( x.Coord( iMaxCoo )) < Abs( x.Z() ))
+      iMaxCoo = 3;
+    return iMaxCoo;
+  }
+  //================================================================================
+  /*!
+   * \brief Store a CutLink
+   */
+  //================================================================================
+
+  const SMDS_MeshNode* Intersector::createNode( const gp_XYZ& p )
+  {
+    const SMDS_MeshNode* n = myMesh->AddNode( p.X(), p.Y(), p.Z() );
+    n->setIsMarked( true ); // cut nodes are marked
+    return n;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Store a CutLink
+   */
+  //================================================================================
+
+  void Intersector::addLink( CutLink& link )
+  {
+    link.myIndex = 0;
+    const CutLink* added = & myCutLinks.Added( link );
+    while ( added->myIntNode.Node() != link.myIntNode.Node() )
+    {
+      if ( !added->myIntNode )
+      {
+        added->myIntNode = link.myIntNode;
+        break;
+      }
+      else
+      {
+        link.myIndex++;
+        added = & myCutLinks.Added( link );
+      }
+    }
+    link.myIndex = 0;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Find a CutLink with an intersection point coincident with that of a given link
+   */
+  //================================================================================
+
+  bool Intersector::findLink( CutLink& link )
+  {
+    link.myIndex = 0;
+    while ( myCutLinks.Contains( link ))
+    {
+      const CutLink* added = & myCutLinks.Added( link );
+      if ( !!added->myIntNode && coincide( added->myIntNode, link.myIntNode, myTol ))
+      {
+        link.myIntNode = added->myIntNode;
+        return true;
+      }
+      link.myIndex++;
+    }
+    return false;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Check if a triangle intersects the plane of another triangle
+   *  \param [in] nodes1 - nodes of triangle 1
+   *  \param [in] n2 - normal of triangle 2
+   *  \param [in] d2 - a constant of the plane equation 2
+   *  \param [out] dist1 - distance of nodes1 from the plane 2
+   *  \param [out] nbOnPlane - number of nodes1 lying on the plane 2
+   *  \return bool - true if the triangle intersects the plane 2
+   */
+  //================================================================================
+
+  bool Intersector::isPlaneIntersected( const gp_XYZ&                       n2,
+                                        const double                        d2,
+                                        const std::vector< SMESH_NodeXYZ >& nodes1,
+                                        std::vector< double > &             dist1,
+                                        int &                               nbOnPlane1,
+                                        int &                               iNotOnPlane1)
+  {
+    iNotOnPlane1 = nbOnPlane1 = 0;
+    dist1.resize( nodes1.size() );
+    for ( size_t i = 0; i < nodes1.size(); ++i )
+    {
+      dist1[i] = n2 * nodes1[i] + d2;
+      if ( Abs( dist1[i] ) < myTol )
+      {
+        ++nbOnPlane1;
+        dist1[i] = 0.;
+      }
+      else
+      {
+        iNotOnPlane1 = i;
+      }
+    }
+    if ( nbOnPlane1 == 0 )
+      for ( size_t i = 0; i < nodes1.size(); ++i )
+        if ( dist1[iNotOnPlane1] * dist1[i] < 0 )
+          return true;
+
+    return nbOnPlane1;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Compute parameters on the plane intersection line of intersections
+   *        of edges of a triangle
+   *  \param [in] nodes - triangle nodes
+   *  \param [in] dist - distance of triangle nodes from the plane of another triangle
+   *  \param [in] nbOnPln - number of nodes lying on the plane of another triangle
+   *  \param [in] iMaxCoo - index of coordinate of max component of the plane intersection line
+   *  \param [out] u - two computed parameters on the plane intersection line
+   *  \param [out] iE - indices of intersected edges
+   */
+  //================================================================================
+
+  void Intersector::computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes,
+                                      const std::vector< double >&        dist,
+                                      const int                           nbOnPln, 
+                                      const int                           iMaxCoo,
+                                      double *                            u,
+                                      int*                                iE)
+  {
+    if ( nbOnPln == 3 )
+    {
+      u[0] = u[1] = 1e+100;
+      return;
+    }
+    int nb = 0;
+    int i1 = 2, i2 = 0;
+    if ( nbOnPln == 1 && ( dist[i1] == 0. || dist[i2] == 0 ))
+    {
+      int i = dist[i1] == 0 ? i1 : i2;
+      u [ 1 ] = nodes[ i ].Coord( iMaxCoo );
+      iE[ 1 ] = i;
+      i1 = i2++;
+    }
+    for ( ; i2 < 3 && nb < 2;  i1 = i2++ )
+    {
+      double dd = dist[i1] - dist[i2];
+      if ( dd != 0. && dist[i2] * dist[i1] <= 0. )
+      {
+        double x1 = nodes[i1].Coord( iMaxCoo );
+        double x2 = nodes[i2].Coord( iMaxCoo );
+        u [ nb ] = x1 + ( x2 - x1 ) * dist[i1] / dd;
+        iE[ nb ] = i1;
+        ++nb;
+      }
+    }
+    if ( u[0] > u[1] )
+    {
+      std::swap( u [0], u [1] );
+      std::swap( iE[0], iE[1] );
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Try to find an intersection node on a link collinear with the plane intersection line
+   */
+  //================================================================================
+
+  void Intersector::findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes,
+                                         const std::vector< double > &       dist,
+                                         CutLink&                            link )
+  {
+    int i1 = ( dist[0] == 0 ? 0 : 1 ), i2 = ( dist[2] == 0 ? 2 : 1 );
+    CutLink link2 = link;
+    link2.Set( nodes[i1].Node(), nodes[i2].Node(), 0 );
+    if ( findLink( link2 ))
+      link.myIntNode = link2.myIntNode;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Compute intersection point of a link1 with a face2
+   */
+  //================================================================================
+
+  void Intersector::intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1,
+                                   const std::vector< double > &       dist1,
+                                   const int                           iEdge1,
+                                   const SMDS_MeshElement*             face2,
+                                   CutLink&                            link1)
+  {
+    const int iEdge2 = ( iEdge1 + 1 ) % nodes1.size();
+    const SMESH_NodeXYZ& p1 = nodes1[ iEdge1 ];
+    const SMESH_NodeXYZ& p2 = nodes1[ iEdge2 ];
+
+    link1.Set( p1.Node(), p2.Node(), face2 );
+    const CutLink* link = & myCutLinks.Added( link1 );
+    if ( !link->IntNode() )
+    {
+      if      ( dist1[ iEdge1 ] == 0. ) link1.myIntNode = p1;
+      else if ( dist1[ iEdge2 ] == 0. ) link1.myIntNode = p2;
+      else
+      {
+        gp_XYZ p = p1 + ( p2 - p1 ) * dist1[ iEdge1 ] / ( dist1[ iEdge1 ] - dist1[ iEdge2 ]);
+        (gp_XYZ&)link1.myIntNode = p;
+      }
+    }
+    else
+    {
+      gp_XYZ p = p1 + ( p2 - p1 ) * dist1[ iEdge1 ] / ( dist1[ iEdge1 ] - dist1[ iEdge2 ]);
+      while ( link->IntNode() )
+      {
+        if ( coincide( p, link->myIntNode, myTol ))
+        {
+          link1.myIntNode = link->myIntNode;
+          break;
+        }
+        link1.myIndex++;
+        link = & myCutLinks.Added( link1 );
+      }
+      if ( !link1.IntNode() )
+      {
+        if      ( dist1[ iEdge1 ] == 0. ) link1.myIntNode = p1;
+        else if ( dist1[ iEdge2 ] == 0. ) link1.myIntNode = p2;
+        else                     (gp_XYZ&)link1.myIntNode = p;
+      }
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Store node replacement in myCutFaces
+   */
+  //================================================================================
+
+  void Intersector::replaceIntNode( const SMDS_MeshNode* nToKeep,
+                                    const SMDS_MeshNode* nToRemove )
+  {
+    if ( nToKeep == nToRemove )
+      return;
+    if ( nToRemove->GetID() < nToKeep->GetID() ) // keep node with lower ID
+      myRemove2KeepNodes.Bind((void*) nToKeep, nToRemove );
+    else
+      myRemove2KeepNodes.Bind((void*) nToRemove, nToKeep );
+  }
+
+  //================================================================================
+  /*!
+   * \brief Compute intersection point on a link of either of faces by choosing
+   *        a link whose parameter on the intersection line in maximal
+   *  \param [in] u1 - parameter on the intersection line of link iE1 of myFace1
+   *  \param [in] u2 - parameter on the intersection line of link iE2 of myFace2
+   *  \param [in] iE1 - index of a link myFace1
+   *  \param [in] iE2 - index of a link myFace2
+   *  \param [out] link - CutLink storing the intersection point
+   *  \param [out] node1 - a node of the 2nd link if two links intersect
+   *  \param [out] node2 - a node of the 2nd link if two links intersect
+   */
+  //================================================================================
+
+  void Intersector::computeIntPoint( const double           u1,
+                                     const double           u2,
+                                     const int              iE1,
+                                     const int              iE2,
+                                     CutLink &              link,
+                                     const SMDS_MeshNode* & node1,
+                                     const SMDS_MeshNode* & node2)
+  {
+    if      ( u1 > u2 + myTol )
+    {
+      intersectLink( myNodes1, myDist1, iE1, myFace2, link );
+      node1 = node2 = 0;
+      if ( myNbOnPlane2 == 2 )
+        findIntPointOnPlane( myNodes2, myDist2, link );
+    }
+    else if ( u2 > u1 + myTol )
+    {
+      intersectLink( myNodes2, myDist2, iE2, myFace1, link );
+      node1 = node2 = 0;
+      if ( myNbOnPlane1 == 2 )
+        findIntPointOnPlane( myNodes1, myDist1, link );
+    }
+    else // edges of two faces intersect the line at the same point
+    {
+      CutLink link2;
+      intersectLink( myNodes1, myDist1, iE1, myFace2, link );
+      intersectLink( myNodes2, myDist2, iE2, myFace1, link2 );
+      node1 = link2.Node1();
+      node2 = link2.Node2();
+
+      if      ( !link.IntNode() && link2.IntNode() )
+        link.myIntNode = link2.myIntNode;
+
+      else if ( !link.IntNode() && !link2.IntNode() )
+        (gp_XYZ&)link.myIntNode = 0.5 * ( link.myIntNode + link2.myIntNode );
+
+      else if ( link.IntNode() && link2.IntNode() )
+        replaceIntNode( link.IntNode(), link2.IntNode() );
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Add intersections to a link collinear with the intersection line
+   */
+  //================================================================================
+
+  void Intersector::cutCollinearLink( const int                           iNotOnPlane1,
+                                      const std::vector< SMESH_NodeXYZ >& nodes1,
+                                      const SMDS_MeshElement*             face2,
+                                      const CutLink&                      link1,
+                                      const CutLink&                      link2)
+
+  {
+    int iN1 = ( iNotOnPlane1 + 1 ) % 3;
+    int iN2 = ( iNotOnPlane1 + 2 ) % 3;
+    CutLink link( nodes1[ iN1 ].Node(), nodes1[ iN2 ].Node(), face2 );
+    if ( link1.myFace != face2 )
+    {
+      link.myIntNode = link1.myIntNode;
+      addLink( link );
+    }
+    if ( link2.myFace != face2 )
+    {
+      link.myIntNode = link2.myIntNode;
+      addLink( link );
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Choose indices on an axis-aligned plane
+   */
+  //================================================================================
+
+  void Intersector::setPlaneIndices( const gp_XYZ& planeNorm )
+  {
+    switch ( MaxIndex( planeNorm )) {
+    case 1: myInd1 = 2; myInd2 = 3; break;
+    case 2: myInd1 = 3; myInd2 = 1; break;
+    case 3: myInd1 = 1; myInd2 = 2; break;
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Intersect two faces
+   */
+  //================================================================================
+
+  void Intersector::Cut( const SMDS_MeshElement* face1,
+                         const SMDS_MeshElement* face2,
+                         const int               nbCommonNodes)
+  {
+    myFace1 = face1;
+    myFace2 = face2;
+    myNodes1.assign( face1->begin_nodes(), face1->end_nodes() );
+    myNodes2.assign( face2->begin_nodes(), face2->end_nodes() );
+
+    const gp_XYZ& n1 = myNormals[ face1->GetID() ];
+    const gp_XYZ& n2 = myNormals[ face2->GetID() ];
+
+    // check if triangles intersect
+    int iNotOnPlane1, iNotOnPlane2;
+    const double d2 = -( n2 * myNodes2[0]);
+    if ( !isPlaneIntersected( n2, d2, myNodes1, myDist1, myNbOnPlane1, iNotOnPlane1 ))
+      return;
+    const double d1 = -( n1 * myNodes1[0]);
+    if ( !isPlaneIntersected( n1, d1, myNodes2, myDist2, myNbOnPlane2, iNotOnPlane2 ))
+      return;
+
+    if ( myNbOnPlane1 == 3 || myNbOnPlane2 == 3 )// triangles are co-planar
+    {
+      setPlaneIndices( myNbOnPlane1 == 3 ? n2 : n1 ); // choose indices on an axis-aligned plane
+      cutCoplanar();
+    }
+    else if ( nbCommonNodes < 2 ) // triangle planes intersect
+    {
+      gp_XYZ lineDir = n1 ^ n2; // intersection line
+
+      // check if intervals of intersections of triangles with lineDir overlap
+
+      double u1[2], u2 [2]; // parameters on lineDir of edge intersection points { minU, maxU }
+      int   iE1[2], iE2[2]; // indices of edges
+      int iMaxCoo = MaxIndex( lineDir );
+      computeIntervals( myNodes1, myDist1, myNbOnPlane1, iMaxCoo, u1, iE1 );
+      computeIntervals( myNodes2, myDist2, myNbOnPlane2, iMaxCoo, u2, iE2 );
+      if ( u1[1] < u2[0] - myTol || u2[1] < u1[0] - myTol )
+        return; // intervals do not overlap
+
+      // make intersection nodes
+
+      const SMDS_MeshNode *l1n1, *l1n2, *l2n1, *l2n2;
+      CutLink link1; // intersection with smaller u on lineDir
+      computeIntPoint( u1[0], u2[0], iE1[0], iE2[0], link1, l1n1, l1n2 );
+      CutLink link2; // intersection with larger u on lineDir
+      computeIntPoint( -u1[1], -u2[1], iE1[1], iE2[1], link2, l2n1, l2n2 );
+
+      const CutFace& cf1 = myCutFaces.Added( CutFace( face1 ));
+      const CutFace& cf2 = myCutFaces.Added( CutFace( face2 ));
+
+      if ( coincide( link1.myIntNode, link2.myIntNode, myTol ))
+      {
+        // intersection is a point
+        if ( link1.IntNode() && link2.IntNode() )
+          replaceIntNode( link1.IntNode(), link2.IntNode() );
+
+        CutLink* link = link2.IntNode() ? &link2 : &link1;
+        if ( !link->IntNode() )
+        {
+          gp_XYZ p = 0.5 * ( link1.myIntNode + link2.myIntNode );
+          link->myIntNode.Set( createNode( p ));
+        }
+        if ( !link1.IntNode() ) link1.myIntNode = link2.myIntNode;
+        if ( !link2.IntNode() ) link2.myIntNode = link1.myIntNode;
+
+        cf1.AddPoint( link1, link2, myTol );
+        cf2.AddPoint( link1, link2, myTol );
+      }
+      else
+      {
+        // intersection is a line segment
+        if ( !link1.IntNode() )
+          link1.myIntNode.Set( createNode( link1.myIntNode ));
+        if ( !link2.IntNode() )
+          link2.myIntNode.Set( createNode( link2.myIntNode ));
+
+        cf1.AddEdge( link1, link2, face2, myNbOnPlane1, iNotOnPlane1 );
+        if ( l1n1 ) link1.Set( l1n1, l1n2, face2 );
+        if ( l2n1 ) link2.Set( l2n1, l2n2, face2 );
+        cf2.AddEdge( link1, link2, face1, myNbOnPlane2, iNotOnPlane2 );
+
+        // add intersections to a link collinear with the intersection line
+        if ( myNbOnPlane1 == 2 && ( link1.myFace != face2 || link2.myFace != face2 ))
+          cutCollinearLink( iNotOnPlane1, myNodes1, face2, link1, link2 );
+
+        if ( myNbOnPlane2 == 2 && ( link1.myFace != face1 || link2.myFace != face1 ))
+          cutCollinearLink( iNotOnPlane2, myNodes2, face1, link1, link2 );
+      }
+
+      addLink( link1 );
+      addLink( link2 );
+
+    } // non co-planar case
+
+    return;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Intersect two 2D line segments
+   */
+  //================================================================================
+
+  bool Intersector::intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1,
+                                       const gp_XY s2p0, const gp_XY s2p1,
+                                       double &    t1,   double &    t2,
+                                       bool &      isCollinear )
+  {
+    gp_XY u = s1p1 - s1p0;
+    gp_XY v = s2p1 - s2p0;
+    gp_XY w = s1p0 - s2p0;
+    double perpDotUV = u * gp_XY( -v.Y(), v.X() );
+    double perpDotVW = v * gp_XY( -w.Y(), w.X() );
+    double perpDotUW = u * gp_XY( -w.Y(), w.X() );
+    double        u2 = u.SquareModulus();
+    double        v2 = v.SquareModulus();
+    if ( u2 < myEps * myEps || v2 < myEps * myEps )
+      return false;
+    if ( perpDotUV * perpDotUV / u2 / v2 < 1e-6 ) // cos ^ 2
+    {
+      if ( !isCollinear )
+        return false; // no need in collinear solution
+      if ( perpDotUW * perpDotUW / u2 > myTol * myTol )
+        return false; // parallel
+
+      // collinear
+      gp_XY w2 = s1p1 - s2p0;
+      if ( Abs( v.X()) + Abs( u.X()) > Abs( v.Y()) + Abs( u.Y())) {
+        t1 = w.X() / v.X();  // params on segment 2
+        t2 = w2.X() / v.X();
+      }
+      else {
+        t1 = w.Y() / v.Y();
+        t2 = w2.Y() / v.Y();
+      }
+      if ( Max( t1,t2 ) <= 0 || Min( t1,t2 ) >= 1 )
+        return false; // no overlap
+      return true;
+    }
+    isCollinear = false;
+
+    t1 = perpDotVW / perpDotUV; // param on segment 1
+    if ( t1 < 0. || t1 > 1. )
+      return false; // intersection not within the segment
+
+    t2 = perpDotUW / perpDotUV; // param on segment 2
+    if ( t2 < 0. || t2 > 1. )
+      return false; // intersection not within the segment
+
+    return true;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Intersect two edges of co-planar triangles
+   *  \param [inout] iE1 - edge index of triangle 1
+   *  \param [inout] iE2 - edge index of triangle 2
+   *  \param [inout] intPoints - intersection points
+   *  \param [inout] nbIntPoints - nb of found intersection points
+   */
+  //================================================================================
+
+  bool Intersector::intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint )
+  {
+    int i01 = iE1, i11 = ( iE1 + 1 ) % 3;
+    int i02 = iE2, i12 = ( iE2 + 1 ) % 3;
+    if (( !intPoint.myIsCollinear ) &&
+        ( myNodes1[ i01 ] == myNodes2[ i02 ] ||
+          myNodes1[ i01 ] == myNodes2[ i12 ] ||
+          myNodes1[ i11 ] == myNodes2[ i02 ] ||
+          myNodes1[ i11 ] == myNodes2[ i12 ] ))
+      return false;
+
+    // segment 1
+    gp_XY s1p0 = p2D( myNodes1[ i01 ]);
+    gp_XY s1p1 = p2D( myNodes1[ i11 ]);
+
+    // segment 2
+    gp_XY s2p0 = p2D( myNodes2[ i02 ]);
+    gp_XY s2p1 = p2D( myNodes2[ i12 ]);
+
+    double t1, t2;
+    if ( !intersectEdgeEdge( s1p0,s1p1, s2p0,s2p1, t1, t2, intPoint.myIsCollinear ))
+      return false;
+
+    intPoint.myEdgeInd[0] = iE1;
+    intPoint.myEdgeInd[1] = iE2;
+    intPoint.myU[0] = t1;
+    intPoint.myU[1] = t2;
+    (gp_XYZ&)intPoint.myNode = myNodes1[i01] * ( 1 - t1 ) + myNodes1[i11] * t1;
+
+    if ( intPoint.myIsCollinear )
+      return true;
+
+    // try to find existing node at intPoint.myNode
+
+    if ( myNodes1[ i01 ] == myNodes2[ i02 ] ||
+         myNodes1[ i01 ] == myNodes2[ i12 ] ||
+         myNodes1[ i11 ] == myNodes2[ i02 ] ||
+         myNodes1[ i11 ] == myNodes2[ i12 ] )
+      return false;
+
+    const double coincTol = myTol * 1e-3;
+
+    CutLink link1( myNodes1[i01].Node(), myNodes1[i11].Node(), myFace2 );
+    CutLink link2( myNodes2[i02].Node(), myNodes2[i12].Node(), myFace1 );
+
+    SMESH_NodeXYZ& n1 = myNodes1[ t1 < 0.5 ? i01 : i11 ];
+    bool same1 = coincide( n1, intPoint.myNode, coincTol );
+    if ( same1 )
+    {
+      link2.myIntNode = intPoint.myNode = n1;
+      addLink( link2 );
+    }
+    SMESH_NodeXYZ& n2 = myNodes2[ t2 < 0.5 ? i02 : i12 ];
+    bool same2 = coincide( n2, intPoint.myNode, coincTol );
+    if ( same2 )
+    {
+      link1.myIntNode = intPoint.myNode = n2;
+      addLink( link1 );
+      if ( same1 )
+      {
+        replaceIntNode( n1.Node(), n2.Node() );
+        return false;
+      }
+      return true;
+    }
+    if ( same1 )
+      return true;
+
+    link1.myIntNode = intPoint.myNode;
+    if ( findLink( link1 ))
+    {
+      intPoint.myNode = link2.myIntNode = link1.myIntNode;
+      addLink( link2 );
+      return true;
+    }
+
+    link2.myIntNode = intPoint.myNode;
+    if ( findLink( link2 ))
+    {
+      intPoint.myNode = link1.myIntNode = link2.myIntNode;
+      addLink( link1 );
+      return true;
+    }
+
+    for ( int is2nd = 0; is2nd < 2; ++is2nd )
+    {
+      const SMDS_MeshElement* f = is2nd ? myFace1 : myFace2;
+      const CutFace&         cf = myCutFaces.Added( CutFace( is2nd ? myFace2 : myFace1 ));
+      for ( size_t i = 0; i < cf.myLinks.size(); ++i )
+        if ( cf.myLinks[i].myFace == f &&
+             //cf.myLinks[i].myIndex != EdgePart::_COPLANAR &&
+             coincide( intPoint.myNode, SMESH_NodeXYZ( cf.myLinks[i].myNode1 ), coincTol ))
+        {
+          intPoint.myNode.Set( cf.myLinks[i].myNode1 );
+          return true;
+        }
+    }
+
+    // make a new node
+
+    intPoint.myNode._node = createNode( intPoint.myNode );
+    link1.myIntNode = link2.myIntNode = intPoint.myNode;
+    addLink( link1 );
+    addLink( link2 );
+
+    return true;
+  }
+
+
+  //================================================================================
+  /*!
+   * \brief Check if a point is contained in a triangle
+   */
+  //================================================================================
+
+  bool Intersector::isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes )
+  {
+    double bc1, bc2;
+    SMESH_MeshAlgos::GetBarycentricCoords( p2D( p ),
+                                           p2D( nodes[0] ), p2D( nodes[1] ), p2D( nodes[2] ),
+                                           bc1, bc2 );
+    return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. );
+  }
+
+  //================================================================================
+  /*!
+   * \brief Intersect two co-planar faces
+   */
+  //================================================================================
+
+  void Intersector::cutCoplanar()
+  {
+    // find intersections of edges
+
+    IntPoint2D intPoints[ 6 ];
+    int      nbIntPoints = 0;
+    for ( int iE1 = 0; iE1 < 3; ++iE1 )
+    {
+      int maxNbIntPoints = nbIntPoints + 2;
+      for ( int iE2 = 0; iE2 < 3 &&  nbIntPoints < maxNbIntPoints; ++iE2 )
+        nbIntPoints += intersectEdgeEdge( iE1, iE2, intPoints[ nbIntPoints ]);
+    }
+    const int minNbOnPlane = Min( myNbOnPlane1, myNbOnPlane2 );
+
+    if ( nbIntPoints == 0 ) // no intersections of edges
+    {
+      bool is1in2;
+      if      ( isPointInTriangle( myNodes1[0], myNodes2 )) // face2 includes face1
+        is1in2 = true;
+      else if ( isPointInTriangle( myNodes2[0], myNodes1 )) // face1 includes face2
+        is1in2 = false;
+      else
+        return;
+
+      // add edges of an inner triangle to an outer one
+
+      const std::vector< SMESH_NodeXYZ >& nodesIn = is1in2 ? myNodes1 : myNodes2;
+      const SMDS_MeshElement*             faceOut = is1in2 ? myFace2  : myFace1;
+      const SMDS_MeshElement*              faceIn = is1in2 ? myFace1  : myFace2;
+
+      const CutFace& outFace = myCutFaces.Added( CutFace( faceOut ));
+      CutLink link1( nodesIn.back().Node(), nodesIn.back().Node(), faceOut );
+      CutLink link2( nodesIn.back().Node(), nodesIn.back().Node(), faceOut );
+
+      link1.myIntNode = nodesIn.back();
+      for ( size_t i = 0; i < nodesIn.size(); ++i )
+      {
+        link2.myIntNode = nodesIn[ i ];
+        outFace.AddEdge( link1, link2, faceIn, minNbOnPlane );
+        link1.myIntNode = link2.myIntNode;
+      }
+    }
+    else
+    {
+      // add parts of edges to a triangle including them
+
+      CutLink link1, link2;
+      IntPoint2D ip0, ip1;
+      ip0.myU[0] = ip0.myU[1] = 0.;
+      ip1.myU[0] = ip1.myU[1] = 1.;
+      ip0.myEdgeInd[0] = ip0.myEdgeInd[1] = ip1.myEdgeInd[0] = ip1.myEdgeInd[1] = 0;
+
+      for ( int isFromFace1 = 0; isFromFace1 < 2; ++isFromFace1 )
+      {
+        const SMDS_MeshElement*                faceTo = isFromFace1 ? myFace2  : myFace1;
+        const SMDS_MeshElement*              faceFrom = isFromFace1 ? myFace1  : myFace2;
+        const std::vector< SMESH_NodeXYZ >&   nodesTo = isFromFace1 ? myNodes2 : myNodes1;
+        const std::vector< SMESH_NodeXYZ >& nodesFrom = isFromFace1 ? myNodes1 : myNodes2;
+        const int                                 iTo = isFromFace1 ? 1 : 0;
+        const int                               iFrom = isFromFace1 ? 0 : 1;
+        //const int                       nbOnPlaneFrom = isFromFace1 ? myNbOnPlane1 : myNbOnPlane2;
+
+        const CutFace* cutFaceTo   = & myCutFaces.Added( CutFace( faceTo ));
+        // const CutFace* cutFaceFrom = 0;
+        // if ( nbOnPlaneFrom > minNbOnPlane )
+        //   cutFaceFrom = & myCutFaces.Added( CutFace( faceTo ));
+
+        link1.myFace = link2.myFace = faceTo;
+
+        IntPoint2DCompare ipCompare( iFrom );
+        TIntPointPtrSet pointsOnEdge( ipCompare ); // IntPoint2D sorted by parameter on edge
+
+        for ( size_t iE = 0; iE < nodesFrom.size(); ++iE )
+        {
+          // get parts of an edge iE
+
+          ip0.myEdgeInd[ iTo ] = iE;
+          ip1.myEdgeInd[ iTo ] = ( iE + 1 ) % nodesFrom.size();
+          ip0.myNode = nodesFrom[ ip0.myEdgeInd[ iTo ]];
+          ip1.myNode = nodesFrom[ ip1.myEdgeInd[ iTo ]];
+
+          pointsOnEdge.clear();
+
+          for ( int iP = 0; iP < nbIntPoints; ++iP )
+            if ( intPoints[ iP ].myEdgeInd[ iFrom ] == iE )
+              pointsOnEdge.insert( & intPoints[ iP ] );
+
+          pointsOnEdge.insert( pointsOnEdge.begin(), & ip0 );
+          pointsOnEdge.insert( pointsOnEdge.end(),   & ip1 );
+
+          // add edge parts to faceTo
+
+          TIntPointPtrSet::iterator ipIt = pointsOnEdge.begin() + 1;
+          for ( ; ipIt != pointsOnEdge.end(); ++ipIt )
+          {
+            const IntPoint2D* p1 = *(ipIt-1);
+            const IntPoint2D* p2 = *ipIt;
+            gp_XYZ middle = 0.5 * ( p1->myNode + p2->myNode );
+            if ( isPointInTriangle( middle, nodesTo ))
+            {
+              p1->InitLink( link1, iTo, ( p1 != & ip0 ) ? nodesTo : nodesFrom );
+              p2->InitLink( link2, iTo, ( p2 != & ip1 ) ? nodesTo : nodesFrom );
+              cutFaceTo->AddEdge( link1, link2, faceFrom, minNbOnPlane );
+
+              // if ( cutFaceFrom )
+              // {
+              //   p1->InitLink( link1, iFrom, nodesFrom );
+              //   p2->InitLink( link2, iFrom, nodesFrom );
+              //   cutFaceTo->AddEdge( link1, link2, faceTo, minNbOnPlane );
+              // }
+            }
+          }
+        }
+      }
+    }
+    return;
+
+  } // Intersector::cutCoplanar()
+
+  //================================================================================
+  /*!
+   * \brief Intersect edges added to myCutFaces
+   */
+  //================================================================================
+
+  void Intersector::intersectNewEdges( const CutFace& cf )
+  {
+    IntPoint2D intPoint;
+
+    if ( cf.NbInternalEdges() < 2 )
+      return;
+
+    const gp_XYZ& faceNorm = myNormals[ cf.myInitFace->GetID() ];
+    setPlaneIndices( faceNorm ); // choose indices on an axis-aligned plane
+
+    size_t limit = cf.myLinks.size() * cf.myLinks.size() * 2;
+
+    for ( size_t i1 = 3; i1 < cf.myLinks.size(); ++i1 )
+    {
+      if ( !cf.myLinks[i1].IsInternal() )
+        continue;
+
+      myIntPointSet.clear();
+      for ( size_t i2 = i1 + 2; i2 < cf.myLinks.size(); ++i2 )
+      {
+        if ( !cf.myLinks[i2].IsInternal() )
+          continue;
+
+        // prepare to intersection
+        myFace1     = cf.myLinks[i1].myFace;
+        myNodes1[0] = cf.myLinks[i1].myNode1;
+        myNodes1[1] = cf.myLinks[i1].myNode2;
+        myFace2     = cf.myLinks[i2].myFace;
+        myNodes2[0] = cf.myLinks[i2].myNode1;
+        myNodes2[1] = cf.myLinks[i2].myNode2;
+
+        // intersect
+        intPoint.myIsCollinear = true; // to find collinear solutions
+        if ( intersectEdgeEdge( 0, 0, intPoint ))
+        {
+          if ( cf.myLinks[i1].IsSame( cf.myLinks[i2] )) // remove i2
+          {
+            cf.myLinks[i1].ReplaceCoplanar( cf.myLinks[i2] );
+            cf.myLinks.erase( cf.myLinks.begin() + i2, cf.myLinks.begin() + i2 + 2 );
+            --i2;
+            continue;
+          }
+          if ( !intPoint.myIsCollinear )
+          {
+            intPoint.myEdgeInd[1] = i2;
+            myIntPointSet.insert( intPoint );
+          }
+          else // if ( intPoint.myIsCollinear ) // overlapping edges
+          {
+            myIntPointSet.clear(); // to recompute
+
+            if ( intPoint.myU[0] > intPoint.myU[1] ) // orient in same direction
+            {
+              std::swap( intPoint.myU[0], intPoint.myU[1] );
+              std::swap( myNodes1[0], myNodes1[1] );
+            }
+            // replace _COPLANAR by _INTERNAL
+            cf.myLinks[i1].ReplaceCoplanar( cf.myLinks[i1+1] );
+            cf.myLinks[i2].ReplaceCoplanar( cf.myLinks[i2+1] );
+
+            if ( coincide( myNodes1[0], myNodes2[0], myTol ) &&
+                 coincide( myNodes1[1], myNodes2[1], myTol ))
+            {
+              cf.myLinks.erase( cf.myLinks.begin() + i2, cf.myLinks.begin() + i2 + 2 );
+              --i2;
+              continue;
+            }
+
+            EdgePart common = cf.myLinks[i1];
+            common.ReplaceCoplanar( cf.myLinks[i2] );
+
+            const SMDS_MeshNode* n1 = myNodes1[0].Node(); // end nodes of an overlapping part
+            const SMDS_MeshNode* n2 = myNodes1[1].Node();
+            size_t i3 = cf.myLinks.size();
+
+            if ( myNodes1[0] != myNodes2[0] ) // a part before the overlapping one
+            {
+              if ( intPoint.myU[0] < 0 )
+                cf.myLinks[i1].Set( myNodes1[0].Node(), myNodes2[0].Node(),
+                                    cf.myLinks[i1].myFace, cf.myLinks[i1].myIndex );
+              else
+                cf.myLinks[i1].Set( myNodes2[0].Node(), myNodes1[0].Node(),
+                                    cf.myLinks[i2].myFace, cf.myLinks[i2].myIndex );
+
+              cf.myLinks[i1+1].Set( cf.myLinks[i1].myNode2,
+                                    cf.myLinks[i1].myNode1,
+                                    cf.myLinks[i1].myFace,
+                                    cf.myLinks[i1].myIndex);
+              n1 = cf.myLinks[i1].myNode2;
+            }
+            else
+              i3 = i1;
+
+            if ( myNodes1[1] != myNodes2[1] ) // a part after the overlapping one
+            {
+              if ( intPoint.myU[1] < 1 )
+                cf.myLinks[i2].Set( myNodes1[1].Node(), myNodes2[1].Node(),
+                                    cf.myLinks[i2].myFace, cf.myLinks[i2].myIndex );
+              else
+                cf.myLinks[i2].Set( myNodes2[1].Node(), myNodes1[1].Node(),
+                                    cf.myLinks[i1].myFace, cf.myLinks[i1].myIndex );
+
+              cf.myLinks[i2+1].Set( cf.myLinks[i2].myNode2,
+                                    cf.myLinks[i2].myNode1,
+                                    cf.myLinks[i2].myFace,
+                                    cf.myLinks[i2].myIndex);
+              n2 = cf.myLinks[i2].myNode1;
+            }
+            else
+              i3 = i2;
+
+            if ( i3 == cf.myLinks.size() )
+              cf.myLinks.resize( i3 + 2 );
+
+            cf.myLinks[i3].Set  ( n1, n2, common.myFace, common.myIndex );
+            cf.myLinks[i3+1].Set( n2, n1, common.myFace, common.myIndex );
+
+            i2 = i1 + 1; // recheck modified i1
+            continue;
+          }
+          //else
+          // {
+          //   // remember a new node
+          //   CutLink link1( myNodes1[0].Node(), myNodes1[1].Node(), cf.myInitFace );
+          //   CutLink link2( myNodes2[0].Node(), myNodes2[1].Node(), cf.myInitFace );
+          //   link2.myIntNode = link1.myIntNode = intPoint.myNode;
+          //   addLink( link1 );
+          //   addLink( link2 );
+
+          //   // split edges
+          //   size_t i = cf.myLinks.size();
+          //   if ( intPoint.myNode != cf.myLinks[ i1 ].myNode1 &&
+          //        intPoint.myNode != cf.myLinks[ i1 ].myNode2 )
+          //   {
+          //     cf.myLinks.push_back( cf.myLinks[ i1 ]);
+          //     cf.myLinks.push_back( cf.myLinks[ i1 + 1 ]);
+          //     cf.myLinks[ i1 ].myNode2 = cf.myLinks[ i1 + 1 ].myNode1 = intPoint.Node();
+          //     cf.myLinks[ i  ].myNode1 = cf.myLinks[ i  + 1 ].myNode2 = intPoint.Node();
+          //   }
+          //   if ( intPoint.myNode != cf.myLinks[ i2 ].myNode1 &&
+          //        intPoint.myNode != cf.myLinks[ i2 ].myNode2 )
+          //   {
+          //     i = cf.myLinks.size();
+          //     cf.myLinks.push_back( cf.myLinks[ i2 ]);
+          //     cf.myLinks.push_back( cf.myLinks[ i2 + 1 ]);
+          //     cf.myLinks[ i2 ].myNode2 = cf.myLinks[ i2 + 1 ].myNode1 = intPoint.Node();
+          //     cf.myLinks[ i  ].myNode1 = cf.myLinks[ i  + 1 ].myNode2 = intPoint.Node();
+          //   }
+          // }
+
+        } // if ( intersectEdgeEdge( 0, 0, intPoint ))
+
+        ++i2;
+        --limit;
+      }
+
+      // split i1 edge and all edges it intersects
+      // don't do it inside intersection loop in order not to loose direction of i1 edge
+      if ( !myIntPointSet.empty() )
+      {
+        cf.myLinks.reserve( cf.myLinks.size() + myIntPointSet.size() * 2 + 2 );
+
+        EdgePart* edge1 = &cf.myLinks[ i1 ];
+        EdgePart* twin1 = &cf.myLinks[ i1 + 1 ];
+
+        TIntPointSet::iterator ipIt = myIntPointSet.begin();
+        for ( ; ipIt != myIntPointSet.end(); ++ipIt ) // int points sorted on i1 edge
+        {
+          size_t i = cf.myLinks.size();
+          if ( ipIt->myNode != edge1->myNode1 &&
+               ipIt->myNode != edge1->myNode2 )
+          {
+            cf.myLinks.push_back( *edge1 );
+            cf.myLinks.push_back( *twin1 );
+            edge1->myNode2          = twin1->myNode1              = ipIt->Node();
+            cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = ipIt->Node();
+            edge1 = & cf.myLinks[ i ];
+            twin1 = & cf.myLinks[ i + 1 ];
+          }
+          size_t i2 = ipIt->myEdgeInd[1];
+          if ( ipIt->myNode != cf.myLinks[ i2 ].myNode1 &&
+               ipIt->myNode != cf.myLinks[ i2 ].myNode2 )
+          {
+            i = cf.myLinks.size();
+            cf.myLinks.push_back( cf.myLinks[ i2 ]);
+            cf.myLinks.push_back( cf.myLinks[ i2 + 1 ]);
+            cf.myLinks[ i2 ].myNode2 = cf.myLinks[ i2 + 1 ].myNode1 = ipIt->Node();
+            cf.myLinks[ i  ].myNode1 = cf.myLinks[ i  + 1 ].myNode2 = ipIt->Node();
+          }
+        }
+        if ( cf.myLinks.size() >= limit )
+          throw SALOME_Exception( "Infinite loop in Intersector::intersectNewEdges()" );
+      }
+      ++i1; // each internal edge encounters twice
+    }
+    return;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Split intersected faces
+   */
+  //================================================================================
+
+  void Intersector::MakeNewFaces( SMESH_MeshAlgos::TEPairVec& theNew2OldFaces,
+                                  SMESH_MeshAlgos::TNPairVec& theNew2OldNodes,
+                                  const double                theSign)
+  {
+    // unmark all nodes except intersection ones
+
+    for ( SMDS_NodeIteratorPtr nIt = myMesh->nodesIterator(); nIt->more(); )
+    {
+      const SMDS_MeshNode* n = nIt->next();
+      if ( n->isMarked() && n->GetID()-1 < (int) theNew2OldNodes.size() )
+        n->setIsMarked( false );
+    }
+    // SMESH_MeshAlgos::MarkElems( myMesh->nodesIterator(), false );
+
+    TCutLinkMap::const_iterator cutLinksIt = myCutLinks.cbegin();
+    // for ( ; cutLinksIt != myCutLinks.cend(); ++cutLinksIt )
+    // {
+    //   const CutLink& link = *cutLinksIt;
+    //   if ( link.IntNode() && link.IntNode()->GetID()-1 < (int) theNew2OldNodes.size() )
+    //     link.IntNode()->setIsMarked( true );
+    // }
+
+    // intersect edges added to myCutFaces
+
+    TCutFaceMap::const_iterator cutFacesIt = myCutFaces.cbegin();
+    for ( ; cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
+    {
+      const CutFace& cf = *cutFacesIt;
+      cf.ReplaceNodes( myRemove2KeepNodes );
+      intersectNewEdges( cf );
+    }
+
+    // make new faces
+
+    EdgeLoopSet                            loopSet;
+    SMESH_MeshAlgos::Triangulate           triangulator;
+    std::vector< EdgePart >                cutOffLinks;
+    TLinkMap                               cutOffCoplanarLinks;
+    std::vector< const CutFace* >          touchedFaces;
+    SMESH_MeshAlgos::TEPairVec::value_type new2OldTria;
+    CutFace                                cutFace(0);
+    std::vector< const SMDS_MeshNode* >    nodes;
+    std::vector<const SMDS_MeshElement *>  faces;
+
+    cutOffLinks.reserve( myCutFaces.Extent() * 2 );
+
+    for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
+    {
+      const CutFace& cf = *cutFacesIt;
+      if ( !cf.IsCut() )
+      {
+        touchedFaces.push_back( & cf );
+        continue;
+      }
+
+      const gp_XYZ& normal = myNormals[ cf.myInitFace->GetID() ];
+
+      // form loops of new faces
+      cf.ReplaceNodes( myRemove2KeepNodes );
+      cf.MakeLoops( loopSet, normal );
+
+      // avoid loops that are not connected to boundary edges of cf.myInitFace
+      if ( cf.RemoveInternalLoops( loopSet ))
+      {
+        intersectNewEdges( cf );
+        cf.MakeLoops( loopSet, normal );
+      }
+      // erase loops that are cut off by face intersections
+      cf.CutOffLoops( loopSet, theSign, myNormals, cutOffLinks, cutOffCoplanarLinks );
+
+      int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
+
+      const SMDS_MeshElement* tria;
+      for ( size_t iL = 0; iL < loopSet.myNbLoops; ++iL )
+      {
+        EdgeLoop& loop = loopSet.myLoops[ iL ];
+        if ( loop.myLinks.size() == 0 )
+          continue;
+
+        int nbTria  = triangulator.GetTriangles( &loop, nodes );
+        int nbNodes = 3 * nbTria;
+        for ( int i = 0; i < nbNodes; i += 3 )
+        {
+          if ( nodes[i] == nodes[i+1] || nodes[i] == nodes[i+2] || nodes[i+1] == nodes[i+2] )
+          {
+#ifdef _DEBUG_
+            std::cerr << "BAD tria" << std::endl;
+            cf.Dump();
+#endif
+            continue;
+          }
+          if (!( tria = myMesh->FindFace( nodes[i], nodes[i+1], nodes[i+2] )))
+            tria = myMesh->AddFace( nodes[i], nodes[i+1], nodes[i+2] );
+          tria->setIsMarked( true ); // not to remove it
+
+          new2OldTria = std::make_pair( tria, theNew2OldFaces[ index ].second );
+          if ( tria->GetID() < (int)theNew2OldFaces.size() )
+            theNew2OldFaces[ tria->GetID() ] = new2OldTria;
+          else
+            theNew2OldFaces.push_back( new2OldTria );
+
+          if ( index == tria->GetID() )
+            index = 0; // do not remove tria
+        }
+      }
+      theNew2OldFaces[ index ].first = 0;
+    }
+
+    // remove split faces
+    for ( size_t id = 1; id < theNew2OldFaces.size(); ++id )
+    {
+      if ( theNew2OldFaces[id].first )
+        continue;
+      if ( const SMDS_MeshElement* f = myMesh->FindElement( id ))
+        myMesh->RemoveFreeElement( f );
+    }
+
+    // remove face connected to cut off parts of cf.myInitFace
+
+    nodes.resize(2);
+    for ( size_t i = 0; i < cutOffLinks.size(); ++i )
+    {
+      //break;
+      nodes[0] = cutOffLinks[i].myNode1;
+      nodes[1] = cutOffLinks[i].myNode2;
+
+      if ( nodes[0] != nodes[1] &&
+           myMesh->GetElementsByNodes( nodes, faces ))
+      {
+        if ( cutOffLinks[i].myFace &&
+             cutOffLinks[i].myIndex != EdgePart::_COPLANAR &&
+             faces.size() == 2 )
+          continue;
+        for ( size_t iF = 0; iF < faces.size(); ++iF )
+        {
+          int index = faces[iF]->GetID();
+          // if ( //faces[iF]->isMarked()         ||  // kept part of cutFace
+          //      !theNew2OldFaces[ index ].first ) // already removed
+          //   continue;
+          cutFace.myInitFace = faces[iF];
+          // if ( myCutFaces.Contains( cutFace )) // keep cutting faces needed in CutOffLoops()
+          // {
+          //   if ( !myCutFaces.Added( cutFace ).IsCut() )
+          //     theNew2OldFaces[ index ].first = 0;
+          //   continue;
+          // }
+          cutFace.myLinks.clear();
+          cutFace.InitLinks();
+          for ( size_t iL = 0; iL < cutFace.myLinks.size(); ++iL )
+            if ( !cutOffLinks[i].IsSame( cutFace.myLinks[ iL ]))
+              cutOffLinks.push_back( cutFace.myLinks[ iL ]);
+
+          theNew2OldFaces[ index ].first = 0;
+          myMesh->RemoveFreeElement( faces[iF] );
+        }
+      }
+    }
+
+    // replace nodes in touched faces
+
+    // treat touched faces
+    for ( size_t i = 0; i < touchedFaces.size(); ++i )
+    {
+      const CutFace& cf = *touchedFaces[i];
+
+      int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
+      if ( !theNew2OldFaces[ index ].first )
+        continue; // already cut off
+
+      if ( !cf.ReplaceNodes( myRemove2KeepNodes ))
+        continue; // just keep as is
+
+      if ( cf.myLinks.size() == 3 )
+      {
+        const SMDS_MeshElement* tria = myMesh->AddFace( cf.myLinks[0].myNode1,
+                                                        cf.myLinks[1].myNode1,
+                                                        cf.myLinks[2].myNode1 );
+        new2OldTria = std::make_pair( tria, theNew2OldFaces[ index ].second );
+        if ( tria->GetID() < (int)theNew2OldFaces.size() )
+          theNew2OldFaces[ tria->GetID() ] = new2OldTria;
+        else
+          theNew2OldFaces.push_back( new2OldTria );
+      }
+      theNew2OldFaces[ index ].first = 0;
+    }
+
+
+    // add used new nodes to theNew2OldNodes
+    SMESH_MeshAlgos::TNPairVec::value_type new2OldNode;
+    new2OldNode.second = NULL;
+    for ( cutLinksIt = myCutLinks.cbegin(); cutLinksIt != myCutLinks.cend(); ++cutLinksIt )
+    {
+      const CutLink& link = *cutLinksIt;
+      if ( link.IntNode() ) // && link.IntNode()->NbInverseElements() > 0 )
+      {
+        new2OldNode.first = link.IntNode();
+        theNew2OldNodes.push_back( new2OldNode );
+      }
+    }
+
+    return;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Debug
+   */
+  //================================================================================
+
+  void CutFace::Dump() const
+  {
+    std::cout << std::endl << "INI F " << myInitFace->GetID() << std::endl;
+    for ( size_t i = 0; i < myLinks.size(); ++i )
+      std::cout << "[" << i << "] ("
+                << char(( myLinks[i].IsInternal() ? 'j' : '0' ) + myLinks[i].myIndex ) << ") "
+                << myLinks[i].myNode1->GetID() << " - " << myLinks[i].myNode2->GetID()
+                << " " << ( myLinks[i].myFace ? 'F' : 'C' )
+                << ( myLinks[i].myFace ? myLinks[i].myFace->GetID() : 0 ) << " " << std::endl;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Add an edge cutting this face
+   *  \param [in] p1 - start point of the edge
+   *  \param [in] p2 - end point of the edge
+   *  \param [in] cutter - a face producing the added cut edge.
+   *  \param [in] nbOnPlane - nb of triangle nodes lying on the plane of the cutter face
+   */
+  //================================================================================
+
+  void CutFace::AddEdge( const CutLink&          p1,
+                         const CutLink&          p2,
+                         const SMDS_MeshElement* cutterFace,
+                         const int               nbOnPlane,
+                         const int               iNotOnPlane) const
+  {
+    int iN[2] = { myInitFace->GetNodeIndex( p1.IntNode() ),
+                  myInitFace->GetNodeIndex( p2.IntNode() ) };
+    if ( iN[0] >= 0 && iN[1] >= 0 )
+    {
+      // the cutting edge is a whole side
+      if ((  cutterFace && nbOnPlane < 3 ) &&
+          !( cutterFace->GetNodeIndex( p1.IntNode() ) >= 0 &&
+             cutterFace->GetNodeIndex( p2.IntNode() ) >= 0 ))
+      {
+        InitLinks();
+        myLinks[ Abs( iN[0] - iN[1] ) == 1 ? Min( iN[0], iN[1] ) : 2 ].myFace = cutterFace;
+      }
+      return;
+    }
+
+    if ( p1.IntNode() == p2.IntNode() )
+    {
+      AddPoint( p1, p2, 1e-10 );
+      return;
+    }
+
+    InitLinks();
+
+    // cut side edges by a new one
+
+    int iEOnPlane = ( nbOnPlane == 2 ) ? ( iNotOnPlane + 1 ) % 3  :  -1;
+
+    double dist[2];
+    for ( int is2nd = 0; is2nd < 2; ++is2nd )
+    {
+      const CutLink& p = is2nd ? p2 : p1;
+      dist[ is2nd ] = 0;
+      if ( iN[ is2nd ] >= 0 )
+        continue;
+
+      int iE = Max( iEOnPlane, myInitFace->GetNodeIndex( p.Node1() ));
+      if ( iE < 0 )
+        continue; // link of other face
+
+      SMESH_NodeXYZ n0 = myLinks[iE].myNode1;
+      dist[ is2nd ]    = ( n0 - p.myIntNode ).SquareModulus();
+
+      for ( size_t i = 0; i < myLinks.size(); ++i )
+        if ( myLinks[i].myIndex == iE )
+        {
+          double d1 = n0.SquareDistance( myLinks[i].myNode1 );
+          if ( d1 < dist[ is2nd ] )
+          {
+            double d2 = n0.SquareDistance( myLinks[i].myNode2 );
+            if ( dist[ is2nd ] < d2 )
+            {
+              myLinks.push_back( myLinks[i] );
+              myLinks.back().myNode1 = myLinks[i].myNode2 = p.IntNode();
+              break;
+            }
+          }
+        }
+    }
+
+    int state = nbOnPlane == 3 ? EdgePart::_COPLANAR : EdgePart::_INTERNAL;
+
+    // look for an existing equal edge
+    if ( nbOnPlane == 2 )
+    {
+      SMESH_NodeXYZ n0 = myLinks[ iEOnPlane ].myNode1;
+      if ( iN[0] >= 0 ) dist[0] = ( n0 - p1.myIntNode ).SquareModulus();
+      if ( iN[1] >= 0 ) dist[1] = ( n0 - p2.myIntNode ).SquareModulus();
+      if ( dist[0] > dist[1] )
+        std::swap( dist[0], dist[1] );
+      for ( size_t i = 0; i < myLinks.size(); ++i )
+      {
+        if ( myLinks[i].myIndex != iEOnPlane )
+          continue;
+        gp_XYZ mid = 0.5 * ( SMESH_NodeXYZ( myLinks[i].myNode1 ) +
+                             SMESH_NodeXYZ( myLinks[i].myNode2 ));
+        double d = ( n0 - mid ).SquareModulus();
+        if ( dist[0] < d && d < dist[1] )
+          myLinks[i].myFace = cutterFace;
+      }
+      return;
+    }
+    else
+    {
+      EdgePart newEdge; newEdge.Set( p1.IntNode(), p2.IntNode(), cutterFace, state );
+      for ( size_t i = 0; i < myLinks.size(); ++i )
+      {
+        if ( myLinks[i].IsSame( newEdge ))
+        {
+          // if ( !myLinks[i].IsInternal() )
+          //   myLinks[ i ].myFace = cutterFace;
+          // else
+          myLinks[ i   ].ReplaceCoplanar( newEdge );
+          myLinks[ i+1 ].ReplaceCoplanar( newEdge );
+          return;
+        }
+        i += myLinks[i].IsInternal();
+      }
+    }
+
+    size_t  i = myLinks.size();
+    myLinks.resize( i + 2 );
+    myLinks[ i   ].Set( p1.IntNode(), p2.IntNode(), cutterFace, state );
+    myLinks[ i+1 ].Set( p2.IntNode(), p1.IntNode(), cutterFace, state );
+  }
+
+  //================================================================================
+  /*!
+   * \brief Add a point cutting this face
+   */
+  //================================================================================
+
+  void CutFace::AddPoint( const CutLink& p1, const CutLink& p2, double tol ) const
+  {
+    if ( myInitFace->GetNodeIndex( p1.IntNode() ) >= 0 ||
+         myInitFace->GetNodeIndex( p2.IntNode() ) >= 0 )
+      return;
+
+    InitLinks();
+
+    const CutLink* link = &p1;
+    int iE = myInitFace->GetNodeIndex( link->Node1() );
+    if ( iE < 0 )
+    {
+      link = &p2;
+      iE = myInitFace->GetNodeIndex( link->Node1() );
+    }
+    if ( iE >= 0 )
+    {
+      // cut an existing edge by the point
+      SMESH_NodeXYZ n0 = link->Node1();
+      double         d = ( n0 - link->myIntNode ).SquareModulus();
+
+      for ( size_t i = 0; i < myLinks.size(); ++i )
+        if ( myLinks[i].myIndex == iE )
+        {
+          double d1 = n0.SquareDistance( myLinks[i].myNode1 );
+          if ( d1 < d )
+          {
+            double d2 = n0.SquareDistance( myLinks[i].myNode2 );
+            if ( d < d2 )
+            {
+              myLinks.push_back( myLinks[i] );
+              myLinks.back().myNode1 = myLinks[i].myNode2 = link->IntNode();
+              return;
+            }
+          }
+        }
+    }
+    else // point is inside the triangle
+    {
+      // // check if a point already added
+      // for ( size_t i = 3; i < myLinks.size(); ++i )
+      //   if ( myLinks[i].myNode1 == p1.IntNode() )
+      //     return;
+
+      // // create a link between the point and the closest corner node
+      // const SMDS_MeshNode* closeNode = myLinks[0].myNode1;
+      // double minDist = p1.myIntNode.SquareDistance( closeNode );
+      // for ( int i = 1; i < 3; ++i )
+      // {
+      //   double dist = p1.myIntNode.SquareDistance( myLinks[i].myNode1 );
+      //   if ( dist < minDist )
+      //   {
+      //     minDist = dist;
+      //     closeNode = myLinks[i].myNode1;
+      //   }
+      // }
+      // if ( minDist > tol * tol )
+      // {
+      //   size_t i = myLinks.size();
+      //   myLinks.resize( i + 2 );
+      //   myLinks[ i   ].Set( p1.IntNode(), closeNode );
+      //   myLinks[ i+1 ].Set( closeNode, p1.IntNode() );
+      // }
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Perform node replacement
+   */
+  //================================================================================
+
+  bool CutFace::ReplaceNodes( const TNNMap& theRm2KeepMap ) const
+  {
+    bool replaced = false;
+    for ( size_t i = 0; i < myLinks.size(); ++i )
+    {
+      while ( theRm2KeepMap.IsBound((Standard_Address) myLinks[i].myNode1 ))
+        replaced = ( myLinks[i].myNode1 = theRm2KeepMap((Standard_Address) myLinks[i].myNode1 ));
+
+      while ( theRm2KeepMap.IsBound((Standard_Address) myLinks[i].myNode2 ))
+        replaced = ( myLinks[i].myNode2 = theRm2KeepMap((Standard_Address) myLinks[i].myNode2 ));
+    }
+
+    //if ( replaced ) // remove equal links
+    {
+      for ( size_t i1 = 0; i1 < myLinks.size(); ++i1 )
+      {
+        if ( myLinks[i1].myNode1 == myLinks[i1].myNode2 )
+        {
+          myLinks.erase( myLinks.begin() + i1,
+                         myLinks.begin() + i1 + 1 + myLinks[i1].IsInternal() );
+          --i1;
+          continue;
+        }
+        size_t i2 = i1 + 1 + myLinks[i1].IsInternal();
+        for ( ; i2 < myLinks.size(); ++i2 )
+        {
+          if ( !myLinks[i2].IsInternal() )
+            continue;
+          if ( myLinks[i1].IsSame( myLinks[i2] ))
+          {
+            myLinks[i1].  ReplaceCoplanar( myLinks[i2] );
+            if ( myLinks[i1].IsInternal() )
+              myLinks[i1+1].ReplaceCoplanar( myLinks[i2+1] );
+            if ( !myLinks[i1].myFace && myLinks[i2].myFace )
+            {
+              myLinks[i1].  myFace = myLinks[i2].myFace;
+              if ( myLinks[i1].IsInternal() )
+                myLinks[i1+1].myFace = myLinks[i2+1].myFace;
+            }
+            myLinks.erase( myLinks.begin() + i2,
+                           myLinks.begin() + i2 + 2 );
+            --i2;
+            continue;
+          }
+          ++i2;
+        }
+        i1 += myLinks[i1].IsInternal();
+      }
+    }
+
+    return replaced;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Initialize myLinks with edges of myInitFace
+   */
+  //================================================================================
+
+  void CutFace::InitLinks() const
+  {
+    if ( !myLinks.empty() ) return;
+
+    int nbNodes = myInitFace->NbNodes();
+    myLinks.reserve( nbNodes * 2 );
+    myLinks.resize( nbNodes );
+
+    for ( int i = 0; i < nbNodes; ++i )
+    {
+      const SMDS_MeshNode* n1 = myInitFace->GetNode( i );
+      const SMDS_MeshNode* n2 = myInitFace->GetNodeWrap( i + 1);
+      myLinks[i].Set( n1, n2, 0, i );
+    }
+  }
+  
+  //================================================================================
+  /*!
+   * \brief Return number of internal edges
+   */
+  //================================================================================
+
+  int CutFace::NbInternalEdges() const
+  {
+    int nb = 0;
+    for ( size_t i = 3; i < myLinks.size(); ++i )
+      nb += myLinks[i].IsInternal();
+
+    return nb / 2; // each internal edge encounters twice
+  }
+
+  //================================================================================
+  /*!
+   * \brief Remove loops that are not connected to boundary edges of myFace by
+   *        adding edges connecting these loops to the boundary
+   */
+  //================================================================================
+
+  bool CutFace::RemoveInternalLoops( EdgeLoopSet& theLoops ) const
+  {
+    size_t nbReachedLoops = 0;
+
+    // count loops including boundary EdgeParts
+    for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
+    {
+      EdgeLoop& loop = theLoops.myLoops[ iL ];
+
+      for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
+        if ( !loop.myLinks[ iE ]->IsInternal() )
+        {
+          nbReachedLoops += loop.SetConnected();
+          break;
+        }
+    }
+    if ( nbReachedLoops == theLoops.myNbLoops )
+      return false; // no unreachable loops
+
+
+    // try to reach all loops by propagating via internal edges shared by loops
+    size_t prevNbReached;
+    do
+    {
+      prevNbReached = nbReachedLoops;
+
+      for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
+      {
+        EdgeLoop& loop = theLoops.myLoops[ iL ];
+        if ( !loop.myIsBndConnected )
+          continue;
+
+        for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
+          if ( loop.myLinks[ iE ]->IsInternal() )
+          {
+            const EdgePart* twinEdge = getTwin( loop.myLinks[ iE ]);
+            EdgeLoop*          loop2 = theLoops.GetLoopOf( twinEdge );
+            if ( loop2->SetConnected() && ++nbReachedLoops == theLoops.myNbLoops )
+              return false; // no unreachable loops
+          }
+      }
+    }
+    while ( prevNbReached < nbReachedLoops );
+
+
+    // add links connecting internal loops with the boundary ones
+
+    for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
+    {
+      EdgeLoop& loop = theLoops.myLoops[ iL ];
+      if ( loop.myIsBndConnected )
+        continue;
+
+      // find a pair of closest nodes
+      const SMDS_MeshNode *closestNode1, *closestNode2;
+      double minDist = 1e100;
+      for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
+      {
+        SMESH_NodeXYZ n1 = loop.myLinks[ iE ]->myNode1;
+
+        for ( size_t i = 0; i < myLinks.size(); ++i )
+        {
+          if ( !loop.Contains( myLinks[i].myNode1 ))
+          {
+            double dist = n1.SquareDistance( myLinks[i].myNode1 );
+            if ( dist < minDist )
+            {
+              minDist = dist;
+              closestNode1 = loop.myLinks[ iE ]->myNode1;
+              closestNode2 = myLinks[i].myNode1;
+            }
+          }
+          if ( myLinks[i].IsInternal() )
+            ++i;
+        }
+      }
+
+      size_t i = myLinks.size();
+      myLinks.resize( i + 2 );
+      myLinks[ i   ].Set( closestNode1, closestNode2 );
+      myLinks[ i+1 ].Set( closestNode2, closestNode1 );
+    }
+
+    return true;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Return equal reversed edge
+   */
+  //================================================================================
+
+  EdgePart* CutFace::getTwin( const EdgePart* edge ) const
+  {
+    size_t i = edge - & myLinks[0];
+
+    if ( i > 2 && myLinks[ i-1 ].IsTwin( *edge ))
+      return & myLinks[ i-1 ];
+
+    if ( i+1 < myLinks.size() &&
+         myLinks[ i+1 ].IsTwin( *edge ))
+      return & myLinks[ i+1 ];
+
+    return 0;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Fill loops of edges
+   */
+  //================================================================================
+
+  void CutFace::MakeLoops( EdgeLoopSet& theLoops, const gp_XYZ& theFaceNorm ) const
+  {
+    theLoops.Init( myLinks );
+
+    if ( myLinks.size() == 3 )
+    {
+      theLoops.AddNewLoop();
+      theLoops.AddEdge( myLinks[0] );
+      theLoops.AddEdge( myLinks[1] );
+      theLoops.AddEdge( myLinks[2] );
+      return;
+    }
+
+    while ( !theLoops.AllEdgesUsed() )
+    {
+      theLoops.AddNewLoop();
+
+      // add 1st edge to a new loop
+      size_t i1;
+      for ( i1 = theLoops.myNbLoops - 1; i1 < myLinks.size(); ++i1 )
+        if ( theLoops.AddEdge( myLinks[i1] ))
+          break;
+
+      EdgePart*             lastEdge = & myLinks[ i1 ];
+      EdgePart*             twinEdge = getTwin( lastEdge );
+      const SMDS_MeshNode* firstNode = lastEdge->myNode1;
+      const SMDS_MeshNode*  lastNode = lastEdge->myNode2;
+
+      do // add the rest edges
+      {
+        theLoops.myCandidates.clear(); // edges starting at lastNode
+        int nbInternal = 0;
+
+        // find candidate edges
+        for ( size_t i = i1 + 1; i < myLinks.size(); ++i )
+          if ( myLinks[ i ].myNode1 == lastNode &&
+               &myLinks[ i ] != twinEdge &&
+               !theLoops.myIsUsedEdge[ i ])
+          {
+            theLoops.myCandidates.push_back( & myLinks[ i ]);
+            nbInternal += myLinks[ i ].IsInternal();
+          }
+
+        // choose among candidates
+        if ( theLoops.myCandidates.size() == 0 )
+        {
+          theLoops.GetLoopOf( lastEdge )->myHasPending = true;
+          lastEdge = twinEdge;
+        }
+        else if ( theLoops.myCandidates.size() == 1 )
+        {
+          lastEdge = theLoops.myCandidates[0];
+        }
+        else if ( nbInternal == 1 && !lastEdge->IsInternal() )
+        {
+          lastEdge = theLoops.myCandidates[ !theLoops.myCandidates[0]->IsInternal() ];
+        }
+        else
+        {
+          gp_Vec  lastVec = *lastEdge;
+          double maxAngle = -2 * M_PI;
+          for ( size_t i = 0; i < theLoops.myCandidates.size(); ++i )
+          {
+            double angle = lastVec.AngleWithRef( *theLoops.myCandidates[i], theFaceNorm );
+            if ( angle > maxAngle )
+            {
+              maxAngle = angle;
+              lastEdge = theLoops.myCandidates[i];
+            }
+          }
+        }
+        theLoops.AddEdge( *lastEdge );
+        lastNode = lastEdge->myNode2;
+        twinEdge = getTwin( lastEdge );
+      }
+      while ( lastNode != firstNode );
+
+    } // while ( !theLoops.AllEdgesUsed() )
+
+    return;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Erase loops that are cut off by face intersections
+   */
+  //================================================================================
+
+  void CutFace::CutOffLoops( EdgeLoopSet&                 theLoops,
+                             const double                 theSign,
+                             const std::vector< gp_XYZ >& theNormals,
+                             std::vector< EdgePart >&     theCutOffLinks,
+                             TLinkMap&                    theCutOffCoplanarLinks) const
+  {
+    EdgePart sideEdge;
+    for ( size_t i = 0; i < myLinks.size(); ++i )
+    {
+      if ( !myLinks[i].myFace )
+        continue;
+
+      EdgeLoop* loop = theLoops.GetLoopOf( & myLinks[i] );
+      if ( !loop || loop->myLinks.empty() || loop->myHasPending )
+        continue;
+
+      bool toErase, isCoplanar = ( myLinks[i].myIndex == EdgePart::_COPLANAR );
+
+      gp_Vec iniNorm = theNormals[ myInitFace->GetID() ];
+      if ( isCoplanar )
+      {
+        toErase = ( myLinks[i].myFace->GetID() > myInitFace->GetID() );
+
+        const EdgePart* twin = getTwin( & myLinks[i] );
+        if ( !twin || twin->myFace == myLinks[i].myFace )
+        {
+          // only one co-planar face includes myLinks[i]
+          gp_Vec inFaceDir = iniNorm ^ myLinks[i];
+          gp_XYZ   edgePnt = SMESH_NodeXYZ( myLinks[i].myNode1 );
+          for ( int iN = 0; iN < 3; ++iN )
+          {
+            gp_Vec inCutFaceDir = ( SMESH_NodeXYZ( myLinks[i].myFace->GetNode( iN )) - edgePnt );
+            if ( inCutFaceDir * inFaceDir < 0 )
+            {
+              toErase = false;
+              break;
+            }
+          }
+        }
+      }
+      else
+      {
+        gp_Vec   cutNorm = theNormals[ myLinks[i].myFace->GetID() ];
+        gp_Vec inFaceDir = iniNorm ^ myLinks[i];
+
+        toErase = inFaceDir * cutNorm * theSign < 0;
+        if ( !toErase )
+        {
+          // erase a neighboring loop
+          loop = 0;
+          if ( const EdgePart* twin = getTwin( & myLinks[i] ))
+            loop = theLoops.GetLoopOf( twin );
+          toErase = ( loop && !loop->myLinks.empty() );
+        }
+      }
+
+      if ( toErase )
+      {
+        if ( !isCoplanar )
+        {
+          // remember whole sides of myInitFace that are cut off
+          for ( size_t iE = 0; iE < loop->myLinks.size(); ++iE )
+          {
+            if ( !loop->myLinks[ iE ]->myFace              &&
+                 !loop->myLinks[ iE ]->IsInternal()     )//   &&
+                 // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked
+                 // !loop->myLinks[ iE ]->myNode2->isMarked() )
+            {
+              int i = loop->myLinks[ iE ]->myIndex;
+              sideEdge.Set( myInitFace->GetNode    ( i   ),
+                            myInitFace->GetNodeWrap( i+1 ));
+              theCutOffLinks.push_back( sideEdge );
+
+              if ( !sideEdge.IsSame( *loop->myLinks[ iE ] )) // nodes replaced
+              {
+                theCutOffLinks.push_back( *loop->myLinks[ iE ] );
+              }
+            }
+            else if ( IsCoplanar( loop->myLinks[ iE ]))
+            {
+              // propagate erasure to a co-planar face
+              theCutOffLinks.push_back( *loop->myLinks[ iE ]);
+            }
+            else if ( loop->myLinks[ iE ]->myFace &&
+                      loop->myLinks[ iE ]->IsInternal() )
+              theCutOffLinks.push_back( *loop->myLinks[ iE ]);
+          }
+
+          // clear the loop
+          theLoops.Erase( loop );
+        }
+      }
+    }
+    return;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Check if the face has cut edges
+   */
+  //================================================================================
+
+  bool CutFace::IsCut() const
+  {
+    if ( myLinks.size() > 3 )
+      return true;
+
+    if ( myLinks.size() == 3 )
+      for ( size_t i = 0; i < 3; ++i )
+        if ( myLinks[i].myFace )
+          return true;
+
+    return false;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Check if an edge is produced by a co-planar cut
+   */
+  //================================================================================
+
+  bool CutFace::IsCoplanar( const EdgePart* edge ) const
+  {
+    if ( edge->myIndex == EdgePart::_COPLANAR )
+    {
+      const EdgePart* twin = getTwin( edge );
+      return ( !twin || twin->myIndex == EdgePart::_COPLANAR );
+    }
+    return false;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Replace _COPLANAR cut edge by _INTERNAL oe vice versa
+   */
+  //================================================================================
+
+  bool EdgePart::ReplaceCoplanar( const EdgePart& e )
+  {
+    if ( myIndex + e.myIndex == _COPLANAR + _INTERNAL )
+    {
+      //check if the faces are connected
+      int nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e.myFace, myFace ).size();
+      bool toReplace = (( myIndex == _INTERNAL && nbCommonNodes > 1 ) ||
+                        ( myIndex == _COPLANAR && nbCommonNodes < 2 ));
+      if ( toReplace )
+      {
+        myIndex = e.myIndex;
+        myFace  = e.myFace;
+        return true;
+      }
+    }
+    return false;
+  }
+
+} // namespace
+
+//================================================================================
+/*!
+ * \brief Create an offsetMesh of given faces
+ *  \param [in] faceIt - the input faces
+ *  \param [out] new2OldFaces - history of faces
+ *  \param [out] new2OldNodes - history of nodes
+ *  \return SMDS_Mesh* - the new offset mesh, a caller should delete
+ */
+//================================================================================
+
+SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
+                                        SMDS_Mesh&           theSrcMesh,
+                                        const double         theOffset,
+                                        const bool           theFixIntersections,
+                                        TEPairVec&           theNew2OldFaces,
+                                        TNPairVec&           theNew2OldNodes)
+{
+  SMDS_Mesh* newMesh = new SMDS_Mesh;
+  theNew2OldFaces.clear();
+  theNew2OldNodes.clear();
+  theNew2OldFaces.push_back
+    ( std::make_pair(( const SMDS_MeshElement*) 0,
+                     ( const SMDS_MeshElement*) 0)); // to have index == face->GetID()
+
+  if ( theSrcMesh.GetMeshInfo().NbFaces( ORDER_QUADRATIC ) > 0 )
+    throw SALOME_Exception( "Offset of quadratic mesh not supported" );
+  if ( theSrcMesh.GetMeshInfo().NbFaces() > theSrcMesh.GetMeshInfo().NbTriangles() )
+    throw SALOME_Exception( "Offset of non-triangular mesh not supported" );
+
+  // copy input faces to the newMesh keeping IDs of nodes
+
+  double minNodeDist = 1e100;
+
+  std::vector< const SMDS_MeshNode* > nodes;
+  while ( theFaceIt->more() )
+  {
+    const SMDS_MeshElement* face = theFaceIt->next();
+    if ( face->GetType() != SMDSAbs_Face ) continue;
+
+    // copy nodes
+    nodes.assign( face->begin_nodes(), face->end_nodes() );
+    for ( size_t i = 0; i < nodes.size(); ++i )
+    {
+      const SMDS_MeshNode* newNode = newMesh->FindNode( nodes[i]->GetID() );
+      if ( !newNode )
+      {
+        SMESH_NodeXYZ xyz( nodes[i] );
+        newNode = newMesh->AddNodeWithID( xyz.X(), xyz.Y(), xyz.Z(), nodes[i]->GetID() );
+        theNew2OldNodes.push_back( std::make_pair( newNode, nodes[i] ));
+        nodes[i] = newNode;
+      }
+    }
+    const SMDS_MeshElement* newFace = 0;
+    switch ( face->GetEntityType() )
+    {
+    case SMDSEntity_Triangle:
+      newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2] );
+      break;
+    case SMDSEntity_Quad_Triangle:
+      newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],
+                                  nodes[3],nodes[4],nodes[5] );
+      break;
+    case SMDSEntity_BiQuad_Triangle:
+      newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],
+                                  nodes[3],nodes[4],nodes[5],nodes[6] );
+      break;
+    case SMDSEntity_Quadrangle:
+      newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3] );
+      break;
+    case SMDSEntity_Quad_Quadrangle:
+      newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3],
+                                  nodes[4],nodes[5],nodes[6],nodes[7] );
+      break;
+    case SMDSEntity_BiQuad_Quadrangle:
+      newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3],nodes[4],
+                                  nodes[5],nodes[6],nodes[7],nodes[8] );
+      break;
+    case SMDSEntity_Polygon:
+      newFace = newMesh->AddPolygonalFace( nodes );
+      break;
+    case SMDSEntity_Quad_Polygon:
+      newFace = newMesh->AddQuadPolygonalFace( nodes );
+      break;
+    default:
+      continue;
+    }
+    theNew2OldFaces.push_back( std::make_pair( newFace, face ));
+
+    SMESH_NodeXYZ pPrev = nodes.back(), p;
+    for ( size_t i = 0; i < nodes.size(); ++i )
+    {
+      p.Set( nodes[i] );
+      double dist = ( pPrev - p ).SquareModulus();
+      if ( dist > std::numeric_limits<double>::min() )
+        minNodeDist = dist;
+      pPrev = p;
+    }
+  } // while ( faceIt->more() )
+
+
+  // compute normals to faces
+  std::vector< gp_XYZ > normals( theNew2OldFaces.size() );
+  for ( size_t i = 1; i < normals.size(); ++i )
+  {
+    if ( !SMESH_MeshAlgos::FaceNormal( theNew2OldFaces[i].second, normals[i] ))
+      normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors
+  }
+
+  const double  tol = 1e-3 * Sqrt( minNodeDist );
+  const double sign = ( theOffset < 0 ? -1 : +1 );
+
+  // translate new nodes by normal to input faces
+  gp_XYZ newXYZ;
+  std::vector< const SMDS_MeshNode* > multiNormalNodes;
+  for ( size_t i = 0; i < theNew2OldNodes.size(); ++i )
+  {
+    const SMDS_MeshNode* newNode = theNew2OldNodes[i].first;
+
+    if ( getTranslatedPosition( newNode, theOffset, tol*10., sign, normals, theSrcMesh, newXYZ ))
+      newMesh->MoveNode( newNode, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
+    else
+      multiNormalNodes.push_back( newNode );
+  }
+  // make multi-normal translation
+  std::vector< SMESH_NodeXYZ > multiPos(10);
+  for ( size_t i = 0; i < multiNormalNodes.size(); ++i )
+  {
+    const SMDS_MeshNode* newNode = multiNormalNodes[i];
+    newNode->setIsMarked( true );
+    SMESH_NodeXYZ oldXYZ = newNode;
+    multiPos.clear();
+    for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
+    {
+      const SMDS_MeshElement* newFace = fIt->next();
+      const int             faceIndex = newFace->GetID();
+      const gp_XYZ&           oldNorm = normals[ faceIndex ];
+      const gp_XYZ             newXYZ = oldXYZ + oldNorm * theOffset;
+      if ( multiPos.empty() )
+      {
+        newMesh->MoveNode( newNode, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
+        multiPos.emplace_back( newNode );
+      }
+      else
+      {
+        newNode = 0;
+        for ( size_t iP = 0; iP < multiPos.size() &&  !newNode; ++iP )
+          if (( multiPos[iP] - newXYZ ).SquareModulus() < tol * tol )
+            newNode = multiPos[iP].Node();
+        if ( !newNode )
+        {
+          newNode = newMesh->AddNode( newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
+          newNode->setIsMarked( true );
+          theNew2OldNodes.push_back( std::make_pair( newNode, theNew2OldNodes[i].second ));
+          multiPos.emplace_back( newNode );
+        }
+      }
+      if ( newNode != oldXYZ.Node() )
+      {
+        nodes.assign( newFace->begin_nodes(), newFace->end_nodes() );
+        nodes[ newFace->GetNodeIndex( oldXYZ.Node() )] = newNode;
+        newMesh->ChangeElementNodes( newFace, & nodes[0], nodes.size() );
+      }
+    }
+  }
+
+  if ( !theFixIntersections )
+    return newMesh;
+
+
+  // remove new faces around concave nodes (they are marked) if the faces are inverted
+  gp_XYZ faceNorm;
+  for ( size_t i = 0; i < theNew2OldNodes.size(); ++i )
+  {
+    const SMDS_MeshNode* newNode = theNew2OldNodes[i].first;
+    //const SMDS_MeshNode* oldNode = theNew2OldNodes[i].second;
+    if ( newNode->isMarked() )
+    {
+      //gp_XYZ moveVec = sign * ( SMESH_NodeXYZ( newNode ) - SMESH_NodeXYZ( oldNode ));
+
+      //bool haveInverseFace = false;
+      for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
+      {
+        const SMDS_MeshElement* newFace = fIt->next();
+        const int             faceIndex = newFace->GetID();
+        const gp_XYZ&           oldNorm = normals[ faceIndex ];
+        if ( !SMESH_MeshAlgos::FaceNormal( newFace, faceNorm, /*normalize=*/false ) ||
+             //faceNorm * moveVec < 0 )
+             faceNorm * oldNorm < 0 )
+        {
+          //haveInverseFace = true;
+          theNew2OldFaces[ faceIndex ].first = 0;
+          newMesh->RemoveFreeElement( newFace );
+          //break;
+        }
+      }
+      // if ( haveInverseFace )
+      // {
+      //   newMesh->MoveNode( newNode, oldNode->X(), oldNode->Y(), oldNode->Z() );
+
+      //   for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
+      //   {
+      //     const SMDS_MeshElement* newFace = fIt->next();
+      //     if ( !SMESH_MeshAlgos::FaceNormal( newFace, normals[ newFace->GetID() ] ))
+      //       normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors
+      //   }
+      // }
+    }
+    // mark all new nodes located closer than theOffset from theSrcMesh
+  }
+
+  // ==================================================
+  // find self-intersections of new faces and fix them
+  // ==================================================
+
+  std::unique_ptr< SMESH_ElementSearcher > fSearcher
+    ( SMESH_MeshAlgos::GetElementSearcher( *newMesh, tol ));
+
+  Intersector intersector( newMesh, tol, normals );
+
+  std::vector< const SMDS_MeshElement* > closeFaces;
+  std::vector< const SMDS_MeshNode* >    faceNodes;
+  Bnd_B3d faceBox;
+  for ( size_t iF = 1; iF < theNew2OldFaces.size(); ++iF )
+  {
+    const SMDS_MeshElement* newFace = theNew2OldFaces[iF].first;
+    if ( !newFace ) continue;
+    faceNodes.assign( newFace->begin_nodes(), newFace->end_nodes() );
+
+    bool isConcaveNode1 = false;
+    for ( size_t iN = 0; iN < faceNodes.size() && !isConcaveNode1; ++iN )
+      isConcaveNode1 = faceNodes[iN]->isMarked();
+
+    // get faces close to a newFace
+    closeFaces.clear();
+    faceBox.Clear();
+    for ( size_t i = 0; i < faceNodes.size(); ++i )
+      faceBox.Add( SMESH_NodeXYZ( faceNodes[i] ));
+    faceBox.Enlarge( tol );
+
+    fSearcher->GetElementsInBox( faceBox, SMDSAbs_Face, closeFaces );
+
+    // intersect the newFace with closeFaces
+
+    for ( size_t i = 0; i < closeFaces.size(); ++i )
+    {
+      const SMDS_MeshElement* closeFace = closeFaces[i];
+      if ( closeFace->GetID() <= newFace->GetID() )
+        continue; // this pair already treated
+
+      // do not intersect connected faces if they have no concave nodes
+      int nbCommonNodes = 0;
+      for ( size_t iN = 0; iN < faceNodes.size(); ++iN )
+        nbCommonNodes += ( closeFace->GetNodeIndex( faceNodes[iN] ) >= 0 );
+
+      if ( !isConcaveNode1 )
+      {
+        bool isConcaveNode2 = false;
+        for ( SMDS_ElemIteratorPtr nIt = closeFace->nodesIterator(); nIt->more(); )
+          if (( isConcaveNode2 = nIt->next()->isMarked() ))
+            break;
+
+        if ( !isConcaveNode2 && nbCommonNodes > 0 )
+          continue;
+      }
+
+      intersector.Cut( newFace, closeFace, nbCommonNodes );
+    }
+  }
+  intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign );
+
+  return newMesh;
+}
diff --git a/src/SMESHUtils/SMESH_Triangulate.cxx b/src/SMESHUtils/SMESH_Triangulate.cxx
new file mode 100644 (file)
index 0000000..6652aea
--- /dev/null
@@ -0,0 +1,341 @@
+// 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_Triangulate.cxx
+// Created   : Thu Jan 18 18:00:13 2018
+// Author    : Edward AGAPOV (eap)
+
+// Extracted from ../DriverSTL/DriverSTL_W_SMDS_Mesh.cxx
+
+#include "SMESH_MeshAlgos.hxx"
+
+#include <Standard_ErrorHandler.hxx>
+#include <Standard_Failure.hxx>
+#include <gp_Ax2.hxx>
+
+using namespace SMESH_MeshAlgos;
+
+//================================================================================
+/*!
+ * \brief Initialization
+ */
+//================================================================================
+
+void Triangulate::PolyVertex::SetNodeAndNext( const SMDS_MeshNode* n,
+                                              PolyVertex&          v )
+{
+  _nxyz.Set( n );
+  _next = &v;
+  v._prev = this;
+}
+//================================================================================
+/*!
+ * \brief Remove self from a polygon
+ */
+//================================================================================
+
+Triangulate::PolyVertex* Triangulate::PolyVertex::Delete()
+{
+  _prev->_next = _next;
+  _next->_prev = _prev;
+  return _next;
+}
+
+//================================================================================
+/*!
+ * \brief Return nodes of a triangle
+ */
+//================================================================================
+
+void Triangulate::PolyVertex::GetTriaNodes( const SMDS_MeshNode** nodes) const
+{
+  nodes[0] = _prev->_nxyz._node;
+  nodes[1] =  this->_nxyz._node;
+  nodes[2] = _next->_nxyz._node;
+}
+
+//================================================================================
+/*!
+ * \brief Compute triangle area
+ */
+//================================================================================
+
+inline static double Area( const gp_XY& xy0, const gp_XY& xy1, const gp_XY& xy2 )
+{
+  gp_XY vPrev = xy0 - xy1;
+  gp_XY vNext = xy2 - xy1;
+  return vNext ^ vPrev;
+}
+
+//================================================================================
+/*!
+ * \brief Compute triangle area
+ */
+//================================================================================
+
+double Triangulate::PolyVertex::TriaArea() const
+{
+  return Area( _prev->_xy, this->_xy, _next->_xy );
+}
+
+//================================================================================
+/*!
+ * \brief Check if a vertex is inside a triangle
+ */
+//================================================================================
+
+bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v )
+{
+  if ( this ->_nxyz == v->_nxyz ||
+       _prev->_nxyz == v->_nxyz ||
+       _next->_nxyz == v->_nxyz )
+    return false;
+
+  gp_XY p = _prev->_xy - v->_xy;
+  gp_XY t =  this->_xy - v->_xy;
+  gp_XY n = _next->_xy - v->_xy;
+  const double tol = -1e-12;
+  return (( p ^ t ) >= tol &&
+          ( t ^ n ) >= tol &&
+          ( n ^ p ) >= tol );
+  // return ( Area( _prev, this, v ) > 0 &&
+  //          Area( this, _next, v ) > 0 &&
+  //          Area( _next, _prev, v ) > 0 );
+}
+
+//================================================================================
+/*!
+ * \brief Triangulate a polygon. Assure correct orientation for concave polygons
+ */
+//================================================================================
+
+bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
+                               const size_t                        nbNodes )
+{
+  // connect nodes into a ring
+  _pv.resize( nbNodes );
+  for ( size_t i = 1; i < nbNodes; ++i )
+    _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i] );
+  _pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], _pv[0] );
+
+  // get a polygon normal
+  gp_XYZ normal(0,0,0), p0,v01,v02;
+  p0  = _pv[0]._nxyz;
+  v01 = _pv[1]._nxyz - p0;
+  for ( size_t i = 2; i < nbNodes; ++i )
+  {
+    v02 = _pv[i]._nxyz - p0;
+    normal += v01 ^ v02;
+    v01 = v02;
+  }
+  // project nodes to the found plane
+  gp_Ax2 axes;
+  try {
+    axes = gp_Ax2( p0, normal, v01 );
+  }
+  catch ( Standard_Failure ) {
+    return false;
+  }
+  for ( size_t i = 0; i < nbNodes; ++i )
+  {
+    gp_XYZ p = _pv[i]._nxyz - p0;
+    _pv[i]._xy.SetX( axes.XDirection().XYZ() * p );
+    _pv[i]._xy.SetY( axes.YDirection().XYZ() * p );
+  }
+
+  // in a loop, find triangles with positive area and having no vertices inside
+  int iN = 0, nbTria = nbNodes - 2;
+  nodes.reserve( nbTria * 3 );
+  const double minArea = 1e-6;
+  PolyVertex* v = &_pv[0], *vi;
+  int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
+  while ( nbBadTria < nbVertices )
+  {
+    if (( isGoodTria = v->TriaArea() > minArea ))
+    {
+      for ( vi = v->_next->_next;
+            vi != v->_prev;
+            vi = vi->_next )
+      {
+        if ( v->IsInsideTria( vi ))
+          break;
+      }
+      isGoodTria = ( vi == v->_prev );
+    }
+    if ( isGoodTria )
+    {
+      v->GetTriaNodes( &nodes[ iN ] );
+      iN += 3;
+      v = v->Delete();
+      if ( --nbVertices == 3 )
+      {
+        // last triangle remains
+        v->GetTriaNodes( &nodes[ iN ] );
+        return true;
+      }
+      nbBadTria = 0;
+    }
+    else
+    {
+      v = v->_next;
+      ++nbBadTria;
+    }
+  }
+
+  // the polygon is invalid; add triangles with positive area
+  nbBadTria = 0;
+  while ( nbBadTria < nbVertices )
+  {
+    isGoodTria = v->TriaArea() > minArea;
+    if ( isGoodTria )
+    {
+      v->GetTriaNodes( &nodes[ iN ] );
+      iN += 3;
+      v = v->Delete();
+      if ( --nbVertices == 3 )
+      {
+        // last triangle remains
+        v->GetTriaNodes( &nodes[ iN ] );
+        return true;
+      }
+      nbBadTria = 0;
+    }
+    else
+    {
+      v = v->_next;
+      ++nbBadTria;
+    }
+  }
+
+  // add all the rest triangles
+  while ( nbVertices >= 3 )
+  {
+    v->GetTriaNodes( &nodes[ iN ] );
+    iN += 3;
+    v = v->Delete();
+    --nbVertices;
+  }
+
+  return true;
+
+} // triangulate()
+
+//================================================================================
+/*!
+ * \brief Return nb triangles in a decomposed mesh face
+ *  \retval int - number of triangles
+ */
+//================================================================================
+
+int Triangulate::GetNbTriangles( const SMDS_MeshElement* face )
+{
+  // WARNING: counting triangles must be coherent with GetTriangles()
+  switch ( face->GetEntityType() )
+  {
+  case SMDSEntity_BiQuad_Triangle:
+  case SMDSEntity_BiQuad_Quadrangle:
+    return face->NbNodes() - 1;
+    // case SMDSEntity_Triangle:
+    // case SMDSEntity_Quad_Triangle:
+    // case SMDSEntity_Quadrangle:
+    // case SMDSEntity_Quad_Quadrangle:
+    // case SMDSEntity_Polygon:
+    // case SMDSEntity_Quad_Polygon:
+  default:
+    return face->NbNodes() - 2;
+  }
+  return 0;
+}
+
+//================================================================================
+/*!
+ * \brief Decompose a mesh face into triangles
+ *  \retval int - number of triangles
+ */
+//================================================================================
+
+int Triangulate::GetTriangles( const SMDS_MeshElement*             face,
+                               std::vector< const SMDS_MeshNode*>& nodes)
+{
+  if ( face->GetType() != SMDSAbs_Face )
+    return 0;
+
+  // WARNING: decomposing into triangles must be coherent with getNbTriangles()
+  int nbTria, i = 0, nbNodes = face->NbNodes();
+  SMDS_NodeIteratorPtr nIt = face->interlacedNodesIterator();
+  nodes.resize( nbNodes * 3 );
+  nodes[ i++ ] = nIt->next();
+  nodes[ i++ ] = nIt->next();
+
+  const SMDSAbs_EntityType type = face->GetEntityType();
+  switch ( type )
+  {
+  case SMDSEntity_BiQuad_Triangle:
+  case SMDSEntity_BiQuad_Quadrangle:
+
+    nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8;
+    nodes[ i++ ] = face->GetNode( nbTria );
+    for ( i = 3; i < 3*(nbTria-1); i += 3 )
+    {
+      nodes[ i+0 ] = nodes[ i-2 ];
+      nodes[ i+1 ] = nIt->next();
+      nodes[ i+2 ] = nodes[ 2 ];
+    }
+    nodes[ i+0 ] = nodes[ i-2 ];
+    nodes[ i+1 ] = nodes[ 0 ];
+    nodes[ i+2 ] = nodes[ 2 ];
+    break;
+
+  case SMDSEntity_Triangle:
+
+    nbTria = 1;
+    nodes[ i++ ] = nIt->next();
+    break;
+
+  default:
+
+    // case SMDSEntity_Quad_Triangle:
+    // case SMDSEntity_Quadrangle:
+    // case SMDSEntity_Quad_Quadrangle:
+    // case SMDSEntity_Polygon:
+    // case SMDSEntity_Quad_Polygon:
+
+    nbTria = nbNodes - 2;
+    while ( nIt->more() )
+      nodes[ i++ ] = nIt->next();
+
+    if ( nbTria > 1 && !triangulate( nodes, nbNodes ))
+    {
+      nIt = face->interlacedNodesIterator();
+      nodes[ 0 ] = nIt->next();
+      nodes[ 1 ] = nIt->next();
+      nodes[ 2 ] = nIt->next();
+      for ( i = 3; i < 3*nbTria; i += 3 )
+      {
+        nodes[ i+0 ] = nodes[ 0 ];
+        nodes[ i+1 ] = nodes[ i-1 ];
+        nodes[ i+2 ] = nIt->next();
+      }
+    }
+  }
+
+  return nbTria;
+}
diff --git a/src/SMESHUtils/SMESH_Triangulate.hxx b/src/SMESHUtils/SMESH_Triangulate.hxx
new file mode 100644 (file)
index 0000000..8734c4b
--- /dev/null
@@ -0,0 +1,51 @@
+// 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_Triangulate.hxx
+// Created   : Thu Jan 18 17:51:34 2018
+// Author    : Edward AGAPOV (eap)
+
+
+#ifndef __SMESH_Triangulate_HXX__
+#define __SMESH_Triangulate_HXX__
+
+/*!
+ * \brief Divide a mesh face into triangles
+ */
+class SMESHUtils_EXPORT SMESH_Triangulate
+{
+ public:
+
+  static int GetNbTriangles( const SMDS_MeshElement* face );
+
+  int GetTriangles( const SMDS_MeshElement*             face,
+                    std::vector< const SMDS_MeshNode*>& nodes);
+
+ private:
+
+  bool triangulate( std::vector< const SMDS_MeshNode*>& nodes, const size_t nbNodes );
+
+  struct PolyVertex;
+  std::vector< PolyVertex > _pv;
+};
+
+
+#endif