Salome HOME
bos #20256: [CEA 18523] Porting SMESH to int 64 bits
[modules/smesh.git] / src / SMESHUtils / SMESH_MeshAlgos.cxx
index 52d7d5415d51943f7b26b5ec4db8e646590a6593..9f0c88e20e9c849f9e1efe4dc7c1c24ea94853f7 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -238,6 +238,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     void getElementsInBox    ( const Bnd_B3d& box,  TElemSeq& foundElems );
     void getElementsInSphere ( const gp_XYZ& center, const double radius, TElemSeq& foundElems );
     ElementBndBoxTree* getLeafAtPoint( const gp_XYZ& point );
+    int  getNbElements();
 
   protected:
     ElementBndBoxTree() {}
@@ -280,7 +281,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
                                        double               tolerance)
     :SMESH_Octree( new LimitAndPool() )
   {
-    int nbElems = mesh.GetMeshInfo().NbElements( elemType );
+    smIdType nbElems = mesh.GetMeshInfo().NbElements( elemType );
     _elements.reserve( nbElems );
 
     TElementBoxPool& elBoPool = getLimitAndPool()->_elBoPool;
@@ -316,7 +317,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 
   //================================================================================
   /*!
-   * \brief Redistrubute element boxes among children
+   * \brief Redistribute element boxes among children
    */
   //================================================================================
 
@@ -466,6 +467,27 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     return 0;
   }
 
+  //================================================================================
+  /*!
+   * \brief Return number of elements
+   */
+  //================================================================================
+
+  int ElementBndBoxTree::getNbElements()
+  {
+    int nb = 0;
+    if ( isLeaf() )
+    {
+      nb = _elements.size();
+    }
+    else
+    {
+      for (int i = 0; i < 8; i++)
+        nb += ((ElementBndBoxTree*) myChildren[i])->getNbElements();
+    }
+    return nb;
+  }
+
   //================================================================================
   /*!
    * \brief Construct the element box
@@ -524,7 +546,8 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
     {
       delete _ebbTree[i]; _ebbTree[i] = NULL;
     }
-    if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0;
+    if ( _nodeSearcher ) delete _nodeSearcher;
+    _nodeSearcher = 0;
   }
   virtual int FindElementsByPoint(const gp_Pnt&                           point,
                                   SMDSAbs_ElementType                     type,
@@ -665,6 +688,9 @@ bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin&           lin
     GC_MakeSegment edge( SMESH_TNodeXYZ( face->GetNode( i )),
                          SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) ));
     anExtCC.Init( lineCurve, edge.Value() );
+    if ( !anExtCC.Extrema().IsDone() ||
+         anExtCC.Extrema().IsParallel() )
+      continue;
     if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol)
     {
       Standard_Real pl, pe;
@@ -848,7 +874,7 @@ FindElementsByPoint(const gp_Pnt&                           point,
 /*!
  * \brief Find an element of given type most close to the given point
  *
- * WARNING: Only face search is implemeneted so far
+ * WARNING: Only edge, face and volume search is implemented so far
  */
 //=======================================================================
 
@@ -1076,7 +1102,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
             const SMDS_MeshElement* prevFace = u_int1->second._face;
             while ( ok && u_int2->second._coincides )
             {
-              if ( SMESH_MeshAlgos::GetCommonNodes(prevFace , u_int2->second._face).empty() )
+              if ( SMESH_MeshAlgos::NbCommonNodes(prevFace , u_int2->second._face) == 0 )
                 ok = false;
               else
               {
@@ -1254,18 +1280,43 @@ gp_XYZ SMESH_ElementSearcherImpl::Project(const gp_Pnt&            point,
   gp_XYZ p = point.XYZ();
   ElementBndBoxTree* ebbLeaf = ebbTree->getLeafAtPoint( p );
   const Bnd_B3d* box = ebbLeaf ? ebbLeaf->getBox() : ebbTree->getBox();
-  double radius = ( box->CornerMax() - box->CornerMin() ).Modulus();
+  gp_XYZ pMin = box->CornerMin(), pMax = box->CornerMax();
+  double radius = Precision::Infinite();
+  if ( ebbLeaf || !box->IsOut( p ))
+  {
+    for ( int i = 1; i <= 3; ++i )
+    {
+      double d = 0.5 * ( pMax.Coord(i) - pMin.Coord(i) );
+      if ( d > Precision::Confusion() )
+        radius = Min( d, radius );
+    }
+    if ( !ebbLeaf )
+      radius /= ebbTree->getHeight( /*full=*/true );
+  }
+  else // p outside of box
+  {
+    for ( int i = 1; i <= 3; ++i )
+    {
+      double d = 0;
+      if ( point.Coord(i) < pMin.Coord(i) )
+        d = pMin.Coord(i) - point.Coord(i);
+      else if ( point.Coord(i) > pMax.Coord(i) )
+        d = point.Coord(i) - pMax.Coord(i);
+      if ( d > Precision::Confusion() )
+        radius = Min( d, radius );
+    }
+  }
 
   ElementBndBoxTree::TElemSeq elems;
   ebbTree->getElementsInSphere( p, radius, elems );
   while ( elems.empty() && radius < 1e100 )
   {
-    radius *= 1.5;
+    radius *= 1.1;
     ebbTree->getElementsInSphere( p, radius, elems );
   }
   gp_XYZ proj, bestProj;
   const SMDS_MeshElement* elem = 0;
-  double minDist = 2 * radius;
+  double minDist = Precision::Infinite();
   ElementBndBoxTree::TElemSeq::iterator e = elems.begin();
   for ( ; e != elems.end(); ++e )
   {
@@ -1277,6 +1328,23 @@ gp_XYZ SMESH_ElementSearcherImpl::Project(const gp_Pnt&            point,
       minDist = d;
     }
   }
+  if ( minDist > radius )
+  {
+    ElementBndBoxTree::TElemSeq elems2;
+    ebbTree->getElementsInSphere( p, minDist, elems2 );
+    for ( e = elems2.begin(); e != elems2.end(); ++e )
+    {
+      if ( elems.count( *e ))
+        continue;
+      double d = SMESH_MeshAlgos::GetDistance( *e, point, &proj );
+      if ( d < minDist )
+      {
+        bestProj = proj;
+        elem = *e;
+        minDist = d;
+      }
+    }
+  }
   if ( closestElem ) *closestElem = elem;
 
   return bestProj;
@@ -1609,7 +1677,7 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face,
   try {
     tgtCS = gp_Ax3( xyz[0], OZ, OX );
   }
-  catch ( Standard_Failure ) {
+  catch ( Standard_Failure& ) {
     return badDistance;
   }
   trsf.SetTransformation( tgtCS );
@@ -1916,7 +1984,11 @@ SMESH_MeshAlgos::FindSharpEdges( SMDS_Mesh* theMesh,
   typedef std::pair< bool, const SMDS_MeshNode* >                            TIsSharpAndMedium;
   typedef NCollection_DataMap< SMESH_TLink, TIsSharpAndMedium, SMESH_TLink > TLinkSharpMap;
 
-  TLinkSharpMap linkIsSharp( theMesh->NbFaces() );
+  TLinkSharpMap linkIsSharp;
+  Standard_Integer nbBuckets = FromSmIdType<Standard_Integer>( theMesh->NbFaces() );
+  if ( nbBuckets > 0 )
+    linkIsSharp.ReSize( nbBuckets );
+
   TIsSharpAndMedium sharpMedium( true, 0 );
   bool                 & isSharp = sharpMedium.first;
   const SMDS_MeshNode* & nMedium = sharpMedium.second;
@@ -2023,7 +2095,10 @@ SMESH_MeshAlgos::SeparateFacesByEdges( SMDS_Mesh* theMesh, const std::vector< Ed
 
   typedef std::vector< const SMDS_MeshElement* >                    TFaceVec;
   typedef NCollection_DataMap< SMESH_TLink, TFaceVec, SMESH_TLink > TFacesByLinks;
-  TFacesByLinks facesByLink( theMesh->NbFaces() );
+  TFacesByLinks facesByLink;
+  Standard_Integer nbBuckets = FromSmIdType<Standard_Integer>( theMesh->NbFaces() );
+  if ( nbBuckets > 0 )
+    facesByLink.ReSize( nbBuckets );
 
   std::vector< const SMDS_MeshNode* > faceNodes;
   for ( SMDS_FaceIteratorPtr faceIt = theMesh->facesIterator(); faceIt->more(); )
@@ -2146,10 +2221,26 @@ bool SMESH_MeshAlgos::FaceNormal(const SMDS_MeshElement* F, gp_XYZ& normal, bool
   return ok;
 }
 
-//=======================================================================
-//function : GetCommonNodes
-//purpose  : Return nodes common to two elements
-//=======================================================================
+//================================================================================
+/*!
+ * \brief Return nodes common to two elements
+ */
+//================================================================================
+
+int SMESH_MeshAlgos::NbCommonNodes(const SMDS_MeshElement* e1,
+                                   const SMDS_MeshElement* e2)
+{
+  int nb = 0;
+  for ( int i = 0 ; i < e1->NbNodes(); ++i )
+    nb += ( e2->GetNodeIndex( e1->GetNode( i )) >= 0 );
+  return nb;
+}
+
+//================================================================================
+/*!
+ * \brief Return nodes common to two elements
+ */
+//================================================================================
 
 std::vector< const SMDS_MeshNode*> SMESH_MeshAlgos::GetCommonNodes(const SMDS_MeshElement* e1,
                                                                    const SMDS_MeshElement* e2)
@@ -2160,6 +2251,7 @@ std::vector< const SMDS_MeshNode*> SMESH_MeshAlgos::GetCommonNodes(const SMDS_Me
       common.push_back( e1->GetNode( i ));
   return common;
 }
+
 //================================================================================
 /*!
  * \brief Return true if node1 encounters first in the face and node2, after
@@ -2347,28 +2439,16 @@ void SMESH_MeshAlgos::Get1DBranches( SMDS_ElemIteratorPtr theEdgeIt,
 
     if ( nbBranches == 2 && !startIsBranchEnd ) // join two branches starting at the same node
     {
-      if ( nodeBranches[0].back() == nodeBranches[1].back() )
-      {
-        // it is a closed branch, keep theStartNode first
-        nodeBranches[0].pop_back();
-        nodeBranches[0].reserve( nodeBranches[0].size() + nodeBranches[1].size() );
-        nodeBranches[0].insert( nodeBranches[0].end(),
-                                nodeBranches[1].rbegin(), nodeBranches[1].rend() );
-        branches[0].reserve( branches[0].size() + branches[1].size() );
-        branches[0].insert( branches[0].end(), branches[1].rbegin(), branches[1].rend() );
-      }
-      else
-      {
-        std::reverse( nodeBranches[0].begin(), nodeBranches[0].end() );
-        nodeBranches[0].pop_back();
-        nodeBranches[0].reserve( nodeBranches[0].size() + nodeBranches[1].size() );
-        nodeBranches[0].insert( nodeBranches[0].end(),
-                                nodeBranches[1].begin(), nodeBranches[1].end() );
-
-        std::reverse( branches[0].begin(), branches[0].end() );
-        branches[0].reserve( branches[0].size() + branches[1].size() );
-        branches[0].insert( branches[0].end(), branches[1].begin(), branches[1].end() );
-      }
+      std::reverse( nodeBranches[0].begin(), nodeBranches[0].end() );
+      nodeBranches[0].pop_back();
+      nodeBranches[0].reserve( nodeBranches[0].size() + nodeBranches[1].size() );
+      nodeBranches[0].insert( nodeBranches[0].end(),
+                              nodeBranches[1].begin(), nodeBranches[1].end() );
+
+      std::reverse( branches[0].begin(), branches[0].end() );
+      branches[0].reserve( branches[0].size() + branches[1].size() );
+      branches[0].insert( branches[0].end(), branches[1].begin(), branches[1].end() );
+
       nodeBranches[1].clear();
       branches[1].clear();
     }