Salome HOME
Windows compatibility.
[modules/smesh.git] / src / SMESHUtils / SMESH_MeshAlgos.cxx
index dbadd82505154cbcbbae5e0fca2c686f205a9697..1eae5f7818e23d205812d4330a7e8437c34349cb 100644 (file)
@@ -441,7 +441,8 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
 {
   SMDS_Mesh*                   _mesh;
   SMDS_ElemIteratorPtr         _meshPartIt;
-  ElementBndBoxTree*           _ebbTree;
+  ElementBndBoxTree*           _ebbTree      [SMDSAbs_NbElementTypes];
+  int                          _ebbTreeHeight[SMDSAbs_NbElementTypes];
   SMESH_NodeSearcherImpl*      _nodeSearcher;
   SMDSAbs_ElementType          _elementType;
   double                       _tolerance;
@@ -451,10 +452,21 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   SMESH_ElementSearcherImpl( SMDS_Mesh&           mesh,
                              double               tol=-1,
                              SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr())
-    : _mesh(&mesh),_meshPartIt(elemIt),_ebbTree(0),_nodeSearcher(0),_tolerance(tol),_outerFacesFound(false) {}
+    : _mesh(&mesh),_meshPartIt(elemIt),_nodeSearcher(0),_tolerance(tol),_outerFacesFound(false)
+  {
+    for ( int i = 0; i < SMDSAbs_NbElementTypes; ++i )
+    {
+      _ebbTree[i] = NULL;
+      _ebbTreeHeight[i] = -1;
+    }
+    _elementType = SMDSAbs_All;
+  }
   virtual ~SMESH_ElementSearcherImpl()
   {
-    if ( _ebbTree )      delete _ebbTree;      _ebbTree      = 0;
+    for ( int i = 0; i < SMDSAbs_NbElementTypes; ++i )
+    {
+      delete _ebbTree[i]; _ebbTree[i] = NULL;
+    }
     if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0;
   }
   virtual int FindElementsByPoint(const gp_Pnt&                      point,
@@ -479,6 +491,12 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   {
     return _outerFaces.empty() || _outerFaces.count(face);
   }
+  int getTreeHeight()
+  {
+    if ( _ebbTreeHeight[ _elementType ] < 0 )
+      _ebbTreeHeight[ _elementType ] = _ebbTree[ _elementType ]->getHeight();
+    return _ebbTreeHeight[ _elementType ];
+  }
 
   struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
   {
@@ -521,9 +539,9 @@ double SMESH_ElementSearcherImpl::getTolerance()
       double boxSize = _nodeSearcher->getTree()->maxSize();
       _tolerance = 1e-8 * boxSize/* / meshInfo.NbNodes()*/;
     }
-    else if ( _ebbTree && meshInfo.NbElements() > 0 )
+    else if ( _ebbTree[_elementType] && meshInfo.NbElements() > 0 )
     {
-      double boxSize = _ebbTree->maxSize();
+      double boxSize = _ebbTree[_elementType]->maxSize();
       _tolerance = 1e-8 * boxSize/* / meshInfo.NbElements()*/;
     }
     if ( _tolerance == 0 )
@@ -544,10 +562,9 @@ double SMESH_ElementSearcherImpl::getTolerance()
       }
       else
       {
-        SMDS_ElemIteratorPtr elemIt =
-            _mesh->elementsIterator( SMDSAbs_ElementType( complexType ));
+        SMDS_ElemIteratorPtr  elemIt = _mesh->elementsIterator( SMDSAbs_ElementType( complexType ));
         const SMDS_MeshElement* elem = elemIt->next();
-        SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+        SMDS_ElemIteratorPtr  nodeIt = elem->nodesIterator();
         SMESH_TNodeXYZ n1( nodeIt->next() );
         elemSize = 0;
         while ( nodeIt->more() )
@@ -672,7 +689,7 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF
           outerFace2 = angle2Face.begin()->second;
       }
     }
-    // store the found outer face and add its links to continue seaching from
+    // store the found outer face and add its links to continue searching from
     if ( outerFace2 )
     {
       _outerFaces.insert( outerFace2 );
@@ -716,6 +733,7 @@ FindElementsByPoint(const gp_Pnt&                      point,
                     vector< const SMDS_MeshElement* >& foundElements)
 {
   foundElements.clear();
+  _elementType = type;
 
   double tolerance = getTolerance();
 
@@ -749,13 +767,12 @@ FindElementsByPoint(const gp_Pnt&                      point,
   // =================================================================================
   else // elements more complex than 0D
   {
-    if ( !_ebbTree || _elementType != type )
+    if ( !_ebbTree[type] )
     {
-      if ( _ebbTree ) delete _ebbTree;
-      _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt, tolerance );
+      _ebbTree[_elementType] = new ElementBndBoxTree( *_mesh, type, _meshPartIt, tolerance );
     }
     TIDSortedElemSet suspectElems;
-    _ebbTree->getElementsNearPoint( point, suspectElems );
+    _ebbTree[ type ]->getElementsNearPoint( point, suspectElems );
     TIDSortedElemSet::iterator elem = suspectElems.begin();
     for ( ; elem != suspectElems.end(); ++elem )
       if ( !SMESH_MeshAlgos::IsOut( *elem, point, tolerance ))
@@ -777,29 +794,29 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
                                           SMDSAbs_ElementType type )
 {
   const SMDS_MeshElement* closestElem = 0;
+  _elementType = type;
 
   if ( type == SMDSAbs_Face || type == SMDSAbs_Volume )
   {
-    if ( !_ebbTree || _elementType != type )
-    {
-      if ( _ebbTree ) delete _ebbTree;
-      _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
-    }
+    ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
+    if ( !ebbTree )
+      ebbTree = new ElementBndBoxTree( *_mesh, type, _meshPartIt );
+
     TIDSortedElemSet suspectElems;
-    _ebbTree->getElementsNearPoint( point, suspectElems );
+    ebbTree->getElementsNearPoint( point, suspectElems );
 
-    if ( suspectElems.empty() && _ebbTree->maxSize() > 0 )
+    if ( suspectElems.empty() && ebbTree->maxSize() > 0 )
     {
-      gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox()->CornerMin() +
-                                 _ebbTree->getBox()->CornerMax() );
+      gp_Pnt boxCenter = 0.5 * ( ebbTree->getBox()->CornerMin() +
+                                 ebbTree->getBox()->CornerMax() );
       double radius = -1;
-      if ( _ebbTree->getBox()->IsOut( point.XYZ() ))
-        radius = point.Distance( boxCenter ) - 0.5 * _ebbTree->maxSize();
+      if ( ebbTree->getBox()->IsOut( point.XYZ() ))
+        radius = point.Distance( boxCenter ) - 0.5 * ebbTree->maxSize();
       if ( radius < 0 )
-        radius = _ebbTree->maxSize() / pow( 2., _ebbTree->getHeight()) / 2;
+        radius = ebbTree->maxSize() / pow( 2., getTreeHeight()) / 2;
       while ( suspectElems.empty() )
       {
-        _ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems );
+        ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems );
         radius *= 1.1;
       }
     }
@@ -862,12 +879,14 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
 
 TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
 {
+  _elementType = SMDSAbs_Face;
+
   double tolerance = getTolerance();
-  if ( !_ebbTree || _elementType != SMDSAbs_Face )
-  {
-    if ( _ebbTree ) delete _ebbTree;
-    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face, _meshPartIt );
-  }
+
+  ElementBndBoxTree*& ebbTree = _ebbTree[ SMDSAbs_Face ];
+  if ( !ebbTree )
+    ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
+
   // Algo: analyse transition of a line starting at the point through mesh boundary;
   // try three lines parallel to axis of the coordinate system and perform rough
   // analysis. If solution is not clear perform thorough analysis.
@@ -883,7 +902,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
     gp_Lin line    ( lineAxis );
 
     TIDSortedElemSet suspectFaces; // faces possibly intersecting the line
-    _ebbTree->getElementsNearLine( lineAxis, suspectFaces );
+    ebbTree->getElementsNearLine( lineAxis, suspectFaces );
 
     // Intersect faces with the line
 
@@ -1091,13 +1110,13 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1&
                                                      SMDSAbs_ElementType                type,
                                                      vector< const SMDS_MeshElement* >& foundElems)
 {
-  if ( !_ebbTree || _elementType != type )
-  {
-    if ( _ebbTree ) delete _ebbTree;
-    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
-  }
+  _elementType = type;
+  ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
+  if ( !ebbTree )
+    ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
+
   TIDSortedElemSet suspectFaces; // elements possibly intersecting the line
-  _ebbTree->getElementsNearLine( line, suspectFaces );
+  ebbTree->getElementsNearLine( line, suspectFaces );
   foundElems.assign( suspectFaces.begin(), suspectFaces.end());
 }
 
@@ -1112,13 +1131,13 @@ void SMESH_ElementSearcherImpl::GetElementsInSphere( const gp_XYZ&
                                                      SMDSAbs_ElementType                type,
                                                      vector< const SMDS_MeshElement* >& foundElems)
 {
-  if ( !_ebbTree || _elementType != type )
-  {
-    if ( _ebbTree ) delete _ebbTree;
-    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
-  }
+  _elementType = type;
+  ElementBndBoxTree*& ebbTree = _ebbTree[ type ];
+  if ( !ebbTree )
+    ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt );
+
   TIDSortedElemSet suspectFaces; // elements possibly intersecting the line
-  _ebbTree->getElementsInSphere( center, radius, suspectFaces );
+  ebbTree->getElementsInSphere( center, radius, suspectFaces );
   foundElems.assign( suspectFaces.begin(), suspectFaces.end() );
 }
 
@@ -1156,8 +1175,8 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
       gp_Vec edge2( xyz[i+1], xyz[(i+2)%nbNodes] );
       faceNorm += edge1 ^ edge2;
     }
-    double normSize = faceNorm.Magnitude();
-    if ( normSize <= tol )
+    double fNormSize = faceNorm.Magnitude();
+    if ( fNormSize <= tol )
     {
       // degenerated face: point is out if it is out of all face edges
       for ( i = 0; i < nbNodes; ++i )
@@ -1168,7 +1187,7 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
       }
       return true;
     }
-    faceNorm /= normSize;
+    faceNorm /= fNormSize;
 
     // check if the point lays on face plane
     gp_Vec n2p( xyz[0], point );
@@ -1197,9 +1216,10 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     // to find intersections of the ray with the boundary.
     gp_Vec ray = n2p;
     gp_Vec plnNorm = ray ^ faceNorm;
-    normSize = plnNorm.Magnitude();
-    if ( normSize <= tol ) return false; // point coincides with the first node
-    plnNorm /= normSize;
+    double n2pSize = plnNorm.Magnitude();
+    if ( n2pSize <= tol ) return false; // point coincides with the first node
+    if ( n2pSize * n2pSize > fNormSize * 100 ) return true; // point is very far
+    plnNorm /= n2pSize;
     // for each node of the face, compute its signed distance to the cutting plane
     vector<double> dist( nbNodes + 1);
     for ( i = 0; i < nbNodes; ++i )
@@ -1245,7 +1265,7 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     if ( rClosest > 0. && rClosest < 1. ) // not node intersection
       return out;
 
-    // ray pass through a face node; analyze transition through an adjacent edge
+    // the ray passes through a face node; analyze transition through an adjacent edge
     gp_Pnt p1 = xyz[ (rClosest == 0.) ? ((iClosest+nbNodes-1) % nbNodes) : (iClosest+1) ];
     gp_Pnt p2 = xyz[ (rClosest == 0.) ? iClosest : ((iClosest+2) % nbNodes) ];
     gp_Vec edgeAdjacent( p1, p2 );