src/SMESHUtils/SMESH_Offset.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 // File      : SMESH_Offset.cxx
  23 // Created   : Mon Dec 25 15:52:38 2017
  24 // Author    : Edward AGAPOV (eap)
  25
  26 #include "SMESH_MeshAlgos.hxx"
  27
  28 #include <SMDS_PolygonalFaceOfNodes.hxx>
  29 #include "SMDS_Mesh.hxx"
  30
  31 #include <Utils_SALOME_Exception.hxx>
  32 #include <Basics_OCCTVersion.hxx>
  33
  34 #include <Bnd_B3d.hxx>
  35 #include <NCollection_Map.hxx>
  36 #include <gp_Lin.hxx>
  37 #include <gp_Pln.hxx>
  38
  39 #include <boost/container/flat_set.hpp>
  40 #include <boost/dynamic_bitset.hpp>
  41
  42 namespace
  43 {
  44   const int theMaxNbFaces = 256; // max number of faces sharing a node
  45
  46   typedef NCollection_DataMap< const SMDS_MeshNode*, const SMDS_MeshNode*, SMESH_Hasher > TNNMap;
  47   typedef NCollection_Map< SMESH_Link, SMESH_TLinkHasher >                                TLinkMap;
  48
  49   //--------------------------------------------------------------------------------
  50   /*!
  51    * \brief Intersected face side storing a node created at this intersection
  52    *        and an intersected face
  53    */
  54   struct CutLink
  55   {
  56     bool                     myReverse;
  57     const SMDS_MeshNode*     myNode[2]; // side nodes. WARNING: don't set them directly, use Set()
  58     mutable SMESH_NodeXYZ    myIntNode; // intersection node
  59     const SMDS_MeshElement*  myFace;    // cutter face
  60     int                      myIndex;   // index of a node on the same link
  61
  62     CutLink(const SMDS_MeshNode*    node1=0,
  63             const SMDS_MeshNode*    node2=0,
  64             const SMDS_MeshElement* face=0,
  65             const int               index=0) { Set ( node1, node2, face, index ); }
  66
  67     void Set( const SMDS_MeshNode*    node1,
  68               const SMDS_MeshNode*    node2,
  69               const SMDS_MeshElement* face,
  70               const int               index=0)
  71     {
  72       myNode[0] = node1; myNode[1] = node2; myFace = face; myIndex = index; myReverse = false;
  73       if ( myNode[0] && ( myReverse = ( myNode[0]->GetID() > myNode[1]->GetID() )))
  74         std::swap( myNode[0], myNode[1] );
  75     }
  76     const SMDS_MeshNode* IntNode() const { return myIntNode.Node(); }
  77     const SMDS_MeshNode* Node1() const { return myNode[ myReverse ]; }
  78     const SMDS_MeshNode* Node2() const { return myNode[ !myReverse ]; }
  79   };
  80
  81   struct CutLinkHasher
  82   {
  83 #if OCC_VERSION_LARGE < 0x07080000
  84     static Standard_Integer HashCode(const CutLink&         link,
  85                                      const Standard_Integer upper)
  86     {
  87       Standard_Integer n = ( link.myNode[0]->GetID() +
  88                              link.myNode[1]->GetID() +
  89                              link.myIndex );
  90       return ::HashCode( n, upper );
  91     }
  92     static Standard_Boolean IsEqual(const CutLink& link1, const CutLink& link2 )
  93     {
  94       return ( link1.myNode[0] == link2.myNode[0] &&
  95                link1.myNode[1] == link2.myNode[1] &&
  96                link1.myIndex == link2.myIndex );
  97     }
  98 #else
  99     size_t operator()(const CutLink& link) const
 100     {
 101       return size_t( link.myNode[0]->GetID() +
 102                      link.myNode[1]->GetID() +
 103                      link.myIndex );
 104     }
 105     bool operator()(const CutLink& link1, const CutLink& link2 ) const
 106     {
 107       return ( link1.myNode[0] == link2.myNode[0] &&
 108                link1.myNode[1] == link2.myNode[1] &&
 109                link1.myIndex == link2.myIndex );
 110     }
 111 #endif
 112   };
 113
 114   typedef NCollection_Map< CutLink, CutLinkHasher > TCutLinkMap;
 115
 116   //--------------------------------------------------------------------------------
 117   /*!
 118    * \brief Part of a divided face edge
 119    */
 120   struct EdgePart
 121   {
 122     const SMDS_MeshNode*    myNode1;
 123     const SMDS_MeshNode*    myNode2;
 124     int                     myIndex; // positive -> side index, negative -> State
 125     const SMDS_MeshElement* myFace;
 126
 127     enum State { _INTERNAL = -1, _COPLANAR = -2, _PENDING = -3 };
 128
 129     void Set( const SMDS_MeshNode*    Node1,
 130               const SMDS_MeshNode*    Node2,
 131               const SMDS_MeshElement* Face = 0,
 132               int                     EdgeIndex = _INTERNAL )
 133     { myNode1 = Node1; myNode2 = Node2; myIndex = EdgeIndex; myFace = Face; }
 134
 135     // bool HasSameNode( const EdgePart& other ) { return ( myNode1 == other.myNode1 ||
 136     //                                                      myNode1 == other.myNode2 ||
 137     //                                                      myNode2 == other.myNode1 ||
 138     //                                                      myNode2 == other.myNode2 );
 139     // }
 140     bool IsInternal() const { return myIndex < 0; }
 141     bool IsTwin( const EdgePart& e ) const { return myNode1 == e.myNode2 && myNode2 == e.myNode1; }
 142     bool IsSame( const EdgePart& e ) const {
 143       return (( myNode1 == e.myNode2 && myNode2 == e.myNode1 ) ||
 144               ( myNode1 == e.myNode1 && myNode2 == e.myNode2 )); }
 145     bool ReplaceCoplanar( const EdgePart& e );
 146     operator SMESH_Link() const { return SMESH_Link( myNode1, myNode2 ); }
 147     operator gp_Vec() const { return SMESH_NodeXYZ( myNode2 ) - SMESH_NodeXYZ( myNode1 ); }
 148   };
 149
 150   //--------------------------------------------------------------------------------
 151   /*!
 152    * \brief Loop of EdgePart's forming a new face which is a part of CutFace
 153    */
 154   struct EdgeLoop : public SMDS_PolygonalFaceOfNodes
 155   {
 156     std::vector< const EdgePart* > myLinks;
 157     bool                           myIsBndConnected; //!< is there a path to CutFace side edges
 158     bool                           myHasPending;     //!< an edge encounters twice
 159
 160     EdgeLoop() : SMDS_PolygonalFaceOfNodes( std::vector<const SMDS_MeshNode *>() ) {}
 161     void Clear() { myLinks.clear(); myIsBndConnected = false; myHasPending = false; }
 162     bool SetConnected() { bool was = myIsBndConnected; myIsBndConnected = true; return !was; }
 163     size_t Contains( const SMDS_MeshNode* n ) const
 164     {
 165       for ( size_t i = 0; i < myLinks.size(); ++i )
 166         if ( myLinks[i]->myNode1 == n ) return i + 1;
 167       return 0;
 168     }
 169     virtual int NbNodes() const { return myLinks.size(); }
 170     virtual SMDS_ElemIteratorPtr nodesIterator() const
 171     {
 172       return setNodes(), SMDS_PolygonalFaceOfNodes::nodesIterator();
 173     }
 174     virtual SMDS_NodeIteratorPtr nodeIterator() const
 175     {
 176       return setNodes(), SMDS_PolygonalFaceOfNodes::nodeIterator();
 177     }
 178     void setNodes() const //!< set nodes to SMDS_PolygonalFaceOfNodes
 179     {
 180       EdgeLoop* me = const_cast<EdgeLoop*>( this );
 181       me->myNodes.resize( NbNodes() );
 182       size_t iMin = 0;
 183       for ( size_t i = 1; i < myNodes.size(); ++i ) {
 184         if ( myLinks[ i ]->myNode1->GetID() < myLinks[ iMin ]->myNode1->GetID() )
 185           iMin = i;
 186       }
 187       for ( size_t i = 0; i < myNodes.size(); ++i )
 188         me->myNodes[ i ] = myLinks[ ( iMin + i ) % myNodes.size() ]->myNode1;
 189     }
 190   };
 191
 192   //--------------------------------------------------------------------------------
 193   /*!
 194    * \brief Set of EdgeLoop's constructed from a CutFace
 195    */
 196   struct EdgeLoopSet
 197   {
 198     std::vector< EdgeLoop >  myLoops;       //!< buffer of EdgeLoop's
 199     size_t                   myNbLoops;     //!< number of constructed loops
 200
 201     const EdgePart*          myEdge0;       //!< & CutFace.myLinks[0]
 202     size_t                   myNbUsedEdges; //!< nb of EdgePart's added to myLoops
 203     boost::dynamic_bitset<>  myIsUsedEdge;  //!< is i-th EdgePart of CutFace is in any EdgeLoop
 204     std::vector< EdgeLoop* > myLoopOfEdge;  //!< EdgeLoop of CutFace.myLinks[i]
 205     std::vector< EdgePart* > myCandidates;  //!< EdgePart's starting at the same node
 206
 207     EdgeLoopSet(): myLoops(100) {}
 208
 209     void Init( const std::vector< EdgePart >& edges )
 210     {
 211       size_t nb = edges.size();
 212       myEdge0 = & edges[0];
 213       myNbLoops = 0;
 214       myNbUsedEdges = 0;
 215       myIsUsedEdge.reset();
 216       myIsUsedEdge.resize( nb, false );
 217       myLoopOfEdge.clear();
 218       myLoopOfEdge.resize( nb, (EdgeLoop*) 0 );
 219     }
 220     EdgeLoop& AddNewLoop()
 221     {
 222       if ( ++myNbLoops >= myLoops.size() )
 223         myLoops.resize( myNbLoops + 10 );
 224       myLoops[ myNbLoops-1 ].Clear();
 225       return myLoops[ myNbLoops-1 ];
 226     }
 227     bool AllEdgesUsed() const { return myNbUsedEdges == myLoopOfEdge.size(); }
 228
 229     bool AddEdge( EdgePart& edge )
 230     {
 231       size_t i = Index( edge );
 232       if ( myIsUsedEdge[ i ])
 233         return false;
 234       myLoops[ myNbLoops-1 ].myLinks.push_back( &edge );
 235       myLoopOfEdge[ i ] = & myLoops[ myNbLoops-1 ];
 236       myIsUsedEdge[ i ] = true;
 237       ++myNbUsedEdges;
 238       return true;
 239     }
 240     void Erase( EdgeLoop* loop )
 241     {
 242       for ( size_t iE = 0; iE < loop->myLinks.size(); ++iE )
 243         myLoopOfEdge[ Index( *loop->myLinks[ iE ] )] = 0;
 244       loop->Clear();
 245     }
 246     void Join( EdgeLoop& loop1, size_t iAfterConcact,
 247                EdgeLoop& loop2, size_t iFromEdge2 )
 248     {
 249       std::vector< const EdgePart* > linksAfterContact( loop1.myLinks.begin() + iAfterConcact,
 250                                                         loop1.myLinks.end() );
 251       loop1.myLinks.reserve( loop2.myLinks.size() + loop1.myLinks.size() );
 252       loop1.myLinks.resize( iAfterConcact );
 253       loop1.myLinks.insert( loop1.myLinks.end(),
 254                             loop2.myLinks.begin() + iFromEdge2, loop2.myLinks.end() );
 255       loop1.myLinks.insert( loop1.myLinks.end(),
 256                             loop2.myLinks.begin(), loop2.myLinks.begin() + iFromEdge2 );
 257       loop1.myLinks.insert( loop1.myLinks.end(),
 258                             linksAfterContact.begin(), linksAfterContact.end() );
 259       loop1.myIsBndConnected = loop2.myIsBndConnected;
 260       loop2.Clear();
 261       for ( size_t iE = 0; iE < loop1.myLinks.size(); ++iE )
 262         myLoopOfEdge[ Index( *loop1.myLinks[ iE ] )] = & loop1;
 263     }
 264     size_t    Index( const EdgePart& edge ) const { return &edge - myEdge0; }
 265     EdgeLoop* GetLoopOf( const EdgePart* edge ) { return myLoopOfEdge[ Index( *edge )]; }
 266   };
 267
 268   //--------------------------------------------------------------------------------
 269   /*!
 270    * \brief Intersections of a face
 271    */
 272   struct CutFace
 273   {
 274     mutable std::vector< EdgePart > myLinks;
 275     const SMDS_MeshElement*         myInitFace;
 276
 277     CutFace( const SMDS_MeshElement* face ): myInitFace( face ) {}
 278     void AddEdge( const CutLink&          p1,
 279                   const CutLink&          p2,
 280                   const SMDS_MeshElement* cutter,
 281                   const int               nbOnPlane,
 282                   const int               iNotOnPlane = -1) const;
 283     void AddPoint( const CutLink& p1, const CutLink& p2, double tol ) const;
 284     bool ReplaceNodes( const TNNMap& theRm2KeepMap ) const;
 285     bool IsCut() const;
 286     int  NbInternalEdges() const;
 287     void MakeLoops( EdgeLoopSet& loops, const gp_XYZ& theFaceNorm ) const;
 288     bool RemoveInternalLoops( EdgeLoopSet& theLoops ) const;
 289     void CutOffLoops( EdgeLoopSet&                 theLoops,
 290                       const double                 theSign,
 291                       const std::vector< gp_XYZ >& theNormals,
 292                       std::vector< EdgePart >&     theCutOffLinks,
 293                       TLinkMap&                    theCutOffCoplanarLinks) const;
 294     void InitLinks() const;
 295     bool IsCoplanar( const EdgePart* edge ) const;
 296     void Dump() const;
 297
 298   private:
 299
 300     EdgePart* getTwin( const EdgePart* edge ) const;
 301   };
 302
 303   struct CutFaceHasher
 304   {
 305 #if OCC_VERSION_LARGE < 0x07080000
 306     static Standard_Integer HashCode(const CutFace& f, const Standard_Integer upper)
 307     {
 308       return ::HashCode( FromSmIdType<int>(f.myInitFace->GetID()), upper );
 309     }
 310     static Standard_Boolean IsEqual(const CutFace& f1, const CutFace& f2 )
 311     {
 312       return f1.myInitFace == f2.myInitFace;
 313     }
 314 #else
 315     size_t operator()(const CutFace& f) const
 316     {
 317       return FromSmIdType<int>(f.myInitFace->GetID());
 318     }
 319
 320     bool operator()(const CutFace& f1, const CutFace& f2) const
 321     {
 322       return f1.myInitFace == f2.myInitFace;
 323     }
 324 #endif
 325   };
 326
 327   typedef NCollection_Map< CutFace, CutFaceHasher > TCutFaceMap;
 328
 329   //--------------------------------------------------------------------------------
 330   /*!
 331    * \brief Intersection point of two edges of co-planar triangles
 332    */
 333   struct IntPoint2D
 334   {
 335     size_t        myEdgeInd[2]; //!< edge indices of triangles
 336     double        myU      [2]; //!< parameter [0,1] on edges of triangles
 337     SMESH_NodeXYZ myNode;       //!< intersection node
 338     bool          myIsCollinear;//!< edges are collinear
 339
 340     IntPoint2D() : myIsCollinear( false ) {}
 341
 342     void InitLink( CutLink& link, int iFace, const std::vector< SMESH_NodeXYZ >& nodes ) const
 343     {
 344       link.Set( nodes[  myEdgeInd[ iFace ]                      ].Node(),
 345                 nodes[( myEdgeInd[ iFace ] + 1 ) % nodes.size() ].Node(),
 346                 link.myFace );
 347       link.myIntNode = myNode;
 348     }
 349     const SMDS_MeshNode* Node() const { return myNode.Node(); }
 350   };
 351   struct IntPoint2DCompare
 352   {
 353     int myI;
 354     IntPoint2DCompare( int iFace=0 ): myI( iFace ) {}
 355     bool operator() ( const IntPoint2D* ip1, const IntPoint2D* ip2 ) const
 356     {
 357       return ip1->myU[ myI ] < ip2->myU[ myI ];
 358     }
 359     bool operator() ( const IntPoint2D& ip1, const IntPoint2D& ip2 ) const
 360     {
 361       return ip1.myU[ myI ] < ip2.myU[ myI ];
 362     }
 363   };
 364   typedef boost::container::flat_set< IntPoint2D, IntPoint2DCompare >  TIntPointSet;
 365   typedef boost::container::flat_set< IntPoint2D*, IntPoint2DCompare > TIntPointPtrSet;
 366
 367   //--------------------------------------------------------------------------------
 368   /*!
 369    * \brief Face used to find translated position of the node
 370    */
 371   struct Face
 372   {
 373     const SMDS_MeshElement* myFace;
 374     SMESH_TNodeXYZ          myNode1; //!< nodes neighboring another node of myFace
 375     SMESH_TNodeXYZ          myNode2;
 376     const gp_XYZ*           myNorm;
 377     bool                    myNodeRightOrder;
 378     void operator=(const SMDS_MeshElement* f) { myFace = f; }
 379     const SMDS_MeshElement* operator->() { return myFace; }
 380     void SetNodes( int i0, int i1 ) //!< set myNode's
 381     {
 382       myNode1.Set( myFace->GetNode( i1 ));
 383       int i2 = ( i0 - 1 + myFace->NbCornerNodes() ) % myFace->NbCornerNodes();
 384       if ( i2 == i1 )
 385         i2 = ( i0 + 1 ) % myFace->NbCornerNodes();
 386       myNode2.Set( myFace->GetNode( i2 ));
 387       myNodeRightOrder = ( Abs( i2-i1 ) == 1 ) ?  i2 > i1  :  i2 < i1;
 388     }
 389     void SetOldNodes( const SMDS_Mesh& theSrcMesh )
 390     {
 391       myNode1.Set( theSrcMesh.FindNode( myNode1->GetID() ));
 392       myNode2.Set( theSrcMesh.FindNode( myNode2->GetID() ));
 393     }
 394     bool SetNormal( const std::vector< gp_XYZ >& faceNormals )
 395     {
 396       myNorm = & faceNormals[ myFace->GetID() ];
 397       return ( myNorm->SquareModulus() > gp::Resolution() * gp::Resolution() );
 398     }
 399     const gp_XYZ& Norm() const { return *myNorm; }
 400   };
 401
 402   //--------------------------------------------------------------------------------
 403   /*!
 404    * \brief Offset plane used to find translated position of the node
 405    */
 406   struct OffsetPlane
 407   {
 408     gp_XYZ myNode;
 409     Face*  myFace;
 410     gp_Pln myPln;
 411     gp_Lin myLines[2]; //!< line of intersection with neighbor OffsetPlane's
 412     bool   myIsLineOk[2];
 413     double myWeight[2];
 414
 415     void   Init( const gp_XYZ& node, Face& tria, double offset )
 416     {
 417       myNode = node;
 418       myFace = & tria;
 419       myPln  = gp_Pln( node + tria.Norm() * offset, tria.Norm() );
 420       myIsLineOk[0] = myIsLineOk[1] = false;
 421       myWeight[0] = myWeight[1] = 0;
 422     }
 423     bool   ComputeIntersectionLine( OffsetPlane& pln );
 424     void   SetSkewLine( const gp_Lin& line );
 425     gp_XYZ GetCommonPoint( int & nbOkPoints, double& sumWeight );
 426     gp_XYZ ProjectNodeOnLine( int & nbOkPoints );
 427     double Weight() const { return myWeight[0] + myWeight[1]; }
 428   };
 429
 430   //================================================================================
 431   /*!
 432    * \brief Set the second line
 433    */
 434   //================================================================================
 435
 436   void OffsetPlane::SetSkewLine( const gp_Lin& line )
 437   {
 438     myLines[1] = line;
 439     gp_XYZ n = myLines[0].Direction().XYZ() ^ myLines[1].Direction().XYZ();
 440     if (( myIsLineOk[1] = n.SquareModulus() > gp::Resolution() ))
 441       myPln = gp_Pln( myPln.Location(), n );
 442   }
 443
 444   //================================================================================
 445   /*!
 446    * \brief Project myNode on myLine[0]
 447    */
 448   //================================================================================
 449
 450   gp_XYZ OffsetPlane::ProjectNodeOnLine( int & nbOkPoints )
 451   {
 452     gp_XYZ p = gp::Origin().XYZ();
 453     if ( myIsLineOk[0] )
 454     {
 455       gp_Vec l2n( myLines[0].Location(), myNode );
 456       double u = l2n * myLines[0].Direction();
 457       p = myLines[0].Location().XYZ() + u * myLines[0].Direction().XYZ();
 458       ++nbOkPoints;
 459     }
 460     return p;
 461   }
 462
 463   //================================================================================
 464   /*!
 465    * \brief Computes intersection point of myLines
 466    */
 467   //================================================================================
 468
 469   gp_XYZ OffsetPlane::GetCommonPoint( int & nbOkPoints, double& sumWeight )
 470   {
 471     if ( !myIsLineOk[0] || !myIsLineOk[1] )
 472     {
 473       // sumWeight += myWeight[0];
 474       // return ProjectNodeOnLine( nbOkPoints ) * myWeight[0];
 475       return gp::Origin().XYZ();
 476     }
 477
 478     gp_XYZ p;
 479
 480     gp_Vec lPerp0 = myLines[0].Direction().XYZ() ^ myPln.Axis().Direction().XYZ();
 481     double  dot01 = lPerp0 * myLines[1].Direction().XYZ();
 482     if ( Abs( dot01 ) > 0.05 )
 483     {
 484       gp_Vec l0l1 = myLines[1].Location().XYZ() - myLines[0].Location().XYZ();
 485       double   u1 = - ( lPerp0 * l0l1 ) / dot01;
 486       p = ( myLines[1].Location().XYZ() + myLines[1].Direction().XYZ() * u1 );
 487     }
 488     else
 489     {
 490       gp_Vec  lv0( myLines[0].Location(), myNode),  lv1(myLines[1].Location(), myNode );
 491       double dot0( lv0 * myLines[0].Direction() ), dot1( lv1 * myLines[1].Direction() );
 492       p  = 0.5 * ( myLines[0].Location().XYZ() + myLines[0].Direction().XYZ() * dot0 );
 493       p += 0.5 * ( myLines[1].Location().XYZ() + myLines[1].Direction().XYZ() * dot1 );
 494     }
 495
 496     sumWeight += Weight();
 497     ++nbOkPoints;
 498
 499     return p * Weight();
 500   }
 501
 502   //================================================================================
 503   /*!
 504    * \brief Compute line of intersection of 2 planes
 505    */
 506   //================================================================================
 507
 508   bool OffsetPlane::ComputeIntersectionLine( OffsetPlane& theNextPln )
 509   {
 510     const gp_XYZ& n1 = myFace->Norm();
 511     const gp_XYZ& n2 = theNextPln.myFace->Norm();
 512
 513     gp_XYZ lineDir = n1 ^ n2;
 514     gp_Pnt linePos;
 515
 516     double x = Abs( lineDir.X() );
 517     double y = Abs( lineDir.Y() );
 518     double z = Abs( lineDir.Z() );
 519
 520     int cooMax; // max coordinate
 521     if (x > y) {
 522       if (x > z) cooMax = 1;
 523       else       cooMax = 3;
 524     }
 525     else {
 526       if (y > z) cooMax = 2;
 527       else       cooMax = 3;
 528     }
 529
 530     bool ok = true;
 531     if ( Abs( lineDir.Coord( cooMax )) < 0.05 )
 532     {
 533       // parallel planes - intersection is an offset of the common edge
 534       linePos  = 0.5 * ( myPln.Location().XYZ() + theNextPln.myPln.Location().XYZ() );
 535       lineDir  = myNode - myFace->myNode2;
 536       ok       = false;
 537       myWeight[0] = 0;
 538     }
 539     else
 540     {
 541       // the constants in the 2 plane equations
 542       double d1 = - ( n1 * myPln.Location().XYZ() );
 543       double d2 = - ( n2 * theNextPln.myPln.Location().XYZ() );
 544
 545       switch ( cooMax ) {
 546       case 1:
 547         linePos.SetX(  0 );
 548         linePos.SetY(( d2*n1.Z() - d1*n2.Z()) / lineDir.X() );
 549         linePos.SetZ(( d1*n2.Y() - d2*n1.Y()) / lineDir.X() );
 550         break;
 551       case 2:
 552         linePos.SetX(( d1*n2.Z() - d2*n1.Z()) / lineDir.Y() );
 553         linePos.SetY(  0 );
 554         linePos.SetZ(( d2*n1.X() - d1*n2.X()) / lineDir.Y() );
 555         break;
 556       case 3:
 557         linePos.SetX(( d2*n1.Y() - d1*n2.Y()) / lineDir.Z() );
 558         linePos.SetY(( d1*n2.X() - d2*n1.X()) / lineDir.Z() );
 559         linePos.SetZ(  0 );
 560       }
 561       myWeight[0] = lineDir.SquareModulus();
 562       if ( n1 * n2 < 0 )
 563         myWeight[0] = 2. - myWeight[0];
 564     }
 565     myLines   [ 0 ].SetDirection( lineDir );
 566     myLines   [ 0 ].SetLocation ( linePos );
 567     myIsLineOk[ 0 ] = ok;
 568
 569     theNextPln.myLines   [ 1 ] = myLines[ 0 ];
 570     theNextPln.myIsLineOk[ 1 ] = ok;
 571     theNextPln.myWeight  [ 1 ] = myWeight[ 0 ];
 572
 573     return ok;
 574   }
 575
 576   //================================================================================
 577   /*!
 578    * \brief Return a translated position of a node
 579    *  \param [in] new2OldNodes - new and old nodes
 580    *  \param [in] faceNormals - normals to input faces
 581    *  \param [in] theSrcMesh - initial mesh
 582    *  \param [in] theNewPos - a computed normal
 583    *  \return bool - true if theNewPos is computed
 584    */
 585   //================================================================================
 586
 587   bool getTranslatedPosition( const SMDS_MeshNode*         theNewNode,
 588                               const double                 theOffset,
 589                               const double                 /*theTol*/,
 590                               const double                 theSign,
 591                               const std::vector< gp_XYZ >& theFaceNormals,
 592                               SMDS_Mesh&                   theSrcMesh,
 593                               gp_XYZ&                      theNewPos)
 594   {
 595     bool useOneNormal = true;
 596
 597     // check if theNewNode needs an average position, i.e. theNewNode is convex
 598     // SMDS_ElemIteratorPtr faceIt = theNewNode->GetInverseElementIterator();
 599     // const SMDS_MeshElement*  f0 = faceIt->next();
 600     // const gp_XYZ&         norm0 = theFaceNormals[ f0->GetID() ];
 601     // const SMESH_NodeXYZ nodePos = theNewNode;
 602     // while ( faceIt->more() )
 603     // {
 604     //   const SMDS_MeshElement* f = faceIt->next();
 605     //   const int         nodeInd = f->GetNodeIndex( theNewNode );
 606     //   SMESH_NodeXYZ    nodePos2 = f->GetWrapNode( nodeInd + 1 );
 607     //   try {
 608     //     const gp_XYZ      nnDir = ( nodePos2 - nodePos ).Normalized();
 609     //   }
 610     //   catch {
 611     //     continue;
 612     //   }
 613     //   const double dot = norm0 * nnDir;
 614     //   bool isConvex =
 615
 616
 617
 618     // get faces surrounding theNewNode and sort them
 619     Face faces[ theMaxNbFaces ];
 620     SMDS_ElemIteratorPtr faceIt = theNewNode->GetInverseElementIterator();
 621     faces[0] = faceIt->next();
 622     while ( !faces[0].SetNormal( theFaceNormals ) && faceIt->more() )
 623       faces[0] = faceIt->next();
 624     int i0 = faces[0]->GetNodeIndex( theNewNode );
 625     int i1 = ( i0 + 1 ) % faces[0]->NbCornerNodes();
 626     faces[0].SetNodes( i0, i1 );
 627     TIDSortedElemSet elemSet, avoidSet;
 628     int iFace = 0;
 629     const SMDS_MeshElement* f;
 630     for ( ; faceIt->more() && iFace < theMaxNbFaces; faceIt->next() )
 631     {
 632       avoidSet.insert( faces[ iFace ].myFace );
 633       f = SMESH_MeshAlgos::FindFaceInSet( theNewNode, faces[ iFace ].myNode2.Node(),
 634                                           elemSet, avoidSet, &i0, &i1 );
 635       if ( !f )
 636       {
 637         std::reverse( &faces[0], &faces[0] + iFace + 1 );
 638         for ( int i = 0; i <= iFace; ++i )
 639         {
 640           std::swap( faces[i].myNode1, faces[i].myNode2 );
 641           faces[i].myNodeRightOrder = !faces[i].myNodeRightOrder;
 642         }
 643         f = SMESH_MeshAlgos::FindFaceInSet( theNewNode, faces[ iFace ].myNode2.Node(),
 644                                             elemSet, avoidSet, &i0, &i1 );
 645         if ( !f )
 646           break;
 647       }
 648       faces[ ++iFace ] = f;
 649       faces[ iFace ].SetNodes( i0, i1 );
 650       faces[ iFace ].SetNormal( theFaceNormals );
 651     }
 652     int nbFaces = iFace + 1;
 653
 654     theNewPos.SetCoord( 0, 0, 0 );
 655     gp_XYZ oldXYZ = SMESH_NodeXYZ( theNewNode );
 656
 657     // check if all faces are co-planar
 658     bool isPlanar = true;
 659     const double tol = 1e-2;
 660     for ( int i = 1; i < nbFaces &&  isPlanar;  ++i )
 661       isPlanar = ( faces[i].Norm() - faces[i-1].Norm() ).SquareModulus() < tol*tol;
 662
 663     if ( isPlanar )
 664     {
 665       theNewPos = oldXYZ + faces[0].Norm() * theOffset;
 666       return useOneNormal;
 667     }
 668
 669     // prepare OffsetPlane's
 670     OffsetPlane pln[ theMaxNbFaces ];
 671     for ( int i = 0; i < nbFaces; ++i )
 672     {
 673       faces[i].SetOldNodes( theSrcMesh );
 674       pln[i].Init( oldXYZ, faces[i], theOffset );
 675     }
 676     // intersect neighboring OffsetPlane's
 677     int nbOkPoints = 0;
 678     for ( int i = 1; i < nbFaces; ++i )
 679       nbOkPoints += pln[ i-1 ].ComputeIntersectionLine( pln[ i ]);
 680     nbOkPoints += pln[ nbFaces-1 ].ComputeIntersectionLine( pln[ 0 ]);
 681
 682     // move intersection lines to over parallel planes
 683     if ( nbOkPoints > 1 )
 684       for ( int i = 0; i < nbFaces; ++i )
 685         if ( pln[i].myIsLineOk[0] && !pln[i].myIsLineOk[1] )
 686           for ( int j = 1; j < nbFaces &&  !pln[i].myIsLineOk[1]; ++j )
 687           {
 688             int i2 = ( i + j ) % nbFaces;
 689             if ( pln[i2].myIsLineOk[0] )
 690               pln[i].SetSkewLine( pln[i2].myLines[0] );
 691           }
 692
 693     // get the translated position
 694     nbOkPoints = 0;
 695     double sumWegith = 0;
 696     const double minWeight = Sin( 30 * M_PI / 180. ) * Sin( 30 * M_PI / 180. );
 697     for ( int i = 0; i < nbFaces; ++i )
 698       if ( pln[ i ].Weight() > minWeight )
 699         theNewPos += pln[ i ].GetCommonPoint( nbOkPoints, sumWegith );
 700
 701     if ( nbOkPoints == 0 )
 702     {
 703       // there is only one feature edge;
 704       // find the theNewPos by projecting oldXYZ to any intersection line
 705       for ( int i = 0; i < nbFaces; ++i )
 706         theNewPos += pln[ i ].ProjectNodeOnLine( nbOkPoints );
 707
 708       if ( nbOkPoints == 0 )
 709       {
 710         theNewPos = oldXYZ + faces[0].Norm() * theOffset;
 711         return useOneNormal;
 712       }
 713       sumWegith = nbOkPoints;
 714     }
 715     theNewPos /= sumWegith;
 716
 717
 718     // mark theNewNode if it is concave
 719     useOneNormal = false;
 720     gp_Vec moveVec( oldXYZ, theNewPos );
 721     for ( int i = 0, iPrev = nbFaces-1; i < nbFaces; iPrev = i++ )
 722     {
 723       gp_Vec nodeVec( oldXYZ, faces[ i ].myNode1 );
 724       double u = ( moveVec * nodeVec ) / nodeVec.SquareMagnitude();
 725       if ( u > 0.5 ) // param [0,1] on nodeVec
 726       {
 727         theNewNode->setIsMarked( true );
 728       }
 729       if ( !useOneNormal )
 730       {
 731         gp_XYZ inFaceVec = faces[ i ].Norm() ^ nodeVec.XYZ();
 732         double       dot = inFaceVec * faces[ iPrev ].Norm();
 733         if ( !faces[ i ].myNodeRightOrder )
 734           dot *= -1;
 735         if ( dot * theSign < 0 )
 736         {
 737           gp_XYZ p1 = oldXYZ + faces[ i ].Norm()     * theOffset;
 738           gp_XYZ p2 = oldXYZ + faces[ iPrev ].Norm() * theOffset;
 739           useOneNormal = ( p1 - p2 ).SquareModulus() > 1e-12;
 740         }
 741       }
 742       if ( useOneNormal && theNewNode->isMarked() )
 743         break;
 744     }
 745
 746     return useOneNormal;
 747   }
 748
 749   //================================================================================
 750   /*!
 751    * \brief Remove small faces
 752    */
 753   //================================================================================
 754
 755   void removeSmallFaces( SMDS_Mesh*                        theMesh,
 756                          SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces,
 757                          const double                      theTol2 )
 758   {
 759     std::vector< SMESH_NodeXYZ > points(3);
 760     std::vector< const SMDS_MeshNode* > nodes(3);
 761     for ( SMDS_ElemIteratorPtr faceIt = theMesh->elementsIterator(); faceIt->more(); )
 762     {
 763       const SMDS_MeshElement* face = faceIt->next();
 764       points.assign( face->begin_nodes(), face->end_nodes() );
 765
 766       SMESH_NodeXYZ* prevN = & points.back();
 767       for ( size_t i = 0; i < points.size(); ++i )
 768       {
 769         double dist2 = ( *prevN - points[ i ]).SquareModulus();
 770         if ( dist2 < theTol2 )
 771         {
 772           const SMDS_MeshNode* nToRemove =
 773             (*prevN)->GetID() > points[ i ]->GetID() ? prevN->Node() : points[ i ].Node();
 774           const SMDS_MeshNode* nToKeep =
 775             nToRemove == points[ i ].Node() ? prevN->Node() : points[ i ].Node();
 776           for ( SMDS_ElemIteratorPtr fIt = nToRemove->GetInverseElementIterator(); fIt->more(); )
 777           {
 778             const SMDS_MeshElement* f = fIt->next();
 779             if ( f == face )
 780               continue;
 781             nodes.assign( f->begin_nodes(), f->end_nodes() );
 782             nodes[ f->GetNodeIndex( nToRemove )] = nToKeep;
 783             theMesh->ChangeElementNodes( f, &nodes[0], nodes.size() );
 784           }
 785           theNew2OldFaces[ face->GetID() ].first = 0;
 786           theMesh->RemoveFreeElement( face );
 787           break;
 788         }
 789         prevN = & points[ i ];
 790       }
 791       continue;
 792     }
 793     return;
 794   }
 795
 796 } // namespace
 797
 798 namespace SMESH_MeshAlgos
 799 {
 800   //--------------------------------------------------------------------------------
 801   /*!
 802    * \brief Intersect faces of a mesh
 803    */
 804   struct Intersector::Algo
 805   {
 806     SMDS_Mesh*                   myMesh;
 807     double                       myTol, myEps;
 808     const std::vector< gp_XYZ >& myNormals;
 809     TCutLinkMap                  myCutLinks; //!< assure sharing of new nodes
 810     TCutFaceMap                  myCutFaces;
 811     TNNMap                       myRemove2KeepNodes; //!< node merge map
 812
 813     // data to intersect 2 faces
 814     const SMDS_MeshElement*      myFace1;
 815     const SMDS_MeshElement*      myFace2;
 816     std::vector< SMESH_NodeXYZ > myNodes1, myNodes2;
 817     std::vector< double >        myDist1,  myDist2;
 818     int                          myInd1, myInd2; // coordinate indices on an axis-aligned plane
 819     int                          myNbOnPlane1, myNbOnPlane2;
 820     TIntPointSet                 myIntPointSet;
 821
 822     Algo( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals )
 823       : myMesh( mesh ),
 824         myTol( tol ),
 825         myEps( 1e-100 ),
 826         //myEps( Sqrt( std::numeric_limits<double>::min() )),
 827         //myEps( gp::Resolution() ),
 828         myNormals( normals )
 829     {}
 830     void Cut( const SMDS_MeshElement* face1,
 831               const SMDS_MeshElement* face2,
 832               const int               nbCommonNodes );
 833     void Cut( const SMDS_MeshElement* face,
 834               SMESH_NodeXYZ&          lineEnd1,
 835               int                     edgeIndex1,
 836               SMESH_NodeXYZ&          lineEnd2,
 837               int                     edgeIndex2 );
 838     void MakeNewFaces( TElemIntPairVec& theNew2OldFaces,
 839                        TNodeIntPairVec& theNew2OldNodes,
 840                        const double     theSign,
 841                        const bool       theOptimize );
 842
 843     void IntersectNewEdges( const CutFace& theCFace );
 844
 845   private:
 846
 847     bool isPlaneIntersected( const gp_XYZ&                       n2,
 848                              const double                        d2,
 849                              const std::vector< SMESH_NodeXYZ >& nodes1,
 850                              std::vector< double > &             dist1,
 851                              int &                               nbOnPlane1,
 852                              int &                               iNotOnPlane1);
 853     void computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes,
 854                            const std::vector< double >&        dist,
 855                            const int                           nbOnPln,
 856                            const int                           iMaxCoo,
 857                            double *                            u,
 858                            int*                                iE);
 859     void cutCoplanar();
 860     void addLink ( CutLink& link );
 861     bool findLink( CutLink& link );
 862     bool coincide( const gp_XYZ& p1, const gp_XYZ& p2, const double tol ) const
 863     {
 864       return ( p1 - p2 ).SquareModulus() < tol * tol;
 865     }
 866     gp_XY p2D( const gp_XYZ& p ) const { return gp_XY( p.Coord( myInd1 ), p.Coord( myInd2 )); }
 867
 868     void intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1,
 869                         const std::vector< double > &       dist1,
 870                         const int                           iEdge1,
 871                         const SMDS_MeshElement*             face2,
 872                         CutLink&                            link1);
 873     void findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes,
 874                               const std::vector< double > &       dist,
 875                               CutLink&                            link );
 876     void replaceIntNode( const SMDS_MeshNode* nToKeep, const SMDS_MeshNode* nToRemove );
 877     void computeIntPoint( const double           u1,
 878                           const double           u2,
 879                           const int              iE1,
 880                           const int              iE2,
 881                           CutLink &              link,
 882                           const SMDS_MeshNode* & node1,
 883                           const SMDS_MeshNode* & node2);
 884     void cutCollinearLink( const int                           iNotOnPlane1,
 885                            const std::vector< SMESH_NodeXYZ >& nodes1,
 886                            const SMDS_MeshElement*             face2,
 887                            const CutLink&                      link1,
 888                            const CutLink&                      link2);
 889     void setPlaneIndices( const gp_XYZ& planeNorm );
 890     bool intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1,
 891                             const gp_XY s2p0, const gp_XY s2p1,
 892                             double &    t1,   double &    t2,
 893                             bool &      isCollinear  );
 894     bool intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint );
 895     bool isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes );
 896     const SMDS_MeshNode* createNode( const gp_XYZ& p );
 897   };
 898
 899   //================================================================================
 900   /*!
 901    * \brief Return coordinate index with maximal abs value
 902    */
 903   //================================================================================
 904
 905   int MaxIndex( const gp_XYZ& x )
 906   {
 907     int iMaxCoo = ( Abs( x.X()) < Abs( x.Y() )) + 1;
 908     if ( Abs( x.Coord( iMaxCoo )) < Abs( x.Z() ))
 909       iMaxCoo = 3;
 910     return iMaxCoo;
 911   }
 912   //================================================================================
 913   /*!
 914    * \brief Store a CutLink
 915    */
 916   //================================================================================
 917
 918   const SMDS_MeshNode* Intersector::Algo::createNode( const gp_XYZ& p )
 919   {
 920     const SMDS_MeshNode* n = myMesh->AddNode( p.X(), p.Y(), p.Z() );
 921     n->setIsMarked( true ); // cut nodes are marked
 922     return n;
 923   }
 924
 925   //================================================================================
 926   /*!
 927    * \brief Store a CutLink
 928    */
 929   //================================================================================
 930
 931   void Intersector::Algo::addLink( CutLink& link )
 932   {
 933     link.myIndex = 0;
 934     const CutLink* added = & myCutLinks.Added( link );
 935     while ( added->myIntNode.Node() != link.myIntNode.Node() )
 936     {
 937       if ( !added->myIntNode )
 938       {
 939         added->myIntNode = link.myIntNode;
 940         break;
 941       }
 942       else
 943       {
 944         link.myIndex++;
 945         added = & myCutLinks.Added( link );
 946       }
 947     }
 948     link.myIndex = 0;
 949   }
 950
 951   //================================================================================
 952   /*!
 953    * \brief Find a CutLink with an intersection point coincident with that of a given link
 954    */
 955   //================================================================================
 956
 957   bool Intersector::Algo::findLink( CutLink& link )
 958   {
 959     link.myIndex = 0;
 960     while ( myCutLinks.Contains( link ))
 961     {
 962       const CutLink* added = & myCutLinks.Added( link );
 963       if ( !!added->myIntNode && coincide( added->myIntNode, link.myIntNode, myTol ))
 964       {
 965         link.myIntNode = added->myIntNode;
 966         return true;
 967       }
 968       link.myIndex++;
 969     }
 970     return false;
 971   }
 972
 973   //================================================================================
 974   /*!
 975    * \brief Check if a triangle intersects the plane of another triangle
 976    *  \param [in] nodes1 - nodes of triangle 1
 977    *  \param [in] n2 - normal of triangle 2
 978    *  \param [in] d2 - a constant of the plane equation 2
 979    *  \param [out] dist1 - distance of nodes1 from the plane 2
 980    *  \param [out] nbOnPlane - number of nodes1 lying on the plane 2
 981    *  \return bool - true if the triangle intersects the plane 2
 982    */
 983   //================================================================================
 984
 985   bool Intersector::Algo::isPlaneIntersected( const gp_XYZ&                       n2,
 986                                               const double                        d2,
 987                                               const std::vector< SMESH_NodeXYZ >& nodes1,
 988                                               std::vector< double > &             dist1,
 989                                               int &                               nbOnPlane1,
 990                                               int &                               iNotOnPlane1)
 991   {
 992     iNotOnPlane1 = nbOnPlane1 = 0;
 993     dist1.resize( nodes1.size() );
 994     for ( size_t i = 0; i < nodes1.size(); ++i )
 995     {
 996       dist1[i] = n2 * nodes1[i] + d2;
 997       if ( Abs( dist1[i] ) < myTol )
 998       {
 999         ++nbOnPlane1;
1000         dist1[i] = 0.;
1001       }
1002       else
1003       {
1004         iNotOnPlane1 = i;
1005       }
1006     }
1007     if ( nbOnPlane1 == 0 )
1008       for ( size_t i = 0; i < nodes1.size(); ++i )
1009         if ( dist1[iNotOnPlane1] * dist1[i] < 0 )
1010           return true;
1011
1012     return nbOnPlane1;
1013   }
1014
1015   //================================================================================
1016   /*!
1017    * \brief Compute parameters on the plane intersection line of intersections
1018    *        of edges of a triangle
1019    *  \param [in] nodes - triangle nodes
1020    *  \param [in] dist - distance of triangle nodes from the plane of another triangle
1021    *  \param [in] nbOnPln - number of nodes lying on the plane of another triangle
1022    *  \param [in] iMaxCoo - index of coordinate of max component of the plane intersection line
1023    *  \param [out] u - two computed parameters on the plane intersection line
1024    *  \param [out] iE - indices of intersected edges
1025    */
1026   //================================================================================
1027
1028   void Intersector::Algo::computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes,
1029                                             const std::vector< double >&        dist,
1030                                             const int                           nbOnPln,
1031                                             const int                           iMaxCoo,
1032                                             double *                            u,
1033                                             int*                                iE)
1034   {
1035     if ( nbOnPln == 3 )
1036     {
1037       u[0] = u[1] = 1e+100;
1038       return;
1039     }
1040     int nb = 0;
1041     int i1 = 2, i2 = 0;
1042     if ( nbOnPln == 1 && ( dist[i1] == 0. || dist[i2] == 0 ))
1043     {
1044       int i = dist[i1] == 0 ? i1 : i2;
1045       u [ 1 ] = nodes[ i ].Coord( iMaxCoo );
1046       iE[ 1 ] = i;
1047       i1 = i2++;
1048     }
1049     for ( ; i2 < 3 && nb < 2;  i1 = i2++ )
1050     {
1051       double dd = dist[i1] - dist[i2];
1052       if ( dd != 0. && dist[i2] * dist[i1] <= 0. )
1053       {
1054         double x1 = nodes[i1].Coord( iMaxCoo );
1055         double x2 = nodes[i2].Coord( iMaxCoo );
1056         u [ nb ] = x1 + ( x2 - x1 ) * dist[i1] / dd;
1057         iE[ nb ] = i1;
1058         ++nb;
1059       }
1060     }
1061     if ( u[0] > u[1] )
1062     {
1063       std::swap( u [0], u [1] );
1064       std::swap( iE[0], iE[1] );
1065     }
1066   }
1067
1068   //================================================================================
1069   /*!
1070    * \brief Try to find an intersection node on a link collinear with the plane intersection line
1071    */
1072   //================================================================================
1073
1074   void Intersector::Algo::findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes,
1075                                                const std::vector< double > &       dist,
1076                                                CutLink&                            link )
1077   {
1078     int i1 = ( dist[0] == 0 ? 0 : 1 ), i2 = ( dist[2] == 0 ? 2 : 1 );
1079     CutLink link2 = link;
1080     link2.Set( nodes[i1].Node(), nodes[i2].Node(), 0 );
1081     if ( findLink( link2 ))
1082       link.myIntNode = link2.myIntNode;
1083   }
1084
1085   //================================================================================
1086   /*!
1087    * \brief Compute intersection point of a link1 with a face2
1088    */
1089   //================================================================================
1090
1091   void Intersector::Algo::intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1,
1092                                          const std::vector< double > &       dist1,
1093                                          const int                           iEdge1,
1094                                          const SMDS_MeshElement*             face2,
1095                                          CutLink&                            link1)
1096   {
1097     const int iEdge2 = ( iEdge1 + 1 ) % nodes1.size();
1098     const SMESH_NodeXYZ& p1 = nodes1[ iEdge1 ];
1099     const SMESH_NodeXYZ& p2 = nodes1[ iEdge2 ];
1100
1101     link1.Set( p1.Node(), p2.Node(), face2 );
1102     const CutLink* link = & myCutLinks.Added( link1 );
1103     if ( !link->IntNode() )
1104     {
1105       if      ( dist1[ iEdge1 ] == 0. ) link1.myIntNode = p1;
1106       else if ( dist1[ iEdge2 ] == 0. ) link1.myIntNode = p2;
1107       else
1108       {
1109         gp_XYZ p = p1 + ( p2 - p1 ) * dist1[ iEdge1 ] / ( dist1[ iEdge1 ] - dist1[ iEdge2 ]);
1110         (gp_XYZ&)link1.myIntNode = p;
1111       }
1112     }
1113     else
1114     {
1115       gp_XYZ p = p1 + ( p2 - p1 ) * dist1[ iEdge1 ] / ( dist1[ iEdge1 ] - dist1[ iEdge2 ]);
1116       while ( link->IntNode() )
1117       {
1118         if ( coincide( p, link->myIntNode, myTol ))
1119         {
1120           link1.myIntNode = link->myIntNode;
1121           break;
1122         }
1123         link1.myIndex++;
1124         link = & myCutLinks.Added( link1 );
1125       }
1126       if ( !link1.IntNode() )
1127       {
1128         if      ( dist1[ iEdge1 ] == 0. ) link1.myIntNode = p1;
1129         else if ( dist1[ iEdge2 ] == 0. ) link1.myIntNode = p2;
1130         else                     (gp_XYZ&)link1.myIntNode = p;
1131       }
1132     }
1133   }
1134
1135   //================================================================================
1136   /*!
1137    * \brief Store node replacement in myCutFaces
1138    */
1139   //================================================================================
1140
1141   void Intersector::Algo::replaceIntNode( const SMDS_MeshNode* nToKeep,
1142                                           const SMDS_MeshNode* nToRemove )
1143   {
1144     if ( nToKeep == nToRemove )
1145       return;
1146     if ( nToRemove->GetID() < nToKeep->GetID() ) // keep node with lower ID
1147       myRemove2KeepNodes.Bind( nToKeep, nToRemove );
1148     else
1149       myRemove2KeepNodes.Bind( nToRemove, nToKeep );
1150   }
1151
1152   //================================================================================
1153   /*!
1154    * \brief Compute intersection point on a link of either of faces by choosing
1155    *        a link whose parameter on the intersection line in maximal
1156    *  \param [in] u1 - parameter on the intersection line of link iE1 of myFace1
1157    *  \param [in] u2 - parameter on the intersection line of link iE2 of myFace2
1158    *  \param [in] iE1 - index of a link myFace1
1159    *  \param [in] iE2 - index of a link myFace2
1160    *  \param [out] link - CutLink storing the intersection point
1161    *  \param [out] node1 - a node of the 2nd link if two links intersect
1162    *  \param [out] node2 - a node of the 2nd link if two links intersect
1163    */
1164   //================================================================================
1165
1166   void Intersector::Algo::computeIntPoint( const double           u1,
1167                                            const double           u2,
1168                                            const int              iE1,
1169                                            const int              iE2,
1170                                            CutLink &              link,
1171                                            const SMDS_MeshNode* & node1,
1172                                            const SMDS_MeshNode* & node2)
1173   {
1174     if      ( u1 > u2 + myTol )
1175     {
1176       intersectLink( myNodes1, myDist1, iE1, myFace2, link );
1177       node1 = node2 = 0;
1178       if ( myNbOnPlane2 == 2 )
1179         findIntPointOnPlane( myNodes2, myDist2, link );
1180     }
1181     else if ( u2 > u1 + myTol )
1182     {
1183       intersectLink( myNodes2, myDist2, iE2, myFace1, link );
1184       node1 = node2 = 0;
1185       if ( myNbOnPlane1 == 2 )
1186         findIntPointOnPlane( myNodes1, myDist1, link );
1187     }
1188     else // edges of two faces intersect the line at the same point
1189     {
1190       CutLink link2;
1191       intersectLink( myNodes1, myDist1, iE1, myFace2, link );
1192       intersectLink( myNodes2, myDist2, iE2, myFace1, link2 );
1193       node1 = link2.Node1();
1194       node2 = link2.Node2();
1195
1196       if      ( !link.IntNode() && link2.IntNode() )
1197         link.myIntNode = link2.myIntNode;
1198
1199       else if ( !link.IntNode() && !link2.IntNode() )
1200         (gp_XYZ&)link.myIntNode = 0.5 * ( link.myIntNode + link2.myIntNode );
1201
1202       else if ( link.IntNode() && link2.IntNode() )
1203         replaceIntNode( link.IntNode(), link2.IntNode() );
1204     }
1205   }
1206
1207   //================================================================================
1208   /*!
1209    * \brief Add intersections to a link collinear with the intersection line
1210    */
1211   //================================================================================
1212
1213   void Intersector::Algo::cutCollinearLink( const int                           iNotOnPlane1,
1214                                             const std::vector< SMESH_NodeXYZ >& nodes1,
1215                                             const SMDS_MeshElement*             face2,
1216                                             const CutLink&                      link1,
1217                                             const CutLink&                      link2)
1218
1219   {
1220     int iN1 = ( iNotOnPlane1 + 1 ) % 3;
1221     int iN2 = ( iNotOnPlane1 + 2 ) % 3;
1222     CutLink link( nodes1[ iN1 ].Node(), nodes1[ iN2 ].Node(), face2 );
1223     if ( link1.myFace != face2 )
1224     {
1225       link.myIntNode = link1.myIntNode;
1226       addLink( link );
1227     }
1228     if ( link2.myFace != face2 )
1229     {
1230       link.myIntNode = link2.myIntNode;
1231       addLink( link );
1232     }
1233   }
1234
1235   //================================================================================
1236   /*!
1237    * \brief Choose indices on an axis-aligned plane
1238    */
1239   //================================================================================
1240
1241   void Intersector::Algo::setPlaneIndices( const gp_XYZ& planeNorm )
1242   {
1243     switch ( MaxIndex( planeNorm )) {
1244     case 1: myInd1 = 2; myInd2 = 3; break;
1245     case 2: myInd1 = 3; myInd2 = 1; break;
1246     case 3: myInd1 = 1; myInd2 = 2; break;
1247     }
1248   }
1249
1250   //================================================================================
1251   /*!
1252    * \brief Intersect two faces
1253    */
1254   //================================================================================
1255
1256   void Intersector::Algo::Cut( const SMDS_MeshElement* face1,
1257                                const SMDS_MeshElement* face2,
1258                                const int               nbCommonNodes)
1259   {
1260     myFace1 = face1;
1261     myFace2 = face2;
1262     myNodes1.assign( face1->begin_nodes(), face1->end_nodes() );
1263     myNodes2.assign( face2->begin_nodes(), face2->end_nodes() );
1264
1265     const gp_XYZ& n1 = myNormals[ face1->GetID() ];
1266     const gp_XYZ& n2 = myNormals[ face2->GetID() ];
1267
1268     // check if triangles intersect
1269     int iNotOnPlane1, iNotOnPlane2;
1270     const double d2 = -( n2 * myNodes2[0]);
1271     if ( !isPlaneIntersected( n2, d2, myNodes1, myDist1, myNbOnPlane1, iNotOnPlane1 ))
1272       return;
1273     const double d1 = -( n1 * myNodes1[0]);
1274     if ( !isPlaneIntersected( n1, d1, myNodes2, myDist2, myNbOnPlane2, iNotOnPlane2 ))
1275       return;
1276
1277     if ( myNbOnPlane1 == 3 || myNbOnPlane2 == 3 )// triangles are co-planar
1278     {
1279       setPlaneIndices( myNbOnPlane1 == 3 ? n2 : n1 ); // choose indices on an axis-aligned plane
1280       cutCoplanar();
1281     }
1282     else if ( nbCommonNodes < 2 ) // triangle planes intersect
1283     {
1284       gp_XYZ lineDir = n1 ^ n2; // intersection line
1285
1286       // check if intervals of intersections of triangles with lineDir overlap
1287
1288       double u1[2], u2 [2]; // parameters on lineDir of edge intersection points { minU, maxU }
1289       int   iE1[2], iE2[2]; // indices of edges
1290       int iMaxCoo = MaxIndex( lineDir );
1291       computeIntervals( myNodes1, myDist1, myNbOnPlane1, iMaxCoo, u1, iE1 );
1292       computeIntervals( myNodes2, myDist2, myNbOnPlane2, iMaxCoo, u2, iE2 );
1293       if ( u1[1] < u2[0] - myTol || u2[1] < u1[0] - myTol )
1294         return; // intervals do not overlap
1295
1296       // make intersection nodes
1297
1298       const SMDS_MeshNode *l1n1, *l1n2, *l2n1, *l2n2;
1299       CutLink link1; // intersection with smaller u on lineDir
1300       computeIntPoint( u1[0], u2[0], iE1[0], iE2[0], link1, l1n1, l1n2 );
1301       CutLink link2; // intersection with larger u on lineDir
1302       computeIntPoint( -u1[1], -u2[1], iE1[1], iE2[1], link2, l2n1, l2n2 );
1303
1304       const CutFace& cf1 = myCutFaces.Added( CutFace( face1 ));
1305       const CutFace& cf2 = myCutFaces.Added( CutFace( face2 ));
1306
1307       if ( coincide( link1.myIntNode, link2.myIntNode, myTol ))
1308       {
1309         // intersection is a point
1310         if ( link1.IntNode() && link2.IntNode() )
1311           replaceIntNode( link1.IntNode(), link2.IntNode() );
1312
1313         CutLink* link = link2.IntNode() ? &link2 : &link1;
1314         if ( !link->IntNode() )
1315         {
1316           gp_XYZ p = 0.5 * ( link1.myIntNode + link2.myIntNode );
1317           link->myIntNode.Set( createNode( p ));
1318         }
1319         if ( !link1.IntNode() ) link1.myIntNode = link2.myIntNode;
1320         if ( !link2.IntNode() ) link2.myIntNode = link1.myIntNode;
1321
1322         cf1.AddPoint( link1, link2, myTol );
1323         if ( l1n1 ) link1.Set( l1n1, l1n2, face2 );
1324         if ( l2n1 ) link2.Set( l2n1, l2n2, face2 );
1325         cf2.AddPoint( link1, link2, myTol );
1326       }
1327       else
1328       {
1329         // intersection is a line segment
1330         if ( !link1.IntNode() )
1331           link1.myIntNode.Set( createNode( link1.myIntNode ));
1332         if ( !link2.IntNode() )
1333           link2.myIntNode.Set( createNode( link2.myIntNode ));
1334
1335         cf1.AddEdge( link1, link2, face2, myNbOnPlane1, iNotOnPlane1 );
1336         if ( l1n1 ) link1.Set( l1n1, l1n2, face2 );
1337         if ( l2n1 ) link2.Set( l2n1, l2n2, face2 );
1338         cf2.AddEdge( link1, link2, face1, myNbOnPlane2, iNotOnPlane2 );
1339
1340         // add intersections to a link collinear with the intersection line
1341         if ( myNbOnPlane1 == 2 && ( link1.myFace != face2 || link2.myFace != face2 ))
1342           cutCollinearLink( iNotOnPlane1, myNodes1, face2, link1, link2 );
1343
1344         if ( myNbOnPlane2 == 2 && ( link1.myFace != face1 || link2.myFace != face1 ))
1345           cutCollinearLink( iNotOnPlane2, myNodes2, face1, link1, link2 );
1346       }
1347
1348       addLink( link1 );
1349       addLink( link2 );
1350
1351     } // non co-planar case
1352
1353     return;
1354   }
1355
1356   //================================================================================
1357   /*!
1358    * \brief Store a face cut by a line given by its ends
1359    *        accompanied by indices of intersected face edges.
1360    *        Edge index is <0 if a line end is inside the face.
1361    *  \param [in] face - a face to cut
1362    *  \param [inout] lineEnd1 - line end coordinates + optional node existing at this point
1363    *  \param [in] edgeIndex1 - index of face edge cut by lineEnd1
1364    *  \param [inout] lineEnd2 - line end coordinates + optional node existing at this point
1365    *  \param [in] edgeIndex2 - index of face edge cut by lineEnd2
1366    */
1367   //================================================================================
1368
1369   void Intersector::Algo::Cut( const SMDS_MeshElement* face,
1370                                SMESH_NodeXYZ&          lineEnd1,
1371                                int                     edgeIndex1,
1372                                SMESH_NodeXYZ&          lineEnd2,
1373                                int                     edgeIndex2 )
1374   {
1375     if ( lineEnd1.Node() && lineEnd2.Node() &&
1376          face->GetNodeIndex( lineEnd1.Node() ) >= 0 &&
1377          face->GetNodeIndex( lineEnd2.Node() ) >= 0 )
1378       return; // intersection at a face node or edge
1379
1380     if ((int) myNormals.size() <= face->GetID() )
1381       const_cast< std::vector< gp_XYZ >& >( myNormals ).resize( face->GetID() + 1 );
1382
1383     const CutFace& cf = myCutFaces.Added( CutFace( face ));
1384     cf.InitLinks();
1385
1386     // look for intersection nodes coincident with line ends
1387     CutLink links[2];
1388     for ( int is2nd = 0; is2nd < 2; ++is2nd )
1389     {
1390       SMESH_NodeXYZ& lineEnd = is2nd ? lineEnd2 : lineEnd1;
1391       int          edgeIndex = is2nd ? edgeIndex2 : edgeIndex1;
1392       CutLink &         link = links[ is2nd ];
1393
1394       link.myIntNode = lineEnd;
1395
1396       for ( size_t i = ( edgeIndex < 0 ? 3 : 0  ); i < cf.myLinks.size(); ++i )
1397         if ( coincide( lineEnd, SMESH_NodeXYZ( cf.myLinks[i].myNode1 ), myTol ))
1398         {
1399           link.myIntNode = cf.myLinks[i].myNode1;
1400           break;
1401         }
1402
1403       if ( edgeIndex >= 0 )
1404       {
1405         link.Set( face->GetNode    ( edgeIndex ),
1406                   face->GetNodeWrap( edgeIndex + 1 ),
1407                   /*cuttingFace=*/0);
1408         findLink( link );
1409       }
1410
1411       if ( !link.myIntNode )
1412         link.myIntNode.Set( createNode( lineEnd ));
1413
1414       lineEnd._node = link.IntNode();
1415
1416       if ( link.myNode[0] )
1417         addLink( link );
1418     }
1419
1420     cf.AddEdge( links[0], links[1], /*face=*/0, /*nbOnPlane=*/0, /*iNotOnPlane=*/-1 );
1421   }
1422
1423   //================================================================================
1424   /*!
1425    * \brief Intersect two 2D line segments
1426    */
1427   //================================================================================
1428
1429   bool Intersector::Algo::intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1,
1430                                              const gp_XY s2p0, const gp_XY s2p1,
1431                                              double &    t1,   double &    t2,
1432                                              bool &      isCollinear )
1433   {
1434     gp_XY u = s1p1 - s1p0;
1435     gp_XY v = s2p1 - s2p0;
1436     gp_XY w = s1p0 - s2p0;
1437     double perpDotUV = u * gp_XY( -v.Y(), v.X() );
1438     double perpDotVW = v * gp_XY( -w.Y(), w.X() );
1439     double perpDotUW = u * gp_XY( -w.Y(), w.X() );
1440     double        u2 = u.SquareModulus();
1441     double        v2 = v.SquareModulus();
1442     if ( u2 < myEps * myEps || v2 < myEps * myEps )
1443       return false;
1444     if ( perpDotUV * perpDotUV / u2 / v2 < 1e-6 ) // cos ^ 2
1445     {
1446       if ( !isCollinear )
1447         return false; // no need in collinear solution
1448       if ( perpDotUW * perpDotUW / u2 > myTol * myTol )
1449         return false; // parallel
1450
1451       // collinear
1452       gp_XY w2 = s1p1 - s2p0;
1453       if ( Abs( v.X()) + Abs( u.X()) > Abs( v.Y()) + Abs( u.Y())) {
1454         t1 = w.X() / v.X();  // params on segment 2
1455         t2 = w2.X() / v.X();
1456       }
1457       else {
1458         t1 = w.Y() / v.Y();
1459         t2 = w2.Y() / v.Y();
1460       }
1461       if ( Max( t1,t2 ) <= 0 || Min( t1,t2 ) >= 1 )
1462         return false; // no overlap
1463       return true;
1464     }
1465     isCollinear = false;
1466
1467     t1 = perpDotVW / perpDotUV; // param on segment 1
1468     if ( t1 < 0. || t1 > 1. )
1469       return false; // intersection not within the segment
1470
1471     t2 = perpDotUW / perpDotUV; // param on segment 2
1472     if ( t2 < 0. || t2 > 1. )
1473       return false; // intersection not within the segment
1474
1475     return true;
1476   }
1477
1478   //================================================================================
1479   /*!
1480    * \brief Intersect two edges of co-planar triangles
1481    *  \param [inout] iE1 - edge index of triangle 1
1482    *  \param [inout] iE2 - edge index of triangle 2
1483    *  \param [inout] intPoints - intersection points
1484    *  \param [inout] nbIntPoints - nb of found intersection points
1485    */
1486   //================================================================================
1487
1488   bool Intersector::Algo::intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint )
1489   {
1490     int i01 = iE1, i11 = ( iE1 + 1 ) % 3;
1491     int i02 = iE2, i12 = ( iE2 + 1 ) % 3;
1492     if (( !intPoint.myIsCollinear ) &&
1493         ( myNodes1[ i01 ] == myNodes2[ i02 ] ||
1494           myNodes1[ i01 ] == myNodes2[ i12 ] ||
1495           myNodes1[ i11 ] == myNodes2[ i02 ] ||
1496           myNodes1[ i11 ] == myNodes2[ i12 ] ))
1497       return false;
1498
1499     // segment 1
1500     gp_XY s1p0 = p2D( myNodes1[ i01 ]);
1501     gp_XY s1p1 = p2D( myNodes1[ i11 ]);
1502
1503     // segment 2
1504     gp_XY s2p0 = p2D( myNodes2[ i02 ]);
1505     gp_XY s2p1 = p2D( myNodes2[ i12 ]);
1506
1507     double t1, t2;
1508     if ( !intersectEdgeEdge( s1p0,s1p1, s2p0,s2p1, t1, t2, intPoint.myIsCollinear ))
1509       return false;
1510
1511     intPoint.myEdgeInd[0] = iE1;
1512     intPoint.myEdgeInd[1] = iE2;
1513     intPoint.myU[0] = t1;
1514     intPoint.myU[1] = t2;
1515     (gp_XYZ&)intPoint.myNode = myNodes1[i01] * ( 1 - t1 ) + myNodes1[i11] * t1;
1516
1517     if ( intPoint.myIsCollinear )
1518       return true;
1519
1520     // try to find existing node at intPoint.myNode
1521
1522     if ( myNodes1[ i01 ] == myNodes2[ i02 ] ||
1523          myNodes1[ i01 ] == myNodes2[ i12 ] ||
1524          myNodes1[ i11 ] == myNodes2[ i02 ] ||
1525          myNodes1[ i11 ] == myNodes2[ i12 ] )
1526       return false;
1527
1528     const double coincTol = myTol * 1e-3;
1529
1530     CutLink link1( myNodes1[i01].Node(), myNodes1[i11].Node(), myFace2 );
1531     CutLink link2( myNodes2[i02].Node(), myNodes2[i12].Node(), myFace1 );
1532
1533     SMESH_NodeXYZ& n1 = myNodes1[ t1 < 0.5 ? i01 : i11 ];
1534     bool same1 = coincide( n1, intPoint.myNode, coincTol );
1535     if ( same1 )
1536     {
1537       link2.myIntNode = intPoint.myNode = n1;
1538       addLink( link2 );
1539     }
1540     SMESH_NodeXYZ& n2 = myNodes2[ t2 < 0.5 ? i02 : i12 ];
1541     bool same2 = coincide( n2, intPoint.myNode, coincTol );
1542     if ( same2 )
1543     {
1544       link1.myIntNode = intPoint.myNode = n2;
1545       addLink( link1 );
1546       if ( same1 )
1547       {
1548         replaceIntNode( n1.Node(), n2.Node() );
1549         return false;
1550       }
1551       return true;
1552     }
1553     if ( same1 )
1554       return true;
1555
1556     link1.myIntNode = intPoint.myNode;
1557     if ( findLink( link1 ))
1558     {
1559       intPoint.myNode = link2.myIntNode = link1.myIntNode;
1560       addLink( link2 );
1561       return true;
1562     }
1563
1564     link2.myIntNode = intPoint.myNode;
1565     if ( findLink( link2 ))
1566     {
1567       intPoint.myNode = link1.myIntNode = link2.myIntNode;
1568       addLink( link1 );
1569       return true;
1570     }
1571
1572     for ( int is2nd = 0; is2nd < 2; ++is2nd )
1573     {
1574       const SMDS_MeshElement* f = is2nd ? myFace1 : myFace2;
1575       if ( !f ) continue;
1576       const CutFace& cf = myCutFaces.Added( CutFace( is2nd ? myFace2 : myFace1 ));
1577       for ( size_t i = 0; i < cf.myLinks.size(); ++i )
1578         if ( cf.myLinks[i].myFace == f &&
1579              //cf.myLinks[i].myIndex != EdgePart::_COPLANAR &&
1580              coincide( intPoint.myNode, SMESH_NodeXYZ( cf.myLinks[i].myNode1 ), coincTol ))
1581         {
1582           intPoint.myNode.Set( cf.myLinks[i].myNode1 );
1583           return true;
1584         }
1585     }
1586
1587     // make a new node
1588
1589     intPoint.myNode._node = createNode( intPoint.myNode );
1590     link1.myIntNode = link2.myIntNode = intPoint.myNode;
1591     addLink( link1 );
1592     addLink( link2 );
1593
1594     return true;
1595   }
1596
1597
1598   //================================================================================
1599   /*!
1600    * \brief Check if a point is contained in a triangle
1601    */
1602   //================================================================================
1603
1604   bool Intersector::Algo::isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes )
1605   {
1606     double bc1, bc2;
1607     SMESH_MeshAlgos::GetBarycentricCoords( p2D( p ),
1608                                            p2D( nodes[0] ), p2D( nodes[1] ), p2D( nodes[2] ),
1609                                            bc1, bc2 );
1610     //return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. );
1611     return ( myTol < bc1 && myTol < bc2 && bc1 + bc2 + myTol < 1. );
1612   }
1613
1614   //================================================================================
1615   /*!
1616    * \brief Intersect two co-planar faces
1617    */
1618   //================================================================================
1619
1620   void Intersector::Algo::cutCoplanar()
1621   {
1622     // find intersections of edges
1623
1624     IntPoint2D intPoints[ 6 ];
1625     int      nbIntPoints = 0;
1626     for ( int iE1 = 0; iE1 < 3; ++iE1 )
1627     {
1628       int maxNbIntPoints = nbIntPoints + 2;
1629       for ( int iE2 = 0; iE2 < 3 &&  nbIntPoints < maxNbIntPoints; ++iE2 )
1630         nbIntPoints += intersectEdgeEdge( iE1, iE2, intPoints[ nbIntPoints ]);
1631     }
1632     const int minNbOnPlane = Min( myNbOnPlane1, myNbOnPlane2 );
1633
1634     if ( nbIntPoints == 0 ) // no intersections of edges
1635     {
1636       bool is1in2;
1637       if      ( isPointInTriangle( myNodes1[0], myNodes2 )) // face2 includes face1
1638         is1in2 = true;
1639       else if ( isPointInTriangle( myNodes2[0], myNodes1 )) // face1 includes face2
1640         is1in2 = false;
1641       else
1642         return;
1643
1644       // add edges of an inner triangle to an outer one
1645
1646       const std::vector< SMESH_NodeXYZ >& nodesIn = is1in2 ? myNodes1 : myNodes2;
1647       const SMDS_MeshElement*             faceOut = is1in2 ? myFace2  : myFace1;
1648       const SMDS_MeshElement*              faceIn = is1in2 ? myFace1  : myFace2;
1649
1650       const CutFace& outFace = myCutFaces.Added( CutFace( faceOut ));
1651       CutLink link1( nodesIn.back().Node(), nodesIn.back().Node(), faceOut );
1652       CutLink link2( nodesIn.back().Node(), nodesIn.back().Node(), faceOut );
1653
1654       link1.myIntNode = nodesIn.back();
1655       for ( size_t i = 0; i < nodesIn.size(); ++i )
1656       {
1657         link2.myIntNode = nodesIn[ i ];
1658         outFace.AddEdge( link1, link2, faceIn, minNbOnPlane );
1659         link1.myIntNode = link2.myIntNode;
1660       }
1661     }
1662     else
1663     {
1664       // add parts of edges to a triangle including them
1665
1666       CutLink link1, link2;
1667       IntPoint2D ip0, ip1;
1668       ip0.myU[0] = ip0.myU[1] = 0.;
1669       ip1.myU[0] = ip1.myU[1] = 1.;
1670       ip0.myEdgeInd[0] = ip0.myEdgeInd[1] = ip1.myEdgeInd[0] = ip1.myEdgeInd[1] = 0;
1671
1672       for ( int isFromFace1 = 0; isFromFace1 < 2; ++isFromFace1 )
1673       {
1674         const SMDS_MeshElement*                faceTo = isFromFace1 ? myFace2  : myFace1;
1675         const SMDS_MeshElement*              faceFrom = isFromFace1 ? myFace1  : myFace2;
1676         const std::vector< SMESH_NodeXYZ >&   nodesTo = isFromFace1 ? myNodes2 : myNodes1;
1677         const std::vector< SMESH_NodeXYZ >& nodesFrom = isFromFace1 ? myNodes1 : myNodes2;
1678         const int                                 iTo = isFromFace1 ? 1 : 0;
1679         const int                               iFrom = isFromFace1 ? 0 : 1;
1680         //const int                       nbOnPlaneFrom = isFromFace1 ? myNbOnPlane1 : myNbOnPlane2;
1681
1682         const CutFace* cutFaceTo   = & myCutFaces.Added( CutFace( faceTo ));
1683         // const CutFace* cutFaceFrom = 0;
1684         // if ( nbOnPlaneFrom > minNbOnPlane )
1685         //   cutFaceFrom = & myCutFaces.Added( CutFace( faceTo ));
1686
1687         link1.myFace = link2.myFace = faceTo;
1688
1689         IntPoint2DCompare ipCompare( iFrom );
1690         TIntPointPtrSet pointsOnEdge( ipCompare ); // IntPoint2D sorted by parameter on edge
1691
1692         for ( size_t iE = 0; iE < nodesFrom.size(); ++iE )
1693         {
1694           // get parts of an edge iE
1695
1696           ip0.myEdgeInd[ iTo ] = iE;
1697           ip1.myEdgeInd[ iTo ] = ( iE + 1 ) % nodesFrom.size();
1698           ip0.myNode = nodesFrom[ ip0.myEdgeInd[ iTo ]];
1699           ip1.myNode = nodesFrom[ ip1.myEdgeInd[ iTo ]];
1700
1701           pointsOnEdge.clear();
1702
1703           for ( int iP = 0; iP < nbIntPoints; ++iP )
1704             if ( intPoints[ iP ].myEdgeInd[ iFrom ] == iE )
1705               pointsOnEdge.insert( & intPoints[ iP ] );
1706
1707           pointsOnEdge.insert( pointsOnEdge.begin(), & ip0 );
1708           pointsOnEdge.insert( pointsOnEdge.end(),   & ip1 );
1709
1710           // add edge parts to faceTo
1711
1712           TIntPointPtrSet::iterator ipIt = pointsOnEdge.begin() + 1;
1713           for ( ; ipIt != pointsOnEdge.end(); ++ipIt )
1714           {
1715             const IntPoint2D* p1 = *(ipIt-1);
1716             const IntPoint2D* p2 = *ipIt;
1717             gp_XYZ middle = 0.5 * ( p1->myNode + p2->myNode );
1718             if ( isPointInTriangle( middle, nodesTo ))
1719             {
1720               p1->InitLink( link1, iTo, ( p1 != & ip0 ) ? nodesTo : nodesFrom );
1721               p2->InitLink( link2, iTo, ( p2 != & ip1 ) ? nodesTo : nodesFrom );
1722               cutFaceTo->AddEdge( link1, link2, faceFrom, minNbOnPlane );
1723
1724               // if ( cutFaceFrom )
1725               // {
1726               //   p1->InitLink( link1, iFrom, nodesFrom );
1727               //   p2->InitLink( link2, iFrom, nodesFrom );
1728               //   cutFaceTo->AddEdge( link1, link2, faceTo, minNbOnPlane );
1729               // }
1730             }
1731           }
1732         }
1733       }
1734     }
1735     return;
1736
1737   } // Intersector::Algo::cutCoplanar()
1738
1739   //================================================================================
1740   /*!
1741    * \brief Intersect edges added to myCutFaces
1742    */
1743   //================================================================================
1744
1745   void Intersector::Algo::IntersectNewEdges( const CutFace& cf )
1746   {
1747     IntPoint2D intPoint;
1748
1749     if ( cf.NbInternalEdges() < 2 )
1750       return;
1751
1752     if ( myNodes1.empty() )
1753     {
1754       myNodes1.resize(2);
1755       myNodes2.resize(2);
1756     }
1757
1758     const gp_XYZ& faceNorm = myNormals[ cf.myInitFace->GetID() ];
1759     setPlaneIndices( faceNorm ); // choose indices on an axis-aligned plane
1760
1761     size_t limit = cf.myLinks.size() * cf.myLinks.size() * 2;
1762
1763     size_t i1 = 3;
1764     while ( cf.myLinks[i1-1].IsInternal() && i1 > 0 )
1765       --i1;
1766
1767     for ( ; i1 < cf.myLinks.size(); ++i1 )
1768     {
1769       if ( !cf.myLinks[i1].IsInternal() )
1770         continue;
1771
1772       myIntPointSet.clear();
1773       for ( size_t i2 = i1 + 2; i2 < cf.myLinks.size(); ++i2 )
1774       {
1775         if ( !cf.myLinks[i2].IsInternal() )
1776           continue;
1777
1778         // prepare to intersection
1779         myFace1     = cf.myLinks[i1].myFace;
1780         myNodes1[0] = cf.myLinks[i1].myNode1;
1781         myNodes1[1] = cf.myLinks[i1].myNode2;
1782         myFace2     = cf.myLinks[i2].myFace;
1783         myNodes2[0] = cf.myLinks[i2].myNode1;
1784         myNodes2[1] = cf.myLinks[i2].myNode2;
1785
1786         // intersect
1787         intPoint.myIsCollinear = true; // to find collinear solutions
1788         if ( intersectEdgeEdge( 0, 0, intPoint ))
1789         {
1790           if ( cf.myLinks[i1].IsSame( cf.myLinks[i2] )) // remove i2
1791           {
1792             cf.myLinks[i1].ReplaceCoplanar( cf.myLinks[i2] );
1793             cf.myLinks.erase( cf.myLinks.begin() + i2, cf.myLinks.begin() + i2 + 2 );
1794             --i2;
1795             continue;
1796           }
1797           if ( !intPoint.myIsCollinear )
1798           {
1799             intPoint.myEdgeInd[1] = i2;
1800             myIntPointSet.insert( intPoint );
1801           }
1802           else // if ( intPoint.myIsCollinear ) // overlapping edges
1803           {
1804             myIntPointSet.clear(); // to recompute
1805
1806             if ( intPoint.myU[0] > intPoint.myU[1] ) // orient in same direction
1807             {
1808               std::swap( intPoint.myU[0], intPoint.myU[1] );
1809               std::swap( myNodes1[0], myNodes1[1] );
1810             }
1811             // replace _COPLANAR by _INTERNAL
1812             cf.myLinks[i1].ReplaceCoplanar( cf.myLinks[i1+1] );
1813             cf.myLinks[i2].ReplaceCoplanar( cf.myLinks[i2+1] );
1814
1815             if ( coincide( myNodes1[0], myNodes2[0], myTol ) &&
1816                  coincide( myNodes1[1], myNodes2[1], myTol ))
1817             {
1818               cf.myLinks.erase( cf.myLinks.begin() + i2, cf.myLinks.begin() + i2 + 2 );
1819               --i2;
1820               continue;
1821             }
1822
1823             EdgePart common = cf.myLinks[i1];
1824             common.ReplaceCoplanar( cf.myLinks[i2] );
1825
1826             const SMDS_MeshNode* n1 = myNodes1[0].Node(); // end nodes of an overlapping part
1827             const SMDS_MeshNode* n2 = myNodes1[1].Node();
1828             size_t i3 = cf.myLinks.size();
1829
1830             if ( myNodes1[0] != myNodes2[0] ) // a part before the overlapping one
1831             {
1832               if ( intPoint.myU[0] < 0 )
1833                 cf.myLinks[i1].Set( myNodes1[0].Node(), myNodes2[0].Node(),
1834                                     cf.myLinks[i1].myFace, cf.myLinks[i1].myIndex );
1835               else
1836                 cf.myLinks[i1].Set( myNodes2[0].Node(), myNodes1[0].Node(),
1837                                     cf.myLinks[i2].myFace, cf.myLinks[i2].myIndex );
1838
1839               cf.myLinks[i1+1].Set( cf.myLinks[i1].myNode2,
1840                                     cf.myLinks[i1].myNode1,
1841                                     cf.myLinks[i1].myFace,
1842                                     cf.myLinks[i1].myIndex);
1843               n1 = cf.myLinks[i1].myNode2;
1844             }
1845             else
1846               i3 = i1;
1847
1848             if ( myNodes1[1] != myNodes2[1] ) // a part after the overlapping one
1849             {
1850               if ( intPoint.myU[1] < 1 )
1851                 cf.myLinks[i2].Set( myNodes1[1].Node(), myNodes2[1].Node(),
1852                                     cf.myLinks[i2].myFace, cf.myLinks[i2].myIndex );
1853               else
1854                 cf.myLinks[i2].Set( myNodes2[1].Node(), myNodes1[1].Node(),
1855                                     cf.myLinks[i1].myFace, cf.myLinks[i1].myIndex );
1856
1857               cf.myLinks[i2+1].Set( cf.myLinks[i2].myNode2,
1858                                     cf.myLinks[i2].myNode1,
1859                                     cf.myLinks[i2].myFace,
1860                                     cf.myLinks[i2].myIndex);
1861               n2 = cf.myLinks[i2].myNode1;
1862             }
1863             else
1864               i3 = i2;
1865
1866             if ( i3 == cf.myLinks.size() )
1867               cf.myLinks.resize( i3 + 2 );
1868
1869             cf.myLinks[i3].Set  ( n1, n2, common.myFace, common.myIndex );
1870             cf.myLinks[i3+1].Set( n2, n1, common.myFace, common.myIndex );
1871
1872             i2 = i1 + 1; // recheck modified i1
1873             continue;
1874           }
1875           //else
1876           // {
1877           //   // remember a new node
1878           //   CutLink link1( myNodes1[0].Node(), myNodes1[1].Node(), cf.myInitFace );
1879           //   CutLink link2( myNodes2[0].Node(), myNodes2[1].Node(), cf.myInitFace );
1880           //   link2.myIntNode = link1.myIntNode = intPoint.myNode;
1881           //   addLink( link1 );
1882           //   addLink( link2 );
1883
1884           //   // split edges
1885           //   size_t i = cf.myLinks.size();
1886           //   if ( intPoint.myNode != cf.myLinks[ i1 ].myNode1 &&
1887           //        intPoint.myNode != cf.myLinks[ i1 ].myNode2 )
1888           //   {
1889           //     cf.myLinks.push_back( cf.myLinks[ i1 ]);
1890           //     cf.myLinks.push_back( cf.myLinks[ i1 + 1 ]);
1891           //     cf.myLinks[ i1 ].myNode2 = cf.myLinks[ i1 + 1 ].myNode1 = intPoint.Node();
1892           //     cf.myLinks[ i  ].myNode1 = cf.myLinks[ i  + 1 ].myNode2 = intPoint.Node();
1893           //   }
1894           //   if ( intPoint.myNode != cf.myLinks[ i2 ].myNode1 &&
1895           //        intPoint.myNode != cf.myLinks[ i2 ].myNode2 )
1896           //   {
1897           //     i = cf.myLinks.size();
1898           //     cf.myLinks.push_back( cf.myLinks[ i2 ]);
1899           //     cf.myLinks.push_back( cf.myLinks[ i2 + 1 ]);
1900           //     cf.myLinks[ i2 ].myNode2 = cf.myLinks[ i2 + 1 ].myNode1 = intPoint.Node();
1901           //     cf.myLinks[ i  ].myNode1 = cf.myLinks[ i  + 1 ].myNode2 = intPoint.Node();
1902           //   }
1903           // }
1904
1905         } // if ( intersectEdgeEdge( 0, 0, intPoint ))
1906
1907         ++i2;
1908         --limit;
1909       }
1910
1911       // split i1 edge and all edges it intersects
1912       // don't do it inside intersection loop in order not to loose direction of i1 edge
1913       if ( !myIntPointSet.empty() )
1914       {
1915         cf.myLinks.reserve( cf.myLinks.size() + myIntPointSet.size() * 2 + 2 );
1916
1917         EdgePart* edge1 = &cf.myLinks[ i1 ];
1918         EdgePart* twin1 = &cf.myLinks[ i1 + 1 ];
1919
1920         TIntPointSet::iterator ipIt = myIntPointSet.begin();
1921         for ( ; ipIt != myIntPointSet.end(); ++ipIt ) // int points sorted on i1 edge
1922         {
1923           size_t i = cf.myLinks.size();
1924           if ( ipIt->myNode != edge1->myNode1 &&
1925                ipIt->myNode != edge1->myNode2 )
1926           {
1927             cf.myLinks.push_back( *edge1 );
1928             cf.myLinks.push_back( *twin1 );
1929             edge1->myNode2          = twin1->myNode1              = ipIt->Node();
1930             cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = ipIt->Node();
1931             edge1 = & cf.myLinks[ i ];
1932             twin1 = & cf.myLinks[ i + 1 ];
1933           }
1934           size_t i2 = ipIt->myEdgeInd[1];
1935           if ( ipIt->myNode != cf.myLinks[ i2 ].myNode1 &&
1936                ipIt->myNode != cf.myLinks[ i2 ].myNode2 )
1937           {
1938             i = cf.myLinks.size();
1939             cf.myLinks.push_back( cf.myLinks[ i2 ]);
1940             cf.myLinks.push_back( cf.myLinks[ i2 + 1 ]);
1941             cf.myLinks[ i2 ].myNode2 = cf.myLinks[ i2 + 1 ].myNode1 = ipIt->Node();
1942             cf.myLinks[ i  ].myNode1 = cf.myLinks[ i  + 1 ].myNode2 = ipIt->Node();
1943           }
1944         }
1945         if ( cf.myLinks.size() >= limit )
1946           throw SALOME_Exception( "Infinite loop in Intersector::Algo::IntersectNewEdges()" );
1947       }
1948       ++i1; // each internal edge encounters twice
1949     }
1950     return;
1951   }
1952
1953   //================================================================================
1954   /*!
1955    * \brief Split intersected faces
1956    */
1957   //================================================================================
1958
1959   void Intersector::Algo::MakeNewFaces( SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces,
1960                                         SMESH_MeshAlgos::TNodeIntPairVec& theNew2OldNodes,
1961                                         const double                      theSign,
1962                                         const bool                        theOptimize)
1963   {
1964     // fill theNew2OldFaces if empty
1965     TCutFaceMap::const_iterator cutFacesIt = myCutFaces.cbegin();
1966     if ( theNew2OldFaces.empty() )
1967       for ( ; cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
1968       {
1969         const CutFace& cf = *cutFacesIt;
1970         smIdType index = cf.myInitFace->GetID(); // index in theNew2OldFaces
1971         if ((smIdType) theNew2OldFaces.size() <= index )
1972           theNew2OldFaces.resize( index + 1 );
1973         theNew2OldFaces[ index ] = std::make_pair( cf.myInitFace, index );
1974       }
1975
1976     // unmark all nodes except intersection ones
1977
1978     for ( SMDS_NodeIteratorPtr nIt = myMesh->nodesIterator(); nIt->more(); )
1979     {
1980       const SMDS_MeshNode* n = nIt->next();
1981       if ( n->isMarked() && n->GetID()-1 < (int) theNew2OldNodes.size() )
1982         n->setIsMarked( false );
1983     }
1984     // SMESH_MeshAlgos::MarkElems( myMesh->nodesIterator(), false );
1985
1986     TCutLinkMap::const_iterator cutLinksIt = myCutLinks.cbegin();
1987     // for ( ; cutLinksIt != myCutLinks.cend(); ++cutLinksIt )
1988     // {
1989     //   const CutLink& link = *cutLinksIt;
1990     //   if ( link.IntNode() && link.IntNode()->GetID()-1 < (int) theNew2OldNodes.size() )
1991     //     link.IntNode()->setIsMarked( true );
1992     // }
1993
1994     // intersect edges added to myCutFaces
1995
1996     for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
1997     {
1998       const CutFace& cf = *cutFacesIt;
1999       cf.ReplaceNodes( myRemove2KeepNodes );
2000       IntersectNewEdges( cf );
2001     }
2002
2003     // make new faces
2004
2005     EdgeLoopSet                            loopSet;
2006     SMESH_MeshAlgos::Triangulate           triangulator( theOptimize );
2007     std::vector< EdgePart >                cutOffLinks;
2008     TLinkMap                               cutOffCoplanarLinks;
2009     std::vector< const CutFace* >          touchedFaces;
2010     SMESH_MeshAlgos::TElemIntPairVec::value_type new2OldTria;
2011     CutFace                                cutFace(0);
2012     std::vector< const SMDS_MeshNode* >    nodes;
2013     std::vector<const SMDS_MeshElement *>  faces;
2014
2015     cutOffLinks.reserve( myCutFaces.Extent() * 2 );
2016
2017     for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
2018     {
2019       const CutFace& cf = *cutFacesIt;
2020       if ( !cf.IsCut() )
2021       {
2022         touchedFaces.push_back( & cf );
2023         continue;
2024       }
2025
2026       const gp_XYZ& normal = myNormals[ cf.myInitFace->GetID() ];
2027
2028       // form loops of new faces
2029       cf.ReplaceNodes( myRemove2KeepNodes );
2030       cf.MakeLoops( loopSet, normal );
2031
2032       // avoid loops that are not connected to boundary edges of cf.myInitFace
2033       if ( cf.RemoveInternalLoops( loopSet ))
2034       {
2035         IntersectNewEdges( cf );
2036         cf.MakeLoops( loopSet, normal );
2037       }
2038       // erase loops that are cut off by face intersections
2039       cf.CutOffLoops( loopSet, theSign, myNormals, cutOffLinks, cutOffCoplanarLinks );
2040
2041       smIdType index = cf.myInitFace->GetID(); // index in theNew2OldFaces
2042
2043       const SMDS_MeshElement* tria;
2044       for ( size_t iL = 0; iL < loopSet.myNbLoops; ++iL )
2045       {
2046         EdgeLoop& loop = loopSet.myLoops[ iL ];
2047         if ( loop.myLinks.size() == 0 )
2048           continue;
2049
2050         int nbTria  = triangulator.GetTriangles( &loop, nodes );
2051         int nbNodes = 3 * nbTria;
2052         for ( int i = 0; i < nbNodes; i += 3 )
2053         {
2054           if ( nodes[i] == nodes[i+1] || nodes[i] == nodes[i+2] || nodes[i+1] == nodes[i+2] )
2055           {
2056             if (SALOME::VerbosityActivated())
2057             {
2058               std::cerr << "BAD tria" << std::endl;
2059               cf.Dump();
2060             }
2061             else
2062             {
2063               if ( i < 0 ) cf.Dump(); // avoid "CutFace::Dump() unused in release mode"
2064             }
2065             continue;
2066           }
2067           if (!( tria = myMesh->FindFace( nodes[i], nodes[i+1], nodes[i+2] )))
2068             tria = myMesh->AddFace( nodes[i], nodes[i+1], nodes[i+2] );
2069           tria->setIsMarked( true ); // not to remove it
2070
2071           new2OldTria = std::make_pair( tria, theNew2OldFaces[ index ].second );
2072           if ( tria->GetID() < (int)theNew2OldFaces.size() )
2073             theNew2OldFaces[ tria->GetID() ] = new2OldTria;
2074           else
2075             theNew2OldFaces.push_back( new2OldTria );
2076
2077           if ( index == tria->GetID() )
2078             index = 0; // do not remove tria
2079         }
2080       }
2081       theNew2OldFaces[ index ].first = 0;
2082     }
2083
2084     // remove split faces
2085     for ( size_t id = 1; id < theNew2OldFaces.size(); ++id )
2086     {
2087       if ( theNew2OldFaces[id].first ||
2088            theNew2OldFaces[id].second == 0 )
2089         continue;
2090       if ( const SMDS_MeshElement* f = myMesh->FindElement( id ))
2091         myMesh->RemoveFreeElement( f );
2092     }
2093
2094     // remove faces that are merged off
2095     for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
2096     {
2097       const CutFace& cf = *cutFacesIt;
2098       if ( !cf.myLinks.empty() || cf.myInitFace->IsNull() )
2099         continue;
2100
2101       nodes.assign( cf.myInitFace->begin_nodes(), cf.myInitFace->end_nodes() );
2102       for ( size_t i = 0; i < nodes.size(); ++i )
2103       {
2104         const SMDS_MeshNode* n = nodes[ i ];
2105         while ( myRemove2KeepNodes.IsBound( n ))
2106           n = myRemove2KeepNodes( n );
2107         if ( n != nodes[ i ] && cf.myInitFace->GetNodeIndex( n ) >= 0 )
2108         {
2109           theNew2OldFaces[ cf.myInitFace->GetID() ].first = 0;
2110           myMesh->RemoveFreeElement( cf.myInitFace );
2111           break;
2112         }
2113       }
2114     }
2115
2116     // remove faces connected to cut off parts of cf.myInitFace
2117
2118     nodes.resize(2);
2119     for ( size_t i = 0; i < cutOffLinks.size(); ++i )
2120     {
2121       //break;
2122       nodes[0] = cutOffLinks[i].myNode1;
2123       nodes[1] = cutOffLinks[i].myNode2;
2124
2125       if ( nodes[0] != nodes[1] &&
2126            myMesh->GetElementsByNodes( nodes, faces ))
2127       {
2128         if ( // cutOffLinks[i].myFace &&
2129             cutOffLinks[i].myIndex != EdgePart::_COPLANAR &&
2130             faces.size() != 1 )
2131           continue;
2132         for ( size_t iF = 0; iF < faces.size(); ++iF )
2133         {
2134           smIdType index = faces[iF]->GetID();
2135           // if ( //faces[iF]->isMarked()         ||  // kept part of cutFace
2136           //      !theNew2OldFaces[ index ].first ) // already removed
2137           //   continue;
2138           cutFace.myInitFace = faces[iF];
2139           // if ( myCutFaces.Contains( cutFace )) // keep cutting faces needed in CutOffLoops()
2140           // {
2141           //   if ( !myCutFaces.Added( cutFace ).IsCut() )
2142           //     theNew2OldFaces[ index ].first = 0;
2143           //   continue;
2144           // }
2145           cutFace.myLinks.clear();
2146           cutFace.InitLinks();
2147           for ( size_t iL = 0; iL < cutFace.myLinks.size(); ++iL )
2148             if ( !cutOffLinks[i].IsSame( cutFace.myLinks[ iL ]))
2149               cutOffLinks.push_back( cutFace.myLinks[ iL ]);
2150
2151           theNew2OldFaces[ index ].first = 0;
2152           myMesh->RemoveFreeElement( faces[iF] );
2153         }
2154       }
2155     }
2156
2157     // replace nodes in touched faces
2158
2159     // treat touched faces
2160     for ( size_t i = 0; i < touchedFaces.size(); ++i )
2161     {
2162       const CutFace& cf = *touchedFaces[i];
2163       if ( cf.myInitFace->IsNull() )
2164         continue;
2165
2166       smIdType index = cf.myInitFace->GetID(); // index in theNew2OldFaces
2167       if ( !theNew2OldFaces[ index ].first )
2168         continue; // already cut off
2169
2170       cf.InitLinks();
2171       if ( !cf.ReplaceNodes( myRemove2KeepNodes ))
2172       {
2173         if ( cf.myLinks.size() == 3 &&
2174              cf.myInitFace->GetNodeIndex( cf.myLinks[0].myNode1 ) >= 0 &&
2175              cf.myInitFace->GetNodeIndex( cf.myLinks[1].myNode1 ) >= 0 &&
2176              cf.myInitFace->GetNodeIndex( cf.myLinks[2].myNode1 ) >= 0 )
2177           continue; // just keep as is
2178       }
2179
2180       if ( cf.myLinks.size() == 3 )
2181       {
2182         const SMDS_MeshElement* tria = myMesh->AddFace( cf.myLinks[0].myNode1,
2183                                                         cf.myLinks[1].myNode1,
2184                                                         cf.myLinks[2].myNode1 );
2185         new2OldTria = std::make_pair( tria, theNew2OldFaces[ index ].second );
2186         if ( tria->GetID() < (int)theNew2OldFaces.size() )
2187           theNew2OldFaces[ tria->GetID() ] = new2OldTria;
2188         else
2189           theNew2OldFaces.push_back( new2OldTria );
2190       }
2191       theNew2OldFaces[ index ].first = 0;
2192     }
2193
2194
2195     // add used new nodes to theNew2OldNodes
2196     SMESH_MeshAlgos::TNodeIntPairVec::value_type new2OldNode;
2197     new2OldNode.second = 0;
2198     for ( cutLinksIt = myCutLinks.cbegin(); cutLinksIt != myCutLinks.cend(); ++cutLinksIt )
2199     {
2200       const CutLink& link = *cutLinksIt;
2201       if ( link.IntNode() ) // && link.IntNode()->NbInverseElements() > 0 )
2202       {
2203         new2OldNode.first = link.IntNode();
2204         theNew2OldNodes.push_back( new2OldNode );
2205       }
2206     }
2207
2208     return;
2209   }
2210
2211   //================================================================================
2212   Intersector::Intersector( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals )
2213   {
2214     myAlgo = new Algo( mesh, tol, normals );
2215   }
2216   //================================================================================
2217   Intersector::~Intersector()
2218   {
2219     delete myAlgo;
2220   }
2221   //================================================================================
2222   //! compute cut of two faces of the mesh
2223   void Intersector::Cut( const SMDS_MeshElement* face1,
2224                          const SMDS_MeshElement* face2,
2225                          const int               nbCommonNodes )
2226   {
2227     myAlgo->Cut( face1, face2, nbCommonNodes );
2228   }
2229   //================================================================================
2230   //! store a face cut by a line given by its ends
2231   //  accompanied by indices of intersected face edges.
2232   //  Edge index is <0 if a line end is inside the face.
2233   void Intersector::Cut( const SMDS_MeshElement* face,
2234                          SMESH_NodeXYZ&          lineEnd1,
2235                          int                     edgeIndex1,
2236                          SMESH_NodeXYZ&          lineEnd2,
2237                          int                     edgeIndex2 )
2238   {
2239     myAlgo->Cut( face, lineEnd1, edgeIndex1, lineEnd2, edgeIndex2 );
2240   }
2241   //================================================================================
2242   //! split all face intersected by Cut() methods
2243   void Intersector::MakeNewFaces( SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces,
2244                                   SMESH_MeshAlgos::TNodeIntPairVec& theNew2OldNodes,
2245                                   const double                      theSign,
2246                                   const bool                        theOptimize )
2247   {
2248     myAlgo->MakeNewFaces( theNew2OldFaces, theNew2OldNodes, theSign, theOptimize );
2249   }
2250   //================================================================================
2251   //! Cut a face by planes, whose normals point to parts to keep
2252   bool Intersector::CutByPlanes(const SMDS_MeshElement*        theFace,
2253                                 const std::vector< gp_Ax1 > &  thePlanes,
2254                                 const double                   theTol,
2255                                 std::vector< TFace > &         theNewFaceConnectivity )
2256   {
2257     theNewFaceConnectivity.clear();
2258
2259     // check if theFace is wholly cut off
2260     std::vector< SMESH_NodeXYZ > facePoints( theFace->begin_nodes(), theFace->end_nodes() );
2261     facePoints.resize( theFace->NbCornerNodes() );
2262     for ( size_t iP = 0; iP < thePlanes.size(); ++iP )
2263     {
2264       size_t nbOut = 0;
2265       const gp_Pnt& O = thePlanes[iP].Location();
2266       for ( size_t i = 0; i < facePoints.size(); ++i )
2267       {
2268         gp_Vec Op( O, facePoints[i] );
2269         nbOut += ( Op * thePlanes[iP].Direction() <= 0 );
2270       }
2271       if ( nbOut == facePoints.size() )
2272         return true;
2273     }
2274
2275     // copy theFace into a temporary mesh
2276     SMDS_Mesh mesh;
2277     Bnd_B3d faceBox;
2278     std::vector< const SMDS_MeshNode* > faceNodes;
2279     faceNodes.resize( facePoints.size() );
2280     for ( size_t i = 0; i < facePoints.size(); ++i )
2281     {
2282       const SMESH_NodeXYZ& n = facePoints[i];
2283       faceNodes[i] = mesh.AddNode( n.X(), n.Y(), n.Z() );
2284       faceBox.Add( n );
2285     }
2286     const SMDS_MeshElement* faceToCut = 0;
2287     switch ( theFace->NbCornerNodes() )
2288     {
2289     case 3:
2290       faceToCut = mesh.AddFace( faceNodes[0], faceNodes[1], faceNodes[2] );
2291       break;
2292     case 4:
2293       faceToCut = mesh.AddFace( faceNodes[0], faceNodes[1], faceNodes[2], faceNodes[3] );
2294       break;
2295     default:
2296       faceToCut = mesh.AddPolygonalFace( faceNodes );
2297     }
2298
2299     std::vector< gp_XYZ > normals( 2 + thePlanes.size() );
2300     SMESH_MeshAlgos::FaceNormal( faceToCut, normals[ faceToCut->GetID() ]);
2301
2302     // add faces corresponding to thePlanes
2303     std::vector< const SMDS_MeshElement* > planeFaces;
2304     double faceSize = Sqrt( faceBox.SquareExtent() );
2305     gp_XYZ   center = 0.5 * ( faceBox.CornerMin() + faceBox.CornerMax() );
2306     for ( size_t i = 0; i < thePlanes.size(); ++i )
2307     {
2308       gp_Ax2 plnAx( thePlanes[i].Location(), thePlanes[i].Direction() );
2309       gp_XYZ O = plnAx.Location().XYZ();
2310       gp_XYZ X = plnAx.XDirection().XYZ();
2311       gp_XYZ Y = plnAx.YDirection().XYZ();
2312       gp_XYZ Z = plnAx.Direction().XYZ();
2313
2314       double dot = ( O - center ) * Z;
2315       gp_XYZ o = center + Z * dot; // center projected to a plane
2316
2317       gp_XYZ p1 = o + X * faceSize * 2;
2318       gp_XYZ p2 = o + Y * faceSize * 2;
2319       gp_XYZ p3 = o - (X + Y ) * faceSize * 2;
2320
2321       const SMDS_MeshNode* n1 = mesh.AddNode( p1.X(), p1.Y(), p1.Z() );
2322       const SMDS_MeshNode* n2 = mesh.AddNode( p2.X(), p2.Y(), p2.Z() );
2323       const SMDS_MeshNode* n3 = mesh.AddNode( p3.X(), p3.Y(), p3.Z() );
2324       planeFaces.push_back( mesh.AddFace( n1, n2, n3 ));
2325
2326       normals[ planeFaces.back()->GetID() ] = thePlanes[i].Direction().XYZ();
2327     }
2328
2329     // cut theFace
2330     Algo algo ( &mesh, theTol, normals );
2331     for ( size_t i = 0; i < planeFaces.size(); ++i )
2332     {
2333       algo.Cut( faceToCut, planeFaces[i], 0 );
2334     }
2335
2336     // retrieve a result
2337     SMESH_MeshAlgos::TElemIntPairVec new2OldFaces;
2338     SMESH_MeshAlgos::TNodeIntPairVec new2OldNodes;
2339     TCutFaceMap::const_iterator cutFacesIt= algo.myCutFaces.cbegin();
2340     for ( ; cutFacesIt != algo.myCutFaces.cend(); ++cutFacesIt )
2341     {
2342       const CutFace& cf = *cutFacesIt;
2343       if ( cf.myInitFace != faceToCut )
2344         continue;
2345
2346       if ( !cf.IsCut() )
2347       {
2348         theNewFaceConnectivity.push_back( facePoints );
2349         break;
2350       }
2351
2352       // intersect cut lines
2353       algo.IntersectNewEdges( cf );
2354
2355       // form loops of new faces
2356       EdgeLoopSet loopSet;
2357       cf.MakeLoops( loopSet, normals[ faceToCut->GetID() ]);
2358
2359       // erase loops that are cut off by thePlanes
2360       const double sign = 1;
2361       std::vector< EdgePart > cutOffLinks;
2362       TLinkMap                cutOffCoplanarLinks;
2363       cf.CutOffLoops( loopSet, sign, normals, cutOffLinks, cutOffCoplanarLinks );
2364
2365       for ( size_t iL = 0; iL < loopSet.myNbLoops; ++iL )
2366       {
2367         EdgeLoop& loop = loopSet.myLoops[ iL ];
2368         if ( loop.myLinks.size() > 0 )
2369         {
2370           facePoints.clear();
2371           for ( SMDS_NodeIteratorPtr nIt = loop.nodeIterator(); nIt->more(); )
2372           {
2373             const SMDS_MeshNode* n = nIt->next();
2374             facePoints.push_back( n );
2375             int iN = faceToCut->GetNodeIndex( n );
2376             if ( iN < 0 )
2377               facePoints.back()._node = 0; // an intersection point
2378             else
2379               facePoints.back()._node = theFace->GetNode( iN );
2380           }
2381           theNewFaceConnectivity.push_back( facePoints );
2382         }
2383       }
2384       break;
2385     }
2386
2387     return theNewFaceConnectivity.empty();
2388   }
2389
2390 } // namespace SMESH_MeshAlgos
2391
2392 namespace
2393 {
2394   //================================================================================
2395   /*!
2396    * \brief Debug
2397    */
2398   //================================================================================
2399
2400   void CutFace::Dump() const
2401   {
2402     std::cout << std::endl << "INI F " << myInitFace->GetID() << std::endl;
2403     for ( size_t i = 0; i < myLinks.size(); ++i )
2404       std::cout << "[" << i << "] ("
2405                 << char(( myLinks[i].IsInternal() ? 'j' : '0' ) + myLinks[i].myIndex ) << ") "
2406                 << myLinks[i].myNode1->GetID() << " - " << myLinks[i].myNode2->GetID()
2407                 << " " << ( myLinks[i].myFace ? 'F' : 'C' )
2408                 << ( myLinks[i].myFace ? myLinks[i].myFace->GetID() : 0 ) << " " << std::endl;
2409   }
2410
2411   //================================================================================
2412   /*!
2413    * \brief Add an edge cutting this face
2414    *  \param [in] p1 - start point of the edge
2415    *  \param [in] p2 - end point of the edge
2416    *  \param [in] cutter - a face producing the added cut edge.
2417    *  \param [in] nbOnPlane - nb of triangle nodes lying on the plane of the cutter face
2418    */
2419   //================================================================================
2420
2421   void CutFace::AddEdge( const CutLink&          p1,
2422                          const CutLink&          p2,
2423                          const SMDS_MeshElement* cutterFace,
2424                          const int               nbOnPlane,
2425                          const int               iNotOnPlane) const
2426   {
2427     int iN[2] = { myInitFace->GetNodeIndex( p1.IntNode() ),
2428                   myInitFace->GetNodeIndex( p2.IntNode() ) };
2429     if ( iN[0] >= 0 && iN[1] >= 0 )
2430     {
2431       // the cutting edge is a whole side
2432       if ((  cutterFace && nbOnPlane < 3 ) &&
2433           !( cutterFace->GetNodeIndex( p1.IntNode() ) >= 0 &&
2434              cutterFace->GetNodeIndex( p2.IntNode() ) >= 0 ))
2435       {
2436         InitLinks();
2437         myLinks[ Abs( iN[0] - iN[1] ) == 1 ? Min( iN[0], iN[1] ) : 2 ].myFace = cutterFace;
2438       }
2439       return;
2440     }
2441
2442     if ( p1.IntNode() == p2.IntNode() )
2443     {
2444       AddPoint( p1, p2, 1e-10 );
2445       return;
2446     }
2447
2448     InitLinks();
2449
2450     // cut side edges by a new one
2451
2452     int iEOnPlane = ( nbOnPlane == 2 ) ? ( iNotOnPlane + 1 ) % 3  :  -1;
2453
2454     double dist[2];
2455     for ( int is2nd = 0; is2nd < 2; ++is2nd )
2456     {
2457       const CutLink& p = is2nd ? p2 : p1;
2458       dist[ is2nd ] = 0;
2459       if ( iN[ is2nd ] >= 0 )
2460         continue;
2461
2462       int iE = Max( iEOnPlane, myInitFace->GetNodeIndex( p.Node1() ));
2463       if ( iE < 0 )
2464         continue; // link of other face
2465
2466       SMESH_NodeXYZ n0 = myLinks[iE].myNode1;
2467       dist[ is2nd ]    = ( n0 - p.myIntNode ).SquareModulus();
2468
2469       for ( size_t i = 0; i < myLinks.size(); ++i )
2470         if ( myLinks[i].myIndex == iE )
2471         {
2472           double d1 = n0.SquareDistance( myLinks[i].myNode1 );
2473           if ( d1 < dist[ is2nd ] )
2474           {
2475             double d2 = n0.SquareDistance( myLinks[i].myNode2 );
2476             if ( dist[ is2nd ] < d2 )
2477             {
2478               myLinks.push_back( myLinks[i] );
2479               myLinks.back().myNode1 = myLinks[i].myNode2 = p.IntNode();
2480               break;
2481             }
2482           }
2483         }
2484     }
2485
2486     int state = nbOnPlane == 3 ? EdgePart::_COPLANAR : EdgePart::_INTERNAL;
2487
2488     // look for an existing equal edge
2489     if ( nbOnPlane == 2 )
2490     {
2491       SMESH_NodeXYZ n0 = myLinks[ iEOnPlane ].myNode1;
2492       if ( iN[0] >= 0 ) dist[0] = ( n0 - p1.myIntNode ).SquareModulus();
2493       if ( iN[1] >= 0 ) dist[1] = ( n0 - p2.myIntNode ).SquareModulus();
2494       if ( dist[0] > dist[1] )
2495         std::swap( dist[0], dist[1] );
2496       for ( size_t i = 0; i < myLinks.size(); ++i )
2497       {
2498         if ( myLinks[i].myIndex != iEOnPlane )
2499           continue;
2500         gp_XYZ mid = 0.5 * ( SMESH_NodeXYZ( myLinks[i].myNode1 ) +
2501                              SMESH_NodeXYZ( myLinks[i].myNode2 ));
2502         double d = ( n0 - mid ).SquareModulus();
2503         if ( dist[0] < d && d < dist[1] )
2504           myLinks[i].myFace = cutterFace;
2505       }
2506       return;
2507     }
2508     else
2509     {
2510       EdgePart newEdge; newEdge.Set( p1.IntNode(), p2.IntNode(), cutterFace, state );
2511       for ( size_t i = 0; i < myLinks.size(); ++i )
2512       {
2513         if ( myLinks[i].IsSame( newEdge ))
2514         {
2515           // if ( !myLinks[i].IsInternal() )
2516           //   myLinks[ i ].myFace = cutterFace;
2517           // else
2518           myLinks[ i ].ReplaceCoplanar( newEdge );
2519           if ( myLinks[i].IsInternal() && i+1 < myLinks.size() )
2520             myLinks[ i+1 ].ReplaceCoplanar( newEdge );
2521           return;
2522         }
2523         i += myLinks[i].IsInternal();
2524       }
2525     }
2526
2527     size_t  i = myLinks.size();
2528     myLinks.resize( i + 2 );
2529     myLinks[ i   ].Set( p1.IntNode(), p2.IntNode(), cutterFace, state );
2530     myLinks[ i+1 ].Set( p2.IntNode(), p1.IntNode(), cutterFace, state );
2531   }
2532
2533   //================================================================================
2534   /*!
2535    * \brief Add a point cutting this face
2536    */
2537   //================================================================================
2538
2539   void CutFace::AddPoint( const CutLink& p1, const CutLink& p2, double /*tol*/ ) const
2540   {
2541     if ( myInitFace->GetNodeIndex( p1.IntNode() ) >= 0 ||
2542          myInitFace->GetNodeIndex( p2.IntNode() ) >= 0 )
2543       return;
2544
2545     InitLinks();
2546
2547     const CutLink* link = &p1;
2548     int iE = myInitFace->GetNodeIndex( link->Node1() );
2549     if ( iE < 0 )
2550     {
2551       link = &p2;
2552       iE = myInitFace->GetNodeIndex( link->Node1() );
2553     }
2554     if ( iE >= 0 )
2555     {
2556       // cut an existing edge by the point
2557       SMESH_NodeXYZ n0 = link->Node1();
2558       double         d = ( n0 - link->myIntNode ).SquareModulus();
2559
2560       for ( size_t i = 0; i < myLinks.size(); ++i )
2561         if ( myLinks[i].myIndex == iE )
2562         {
2563           double d1 = n0.SquareDistance( myLinks[i].myNode1 );
2564           if ( d1 < d )
2565           {
2566             double d2 = n0.SquareDistance( myLinks[i].myNode2 );
2567             if ( d < d2 )
2568             {
2569               myLinks.push_back( myLinks[i] );
2570               myLinks.back().myNode1 = myLinks[i].myNode2 = link->IntNode();
2571               return;
2572             }
2573           }
2574         }
2575     }
2576     else // point is inside the triangle
2577     {
2578       // // check if a point already added
2579       // for ( size_t i = 3; i < myLinks.size(); ++i )
2580       //   if ( myLinks[i].myNode1 == p1.IntNode() )
2581       //     return;
2582
2583       // // create a link between the point and the closest corner node
2584       // const SMDS_MeshNode* closeNode = myLinks[0].myNode1;
2585       // double minDist = p1.myIntNode.SquareDistance( closeNode );
2586       // for ( int i = 1; i < 3; ++i )
2587       // {
2588       //   double dist = p1.myIntNode.SquareDistance( myLinks[i].myNode1 );
2589       //   if ( dist < minDist )
2590       //   {
2591       //     minDist = dist;
2592       //     closeNode = myLinks[i].myNode1;
2593       //   }
2594       // }
2595       // if ( minDist > tol * tol )
2596       // {
2597       //   size_t i = myLinks.size();
2598       //   myLinks.resize( i + 2 );
2599       //   myLinks[ i   ].Set( p1.IntNode(), closeNode );
2600       //   myLinks[ i+1 ].Set( closeNode, p1.IntNode() );
2601       // }
2602     }
2603   }
2604
2605   //================================================================================
2606   /*!
2607    * \brief Perform node replacement
2608    */
2609   //================================================================================
2610
2611   bool CutFace::ReplaceNodes( const TNNMap& theRm2KeepMap ) const
2612   {
2613     bool replaced = false;
2614     for ( size_t i = 0; i < myLinks.size(); ++i )
2615     {
2616       while ( theRm2KeepMap.IsBound( myLinks[i].myNode1 ))
2617         replaced = ( myLinks[i].myNode1 = theRm2KeepMap( myLinks[i].myNode1 ));
2618
2619       while ( theRm2KeepMap.IsBound( myLinks[i].myNode2 ))
2620         replaced = ( myLinks[i].myNode2 = theRm2KeepMap( myLinks[i].myNode2 ));
2621     }
2622
2623     //if ( replaced ) // remove equal links
2624     {
2625       for ( size_t i1 = 0; i1 < myLinks.size(); ++i1 )
2626       {
2627         if ( myLinks[i1].myNode1 == myLinks[i1].myNode2 )
2628         {
2629           myLinks.erase( myLinks.begin() + i1,
2630                          myLinks.begin() + i1 + 1 + myLinks[i1].IsInternal() );
2631           --i1;
2632           continue;
2633         }
2634         size_t i2 = i1 + 1 + myLinks[i1].IsInternal();
2635         for ( ; i2 < myLinks.size(); ++i2 )
2636         {
2637           if ( !myLinks[i2].IsInternal() )
2638             continue;
2639           if ( myLinks[i1].IsSame( myLinks[i2] ))
2640           {
2641             myLinks[i1].  ReplaceCoplanar( myLinks[i2] );
2642             if ( myLinks[i1].IsInternal() )
2643               myLinks[i1+1].ReplaceCoplanar( myLinks[i2+1] );
2644             if ( !myLinks[i1].myFace && myLinks[i2].myFace )
2645             {
2646               myLinks[i1].  myFace = myLinks[i2].myFace;
2647               if ( myLinks[i1].IsInternal() )
2648                 myLinks[i1+1].myFace = myLinks[i2+1].myFace;
2649             }
2650             myLinks.erase( myLinks.begin() + i2,
2651                            myLinks.begin() + i2 + 2 );
2652             --i2;
2653             continue;
2654           }
2655           ++i2;
2656         }
2657         i1 += myLinks[i1].IsInternal();
2658       }
2659     }
2660
2661     return replaced;
2662   }
2663
2664   //================================================================================
2665   /*!
2666    * \brief Initialize myLinks with edges of myInitFace
2667    */
2668   //================================================================================
2669
2670   void CutFace::InitLinks() const
2671   {
2672     if ( !myLinks.empty() ) return;
2673
2674     int nbNodes = myInitFace->NbNodes();
2675     myLinks.reserve( nbNodes * 2 );
2676     myLinks.resize( nbNodes );
2677
2678     for ( int i = 0; i < nbNodes; ++i )
2679     {
2680       const SMDS_MeshNode* n1 = myInitFace->GetNode( i );
2681       const SMDS_MeshNode* n2 = myInitFace->GetNodeWrap( i + 1);
2682       myLinks[i].Set( n1, n2, 0, i );
2683     }
2684   }
2685
2686   //================================================================================
2687   /*!
2688    * \brief Return number of internal edges
2689    */
2690   //================================================================================
2691
2692   int CutFace::NbInternalEdges() const
2693   {
2694     int nb = 0;
2695     for ( size_t i = 3; i < myLinks.size(); ++i )
2696       nb += myLinks[i].IsInternal();
2697
2698     return nb / 2; // each internal edge encounters twice
2699   }
2700
2701   //================================================================================
2702   /*!
2703    * \brief Remove loops that are not connected to boundary edges of myFace by
2704    *        adding edges connecting these loops to the boundary.
2705    *        Such loops must be removed as they form polygons with ring topology.
2706    */
2707   //================================================================================
2708
2709   bool CutFace::RemoveInternalLoops( EdgeLoopSet& theLoops ) const
2710   {
2711     size_t nbReachedLoops = 0;
2712
2713     // count loops including boundary EdgeParts
2714     for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
2715     {
2716       EdgeLoop& loop = theLoops.myLoops[ iL ];
2717
2718       for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
2719         if ( !loop.myLinks[ iE ]->IsInternal() )
2720         {
2721           nbReachedLoops += loop.SetConnected();
2722           break;
2723         }
2724     }
2725     if ( nbReachedLoops == theLoops.myNbLoops )
2726       return false; // no unreachable loops
2727
2728
2729     // try to reach all loops by propagating via internal edges shared by loops
2730     size_t prevNbReached;
2731     do
2732     {
2733       prevNbReached = nbReachedLoops;
2734
2735       for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
2736       {
2737         EdgeLoop& loop = theLoops.myLoops[ iL ];
2738         if ( !loop.myIsBndConnected )
2739           continue;
2740
2741         for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
2742           if ( loop.myLinks[ iE ]->IsInternal() )
2743           {
2744             const EdgePart* twinEdge = getTwin( loop.myLinks[ iE ]);
2745             EdgeLoop*          loop2 = theLoops.GetLoopOf( twinEdge );
2746             if ( loop2->SetConnected() && ++nbReachedLoops == theLoops.myNbLoops )
2747               return false; // no unreachable loops
2748           }
2749       }
2750     }
2751     while ( prevNbReached < nbReachedLoops );
2752
2753
2754
2755     for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
2756     {
2757       EdgeLoop& loop = theLoops.myLoops[ iL ];
2758       if ( loop.myIsBndConnected || loop.myLinks.size() == 0 )
2759         continue;
2760
2761       if ( loop.myHasPending )
2762       {
2763         // try to join the loop to another one, with which it contacts at a node
2764
2765         // look for a node where the loop reverses
2766         const EdgePart* edgePrev = loop.myLinks.back();
2767         for ( size_t iE = 0; iE < loop.myLinks.size(); edgePrev = loop.myLinks[ iE++ ] )
2768         {
2769           if ( !edgePrev->IsTwin( *loop.myLinks[ iE ]))
2770             continue;
2771           const SMDS_MeshNode* reverseNode = edgePrev->myNode2;
2772
2773           // look for a loop including reverseNode
2774           size_t iContactEdge2; // index(+1) of edge starting at reverseNode
2775           for ( size_t iL2 = 0; iL2 < theLoops.myNbLoops; ++iL2 )
2776           {
2777             if ( iL == iL2 )
2778               continue;
2779             EdgeLoop& loop2 = theLoops.myLoops[ iL2 ];
2780             if ( ! ( iContactEdge2 = loop2.Contains( reverseNode )))
2781               continue;
2782
2783             // insert loop2 into the loop
2784             theLoops.Join( loop, iE, loop2, iContactEdge2 - 1 );
2785             break;
2786           }
2787           if ( loop.myIsBndConnected )
2788             break;
2789         }
2790
2791         if ( loop.myIsBndConnected )
2792           continue;
2793       }
2794
2795       // add links connecting internal loops with the boundary ones
2796
2797       // find a pair of closest nodes
2798       const SMDS_MeshNode *closestNode1 = 0, *closestNode2 = 0;
2799       double minDist = 1e100;
2800       for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
2801       {
2802         SMESH_NodeXYZ n1 = loop.myLinks[ iE ]->myNode1;
2803
2804         for ( size_t i = 0; i < myLinks.size(); ++i )
2805         {
2806           if ( !loop.Contains( myLinks[i].myNode1 ))
2807           {
2808             double dist = n1.SquareDistance( myLinks[i].myNode1 );
2809             if ( dist < minDist )
2810             {
2811               minDist = dist;
2812               closestNode1 = loop.myLinks[ iE ]->myNode1;
2813               closestNode2 = myLinks[i].myNode1;
2814             }
2815           }
2816           if ( myLinks[i].IsInternal() )
2817             ++i;
2818         }
2819       }
2820
2821       size_t i = myLinks.size();
2822       myLinks.resize( i + 2 );
2823       myLinks[ i   ].Set( closestNode1, closestNode2 );
2824       myLinks[ i+1 ].Set( closestNode2, closestNode1 );
2825     }
2826
2827     return true;
2828   }
2829
2830   //================================================================================
2831   /*!
2832    * \brief Return equal reversed edge
2833    */
2834   //================================================================================
2835
2836   EdgePart* CutFace::getTwin( const EdgePart* edge ) const
2837   {
2838     size_t i = edge - & myLinks[0];
2839
2840     if ( i > 2 && myLinks[ i-1 ].IsTwin( *edge ))
2841       return & myLinks[ i-1 ];
2842
2843     if ( i+1 < myLinks.size() &&
2844          myLinks[ i+1 ].IsTwin( *edge ))
2845       return & myLinks[ i+1 ];
2846
2847     return 0;
2848   }
2849
2850   //================================================================================
2851   /*!
2852    * \brief Fill loops of edges
2853    */
2854   //================================================================================
2855
2856   void CutFace::MakeLoops( EdgeLoopSet& theLoops, const gp_XYZ& theFaceNorm ) const
2857   {
2858     theLoops.Init( myLinks );
2859
2860     if ( myLinks.size() == 3 )
2861     {
2862       theLoops.AddNewLoop();
2863       theLoops.AddEdge( myLinks[0] );
2864       if ( myLinks[0].myNode2 == myLinks[1].myNode1 )
2865       {
2866         theLoops.AddEdge( myLinks[1] );
2867         theLoops.AddEdge( myLinks[2] );
2868       }
2869       else
2870       {
2871         theLoops.AddEdge( myLinks[2] );
2872         theLoops.AddEdge( myLinks[1] );
2873       }
2874       return;
2875     }
2876
2877     while ( !theLoops.AllEdgesUsed() )
2878     {
2879       EdgeLoop& loop = theLoops.AddNewLoop();
2880
2881       // add 1st edge to a new loop
2882       size_t i1;
2883       for ( i1 = theLoops.myNbLoops - 1; i1 < myLinks.size(); ++i1 )
2884         if ( theLoops.AddEdge( myLinks[i1] ))
2885           break;
2886
2887       EdgePart*             lastEdge = & myLinks[ i1 ];
2888       EdgePart*             twinEdge = getTwin( lastEdge );
2889       const SMDS_MeshNode* firstNode = lastEdge->myNode1;
2890       const SMDS_MeshNode*  lastNode = lastEdge->myNode2;
2891
2892       do // add the rest edges
2893       {
2894         theLoops.myCandidates.clear(); // edges starting at lastNode
2895         int nbInternal = 0;
2896
2897         // find candidate edges
2898         for ( size_t i = i1 + 1; i < myLinks.size(); ++i )
2899           if ( myLinks[ i ].myNode1 == lastNode &&
2900                &myLinks[ i ] != twinEdge &&
2901                !theLoops.myIsUsedEdge[ i ])
2902           {
2903             theLoops.myCandidates.push_back( & myLinks[ i ]);
2904             nbInternal += myLinks[ i ].IsInternal();
2905           }
2906
2907         // choose among candidates
2908         if ( theLoops.myCandidates.size() == 0 )
2909         {
2910           loop.myHasPending = bool( twinEdge );
2911           lastEdge = twinEdge;
2912         }
2913         else if ( theLoops.myCandidates.size() == 1 )
2914         {
2915           lastEdge = theLoops.myCandidates[0];
2916         }
2917         else if ( nbInternal == 1 && !lastEdge->IsInternal() )
2918         {
2919           lastEdge = theLoops.myCandidates[ !theLoops.myCandidates[0]->IsInternal() ];
2920         }
2921         else
2922         {
2923           gp_Vec  lastVec = *lastEdge;
2924           double maxAngle = -2 * M_PI;
2925           for ( size_t i = 0; i < theLoops.myCandidates.size(); ++i )
2926           {
2927             double angle = lastVec.AngleWithRef( *theLoops.myCandidates[i], theFaceNorm );
2928             if ( angle > maxAngle )
2929             {
2930               maxAngle = angle;
2931               lastEdge = theLoops.myCandidates[i];
2932             }
2933           }
2934         }
2935         theLoops.AddEdge( *lastEdge );
2936         lastNode = lastEdge->myNode2;
2937         twinEdge = getTwin( lastEdge );
2938       }
2939       while ( lastNode != firstNode );
2940
2941
2942       if ( twinEdge == & myLinks[ i1 ])
2943         loop.myHasPending = true;
2944
2945     } // while ( !theLoops.AllEdgesUsed() )
2946
2947     return;
2948   }
2949
2950   //================================================================================
2951   /*!
2952    * \brief Erase loops that are cut off by face intersections
2953    */
2954   //================================================================================
2955
2956   void CutFace::CutOffLoops( EdgeLoopSet&                 theLoops,
2957                              const double                 theSign,
2958                              const std::vector< gp_XYZ >& theNormals,
2959                              std::vector< EdgePart >&     theCutOffLinks,
2960                              TLinkMap&                    /*theCutOffCoplanarLinks*/) const
2961   {
2962     EdgePart sideEdge;
2963     boost::container::flat_set< const SMDS_MeshElement* > checkedCoplanar;
2964
2965     for ( size_t i = 0; i < myLinks.size(); ++i )
2966     {
2967       if ( !myLinks[i].myFace )
2968         continue;
2969
2970       EdgeLoop* loop = theLoops.GetLoopOf( & myLinks[i] );
2971       if ( !loop || loop->myLinks.empty() || loop->myHasPending )
2972         continue;
2973
2974       bool toErase, isCoplanar = ( myLinks[i].myIndex == EdgePart::_COPLANAR );
2975
2976       gp_Vec iniNorm = theNormals[ myInitFace->GetID() ];
2977       if ( isCoplanar )
2978       {
2979         toErase = ( myLinks[i].myFace->GetID() > myInitFace->GetID() );
2980
2981         const EdgePart* twin = getTwin( & myLinks[i] );
2982         if ( !twin || twin->myFace == myLinks[i].myFace )
2983         {
2984           // only one co-planar face includes myLinks[i]
2985           gp_Vec inFaceDir = iniNorm ^ myLinks[i];
2986           gp_XYZ   edgePnt = SMESH_NodeXYZ( myLinks[i].myNode1 );
2987           for ( int iN = 0; iN < 3; ++iN )
2988           {
2989             gp_Vec inCutFaceDir = ( SMESH_NodeXYZ( myLinks[i].myFace->GetNode( iN )) - edgePnt );
2990             if ( inCutFaceDir * inFaceDir < 0 )
2991             {
2992               toErase = false;
2993               break;
2994             }
2995           }
2996         }
2997       }
2998       else
2999       {
3000         gp_Vec   cutNorm = theNormals[ myLinks[i].myFace->GetID() ];
3001         gp_Vec inFaceDir = iniNorm ^ myLinks[i];
3002
3003         toErase = inFaceDir * cutNorm * theSign < 0;
3004         if ( !toErase )
3005         {
3006           // erase a neighboring loop
3007           loop = 0;
3008           if ( const EdgePart* twin = getTwin( & myLinks[i] ))
3009             loop = theLoops.GetLoopOf( twin );
3010           toErase = ( loop && !loop->myLinks.empty() );
3011         }
3012
3013         if ( toErase ) // do not erase if cutFace is connected to a co-planar cutFace
3014         {
3015           checkedCoplanar.clear();
3016           for ( size_t iE = 0; iE < myLinks.size() &&  toErase; ++iE )
3017           {
3018             if ( !myLinks[iE].myFace || myLinks[iE].myIndex != EdgePart::_COPLANAR )
3019               continue;
3020             bool isAdded = checkedCoplanar.insert( myLinks[iE].myFace ).second;
3021             if ( !isAdded )
3022               continue;
3023             toErase = ( SMESH_MeshAlgos::NbCommonNodes( myLinks[i ].myFace,
3024                                                         myLinks[iE].myFace ) < 1 );
3025           }
3026         }
3027       }
3028
3029       if ( toErase )
3030       {
3031         if ( !isCoplanar )
3032         {
3033           // remember whole sides of myInitFace that are cut off
3034           for ( size_t iE = 0; iE < loop->myLinks.size(); ++iE )
3035           {
3036             if ( !loop->myLinks[ iE ]->myFace              &&
3037                  !loop->myLinks[ iE ]->IsInternal()     )//   &&
3038               // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked
3039               // !loop->myLinks[ iE ]->myNode2->isMarked() )
3040             {
3041               int i = loop->myLinks[ iE ]->myIndex;
3042               sideEdge.Set( myInitFace->GetNode    ( i   ),
3043                             myInitFace->GetNodeWrap( i+1 ));
3044               theCutOffLinks.push_back( sideEdge );
3045
3046               if ( !sideEdge.IsSame( *loop->myLinks[ iE ] )) // nodes replaced
3047               {
3048                 theCutOffLinks.push_back( *loop->myLinks[ iE ] );
3049               }
3050             }
3051             else if ( IsCoplanar( loop->myLinks[ iE ]))
3052             {
3053               // propagate erasure to a co-planar face
3054               theCutOffLinks.push_back( *loop->myLinks[ iE ]);
3055             }
3056             else if ( loop->myLinks[ iE ]->myFace &&
3057                       loop->myLinks[ iE ]->IsInternal() )
3058               theCutOffLinks.push_back( *loop->myLinks[ iE ]);
3059           }
3060
3061           // clear the loop
3062           theLoops.Erase( loop );
3063         }
3064       }
3065     }
3066     return;
3067   }
3068
3069   //================================================================================
3070   /*!
3071    * \brief Check if the face has cut edges
3072    */
3073   //================================================================================
3074
3075   bool CutFace::IsCut() const
3076   {
3077     if ( myLinks.size() > 3 )
3078       return true;
3079
3080     if ( myLinks.size() == 3 )
3081       for ( size_t i = 0; i < 3; ++i )
3082         if ( myLinks[i].myFace )
3083           return true;
3084
3085     return false;
3086   }
3087
3088   //================================================================================
3089   /*!
3090    * \brief Check if an edge is produced by a co-planar cut
3091    */
3092   //================================================================================
3093
3094   bool CutFace::IsCoplanar( const EdgePart* edge ) const
3095   {
3096     if ( edge->myIndex == EdgePart::_COPLANAR )
3097     {
3098       const EdgePart* twin = getTwin( edge );
3099       return ( !twin || twin->myIndex == EdgePart::_COPLANAR );
3100     }
3101     return false;
3102   }
3103
3104   //================================================================================
3105   /*!
3106    * \brief Replace _COPLANAR cut edge by _INTERNAL or vice versa
3107    */
3108   //================================================================================
3109
3110   bool EdgePart::ReplaceCoplanar( const EdgePart& e )
3111   {
3112     if ( myIndex + e.myIndex == _COPLANAR + _INTERNAL )
3113     {
3114       //check if the faces are connected
3115       int nbCommonNodes = 0;
3116       if ( e.myFace && myFace )
3117         nbCommonNodes = SMESH_MeshAlgos::NbCommonNodes( e.myFace, myFace );
3118       bool toReplace = (( myIndex == _INTERNAL && nbCommonNodes > 1 ) ||
3119                         ( myIndex == _COPLANAR && nbCommonNodes < 2 ));
3120       if ( toReplace )
3121       {
3122         myIndex = e.myIndex;
3123         myFace  = e.myFace;
3124         return true;
3125       }
3126     }
3127     return false;
3128   }
3129
3130 } // namespace
3131
3132 //================================================================================
3133 /*!
3134  * \brief Create an offsetMesh of given faces
3135  *  \param [in] faceIt - the input faces
3136  *  \param [out] new2OldFaces - history of faces (new face -> old face ID)
3137  *  \param [out] new2OldNodes - history of nodes (new node -> old node ID)
3138  *  \return SMDS_Mesh* - the new offset mesh, a caller should delete
3139  */
3140 //================================================================================
3141
3142 SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
3143                                         SMDS_Mesh&           theSrcMesh,
3144                                         const double         theOffset,
3145                                         const bool           theFixIntersections,
3146                                         TElemIntPairVec&     theNew2OldFaces,
3147                                         TNodeIntPairVec&     theNew2OldNodes)
3148 {
3149   if ( theSrcMesh.GetMeshInfo().NbFaces( ORDER_QUADRATIC ) > 0 )
3150     throw SALOME_Exception( "Offset of quadratic mesh not supported" );
3151   if ( theSrcMesh.GetMeshInfo().NbFaces() > theSrcMesh.GetMeshInfo().NbTriangles() )
3152     throw SALOME_Exception( "Offset of non-triangular mesh not supported" );
3153
3154   SMDS_Mesh* newMesh = new SMDS_Mesh;
3155   theNew2OldFaces.clear();
3156   theNew2OldNodes.clear();
3157   theNew2OldFaces.push_back
3158     ( std::make_pair(( const SMDS_MeshElement*) 0, 0)); // to have index == face->GetID()
3159
3160   // copy input faces to the newMesh keeping IDs of nodes
3161
3162   double minNodeDist = 1e100;
3163
3164   std::vector< const SMDS_MeshNode* > nodes;
3165   while ( theFaceIt->more() )
3166   {
3167     const SMDS_MeshElement* face = theFaceIt->next();
3168     if ( face->GetType() != SMDSAbs_Face ) continue;
3169
3170     // copy nodes
3171     nodes.assign( face->begin_nodes(), face->end_nodes() );
3172     for ( size_t i = 0; i < nodes.size(); ++i )
3173     {
3174       const SMDS_MeshNode* newNode = newMesh->FindNode( nodes[i]->GetID() );
3175       if ( !newNode )
3176       {
3177         SMESH_NodeXYZ xyz( nodes[i] );
3178         newNode = newMesh->AddNodeWithID( xyz.X(), xyz.Y(), xyz.Z(), nodes[i]->GetID() );
3179         theNew2OldNodes.push_back( std::make_pair( newNode, nodes[i]->GetID() ));
3180         nodes[i] = newNode;
3181       }
3182     }
3183     const SMDS_MeshElement* newFace = 0;
3184     switch ( face->GetEntityType() )
3185     {
3186     case SMDSEntity_Triangle:
3187       newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2] );
3188       break;
3189     case SMDSEntity_Quad_Triangle:
3190       newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],
3191                                   nodes[3],nodes[4],nodes[5] );
3192       break;
3193     case SMDSEntity_BiQuad_Triangle:
3194       newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],
3195                                   nodes[3],nodes[4],nodes[5],nodes[6] );
3196       break;
3197     case SMDSEntity_Quadrangle:
3198       newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3] );
3199       break;
3200     case SMDSEntity_Quad_Quadrangle:
3201       newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3],
3202                                   nodes[4],nodes[5],nodes[6],nodes[7] );
3203       break;
3204     case SMDSEntity_BiQuad_Quadrangle:
3205       newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3],nodes[4],
3206                                   nodes[5],nodes[6],nodes[7],nodes[8] );
3207       break;
3208     case SMDSEntity_Polygon:
3209       newFace = newMesh->AddPolygonalFace( nodes );
3210       break;
3211     case SMDSEntity_Quad_Polygon:
3212       newFace = newMesh->AddQuadPolygonalFace( nodes );
3213       break;
3214     default:
3215       continue;
3216     }
3217     theNew2OldFaces.push_back( std::make_pair( newFace, face->GetID() ));
3218
3219     SMESH_NodeXYZ pPrev = nodes.back(), p;
3220     for ( size_t i = 0; i < nodes.size(); ++i )
3221     {
3222       p.Set( nodes[i] );
3223       double dist = ( pPrev - p ).SquareModulus();
3224       if ( dist < minNodeDist && dist > std::numeric_limits<double>::min() )
3225         minNodeDist = dist;
3226       pPrev = p;
3227     }
3228   } // while ( faceIt->more() )
3229
3230
3231   // compute normals to faces
3232   std::vector< gp_XYZ > normals( theNew2OldFaces.size() );
3233   for ( size_t i = 1; i < normals.size(); ++i )
3234   {
3235     if ( !SMESH_MeshAlgos::FaceNormal( theNew2OldFaces[i].first, normals[i] ))
3236       normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors
3237   }
3238
3239   const double sign = ( theOffset < 0 ? -1 : +1 );
3240   const double  tol = Min( 1e-3 * Sqrt( minNodeDist ),
3241                            1e-2 * theOffset * sign );
3242
3243   // translate new nodes by normal to input faces
3244   gp_XYZ newXYZ;
3245   std::vector< const SMDS_MeshNode* > multiNormalNodes;
3246   for ( size_t i = 0; i < theNew2OldNodes.size(); ++i )
3247   {
3248     const SMDS_MeshNode* newNode = theNew2OldNodes[i].first;
3249
3250     if ( getTranslatedPosition( newNode, theOffset, tol*10., sign, normals, theSrcMesh, newXYZ ))
3251       newMesh->MoveNode( newNode, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
3252     else
3253       multiNormalNodes.push_back( newNode );
3254   }
3255   // make multi-normal translation
3256   std::vector< SMESH_NodeXYZ > multiPos(10);
3257   for ( size_t i = 0; i < multiNormalNodes.size(); ++i )
3258   {
3259     const SMDS_MeshNode* newNode = multiNormalNodes[i];
3260     newNode->setIsMarked( true );
3261     SMESH_NodeXYZ oldXYZ = newNode;
3262     multiPos.clear();
3263     for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
3264     {
3265       const SMDS_MeshElement* newFace = fIt->next();
3266       const smIdType        faceIndex = newFace->GetID();
3267       const gp_XYZ&           oldNorm = normals[ faceIndex ];
3268       const gp_XYZ             newXYZ = oldXYZ + oldNorm * theOffset;
3269       if ( multiPos.empty() )
3270       {
3271         newMesh->MoveNode( newNode, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
3272         multiPos.emplace_back( newNode );
3273       }
3274       else
3275       {
3276         newNode = 0;
3277         for ( size_t iP = 0; iP < multiPos.size() &&  !newNode; ++iP )
3278           if (( multiPos[iP] - newXYZ ).SquareModulus() < tol * tol )
3279             newNode = multiPos[iP].Node();
3280         if ( !newNode )
3281         {
3282           newNode = newMesh->AddNode( newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
3283           newNode->setIsMarked( true );
3284           theNew2OldNodes.push_back( std::make_pair( newNode, 0 ));
3285           multiPos.emplace_back( newNode );
3286         }
3287       }
3288       if ( newNode != oldXYZ.Node() )
3289       {
3290         nodes.assign( newFace->begin_nodes(), newFace->end_nodes() );
3291         nodes[ newFace->GetNodeIndex( oldXYZ.Node() )] = newNode;
3292         newMesh->ChangeElementNodes( newFace, & nodes[0], nodes.size() );
3293       }
3294     }
3295   }
3296
3297   if ( !theFixIntersections )
3298     return newMesh;
3299
3300
3301   // remove new faces around concave nodes (they are marked) if the faces are inverted
3302   gp_XYZ faceNorm;
3303   for ( size_t i = 0; i < theNew2OldNodes.size(); ++i )
3304   {
3305     const SMDS_MeshNode* newNode = theNew2OldNodes[i].first;
3306     //const SMDS_MeshNode* oldNode = theNew2OldNodes[i].second;
3307     if ( newNode->isMarked() )
3308     {
3309       //gp_XYZ moveVec = sign * ( SMESH_NodeXYZ( newNode ) - SMESH_NodeXYZ( oldNode ));
3310
3311       //bool haveInverseFace = false;
3312       for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
3313       {
3314         const SMDS_MeshElement* newFace = fIt->next();
3315         const smIdType        faceIndex = newFace->GetID();
3316         const gp_XYZ&           oldNorm = normals[ faceIndex ];
3317         if ( !SMESH_MeshAlgos::FaceNormal( newFace, faceNorm, /*normalize=*/false ) ||
3318              //faceNorm * moveVec < 0 )
3319              faceNorm * oldNorm < 0 )
3320         {
3321           //haveInverseFace = true;
3322           theNew2OldFaces[ faceIndex ].first = 0;
3323           newMesh->RemoveFreeElement( newFace );
3324           //break;
3325         }
3326       }
3327       // if ( haveInverseFace )
3328       // {
3329       //   newMesh->MoveNode( newNode, oldNode->X(), oldNode->Y(), oldNode->Z() );
3330
3331       //   for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
3332       //   {
3333       //     const SMDS_MeshElement* newFace = fIt->next();
3334       //     if ( !SMESH_MeshAlgos::FaceNormal( newFace, normals[ newFace->GetID() ] ))
3335       //       normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors
3336       //   }
3337       // }
3338     }
3339     // mark all new nodes located closer than theOffset from theSrcMesh
3340   }
3341
3342   removeSmallFaces( newMesh, theNew2OldFaces, tol*tol );
3343
3344   // ==================================================
3345   // find self-intersections of new faces and fix them
3346   // ==================================================
3347
3348   std::unique_ptr< SMESH_ElementSearcher > fSearcher
3349     ( SMESH_MeshAlgos::GetElementSearcher( *newMesh, tol ));
3350
3351   Intersector intersector( newMesh, tol, normals );
3352
3353   std::vector< const SMDS_MeshElement* > closeFaces;
3354   std::vector< SMESH_NodeXYZ >           faceNodes;
3355   Bnd_B3d faceBox;
3356
3357   for ( size_t iF = 1; iF < theNew2OldFaces.size(); ++iF )
3358   {
3359     const SMDS_MeshElement* newFace = theNew2OldFaces[iF].first;
3360     if ( !newFace ) continue;
3361     faceNodes.assign( newFace->begin_nodes(), newFace->end_nodes() );
3362
3363     bool isConcaveNode1 = false;
3364     for ( size_t iN = 0; iN < faceNodes.size() && !isConcaveNode1; ++iN )
3365       isConcaveNode1 = faceNodes[iN]->isMarked();
3366
3367     // get faces close to a newFace
3368     closeFaces.clear();
3369     faceBox.Clear();
3370     for ( size_t i = 0; i < faceNodes.size(); ++i )
3371       faceBox.Add( faceNodes[i] );
3372     faceBox.Enlarge( tol );
3373
3374     fSearcher->GetElementsInBox( faceBox, SMDSAbs_Face, closeFaces );
3375
3376     // intersect the newFace with closeFaces
3377
3378     for ( size_t i = 0; i < closeFaces.size(); ++i )
3379     {
3380       const SMDS_MeshElement* closeFace = closeFaces[i];
3381       if ( closeFace->GetID() <= newFace->GetID() )
3382         continue; // this pair already treated
3383
3384       // do not intersect connected faces if they have no concave nodes
3385       int nbCommonNodes = 0;
3386       for ( size_t iN = 0; iN < faceNodes.size(); ++iN )
3387         nbCommonNodes += ( closeFace->GetNodeIndex( faceNodes[iN].Node() ) >= 0 );
3388
3389       if ( !isConcaveNode1 )
3390       {
3391         bool isConcaveNode2 = false;
3392         for ( SMDS_ElemIteratorPtr nIt = closeFace->nodesIterator(); nIt->more(); )
3393           if (( isConcaveNode2 = nIt->next()->isMarked() ))
3394             break;
3395
3396         if ( !isConcaveNode2 && nbCommonNodes > 0 )
3397         {
3398           if ( normals[ newFace->GetID() ] * normals[ closeFace->GetID() ] < 1.0 )
3399             continue; // not co-planar
3400         }
3401       }
3402
3403       intersector.Cut( newFace, closeFace, nbCommonNodes );
3404     }
3405   }
3406   intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign, /*optimize=*/true );
3407
3408   return newMesh;
3409 }