src/SMESHUtils/SMESH_Slot.cxx

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