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