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