src/SMESHUtils/SMESH_FreeBorders.cxx

   1 // Copyright (C) 2007-2021  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
 135         for ( BEdge* be = myNext; be && be->myID < 0; be = be->myNext )
 136         {
 137           be->myID = ++id;
 138         }
 139       }
 140     }
 141     //================================================================================
 142     /*!
 143      * \brief Checks if a point is closer to this BEdge than tol
 144      */
 145     //================================================================================
 146
 147     bool IsOut( const gp_XYZ& point, const double tol, double& u ) const
 148     {
 149       gp_XYZ  me = *myBNode2 - *myBNode1;
 150       gp_XYZ n1p = point     - *myBNode1;
 151       u = ( me * n1p ) / me.SquareModulus(); // param [0,1] on this
 152       if ( u < 0. ) return ( n1p.SquareModulus() > tol * tol );
 153       if ( u > 1. ) return ( ( point - *myBNode2 ).SquareModulus() > tol * tol );
 154
 155       gp_XYZ  proj = ( 1. - u ) * *myBNode1 + u * *myBNode2; // projection of the point on this
 156       double dist2 = ( point - proj ).SquareModulus();
 157       return ( dist2 > tol * tol );
 158     }
 159     //================================================================================
 160     /*!
 161      * \brief Checks if two BEdges can be considered as overlapping
 162      */
 163     //================================================================================
 164
 165     bool IsOverlappingProjection( const BEdge* toE, const double u, bool is1st ) const
 166     {
 167       // is1st shows which end of toE is projected on this at u
 168       double u2;
 169       const double eps = 0.1;
 170       if ( myBNode1->IsCloseEdge( toE, &u2 ) ||
 171            myBNode2->IsCloseEdge( toE, &u2 ))
 172         return (( 0 < u2 && u2 < 1 ) &&     // u2 is proj param of myBNode's on toE
 173                 ( Abs( u2 - int( !is1st )) > eps ));
 174
 175       const BNode* n = is1st ? toE->myBNode2 : toE->myBNode1;
 176       if ( this == n->GetCloseEdgeOfBorder( this->myBorderID, &u2 ))
 177         return Abs( u - u2 ) > eps;
 178       return false;
 179     }
 180     //================================================================================
 181     /*!
 182      * \brief Finds all neighbor BEdge's having the same close borders
 183      */
 184     //================================================================================
 185
 186     bool GetRangeOfSameCloseBorders(BEdge* eRange[2], const std::set< int >& bordIDs)
 187     {
 188       if ( this->myCloseBorders != bordIDs )
 189         return false;
 190
 191       if ( bordIDs.size() == 1 && bordIDs.count( myBorderID )) // border close to self
 192       {
 193         double u;
 194         eRange[0] = this;
 195         while ( eRange[0]->myBNode1->GetCloseEdgeOfBorder( myBorderID, &u ))
 196         {
 197           if ( eRange[0]->myPrev == this || u < 0 || u > 1 )
 198             break;
 199           eRange[0] = eRange[0]->myPrev;
 200         }
 201         eRange[1] = this;
 202         while ( eRange[1]->myBNode2->GetCloseEdgeOfBorder( myBorderID, &u ))
 203         {
 204           if ( eRange[1]->myNext == this || u < 0 || u > 1 )
 205             break;
 206           eRange[1] = eRange[1]->myNext;
 207         }
 208       }
 209       else
 210       {
 211         eRange[0] = this;
 212         while ( eRange[0]->myPrev && eRange[0]->myPrev->myCloseBorders == bordIDs )
 213         {
 214           if ( eRange[0]->myPrev == this )
 215             break;
 216           eRange[0] = eRange[0]->myPrev;
 217         }
 218
 219         eRange[1] = this;
 220         if ( eRange[0]->myPrev != this ) // not closed border
 221           while ( eRange[1]->myNext && eRange[1]->myNext->myCloseBorders == bordIDs )
 222           {
 223             if ( eRange[1]->myNext == this )
 224               break;
 225             eRange[1] = eRange[1]->myNext;
 226           }
 227       }
 228
 229       if ( eRange[0] == eRange[1] )
 230       {
 231         std::set<int>::iterator closeBord = eRange[0]->myCloseBorders.begin();
 232         for ( ; closeBord != eRange[0]->myCloseBorders.end(); ++closeBord )
 233         {
 234           if ( BEdge* be = eRange[0]->myBNode1->GetCloseEdgeOfBorder( *closeBord ))
 235             if ( be->myCloseBorders == eRange[0]->myCloseBorders )
 236               return true;
 237           if ( BEdge* be = eRange[0]->myBNode2->GetCloseEdgeOfBorder( *closeBord ))
 238             if ( be->myCloseBorders == eRange[0]->myCloseBorders )
 239               return true;
 240         }
 241         return false;
 242       }
 243       return true;
 244     }
 245   }; // class BEdge
 246
 247   //================================================================================
 248   /*!
 249    * \brief Checks if all border parts include the whole closed border, and if so
 250    *        returns \c true and choose starting BEdge's with most coincident nodes
 251    */
 252   //================================================================================
 253
 254   bool chooseStartOfClosedBorders( std::vector< BEdge* >& ranges ) // PAL23078#c21002
 255   {
 256     bool allClosed = true;
 257     for ( size_t iR = 1; iR < ranges.size() && allClosed; iR += 2 )
 258       allClosed = ( ranges[ iR-1 ]->myPrev == ranges[ iR ] );
 259     if ( !allClosed )
 260       return allClosed;
 261
 262     double u, minDiff = Precision::Infinite();
 263     std::vector< BEdge* > start( ranges.size() / 2 );
 264     BEdge* range0 = start[0] = ranges[0];
 265     do
 266     {
 267       double maxDiffU = 0;
 268       double  maxDiff = 0;
 269       for ( size_t iR = 3; iR < ranges.size(); iR += 2 )
 270       {
 271         int borderID = ranges[iR]->myBorderID;
 272         if ( BEdge* e = start[0]->myBNode1->GetCloseEdgeOfBorder( borderID, & u ))
 273         {
 274           start[ iR / 2 ] = e;
 275           double diffU = Min( Abs( u ), Abs( 1.-u ));
 276           double  diff = e->myBNode1->SquareDistance( *e->myBNode2 ) * diffU * diffU;
 277           maxDiffU = Max( diffU, maxDiffU );
 278           maxDiff  = Max( diff,  maxDiff );
 279         }
 280       }
 281       if ( maxDiff < minDiff )
 282       {
 283         minDiff = maxDiff;
 284         for ( size_t iR = 1; iR < ranges.size(); iR += 2 )
 285         {
 286           ranges[ iR-1 ] = start[ iR/2 ];
 287           ranges[ iR   ] = ranges[ iR-1]->myPrev;
 288         }
 289       }
 290       if ( maxDiffU < 1e-6 )
 291         break;
 292       start[0] = start[0]->myNext;
 293     }
 294     while ( start[0] != range0 );
 295
 296     return allClosed;
 297   }
 298
 299   //================================================================================
 300   /*!
 301    * \brief Tries to include neighbor BEdge's into a border part
 302    */
 303   //================================================================================
 304
 305   void extendPart( BEdge* & e1, BEdge* & e2, const std::set< int >& bordIDs, int groupID )
 306   {
 307     if (( e1->myPrev == e2 ) ||
 308         ( e1 == e2 && e1->myPrev && e1->myPrev->myInGroup == groupID ))
 309       return; // full free border already
 310
 311     double u;
 312     BEdge* be;
 313     std::set<int>::const_iterator bord;
 314     if ( e1->myPrev )
 315     {
 316       for ( bord = bordIDs.begin(); bord != bordIDs.end(); ++bord )
 317         if ((( be = e1->myBNode1->GetCloseEdgeOfBorder( *bord, &u ))) &&
 318             (  be->myInGroup == groupID ) &&
 319             (  0 < u && u < 1 ) &&
 320             (  be->IsOverlappingProjection( e1->myPrev, u, false )))
 321         {
 322           e1 = e1->myPrev;
 323           break;
 324         }
 325       if ( bord == bordIDs.end() && // not extended
 326            e1->myBNode1->HasCloseEdgeWithNode( e1->myPrev->myBNode1 ))
 327       {
 328         e1 = e1->myPrev;
 329       }
 330       e1->myInGroup = groupID;
 331     }
 332     if ( e2->myNext )
 333     {
 334       for ( bord = bordIDs.begin(); bord != bordIDs.end(); ++bord )
 335         if ((( be = e2->myBNode2->GetCloseEdgeOfBorder( *bord, &u ))) &&
 336             (  be->myInGroup == groupID ) &&
 337             (  0 < u && u < 1 ) &&
 338             (  be->IsOverlappingProjection( e2->myNext, u, true )))
 339         {
 340           e2 = e2->myNext;
 341           break;
 342         }
 343       if ( bord == bordIDs.end() && // not extended
 344            e2->myBNode2->HasCloseEdgeWithNode( e2->myNext->myBNode2 ))
 345       {
 346         e2 = e2->myNext;
 347       }
 348       e2->myInGroup = groupID;
 349     }
 350   }
 351
 352   //================================================================================
 353   /*!
 354    * \brief Connect BEdge's incident at this node
 355    */
 356   //================================================================================
 357
 358   void BNode::AddLinked( BEdge* e ) const
 359   {
 360     myLinkedEdges.reserve(2);
 361     myLinkedEdges.push_back( e );
 362     if ( myLinkedEdges.size() < 2 ) return;
 363
 364     if ( myLinkedEdges.size() == 2 )
 365     {
 366       myLinkedEdges[0]->AddLinked( myLinkedEdges[1] );
 367       myLinkedEdges[1]->AddLinked( myLinkedEdges[0] );
 368     }
 369     else
 370     {
 371       for ( size_t i = 0; i < myLinkedEdges.size(); ++i )
 372         for ( size_t j = 0; j < myLinkedEdges.size(); ++j )
 373           if ( i != j )
 374             myLinkedEdges[i]->RemoveLinked( myLinkedEdges[j] );
 375     }
 376   }
 377   void BNode::AddClose ( const BEdge* e, double u ) const
 378   {
 379     if ( ! e->Contains( this ))
 380       myCloseEdges.push_back( std::make_pair( const_cast< BEdge* >( e ), u ));
 381   }
 382   BEdge* BNode::GetCloseEdgeOfBorder( int borderID, double * uPtr ) const
 383   {
 384     BEdge* e = 0;
 385     double u = 0;
 386     for ( size_t i = 0; i < myCloseEdges.size(); ++i )
 387       if ( borderID == GetCloseEdge( i )->myBorderID )
 388       {
 389         if ( e && Abs( u - 0.5 ) < Abs( GetCloseU( i ) - 0.5 ))
 390           continue;
 391         u = GetCloseU( i );
 392         e = GetCloseEdge ( i );
 393       }
 394     if ( uPtr ) *uPtr = u;
 395     return e;
 396   }
 397   bool BNode::HasCloseEdgeWithNode( const BNode* n ) const
 398   {
 399     for ( size_t i = 0; i < myCloseEdges.size(); ++i )
 400       if ( GetCloseEdge( i )->Contains( n ) &&
 401            0 < GetCloseU( i ) && GetCloseU( i ) < 1 )
 402         return true;
 403     return false;
 404   }
 405   bool BNode::IsCloseEdge( const BEdge* e, double * uPtr ) const
 406   {
 407     for ( size_t i = 0; i < myCloseEdges.size(); ++i )
 408       if ( e == GetCloseEdge( i ) )
 409       {
 410         if ( uPtr ) *uPtr = GetCloseU( i );
 411         return true;
 412       }
 413     return false;
 414   }
 415
 416   /// Accessor to SMDS_MeshElement* inherited by BEdge
 417   struct ElemAcess
 418   {
 419     static const SMDS_MeshElement* value( std::vector< BEdge >::const_iterator it)
 420     {
 421       return & (*it);
 422     }
 423   };
 424   /// Iterator over a vector of BEdge's
 425   static SMDS_ElemIteratorPtr getElemIterator( const std::vector< BEdge > & bedges )
 426   {
 427     typedef SMDS_SetIterator
 428       < const SMDS_MeshElement*, std::vector< BEdge >::const_iterator, ElemAcess > BEIter;
 429     return SMDS_ElemIteratorPtr( new BEIter( bedges.begin(), bedges.end() ));
 430   }
 431
 432 } // namespace
 433
 434 //================================================================================
 435 /*
 436  * Returns groups of TFreeBorder's coincident within the given tolerance.
 437  * If the tolerance <= 0.0 then one tenth of an average size of elements adjacent
 438  * to free borders being compared is used.
 439  */
 440 //================================================================================
 441
 442 void SMESH_MeshAlgos::FindCoincidentFreeBorders(SMDS_Mesh&              mesh,
 443                                                 double                  tolerance,
 444                                                 CoincidentFreeBorders & foundFreeBordes)
 445 {
 446   // find free links
 447   typedef NCollection_DataMap<SMESH_TLink, const SMDS_MeshElement*, SMESH_TLink > TLink2FaceMap;
 448   TLink2FaceMap linkMap;
 449   int nbSharedLinks = 0;
 450   SMDS_FaceIteratorPtr faceIt = mesh.facesIterator();
 451   while ( faceIt->more() )
 452   {
 453     const SMDS_MeshElement* face = faceIt->next();
 454     if ( !face ) continue;
 455
 456     const SMDS_MeshNode*     n0 = face->GetNode( face->NbNodes() - 1 );
 457     SMDS_NodeIteratorPtr nodeIt = face->interlacedNodesIterator();
 458     while ( nodeIt->more() )
 459     {
 460       const SMDS_MeshNode* n1 = nodeIt->next();
 461       SMESH_TLink link( n0, n1 );
 462       if ( const SMDS_MeshElement** faceInMap = linkMap.ChangeSeek( link ))
 463       {
 464         nbSharedLinks += bool( *faceInMap );
 465         *faceInMap = 0;
 466       }
 467       else
 468       {
 469         linkMap.Bind( link, face );
 470       }
 471       n0 = n1;
 472     }
 473   }
 474   if ( linkMap.Extent() == nbSharedLinks )
 475     return;
 476
 477   // form free borders
 478   std::set   < BNode > bNodes;
 479   std::vector< BEdge > bEdges( linkMap.Extent() - nbSharedLinks );
 480
 481   TLink2FaceMap::Iterator linkIt( linkMap );
 482   for ( int iEdge = 0; linkIt.More(); linkIt.Next() )
 483   {
 484     if ( !linkIt.Value() ) continue;
 485     const SMESH_TLink & link = linkIt.Key();
 486     std::set< BNode >::iterator n1 = bNodes.insert( BNode( link.node1() )).first;
 487     std::set< BNode >::iterator n2 = bNodes.insert( BNode( link.node2() )).first;
 488     bEdges[ iEdge ].Set( &*n1, &*n2, linkIt.Value(), iEdge+1 );
 489     n1->AddLinked( & bEdges[ iEdge ] );
 490     n2->AddLinked( & bEdges[ iEdge ] );
 491     ++iEdge;
 492   }
 493   linkMap.Clear();
 494
 495   // assign IDs to borders
 496   std::vector< BEdge* > borders; // 1st of connected (via myPrev and myNext) edges
 497   std::set< BNode >::iterator bn = bNodes.begin();
 498   for ( ; bn != bNodes.end(); ++bn )
 499   {
 500     for ( size_t i = 0; i < bn->myLinkedEdges.size(); ++i )
 501     {
 502       if ( bn->myLinkedEdges[i]->myBorderID < 0 )
 503       {
 504         BEdge* be = bn->myLinkedEdges[i];
 505         int borderID = borders.size();
 506         borders.push_back( be );
 507         for ( ; be && be->myBorderID < 0; be = be->myNext )
 508         {
 509           be->myBorderID = borderID;
 510           be->Orient();
 511         }
 512         bool isClosed = ( be == bn->myLinkedEdges[i] );
 513         be = bn->myLinkedEdges[i]->myPrev;
 514         for ( ; be && be->myBorderID < 0; be = be->myPrev )
 515         {
 516           be->myBorderID = borderID;
 517           be->Orient();
 518         }
 519         if ( !isClosed )
 520           while ( borders.back()->myPrev )
 521             borders.back() = borders.back()->myPrev;
 522
 523         borders.back()->SetID( 0 ); // set IDs to all edges of the border
 524       }
 525     }
 526   }
 527
 528   // compute tolerance of each border
 529   double maxTolerance = tolerance;
 530   std::vector< double > bordToler( borders.size(), tolerance );
 531   if ( maxTolerance < std::numeric_limits< double >::min() )
 532   {
 533     // no tolerance provided by the user; compute tolerance of each border
 534     // as one tenth of an average size of faces adjacent to a border
 535     for ( size_t i = 0; i < borders.size(); ++i )
 536     {
 537       double avgFaceSize = 0;
 538       int    nbFaces     = 0;
 539       BEdge* be = borders[ i ];
 540       do {
 541         double facePerimeter = 0;
 542         gp_Pnt p0 = SMESH_TNodeXYZ( be->myFace->GetNode( be->myFace->NbNodes() - 1 ));
 543         SMDS_NodeIteratorPtr nodeIt = be->myFace->interlacedNodesIterator();
 544         while ( nodeIt->more() )
 545         {
 546           gp_Pnt p1 = SMESH_TNodeXYZ( nodeIt->next() );
 547           facePerimeter += p0.Distance( p1 );
 548           p0 = p1;
 549         }
 550         avgFaceSize += ( facePerimeter / be->myFace->NbCornerNodes() );
 551         nbFaces++;
 552
 553         be = be->myNext;
 554       }
 555       while ( be && be != borders[i] );
 556
 557       bordToler[ i ] = 0.1 * avgFaceSize / nbFaces;
 558       maxTolerance = Max( maxTolerance, bordToler[ i ]);
 559     }
 560   }
 561
 562   // for every border node find close border edges
 563   SMESH_ElementSearcher* searcher =
 564     GetElementSearcher( mesh, getElemIterator( bEdges ), maxTolerance );
 565   SMESHUtils::Deleter< SMESH_ElementSearcher > searcherDeleter( searcher );
 566   std::vector< const SMDS_MeshElement* > candidateEdges;
 567   for ( bn = bNodes.begin(); bn != bNodes.end(); ++bn )
 568   {
 569     searcher->FindElementsByPoint( *bn, SMDSAbs_Edge, candidateEdges );
 570     if ( candidateEdges.size() <= bn->myLinkedEdges.size() )
 571       continue;
 572
 573     double nodeTol = 0, u;
 574     for ( size_t i = 0; i < bn->myLinkedEdges.size(); ++i )
 575       nodeTol = Max( nodeTol, bordToler[ bn->myLinkedEdges[ i ]->myBorderID ]);
 576
 577     for ( size_t i = 0; i < candidateEdges.size(); ++i )
 578     {
 579       const BEdge* be = static_cast< const BEdge* >( candidateEdges[ i ]);
 580       double      tol = Max( nodeTol, bordToler[ be->myBorderID ]);
 581       if ( !be->IsOut( *bn, tol, u ))
 582         bn->AddClose( be, u );
 583     }
 584   }
 585
 586   // for every border edge find close borders
 587
 588   std::vector< BEdge* > closeEdges;
 589   for ( size_t i = 0; i < bEdges.size(); ++i )
 590   {
 591     BEdge& be = bEdges[i];
 592     if ( be.myBNode1->myCloseEdges.empty() ||
 593          be.myBNode2->myCloseEdges.empty() )
 594       continue;
 595
 596     closeEdges.clear();
 597     for ( size_t iE1 = 0; iE1 < be.myBNode1->myCloseEdges.size(); ++iE1 )
 598     {
 599       // find edges of the same border close to both nodes of the edge
 600       BEdge* closeE1 = be.myBNode1->GetCloseEdge( iE1 );
 601       BEdge* closeE2 = be.myBNode2->GetCloseEdgeOfBorder( closeE1->myBorderID );
 602       if ( !closeE2 )
 603         continue;
 604       // check that edges connecting closeE1 and closeE2 (if any) are also close to 'be'
 605       if ( closeE1 != closeE2 )
 606       {
 607         bool coincide;
 608         for ( int j = 0; j < 2; ++j ) // move closeE1 -> closeE2 or inversely
 609         {
 610           BEdge* ce = closeE1;
 611           do {
 612             coincide = ( ce->myBNode2->GetCloseEdgeOfBorder( be.myBorderID ));
 613             ce       = ce->myNext;
 614           } while ( coincide && ce && ce != closeE2 );
 615
 616           if ( coincide && ce == closeE2 )
 617             break;
 618           if ( j == 0 )
 619             std::swap( closeE1, closeE2 );
 620           coincide = false;
 621         }
 622         if ( !coincide )
 623           continue;
 624         closeEdges.push_back( closeE1 );
 625         closeEdges.push_back( closeE2 );
 626       }
 627       else
 628       {
 629         closeEdges.push_back( closeE1 );
 630       }
 631       be.myCloseBorders.insert( closeE1->myBorderID );
 632     }
 633     if ( !closeEdges.empty() )
 634     {
 635       be.myCloseBorders.insert( be.myBorderID );
 636       // for ( size_t iB = 0; iB < closeEdges.size(); ++iB )
 637       //   closeEdges[ iB ]->myCloseBorders.insert( be.myCloseBorders.begin(),
 638       //                                            be.myCloseBorders.end() );
 639     }
 640   }
 641
 642   // Fill in CoincidentFreeBorders
 643
 644   // save nodes of free borders
 645   foundFreeBordes._borders.resize( borders.size() );
 646   for ( size_t i = 0; i < borders.size(); ++i )
 647   {
 648     BEdge* be = borders[i];
 649     foundFreeBordes._borders[i].push_back( be->myBNode1->Node() );
 650     do {
 651       foundFreeBordes._borders[i].push_back( be->myBNode2->Node() );
 652       be = be->myNext;
 653     }
 654     while ( be && be != borders[i] );
 655   }
 656
 657   // form groups of coincident parts of free borders
 658
 659   TFreeBorderPart       part;
 660   TCoincidentGroup      group;
 661   std::vector< BEdge* > ranges; // couples of edges delimiting parts
 662   BEdge* be = 0; // a current edge
 663   int skipGroup = bEdges.size(); // a group ID used to avoid repeating treatment of edges
 664
 665   for ( int i = 0, nbBords = borders.size(); i < nbBords; i += bool(!be) )
 666   {
 667     if ( !be )
 668       be = borders[i];
 669
 670     // look for an edge close to other borders
 671     do {
 672       if ( !be->IsInGroup() && !be->myCloseBorders.empty() )
 673         break;
 674       be = be->myNext;
 675     } while ( be && be != borders[i] );
 676
 677     if ( !be || be->IsInGroup() || be->myCloseBorders.empty() )
 678     {
 679       be = 0;
 680       continue; // all edges of a border are treated or non-coincident
 681     }
 682     group.clear();
 683     ranges.clear();
 684
 685     // look for the 1st and last edge of a coincident group
 686     BEdge* beRange[2];
 687     if ( !be->GetRangeOfSameCloseBorders( beRange, be->myCloseBorders ))
 688     {
 689       be->myInGroup = skipGroup;
 690       be = be->myNext;
 691       continue;
 692     }
 693
 694     ranges.push_back( beRange[0] );
 695     ranges.push_back( beRange[1] );
 696
 697     int groupID = foundFreeBordes._coincidentGroups.size();
 698     be = beRange[0];
 699     be->myInGroup = groupID;
 700     while ( be != beRange[1] )
 701     {
 702       be->myInGroup = groupID;
 703       be = be->myNext;
 704     }
 705     beRange[1]->myInGroup = groupID;
 706
 707     // get starting edge of each close border
 708     closeEdges.clear();
 709     be = beRange[0];
 710     if ( be->myCloseBorders.empty() )
 711       be = beRange[0]->myNext;
 712     std::set<int>::iterator closeBord = be->myCloseBorders.begin();
 713     for ( ; closeBord != be->myCloseBorders.end(); ++closeBord )
 714       if ( BEdge* e = be->myBNode2->GetCloseEdgeOfBorder( *closeBord ))
 715         closeEdges.push_back( e );
 716
 717     for ( size_t iE = 0; iE < closeEdges.size(); ++iE )
 718       if ( be->myCloseBorders != closeEdges[iE]->myCloseBorders )
 719       {
 720         closeBord = closeEdges[iE]->myCloseBorders.begin();
 721         for ( ; closeBord != closeEdges[iE]->myCloseBorders.end(); ++closeBord )
 722           if ( !be->myCloseBorders.count( *closeBord ))
 723             if ( BEdge* e = closeEdges[iE]->myBNode2->GetCloseEdgeOfBorder( *closeBord ))
 724               if ( std::find( closeEdges.begin(), closeEdges.end(), e ) == closeEdges.end() )
 725                 closeEdges.push_back( e );
 726       }
 727
 728     // add parts of other borders
 729
 730     BEdge* be1st = beRange[0];
 731     for ( size_t iE = 0; iE < closeEdges.size(); ++iE )
 732     {
 733       be = closeEdges[ iE ];
 734       if ( !be ) continue;
 735
 736       bool ok = be->GetRangeOfSameCloseBorders( beRange, be->myCloseBorders );
 737       // if ( !ok && be->myPrev )
 738       //   ok = be->myPrev->GetRangeOfSameCloseBorders( beRange, be1st->myCloseBorders );
 739       // if ( !ok && be->myNext )
 740       //   ok = be->myNext->GetRangeOfSameCloseBorders( beRange, be1st->myCloseBorders );
 741       if ( !ok )
 742         continue;
 743
 744       be = beRange[0];
 745
 746       ranges.push_back( beRange[0] );
 747       ranges.push_back( beRange[1] );
 748
 749       be->myInGroup = groupID;
 750       while ( be != beRange[1] )
 751       {
 752         be->myInGroup = groupID;
 753         be = be->myNext;
 754       }
 755       beRange[1]->myInGroup = groupID;
 756     }
 757
 758     if ( ranges.size() > 2 )
 759     {
 760       if ( !chooseStartOfClosedBorders( ranges ))
 761         for ( size_t iR = 1; iR < ranges.size(); iR += 2 )
 762           extendPart( ranges[ iR-1 ], ranges[ iR ], be1st->myCloseBorders, groupID );
 763
 764       // fill in a group
 765       beRange[0] = ranges[0];
 766       beRange[1] = ranges[1];
 767
 768       part._border   = i;
 769       part._node1    = beRange[0]->myID;
 770       part._node2    = beRange[0]->myID + 1;
 771       part._nodeLast = beRange[1]->myID + 1;
 772       group.push_back( part );
 773
 774       be1st = beRange[0];
 775       for ( size_t iR = 3; iR < ranges.size(); iR += 2 )
 776       {
 777         beRange[0] = ranges[iR-1];
 778         beRange[1] = ranges[iR-0];
 779
 780         // find out mutual orientation of borders
 781         double u1, u2;
 782         be1st       ->IsOut( *beRange[ 0 ]->myBNode1, maxTolerance, u1 );
 783         beRange[ 0 ]->IsOut( *be1st->myBNode1,        maxTolerance, u2 );
 784         bool reverse = (( u1 < 0 || u1 > 1 ) && ( u2 < 0 || u2 > 1 ));
 785
 786         // fill in a group
 787         part._border   = beRange[0]->myBorderID;
 788         if ( reverse ) {
 789           part._node1    = beRange[1]->myID + 1;
 790           part._node2    = beRange[1]->myID;
 791           part._nodeLast = beRange[0]->myID;
 792         }
 793         else  {
 794           part._node1    = beRange[0]->myID;
 795           part._node2    = beRange[0]->myID + 1;
 796           part._nodeLast = beRange[1]->myID + 1;
 797         }
 798         // if ( group[0]._node2 != part._node2 )
 799         group.push_back( part );
 800       }
 801       //if ( group.size() > 1 )
 802         foundFreeBordes._coincidentGroups.push_back( group );
 803     }
 804     else
 805     {
 806       beRange[0] = ranges[0];
 807       beRange[1] = ranges[1];
 808
 809       be = beRange[0];
 810       be->myInGroup = skipGroup;
 811       while ( be != beRange[1] )
 812       {
 813         be->myInGroup = skipGroup;
 814         be = be->myNext;
 815       }
 816       beRange[1]->myInGroup = skipGroup;
 817     }
 818
 819     be = ranges[1];
 820
 821   } // loop on free borders
 822
 823   return;
 824
 825 } // SMESH_MeshAlgos::FindCoincidentFreeBorders()
 826
 827 //================================================================================
 828 /*
 829  * Returns all TFreeBorder's. Optionally check if the mesh is manifold
 830  * and if faces are correctly oriented.
 831  */
 832 //================================================================================
 833
 834 void SMESH_MeshAlgos::FindFreeBorders(SMDS_Mesh&       theMesh,
 835                                       TFreeBorderVec & theFoundFreeBordes,
 836                                       const bool       theClosedOnly,
 837                                       bool*            theIsManifold,
 838                                       bool*            theIsGoodOri)
 839 {
 840   bool isManifold = true;
 841
 842   // find free links
 843   typedef NCollection_DataMap<SMESH_TLink, const SMDS_MeshElement*, SMESH_TLink > TLink2FaceMap;
 844   TLink2FaceMap linkMap;
 845   int nbSharedLinks = 0;
 846   SMDS_FaceIteratorPtr faceIt = theMesh.facesIterator();
 847   while ( faceIt->more() )
 848   {
 849     const SMDS_MeshElement* face = faceIt->next();
 850     if ( !face ) continue;
 851
 852     const SMDS_MeshNode*     n0 = face->GetNode( face->NbNodes() - 1 );
 853     SMDS_NodeIteratorPtr nodeIt = face->interlacedNodesIterator();
 854     while ( nodeIt->more() )
 855     {
 856       const SMDS_MeshNode* n1 = nodeIt->next();
 857       SMESH_TLink link( n0, n1 );
 858       if ( const SMDS_MeshElement** faceInMap = linkMap.ChangeSeek( link ))
 859       {
 860         if ( *faceInMap )
 861         {
 862           if ( theIsGoodOri && *theIsGoodOri && !IsRightOrder( *faceInMap, n1, n0 ))
 863             *theIsGoodOri = false;
 864         }
 865         else
 866         {
 867           isManifold = false;
 868         }
 869         nbSharedLinks += bool( *faceInMap );
 870         *faceInMap = 0;
 871       }
 872       else
 873       {
 874         linkMap.Bind( link, face );
 875       }
 876       n0 = n1;
 877     }
 878   }
 879   if ( theIsManifold )
 880     *theIsManifold = isManifold;
 881
 882   if ( linkMap.Extent() == nbSharedLinks )
 883     return;
 884
 885   // form free borders
 886   std::set   < BNode > bNodes;
 887   std::vector< BEdge > bEdges( linkMap.Extent() - nbSharedLinks );
 888
 889   TLink2FaceMap::Iterator linkIt( linkMap );
 890   for ( int iEdge = 0; linkIt.More(); linkIt.Next() )
 891   {
 892     if ( !linkIt.Value() ) continue;
 893     const SMESH_TLink & link = linkIt.Key();
 894     std::set< BNode >::iterator n1 = bNodes.insert( BNode( link.node1() )).first;
 895     std::set< BNode >::iterator n2 = bNodes.insert( BNode( link.node2() )).first;
 896     bEdges[ iEdge ].Set( &*n1, &*n2, linkIt.Value(), iEdge+1 );
 897     n1->AddLinked( & bEdges[ iEdge ] );
 898     n2->AddLinked( & bEdges[ iEdge ] );
 899     ++iEdge;
 900   }
 901   linkMap.Clear();
 902
 903   // assign IDs to borders
 904   std::vector< BEdge* > borders; // 1st of connected (via myPrev and myNext) edges
 905   std::set< BNode >::iterator bn = bNodes.begin();
 906   for ( ; bn != bNodes.end(); ++bn )
 907   {
 908     for ( size_t i = 0; i < bn->myLinkedEdges.size(); ++i )
 909     {
 910       if ( bn->myLinkedEdges[i]->myBorderID < 0 )
 911       {
 912         BEdge* be = bn->myLinkedEdges[i];
 913         int borderID = borders.size();
 914         borders.push_back( be );
 915         for ( ; be && be->myBorderID < 0; be = be->myNext )
 916         {
 917           be->myBorderID = borderID;
 918           be->Orient();
 919         }
 920         bool isClosed = ( be == bn->myLinkedEdges[i] );
 921         if ( !isClosed && theClosedOnly )
 922         {
 923           borders.pop_back();
 924           continue;
 925         }
 926         be = bn->myLinkedEdges[i]->myPrev;
 927         for ( ; be && be->myBorderID < 0; be = be->myPrev )
 928         {
 929           be->myBorderID = borderID;
 930           be->Orient();
 931         }
 932         if ( !isClosed )
 933           while ( borders.back()->myPrev )
 934             borders.back() = borders.back()->myPrev;
 935       }
 936     }
 937   }
 938   theFoundFreeBordes.resize( borders.size() );
 939   for ( size_t i = 0; i < borders.size(); ++i )
 940   {
 941     TFreeBorder & bordNodes = theFoundFreeBordes[ i ];
 942     BEdge*               be = borders[i];
 943
 944     size_t cnt = 1;
 945     for ( be = be->myNext; be && be != borders[i]; be = be->myNext )
 946       ++cnt;
 947     bordNodes.resize( cnt + 1 );
 948
 949     BEdge* beLast = 0;
 950     for ( be = borders[i], cnt = 0;
 951           be && cnt < bordNodes.size()-1;
 952           be = be->myNext, ++cnt )
 953     {
 954       bordNodes[ cnt ] = be->myBNode1->Node();
 955       beLast = be;
 956     }
 957     if ( beLast )
 958       bordNodes.back() = beLast->myBNode2->Node();
 959   }
 960 }