Salome HOME
Minor doc changes
[modules/smesh.git] / src / SMESHUtils / SMESH_MeshAlgos.cxx
index c42f020236e096a1744df1235270a1d0f3d4661b..e65e20114c42723fd8ec3ff435c115a73205332c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  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
@@ -60,7 +60,8 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
   /*!
    * \brief Constructor
    */
-  SMESH_NodeSearcherImpl( const SMDS_Mesh* theMesh )
+  SMESH_NodeSearcherImpl( const SMDS_Mesh*     theMesh   = 0,
+                          SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr() )
   {
     myMesh = ( SMDS_Mesh* ) theMesh;
 
@@ -70,6 +71,14 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
       while ( nIt->more() )
         nodes.insert( nodes.end(), nIt->next() );
     }
+    else if ( theElemIt )
+    {
+      while ( theElemIt->more() )
+      {
+        const SMDS_MeshElement* e = theElemIt->next();
+        nodes.insert( e->begin_nodes(), e->end_nodes() );
+      }
+    }
     myOctreeNode = new SMESH_OctreeNode(nodes) ;
 
     // get max size of a leaf box
@@ -458,6 +467,10 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   void GetElementsNearLine( const gp_Ax1&                      line,
                             SMDSAbs_ElementType                type,
                             vector< const SMDS_MeshElement* >& foundElems);
+  void GetElementsInSphere( const gp_XYZ&                      center,
+                            const double                       radius,
+                            SMDSAbs_ElementType                type,
+                            vector< const SMDS_MeshElement* >& foundElems);
   double getTolerance();
   bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
                             const double tolerance, double & param);
@@ -466,6 +479,7 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   {
     return _outerFaces.empty() || _outerFaces.count(face);
   }
+
   struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
   {
     const SMDS_MeshElement* _face;
@@ -570,7 +584,7 @@ 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);
+    anExtCC.Init( lineCurve, edge.Value() );
     if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol)
     {
       Quantity_Parameter pl, pe;
@@ -646,6 +660,7 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF
         set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin();
         for ( ; face != faces.end(); ++face )
         {
+          if ( *face == outerFace ) continue;
           if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false ))
             continue;
           gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2;
@@ -660,8 +675,8 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF
     // store the found outer face and add its links to continue seaching from
     if ( outerFace2 )
     {
-      _outerFaces.insert( outerFace );
-      int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 );
+      _outerFaces.insert( outerFace2 );
+      int nbNodes = outerFace2->NbCornerNodes();
       for ( int i = 0; i < nbNodes; ++i )
       {
         SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes));
@@ -708,8 +723,12 @@ FindElementsByPoint(const gp_Pnt&                      point,
   if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement || type == SMDSAbs_Ball)
   {
     if ( !_nodeSearcher )
-      _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
-
+    {
+      if ( _meshPartIt )
+        _nodeSearcher = new SMESH_NodeSearcherImpl( 0, _meshPartIt );
+      else
+        _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
+    }
     std::vector< const SMDS_MeshNode* > foundNodes;
     _nodeSearcher->FindNearPoint( point, tolerance, foundNodes );
 
@@ -887,8 +906,9 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
       }
       else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 )
       {
+        double tol = 1e-4 * Sqrt( fNorm.Modulus() );
         gp_Pnt intersectionPoint = intersection.Point(1);
-        if ( !SMESH_MeshAlgos::IsOut( *face, intersectionPoint, tolerance ))
+        if ( !SMESH_MeshAlgos::IsOut( *face, intersectionPoint, tol ))
           u2inters.insert(make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm )));
       }
     }
@@ -1078,6 +1098,27 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1&
   foundElems.assign( suspectFaces.begin(), suspectFaces.end());
 }
 
+//=======================================================================
+/*
+ * Return elements whose bounding box intersects a sphere
+ */
+//=======================================================================
+
+void SMESH_ElementSearcherImpl::GetElementsInSphere( const gp_XYZ&                      center,
+                                                     const double                       radius,
+                                                     SMDSAbs_ElementType                type,
+                                                     vector< const SMDS_MeshElement* >& foundElems)
+{
+  if ( !_ebbTree || _elementType != type )
+  {
+    if ( _ebbTree ) delete _ebbTree;
+    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
+  }
+  TIDSortedElemSet suspectFaces; // elements possibly intersecting the line
+  _ebbTree->getElementsInSphere( center, radius, suspectFaces );
+  foundElems.assign( suspectFaces.begin(), suspectFaces.end() );
+}
+
 //=======================================================================
 /*!
  * \brief Return true if the point is IN or ON of the element
@@ -1093,14 +1134,11 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
 
   // get ordered nodes
 
-  vector< SMESH_TNodeXYZ > xyz;
+  vector< SMESH_TNodeXYZ > xyz; xyz.reserve( element->NbNodes()+1 );
 
   SMDS_ElemIteratorPtr nodeIt = element->interlacedNodesElemIterator();
-  while ( nodeIt->more() )
-  {
-    SMESH_TNodeXYZ node = nodeIt->next();
-    xyz.push_back( node );
-  }
+  for ( int i = 0; nodeIt->more(); ++i )
+    xyz.push_back( SMESH_TNodeXYZ( nodeIt->next() ));
 
   int i, nbNodes = (int) xyz.size(); // central node of biquadratic is missing
 
@@ -1131,8 +1169,23 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
 
     // check if the point lays on face plane
     gp_Vec n2p( xyz[0], point );
-    if ( fabs( n2p * faceNorm ) > tol )
-      return true; // not on face plane
+    double dot = n2p * faceNorm;
+    if ( Abs( dot ) > tol ) // not on face plane
+    {
+      bool isOut = true;
+      if ( nbNodes > 3 ) // maybe the face is not planar
+      {
+        double elemThick = 0;
+        for ( i = 1; i < nbNodes; ++i )
+        {
+          gp_Vec n2n( xyz[0], xyz[i] );
+          elemThick = Max( elemThick, Abs( n2n * faceNorm ));
+        }
+        isOut = Abs( dot ) > elemThick + tol;
+      }
+      if ( isOut )
+        return isOut;
+    }
 
     // check if point is out of face boundary:
     // define it by closest transition of a ray point->infinity through face boundary
@@ -1144,7 +1197,7 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     normSize = plnNorm.Magnitude();
     if ( normSize <= tol ) return false; // point coincides with the first node
     plnNorm /= normSize;
-    // for each node of the face, compute its signed distance to the plane
+    // 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 )
     {
@@ -1154,12 +1207,12 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     dist.back() = dist.front();
     // find the closest intersection
     int    iClosest = -1;
-    double rClosest, distClosest = 1e100;;
+    double rClosest = 0, distClosest = 1e100;
     gp_Pnt pClosest;
     for ( i = 0; i < nbNodes; ++i )
     {
       double r;
-      if ( fabs( dist[i]) < tol )
+      if ( fabs( dist[i] ) < tol )
         r = 0.;
       else if ( fabs( dist[i+1]) < tol )
         r = 1.;
@@ -1168,17 +1221,14 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
       else
         continue; // no intersection
       gp_Pnt pInt = xyz[i] * (1.-r) + xyz[i+1] * r;
-      gp_Vec p2int ( point, pInt);
-      if ( p2int * ray > -tol ) // right half-space
+      gp_Vec p2int( point, pInt);
+      double intDist = p2int.SquareMagnitude();
+      if ( intDist < distClosest )
       {
-        double intDist = p2int.SquareMagnitude();
-        if ( intDist < distClosest )
-        {
-          iClosest = i;
-          rClosest = r;
-          pClosest = pInt;
-          distClosest = intDist;
-        }
+        iClosest = i;
+        rClosest = r;
+        pClosest = pInt;
+        distClosest = intDist;
       }
     }
     if ( iClosest < 0 )
@@ -1202,6 +1252,7 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     bool covexCorner = ( edgeNorm * edgeAdjacent * (rClosest==1. ? 1. : -1.)) < 0;
     return covexCorner ? (out || out2) : (out && out2);
   }
+
   if ( element->GetType() == SMDSAbs_Edge ) // --------------------------------------------------
   {
     // point is out of edge if it is NOT ON any straight part of edge
@@ -1229,10 +1280,11 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     }
     return true;
   }
+
   // Node or 0D element -------------------------------------------------------------------------
   {
     gp_Vec n2p ( xyz[0], point );
-    return n2p.SquareMagnitude() <= tol * tol;
+    return n2p.SquareMagnitude() > tol * tol;
   }
   return true;
 }
@@ -1668,6 +1720,17 @@ SMESH_NodeSearcher* SMESH_MeshAlgos::GetNodeSearcher(SMDS_Mesh& mesh)
   return new SMESH_NodeSearcherImpl( &mesh );
 }
 
+//=======================================================================
+/*!
+ * \brief Return SMESH_NodeSearcher
+ */
+//=======================================================================
+
+SMESH_NodeSearcher* SMESH_MeshAlgos::GetNodeSearcher(SMDS_ElemIteratorPtr elemIt)
+{
+  return new SMESH_NodeSearcherImpl( 0, elemIt );
+}
+
 //=======================================================================
 /*!
  * \brief Return SMESH_ElementSearcher