]> SALOME platform Git repositories - modules/smesh.git/blob - src/SMESHUtils/SMESH_Slot.cxx
Salome HOME
Merge branch 'master' into gni/adaptation
[modules/smesh.git] / src / SMESHUtils / SMESH_Slot.cxx
1 // Copyright (C) 2018-2020  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // File      : SMESH_Slot.cxx
20 // Created   : Fri Nov 30 15:58:37 2018
21 // Author    : Edward AGAPOV (eap)
22
23 #include "SMESH_MeshAlgos.hxx"
24
25 #include "ObjectPool.hxx"
26 #include "SMDS_LinearEdge.hxx"
27 #include "SMDS_Mesh.hxx"
28 #include "SMDS_MeshGroup.hxx"
29
30 #include <IntAna_IntConicQuad.hxx>
31 #include <IntAna_Quadric.hxx>
32 #include <NCollection_DataMap.hxx>
33 #include <NCollection_Map.hxx>
34 #include <Precision.hxx>
35 #include <gp_Ax1.hxx>
36 #include <gp_Cylinder.hxx>
37 #include <gp_Dir.hxx>
38 #include <gp_Lin.hxx>
39 #include <gp_Pln.hxx>
40 #include <gp_Pnt.hxx>
41 #include <gp_Vec.hxx>
42
43 #include <Utils_SALOME_Exception.hxx>
44
45 namespace
46 {
47   typedef SMESH_MeshAlgos::Edge TEdge;
48   
49   //================================================================================
50   //! point of intersection of a face edge with the cylinder
51   struct IntPoint
52   {
53     SMESH_NodeXYZ myNode;        // point and a node
54     int           myEdgeIndex;   // face edge index
55     bool          myIsOutPln[2]; // isOut of two planes
56
57     double SquareDistance( const IntPoint& p ) const { return ( myNode-p.myNode ).SquareModulus(); }
58   };
59
60   //================================================================================
61   //! cut of a face
62   struct Cut
63   {
64     IntPoint myIntPnt1, myIntPnt2;
65     const SMDS_MeshElement* myFace;
66
67     const IntPoint& operator[]( size_t i ) const { return i ? myIntPnt2 : myIntPnt1; }
68
69     double SquareDistance( const gp_Pnt& p, gp_XYZ & pClosest ) const
70     {
71       gp_Vec edge( myIntPnt1.myNode, myIntPnt2.myNode );
72       gp_Vec n1p ( myIntPnt1.myNode, p  );
73       double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
74       if ( u <= 0. )
75       {
76         pClosest = myIntPnt1.myNode;
77         return n1p.SquareMagnitude();
78       }
79       if ( u >= 1. )
80       {
81         pClosest = myIntPnt2.myNode;
82         return p.SquareDistance( myIntPnt2.myNode );
83       }
84       pClosest = myIntPnt1.myNode + u * edge.XYZ(); // projection of the point on the edge
85       return p.SquareDistance( pClosest );
86     }
87   };
88
89   //================================================================================
90   //! poly-line segment
91   struct Segment
92   {
93     typedef std::vector< Cut > TCutList;
94
95     const SMDS_MeshElement*        myEdge;
96     TCutList                       myCuts;
97     std::vector< const IntPoint* > myFreeEnds; // ends of cut edges
98
99     Segment( const SMDS_MeshElement* e = 0 ): myEdge(e) { myCuts.reserve( 4 ); }
100
101     // return its axis
102     gp_Ax1 Ax1( bool reversed = false ) const
103     {
104       SMESH_NodeXYZ n1 = myEdge->GetNode(  reversed );
105       SMESH_NodeXYZ n2 = myEdge->GetNode( !reversed );
106       return gp_Ax1( n1, gp_Dir( n2 - n1 ));
107     }
108
109     // return a node
110     const SMDS_MeshNode* Node(int i) const
111     {
112       return myEdge->GetNode( i % 2 );
113     }
114
115     // store an intersection edge forming the slot border
116     void AddCutEdge( const IntPoint& p1,
117                      const IntPoint& p2,
118                      const SMDS_MeshElement* myFace )
119     {
120       myCuts.push_back( Cut({ p1, p2, myFace }));
121     }
122
123     // return number of not shared IntPoint's
124     int NbFreeEnds( double tol )
125     {
126       if ( myCuts.empty() )
127         return 0;
128       if ( myFreeEnds.empty() )
129       {
130         // remove degenerated cuts
131         // for ( size_t iC1 = 0; iC1 < myCuts.size(); ++iC1 )
132         //   if ( myCuts[ iC1 ][ 0 ].myNode == myCuts[ iC1 ][ 1 ].myNode )
133         //   {
134         //     if ( iC1 < myCuts.size() - 1 )
135         //       myCuts[ iC1 ] = myCuts.back();
136         //     myCuts.pop_back();
137         //   }
138
139         int nbShared = 0;
140         std::vector< bool > isSharedPnt( myCuts.size() * 2, false );
141         for ( size_t iC1 = 0; iC1 < myCuts.size() - 1; ++iC1 )
142           for ( size_t iP1 = 0; iP1 < 2; ++iP1 )
143           {
144             size_t i1 = iC1 * 2 + iP1;
145             if ( isSharedPnt[ i1 ])
146               continue;
147             for ( size_t iC2 = iC1 + 1; iC2 < myCuts.size(); ++iC2 )
148               for ( size_t iP2 = 0; iP2 < 2; ++iP2 )
149               {
150                 size_t i2 = iC2 * 2 + iP2;
151                 if ( isSharedPnt[ i2 ])
152                   continue;
153                 if ( myCuts[ iC1 ][ iP1 ].SquareDistance( myCuts[ iC2 ][ iP2 ]) < tol * tol )
154                 {
155                   nbShared += 2;
156                   if ( myCuts[ iC1 ][ 0 ].SquareDistance( myCuts[ iC1 ][ 1 ]) < tol * tol )
157                     isSharedPnt[ iC1 * 2 ] = isSharedPnt[ iC1 * 2 + 1 ] = true;
158                   else if ( myCuts[ iC2 ][ 0 ].SquareDistance( myCuts[ iC2 ][ 1 ]) < tol * tol )
159                     isSharedPnt[ iC2 * 2 ] = isSharedPnt[ iC2 * 2 + 1 ] = true;
160                   else
161                     isSharedPnt[ i1 ] = isSharedPnt[ i2 ] = true;
162                 }
163               }
164           }
165         myFreeEnds.reserve( isSharedPnt.size() - nbShared );
166         for ( size_t i = 0; i < isSharedPnt.size(); ++i )
167           if ( !isSharedPnt[ i ] )
168           {
169             int iP = i % 2;
170             int iC = i / 2;
171             myFreeEnds.push_back( & myCuts[ iC ][ iP ]);
172           }
173       }
174       return myFreeEnds.size();
175     }
176   };
177   typedef ObjectPoolIterator<Segment> TSegmentIterator;
178
179
180   //================================================================================
181   //! Segments and plane separating domains of segments, at common node
182   struct NodeData
183   {
184     std::vector< Segment* > mySegments;
185     gp_Ax1                  myPlane; // oriented OK for mySegments[0]
186
187     void AddSegment( Segment* seg, const SMDS_MeshNode* n )
188     {
189       mySegments.reserve(2);
190       mySegments.push_back( seg );
191       if ( mySegments.size() == 1 )
192       {
193         myPlane = mySegments[0]->Ax1( mySegments[0]->myEdge->GetNodeIndex( n ));
194       }
195       else
196       {
197         gp_Ax1 axis2 = mySegments[1]->Ax1( mySegments[1]->myEdge->GetNodeIndex( n ));
198         myPlane.SetDirection( myPlane.Direction().XYZ() - axis2.Direction().XYZ() );
199       }
200     }
201     gp_Ax1 Plane( const Segment* seg )
202     {
203       return ( seg == mySegments[0] ) ? myPlane : myPlane.Reversed();
204     }
205   };
206   typedef NCollection_DataMap< const SMDS_MeshNode*, NodeData, SMESH_Hasher > TSegmentsOfNode;
207
208
209   //================================================================================
210   /*!
211    * \brief Intersect a face edge given by its nodes with a cylinder.
212    */
213   //================================================================================
214
215   bool intersectEdge( const gp_Cylinder&      cyl,
216                       const SMESH_NodeXYZ&    n1,
217                       const SMESH_NodeXYZ&    n2,
218                       const double            tol,
219                       std::vector< IntPoint >& intPoints )
220   {
221     gp_Lin line( gp_Ax1( n1, gp_Dir( n2 - n1 )));
222     IntAna_IntConicQuad intersection( line, IntAna_Quadric( cyl ));
223
224     if ( !intersection.IsDone()     ||
225          intersection.IsParallel()  ||
226          intersection.IsInQuadric() ||
227          intersection.NbPoints() == 0 )
228       return false;
229
230     gp_Vec edge( n1, n2 );
231
232     size_t oldNbPnts = intPoints.size();
233     for ( int iP = 1; iP <= intersection.NbPoints(); ++iP )
234     {
235       const gp_Pnt& p = intersection.Point( iP );
236
237       gp_Vec n1p ( n1, p );
238       const SMDS_MeshNode* n = 0;
239
240       double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
241       if ( u <= 0. ) {
242         if ( p.SquareDistance( n1 ) < tol * tol )
243           n = n1.Node();
244         else
245           continue;
246       }
247       else if ( u >= 1. ) {
248         if ( p.SquareDistance( n2 ) < tol * tol )
249           n = n2.Node();
250         else
251           continue;
252       }
253       else {
254         if      ( p.SquareDistance( n1 ) < tol * tol )
255           n = n1.Node();
256         else if ( p.SquareDistance( n2 ) < tol * tol )
257           n = n2.Node();
258       }
259
260       intPoints.push_back( IntPoint() );
261       if ( n )
262         intPoints.back().myNode.Set( n );
263       else
264         intPoints.back().myNode.SetCoord( p.X(),p.Y(),p.Z() );
265     }
266
267     // set points order along an edge
268     if ( intPoints.size() - oldNbPnts == 2 &&
269          intersection.ParamOnConic( 1 ) > intersection.ParamOnConic( 2 ))
270     {
271       int i = intPoints.size() - 1;
272       std::swap( intPoints[ i ], intPoints[ i - 1 ]);
273     }
274
275     return intPoints.size() - oldNbPnts > 0;
276   }
277
278   //================================================================================
279   /*!
280    * \brief Return signed distance between a point and a plane
281    */
282   //================================================================================
283
284   double signedDist( const gp_Pnt& p, const gp_Ax1& planeNormal )
285   {
286     const gp_Pnt& O = planeNormal.Location();
287     gp_Vec Op( O, p );
288     return Op * planeNormal.Direction();
289   }
290
291   //================================================================================
292   /*!
293    * \brief Check if a point is outside a segment domain bound by two planes
294    */
295   //================================================================================
296
297   bool isOut( const gp_Pnt& p, const gp_Ax1* planeNormal, bool* isOutPtr, int nbPln = 2 )
298   {
299     isOutPtr[0] = isOutPtr[1] = false;
300
301     for ( int i = 0; i < nbPln; ++i )
302     {
303       isOutPtr[i] = ( signedDist( p, planeNormal[i] ) <= 0. );
304     }
305     return ( isOutPtr[0] && isOutPtr[1] );
306   }
307
308   //================================================================================
309   /*!
310    * \brief Check if a segment between two points is outside a segment domain bound by two planes
311    */
312   //================================================================================
313
314   bool isSegmentOut( bool* isOutPtr1, bool* isOutPtr2 )
315   {
316     return (( isOutPtr1[0] && isOutPtr2[0] ) ||
317             ( isOutPtr1[1] && isOutPtr2[1] ));
318   }
319
320   //================================================================================
321   /*!
322    * \brief cut off ip1 from edge (ip1 - ip2) by a plane
323    */
324   //================================================================================
325
326   void cutOff( IntPoint & ip1, const IntPoint & ip2, const gp_Ax1& planeNormal, double tol )
327   {
328     gp_Lin lin( ip1.myNode, ( ip2.myNode - ip1.myNode ));
329     gp_Pln pln( planeNormal.Location(), planeNormal.Direction() );
330
331     IntAna_IntConicQuad intersection( lin, pln, Precision::Angular/*Tolerance*/() );
332     if ( intersection.IsDone()      &&
333          !intersection.IsParallel()  &&
334          !intersection.IsInQuadric() &&
335          intersection.NbPoints() == 1 )
336     {
337       if ( intersection.Point( 1 ).SquareDistance( ip1.myNode ) > tol * tol )
338       {
339         static_cast< gp_XYZ& >( ip1.myNode ) = intersection.Point( 1 ).XYZ();
340         ip1.myNode._node = 0;
341         ip1.myEdgeIndex = -1;
342       }
343     }
344   }
345
346   //================================================================================
347   /*!
348    * \brief Assure that face normal is computed in faceNormals vector 
349    */
350   //================================================================================
351
352   const gp_XYZ& computeNormal( const SMDS_MeshElement* face,
353                                std::vector< gp_XYZ >&  faceNormals )
354   {
355     bool toCompute;
356     if ((int) faceNormals.size() <= face->GetID() )
357     {
358       toCompute = true;
359       faceNormals.resize( face->GetID() + 1 );
360     }
361     else
362     {
363       toCompute = faceNormals[ face->GetID() ].SquareModulus() == 0.;
364     }
365     if ( toCompute )
366       SMESH_MeshAlgos::FaceNormal( face, faceNormals[ face->GetID() ], /*normalized=*/false );
367
368     return faceNormals[ face->GetID() ];
369   }
370
371   typedef std::vector< SMDS_MeshGroup* > TGroupVec;
372
373   //================================================================================
374   /*!
375    * \brief Fill theFaceID2Groups map for a given face
376    *  \param [in] theFace - the face
377    *  \param [in] theGroupsToUpdate - list of groups to treat
378    *  \param [out] theFaceID2Groups - the map to fill in
379    *  \param [out] theWorkGroups - a working buffer of groups
380    */
381   //================================================================================
382
383   void findGroups( const SMDS_MeshElement *                theFace,
384                    TGroupVec &                             theGroupsToUpdate,
385                    NCollection_DataMap< int, TGroupVec > & theFaceID2Groups,
386                    TGroupVec &                             theWorkGroups )
387   {
388     theWorkGroups.clear();
389     for ( size_t i = 0; i < theGroupsToUpdate.size(); ++i )
390       if ( theGroupsToUpdate[i]->Contains( theFace ))
391         theWorkGroups.push_back( theGroupsToUpdate[i] );
392
393     if ( !theWorkGroups.empty() )
394       theFaceID2Groups.Bind( theFace->GetID(), theWorkGroups );
395   }
396
397   //================================================================================
398   /*!
399    * \brief Check distance between a point and an edge defined by a couple of nodes
400    */
401   //================================================================================
402
403   bool isOnEdge( const SMDS_MeshNode* n1,
404                  const SMDS_MeshNode* n2,
405                  const gp_Pnt&        p,
406                  const double         tol )
407   {
408     SMDS_LinearEdge edge( n1, n2 );
409     return ( SMESH_MeshAlgos::GetDistance( &edge, p ) < tol );
410   }
411
412   //================================================================================
413   /*!
414    * \return Index of intersection point detected on a triangle cut by planes
415    *  \param [in] i - index of a cut triangle side
416    *  \param [in] n1 - 1st point of a cut triangle side
417    *  \param [in] n2 - 2nd point of a cut triangle side
418    *  \param [in] face - a not cut triangle
419    *  \param [in] intPoint - the intersection point
420    *  \param [in] faceNodes - nodes of not cut triangle
421    *  \param [in] tol - tolerance
422    */
423   //================================================================================
424
425   int edgeIndex( const int                                  i,
426                  const SMESH_NodeXYZ&                       n1,
427                  const SMESH_NodeXYZ&                       n2,
428                  const SMDS_MeshElement*                    face,
429                  const IntPoint&                            intPoint,
430                  const std::vector< const SMDS_MeshNode* >& faceNodes,
431                  const double                               tol )
432   {
433     if ( n1.Node() && n2.Node() )
434       return face->GetNodeIndex( n1.Node() );
435
436     // project intPoint to sides of face
437     for ( size_t i = 1; i < faceNodes.size(); ++i )
438       if ( isOnEdge( faceNodes[ i-1 ], faceNodes[ i ], intPoint.myNode, tol ))
439         return i - 1;
440
441     return -(i+1);
442   }
443
444   //================================================================================
445   /*!
446    * \brief Find a neighboring segment and its next node
447    *  \param [in] curSegment - a current segment
448    *  \param [in,out] curNode - a current node to update
449    *  \param [in] segmentsOfNode - map of segments of nodes
450    *  \return Segment* - the found segment
451    */
452   //================================================================================
453
454   Segment* nextSegment( const Segment*         curSegment,
455                         const SMDS_MeshNode* & curNode,
456                         const TSegmentsOfNode& segmentsOfNode )
457   {
458     Segment* neighborSeg = 0;
459     const NodeData& noData = segmentsOfNode( curNode );
460     for ( size_t iS = 0; iS < noData.mySegments.size()  && !neighborSeg; ++iS )
461       if ( noData.mySegments[ iS ] != curSegment )
462         neighborSeg = noData.mySegments[ iS ];
463
464     if ( neighborSeg )
465     {
466       int iN = ( neighborSeg->Node(0) == curNode );
467       curNode = neighborSeg->Node( iN );
468     }
469     return neighborSeg;
470   }
471
472   //================================================================================
473   /*!
474    * \brief Tries to find a segment to which a given point is too close
475    *  \param [in] p - the point
476    *  \param [in] minDist - minimal allowed distance from segment
477    *  \param [in] curSegment - start segment
478    *  \param [in] curNode - start node
479    *  \param [in] segmentsOfNode - map of segments of nodes
480    *  \return bool - true if a too close segment found
481    */
482   //================================================================================
483
484   const Segment* findTooCloseSegment( const IntPoint&        p,
485                                       const double           minDist,
486                                       const double           tol,
487                                       const Segment*         curSegment,
488                                       const SMDS_MeshNode*   curNode,
489                                       const TSegmentsOfNode& segmentsOfNode )
490   {
491     double prevDist = Precision::Infinite();
492     while ( curSegment )
493     {
494       double dist = SMESH_MeshAlgos::GetDistance( curSegment->myEdge, p.myNode );
495       if ( dist < minDist )
496       {
497         // check if dist is less than distance of curSegment to its cuts
498         // double minCutDist = prevDist;
499         // bool     coincide = false;
500         // for ( size_t iC = 0; iC < curSegment->myCuts.size(); ++iC )
501         // {
502         //   if (( coincide = ( curSegment->myCuts[iC].SquareDistance( p.myNode ) < tol * tol )))
503         //     break;
504         //   for ( size_t iP = 0; iP < 2; ++iP )
505         //   {
506         //     double cutDist = SMESH_MeshAlgos::GetDistance( curSegment->myEdge,
507         //                                                    curSegment->myCuts[iC][iP].myNode );
508         //     minCutDist = std::min( minCutDist, cutDist );
509         //   }
510         // }
511         // if ( !coincide && minCutDist > dist )
512         return curSegment;
513       }
514       if ( dist > prevDist )
515         break;
516       prevDist   = dist;
517       curSegment = nextSegment( curSegment, curNode, segmentsOfNode );
518     }
519     return 0;
520   }
521 }
522
523 //================================================================================
524 /*!
525  * \brief Create a slot of given width around given 1D elements lying on a triangle mesh.
526  * The slot is constructed by cutting faces by cylindrical surfaces made around each segment.
527  * \return Edges located at the slot boundary
528  */
529 //================================================================================
530
531 std::vector< SMESH_MeshAlgos::Edge >
532 SMESH_MeshAlgos::MakeSlot( SMDS_ElemIteratorPtr             theSegmentIt,
533                            double                           theWidth,
534                            SMDS_Mesh*                       theMesh,
535                            std::vector< SMDS_MeshGroup* > & theGroupsToUpdate)
536 {
537   std::vector< Edge > bndEdges;
538
539   if ( !theSegmentIt || !theSegmentIt->more() || !theMesh || theWidth == 0.)
540     return bndEdges;
541
542   // ----------------------------------------------------------------------------------
543   // put the input segments to a data map in order to be able finding neighboring ones
544   // ----------------------------------------------------------------------------------
545
546   TSegmentsOfNode segmentsOfNode;
547   ObjectPool< Segment > segmentPool;
548
549   while( theSegmentIt->more() )
550   {
551     const SMDS_MeshElement* edge = theSegmentIt->next();
552     if ( edge->GetType() != SMDSAbs_Edge )
553       throw SALOME_Exception( "A segment is not a mesh edge");
554
555     Segment* segment = segmentPool.getNew();
556     segment->myEdge = edge;
557
558     for ( SMDS_NodeIteratorPtr nIt = edge->nodeIterator(); nIt->more(); )
559     {
560       const SMDS_MeshNode* n = nIt->next();
561       NodeData* noData = segmentsOfNode.ChangeSeek( n );
562       if ( !noData )
563         noData = segmentsOfNode.Bound( n, NodeData() );
564       noData->AddSegment( segment, n );
565     }
566   }
567
568   // ---------------------------------
569   // Cut the mesh around the segments
570   // ---------------------------------
571
572   const double tol = Precision::Confusion();
573   const double angularTol = 1e-5;
574   std::vector< gp_XYZ > faceNormals;
575   SMESH_MeshAlgos::Intersector meshIntersector( theMesh, tol, faceNormals );
576   std::unique_ptr< SMESH_ElementSearcher> faceSearcher;
577
578   std::vector< NLink > startEdges;
579   std::vector< const SMDS_MeshNode* > faceNodes(4), edgeNodes(2);
580   std::vector<const SMDS_MeshElement *> faces(2);
581   NCollection_Map<const SMDS_MeshElement*, SMESH_Hasher > checkedFaces;
582   std::vector< IntPoint > intPoints, p(2);
583   std::vector< SMESH_NodeXYZ > facePoints(4);
584   std::vector< Intersector::TFace > cutFacePoints;
585
586   NCollection_DataMap< int, TGroupVec > faceID2Groups;
587   TGroupVec groupVec;
588
589   std::vector< gp_Ax1 > planeNormalVec(2);
590   gp_Ax1 * planeNormal = & planeNormalVec[0];
591   
592   for ( TSegmentIterator segIt( segmentPool ); segIt.more(); ) // loop on all segments
593   {
594     Segment* segment = const_cast< Segment* >( segIt.next() );
595
596     gp_Lin      segLine( segment->Ax1() );
597     gp_Ax3      cylAxis( segLine.Location(), segLine.Direction() );
598     gp_Cylinder segCylinder( cylAxis, 0.5 * theWidth );
599     double      radius2( segCylinder.Radius() * segCylinder.Radius() );
600
601     // get normals of planes separating domains of neighboring segments
602     for ( int i = 0; i < 2; ++i ) // loop on 2 segment ends
603     {
604       const SMDS_MeshNode* n = segment->Node( i );
605       planeNormal[i] = segmentsOfNode( n ).Plane( segment );
606     }
607
608     // we explore faces around a segment starting from face edges;
609     // initialize a list of starting edges
610     startEdges.clear();
611     {
612       // get a face to start searching intersected faces from
613       const SMDS_MeshNode*      n0 = segment->Node( 0 );
614       SMDS_ElemIteratorPtr     fIt = n0->GetInverseElementIterator( SMDSAbs_Face );
615       const SMDS_MeshElement* face = ( fIt->more() ) ? fIt->next() : 0;
616       if ( !theMesh->Contains( face ))
617       {
618         if ( !faceSearcher )
619           faceSearcher.reset( SMESH_MeshAlgos::GetElementSearcher( *theMesh ));
620         face = faceSearcher->FindClosestTo( SMESH_NodeXYZ( n0 ), SMDSAbs_Face );
621       }
622       // collect face edges
623       int nbNodes = face->NbCornerNodes();
624       faceNodes.assign( face->begin_nodes(), face->end_nodes() );
625       faceNodes.resize( nbNodes + 1 );
626       faceNodes[ nbNodes ] = faceNodes[ 0 ];
627       for ( int i = 0; i < nbNodes; ++i )
628         startEdges.push_back( NLink( faceNodes[i], faceNodes[i+1] ));
629     }
630
631     // intersect faces located around a segment
632     checkedFaces.Clear();
633     while ( !startEdges.empty() )
634     {
635       edgeNodes[0] = startEdges[0].first;
636       edgeNodes[1] = startEdges[0].second;
637
638       theMesh->GetElementsByNodes( edgeNodes, faces, SMDSAbs_Face );
639       for ( size_t iF = 0; iF < faces.size(); ++iF ) // loop on faces sharing a start edge
640       {
641         const SMDS_MeshElement* face = faces[iF];
642         if ( !checkedFaces.Add( face ))
643           continue;
644
645         int nbNodes = face->NbCornerNodes();
646         if ( nbNodes != 3 )
647           throw SALOME_Exception( "MakeSlot() accepts triangles only" );
648         faceNodes.assign( face->begin_nodes(), face->end_nodes() );
649         faceNodes.resize( nbNodes + 1 );
650         faceNodes[ nbNodes ] = faceNodes[ 0 ];
651         facePoints.assign( faceNodes.begin(), faceNodes.end() );
652
653         // check if cylinder axis || face
654         const gp_XYZ& faceNorm = computeNormal( face, faceNormals );
655         bool isCylinderOnFace  = ( Abs( faceNorm * cylAxis.Direction().XYZ() ) < tol );
656
657         if ( !isCylinderOnFace )
658         {
659           if ( Intersector::CutByPlanes( face, planeNormalVec, tol, cutFacePoints ))
660             continue; // whole face cut off
661           facePoints.swap( cutFacePoints[0] );
662           facePoints.push_back( facePoints[0] );
663         }
664
665         // find intersection points on face edges
666         intPoints.clear();
667         int nbPoints = facePoints.size()-1;
668         int nbFarPoints = 0;
669         for ( int i = 0; i < nbPoints; ++i )
670         {
671           const SMESH_NodeXYZ& n1 = facePoints[i];
672           const SMESH_NodeXYZ& n2 = facePoints[i+1];
673
674           size_t iP = intPoints.size();
675           intersectEdge( segCylinder, n1, n2, tol, intPoints );
676
677           // save edge index
678           if ( isCylinderOnFace )
679             for ( ; iP < intPoints.size(); ++iP )
680               intPoints[ iP ].myEdgeIndex = i;
681           else
682             for ( ; iP < intPoints.size(); ++iP )
683               intPoints[ iP ].myEdgeIndex = edgeIndex( i, n1, n2, face,
684                                                        intPoints[ iP ], faceNodes, tol );
685
686           nbFarPoints += ( segLine.SquareDistance( n1 ) > radius2 );
687         }
688
689         // feed startEdges
690         if ( nbFarPoints < nbPoints || !intPoints.empty() )
691           for ( size_t i = 1; i < faceNodes.size(); ++i )
692           {
693             const SMESH_NodeXYZ& n1 = faceNodes[i];
694             const SMESH_NodeXYZ& n2 = faceNodes[i-1];
695             isOut( n1, planeNormal, p[0].myIsOutPln );
696             isOut( n2, planeNormal, p[1].myIsOutPln );
697             if ( !isSegmentOut( p[0].myIsOutPln, p[1].myIsOutPln ))
698             {
699               startEdges.push_back( NLink( n1.Node(), n2.Node() ));
700             }
701           }
702
703         if ( intPoints.size() < 2 )
704           continue;
705
706         // classify intPoints by planes
707         for ( size_t i = 0; i < intPoints.size(); ++i )
708           isOut( intPoints[i].myNode, planeNormal, intPoints[i].myIsOutPln );
709
710         // cut the face
711
712         if ( intPoints.size() > 2 )
713           intPoints.push_back( intPoints[0] );
714
715         for ( size_t iE = 1; iE < intPoints.size(); ++iE ) // 2 <= intPoints.size() <= 5
716         {
717           if (( intPoints[iE].myIsOutPln[0] && intPoints[iE].myIsOutPln[1]   ) ||
718               ( isSegmentOut( intPoints[iE].myIsOutPln, intPoints[iE-1].myIsOutPln )))
719             continue; // intPoint is out of domain
720
721           // check if a cutting edge connecting two intPoints is on cylinder surface
722           if ( intPoints[iE].myEdgeIndex == intPoints[iE-1].myEdgeIndex )
723             continue; // on same edge
724           if ( intPoints[iE].myNode.Node() &&
725                intPoints[iE].myNode == intPoints[iE-1].myNode ) // coincide
726             continue;
727
728           gp_XYZ edegDir = intPoints[iE].myNode - intPoints[iE-1].myNode;
729
730           bool toCut; // = edegDir.SquareModulus() > tol * tol;
731           if ( intPoints.size() == 2 )
732             toCut = true;
733           else if ( isCylinderOnFace )
734             toCut = cylAxis.Direction().IsParallel( edegDir, angularTol );
735           else
736           {
737             SMESH_NodeXYZ nBetween;
738             int eInd = intPoints[iE-1].myEdgeIndex;
739             if ( eInd < 0 )
740               nBetween = facePoints[( 1 - (eInd-1)) % nbPoints ];
741             else
742               nBetween = faceNodes[( 1 + eInd ) % nbNodes ];
743             toCut = ( segLine.SquareDistance( nBetween ) > radius2 );
744           }
745           if ( !toCut )
746             continue;
747
748           // limit the edge by planes
749           if ( intPoints[iE].myIsOutPln[0] ||
750                intPoints[iE].myIsOutPln[1] )
751             cutOff( intPoints[iE], intPoints[iE-1],
752                     planeNormal[ intPoints[iE].myIsOutPln[1] ], tol );
753
754           if ( intPoints[iE-1].myIsOutPln[0] ||
755                intPoints[iE-1].myIsOutPln[1] )
756             cutOff( intPoints[iE-1], intPoints[iE],
757                     planeNormal[ intPoints[iE-1].myIsOutPln[1] ], tol );
758
759           gp_XYZ edegDirNew = intPoints[iE].myNode - intPoints[iE-1].myNode;
760           if ( edegDir * edegDirNew < 0 ||
761                edegDir.SquareModulus() < tol * tol )
762             continue; // fully cut off
763
764           segment->AddCutEdge( intPoints[iE], intPoints[iE-1], face );
765
766         }
767       }  // loop on faces sharing an edge
768
769       startEdges[0] = startEdges.back();
770       startEdges.pop_back();
771
772     } // loop on startEdges
773   } // loop on all input segments
774
775
776   // ----------------------------------------------------------
777   // If a plane fully cuts off edges of one side of a segment,
778   // it also may cut edges of adjacent segments
779   // ----------------------------------------------------------
780
781   for ( TSegmentIterator segIt( segmentPool ); segIt.more(); ) // loop on all segments
782   {
783     Segment* segment = const_cast< Segment* >( segIt.next() );
784     if ( segment->NbFreeEnds( tol ) >= 4 )
785       continue;
786
787     for ( int iE = 0; iE < 2; ++iE ) // loop on 2 segment ends
788     {
789       const SMDS_MeshNode* n1 = segment->Node( iE );
790       const SMDS_MeshNode* n2 = segment->Node( 1 - iE );
791       planeNormal[0] = segmentsOfNode( n1 ).Plane( segment );
792
793       bool isNeighborCut;
794       Segment* neighborSeg = segment;
795       do // check segments connected to the segment via n2
796       {
797         neighborSeg = nextSegment( neighborSeg, n2, segmentsOfNode );
798         if ( !neighborSeg )
799           break;
800
801         isNeighborCut = false;
802         for ( size_t iC = 0; iC < neighborSeg->myCuts.size(); ++iC ) // check cut edges
803         {
804           IntPoint* intPnt = &( neighborSeg->myCuts[iC].myIntPnt1 );
805           isOut( intPnt[0].myNode, planeNormal, intPnt[0].myIsOutPln, 1 );
806           isOut( intPnt[1].myNode, planeNormal, intPnt[1].myIsOutPln, 1 );
807           const Segment * closeSeg[2] = { 0, 0 };
808           if ( intPnt[0].myIsOutPln[0] )
809             closeSeg[0] = findTooCloseSegment( intPnt[0], 0.5 * theWidth - 1e-3*tol, tol,
810                                                segment, n1, segmentsOfNode );
811           if ( intPnt[1].myIsOutPln[0] )
812             closeSeg[1] = findTooCloseSegment( intPnt[1], 0.5 * theWidth - 1e-3*tol, tol,
813                                                segment, n1, segmentsOfNode );
814           int nbCut = bool( closeSeg[0] ) + bool( closeSeg[1] );
815           if ( nbCut == 0 )
816             continue;
817           isNeighborCut = true;
818           if ( nbCut == 2 ) // remove a cut
819           {
820             if ( iC < neighborSeg->myCuts.size() - 1 )
821               neighborSeg->myCuts[iC] = neighborSeg->myCuts.back();
822             neighborSeg->myCuts.pop_back();
823           }
824           else // shorten cuts of 1) neighborSeg and 2) closeSeg
825           {
826             // 1)
827             int iP = bool( closeSeg[1] );
828             gp_Lin      segLine( closeSeg[iP]->Ax1() );
829             gp_Ax3      cylAxis( segLine.Location(), segLine.Direction() );
830             gp_Cylinder cyl( cylAxis, 0.5 * theWidth );
831             intPoints.clear();
832             if ( intersectEdge( cyl, intPnt[iP].myNode, intPnt[1-iP].myNode, tol, intPoints ) &&
833                  intPoints[0].SquareDistance( intPnt[iP] ) > tol * tol )
834               intPnt[iP].myNode = intPoints[0].myNode;
835
836             // 2)
837             double minCutDist = theWidth;
838             gp_XYZ projection, closestProj;
839             int    iCut = -1;
840             for ( size_t iC2 = 0; iC2 < closeSeg[iP]->myCuts.size(); ++iC2 )
841             {
842               double cutDist = closeSeg[iP]->myCuts[iC2].SquareDistance( intPnt[iP].myNode,
843                                                                         projection );
844               if ( cutDist < minCutDist )
845               {
846                 closestProj = projection;
847                 minCutDist  = cutDist;
848                 iCut        = iC2;
849                 if ( minCutDist < tol * tol )
850                   break;
851               }
852             }
853             if ( iCut < 0 )
854               continue; // ???
855             double d1 = SMESH_MeshAlgos::GetDistance( neighborSeg->myEdge,
856                                                       closeSeg[iP]->myCuts[iCut][0].myNode );
857             double d2 = SMESH_MeshAlgos::GetDistance( neighborSeg->myEdge,
858                                                       closeSeg[iP]->myCuts[iCut][1].myNode );
859             int iP2 = ( d2 < d1 );
860             IntPoint& ip = const_cast< IntPoint& >( closeSeg[iP]->myCuts[iCut][iP2] );
861             ip = intPnt[iP];
862           }
863           // update myFreeEnds
864           neighborSeg->myFreeEnds.clear();
865           neighborSeg->NbFreeEnds( tol );
866         }
867       }
868       while ( isNeighborCut );
869     }
870   }
871
872   // -----------------------
873   // Cut faces by cut edges
874   // -----------------------
875
876   for ( TSegmentIterator segIt( segmentPool ); segIt.more(); ) // loop on all segments
877   {
878     Segment* segment = const_cast< Segment* >( segIt.next() );
879     for ( size_t iC = 0; iC < segment->myCuts.size(); ++iC )
880     {
881       Cut & cut = segment->myCuts[ iC ];
882       computeNormal( cut.myFace, faceNormals );
883       meshIntersector.Cut( cut.myFace,
884                            cut.myIntPnt1.myNode, cut.myIntPnt1.myEdgeIndex,
885                            cut.myIntPnt2.myNode, cut.myIntPnt2.myEdgeIndex );
886
887       Edge e = { cut.myIntPnt1.myNode.Node(), cut.myIntPnt2.myNode.Node(), 0 };
888       bndEdges.push_back( e );
889
890       findGroups( cut.myFace, theGroupsToUpdate, faceID2Groups, groupVec );
891     }
892   }
893
894   // -----------------------------------------
895   // Make cut at the end of group of segments
896   // -----------------------------------------
897
898   std::vector<const SMDS_MeshElement*> polySegments;
899
900   for ( TSegmentsOfNode::Iterator nSegsIt( segmentsOfNode ); nSegsIt.More(); nSegsIt.Next() )
901   {
902     const NodeData& noData = nSegsIt.Value();
903     if ( noData.mySegments.size() != 1 )
904       continue;
905
906     const Segment* segment = noData.mySegments[0];
907
908     // find two IntPoint's of cut edges to make a cut between
909     if ( segment->myFreeEnds.size() != 4 )
910       throw SALOME_Exception( "MakeSlot(): too short end edge?" );
911     std::multimap< double, const IntPoint* > dist2IntPntMap;
912     for ( size_t iE = 0; iE < segment->myFreeEnds.size(); ++iE )
913     {
914       const SMESH_NodeXYZ& n = segment->myFreeEnds[ iE ]->myNode;
915       double d = Abs( signedDist( n, noData.myPlane ));
916       dist2IntPntMap.insert( std::make_pair( d, segment->myFreeEnds[ iE ]));
917     }
918     std::multimap< double, const IntPoint* >::iterator d2ip = dist2IntPntMap.begin();
919     SMESH_MeshAlgos::PolySegment linkNodes;
920     linkNodes.myXYZ[0] = d2ip->second->myNode;
921     linkNodes.myXYZ[1] = (++d2ip)->second->myNode;
922     linkNodes.myVector = noData.myPlane.Direction() ^ (linkNodes.myXYZ[0] - linkNodes.myXYZ[1]);
923     linkNodes.myNode1[ 0 ] = linkNodes.myNode2[ 0 ] = 0;
924     linkNodes.myNode1[ 1 ] = linkNodes.myNode2[ 1 ] = 0;
925
926     // create segments connecting linkNodes
927     std::vector<const SMDS_MeshElement*> newSegments;
928     std::vector<const SMDS_MeshNode*>    newNodes;
929     SMESH_MeshAlgos::TListOfPolySegments polySegs(1, linkNodes);
930     SMESH_MeshAlgos::MakePolyLine( theMesh, polySegs, newSegments, newNodes,
931                                    /*group=*/0, faceSearcher.get() );
932     // cut faces by newSegments
933     intPoints.resize(2);
934     for ( size_t i = 0; i < newSegments.size(); ++i )
935     {
936       intPoints[0].myNode = edgeNodes[0] = newSegments[i]->GetNode(0);
937       intPoints[1].myNode = edgeNodes[1] = newSegments[i]->GetNode(1);
938
939       // find an underlying face
940       gp_XYZ                middle = 0.5 * ( intPoints[0].myNode + intPoints[1].myNode );
941       const SMDS_MeshElement* face = faceSearcher->FindClosestTo( middle, SMDSAbs_Face );
942
943       // find intersected edges of the face
944       int nbNodes = face->NbCornerNodes();
945       faceNodes.assign( face->begin_nodes(), face->end_nodes() );
946       faceNodes.resize( nbNodes + 1 );
947       faceNodes[ nbNodes ] = faceNodes[ 0 ];
948       for ( int iP = 0; iP < 2; ++iP )
949       {
950         intPoints[iP].myEdgeIndex = -1;
951         for ( int iN = 0; iN < nbNodes &&  intPoints[iP].myEdgeIndex < 0; ++iN )
952         {
953           if ( isOnEdge( faceNodes[iN], faceNodes[iN+1], intPoints[iP].myNode, tol ))
954             intPoints[iP].myEdgeIndex = iN;
955         }
956       }
957
958
959       // face cut
960       computeNormal( face, faceNormals );
961       meshIntersector.Cut( face,
962                            intPoints[0].myNode, intPoints[0].myEdgeIndex,
963                            intPoints[1].myNode, intPoints[1].myEdgeIndex );
964
965       Edge e = { intPoints[0].myNode.Node(), intPoints[1].myNode.Node(), 0 };
966       bndEdges.push_back( e );
967
968       findGroups( face, theGroupsToUpdate, faceID2Groups, groupVec );
969
970       // add cut points to an adjacent face at ends of poly-line
971       // if they fall onto face edges
972       if (( i == 0                       && intPoints[0].myEdgeIndex >= 0 ) ||
973           ( i == newSegments.size() - 1  && intPoints[1].myEdgeIndex >= 0 ))
974       {
975         for ( int iE = 0; iE < 2; ++iE ) // loop on ends of a new segment
976         {
977           if ( iE ? ( i != newSegments.size() - 1 ) : ( i != 0 ))
978             continue;
979           int iEdge = intPoints[ iE ].myEdgeIndex;
980           edgeNodes[0] = faceNodes[ iEdge ];
981           edgeNodes[1] = faceNodes[ iEdge+1 ];
982           theMesh->GetElementsByNodes( edgeNodes, faces, SMDSAbs_Face );
983           for ( size_t iF = 0; iF < faces.size(); ++iF )
984             if ( faces[iF] != face )
985             {
986               int iN1 = faces[iF]->GetNodeIndex( edgeNodes[0] );
987               int iN2 = faces[iF]->GetNodeIndex( edgeNodes[1] );
988               intPoints[ iE ].myEdgeIndex = Abs( iN1 - iN2 ) == 1 ? Min( iN1, iN2 ) : 2;
989               computeNormal( faces[iF], faceNormals );
990               meshIntersector.Cut( faces[iF],
991                                    intPoints[iE].myNode, intPoints[iE].myEdgeIndex,
992                                    intPoints[iE].myNode, intPoints[iE].myEdgeIndex );
993
994               findGroups( faces[iF], theGroupsToUpdate, faceID2Groups, groupVec );
995             }
996         }
997       }
998
999     } // loop on newSegments
1000
1001     polySegments.insert( polySegments.end(), newSegments.begin(), newSegments.end() );
1002
1003   } // loop on map of input segments
1004
1005   // actual mesh splitting
1006   TElemIntPairVec new2OldFaces;
1007   TNodeIntPairVec new2OldNodes;
1008   meshIntersector.MakeNewFaces( new2OldFaces, new2OldNodes, /*sign=*/1, /*optimize=*/true );
1009
1010   // add new faces to theGroupsToUpdate
1011   for ( size_t i = 0; i < new2OldFaces.size(); ++i )
1012   {
1013     const SMDS_MeshElement* newFace = new2OldFaces[i].first;
1014     const int             oldFaceID = new2OldFaces[i].second;
1015     if ( !newFace ) continue;
1016
1017     if ( TGroupVec* groups = const_cast< TGroupVec* >( faceID2Groups.Seek( oldFaceID )))
1018       for ( size_t iG = 0; iG < groups->size(); ++iG )
1019         (*groups)[ iG ]->Add( newFace );
1020   }
1021
1022   // remove poly-line edges
1023   for ( size_t i = 0; i < polySegments.size(); ++i )
1024   {
1025     edgeNodes[0] = polySegments[i]->GetNode(0);
1026     edgeNodes[1] = polySegments[i]->GetNode(1);
1027
1028     theMesh->RemoveFreeElement( polySegments[i] );
1029
1030     if ( edgeNodes[0]->NbInverseElements() == 0 )
1031       theMesh->RemoveNode( edgeNodes[0] );
1032     if ( edgeNodes[1]->NbInverseElements() == 0 )
1033       theMesh->RemoveNode( edgeNodes[1] );
1034   }
1035
1036   return bndEdges;
1037 }