src/SMESHUtils/SMESH_FreeBorders.cxx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : SMESH_FreeBorders.cxx
  20 // Created   : Tue Sep  8 17:08:39 2015
  21 // Author    : Edward AGAPOV (eap)
  22
  23 //================================================================================
  24 // Implementation of SMESH_MeshAlgos::FindCoincidentFreeBorders()
  25 //================================================================================
  26
  27 #include "SMESH_MeshAlgos.hxx"
  28
  29 #include "SMDS_LinearEdge.hxx"
  30 #include "SMDS_Mesh.hxx"
  31 #include "SMDS_SetIterator.hxx"
  32
  33 #include <algorithm>
  34 #include <limits>
  35 #include <set>
  36 #include <vector>
  37
  38 #include <NCollection_DataMap.hxx>
  39 #include <Precision.hxx>
  40 #include <gp_Pnt.hxx>
  41
  42 using namespace SMESH_MeshAlgos;
  43
  44 namespace
  45 {
  46   struct BEdge;
  47
  48   /*!
  49    * \brief Node on a free border
  50    */
  51   struct BNode : public SMESH_TNodeXYZ
  52   {
  53     mutable std::vector< BEdge* > myLinkedEdges;
  54     mutable std::vector< std::pair < BEdge*, double > > myCloseEdges; // edge & U
  55
  56     BNode(const SMDS_MeshNode * node): SMESH_TNodeXYZ( node ) {}
  57     const SMDS_MeshNode * Node() const { return _node; }
  58     void   AddLinked( BEdge* e ) const;
  59     void   AddClose ( const BEdge* e, double u ) const;
  60     BEdge* GetCloseEdge( size_t i ) const { return myCloseEdges[i].first; }
  61     double GetCloseU( size_t i ) const { return myCloseEdges[i].second; }
  62     BEdge* GetCloseEdgeOfBorder( int borderID, double * u = 0 ) const;
  63     bool   HasCloseEdgeWithNode( const BNode* n ) const;
  64     bool   IsCloseEdge( const BEdge*, double * u = 0 ) const;
  65     bool operator<(const BNode& other) const { return Node()->GetID() < other.Node()->GetID(); }
  66     double SquareDistance(const BNode& e2) const { return ( e2 - *this ).SquareModulus(); }
  67   };
  68   /*!
  69    * \brief Edge of a free border
  70    */
  71   struct BEdge : public SMDS_LinearEdge
  72   {
  73     const BNode*            myBNode1;
  74     const BNode*            myBNode2;
  75     int                     myBorderID;
  76     int                     myID; // within a border
  77     BEdge*                  myPrev;
  78     BEdge*                  myNext;
  79     const SMDS_MeshElement* myFace;
  80     std::set< int >         myCloseBorders;
  81     int                     myInGroup;
  82
  83     BEdge():SMDS_LinearEdge( 0, 0 ), myBorderID(-1), myID(-1), myPrev(0), myNext(0), myInGroup(-1) {}
  84
  85     void Set( const BNode *           node1,
  86               const BNode *           node2,
  87               const SMDS_MeshElement* face,
  88               const int               ID)
  89     {
  90       myBNode1   = node1;
  91       myBNode2   = node2;
  92       myNodes[0] = node1->Node();
  93       myNodes[1] = node2->Node();
  94       myFace     = face;
  95       setId( ID ); // mesh element ID
  96     }
  97     bool IsInGroup() const
  98     {
  99       return myInGroup >= 0;
 100     }
 101     bool Contains( const BNode* n ) const
 102     {
 103       return ( n == myBNode1 || n == myBNode2 );
 104     }
 105     void AddLinked( BEdge* e )
 106     {
 107       if ( e->Contains( myBNode1 )) myPrev = e;
 108       else                          myNext = e;
 109     }
 110     void RemoveLinked( BEdge* e )
 111     {
 112       if ( myPrev == e ) myPrev = 0;
 113       if ( myNext == e ) myNext = 0;
 114     }
 115     void Reverse()
 116     {
 117       std::swap( myBNode1, myBNode2 );
 118       myNodes[0] = myBNode1->Node();
 119       myNodes[1] = myBNode2->Node();
 120     }
 121     void Orient()
 122     {
 123       if (( myPrev && !myPrev->Contains( myBNode1 )) ||
 124           ( myNext && !myNext->Contains( myBNode2 )))
 125         std::swap( myPrev, myNext );
 126       if ( myPrev && myPrev->myBNode2 != myBNode1 ) myPrev->Reverse();
 127       if ( myNext && myNext->myBNode1 != myBNode2 ) myNext->Reverse();
 128     }
 129     void SetID( int id )
 130     {
 131       if ( myID < 0 )
 132       {
 133         myID = id;
 134         if ( myNext )
 135           myNext->SetID( id + 1 );
 136       }
 137     }
 138     //================================================================================
 139     /*!
 140      * \brief Checks if a point is closer to this BEdge than tol
 141      */
 142     //================================================================================
 143
 144     bool IsOut( const gp_XYZ& point, const double tol, double& u ) const
 145     {
 146       gp_XYZ  me = *myBNode2 - *myBNode1;
 147       gp_XYZ n1p = point     - *myBNode1;
 148       u = ( me * n1p ) / me.SquareModulus(); // param [0,1] on this
 149       if ( u < 0. ) return ( n1p.SquareModulus() > tol * tol );
 150       if ( u > 1. ) return ( ( point - *myBNode2 ).SquareModulus() > tol * tol );
 151
 152       gp_XYZ  proj = ( 1. - u ) * *myBNode1 + u * *myBNode2; // projection of the point on this
 153       double dist2 = ( point - proj ).SquareModulus();
 154       return ( dist2 > tol * tol );
 155     }
 156     //================================================================================
 157     /*!
 158      * \brief Checks if two BEdges can be considered as overlapping
 159      */
 160     //================================================================================
 161
 162     bool IsOverlappingProjection( const BEdge* toE, const double u, bool is1st ) const
 163     {
 164       // is1st shows which end of toE is projected on this at u
 165       double u2;
 166       const double eps = 0.1;
 167       if ( myBNode1->IsCloseEdge( toE, &u2 ) ||
 168            myBNode2->IsCloseEdge( toE, &u2 ))
 169         return (( 0 < u2 && u2 < 1 ) &&     // u2 is proj param of myBNode's on toE
 170                 ( Abs( u2 - int( !is1st )) > eps ));
 171
 172       const BNode* n = is1st ? toE->myBNode2 : toE->myBNode1;
 173       if ( this == n->GetCloseEdgeOfBorder( this->myBorderID, &u2 ))
 174         return Abs( u - u2 ) > eps;
 175       return false;
 176     }
 177     //================================================================================
 178     /*!
 179      * \brief Finds all neighbor BEdge's having the same close borders
 180      */
 181     //================================================================================
 182
 183     bool GetRangeOfSameCloseBorders(BEdge* eRange[2], const std::set< int >& bordIDs)
 184     {
 185       if ( this->myCloseBorders != bordIDs )
 186         return false;
 187
 188       if ( bordIDs.size() == 1 && bordIDs.count( myBorderID )) // border close to self
 189       {
 190         double u;
 191         eRange[0] = this;
 192         while ( eRange[0]->myBNode1->GetCloseEdgeOfBorder( myBorderID, &u ))
 193         {
 194           if ( eRange[0]->myPrev == this || u < 0 || u > 1 )
 195             break;
 196           eRange[0] = eRange[0]->myPrev;
 197         }
 198         eRange[1] = this;
 199         while ( eRange[1]->myBNode2->GetCloseEdgeOfBorder( myBorderID, &u ))
 200         {
 201           if ( eRange[1]->myNext == this || u < 0 || u > 1 )
 202             break;
 203           eRange[1] = eRange[1]->myNext;
 204         }
 205       }
 206       else
 207       {
 208         eRange[0] = this;
 209         while ( eRange[0]->myPrev && eRange[0]->myPrev->myCloseBorders == bordIDs )
 210         {
 211           if ( eRange[0]->myPrev == this )
 212             break;
 213           eRange[0] = eRange[0]->myPrev;
 214         }
 215
 216         eRange[1] = this;
 217         if ( eRange[0]->myPrev != this ) // not closed border
 218           while ( eRange[1]->myNext && eRange[1]->myNext->myCloseBorders == bordIDs )
 219           {
 220             if ( eRange[1]->myNext == this )
 221               break;
 222             eRange[1] = eRange[1]->myNext;
 223           }
 224       }
 225
 226       if ( eRange[0] == eRange[1] )
 227       {
 228         std::set<int>::iterator closeBord = eRange[0]->myCloseBorders.begin();
 229         for ( ; closeBord != eRange[0]->myCloseBorders.end(); ++closeBord )
 230         {
 231           if ( BEdge* be = eRange[0]->myBNode1->GetCloseEdgeOfBorder( *closeBord ))
 232             if ( be->myCloseBorders == eRange[0]->myCloseBorders )
 233               return true;
 234           if ( BEdge* be = eRange[0]->myBNode2->GetCloseEdgeOfBorder( *closeBord ))
 235             if ( be->myCloseBorders == eRange[0]->myCloseBorders )
 236               return true;
 237         }
 238         return false;
 239       }
 240       return true;
 241     }
 242   }; // class BEdge
 243
 244   //================================================================================
 245   /*!
 246    * \brief Checks if all border parts include the whole closed border, and if so
 247    *        returns \c true and choose starting BEdge's with most coincident nodes
 248    */
 249   //================================================================================
 250
 251   bool chooseStartOfClosedBorders( std::vector< BEdge* >& ranges ) // PAL23078#c21002
 252   {
 253     bool allClosed = true;
 254     for ( size_t iR = 1; iR < ranges.size() && allClosed; iR += 2 )
 255       allClosed = ( ranges[ iR-1 ]->myPrev == ranges[ iR ] );
 256     if ( !allClosed )
 257       return allClosed;
 258
 259     double u, minDiff = Precision::Infinite();
 260     std::vector< BEdge* > start( ranges.size() / 2 );
 261     BEdge* range0 = start[0] = ranges[0];
 262     do
 263     {
 264       double maxDiffU = 0;
 265       double  maxDiff = 0;
 266       for ( size_t iR = 3; iR < ranges.size(); iR += 2 )
 267       {
 268         int borderID = ranges[iR]->myBorderID;
 269         if ( BEdge* e = start[0]->myBNode1->GetCloseEdgeOfBorder( borderID, & u ))
 270         {
 271           start[ iR / 2 ] = e;
 272           double diffU = Min( Abs( u ), Abs( 1.-u ));
 273           double  diff = e->myBNode1->SquareDistance( *e->myBNode2 ) * diffU * diffU;
 274           maxDiffU = Max( diffU, maxDiffU );
 275           maxDiff  = Max( diff,  maxDiff );
 276         }
 277       }
 278       if ( maxDiff < minDiff )
 279       {
 280         minDiff = maxDiff;
 281         for ( size_t iR = 1; iR < ranges.size(); iR += 2 )
 282         {
 283           ranges[ iR-1 ] = start[ iR/2 ];
 284           ranges[ iR   ] = ranges[ iR-1]->myPrev;
 285         }
 286       }
 287       if ( maxDiffU < 1e-6 )
 288         break;
 289       start[0] = start[0]->myNext;
 290     }
 291     while ( start[0] != range0 );
 292
 293     return allClosed;
 294   }
 295
 296   //================================================================================
 297   /*!
 298    * \brief Tries to include neighbor BEdge's into a border part
 299    */
 300   //================================================================================
 301
 302   void extendPart( BEdge* & e1, BEdge* & e2, const std::set< int >& bordIDs, int groupID )
 303   {
 304     if (( e1->myPrev == e2 ) ||
 305         ( e1 == e2 && e1->myPrev && e1->myPrev->myInGroup == groupID ))
 306       return; // full free border already
 307
 308     double u;
 309     BEdge* be;
 310     std::set<int>::const_iterator bord;
 311     if ( e1->myPrev )
 312     {
 313       for ( bord = bordIDs.begin(); bord != bordIDs.end(); ++bord )
 314         if ((( be = e1->myBNode1->GetCloseEdgeOfBorder( *bord, &u ))) &&
 315             (  be->myInGroup == groupID ) &&
 316             (  0 < u && u < 1 ) &&
 317             (  be->IsOverlappingProjection( e1->myPrev, u, false )))
 318         {
 319           e1 = e1->myPrev;
 320           break;
 321         }
 322       if ( bord == bordIDs.end() && // not extended
 323            e1->myBNode1->HasCloseEdgeWithNode( e1->myPrev->myBNode1 ))
 324       {
 325         e1 = e1->myPrev;
 326       }
 327       e1->myInGroup = groupID;
 328     }
 329     if ( e2->myNext )
 330     {
 331       for ( bord = bordIDs.begin(); bord != bordIDs.end(); ++bord )
 332         if ((( be = e2->myBNode2->GetCloseEdgeOfBorder( *bord, &u ))) &&
 333             (  be->myInGroup == groupID ) &&
 334             (  0 < u && u < 1 ) &&
 335             (  be->IsOverlappingProjection( e2->myNext, u, true )))
 336         {
 337           e2 = e2->myNext;
 338           break;
 339         }
 340       if ( bord == bordIDs.end() && // not extended
 341            e2->myBNode2->HasCloseEdgeWithNode( e2->myNext->myBNode2 ))
 342       {
 343         e2 = e2->myNext;
 344       }
 345       e2->myInGroup = groupID;
 346     }
 347   }
 348
 349   //================================================================================
 350   /*!
 351    * \brief Connect BEdge's incident at this node
 352    */
 353   //================================================================================
 354
 355   void BNode::AddLinked( BEdge* e ) const
 356   {
 357     myLinkedEdges.reserve(2);
 358     myLinkedEdges.push_back( e );
 359     if ( myLinkedEdges.size() < 2 ) return;
 360
 361     if ( myLinkedEdges.size() == 2 )
 362     {
 363       myLinkedEdges[0]->AddLinked( myLinkedEdges[1] );
 364       myLinkedEdges[1]->AddLinked( myLinkedEdges[0] );
 365     }
 366     else
 367     {
 368       for ( size_t i = 0; i < myLinkedEdges.size(); ++i )
 369         for ( size_t j = 0; j < myLinkedEdges.size(); ++j )
 370           if ( i != j )
 371             myLinkedEdges[i]->RemoveLinked( myLinkedEdges[j] );
 372     }
 373   }
 374   void BNode::AddClose ( const BEdge* e, double u ) const
 375   {
 376     if ( ! e->Contains( this ))
 377       myCloseEdges.push_back( std::make_pair( const_cast< BEdge* >( e ), u ));
 378   }
 379   BEdge* BNode::GetCloseEdgeOfBorder( int borderID, double * uPtr ) const
 380   {
 381     BEdge* e = 0;
 382     double u = 0;
 383     for ( size_t i = 0; i < myCloseEdges.size(); ++i )
 384       if ( borderID == GetCloseEdge( i )->myBorderID )
 385       {
 386         if ( e && Abs( u - 0.5 ) < Abs( GetCloseU( i ) - 0.5 ))
 387           continue;
 388         u = GetCloseU( i );
 389         e = GetCloseEdge ( i );
 390       }
 391     if ( uPtr ) *uPtr = u;
 392     return e;
 393   }
 394   bool BNode::HasCloseEdgeWithNode( const BNode* n ) const
 395   {
 396     for ( size_t i = 0; i < myCloseEdges.size(); ++i )
 397       if ( GetCloseEdge( i )->Contains( n ) &&
 398            0 < GetCloseU( i ) && GetCloseU( i ) < 1 )
 399         return true;
 400     return false;
 401   }
 402   bool BNode::IsCloseEdge( const BEdge* e, double * uPtr ) const
 403   {
 404     for ( size_t i = 0; i < myCloseEdges.size(); ++i )
 405       if ( e == GetCloseEdge( i ) )
 406       {
 407         if ( uPtr ) *uPtr = GetCloseU( i );
 408         return true;
 409       }
 410     return false;
 411   }
 412
 413   /// Accessor to SMDS_MeshElement* inherited by BEdge
 414   struct ElemAcess
 415   {
 416     static const SMDS_MeshElement* value( std::vector< BEdge >::const_iterator it)
 417     {
 418       return & (*it);
 419     }
 420   };
 421   /// Iterator over a vector of BEdge's
 422   static SMDS_ElemIteratorPtr getElemIterator( const std::vector< BEdge > & bedges )
 423   {
 424     typedef SMDS_SetIterator
 425       < const SMDS_MeshElement*, std::vector< BEdge >::const_iterator, ElemAcess > BEIter;
 426     return SMDS_ElemIteratorPtr( new BEIter( bedges.begin(), bedges.end() ));
 427   }
 428
 429 } // namespace
 430
 431 //================================================================================
 432 /*
 433  * Returns groups of TFreeBorder's coincident within the given tolerance.
 434  * If the tolerance <= 0.0 then one tenth of an average size of elements adjacent
 435  * to free borders being compared is used.
 436  */
 437 //================================================================================
 438
 439 void SMESH_MeshAlgos::FindCoincidentFreeBorders(SMDS_Mesh&              mesh,
 440                                                 double                  tolerance,
 441                                                 CoincidentFreeBorders & foundFreeBordes)
 442 {
 443   // find free links
 444   typedef NCollection_DataMap<SMESH_TLink, const SMDS_MeshElement*, SMESH_TLink > TLink2FaceMap;
 445   TLink2FaceMap linkMap;
 446   int nbSharedLinks = 0;
 447   SMDS_FaceIteratorPtr faceIt = mesh.facesIterator();
 448   while ( faceIt->more() )
 449   {
 450     const SMDS_MeshElement* face = faceIt->next();
 451     if ( !face ) continue;
 452
 453     const SMDS_MeshNode*     n0 = face->GetNode( face->NbNodes() - 1 );
 454     SMDS_NodeIteratorPtr nodeIt = face->interlacedNodesIterator();
 455     while ( nodeIt->more() )
 456     {
 457       const SMDS_MeshNode* n1 = nodeIt->next();
 458       SMESH_TLink link( n0, n1 );
 459       if ( const SMDS_MeshElement** faceInMap = linkMap.ChangeSeek( link ))
 460       {
 461         nbSharedLinks += bool( *faceInMap );
 462         *faceInMap = 0;
 463       }
 464       else
 465       {
 466         linkMap.Bind( link, face );
 467       }
 468       n0 = n1;
 469     }
 470   }
 471   if ( linkMap.Extent() == nbSharedLinks )
 472     return;
 473
 474   // form free borders
 475   std::set   < BNode > bNodes;
 476   std::vector< BEdge > bEdges( linkMap.Extent() - nbSharedLinks );
 477
 478   TLink2FaceMap::Iterator linkIt( linkMap );
 479   for ( int iEdge = 0; linkIt.More(); linkIt.Next() )
 480   {
 481     if ( !linkIt.Value() ) continue;
 482     const SMESH_TLink & link = linkIt.Key();
 483     std::set< BNode >::iterator n1 = bNodes.insert( BNode( link.node1() )).first;
 484     std::set< BNode >::iterator n2 = bNodes.insert( BNode( link.node2() )).first;
 485     bEdges[ iEdge ].Set( &*n1, &*n2, linkIt.Value(), iEdge+1 );
 486     n1->AddLinked( & bEdges[ iEdge ] );
 487     n2->AddLinked( & bEdges[ iEdge ] );
 488     ++iEdge;
 489   }
 490   linkMap.Clear();
 491
 492   // assign IDs to borders
 493   std::vector< BEdge* > borders; // 1st of connected (via myPrev and myNext) edges
 494   std::set< BNode >::iterator bn = bNodes.begin();
 495   for ( ; bn != bNodes.end(); ++bn )
 496   {
 497     for ( size_t i = 0; i < bn->myLinkedEdges.size(); ++i )
 498     {
 499       if ( bn->myLinkedEdges[i]->myBorderID < 0 )
 500       {
 501         BEdge* be = bn->myLinkedEdges[i];
 502         int borderID = borders.size();
 503         borders.push_back( be );
 504         for ( ; be && be->myBorderID < 0; be = be->myNext )
 505         {
 506           be->myBorderID = borderID;
 507           be->Orient();
 508         }
 509         bool isClosed = ( be == bn->myLinkedEdges[i] );
 510         be = bn->myLinkedEdges[i]->myPrev;
 511         for ( ; be && be->myBorderID < 0; be = be->myPrev )
 512         {
 513           be->myBorderID = borderID;
 514           be->Orient();
 515         }
 516         if ( !isClosed )
 517           while ( borders.back()->myPrev )
 518             borders.back() = borders.back()->myPrev;
 519
 520         borders.back()->SetID( 0 ); // set IDs to all edges of the border
 521       }
 522     }
 523   }
 524
 525   // compute tolerance of each border
 526   double maxTolerance = tolerance;
 527   std::vector< double > bordToler( borders.size(), tolerance );
 528   if ( maxTolerance < std::numeric_limits< double >::min() )
 529   {
 530     // no tolerance provided by the user; compute tolerance of each border
 531     // as one tenth of an average size of faces adjacent to a border
 532     for ( size_t i = 0; i < borders.size(); ++i )
 533     {
 534       double avgFaceSize = 0;
 535       int    nbFaces     = 0;
 536       BEdge* be = borders[ i ];
 537       do {
 538         double facePerimeter = 0;
 539         gp_Pnt p0 = SMESH_TNodeXYZ( be->myFace->GetNode( be->myFace->NbNodes() - 1 ));
 540         SMDS_NodeIteratorPtr nodeIt = be->myFace->interlacedNodesIterator();
 541         while ( nodeIt->more() )
 542         {
 543           gp_Pnt p1 = SMESH_TNodeXYZ( nodeIt->next() );
 544           facePerimeter += p0.Distance( p1 );
 545           p0 = p1;
 546         }
 547         avgFaceSize += ( facePerimeter / be->myFace->NbCornerNodes() );
 548         nbFaces++;
 549
 550         be = be->myNext;
 551       }
 552       while ( be && be != borders[i] );
 553
 554       bordToler[ i ] = 0.1 * avgFaceSize / nbFaces;
 555       maxTolerance = Max( maxTolerance, bordToler[ i ]);
 556     }
 557   }
 558
 559   // for every border node find close border edges
 560   SMESH_ElementSearcher* searcher =
 561     GetElementSearcher( mesh, getElemIterator( bEdges ), maxTolerance );
 562   SMESHUtils::Deleter< SMESH_ElementSearcher > searcherDeleter( searcher );
 563   std::vector< const SMDS_MeshElement* > candidateEdges;
 564   for ( bn = bNodes.begin(); bn != bNodes.end(); ++bn )
 565   {
 566     searcher->FindElementsByPoint( *bn, SMDSAbs_Edge, candidateEdges );
 567     if ( candidateEdges.size() <= bn->myLinkedEdges.size() )
 568       continue;
 569
 570     double nodeTol = 0, u;
 571     for ( size_t i = 0; i < bn->myLinkedEdges.size(); ++i )
 572       nodeTol = Max( nodeTol, bordToler[ bn->myLinkedEdges[ i ]->myBorderID ]);
 573
 574     for ( size_t i = 0; i < candidateEdges.size(); ++i )
 575     {
 576       const BEdge* be = static_cast< const BEdge* >( candidateEdges[ i ]);
 577       double      tol = Max( nodeTol, bordToler[ be->myBorderID ]);
 578       if ( !be->IsOut( *bn, tol, u ))
 579         bn->AddClose( be, u );
 580     }
 581   }
 582
 583   // for every border edge find close borders
 584
 585   std::vector< BEdge* > closeEdges;
 586   for ( size_t i = 0; i < bEdges.size(); ++i )
 587   {
 588     BEdge& be = bEdges[i];
 589     if ( be.myBNode1->myCloseEdges.empty() ||
 590          be.myBNode2->myCloseEdges.empty() )
 591       continue;
 592
 593     closeEdges.clear();
 594     for ( size_t iE1 = 0; iE1 < be.myBNode1->myCloseEdges.size(); ++iE1 )
 595     {
 596       // find edges of the same border close to both nodes of the edge
 597       BEdge* closeE1 = be.myBNode1->GetCloseEdge( iE1 );
 598       BEdge* closeE2 = be.myBNode2->GetCloseEdgeOfBorder( closeE1->myBorderID );
 599       if ( !closeE2 )
 600         continue;
 601       // check that edges connecting closeE1 and closeE2 (if any) are also close to 'be'
 602       if ( closeE1 != closeE2 )
 603       {
 604         bool coincide;
 605         for ( int j = 0; j < 2; ++j ) // move closeE1 -> closeE2 or inversely
 606         {
 607           BEdge* ce = closeE1;
 608           do {
 609             coincide = ( ce->myBNode2->GetCloseEdgeOfBorder( be.myBorderID ));
 610             ce       = ce->myNext;
 611           } while ( coincide && ce && ce != closeE2 );
 612
 613           if ( coincide && ce == closeE2 )
 614             break;
 615           if ( j == 0 )
 616             std::swap( closeE1, closeE2 );
 617           coincide = false;
 618         }
 619         if ( !coincide )
 620           continue;
 621         closeEdges.push_back( closeE1 );
 622         closeEdges.push_back( closeE2 );
 623       }
 624       else
 625       {
 626         closeEdges.push_back( closeE1 );
 627       }
 628       be.myCloseBorders.insert( closeE1->myBorderID );
 629     }
 630     if ( !closeEdges.empty() )
 631     {
 632       be.myCloseBorders.insert( be.myBorderID );
 633       // for ( size_t iB = 0; iB < closeEdges.size(); ++iB )
 634       //   closeEdges[ iB ]->myCloseBorders.insert( be.myCloseBorders.begin(),
 635       //                                            be.myCloseBorders.end() );
 636     }
 637   }
 638
 639   // Fill in CoincidentFreeBorders
 640
 641   // save nodes of free borders
 642   foundFreeBordes._borders.resize( borders.size() );
 643   for ( size_t i = 0; i < borders.size(); ++i )
 644   {
 645     BEdge* be = borders[i];
 646     foundFreeBordes._borders[i].push_back( be->myBNode1->Node() );
 647     do {
 648       foundFreeBordes._borders[i].push_back( be->myBNode2->Node() );
 649       be = be->myNext;
 650     }
 651     while ( be && be != borders[i] );
 652   }
 653
 654   // form groups of coincident parts of free borders
 655
 656   TFreeBorderPart       part;
 657   TCoincidentGroup      group;
 658   std::vector< BEdge* > ranges; // couples of edges delimiting parts
 659   BEdge* be = 0; // a current edge
 660   int skipGroup = bEdges.size(); // a group ID used to avoid repeating treatment of edges
 661
 662   for ( int i = 0, nbBords = borders.size(); i < nbBords; i += bool(!be) )
 663   {
 664     if ( !be )
 665       be = borders[i];
 666
 667     // look for an edge close to other borders
 668     do {
 669       if ( !be->IsInGroup() && !be->myCloseBorders.empty() )
 670         break;
 671       be = be->myNext;
 672     } while ( be && be != borders[i] );
 673
 674     if ( !be || be->IsInGroup() || be->myCloseBorders.empty() )
 675     {
 676       be = 0;
 677       continue; // all edges of a border are treated or non-coincident
 678     }
 679     group.clear();
 680     ranges.clear();
 681
 682     // look for the 1st and last edge of a coincident group
 683     BEdge* beRange[2];
 684     if ( !be->GetRangeOfSameCloseBorders( beRange, be->myCloseBorders ))
 685     {
 686       be->myInGroup = skipGroup;
 687       be = be->myNext;
 688       continue;
 689     }
 690
 691     ranges.push_back( beRange[0] );
 692     ranges.push_back( beRange[1] );
 693
 694     int groupID = foundFreeBordes._coincidentGroups.size();
 695     be = beRange[0];
 696     be->myInGroup = groupID;
 697     while ( be != beRange[1] )
 698     {
 699       be->myInGroup = groupID;
 700       be = be->myNext;
 701     }
 702     beRange[1]->myInGroup = groupID;
 703
 704     // get starting edge of each close border
 705     closeEdges.clear();
 706     be = beRange[0];
 707     if ( be->myCloseBorders.empty() )
 708       be = beRange[0]->myNext;
 709     std::set<int>::iterator closeBord = be->myCloseBorders.begin();
 710     for ( ; closeBord != be->myCloseBorders.end(); ++closeBord )
 711       if ( BEdge* e = be->myBNode2->GetCloseEdgeOfBorder( *closeBord ))
 712         closeEdges.push_back( e );
 713
 714     for ( size_t iE = 0; iE < closeEdges.size(); ++iE )
 715       if ( be->myCloseBorders != closeEdges[iE]->myCloseBorders )
 716       {
 717         closeBord = closeEdges[iE]->myCloseBorders.begin();
 718         for ( ; closeBord != closeEdges[iE]->myCloseBorders.end(); ++closeBord )
 719           if ( !be->myCloseBorders.count( *closeBord ))
 720             if ( BEdge* e = closeEdges[iE]->myBNode2->GetCloseEdgeOfBorder( *closeBord ))
 721               if ( std::find( closeEdges.begin(), closeEdges.end(), e ) == closeEdges.end() )
 722                 closeEdges.push_back( e );
 723       }
 724
 725     // add parts of other borders
 726
 727     BEdge* be1st = beRange[0];
 728     for ( size_t iE = 0; iE < closeEdges.size(); ++iE )
 729     {
 730       be = closeEdges[ iE ];
 731       if ( !be ) continue;
 732
 733       bool ok = be->GetRangeOfSameCloseBorders( beRange, be->myCloseBorders );
 734       // if ( !ok && be->myPrev )
 735       //   ok = be->myPrev->GetRangeOfSameCloseBorders( beRange, be1st->myCloseBorders );
 736       // if ( !ok && be->myNext )
 737       //   ok = be->myNext->GetRangeOfSameCloseBorders( beRange, be1st->myCloseBorders );
 738       if ( !ok )
 739         continue;
 740
 741       be = beRange[0];
 742
 743       ranges.push_back( beRange[0] );
 744       ranges.push_back( beRange[1] );
 745
 746       be->myInGroup = groupID;
 747       while ( be != beRange[1] )
 748       {
 749         be->myInGroup = groupID;
 750         be = be->myNext;
 751       }
 752       beRange[1]->myInGroup = groupID;
 753     }
 754
 755     if ( ranges.size() > 2 )
 756     {
 757       if ( !chooseStartOfClosedBorders( ranges ))
 758         for ( size_t iR = 1; iR < ranges.size(); iR += 2 )
 759           extendPart( ranges[ iR-1 ], ranges[ iR ], be1st->myCloseBorders, groupID );
 760
 761       // fill in a group
 762       beRange[0] = ranges[0];
 763       beRange[1] = ranges[1];
 764
 765       part._border   = i;
 766       part._node1    = beRange[0]->myID;
 767       part._node2    = beRange[0]->myID + 1;
 768       part._nodeLast = beRange[1]->myID + 1;
 769       group.push_back( part );
 770
 771       be1st = beRange[0];
 772       for ( size_t iR = 3; iR < ranges.size(); iR += 2 )
 773       {
 774         beRange[0] = ranges[iR-1];
 775         beRange[1] = ranges[iR-0];
 776
 777         // find out mutual orientation of borders
 778         double u1, u2;
 779         be1st       ->IsOut( *beRange[ 0 ]->myBNode1, maxTolerance, u1 );
 780         beRange[ 0 ]->IsOut( *be1st->myBNode1,        maxTolerance, u2 );
 781         bool reverse = (( u1 < 0 || u1 > 1 ) && ( u2 < 0 || u2 > 1 ));
 782
 783         // fill in a group
 784         part._border   = beRange[0]->myBorderID;
 785         if ( reverse ) {
 786           part._node1    = beRange[1]->myID + 1;
 787           part._node2    = beRange[1]->myID;
 788           part._nodeLast = beRange[0]->myID;
 789         }
 790         else  {
 791           part._node1    = beRange[0]->myID;
 792           part._node2    = beRange[0]->myID + 1;
 793           part._nodeLast = beRange[1]->myID + 1;
 794         }
 795         // if ( group[0]._node2 != part._node2 )
 796         group.push_back( part );
 797       }
 798       //if ( group.size() > 1 )
 799         foundFreeBordes._coincidentGroups.push_back( group );
 800     }
 801     else
 802     {
 803       beRange[0] = ranges[0];
 804       beRange[1] = ranges[1];
 805
 806       be = beRange[0];
 807       be->myInGroup = skipGroup;
 808       while ( be != beRange[1] )
 809       {
 810         be->myInGroup = skipGroup;
 811         be = be->myNext;
 812       }
 813       beRange[1]->myInGroup = skipGroup;
 814     }
 815
 816     be = ranges[1];
 817
 818   } // loop on free borders
 819
 820   return;
 821
 822 } // SMESH_MeshAlgos::FindCoincidentFreeBorders()
 823