src/SMESHUtils/SMESH_Offset.cxx

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