src/SMESHUtils/SMESH_Offset.cxx

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