Salome HOME
Merge commit '6600bcec782fc8b6c72871fe6e08bd19a34a4e2b'
[modules/smesh.git] / src / SMESHUtils / SMESH_MeshAlgos.cxx
index db6f5e690b5eb2d33490c7ddcbe316e7da973ef3..e8b1187986c2b9d48745e3cd29fe0fa896521c6c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2020  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2023  CEA, EDF, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -253,6 +253,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
       void init(const SMDS_MeshElement* elem, double tolerance);
     };
     std::vector< ElementBox* > _elements;
+    //std::string _name;
 
     typedef ObjectPool< ElementBox > TElementBoxPool;
 
@@ -281,15 +282,13 @@ 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;
 
-#ifdef _DEBUG_
-    if ( theElemIt && !theElemIt->more() )
+    if (SALOME::VerbosityActivated() && theElemIt && !theElemIt->more() )
       std::cout << "WARNING: ElementBndBoxTree constructed on empty iterator!" << std::endl;
-#endif
 
     SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType );
     while ( elemIt->more() )
@@ -342,6 +341,15 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 
       if ( child->isLeaf() && child->_elements.capacity() > child->_elements.size() )
         SMESHUtils::CompactVector( child->_elements );
+
+      // child->_name = _name + char('0' + j);
+      // if ( child->isLeaf() && child->_elements.size() )
+      // {
+      //   cout << child->_name << " ";
+      //   for ( size_t i = 0; i < child->_elements.size(); ++i )
+      //     cout << child->_elements[i]->_element->GetID() << " ";
+      //   cout << endl;
+      // }
     }
   }
 
@@ -546,7 +554,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,
@@ -687,6 +696,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;
@@ -1673,7 +1685,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 );
@@ -1980,7 +1992,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;
@@ -2087,7 +2103,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(); )
@@ -2241,6 +2260,56 @@ std::vector< const SMDS_MeshNode*> SMESH_MeshAlgos::GetCommonNodes(const SMDS_Me
   return common;
 }
 
+//================================================================================
+/*!
+ * \brief Return true if a node is on a boundary of 2D mesh.
+ *        Optionally returns two neighboring boundary nodes (or more in non-manifold mesh)
+ */
+//================================================================================
+
+bool SMESH_MeshAlgos::IsOn2DBoundary( const SMDS_MeshNode*                 theNode,
+                                      std::vector< const SMDS_MeshNode*> * theNeibors )
+{
+  typedef NCollection_DataMap< SMESH_TLink, int, SMESH_TLink > TLinkCountMap;
+  TLinkCountMap linkCountMap( 10 );
+
+  int nbFreeLinks = 0;
+  for ( SMDS_ElemIteratorPtr fIt = theNode->GetInverseElementIterator(SMDSAbs_Face); fIt->more(); )
+  {
+    const SMDS_MeshElement* face = fIt->next();
+    const int          nbCorners = face->NbCornerNodes();
+
+    int    iN = face->GetNodeIndex( theNode );
+    int iPrev = ( iN - 1 + nbCorners ) % nbCorners;
+    int iNext = ( iN + 1 ) % nbCorners;
+
+    for ( int i : { iPrev, iNext } )
+    {
+      SMESH_TLink link( theNode, face->GetNode( i ));
+      int* count = linkCountMap.ChangeSeek( link );
+      if ( count )  ++( *count );
+      else          linkCountMap.Bind( link, 1 );
+
+      if ( !count ) ++nbFreeLinks;
+      else          --nbFreeLinks;
+    }
+  }
+
+  if ( theNeibors )
+  {
+    theNeibors->clear();
+    theNeibors->reserve( nbFreeLinks );
+    for ( TLinkCountMap::Iterator linkIt( linkCountMap ); linkIt.More(); linkIt.Next() )
+      if ( linkIt.Value() == 1 )
+      {
+        theNeibors->push_back( linkIt.Key().node1() );
+        if ( theNeibors->back() == theNode )
+          theNeibors->back() = linkIt.Key().node2();
+      }
+  }
+  return nbFreeLinks > 0;
+}
+
 //================================================================================
 /*!
  * \brief Return true if node1 encounters first in the face and node2, after
@@ -2505,3 +2574,111 @@ SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh&
 {
   return new SMESH_ElementSearcherImpl( mesh, tolerance, elemIt );
 }
+
+
+//================================================================================
+/*!
+ * \brief Intersect a ray with a convex volume
+ *  \param [in] ray - the ray
+ *  \param [in] rayLen - ray length
+ *  \param [in] vol - the volume
+ *  \param [out] tMin - return a ray parameter where the ray enters the volume
+ *  \param [out] tMax - return a ray parameter where the ray exit the volume
+ *  \param [out] iFacetMin - facet index where the ray enters the volume
+ *  \param [out] iFacetMax - facet index  where the ray exit the volume
+ *  \return bool - true if the ray intersects the volume
+ */
+//================================================================================
+
+bool SMESH_MeshAlgos::IntersectRayVolume( const gp_Ax1& ray,
+                                          const double rayLen,
+                                          const SMDS_MeshElement* vol,
+                                          double & tMin,
+                                          double & tMax,
+                                          int & iFacetMin,
+                                          int & iFacetMax)
+{
+  /* Ray-Convex Polyhedron Intersection Test by Eric Haines, erich@eye.com
+   *
+   * This test checks the ray against each face of a polyhedron, checking whether
+   * the set of intersection points found for each ray-plane intersection
+   * overlaps the previous intersection results.  If there is no overlap (i.e.
+   * no line segment along the ray that is inside the polyhedron), then the
+   * ray misses and returns false; else true.
+   */
+  SMDS_VolumeTool vTool;
+  if ( !vTool.Set( vol ))
+    return false;
+
+  tMin = -Precision::Infinite() ;
+  tMax = rayLen ;
+
+  /* Test each plane in polyhedron */
+  for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+  {
+    const SMDS_MeshNode** fNodes = vTool.GetFaceNodes( iF );
+    gp_XYZ normal;
+    vTool.GetFaceNormal( iF,
+                         normal.ChangeCoord(1),
+                         normal.ChangeCoord(2),
+                         normal.ChangeCoord(3) );
+    double D = - ( normal * SMESH_NodeXYZ( fNodes[0] ));
+
+    /* Compute intersection point T and sidedness */
+    double vd = ray.Direction().XYZ() * normal;
+    double vn = ray.Location().XYZ() * normal + D;
+    if ( vd == 0.0 ) {
+      /* ray is parallel to plane - check if ray origin is inside plane's
+         half-space */
+      if ( vn > 0.0 )
+        /* ray origin is outside half-space */
+        return false;
+    }
+    else
+    {
+      /* ray not parallel - get distance to plane */
+      double t = -vn / vd ;
+      if ( vd < 0.0 )
+      {
+        /* front face - T is a near point */
+        if ( t > tMax ) return false;
+        if ( t > tMin ) {
+          /* hit near face */
+          tMin = t ;
+          iFacetMin = iF;
+        }
+      }
+      else
+      {
+        /* back face - T is a far point */
+        if ( t < tMin ) return false;
+        if ( t < tMax ) {
+          /* hit far face */
+          tMax = t ;
+          iFacetMax = iF;
+        }
+      }
+    }
+  }
+
+  /* survived all tests */
+  /* Note: if ray originates on polyhedron, may want to change 0.0 to some
+   * epsilon to avoid intersecting the originating face.
+   */
+  if ( tMin >= 0.0 ) {
+    /* outside, hitting front face */
+    return true;
+  }
+  else
+  {
+    if ( tMax < rayLen ) {
+      /* inside, hitting back face */
+      return true;
+    }
+    else
+    {
+      /* inside, but back face beyond tmax */
+      return false;
+    }
+  }
+}