src/SMESHUtils/SMESH_MeshAlgos.cxx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 // File      : SMESH_MeshAlgos.hxx
  23 // Created   : Tue Apr 30 18:00:36 2013
  24 // Author    : Edward AGAPOV (eap)
  25
  26 // This file holds some low level algorithms extracted from SMESH_MeshEditor
  27 // to make them accessible from Controls package
  28
  29 #include "SMESH_MeshAlgos.hxx"
  30
  31 #include "ObjectPool.hxx"
  32 #include "SMDS_FaceOfNodes.hxx"
  33 #include "SMDS_LinearEdge.hxx"
  34 #include "SMDS_Mesh.hxx"
  35 #include "SMDS_PolygonalFaceOfNodes.hxx"
  36 #include "SMDS_VolumeTool.hxx"
  37 #include "SMESH_OctreeNode.hxx"
  38
  39 #include <Utils_SALOME_Exception.hxx>
  40
  41 #include <GC_MakeSegment.hxx>
  42 #include <GeomAPI_ExtremaCurveCurve.hxx>
  43 #include <Geom_Line.hxx>
  44 #include <IntAna_IntConicQuad.hxx>
  45 #include <IntAna_Quadric.hxx>
  46 #include <gp_Lin.hxx>
  47 #include <gp_Pln.hxx>
  48 #include <NCollection_DataMap.hxx>
  49
  50 #include <limits>
  51 #include <numeric>
  52
  53 #include <boost/container/flat_set.hpp>
  54
  55 //=======================================================================
  56 /*!
  57  * \brief Implementation of search for the node closest to point
  58  */
  59 //=======================================================================
  60
  61 struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
  62 {
  63   //---------------------------------------------------------------------
  64   /*!
  65    * \brief Constructor
  66    */
  67   SMESH_NodeSearcherImpl( const SMDS_Mesh*     theMesh   = 0,
  68                           SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr() )
  69   {
  70     myMesh = ( SMDS_Mesh* ) theMesh;
  71
  72     TIDSortedNodeSet nodes;
  73     if ( theMesh ) {
  74       SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
  75       while ( nIt->more() )
  76         nodes.insert( nodes.end(), nIt->next() );
  77     }
  78     else if ( theElemIt )
  79     {
  80       while ( theElemIt->more() )
  81       {
  82         const SMDS_MeshElement* e = theElemIt->next();
  83         nodes.insert( e->begin_nodes(), e->end_nodes() );
  84       }
  85     }
  86     myOctreeNode = new SMESH_OctreeNode(nodes) ;
  87
  88     // get max size of a leaf box
  89     SMESH_OctreeNode* tree = myOctreeNode;
  90     while ( !tree->isLeaf() )
  91     {
  92       SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
  93       if ( cIt->more() )
  94         tree = cIt->next();
  95     }
  96     myHalfLeafSize = tree->maxSize() / 2.;
  97   }
  98
  99   //---------------------------------------------------------------------
 100   /*!
 101    * \brief Move node and update myOctreeNode accordingly
 102    */
 103   void MoveNode( const SMDS_MeshNode* node, const gp_Pnt& toPnt )
 104   {
 105     myOctreeNode->UpdateByMoveNode( node, toPnt );
 106     myMesh->MoveNode( node, toPnt.X(), toPnt.Y(), toPnt.Z() );
 107   }
 108
 109   //---------------------------------------------------------------------
 110   /*!
 111    * \brief Do it's job
 112    */
 113   const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt )
 114   {
 115     std::map<double, const SMDS_MeshNode*> dist2Nodes;
 116     myOctreeNode->NodesAround( thePnt.Coord(), dist2Nodes, myHalfLeafSize );
 117     if ( !dist2Nodes.empty() )
 118       return dist2Nodes.begin()->second;
 119
 120     std::vector<const SMDS_MeshNode*> nodes;
 121     //myOctreeNode->NodesAround( &tgtNode, &nodes, myHalfLeafSize );
 122
 123     double minSqDist = DBL_MAX;
 124     if ( nodes.empty() )  // get all nodes of OctreeNode's closest to thePnt
 125     {
 126       // sort leafs by their distance from thePnt
 127       typedef std::multimap< double, SMESH_OctreeNode* > TDistTreeMap;
 128       TDistTreeMap treeMap;
 129       std::list< SMESH_OctreeNode* > treeList;
 130       std::list< SMESH_OctreeNode* >::iterator trIt;
 131       treeList.push_back( myOctreeNode );
 132
 133       gp_XYZ pointNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
 134       bool pointInside = myOctreeNode->isInside( pointNode, myHalfLeafSize );
 135       for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt)
 136       {
 137         SMESH_OctreeNode* tree = *trIt;
 138         if ( !tree->isLeaf() ) // put children to the queue
 139         {
 140           if ( pointInside && !tree->isInside( pointNode, myHalfLeafSize )) continue;
 141           SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
 142           while ( cIt->more() )
 143             treeList.push_back( cIt->next() );
 144         }
 145         else if ( tree->NbNodes() ) // put a tree to the treeMap
 146         {
 147           const Bnd_B3d& box = *tree->getBox();
 148           double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
 149           treeMap.insert( std::make_pair( sqDist, tree ));
 150         }
 151       }
 152       // find distance after which there is no sense to check tree's
 153       double sqLimit = DBL_MAX;
 154       TDistTreeMap::iterator sqDist_tree = treeMap.begin();
 155       if ( treeMap.size() > 5 ) {
 156         SMESH_OctreeNode* closestTree = sqDist_tree->second;
 157         const Bnd_B3d& box = *closestTree->getBox();
 158         double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
 159         sqLimit = limit * limit;
 160       }
 161       // get all nodes from trees
 162       for ( ; sqDist_tree != treeMap.end(); ++sqDist_tree) {
 163         if ( sqDist_tree->first > sqLimit )
 164           break;
 165         SMESH_OctreeNode* tree = sqDist_tree->second;
 166         tree->AllNodesAround( tree->GetNodeIterator()->next(), &nodes );
 167       }
 168     }
 169     // find closest among nodes
 170     minSqDist = DBL_MAX;
 171     const SMDS_MeshNode* closestNode = 0;
 172     for ( size_t i = 0; i < nodes.size(); ++i )
 173     {
 174       double sqDist = thePnt.SquareDistance( SMESH_NodeXYZ( nodes[ i ]));
 175       if ( minSqDist > sqDist ) {
 176         closestNode = nodes[ i ];
 177         minSqDist = sqDist;
 178       }
 179     }
 180     return closestNode;
 181   }
 182
 183   //---------------------------------------------------------------------
 184   /*!
 185    * \brief Finds nodes located within a tolerance near a point
 186    */
 187   int FindNearPoint(const gp_Pnt&                        point,
 188                     const double                         tolerance,
 189                     std::vector< const SMDS_MeshNode* >& foundNodes)
 190   {
 191     myOctreeNode->NodesAround( point.Coord(), foundNodes, tolerance );
 192     return foundNodes.size();
 193   }
 194
 195   //---------------------------------------------------------------------
 196   /*!
 197    * \brief Destructor
 198    */
 199   ~SMESH_NodeSearcherImpl() { delete myOctreeNode; }
 200
 201   //---------------------------------------------------------------------
 202   /*!
 203    * \brief Return the node tree
 204    */
 205   const SMESH_OctreeNode* getTree() const { return myOctreeNode; }
 206
 207 private:
 208   SMESH_OctreeNode* myOctreeNode;
 209   SMDS_Mesh*        myMesh;
 210   double            myHalfLeafSize; // max size of a leaf box
 211 };
 212
 213 // ========================================================================
 214 namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 215 {
 216   const int MaxNbElemsInLeaf = 10; // maximal number of elements in a leaf of tree
 217   const int MaxLevel         = 7;  // maximal tree height -> nb terminal boxes: 8^7 = 2097152
 218   const double NodeRadius = 1e-9;  // to enlarge bnd box of element
 219
 220   //=======================================================================
 221   /*!
 222    * \brief Octal tree of bounding boxes of elements
 223    */
 224   //=======================================================================
 225
 226   class ElementBndBoxTree : public SMESH_Octree
 227   {
 228   public:
 229
 230     typedef boost::container::flat_set< const SMDS_MeshElement*, TIDCompare > TElemSeq;
 231
 232     ElementBndBoxTree(const SMDS_Mesh&     mesh,
 233                       SMDSAbs_ElementType  elemType,
 234                       SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(),
 235                       double               tolerance = NodeRadius );
 236     void getElementsNearPoint( const gp_Pnt& point, TElemSeq& foundElems );
 237     void getElementsNearLine ( const gp_Ax1& line,  TElemSeq& foundElems );
 238     void getElementsInBox    ( const Bnd_B3d& box,  TElemSeq& foundElems );
 239     void getElementsInSphere ( const gp_XYZ& center, const double radius, TElemSeq& foundElems );
 240     ElementBndBoxTree* getLeafAtPoint( const gp_XYZ& point );
 241
 242   protected:
 243     ElementBndBoxTree() {}
 244     SMESH_Octree* newChild() const { return new ElementBndBoxTree; }
 245     void          buildChildrenData();
 246     Bnd_B3d*      buildRootBox();
 247   private:
 248     //!< Bounding box of element
 249     struct ElementBox : public Bnd_B3d
 250     {
 251       const SMDS_MeshElement* _element;
 252       void init(const SMDS_MeshElement* elem, double tolerance);
 253     };
 254     std::vector< ElementBox* > _elements;
 255
 256     typedef ObjectPool< ElementBox > TElementBoxPool;
 257
 258     //!< allocator of ElementBox's and SMESH_TreeLimit
 259     struct LimitAndPool : public SMESH_TreeLimit
 260     {
 261       TElementBoxPool _elBoPool;
 262       LimitAndPool():SMESH_TreeLimit( MaxLevel, /*minSize=*/0. ) {}
 263     };
 264     LimitAndPool* getLimitAndPool() const
 265     {
 266       SMESH_TreeLimit* limitAndPool = const_cast< SMESH_TreeLimit* >( myLimit );
 267       return static_cast< LimitAndPool* >( limitAndPool );
 268     }
 269   };
 270
 271   //================================================================================
 272   /*!
 273    * \brief ElementBndBoxTree creation
 274    */
 275   //================================================================================
 276
 277   ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh&     mesh,
 278                                        SMDSAbs_ElementType  elemType,
 279                                        SMDS_ElemIteratorPtr theElemIt,
 280                                        double               tolerance)
 281     :SMESH_Octree( new LimitAndPool() )
 282   {
 283     int nbElems = mesh.GetMeshInfo().NbElements( elemType );
 284     _elements.reserve( nbElems );
 285
 286     TElementBoxPool& elBoPool = getLimitAndPool()->_elBoPool;
 287
 288     SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType );
 289     while ( elemIt->more() )
 290     {
 291       ElementBox* eb = elBoPool.getNew();
 292       eb->init( elemIt->next(), tolerance );
 293       _elements.push_back( eb );
 294     }
 295     compute();
 296   }
 297
 298   //================================================================================
 299   /*!
 300    * \brief Return the maximal box
 301    */
 302   //================================================================================
 303
 304   Bnd_B3d* ElementBndBoxTree::buildRootBox()
 305   {
 306     Bnd_B3d* box = new Bnd_B3d;
 307     for ( size_t i = 0; i < _elements.size(); ++i )
 308       box->Add( *_elements[i] );
 309     return box;
 310   }
 311
 312   //================================================================================
 313   /*!
 314    * \brief Redistrubute element boxes among children
 315    */
 316   //================================================================================
 317
 318   void ElementBndBoxTree::buildChildrenData()
 319   {
 320     for ( size_t i = 0; i < _elements.size(); ++i )
 321     {
 322       for (int j = 0; j < 8; j++)
 323       {
 324         if ( !_elements[i]->IsOut( *myChildren[j]->getBox() ))
 325           ((ElementBndBoxTree*)myChildren[j])->_elements.push_back( _elements[i]);
 326       }
 327     }
 328     //_size = _elements.size();
 329     SMESHUtils::FreeVector( _elements ); // = _elements.clear() + free memory
 330
 331     for (int j = 0; j < 8; j++)
 332     {
 333       ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j]);
 334       if ((int) child->_elements.size() <= MaxNbElemsInLeaf )
 335         child->myIsLeaf = true;
 336
 337       if ( child->isLeaf() && child->_elements.capacity() > child->_elements.size() )
 338         SMESHUtils::CompactVector( child->_elements );
 339     }
 340   }
 341
 342   //================================================================================
 343   /*!
 344    * \brief Return elements which can include the point
 345    */
 346   //================================================================================
 347
 348   void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt& point, TElemSeq& foundElems)
 349   {
 350     if ( getBox()->IsOut( point.XYZ() ))
 351       return;
 352
 353     if ( isLeaf() )
 354     {
 355       for ( size_t i = 0; i < _elements.size(); ++i )
 356         if ( !_elements[i]->IsOut( point.XYZ() ))
 357           foundElems.insert( _elements[i]->_element );
 358     }
 359     else
 360     {
 361       for (int i = 0; i < 8; i++)
 362         ((ElementBndBoxTree*) myChildren[i])->getElementsNearPoint( point, foundElems );
 363     }
 364   }
 365
 366   //================================================================================
 367   /*!
 368    * \brief Return elements which can be intersected by the line
 369    */
 370   //================================================================================
 371
 372   void ElementBndBoxTree::getElementsNearLine( const gp_Ax1& line, TElemSeq& foundElems )
 373   {
 374     if ( getBox()->IsOut( line ))
 375       return;
 376
 377     if ( isLeaf() )
 378     {
 379       for ( size_t i = 0; i < _elements.size(); ++i )
 380         if ( !_elements[i]->IsOut( line ) )
 381           foundElems.insert( _elements[i]->_element );
 382     }
 383     else
 384     {
 385       for (int i = 0; i < 8; i++)
 386         ((ElementBndBoxTree*) myChildren[i])->getElementsNearLine( line, foundElems );
 387     }
 388   }
 389
 390   //================================================================================
 391   /*!
 392    * \brief Return elements from leaves intersecting the sphere
 393    */
 394   //================================================================================
 395
 396   void ElementBndBoxTree::getElementsInSphere ( const gp_XYZ& center,
 397                                                 const double  radius,
 398                                                 TElemSeq&     foundElems)
 399   {
 400     if ( getBox()->IsOut( center, radius ))
 401       return;
 402
 403     if ( isLeaf() )
 404     {
 405       for ( size_t i = 0; i < _elements.size(); ++i )
 406         if ( !_elements[i]->IsOut( center, radius ))
 407           foundElems.insert( _elements[i]->_element );
 408     }
 409     else
 410     {
 411       for (int i = 0; i < 8; i++)
 412         ((ElementBndBoxTree*) myChildren[i])->getElementsInSphere( center, radius, foundElems );
 413     }
 414   }
 415
 416   //================================================================================
 417   /*!
 418    * \brief Return elements from leaves intersecting the box
 419    */
 420   //================================================================================
 421
 422   void ElementBndBoxTree::getElementsInBox( const Bnd_B3d& box,  TElemSeq& foundElems )
 423   {
 424     if ( getBox()->IsOut( box ))
 425       return;
 426
 427     if ( isLeaf() )
 428     {
 429       for ( size_t i = 0; i < _elements.size(); ++i )
 430         if ( !_elements[i]->IsOut( box ))
 431           foundElems.insert( _elements[i]->_element );
 432     }
 433     else
 434     {
 435       for (int i = 0; i < 8; i++)
 436         ((ElementBndBoxTree*) myChildren[i])->getElementsInBox( box, foundElems );
 437     }
 438   }
 439
 440   //================================================================================
 441   /*!
 442    * \brief Return a leaf including a point
 443    */
 444   //================================================================================
 445
 446   ElementBndBoxTree* ElementBndBoxTree::getLeafAtPoint( const gp_XYZ& point )
 447   {
 448     if ( getBox()->IsOut( point ))
 449       return 0;
 450
 451     if ( isLeaf() )
 452     {
 453       return this;
 454     }
 455     else
 456     {
 457       for (int i = 0; i < 8; i++)
 458         if ( ElementBndBoxTree* l = ((ElementBndBoxTree*) myChildren[i])->getLeafAtPoint( point ))
 459           return l;
 460     }
 461     return 0;
 462   }
 463
 464   //================================================================================
 465   /*!
 466    * \brief Construct the element box
 467    */
 468   //================================================================================
 469
 470   void ElementBndBoxTree::ElementBox::init(const SMDS_MeshElement* elem, double tolerance)
 471   {
 472     _element  = elem;
 473     SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
 474     while ( nIt->more() )
 475       Add( SMESH_NodeXYZ( nIt->next() ));
 476     Enlarge( tolerance );
 477   }
 478
 479 } // namespace
 480
 481 //=======================================================================
 482 /*!
 483  * \brief Implementation of search for the elements by point and
 484  *        of classification of point in 2D mesh
 485  */
 486 //=======================================================================
 487
 488 SMESH_ElementSearcher::~SMESH_ElementSearcher()
 489 {
 490 }
 491
 492 struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
 493 {
 494   SMDS_Mesh*                        _mesh;
 495   SMDS_ElemIteratorPtr              _meshPartIt;
 496   ElementBndBoxTree*                _ebbTree      [SMDSAbs_NbElementTypes];
 497   int                               _ebbTreeHeight[SMDSAbs_NbElementTypes];
 498   SMESH_NodeSearcherImpl*           _nodeSearcher;
 499   SMDSAbs_ElementType               _elementType;
 500   double                            _tolerance;
 501   bool                              _outerFacesFound;
 502   std::set<const SMDS_MeshElement*> _outerFaces; // empty means "no internal faces at all"
 503
 504   SMESH_ElementSearcherImpl( SMDS_Mesh&           mesh,
 505                              double               tol=-1,
 506                              SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr())
 507     : _mesh(&mesh),_meshPartIt(elemIt),_nodeSearcher(0),_tolerance(tol),_outerFacesFound(false)
 508   {
 509     for ( int i = 0; i < SMDSAbs_NbElementTypes; ++i )
 510     {
 511       _ebbTree[i] = NULL;
 512       _ebbTreeHeight[i] = -1;
 513     }
 514     _elementType = SMDSAbs_All;
 515   }
 516   virtual ~SMESH_ElementSearcherImpl()
 517   {
 518     for ( int i = 0; i < SMDSAbs_NbElementTypes; ++i )
 519     {
 520       delete _ebbTree[i]; _ebbTree[i] = NULL;
 521     }
 522     if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0;
 523   }
 524   virtual int FindElementsByPoint(const gp_Pnt&                           point,
 525                                   SMDSAbs_ElementType                     type,
 526                                   std::vector< const SMDS_MeshElement* >& foundElements);
 527   virtual TopAbs_State GetPointState(const gp_Pnt& point);
 528   virtual const SMDS_MeshElement* FindClosestTo( const gp_Pnt&       point,
 529                                                  SMDSAbs_ElementType type );
 530
 531   virtual void GetElementsNearLine( const gp_Ax1&                           line,
 532                                     SMDSAbs_ElementType                     type,
 533                                     std::vector< const SMDS_MeshElement* >& foundElems);
 534   virtual void GetElementsInSphere( const gp_XYZ&                           center,
 535                                     const double                            radius,
 536                                     SMDSAbs_ElementType                     type,
 537                                     std::vector< const SMDS_MeshElement* >& foundElems);
 538   virtual void GetElementsInBox( const Bnd_B3d&                          box,
 539                                  SMDSAbs_ElementType                     type,
 540                                  std::vector< const SMDS_MeshElement* >& foundElems);
 541   virtual gp_XYZ Project(const gp_Pnt&            point,
 542                          SMDSAbs_ElementType      type,
 543                          const SMDS_MeshElement** closestElem);
 544   double getTolerance();
 545   bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
 546                             const double tolerance, double & param);
 547   void findOuterBoundary(const SMDS_MeshElement* anyOuterFace);
 548   bool isOuterBoundary(const SMDS_MeshElement* face) const
 549   {
 550     return _outerFaces.empty() || _outerFaces.count(face);
 551   }
 552   int getTreeHeight()
 553   {
 554     if ( _ebbTreeHeight[ _elementType ] < 0 )
 555       _ebbTreeHeight[ _elementType ] = _ebbTree[ _elementType ]->getHeight();
 556     return _ebbTreeHeight[ _elementType ];
 557   }
 558
 559   struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
 560   {
 561     const SMDS_MeshElement* _face;
 562     gp_Vec                  _faceNorm;
 563     bool                    _coincides; //!< the line lays in face plane
 564     TInters(const SMDS_MeshElement* face, const gp_Vec& faceNorm, bool coinc=false)
 565       : _face(face), _faceNorm( faceNorm ), _coincides( coinc ) {}
 566   };
 567   struct TFaceLink //!< link and faces sharing it (used in findOuterBoundary())
 568   {
 569     SMESH_TLink      _link;
 570     TIDSortedElemSet _faces;
 571     TFaceLink( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshElement* face)
 572       : _link( n1, n2 ), _faces( &face, &face + 1) {}
 573   };
 574 };
 575
 576 ostream& operator<< (ostream& out, const SMESH_ElementSearcherImpl::TInters& i)
 577 {
 578   return out << "TInters(face=" << ( i._face ? i._face->GetID() : 0)
 579              << ", _coincides="<<i._coincides << ")";
 580 }
 581
 582 //=======================================================================
 583 /*!
 584  * \brief define tolerance for search
 585  */
 586 //=======================================================================
 587
 588 double SMESH_ElementSearcherImpl::getTolerance()
 589 {
 590   if ( _tolerance < 0 )
 591   {
 592     const SMDS_MeshInfo& meshInfo = _mesh->GetMeshInfo();
 593
 594     _tolerance = 0;
 595     if ( _nodeSearcher && meshInfo.NbNodes() > 1 )
 596     {
 597       double boxSize = _nodeSearcher->getTree()->maxSize();
 598       _tolerance = 1e-8 * boxSize/* / meshInfo.NbNodes()*/;
 599     }
 600     else if ( _ebbTree[_elementType] && meshInfo.NbElements() > 0 )
 601     {
 602       double boxSize = _ebbTree[_elementType]->maxSize();
 603       _tolerance = 1e-8 * boxSize/* / meshInfo.NbElements()*/;
 604     }
 605     if ( _tolerance == 0 )
 606     {
 607       // define tolerance by size of a most complex element
 608       int complexType = SMDSAbs_Volume;
 609       while ( complexType > SMDSAbs_All &&
 610               meshInfo.NbElements( SMDSAbs_ElementType( complexType )) < 1 )
 611         --complexType;
 612       if ( complexType == SMDSAbs_All ) return 0; // empty mesh
 613       double elemSize;
 614       if ( complexType == int( SMDSAbs_Node ))
 615       {
 616         SMDS_NodeIteratorPtr nodeIt = _mesh->nodesIterator();
 617         elemSize = 1;
 618         if ( meshInfo.NbNodes() > 2 )
 619           elemSize = SMESH_TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() );
 620       }
 621       else
 622       {
 623         SMDS_ElemIteratorPtr  elemIt = _mesh->elementsIterator( SMDSAbs_ElementType( complexType ));
 624         const SMDS_MeshElement* elem = elemIt->next();
 625         SMDS_ElemIteratorPtr  nodeIt = elem->nodesIterator();
 626         SMESH_TNodeXYZ n1( nodeIt->next() );
 627         elemSize = 0;
 628         while ( nodeIt->more() )
 629         {
 630           double dist = n1.Distance( static_cast<const SMDS_MeshNode*>( nodeIt->next() ));
 631           elemSize = std::max( dist, elemSize );
 632         }
 633       }
 634       _tolerance = 1e-4 * elemSize;
 635     }
 636   }
 637   return _tolerance;
 638 }
 639
 640 //================================================================================
 641 /*!
 642  * \brief Find intersection of the line and an edge of face and return parameter on line
 643  */
 644 //================================================================================
 645
 646 bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin&           line,
 647                                                      const SMDS_MeshElement* face,
 648                                                      const double            tol,
 649                                                      double &                param)
 650 {
 651   int nbInts = 0;
 652   param = 0;
 653
 654   GeomAPI_ExtremaCurveCurve anExtCC;
 655   Handle(Geom_Curve) lineCurve = new Geom_Line( line );
 656
 657   int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
 658   for ( int i = 0; i < nbNodes && nbInts < 2; ++i )
 659   {
 660     GC_MakeSegment edge( SMESH_TNodeXYZ( face->GetNode( i )),
 661                          SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) ));
 662     anExtCC.Init( lineCurve, edge.Value() );
 663     if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol)
 664     {
 665       Standard_Real pl, pe;
 666       anExtCC.LowerDistanceParameters( pl, pe );
 667       param += pl;
 668       if ( ++nbInts == 2 )
 669         break;
 670     }
 671   }
 672   if ( nbInts > 0 ) param /= nbInts;
 673   return nbInts > 0;
 674 }
 675 //================================================================================
 676 /*!
 677  * \brief Find all faces belonging to the outer boundary of mesh
 678  */
 679 //================================================================================
 680
 681 void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerFace)
 682 {
 683   if ( _outerFacesFound ) return;
 684
 685   // Collect all outer faces by passing from one outer face to another via their links
 686   // and BTW find out if there are internal faces at all.
 687
 688   // checked links and links where outer boundary meets internal one
 689   std::set< SMESH_TLink > visitedLinks, seamLinks;
 690
 691   // links to treat with already visited faces sharing them
 692   std::list < TFaceLink > startLinks;
 693
 694   // load startLinks with the first outerFace
 695   startLinks.push_back( TFaceLink( outerFace->GetNode(0), outerFace->GetNode(1), outerFace));
 696   _outerFaces.insert( outerFace );
 697
 698   TIDSortedElemSet emptySet;
 699   while ( !startLinks.empty() )
 700   {
 701     const SMESH_TLink& link  = startLinks.front()._link;
 702     TIDSortedElemSet&  faces = startLinks.front()._faces;
 703
 704     outerFace = *faces.begin();
 705     // find other faces sharing the link
 706     const SMDS_MeshElement* f;
 707     while (( f = SMESH_MeshAlgos::FindFaceInSet(link.node1(), link.node2(), emptySet, faces )))
 708       faces.insert( f );
 709
 710     // select another outer face among the found
 711     const SMDS_MeshElement* outerFace2 = 0;
 712     if ( faces.size() == 2 )
 713     {
 714       outerFace2 = (outerFace == *faces.begin() ? *faces.rbegin() : *faces.begin());
 715     }
 716     else if ( faces.size() > 2 )
 717     {
 718       seamLinks.insert( link );
 719
 720       // link direction within the outerFace
 721       gp_Vec n1n2( SMESH_TNodeXYZ( link.node1()),
 722                    SMESH_TNodeXYZ( link.node2()));
 723       int i1 = outerFace->GetNodeIndex( link.node1() );
 724       int i2 = outerFace->GetNodeIndex( link.node2() );
 725       bool rev = ( abs(i2-i1) == 1 ? i1 > i2 : i2 > i1 );
 726       if ( rev ) n1n2.Reverse();
 727       // outerFace normal
 728       gp_XYZ ofNorm, fNorm;
 729       if ( SMESH_MeshAlgos::FaceNormal( outerFace, ofNorm, /*normalized=*/false ))
 730       {
 731         // direction from the link inside outerFace
 732         gp_Vec dirInOF = gp_Vec( ofNorm ) ^ n1n2;
 733         // sort all other faces by angle with the dirInOF
 734         std::map< double, const SMDS_MeshElement* > angle2Face;
 735         std::set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin();
 736         for ( ; face != faces.end(); ++face )
 737         {
 738           if ( *face == outerFace ) continue;
 739           if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false ))
 740             continue;
 741           gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2;
 742           double angle = dirInOF.AngleWithRef( dirInF, n1n2 );
 743           if ( angle < 0 ) angle += 2. * M_PI;
 744           angle2Face.insert( std::make_pair( angle, *face ));
 745         }
 746         if ( !angle2Face.empty() )
 747           outerFace2 = angle2Face.begin()->second;
 748       }
 749     }
 750     // store the found outer face and add its links to continue searching from
 751     if ( outerFace2 )
 752     {
 753       _outerFaces.insert( outerFace2 );
 754       int nbNodes = outerFace2->NbCornerNodes();
 755       for ( int i = 0; i < nbNodes; ++i )
 756       {
 757         SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes));
 758         if ( visitedLinks.insert( link2 ).second )
 759           startLinks.push_back( TFaceLink( link2.node1(), link2.node2(), outerFace2 ));
 760       }
 761     }
 762     startLinks.pop_front();
 763   }
 764   _outerFacesFound = true;
 765
 766   if ( !seamLinks.empty() )
 767   {
 768     // There are internal boundaries touching the outher one,
 769     // find all faces of internal boundaries in order to find
 770     // faces of boundaries of holes, if any.
 771
 772   }
 773   else
 774   {
 775     _outerFaces.clear();
 776   }
 777 }
 778
 779 //=======================================================================
 780 /*!
 781  * \brief Find elements of given type where the given point is IN or ON.
 782  *        Returns nb of found elements and elements them-selves.
 783  *
 784  * 'ALL' type means elements of any type excluding nodes, balls and 0D elements
 785  */
 786 //=======================================================================
 787
 788 int SMESH_ElementSearcherImpl::
 789 FindElementsByPoint(const gp_Pnt&                           point,
 790                     SMDSAbs_ElementType                     type,
 791                     std::vector< const SMDS_MeshElement* >& foundElements)
 792 {
 793   foundElements.clear();
 794   _elementType = type;
 795
 796   double tolerance = getTolerance();
 797
 798   // =================================================================================
 799   if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement || type == SMDSAbs_Ball)
 800   {
 801     if ( !_nodeSearcher )
 802     {
 803       if ( _meshPartIt )
 804         _nodeSearcher = new SMESH_NodeSearcherImpl( 0, _meshPartIt );
 805       else
 806         _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
 807     }
 808     std::vector< const SMDS_MeshNode* > foundNodes;
 809     _nodeSearcher->FindNearPoint( point, tolerance, foundNodes );
 810
 811     if ( type == SMDSAbs_Node )
 812     {
 813       foundElements.assign( foundNodes.begin(), foundNodes.end() );
 814     }
 815     else
 816     {
 817       for ( size_t i = 0; i < foundNodes.size(); ++i )
 818       {
 819         SMDS_ElemIteratorPtr elemIt = foundNodes[i]->GetInverseElementIterator( type );
 820         while ( elemIt->more() )
 821           foundElements.push_back( elemIt->next() );
 822       }
 823     }
 824   }
 825   // =================================================================================
 826   else // elements more complex than 0D
 827   {
 828     if ( !_ebbTree[type] )
 829     {
 830       _ebbTree[_elementType] = new ElementBndBoxTree( *_mesh, type, _meshPartIt, tolerance );
 831     }
 832     ElementBndBoxTree::TElemSeq suspectElems;
 833     _ebbTree[ type ]->getElementsNearPoint( point, suspectElems );
 834     ElementBndBoxTree::TElemSeq::iterator elem = suspectElems.begin();
 835     for ( ; elem != suspectElems.end(); ++elem )
 836       if ( !SMESH_MeshAlgos::IsOut( *elem, point, tolerance ))
 837         foundElements.push_back( *elem );
 838   }
 839   return foundElements.size();
 840 }
 841
 842 //=======================================================================
 843 /*!
 844  * \brief Find an element of given type most close to the given point
 845  *
 846  * WARNING: Only face search is implemeneted so far
 847  */
 848 //=======================================================================
 849
 850 const SMDS_MeshElement*
 851 SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
 852                                           SMDSAbs_ElementType type )
 853 {
 854   const SMDS_MeshElement* closestElem = 0;
 855   _elementType = type;
 856
 857   if ( type == SMDSAbs_Face ||
 858        type == SMDSAbs_Volume ||
 859        type == SMDSAbs_Edge )
 860   {
 861     ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
 862     if ( !ebbTree )
 863       ebbTree = new ElementBndBoxTree( *_mesh, type, _meshPartIt );
 864
 865     ElementBndBoxTree::TElemSeq suspectElems;
 866     ebbTree->getElementsNearPoint( point, suspectElems );
 867
 868     if ( suspectElems.empty() && ebbTree->maxSize() > 0 )
 869     {
 870       gp_Pnt boxCenter = 0.5 * ( ebbTree->getBox()->CornerMin() +
 871                                  ebbTree->getBox()->CornerMax() );
 872       double radius = -1;
 873       if ( ebbTree->getBox()->IsOut( point.XYZ() ))
 874         radius = point.Distance( boxCenter ) - 0.5 * ebbTree->maxSize();
 875       if ( radius < 0 )
 876         radius = ebbTree->maxSize() / pow( 2., getTreeHeight()) / 2;
 877       while ( suspectElems.empty() )
 878       {
 879         ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems );
 880         radius *= 1.1;
 881       }
 882     }
 883     double minDist = std::numeric_limits<double>::max();
 884     std::multimap< double, const SMDS_MeshElement* > dist2face;
 885     ElementBndBoxTree::TElemSeq::iterator elem = suspectElems.begin();
 886     for ( ; elem != suspectElems.end(); ++elem )
 887     {
 888       double dist = SMESH_MeshAlgos::GetDistance( *elem, point );
 889       if ( dist < minDist + 1e-10)
 890       {
 891         minDist = dist;
 892         dist2face.insert( dist2face.begin(), std::make_pair( dist, *elem ));
 893       }
 894     }
 895     if ( !dist2face.empty() )
 896     {
 897       std::multimap< double, const SMDS_MeshElement* >::iterator d2f = dist2face.begin();
 898       closestElem = d2f->second;
 899       // if there are several elements at the same distance, select one
 900       // with GC closest to the point
 901       typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
 902       double minDistToGC = 0;
 903       for ( ++d2f; d2f != dist2face.end() && fabs( d2f->first - minDist ) < 1e-10; ++d2f )
 904       {
 905         if ( minDistToGC == 0 )
 906         {
 907           gp_XYZ gc(0,0,0);
 908           gc = accumulate( TXyzIterator(closestElem->nodesIterator()),
 909                            TXyzIterator(), gc ) / closestElem->NbNodes();
 910           minDistToGC = point.SquareDistance( gc );
 911         }
 912         gp_XYZ gc(0,0,0);
 913         gc = accumulate( TXyzIterator( d2f->second->nodesIterator()),
 914                          TXyzIterator(), gc ) / d2f->second->NbNodes();
 915         double d = point.SquareDistance( gc );
 916         if ( d < minDistToGC )
 917         {
 918           minDistToGC = d;
 919           closestElem = d2f->second;
 920         }
 921       }
 922       // cout << "FindClosestTo( " <<point.X()<<", "<<point.Y()<<", "<<point.Z()<<" ) FACE "
 923       //      <<closestElem->GetID() << " DIST " << minDist << endl;
 924     }
 925   }
 926   else
 927   {
 928     // NOT IMPLEMENTED SO FAR
 929   }
 930   return closestElem;
 931 }
 932
 933
 934 //================================================================================
 935 /*!
 936  * \brief Classify the given point in the closed 2D mesh
 937  */
 938 //================================================================================
 939
 940 TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
 941 {
 942   _elementType = SMDSAbs_Face;
 943
 944   double tolerance = getTolerance();
 945
 946   ElementBndBoxTree*& ebbTree = _ebbTree[ SMDSAbs_Face ];
 947   if ( !ebbTree )
 948     ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
 949
 950   // Algo: analyse transition of a line starting at the point through mesh boundary;
 951   // try three lines parallel to axis of the coordinate system and perform rough
 952   // analysis. If solution is not clear perform thorough analysis.
 953
 954   const int nbAxes = 3;
 955   gp_Dir axisDir[ nbAxes ] = { gp::DX(), gp::DY(), gp::DZ() };
 956   std::map< double, TInters >   paramOnLine2TInters[ nbAxes ];
 957   std::list< TInters > tangentInters[ nbAxes ]; // of faces whose plane includes the line
 958   std::multimap< int, int > nbInt2Axis; // to find the simplest case
 959   for ( int axis = 0; axis < nbAxes; ++axis )
 960   {
 961     gp_Ax1 lineAxis( point, axisDir[axis]);
 962     gp_Lin line    ( lineAxis );
 963
 964     ElementBndBoxTree::TElemSeq suspectFaces; // faces possibly intersecting the line
 965     ebbTree->getElementsNearLine( lineAxis, suspectFaces );
 966
 967     // Intersect faces with the line
 968
 969     std::map< double, TInters > & u2inters = paramOnLine2TInters[ axis ];
 970     ElementBndBoxTree::TElemSeq::iterator face = suspectFaces.begin();
 971     for ( ; face != suspectFaces.end(); ++face )
 972     {
 973       // get face plane
 974       gp_XYZ fNorm;
 975       if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false)) continue;
 976       gp_Pln facePlane( SMESH_TNodeXYZ( (*face)->GetNode(0)), fNorm );
 977
 978       // perform intersection
 979       IntAna_IntConicQuad intersection( line, IntAna_Quadric( facePlane ));
 980       if ( !intersection.IsDone() )
 981         continue;
 982       if ( intersection.IsInQuadric() )
 983       {
 984         tangentInters[ axis ].push_back( TInters( *face, fNorm, true ));
 985       }
 986       else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 )
 987       {
 988         double tol = 1e-4 * Sqrt( fNorm.Modulus() );
 989         gp_Pnt intersectionPoint = intersection.Point(1);
 990         if ( !SMESH_MeshAlgos::IsOut( *face, intersectionPoint, tol ))
 991           u2inters.insert( std::make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm )));
 992       }
 993     }
 994     // Analyse intersections roughly
 995
 996     int nbInter = u2inters.size();
 997     if ( nbInter == 0 )
 998       return TopAbs_OUT;
 999
1000     double f = u2inters.begin()->first, l = u2inters.rbegin()->first;
1001     if ( nbInter == 1 ) // not closed mesh
1002       return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
1003
1004     if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
1005       return TopAbs_ON;
1006
1007     if ( (f<0) == (l<0) )
1008       return TopAbs_OUT;
1009
1010     int nbIntBeforePoint = std::distance( u2inters.begin(), u2inters.lower_bound(0));
1011     int nbIntAfterPoint  = nbInter - nbIntBeforePoint;
1012     if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 )
1013       return TopAbs_IN;
1014
1015     nbInt2Axis.insert( std::make_pair( std::min( nbIntBeforePoint, nbIntAfterPoint ), axis ));
1016
1017     if ( _outerFacesFound ) break; // pass to thorough analysis
1018
1019   } // three attempts - loop on CS axes
1020
1021   // Analyse intersections thoroughly.
1022   // We make two loops maximum, on the first one we only exclude touching intersections,
1023   // on the second, if situation is still unclear, we gather and use information on
1024   // position of faces (internal or outer). If faces position is already gathered,
1025   // we make the second loop right away.
1026
1027   for ( int hasPositionInfo = _outerFacesFound; hasPositionInfo < 2; ++hasPositionInfo )
1028   {
1029     std::multimap< int, int >::const_iterator nb_axis = nbInt2Axis.begin();
1030     for ( ; nb_axis != nbInt2Axis.end(); ++nb_axis )
1031     {
1032       int axis = nb_axis->second;
1033       std::map< double, TInters > & u2inters = paramOnLine2TInters[ axis ];
1034
1035       gp_Ax1 lineAxis( point, axisDir[axis]);
1036       gp_Lin line    ( lineAxis );
1037
1038       // add tangent intersections to u2inters
1039       double param;
1040       std::list< TInters >::const_iterator tgtInt = tangentInters[ axis ].begin();
1041       for ( ; tgtInt != tangentInters[ axis ].end(); ++tgtInt )
1042         if ( getIntersParamOnLine( line, tgtInt->_face, tolerance, param ))
1043           u2inters.insert( std::make_pair( param, *tgtInt ));
1044       tangentInters[ axis ].clear();
1045
1046       // Count intersections before and after the point excluding touching ones.
1047       // If hasPositionInfo we count intersections of outer boundary only
1048
1049       int nbIntBeforePoint = 0, nbIntAfterPoint = 0;
1050       double f = std::numeric_limits<double>::max(), l = -std::numeric_limits<double>::max();
1051       std::map< double, TInters >::iterator u_int1 = u2inters.begin(), u_int2 = u_int1;
1052       bool ok = ! u_int1->second._coincides;
1053       while ( ok && u_int1 != u2inters.end() )
1054       {
1055         double u = u_int1->first;
1056         bool touchingInt = false;
1057         if ( ++u_int2 != u2inters.end() )
1058         {
1059           // skip intersections at the same point (if the line passes through edge or node)
1060           int nbSamePnt = 0;
1061           while ( u_int2 != u2inters.end() && fabs( u_int2->first - u ) < tolerance )
1062           {
1063             ++nbSamePnt;
1064             ++u_int2;
1065           }
1066
1067           // skip tangent intersections
1068           int nbTgt = 0;
1069           if ( u_int2 != u2inters.end() )
1070           {
1071             const SMDS_MeshElement* prevFace = u_int1->second._face;
1072             while ( ok && u_int2->second._coincides )
1073             {
1074               if ( SMESH_MeshAlgos::GetCommonNodes(prevFace , u_int2->second._face).empty() )
1075                 ok = false;
1076               else
1077               {
1078                 nbTgt++;
1079                 u_int2++;
1080                 ok = ( u_int2 != u2inters.end() );
1081               }
1082             }
1083           }
1084           if ( !ok ) break;
1085
1086           // skip intersections at the same point after tangent intersections
1087           if ( nbTgt > 0 )
1088           {
1089             double u2 = u_int2->first;
1090             ++u_int2;
1091             while ( u_int2 != u2inters.end() && fabs( u_int2->first - u2 ) < tolerance )
1092             {
1093               ++nbSamePnt;
1094               ++u_int2;
1095             }
1096           }
1097           // decide if we skipped a touching intersection
1098           if ( nbSamePnt + nbTgt > 0 )
1099           {
1100             double minDot = std::numeric_limits<double>::max(), maxDot = -minDot;
1101             std::map< double, TInters >::iterator u_int = u_int1;
1102             for ( ; u_int != u_int2; ++u_int )
1103             {
1104               if ( u_int->second._coincides ) continue;
1105               double dot = u_int->second._faceNorm * line.Direction();
1106               if ( dot > maxDot ) maxDot = dot;
1107               if ( dot < minDot ) minDot = dot;
1108             }
1109             touchingInt = ( minDot*maxDot < 0 );
1110           }
1111         }
1112         if ( !touchingInt )
1113         {
1114           if ( !hasPositionInfo || isOuterBoundary( u_int1->second._face ))
1115           {
1116             if ( u < 0 )
1117               ++nbIntBeforePoint;
1118             else
1119               ++nbIntAfterPoint;
1120           }
1121           if ( u < f ) f = u;
1122           if ( u > l ) l = u;
1123         }
1124
1125         u_int1 = u_int2; // to next intersection
1126
1127       } // loop on intersections with one line
1128
1129       if ( ok )
1130       {
1131         if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
1132           return TopAbs_ON;
1133
1134         if ( nbIntBeforePoint == 0  || nbIntAfterPoint == 0)
1135           return TopAbs_OUT;
1136
1137         if ( nbIntBeforePoint + nbIntAfterPoint == 1 ) // not closed mesh
1138           return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
1139
1140         if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 )
1141           return TopAbs_IN;
1142
1143         if ( (f<0) == (l<0) )
1144           return TopAbs_OUT;
1145
1146         if ( hasPositionInfo )
1147           return nbIntBeforePoint % 2 ? TopAbs_IN : TopAbs_OUT;
1148       }
1149     } // loop on intersections of the tree lines - thorough analysis
1150
1151     if ( !hasPositionInfo )
1152     {
1153       // gather info on faces position - is face in the outer boundary or not
1154       std::map< double, TInters > & u2inters = paramOnLine2TInters[ 0 ];
1155       findOuterBoundary( u2inters.begin()->second._face );
1156     }
1157
1158   } // two attempts - with and w/o faces position info in the mesh
1159
1160   return TopAbs_UNKNOWN;
1161 }
1162
1163 //=======================================================================
1164 /*!
1165  * \brief Return elements possibly intersecting the line
1166  */
1167 //=======================================================================
1168
1169 void SMESH_ElementSearcherImpl::
1170 GetElementsNearLine( const gp_Ax1&                           line,
1171                      SMDSAbs_ElementType                     type,
1172                      std::vector< const SMDS_MeshElement* >& foundElems)
1173 {
1174   _elementType = type;
1175   ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
1176   if ( !ebbTree )
1177     ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
1178
1179   ElementBndBoxTree::TElemSeq elems;
1180   ebbTree->getElementsNearLine( line, elems );
1181
1182   foundElems.insert( foundElems.end(), elems.begin(), elems.end() );
1183 }
1184
1185 //=======================================================================
1186 /*
1187  * Return elements whose bounding box intersects a sphere
1188  */
1189 //=======================================================================
1190
1191 void SMESH_ElementSearcherImpl::
1192 GetElementsInSphere( const gp_XYZ&                           center,
1193                      const double                            radius,
1194                      SMDSAbs_ElementType                     type,
1195                      std::vector< const SMDS_MeshElement* >& foundElems)
1196 {
1197   _elementType = type;
1198   ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
1199   if ( !ebbTree )
1200     ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
1201
1202   ElementBndBoxTree::TElemSeq elems;
1203   ebbTree->getElementsInSphere( center, radius, elems );
1204
1205   foundElems.insert( foundElems.end(), elems.begin(), elems.end() );
1206 }
1207
1208 //=======================================================================
1209 /*
1210  * Return elements whose bounding box intersects a given bounding box
1211  */
1212 //=======================================================================
1213
1214 void SMESH_ElementSearcherImpl::
1215 GetElementsInBox( const Bnd_B3d&                          box,
1216                   SMDSAbs_ElementType                     type,
1217                   std::vector< const SMDS_MeshElement* >& foundElems)
1218 {
1219   _elementType = type;
1220   ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
1221   if ( !ebbTree )
1222     ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt, getTolerance() );
1223
1224   ElementBndBoxTree::TElemSeq elems;
1225   ebbTree->getElementsInBox( box, elems );
1226
1227   foundElems.insert( foundElems.end(), elems.begin(), elems.end() );
1228 }
1229
1230 //=======================================================================
1231 /*
1232  * \brief Return a projection of a given point to a mesh.
1233  *        Optionally return the closest element
1234  */
1235 //=======================================================================
1236
1237 gp_XYZ SMESH_ElementSearcherImpl::Project(const gp_Pnt&            point,
1238                                           SMDSAbs_ElementType      type,
1239                                           const SMDS_MeshElement** closestElem)
1240 {
1241   _elementType = type;
1242   if ( _mesh->GetMeshInfo().NbElements( _elementType ) == 0 )
1243     throw SALOME_Exception( LOCALIZED( "No elements of given type in the mesh" ));
1244
1245   ElementBndBoxTree*& ebbTree = _ebbTree[ _elementType ];
1246   if ( !ebbTree )
1247     ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
1248
1249   gp_XYZ p = point.XYZ();
1250   ElementBndBoxTree* ebbLeaf = ebbTree->getLeafAtPoint( p );
1251   const Bnd_B3d* box = ebbLeaf ? ebbLeaf->getBox() : ebbTree->getBox();
1252   double radius = ( box->CornerMax() - box->CornerMin() ).Modulus();
1253
1254   ElementBndBoxTree::TElemSeq elems;
1255   ebbTree->getElementsInSphere( p, radius, elems );
1256   while ( elems.empty() )
1257   {
1258     radius *= 1.5;
1259     ebbTree->getElementsInSphere( p, radius, elems );
1260   }
1261   gp_XYZ proj, bestProj;
1262   const SMDS_MeshElement* elem = 0;
1263   double minDist = 2 * radius;
1264   ElementBndBoxTree::TElemSeq::iterator e = elems.begin();
1265   for ( ; e != elems.end(); ++e )
1266   {
1267     double d = SMESH_MeshAlgos::GetDistance( *e, p, &proj );
1268     if ( d < minDist )
1269     {
1270       bestProj = proj;
1271       elem = *e;
1272       minDist = d;
1273     }
1274   }
1275   if ( closestElem ) *closestElem = elem;
1276
1277   return bestProj;
1278 }
1279
1280 //=======================================================================
1281 /*!
1282  * \brief Return true if the point is IN or ON of the element
1283  */
1284 //=======================================================================
1285
1286 bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol )
1287 {
1288   if ( element->GetType() == SMDSAbs_Volume)
1289   {
1290     return SMDS_VolumeTool( element ).IsOut( point.X(), point.Y(), point.Z(), tol );
1291   }
1292
1293   // get ordered nodes
1294
1295   std::vector< SMESH_TNodeXYZ > xyz; xyz.reserve( element->NbNodes()+1 );
1296
1297   SMDS_NodeIteratorPtr nodeIt = element->interlacedNodesIterator();
1298   for ( int i = 0; nodeIt->more(); ++i )
1299     xyz.push_back( SMESH_TNodeXYZ( nodeIt->next() ));
1300
1301   int i, nbNodes = (int) xyz.size(); // central node of biquadratic is missing
1302
1303   if ( element->GetType() == SMDSAbs_Face ) // --------------------------------------------------
1304   {
1305     // compute face normal
1306     gp_Vec faceNorm(0,0,0);
1307     xyz.push_back( xyz.front() );
1308     for ( i = 0; i < nbNodes; ++i )
1309     {
1310       gp_Vec edge1( xyz[i+1], xyz[i]);
1311       gp_Vec edge2( xyz[i+1], xyz[(i+2)%nbNodes] );
1312       faceNorm += edge1 ^ edge2;
1313     }
1314     double fNormSize = faceNorm.Magnitude();
1315     if ( fNormSize <= tol )
1316     {
1317       // degenerated face: point is out if it is out of all face edges
1318       for ( i = 0; i < nbNodes; ++i )
1319       {
1320         SMDS_LinearEdge edge( xyz[i]._node, xyz[i+1]._node );
1321         if ( !IsOut( &edge, point, tol ))
1322           return false;
1323       }
1324       return true;
1325     }
1326     faceNorm /= fNormSize;
1327
1328     // check if the point lays on face plane
1329     gp_Vec n2p( xyz[0], point );
1330     double dot = n2p * faceNorm;
1331     if ( Abs( dot ) > tol ) // not on face plane
1332     {
1333       bool isOut = true;
1334       if ( nbNodes > 3 ) // maybe the face is not planar
1335       {
1336         double elemThick = 0;
1337         for ( i = 1; i < nbNodes; ++i )
1338         {
1339           gp_Vec n2n( xyz[0], xyz[i] );
1340           elemThick = Max( elemThick, Abs( n2n * faceNorm ));
1341         }
1342         isOut = Abs( dot ) > elemThick + tol;
1343       }
1344       if ( isOut )
1345         return isOut;
1346     }
1347
1348     // check if point is out of face boundary:
1349     // define it by closest transition of a ray point->infinity through face boundary
1350     // on the face plane.
1351     // First, find normal of a plane perpendicular to face plane, to be used as a cutting tool
1352     // to find intersections of the ray with the boundary.
1353     gp_Vec ray = n2p;
1354     gp_Vec plnNorm = ray ^ faceNorm;
1355     double n2pSize = plnNorm.Magnitude();
1356     if ( n2pSize <= tol ) return false; // point coincides with the first node
1357     if ( n2pSize * n2pSize > fNormSize * 100 ) return true; // point is very far
1358     plnNorm /= n2pSize;
1359     // for each node of the face, compute its signed distance to the cutting plane
1360     std::vector<double> dist( nbNodes + 1);
1361     for ( i = 0; i < nbNodes; ++i )
1362     {
1363       gp_Vec n2p( xyz[i], point );
1364       dist[i] = n2p * plnNorm;
1365     }
1366     dist.back() = dist.front();
1367     // find the closest intersection
1368     int    iClosest = -1;
1369     double rClosest = 0, distClosest = 1e100;
1370     gp_Pnt pClosest;
1371     for ( i = 0; i < nbNodes; ++i )
1372     {
1373       double r;
1374       if ( fabs( dist[i] ) < tol )
1375         r = 0.;
1376       else if ( fabs( dist[i+1]) < tol )
1377         r = 1.;
1378       else if ( dist[i] * dist[i+1] < 0 )
1379         r = dist[i] / ( dist[i] - dist[i+1] );
1380       else
1381         continue; // no intersection
1382       gp_Pnt pInt = xyz[i] * (1.-r) + xyz[i+1] * r;
1383       gp_Vec p2int( point, pInt);
1384       double intDist = p2int.SquareMagnitude();
1385       if ( intDist < distClosest )
1386       {
1387         iClosest = i;
1388         rClosest = r;
1389         pClosest = pInt;
1390         distClosest = intDist;
1391       }
1392     }
1393     if ( iClosest < 0 )
1394       return true; // no intesections - out
1395
1396     // analyse transition
1397     gp_Vec edge( xyz[iClosest], xyz[iClosest+1] );
1398     gp_Vec edgeNorm = -( edge ^ faceNorm ); // normal to intersected edge pointing out of face
1399     gp_Vec p2int ( point, pClosest );
1400     bool out = (edgeNorm * p2int) < -tol;
1401     if ( rClosest > 0. && rClosest < 1. ) // not node intersection
1402       return out;
1403
1404     // the ray passes through a face node; analyze transition through an adjacent edge
1405     gp_Pnt p1 = xyz[ (rClosest == 0.) ? ((iClosest+nbNodes-1) % nbNodes) : (iClosest+1) ];
1406     gp_Pnt p2 = xyz[ (rClosest == 0.) ? iClosest : ((iClosest+2) % nbNodes) ];
1407     gp_Vec edgeAdjacent( p1, p2 );
1408     gp_Vec edgeNorm2 = -( edgeAdjacent ^ faceNorm );
1409     bool out2 = (edgeNorm2 * p2int) < -tol;
1410
1411     bool covexCorner = ( edgeNorm * edgeAdjacent * (rClosest==1. ? 1. : -1.)) < 0;
1412     return covexCorner ? (out || out2) : (out && out2);
1413   }
1414
1415   if ( element->GetType() == SMDSAbs_Edge ) // --------------------------------------------------
1416   {
1417     // point is out of edge if it is NOT ON any straight part of edge
1418     // (we consider quadratic edge as being composed of two straight parts)
1419     for ( i = 1; i < nbNodes; ++i )
1420     {
1421       gp_Vec edge( xyz[i-1], xyz[i] );
1422       gp_Vec n1p ( xyz[i-1], point  );
1423       double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
1424       if ( u <= 0. ) {
1425         if ( n1p.SquareMagnitude() < tol * tol )
1426           return false;
1427         continue;
1428       }
1429       if ( u >= 1. ) {
1430         if ( point.SquareDistance( xyz[i] ) < tol * tol )
1431           return false;
1432         continue;
1433       }
1434       gp_XYZ proj = ( 1. - u ) * xyz[i-1] + u * xyz[i]; // projection of the point on the edge
1435       double dist2 = point.SquareDistance( proj );
1436       if ( dist2 > tol * tol )
1437         continue;
1438       return false; // point is ON this part
1439     }
1440     return true;
1441   }
1442
1443   // Node or 0D element -------------------------------------------------------------------------
1444   {
1445     gp_Vec n2p ( xyz[0], point );
1446     return n2p.SquareMagnitude() > tol * tol;
1447   }
1448   return true;
1449 }
1450
1451 //=======================================================================
1452 namespace
1453 {
1454   // Position of a point relative to a segment
1455   //            .           .
1456   //            .  LEFT     .
1457   //            .           .
1458   //  VERTEX 1  o----ON----->  VERTEX 2
1459   //            .           .
1460   //            .  RIGHT    .
1461   //            .           .
1462   enum PositionName { POS_LEFT = 1, POS_VERTEX = 2, POS_RIGHT = 4, //POS_ON = 8,
1463                       POS_ALL = POS_LEFT | POS_RIGHT | POS_VERTEX };
1464   struct PointPos
1465   {
1466     PositionName _name;
1467     int          _index; // index of vertex or segment
1468
1469     PointPos( PositionName n, int i=-1 ): _name(n), _index(i) {}
1470     bool operator < (const PointPos& other ) const
1471     {
1472       if ( _name == other._name )
1473         return  ( _index < 0 || other._index < 0 ) ? false : _index < other._index;
1474       return _name < other._name;
1475     }
1476   };
1477
1478   //================================================================================
1479   /*!
1480    * \brief Return position of a point relative to a segment
1481    *  \param point2D      - the point to analyze position of
1482    *  \param segEnds      - end points of segments
1483    *  \param index0       - 0-based index of the first point of segment
1484    *  \param posToFindOut - flags of positions to detect
1485    *  \retval PointPos - point position
1486    */
1487   //================================================================================
1488
1489   PointPos getPointPosition( const gp_XY& point2D,
1490                              const gp_XY* segEnds,
1491                              const int    index0 = 0,
1492                              const int    posToFindOut = POS_ALL)
1493   {
1494     const gp_XY& p1 = segEnds[ index0   ];
1495     const gp_XY& p2 = segEnds[ index0+1 ];
1496     const gp_XY grad = p2 - p1;
1497
1498     if ( posToFindOut & POS_VERTEX )
1499     {
1500       // check if the point2D is at "vertex 1" zone
1501       gp_XY pp1[2] = { p1, gp_XY( p1.X() - grad.Y(),
1502                                   p1.Y() + grad.X() ) };
1503       if ( getPointPosition( point2D, pp1, 0, POS_LEFT|POS_RIGHT )._name == POS_LEFT )
1504         return PointPos( POS_VERTEX, index0 );
1505
1506       // check if the point2D is at "vertex 2" zone
1507       gp_XY pp2[2] = { p2, gp_XY( p2.X() - grad.Y(),
1508                                   p2.Y() + grad.X() ) };
1509       if ( getPointPosition( point2D, pp2, 0, POS_LEFT|POS_RIGHT )._name == POS_RIGHT )
1510         return PointPos( POS_VERTEX, index0 + 1);
1511     }
1512     double edgeEquation =
1513       ( point2D.X() - p1.X() ) * grad.Y() - ( point2D.Y() - p1.Y() ) * grad.X();
1514     return PointPos( edgeEquation < 0 ? POS_LEFT : POS_RIGHT, index0 );
1515   }
1516 }
1517
1518 //=======================================================================
1519 /*!
1520  * \brief Return minimal distance from a point to an element
1521  *
1522  * Currently we ignore non-planarity and 2nd order of face
1523  */
1524 //=======================================================================
1525
1526 double SMESH_MeshAlgos::GetDistance( const SMDS_MeshElement* elem,
1527                                      const gp_Pnt&           point,
1528                                      gp_XYZ*                 closestPnt )
1529 {
1530   switch ( elem->GetType() )
1531   {
1532   case SMDSAbs_Volume:
1533     return GetDistance( static_cast<const SMDS_MeshVolume*>( elem ), point, closestPnt );
1534   case SMDSAbs_Face:
1535     return GetDistance( static_cast<const SMDS_MeshFace*>( elem ), point, closestPnt );
1536   case SMDSAbs_Edge:
1537     return GetDistance( static_cast<const SMDS_MeshEdge*>( elem ), point, closestPnt );
1538   case SMDSAbs_Node:
1539     if ( closestPnt ) *closestPnt = SMESH_TNodeXYZ( elem );
1540     return point.Distance( SMESH_TNodeXYZ( elem ));
1541   default:;
1542   }
1543   return -1;
1544 }
1545
1546 //=======================================================================
1547 /*!
1548  * \brief Return minimal distance from a point to a face
1549  *
1550  * Currently we ignore non-planarity and 2nd order of face
1551  */
1552 //=======================================================================
1553
1554 double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face,
1555                                      const gp_Pnt&        point,
1556                                      gp_XYZ*              closestPnt )
1557 {
1558   const double badDistance = -1;
1559   if ( !face ) return badDistance;
1560
1561   // coordinates of nodes (medium nodes, if any, ignored)
1562   typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
1563   std::vector<gp_XYZ> xyz( TXyzIterator( face->nodesIterator()), TXyzIterator() );
1564   xyz.resize( face->NbCornerNodes()+1 );
1565
1566   // transformation to get xyz[0] lies on the origin, xyz[1] lies on the Z axis,
1567   // and xyz[2] lies in the XZ plane. This is to pass to 2D space on XZ plane.
1568   gp_Trsf trsf;
1569   gp_Vec OZ ( xyz[0], xyz[1] );
1570   gp_Vec OX ( xyz[0], xyz[2] );
1571   if ( OZ.Magnitude() < std::numeric_limits<double>::min() )
1572   {
1573     if ( xyz.size() < 4 ) return badDistance;
1574     OZ = gp_Vec ( xyz[0], xyz[2] );
1575     OX = gp_Vec ( xyz[0], xyz[3] );
1576   }
1577   gp_Ax3 tgtCS;
1578   try {
1579     tgtCS = gp_Ax3( xyz[0], OZ, OX );
1580   }
1581   catch ( Standard_Failure ) {
1582     return badDistance;
1583   }
1584   trsf.SetTransformation( tgtCS );
1585
1586   // move all the nodes to 2D
1587   std::vector<gp_XY> xy( xyz.size() );
1588   for ( size_t i = 0;i < xyz.size()-1; ++i )
1589   {
1590     gp_XYZ p3d = xyz[i];
1591     trsf.Transforms( p3d );
1592     xy[i].SetCoord( p3d.X(), p3d.Z() );
1593   }
1594   xyz.back() = xyz.front();
1595   xy.back() = xy.front();
1596
1597   // // move the point in 2D
1598   gp_XYZ tmpPnt = point.XYZ();
1599   trsf.Transforms( tmpPnt );
1600   gp_XY point2D( tmpPnt.X(), tmpPnt.Z() );
1601
1602   // loop on edges of the face to analyze point position ralative to the face
1603   std::set< PointPos > pntPosSet;
1604   for ( size_t i = 1; i < xy.size(); ++i )
1605   {
1606     PointPos pos = getPointPosition( point2D, &xy[0], i-1 );
1607     pntPosSet.insert( pos );
1608   }
1609
1610   // compute distance
1611
1612   double minDist2 = Precision::Infinite();
1613   for ( std::set< PointPos >::iterator posIt = pntPosSet.begin(); posIt != pntPosSet.end(); ++posIt)
1614   {
1615     PointPos pos = *posIt;
1616     if ( pos._name != pntPosSet.begin()->_name )
1617       break;
1618     switch ( pos._name )
1619     {
1620     case POS_LEFT: // point is most close to an edge
1621     {
1622       gp_Vec edge( xyz[ pos._index ], xyz[ pos._index+1 ]);
1623       gp_Vec n1p ( xyz[ pos._index ], point  );
1624       double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
1625       // projection of the point on the edge
1626       gp_XYZ proj = xyz[ pos._index ] + u * edge.XYZ();
1627       double dist2 = point.SquareDistance( proj );
1628       if ( dist2 < minDist2 )
1629       {
1630         if ( closestPnt ) *closestPnt = proj;
1631         minDist2 = dist2;
1632       }
1633       break;
1634     }
1635
1636     case POS_RIGHT: // point is inside the face
1637     {
1638       double distToFacePlane = Abs( tmpPnt.Y() );
1639       if ( closestPnt )
1640       {
1641         if ( distToFacePlane < std::numeric_limits<double>::min() ) {
1642           *closestPnt = point.XYZ();
1643         }
1644         else {
1645           tmpPnt.SetY( 0 );
1646           trsf.Inverted().Transforms( tmpPnt );
1647           *closestPnt = tmpPnt;
1648         }
1649       }
1650       return distToFacePlane;
1651     }
1652
1653     case POS_VERTEX: // point is most close to a node
1654     {
1655       double dist2 = point.SquareDistance( xyz[ pos._index ]);
1656       if ( dist2 < minDist2 )
1657       {
1658         if ( closestPnt ) *closestPnt = xyz[ pos._index ];
1659         minDist2 = dist2;
1660       }
1661       break;
1662     }
1663     default:;
1664       return badDistance;
1665     }
1666   }
1667   return Sqrt( minDist2 );
1668 }
1669
1670 //=======================================================================
1671 /*!
1672  * \brief Return minimal distance from a point to an edge
1673  */
1674 //=======================================================================
1675
1676 double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg,
1677                                      const gp_Pnt&        point,
1678                                      gp_XYZ*              closestPnt )
1679 {
1680   double dist = Precision::Infinite();
1681   if ( !seg ) return dist;
1682
1683   int i = 0, nbNodes = seg->NbNodes();
1684
1685   std::vector< SMESH_TNodeXYZ > xyz( nbNodes );
1686   for ( SMDS_NodeIteratorPtr nodeIt = seg->interlacedNodesIterator(); nodeIt->more(); i++ )
1687     xyz[ i ].Set( nodeIt->next() );
1688
1689   for ( i = 1; i < nbNodes; ++i )
1690   {
1691     gp_Vec edge( xyz[i-1], xyz[i] );
1692     gp_Vec n1p ( xyz[i-1], point  );
1693     double d, u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
1694     if ( u <= 0. ) {
1695       if (( d = n1p.SquareMagnitude() ) < dist ) {
1696         dist = d;
1697         if ( closestPnt ) *closestPnt = xyz[i-1];
1698       }
1699     }
1700     else if ( u >= 1. ) {
1701       if (( d = point.SquareDistance( xyz[i] )) < dist ) {
1702         dist = d;
1703         if ( closestPnt ) *closestPnt = xyz[i];
1704       }
1705     }
1706     else {
1707       gp_XYZ proj = xyz[i-1] + u * edge.XYZ(); // projection of the point on the edge
1708       if (( d = point.SquareDistance( proj )) < dist ) {
1709         dist = d;
1710         if ( closestPnt ) *closestPnt = proj;
1711       }
1712     }
1713   }
1714   return Sqrt( dist );
1715 }
1716
1717 //=======================================================================
1718 /*!
1719  * \brief Return minimal distance from a point to a volume
1720  *
1721  * Currently we ignore non-planarity and 2nd order
1722  */
1723 //=======================================================================
1724
1725 double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume,
1726                                      const gp_Pnt&          point,
1727                                      gp_XYZ*                closestPnt )
1728 {
1729   SMDS_VolumeTool vTool( volume );
1730   vTool.SetExternalNormal();
1731   const int iQ = volume->IsQuadratic() ? 2 : 1;
1732
1733   double n[3], bc[3];
1734   double minDist = 1e100, dist;
1735   gp_XYZ closeP = point.XYZ();
1736   bool isOut = false;
1737   for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
1738   {
1739     // skip a facet with normal not "looking at" the point
1740     if ( !vTool.GetFaceNormal( iF, n[0], n[1], n[2] ) ||
1741          !vTool.GetFaceBaryCenter( iF, bc[0], bc[1], bc[2] ))
1742       continue;
1743     gp_XYZ bcp = point.XYZ() - gp_XYZ( bc[0], bc[1], bc[2] );
1744     if ( gp_XYZ( n[0], n[1], n[2] ) * bcp < 1e-6 )
1745       continue;
1746
1747     // find distance to a facet
1748     const SMDS_MeshNode** nodes = vTool.GetFaceNodes( iF );
1749     switch ( vTool.NbFaceNodes( iF ) / iQ ) {
1750     case 3:
1751     {
1752       SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ] );
1753       dist = GetDistance( &tmpFace, point, closestPnt );
1754       break;
1755     }
1756     case 4:
1757     {
1758       SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ], nodes[ 3*iQ ]);
1759       dist = GetDistance( &tmpFace, point, closestPnt );
1760       break;
1761     }
1762     default:
1763       std::vector<const SMDS_MeshNode *> nvec( nodes, nodes + vTool.NbFaceNodes( iF ));
1764       SMDS_PolygonalFaceOfNodes tmpFace( nvec );
1765       dist = GetDistance( &tmpFace, point, closestPnt );
1766     }
1767     if ( dist < minDist )
1768     {
1769       minDist = dist;
1770       isOut = true;
1771       if ( closestPnt ) closeP = *closestPnt;
1772     }
1773   }
1774   if ( isOut )
1775   {
1776     if ( closestPnt ) *closestPnt = closeP;
1777     return minDist;
1778   }
1779
1780   return 0; // point is inside the volume
1781 }
1782
1783 //================================================================================
1784 /*!
1785  * \brief Returns barycentric coordinates of a point within a triangle.
1786  *        A not returned bc2 = 1. - bc0 - bc1.
1787  *        The point lies within the triangle if ( bc0 >= 0 && bc1 >= 0 && bc0+bc1 <= 1 )
1788  */
1789 //================================================================================
1790
1791 void SMESH_MeshAlgos::GetBarycentricCoords( const gp_XY& p,
1792                                             const gp_XY& t0,
1793                                             const gp_XY& t1,
1794                                             const gp_XY& t2,
1795                                             double &     bc0,
1796                                             double &     bc1)
1797 {
1798   const double // matrix 2x2
1799     T11 = t0.X()-t2.X(), T12 = t1.X()-t2.X(),
1800     T21 = t0.Y()-t2.Y(), T22 = t1.Y()-t2.Y();
1801   const double Tdet = T11*T22 - T12*T21; // matrix determinant
1802   if ( Abs( Tdet ) < std::numeric_limits<double>::min() )
1803   {
1804     bc0 = bc1 = 2.;
1805     return;
1806   }
1807   // matrix inverse
1808   const double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
1809   // vector
1810   const double r11 = p.X()-t2.X(), r12 = p.Y()-t2.Y();
1811   // barycentric coordinates: multiply matrix by vector
1812   bc0 = (t11 * r11 + t12 * r12)/Tdet;
1813   bc1 = (t21 * r11 + t22 * r12)/Tdet;
1814 }
1815
1816 //=======================================================================
1817 //function : FindFaceInSet
1818 //purpose  : Return a face having linked nodes n1 and n2 and which is
1819 //           - not in avoidSet,
1820 //           - in elemSet provided that !elemSet.empty()
1821 //           i1 and i2 optionally returns indices of n1 and n2
1822 //=======================================================================
1823
1824 const SMDS_MeshElement*
1825 SMESH_MeshAlgos::FindFaceInSet(const SMDS_MeshNode*    n1,
1826                                const SMDS_MeshNode*    n2,
1827                                const TIDSortedElemSet& elemSet,
1828                                const TIDSortedElemSet& avoidSet,
1829                                int*                    n1ind,
1830                                int*                    n2ind)
1831
1832 {
1833   int i1 = 0, i2 = 0;
1834   const SMDS_MeshElement* face = 0;
1835
1836   SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
1837   while ( invElemIt->more() && !face ) // loop on inverse faces of n1
1838   {
1839     const SMDS_MeshElement* elem = invElemIt->next();
1840     if (avoidSet.count( elem ))
1841       continue;
1842     if ( !elemSet.empty() && !elemSet.count( elem ))
1843       continue;
1844     // index of n1
1845     i1 = elem->GetNodeIndex( n1 );
1846     // find a n2 linked to n1
1847     int nbN = elem->IsQuadratic() ? elem->NbNodes()/2 : elem->NbNodes();
1848     for ( int di = -1; di < 2 && !face; di += 2 )
1849     {
1850       i2 = (i1+di+nbN) % nbN;
1851       if ( elem->GetNode( i2 ) == n2 )
1852         face = elem;
1853     }
1854     if ( !face && elem->IsQuadratic())
1855     {
1856       // analysis for quadratic elements using all nodes
1857       SMDS_NodeIteratorPtr anIter = elem->interlacedNodesIterator();
1858       const SMDS_MeshNode* prevN = static_cast<const SMDS_MeshNode*>( anIter->next() );
1859       for ( i1 = -1, i2 = 0; anIter->more() && !face; i1++, i2++ )
1860       {
1861         const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( anIter->next() );
1862         if ( n1 == prevN && n2 == n )
1863         {
1864           face = elem;
1865         }
1866         else if ( n2 == prevN && n1 == n )
1867         {
1868           face = elem; std::swap( i1, i2 );
1869         }
1870         prevN = n;
1871       }
1872     }
1873   }
1874   if ( n1ind ) *n1ind = i1;
1875   if ( n2ind ) *n2ind = i2;
1876   return face;
1877 }
1878
1879 //================================================================================
1880 /*!
1881  * Return sharp edges of faces and non-manifold ones. Optionally adds existing edges.
1882  */
1883 //================================================================================
1884
1885 std::vector< SMESH_MeshAlgos::Edge >
1886 SMESH_MeshAlgos::FindSharpEdges( SMDS_Mesh* theMesh,
1887                                  double     theAngle,
1888                                  bool       theAddExisting )
1889 {
1890   std::vector< Edge > resultEdges;
1891   if ( !theMesh ) return resultEdges;
1892
1893   typedef std::pair< bool, const SMDS_MeshNode* >                            TIsSharpAndMedium;
1894   typedef NCollection_DataMap< SMESH_TLink, TIsSharpAndMedium, SMESH_TLink > TLinkSharpMap;
1895
1896   TLinkSharpMap linkIsSharp( theMesh->NbFaces() );
1897   TIsSharpAndMedium sharpMedium( true, 0 );
1898   bool                 & isSharp = sharpMedium.first;
1899   const SMDS_MeshNode* & nMedium = sharpMedium.second;
1900
1901   if ( theAddExisting )
1902   {
1903     for ( SMDS_EdgeIteratorPtr edgeIt = theMesh->edgesIterator(); edgeIt->more(); )
1904     {
1905       const SMDS_MeshElement* edge = edgeIt->next();
1906       nMedium = ( edge->IsQuadratic() ) ? edge->GetNode(2) : 0;
1907       linkIsSharp.Bind( SMESH_TLink( edge->GetNode(0), edge->GetNode(1)), sharpMedium );
1908     }
1909   }
1910
1911   // check angles between face normals
1912
1913   const double angleCos = Cos( theAngle * M_PI / 180. ), angleCos2 = angleCos * angleCos;
1914   gp_XYZ norm1, norm2;
1915   std::vector< const SMDS_MeshNode* > faceNodes, linkNodes(2);
1916   std::vector<const SMDS_MeshElement *> linkFaces;
1917
1918   int nbSharp = linkIsSharp.Extent();
1919   for ( SMDS_FaceIteratorPtr faceIt = theMesh->facesIterator(); faceIt->more(); )
1920   {
1921     const SMDS_MeshElement* face = faceIt->next();
1922     size_t             nbCorners = face->NbCornerNodes();
1923
1924     faceNodes.assign( face->begin_nodes(), face->end_nodes() );
1925     if ( faceNodes.size() == nbCorners )
1926       faceNodes.resize( nbCorners * 2, 0 );
1927
1928     const SMDS_MeshNode* nPrev = faceNodes[ nbCorners-1 ];
1929     for ( size_t i = 0; i < nbCorners; ++i )
1930     {
1931       SMESH_TLink link( nPrev, faceNodes[i] );
1932       if ( !linkIsSharp.IsBound( link ))
1933       {
1934         linkNodes[0] = link.node1();
1935         linkNodes[1] = link.node2();
1936         linkFaces.clear();
1937         theMesh->GetElementsByNodes( linkNodes, linkFaces, SMDSAbs_Face );
1938
1939         isSharp = false;
1940         if ( linkFaces.size() > 2 )
1941         {
1942           isSharp = true;
1943         }
1944         else if ( linkFaces.size() == 2 &&
1945                   FaceNormal( linkFaces[0], norm1, /*normalize=*/false ) &&
1946                   FaceNormal( linkFaces[1], norm2, /*normalize=*/false ))
1947         {
1948           double dot = norm1 * norm2; // == cos * |norm1| * |norm2|
1949           if (( dot < 0 ) == ( angleCos < 0 ))
1950           {
1951             double cos2 = dot * dot / norm1.SquareModulus() / norm2.SquareModulus();
1952             isSharp = ( angleCos < 0 ) ? ( cos2 > angleCos2 ) : ( cos2 < angleCos2 );
1953           }
1954           else
1955           {
1956             isSharp = ( angleCos > 0 );
1957           }
1958         }
1959         nMedium = faceNodes[( i-1+nbCorners ) % nbCorners + nbCorners ];
1960
1961         linkIsSharp.Bind( link, sharpMedium );
1962         nbSharp += isSharp;
1963       }
1964
1965       nPrev = faceNodes[i];
1966     }
1967   }
1968
1969   resultEdges.resize( nbSharp );
1970   TLinkSharpMap::Iterator linkIsSharpIter( linkIsSharp );
1971   for ( int i = 0; linkIsSharpIter.More() && i < nbSharp; linkIsSharpIter.Next() )
1972   {
1973     const SMESH_TLink&                link = linkIsSharpIter.Key();
1974     const TIsSharpAndMedium& isSharpMedium = linkIsSharpIter.Value();
1975     if ( isSharpMedium.first )
1976     {
1977       Edge & edge  = resultEdges[ i++ ];
1978       edge._node1  = link.node1();
1979       edge._node2  = link.node2();
1980       edge._medium = isSharpMedium.second;
1981     }
1982   }
1983
1984   return resultEdges;
1985 }
1986
1987 //================================================================================
1988 /*!
1989  * Distribute all faces of the mesh between groups using given edges as group boundaries
1990  */
1991 //================================================================================
1992
1993 std::vector< std::vector< const SMDS_MeshElement* > >
1994 SMESH_MeshAlgos::SeparateFacesByEdges( SMDS_Mesh* theMesh, const std::vector< Edge >& theEdges )
1995 {
1996   std::vector< std::vector< const SMDS_MeshElement* > > groups;
1997   if ( !theMesh ) return groups;
1998
1999   // build map of face edges (SMESH_TLink) and their faces
2000
2001   typedef std::vector< const SMDS_MeshElement* >                    TFaceVec;
2002   typedef NCollection_DataMap< SMESH_TLink, TFaceVec, SMESH_TLink > TFacesByLinks;
2003   TFacesByLinks facesByLink( theMesh->NbFaces() );
2004
2005   std::vector< const SMDS_MeshNode* > faceNodes;
2006   for ( SMDS_FaceIteratorPtr faceIt = theMesh->facesIterator(); faceIt->more(); )
2007   {
2008     const SMDS_MeshElement* face = faceIt->next();
2009     size_t             nbCorners = face->NbCornerNodes();
2010
2011     faceNodes.assign( face->begin_nodes(), face->end_nodes() );
2012     faceNodes.resize( nbCorners + 1 );
2013     faceNodes[ nbCorners ] = faceNodes[0];
2014
2015     face->setIsMarked( false );
2016
2017     for ( size_t i = 0; i < nbCorners; ++i )
2018     {
2019       SMESH_TLink link( faceNodes[i], faceNodes[i+1] );
2020       TFaceVec* linkFaces = facesByLink.ChangeSeek( link );
2021       if ( !linkFaces )
2022       {
2023         linkFaces = facesByLink.Bound( link, TFaceVec() );
2024         linkFaces->reserve(2);
2025       }
2026       linkFaces->push_back( face );
2027     }
2028   }
2029
2030   // remove the given edges from facesByLink map
2031
2032   for ( size_t i = 0; i < theEdges.size(); ++i )
2033   {
2034     SMESH_TLink link( theEdges[i]._node1, theEdges[i]._node2 );
2035     facesByLink.UnBind( link );
2036   }
2037
2038   // faces connected via links of facesByLink map form a group
2039
2040   while ( !facesByLink.IsEmpty() )
2041   {
2042     groups.push_back( TFaceVec() );
2043     TFaceVec & group = groups.back();
2044
2045     group.push_back( TFacesByLinks::Iterator( facesByLink ).Value()[0] );
2046     group.back()->setIsMarked( true );
2047
2048     for ( size_t iF = 0; iF < group.size(); ++iF )
2049     {
2050       const SMDS_MeshElement* face = group[iF];
2051       size_t             nbCorners = face->NbCornerNodes();
2052       faceNodes.assign( face->begin_nodes(), face->end_nodes() );
2053       faceNodes.resize( nbCorners + 1 );
2054       faceNodes[ nbCorners ] = faceNodes[0];
2055
2056       for ( size_t iN = 0; iN < nbCorners; ++iN )
2057       {
2058         SMESH_TLink link( faceNodes[iN], faceNodes[iN+1] );
2059         if ( const TFaceVec* faces = facesByLink.Seek( link ))
2060         {
2061           const TFaceVec& faceNeighbors = *faces;
2062           for ( size_t i = 0; i < faceNeighbors.size(); ++i )
2063             if ( !faceNeighbors[i]->isMarked() )
2064             {
2065               group.push_back( faceNeighbors[i] );
2066               faceNeighbors[i]->setIsMarked( true );
2067             }
2068           facesByLink.UnBind( link );
2069         }
2070       }
2071     }
2072   }
2073
2074   // find faces that are alone in its group; they were not in facesByLink
2075
2076   int nbInGroups = 0;
2077   for ( size_t i = 0; i < groups.size(); ++i )
2078     nbInGroups += groups[i].size();
2079   if ( nbInGroups < theMesh->NbFaces() )
2080   {
2081     for ( SMDS_FaceIteratorPtr faceIt = theMesh->facesIterator(); faceIt->more(); )
2082     {
2083       const SMDS_MeshElement* face = faceIt->next();
2084       if ( !face->isMarked() )
2085       {
2086         groups.push_back( TFaceVec() );
2087         groups.back().push_back( face );
2088       }
2089     }
2090   }
2091
2092   return groups;
2093 }
2094
2095 //================================================================================
2096 /*!
2097  * \brief Calculate normal of a mesh face
2098  */
2099 //================================================================================
2100
2101 bool SMESH_MeshAlgos::FaceNormal(const SMDS_MeshElement* F, gp_XYZ& normal, bool normalized)
2102 {
2103   if ( !F || F->GetType() != SMDSAbs_Face )
2104     return false;
2105
2106   normal.SetCoord(0,0,0);
2107   int nbNodes = F->NbCornerNodes();
2108   for ( int i = 0; i < nbNodes-2; ++i )
2109   {
2110     gp_XYZ p[3];
2111     for ( int n = 0; n < 3; ++n )
2112     {
2113       const SMDS_MeshNode* node = F->GetNode( i + n );
2114       p[n].SetCoord( node->X(), node->Y(), node->Z() );
2115     }
2116     normal += ( p[2] - p[1] ) ^ ( p[0] - p[1] );
2117   }
2118   double size2 = normal.SquareModulus();
2119   bool ok = ( size2 > std::numeric_limits<double>::min() * std::numeric_limits<double>::min());
2120   if ( normalized && ok )
2121     normal /= sqrt( size2 );
2122
2123   return ok;
2124 }
2125
2126 //=======================================================================
2127 //function : GetCommonNodes
2128 //purpose  : Return nodes common to two elements
2129 //=======================================================================
2130
2131 std::vector< const SMDS_MeshNode*> SMESH_MeshAlgos::GetCommonNodes(const SMDS_MeshElement* e1,
2132                                                                    const SMDS_MeshElement* e2)
2133 {
2134   std::vector< const SMDS_MeshNode*> common;
2135   for ( int i = 0 ; i < e1->NbNodes(); ++i )
2136     if ( e2->GetNodeIndex( e1->GetNode( i )) >= 0 )
2137       common.push_back( e1->GetNode( i ));
2138   return common;
2139 }
2140 //================================================================================
2141 /*!
2142  * \brief Return true if node1 encounters first in the face and node2, after
2143  */
2144 //================================================================================
2145
2146 bool SMESH_MeshAlgos::IsRightOrder( const SMDS_MeshElement* face,
2147                                     const SMDS_MeshNode*    node0,
2148                                     const SMDS_MeshNode*    node1 )
2149 {
2150   int i0 = face->GetNodeIndex( node0 );
2151   int i1 = face->GetNodeIndex( node1 );
2152   if ( face->IsQuadratic() )
2153   {
2154     if ( face->IsMediumNode( node0 ))
2155     {
2156       i0 -= ( face->NbNodes()/2 - 1 );
2157       i1 *= 2;
2158     }
2159     else
2160     {
2161       i1 -= ( face->NbNodes()/2 - 1 );
2162       i0 *= 2;
2163     }
2164   }
2165   int diff = i1 - i0;
2166   return ( diff == 1 ) || ( diff == -face->NbNodes()+1 );
2167 }
2168
2169 //=======================================================================
2170 /*!
2171  * \brief Return SMESH_NodeSearcher
2172  */
2173 //=======================================================================
2174
2175 SMESH_NodeSearcher* SMESH_MeshAlgos::GetNodeSearcher(SMDS_Mesh& mesh)
2176 {
2177   return new SMESH_NodeSearcherImpl( &mesh );
2178 }
2179
2180 //=======================================================================
2181 /*!
2182  * \brief Return SMESH_NodeSearcher
2183  */
2184 //=======================================================================
2185
2186 SMESH_NodeSearcher* SMESH_MeshAlgos::GetNodeSearcher(SMDS_ElemIteratorPtr elemIt)
2187 {
2188   return new SMESH_NodeSearcherImpl( 0, elemIt );
2189 }
2190
2191 //=======================================================================
2192 /*!
2193  * \brief Return SMESH_ElementSearcher
2194  */
2195 //=======================================================================
2196
2197 SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh& mesh,
2198                                                            double     tolerance)
2199 {
2200   return new SMESH_ElementSearcherImpl( mesh, tolerance );
2201 }
2202
2203 //=======================================================================
2204 /*!
2205  * \brief Return SMESH_ElementSearcher acting on a sub-set of elements
2206  */
2207 //=======================================================================
2208
2209 SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh&           mesh,
2210                                                            SMDS_ElemIteratorPtr elemIt,
2211                                                            double               tolerance)
2212 {
2213   return new SMESH_ElementSearcherImpl( mesh, tolerance, elemIt );
2214 }