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
  29 #include <IntAna_IntConicQuad.hxx>
  30 #include <IntAna_Quadric.hxx>
  31 #include <NCollection_DataMap.hxx>
  32 #include <NCollection_Map.hxx>
  33 #include <Precision.hxx>
  34 #include <gp_Ax1.hxx>
  35 #include <gp_Cylinder.hxx>
  36 #include <gp_Dir.hxx>
  37 #include <gp_Lin.hxx>
  38 #include <gp_Pln.hxx>
  39 #include <gp_Pnt.hxx>
  40 #include <gp_Vec.hxx>
  41
  42 #include <Utils_SALOME_Exception.hxx>
  43
  44 #include <boost/container/flat_set.hpp>
  45
  46 namespace
  47 {
  48   typedef SMESH_MeshAlgos::Edge TEdge;
  49
  50   //================================================================================
  51   //! point of intersection of a face edge with the cylinder
  52   struct IntPoint
  53   {
  54     SMESH_NodeXYZ myNode;        // point and a node
  55     int           myEdgeIndex;   // face edge index
  56     bool          myIsOutPln[2]; // isOut of two planes
  57   };
  58
  59   //================================================================================
  60   //! poly-line segment
  61   struct Segment
  62   {
  63     typedef boost::container::flat_set< const SMDS_MeshNode* > TNodeSet;
  64     //typedef std::list< TEdge > TEdgeList;
  65
  66     const SMDS_MeshElement* myEdge;
  67     TNodeSet                myEndNodes; // ends of cut edges
  68     //TEdgeList               myCutEdges[2];
  69
  70
  71     // return its axis
  72     gp_Ax1 Ax1( bool reversed = false ) const
  73     {
  74       SMESH_NodeXYZ n1 = myEdge->GetNode(  reversed );
  75       SMESH_NodeXYZ n2 = myEdge->GetNode( !reversed );
  76       return gp_Ax1( n1, gp_Dir( n2 - n1 ));
  77     }
  78     // return a node
  79     const SMDS_MeshNode* Node(int i) const
  80     {
  81       return myEdge->GetNode( i % 2 );
  82     }
  83     // store an intersection edge forming the slot border
  84     void AddEdge( TEdge& e, double tol )
  85     {
  86       const SMDS_MeshNode** nodes = & e._node1;
  87       for ( int i = 0; i < 2; ++i )
  88       {
  89         std::pair< TNodeSet::iterator, bool > nItAdded = myEndNodes.insert( nodes[ i ]);
  90         if ( !nItAdded.second )
  91           myEndNodes.erase( nItAdded.first );
  92       }
  93     }
  94     // { -- PREV version
  95     //   int i = myCutEdges[0].empty() ? 0 : 1;
  96     //   std::insert_iterator< TEdgeList > where = inserter( myCutEdges[i], myCutEdges[i].begin() );
  97
  98     //   //double minDist = 1e100;
  99     //   SMESH_NodeXYZ nNew[2] = { e._node1, e._node2 };
 100     //   int iNewMin = 0, iCurMin = 1;
 101     //   for ( i = 0; i < 2; ++i )
 102     //   {
 103     //     if ( myCutEdges[i].empty() )
 104     //       continue;
 105     //     SMESH_NodeXYZ nCur[2] = { myCutEdges[i].front()._node1,
 106     //                               myCutEdges[i].back()._node2 };
 107     //     for   ( int iN = 0; iN < 2; ++iN )
 108     //       for ( int iC = 0; iC < 2; ++iC )
 109     //       {
 110     //         if (( nCur[iC].Node() && nCur[iC] == nNew[iN] ) ||
 111     //             ( nCur[iC] - nNew[iN] ).SquareModulus() < tol * tol )
 112     //         {
 113     //           where   = inserter( myCutEdges[i], iC ? myCutEdges[i].end() : myCutEdges[i].begin() );
 114     //           iNewMin = iN;
 115     //           iCurMin = iC;
 116     //           //minDist = dist;
 117     //           iN = 2;
 118     //           break;
 119     //         }
 120     //       }
 121     //   }
 122     //   if ( iNewMin == iCurMin )
 123     //     std::swap( e._node1, e._node2 );
 124
 125     //   where = e;
 126     // }
 127     Segment( const SMDS_MeshElement* e = 0 ): myEdge(e) { myEndNodes.reserve( 4 ); }
 128   };
 129   typedef ObjectPoolIterator<Segment> TSegmentIterator;
 130
 131
 132   //================================================================================
 133   /*!
 134    * \brief Intersect a face edge given by its nodes with a cylinder.
 135    */
 136   //================================================================================
 137
 138   void intersectEdge( const gp_Cylinder&      cyl,
 139                       const SMESH_NodeXYZ&    n1,
 140                       const SMESH_NodeXYZ&    n2,
 141                       const double            tol,
 142                       std::vector< IntPoint >& intPoints )
 143   {
 144     gp_Lin line( gp_Ax1( n1, gp_Dir( n2 - n1 )));
 145     IntAna_IntConicQuad intersection( line, IntAna_Quadric( cyl ));
 146
 147     if ( !intersection.IsDone()     ||
 148          intersection.IsParallel()  ||
 149          intersection.IsInQuadric() ||
 150          intersection.NbPoints() == 0 )
 151       return;
 152
 153     gp_Vec edge( n1, n2 );
 154
 155     size_t oldNbPnts = intPoints.size();
 156     for ( int iP = 1; iP <= intersection.NbPoints(); ++iP )
 157     {
 158       const gp_Pnt& p = intersection.Point( iP );
 159
 160       gp_Vec n1p ( n1, p );
 161       const SMDS_MeshNode* n = 0;
 162
 163       double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
 164       if ( u <= 0. ) {
 165         if ( p.SquareDistance( n1 ) < tol * tol )
 166           n = n1.Node();
 167         else
 168           continue;
 169       }
 170       else if ( u >= 1. ) {
 171         if ( p.SquareDistance( n2 ) < tol * tol )
 172           n = n2.Node();
 173         else
 174           continue;
 175       }
 176       else {
 177         if      ( p.SquareDistance( n1 ) < tol * tol )
 178           n = n1.Node();
 179         else if ( p.SquareDistance( n2 ) < tol * tol )
 180           n = n2.Node();
 181       }
 182
 183       intPoints.push_back( IntPoint() );
 184       if ( n )
 185         intPoints.back().myNode.Set( n );
 186       else
 187         intPoints.back().myNode.SetCoord( p.X(),p.Y(),p.Z() );
 188     }
 189
 190     // set points order along an edge
 191     if ( intPoints.size() - oldNbPnts == 2 &&
 192          intersection.ParamOnConic( 1 ) > intersection.ParamOnConic( 2 ))
 193     {
 194       int i = intPoints.size() - 1;
 195       std::swap( intPoints[ i ], intPoints[ i - 1 ]);
 196     }
 197
 198     return;
 199   }
 200
 201   //================================================================================
 202   /*!
 203    * \brief Return signed distance between a point and a plane
 204    */
 205   //================================================================================
 206
 207   double signedDist( const gp_Pnt& p, const gp_Ax1& planeNormal )
 208   {
 209     const gp_Pnt& O = planeNormal.Location();
 210     gp_Vec Op( O, p );
 211     return Op * planeNormal.Direction();
 212   }
 213
 214   //================================================================================
 215   /*!
 216    * \brief Check if a point is outside a segment domain bound by two planes
 217    */
 218   //================================================================================
 219
 220   bool isOut( const gp_Pnt& p, const gp_Ax1* planeNormal, bool* isOutPtr )
 221   {
 222     isOutPtr[0] = isOutPtr[1] = false;
 223
 224     for ( int i = 0; i < 2; ++i )
 225     {
 226       isOutPtr[i] = ( signedDist( p, planeNormal[i] ) <= 0. );
 227     }
 228     return ( isOutPtr[0] && isOutPtr[1] );
 229   }
 230
 231   //================================================================================
 232   /*!
 233    * \brief Check if a segment between two points is outside a segment domain bound by two planes
 234    */
 235   //================================================================================
 236
 237   bool isSegmentOut( bool* isOutPtr1, bool* isOutPtr2 )
 238   {
 239     return (( isOutPtr1[0] && isOutPtr2[0] ) ||
 240             ( isOutPtr1[1] && isOutPtr2[1] ));
 241   }
 242
 243   //================================================================================
 244   /*!
 245    * \brief cut off ip1 from edge (ip1 - ip2) by a plane
 246    */
 247   //================================================================================
 248
 249   void cutOff( IntPoint & ip1, const IntPoint & ip2, const gp_Ax1& planeNormal, double tol )
 250   {
 251     gp_Lin lin( ip1.myNode, ( ip2.myNode - ip1.myNode ));
 252     gp_Pln pln( planeNormal.Location(), planeNormal.Direction() );
 253
 254     IntAna_IntConicQuad intersection( lin, pln, Precision::Angular/*Tolerance*/() );
 255     if ( intersection.IsDone()      &&
 256          !intersection.IsParallel()  &&
 257          !intersection.IsInQuadric() &&
 258          intersection.NbPoints() == 1 )
 259     {
 260       if ( intersection.Point( 1 ).SquareDistance( ip1.myNode ) > tol * tol )
 261       {
 262         static_cast< gp_XYZ& >( ip1.myNode ) = intersection.Point( 1 ).XYZ();
 263         ip1.myNode._node = 0;
 264         ip1.myEdgeIndex = -1;
 265       }
 266     }
 267   }
 268
 269   //================================================================================
 270   /*!
 271    * \brief Assure that face normal is computed in faceNormals vector
 272    */
 273   //================================================================================
 274
 275   const gp_XYZ& computeNormal( const SMDS_MeshElement* face,
 276                                std::vector< gp_XYZ >&  faceNormals )
 277   {
 278     bool toCompute;
 279     if ((int) faceNormals.size() <= face->GetID() )
 280     {
 281       toCompute = true;
 282       faceNormals.resize( face->GetID() + 1 );
 283     }
 284     else
 285     {
 286       toCompute = faceNormals[ face->GetID() ].SquareModulus() == 0.;
 287     }
 288     if ( toCompute )
 289       SMESH_MeshAlgos::FaceNormal( face, faceNormals[ face->GetID() ], /*normalized=*/false );
 290
 291     return faceNormals[ face->GetID() ];
 292   }
 293 }
 294
 295 //================================================================================
 296 /*!
 297  * \brief Create a slot of given width around given 1D elements lying on a triangle mesh.
 298  * The slot is consrtucted by cutting faces by cylindrical surfaces made around each segment.
 299  * \return Edges located at the slot boundary
 300  */
 301 //================================================================================
 302
 303 std::vector< SMESH_MeshAlgos::Edge >
 304 SMESH_MeshAlgos::MakeSlot( SMDS_ElemIteratorPtr theSegmentIt,
 305                            double               theWidth,
 306                            SMDS_Mesh*           theMesh)
 307 {
 308   std::vector< Edge > bndEdges;
 309
 310   if ( !theSegmentIt || !theSegmentIt->more() || !theMesh || theWidth == 0.)
 311     return bndEdges;
 312
 313   // put the input segments to a data map in order to be able finding neighboring ones
 314
 315   typedef std::vector< Segment* >                                                TSegmentVec;
 316   typedef NCollection_DataMap< const SMDS_MeshNode*, TSegmentVec, SMESH_Hasher > TSegmentsOfNode;
 317   TSegmentsOfNode segmentsOfNode;
 318   ObjectPool< Segment > segmentPool;
 319
 320   while( theSegmentIt->more() )
 321   {
 322     const SMDS_MeshElement* edge = theSegmentIt->next();
 323     if ( edge->GetType() != SMDSAbs_Edge )
 324       throw SALOME_Exception( "A segment is not a mesh edge");
 325
 326     Segment* segment = segmentPool.getNew();
 327     segment->myEdge = edge;
 328
 329     for ( SMDS_NodeIteratorPtr nIt = edge->nodeIterator(); nIt->more(); )
 330     {
 331       const SMDS_MeshNode* n = nIt->next();
 332       TSegmentVec* segVec = segmentsOfNode.ChangeSeek( n );
 333       if ( !segVec )
 334         segVec = segmentsOfNode.Bound( n, TSegmentVec() );
 335       segVec->reserve(2);
 336       segVec->push_back( segment );
 337     }
 338   }
 339
 340   // Cut the mesh around the segments
 341
 342   const double tol = Precision::Confusion();
 343   std::vector< gp_XYZ > faceNormals;
 344   SMESH_MeshAlgos::Intersector meshIntersector( theMesh, tol, faceNormals );
 345   std::unique_ptr< SMESH_ElementSearcher> faceSearcher;
 346
 347   std::vector< NLink > startEdges;
 348   std::vector< const SMDS_MeshNode* > faceNodes(4), edgeNodes(2);
 349   std::vector<const SMDS_MeshElement *> faces(2);
 350   NCollection_Map<const SMDS_MeshElement*, SMESH_Hasher > checkedFaces;
 351   std::vector< IntPoint > intPoints, p(2);
 352   std::vector< SMESH_NodeXYZ > facePoints(4);
 353   std::vector< Intersector::TFace > cutFacePoints;
 354
 355   std::vector< gp_Ax1 > planeNormalVec(2);
 356   gp_Ax1 * planeNormal = & planeNormalVec[0];
 357
 358   for ( TSegmentIterator segIt( segmentPool ); segIt.more(); ) // loop on all segments
 359   {
 360     Segment* segment = const_cast< Segment* >( segIt.next() );
 361
 362     gp_Lin      segLine( segment->Ax1() );
 363     gp_Ax3      cylAxis( segLine.Location(), segLine.Direction() );
 364     gp_Cylinder segCylinder( cylAxis, 0.5 * theWidth );
 365     double      radius2( segCylinder.Radius() * segCylinder.Radius() );
 366
 367     // get normals of planes separating domains of neighboring segments
 368     for ( int i = 0; i < 2; ++i ) // loop on 2 segment ends
 369     {
 370       planeNormal[i] = segment->Ax1(i);
 371
 372       const SMDS_MeshNode*    n = segment->Node( i );
 373       const TSegmentVec& segVec = segmentsOfNode( n );
 374       for ( size_t iS = 0; iS < segVec.size(); ++iS )
 375       {
 376         if ( segVec[iS] == segment )
 377           continue;
 378
 379         gp_Ax1 axis2 = segVec[iS]->Ax1();
 380         if ( n != segVec[iS]->Node( 1 ))
 381           axis2.Reverse(); // along a wire
 382
 383         planeNormal[i].SetDirection( planeNormal[i].Direction().XYZ() + axis2.Direction().XYZ() );
 384       }
 385     }
 386
 387     // we explore faces around a segment starting from face edges;
 388     // initialize a list of starting edges
 389     startEdges.clear();
 390     {
 391       // get a face to start searching intersected faces from
 392       const SMDS_MeshNode*      n0 = segment->Node( 0 );
 393       SMDS_ElemIteratorPtr     fIt = n0->GetInverseElementIterator( SMDSAbs_Face );
 394       const SMDS_MeshElement* face = ( fIt->more() ) ? fIt->next() : 0;
 395       if ( !theMesh->Contains( face ))
 396       {
 397         if ( !faceSearcher )
 398           faceSearcher.reset( SMESH_MeshAlgos::GetElementSearcher( *theMesh ));
 399         face = faceSearcher->FindClosestTo( SMESH_NodeXYZ( n0 ), SMDSAbs_Face );
 400       }
 401       // collect face edges
 402       int nbNodes = face->NbCornerNodes();
 403       faceNodes.assign( face->begin_nodes(), face->end_nodes() );
 404       faceNodes.resize( nbNodes + 1 );
 405       faceNodes[ nbNodes ] = faceNodes[ 0 ];
 406       for ( int i = 0; i < nbNodes; ++i )
 407         startEdges.push_back( NLink( faceNodes[i], faceNodes[i+1] ));
 408     }
 409
 410     // intersect faces located around a segment
 411     checkedFaces.Clear();
 412     while ( !startEdges.empty() )
 413     {
 414       edgeNodes[0] = startEdges[0].first;
 415       edgeNodes[1] = startEdges[0].second;
 416
 417       theMesh->GetElementsByNodes( edgeNodes, faces, SMDSAbs_Face );
 418       for ( size_t iF = 0; iF < faces.size(); ++iF ) // loop on faces sharing a start edge
 419       {
 420         const SMDS_MeshElement* face = faces[iF];
 421         if ( !checkedFaces.Add( face ))
 422           continue;
 423
 424         int nbNodes = face->NbCornerNodes();
 425         if ( nbNodes != 3 )
 426           throw SALOME_Exception( "MakeSlot() accepts triangles only" );
 427         facePoints.assign( face->begin_nodes(), face->end_nodes() );
 428         facePoints.resize( nbNodes + 1 );
 429         facePoints[ nbNodes ] = facePoints[ 0 ];
 430
 431         // check if cylinder axis || face
 432         const gp_XYZ& faceNorm = computeNormal( face, faceNormals );
 433         bool isCylinderOnFace  = ( Abs( faceNorm * cylAxis.Direction().XYZ() ) < tol );
 434
 435         if ( !isCylinderOnFace )
 436         {
 437           if ( Intersector::CutByPlanes( face, planeNormalVec, tol, cutFacePoints ))
 438             continue; // whole face cut off
 439           facePoints.swap( cutFacePoints[0] );
 440           facePoints.push_back( facePoints[0] );
 441         }
 442
 443         // find intersection points on face edges
 444         intPoints.clear();
 445         int nbPoints = facePoints.size()-1;
 446         int nbFarPoints = 0;
 447         for ( int i = 0; i < nbPoints; ++i )
 448         {
 449           const SMESH_NodeXYZ& n1 = facePoints[i];
 450           const SMESH_NodeXYZ& n2 = facePoints[i+1];
 451
 452           size_t iP = intPoints.size();
 453           intersectEdge( segCylinder, n1, n2, tol, intPoints );
 454
 455           // save edge index
 456           if ( isCylinderOnFace )
 457             for ( ; iP < intPoints.size(); ++iP )
 458               intPoints[ iP ].myEdgeIndex = i;
 459           else
 460             for ( ; iP < intPoints.size(); ++iP )
 461               if ( n1.Node() && n2.Node() )
 462                 intPoints[ iP ].myEdgeIndex = face->GetNodeIndex( n1.Node() );
 463               else
 464                 intPoints[ iP ].myEdgeIndex = -(i+1);
 465
 466           nbFarPoints += ( segLine.SquareDistance( n1 ) > radius2 );
 467         }
 468
 469         // feed startEdges
 470         if ( nbFarPoints < nbPoints || !intPoints.empty() )
 471           for ( int i = 0; i < nbPoints; ++i )
 472           {
 473             const SMESH_NodeXYZ& n1 = facePoints[i];
 474             const SMESH_NodeXYZ& n2 = facePoints[i+1];
 475             if ( n1.Node() && n2.Node() )
 476             {
 477               isOut( n1, planeNormal, p[0].myIsOutPln );
 478               isOut( n2, planeNormal, p[1].myIsOutPln );
 479               if ( !isSegmentOut( p[0].myIsOutPln, p[1].myIsOutPln ))
 480               {
 481                 startEdges.push_back( NLink( n1.Node(), n2.Node() ));
 482               }
 483             }
 484           }
 485
 486         if ( intPoints.size() < 2 )
 487           continue;
 488
 489         // classify intPoints by planes
 490         for ( size_t i = 0; i < intPoints.size(); ++i )
 491           isOut( intPoints[i].myNode, planeNormal, intPoints[i].myIsOutPln );
 492
 493         // cut the face
 494
 495         if ( intPoints.size() > 2 )
 496           intPoints.push_back( intPoints[0] );
 497
 498         for ( size_t iE = 1; iE < intPoints.size(); ++iE ) // 2 <= intPoints.size() <= 5
 499         {
 500           if (( intPoints[iE].myIsOutPln[0] && intPoints[iE].myIsOutPln[1]   ) ||
 501               ( isSegmentOut( intPoints[iE].myIsOutPln, intPoints[iE-1].myIsOutPln )))
 502             continue; // intPoint is out of domain
 503
 504           // check if a cutting edge connecting two intPoints is on cylinder surface
 505           if ( intPoints[iE].myEdgeIndex == intPoints[iE-1].myEdgeIndex )
 506             continue; // on same edge
 507           if ( intPoints[iE].myNode.Node() &&
 508                intPoints[iE].myNode == intPoints[iE-1].myNode ) // coincide
 509             continue;
 510
 511           gp_XYZ edegDir = intPoints[iE].myNode - intPoints[iE-1].myNode;
 512
 513           bool toCut; // = edegDir.SquareModulus() > tol * tol;
 514           if ( intPoints.size() == 2 )
 515             toCut = true;
 516           else if ( isCylinderOnFace )
 517             toCut = cylAxis.Direction().IsParallel( edegDir, tol );
 518           else
 519           {
 520             SMESH_NodeXYZ nBetween;
 521             int eInd = intPoints[iE-1].myEdgeIndex;
 522             if ( eInd < 0 )
 523               nBetween = facePoints[( 1 - (eInd-1)) % nbPoints ];
 524             else
 525               nBetween = faceNodes[( 1 + eInd ) % nbNodes ];
 526             toCut = ( segLine.SquareDistance( nBetween ) > radius2 );
 527           }
 528           if ( !toCut )
 529             continue;
 530
 531           // limit the edge by planes
 532           if ( intPoints[iE].myIsOutPln[0] ||
 533                intPoints[iE].myIsOutPln[1] )
 534             cutOff( intPoints[iE], intPoints[iE-1],
 535                     planeNormal[ intPoints[iE].myIsOutPln[1] ], tol );
 536
 537           if ( intPoints[iE-1].myIsOutPln[0] ||
 538                intPoints[iE-1].myIsOutPln[1] )
 539             cutOff( intPoints[iE-1], intPoints[iE],
 540                     planeNormal[ intPoints[iE-1].myIsOutPln[1] ], tol );
 541
 542           edegDir = intPoints[iE].myNode - intPoints[iE-1].myNode;
 543           if ( edegDir.SquareModulus() < tol * tol )
 544             continue; // fully cut off
 545
 546           // face cut
 547           meshIntersector.Cut( face,
 548                                intPoints[iE-1].myNode, intPoints[iE-1].myEdgeIndex,
 549                                intPoints[iE  ].myNode, intPoints[iE  ].myEdgeIndex );
 550
 551           Edge e = { intPoints[iE].myNode.Node(), intPoints[iE-1].myNode.Node(), 0 };
 552           segment->AddEdge( e, tol );
 553           bndEdges.push_back( e );
 554         }
 555       }  // loop on faces sharing an edge
 556
 557       startEdges[0] = startEdges.back();
 558       startEdges.pop_back();
 559
 560     } // loop on startEdges
 561   } // loop on all input segments
 562
 563
 564   // Make cut at the end of group of segments
 565
 566   std::vector<const SMDS_MeshElement*> polySegments;
 567
 568   for ( TSegmentsOfNode::Iterator nSegsIt( segmentsOfNode ); nSegsIt.More(); nSegsIt.Next() )
 569   {
 570     const TSegmentVec& segVec = nSegsIt.Value();
 571     if ( segVec.size() != 1 )
 572       continue;
 573
 574     const Segment*       segment = segVec[0];
 575     const SMDS_MeshNode* segNode = nSegsIt.Key();
 576
 577     // find two end nodes of cut edges to make a cut between
 578     if ( segment->myEndNodes.size() != 4 )
 579       throw SALOME_Exception( "MakeSlot(): too short end edge?" );
 580     SMESH_MeshAlgos::PolySegment linkNodes;
 581     gp_Ax1 planeNorm = segment->Ax1( segNode != segment->Node(0) );
 582     double minDist[2] = { 1e100, 1e100 };
 583     Segment::TNodeSet::const_iterator nIt = segment->myEndNodes.begin();
 584     for ( ; nIt != segment->myEndNodes.end(); ++nIt )
 585     {
 586       SMESH_NodeXYZ n = *nIt;
 587       double d = Abs( signedDist( n, planeNorm ));
 588       double diff1 = minDist[0] - d, diff2 = minDist[1] - d;
 589       int i;
 590       if ( diff1 > 0 && diff2 > 0 )
 591       {
 592         i = ( diff1 < diff2 );
 593       }
 594       else if ( diff1 > 0 )
 595       {
 596         i = 0;
 597       }
 598       else if ( diff2 > 0 )
 599       {
 600         i = 1;
 601       }
 602       else
 603       {
 604         continue;
 605       }
 606       linkNodes.myXYZ[ i ] = n;
 607       minDist        [ i ] = d;
 608     }
 609     // for ( int iSide = 0; iSide < 2; ++iSide )
 610     // {
 611     //   if ( segment->myCutEdges[ iSide ].empty() )
 612     //     throw SALOME_Exception( "MakeSlot(): too short end edge?" );
 613     //   SMESH_NodeXYZ n1 = segment->myCutEdges[ iSide ].front()._node1;
 614     //   SMESH_NodeXYZ n2 = segment->myCutEdges[ iSide ].back ()._node2;
 615     //   double d1 = Abs( signedDist( n1, planeNorm ));
 616     //   double d2 = Abs( signedDist( n2, planeNorm ));
 617     //   linkNodes.myXYZ  [ iSide ] = ( d1 < d2 ) ? n1 : n2;
 618     //   linkNodes.myNode1[ iSide ] = linkNodes.myNode2[ iSide ] = 0;
 619     // }
 620     linkNodes.myVector = planeNorm.Direction() ^ (linkNodes.myXYZ[0] - linkNodes.myXYZ[1]);
 621     linkNodes.myNode1[ 0 ] = linkNodes.myNode2[ 0 ] = 0;
 622     linkNodes.myNode1[ 1 ] = linkNodes.myNode2[ 1 ] = 0;
 623
 624     // create segments connecting linkNodes
 625     std::vector<const SMDS_MeshElement*> newSegments;
 626     std::vector<const SMDS_MeshNode*>    newNodes;
 627     SMESH_MeshAlgos::TListOfPolySegments polySegs(1, linkNodes);
 628     SMESH_MeshAlgos::MakePolyLine( theMesh, polySegs, newSegments, newNodes,
 629                                    /*group=*/0, faceSearcher.get() );
 630     // cut faces by newSegments
 631     intPoints.resize(2);
 632     for ( size_t i = 0; i < newSegments.size(); ++i )
 633     {
 634       intPoints[0].myNode = edgeNodes[0] = newSegments[i]->GetNode(0);
 635       intPoints[1].myNode = edgeNodes[1] = newSegments[i]->GetNode(1);
 636
 637       // find an underlying face
 638       gp_XYZ                middle = 0.5 * ( intPoints[0].myNode + intPoints[1].myNode );
 639       const SMDS_MeshElement* face = faceSearcher->FindClosestTo( middle, SMDSAbs_Face );
 640
 641       // find intersected edges of the face
 642       int nbNodes = face->NbCornerNodes();
 643       faceNodes.assign( face->begin_nodes(), face->end_nodes() );
 644       faceNodes.resize( nbNodes + 1 );
 645       faceNodes[ nbNodes ] = faceNodes[ 0 ];
 646       for ( int iP = 0; iP < 2; ++iP )
 647       {
 648         intPoints[iP].myEdgeIndex = -1;
 649         for ( int iN = 0; iN < nbNodes &&  intPoints[iP].myEdgeIndex < 0; ++iN )
 650         {
 651           SMDS_LinearEdge edge( faceNodes[iN], faceNodes[iN+1] );
 652           if ( SMESH_MeshAlgos::GetDistance( &edge, intPoints[iP].myNode) < tol )
 653             intPoints[iP].myEdgeIndex = iN;
 654         }
 655       }
 656
 657
 658       // face cut
 659       computeNormal( face, faceNormals );
 660       meshIntersector.Cut( face,
 661                            intPoints[0].myNode, intPoints[0].myEdgeIndex,
 662                            intPoints[1].myNode, intPoints[1].myEdgeIndex );
 663
 664       Edge e = { intPoints[0].myNode.Node(), intPoints[1].myNode.Node(), 0 };
 665       bndEdges.push_back( e );
 666
 667       // add cut points to an adjacent face at ends of poly-line
 668       // if they fall onto face edges
 669       if (( i == 0                       && intPoints[0].myEdgeIndex >= 0 ) ||
 670           ( i == newSegments.size() - 1  && intPoints[1].myEdgeIndex >= 0 ))
 671       {
 672         for ( int iE = 0; iE < 2; ++iE ) // loop on ends of a new segment
 673         {
 674           if ( iE ? ( i != newSegments.size() - 1 ) : ( i != 0 ))
 675             continue;
 676           int iEdge = intPoints[ iE ].myEdgeIndex;
 677           edgeNodes[0] = faceNodes[ iEdge ];
 678           edgeNodes[1] = faceNodes[ iEdge+1 ];
 679           theMesh->GetElementsByNodes( edgeNodes, faces, SMDSAbs_Face );
 680           for ( size_t iF = 0; iF < faces.size(); ++iF )
 681             if ( faces[iF] != face )
 682             {
 683               int iN1 = faces[iF]->GetNodeIndex( edgeNodes[0] );
 684               int iN2 = faces[iF]->GetNodeIndex( edgeNodes[1] );
 685               intPoints[ iE ].myEdgeIndex = Abs( iN1 - iN2 ) == 1 ? Min( iN1, iN2 ) : 2;
 686               computeNormal( faces[iF], faceNormals );
 687               meshIntersector.Cut( faces[iF],
 688                                    intPoints[iE].myNode, intPoints[iE].myEdgeIndex,
 689                                    intPoints[iE].myNode, intPoints[iE].myEdgeIndex );
 690             }
 691         }
 692       }
 693
 694     } // loop on newSegments
 695
 696     polySegments.insert( polySegments.end(), newSegments.begin(), newSegments.end() );
 697
 698   } // loop on map of input segments
 699
 700   // actual mesh splitting
 701   TElemIntPairVec new2OldFaces;
 702   TNodeIntPairVec new2OldNodes;
 703   meshIntersector.MakeNewFaces( new2OldFaces, new2OldNodes, /*sign=*/1, /*optimize=*/true );
 704
 705   // remove poly-line edges
 706   for ( size_t i = 0; i < polySegments.size(); ++i )
 707   {
 708     edgeNodes[0] = polySegments[i]->GetNode(0);
 709     edgeNodes[1] = polySegments[i]->GetNode(1);
 710
 711     theMesh->RemoveFreeElement( polySegments[i] );
 712
 713     if ( edgeNodes[0]->NbInverseElements() == 0 )
 714       theMesh->RemoveNode( edgeNodes[0] );
 715     if ( edgeNodes[1]->NbInverseElements() == 0 )
 716       theMesh->RemoveNode( edgeNodes[1] );
 717   }
 718
 719   return bndEdges;
 720 }