Salome HOME
23173: EDF 11552 - Problem using Add 0D element function
[modules/smesh.git] / src / SMESHUtils / SMESH_MeshAlgos.cxx
index 2ecd29a9044f5b483c2645a8b13e7be1a4b700cf..a5171703b8f692055813d3ef2a70c29d2d364510 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2015  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
@@ -6,7 +6,7 @@
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 
 #include "SMESH_MeshAlgos.hxx"
 
+#include "SMDS_FaceOfNodes.hxx"
 #include "SMDS_LinearEdge.hxx"
-#include "SMDS_VolumeTool.hxx"
 #include "SMDS_Mesh.hxx"
+#include "SMDS_PolygonalFaceOfNodes.hxx"
+#include "SMDS_VolumeTool.hxx"
 #include "SMESH_OctreeNode.hxx"
 
 #include <GC_MakeSegment.hxx>
@@ -166,6 +168,18 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
     return closestNode;
   }
 
+  //---------------------------------------------------------------------
+  /*!
+   * \brief Finds nodes located within a tolerance near a point 
+   */
+  int FindNearPoint(const gp_Pnt&                        point,
+                    const double                         tolerance,
+                    std::vector< const SMDS_MeshNode* >& foundNodes)
+  {
+    myOctreeNode->NodesAround( point.Coord(), foundNodes, tolerance );
+    return foundNodes.size();
+  }
+
   //---------------------------------------------------------------------
   /*!
    * \brief Destructor
@@ -425,8 +439,10 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   bool                         _outerFacesFound;
   set<const SMDS_MeshElement*> _outerFaces; // empty means "no internal faces at all"
 
-  SMESH_ElementSearcherImpl( SMDS_Mesh& mesh, SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr())
-    : _mesh(&mesh),_meshPartIt(elemIt),_ebbTree(0),_nodeSearcher(0),_tolerance(-1),_outerFacesFound(false) {}
+  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) {}
   virtual ~SMESH_ElementSearcherImpl()
   {
     if ( _ebbTree )      delete _ebbTree;      _ebbTree      = 0;
@@ -694,21 +710,21 @@ FindElementsByPoint(const gp_Pnt&                      point,
     if ( !_nodeSearcher )
       _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
 
-    const SMDS_MeshNode* closeNode = _nodeSearcher->FindClosestTo( point );
-    if ( !closeNode ) return foundElements.size();
-
-    if ( point.Distance( SMESH_TNodeXYZ( closeNode )) > tolerance )
-      return foundElements.size(); // to far from any node
+    std::vector< const SMDS_MeshNode* > foundNodes;
+    _nodeSearcher->FindNearPoint( point, tolerance, foundNodes );
 
     if ( type == SMDSAbs_Node )
     {
-      foundElements.push_back( closeNode );
+      foundElements.assign( foundNodes.begin(), foundNodes.end() );
     }
     else
     {
-      SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( type );
-      while ( elemIt->more() )
-        foundElements.push_back( elemIt->next() );
+      for ( size_t i = 0; i < foundNodes.size(); ++i )
+      {
+        SMDS_ElemIteratorPtr elemIt = foundNodes[i]->GetInverseElementIterator( type );
+        while ( elemIt->more() )
+          foundElements.push_back( elemIt->next() );
+      }
     }
   }
   // =================================================================================
@@ -743,7 +759,7 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
 {
   const SMDS_MeshElement* closestElem = 0;
 
-  if ( type == SMDSAbs_Face )
+  if ( type == SMDSAbs_Face || type == SMDSAbs_Volume )
   {
     if ( !_ebbTree || _elementType != type )
     {
@@ -757,10 +773,10 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
     {
       gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox()->CornerMin() +
                                  _ebbTree->getBox()->CornerMax() );
-      double radius;
+      double radius = -1;
       if ( _ebbTree->getBox()->IsOut( point.XYZ() ))
         radius = point.Distance( boxCenter ) - 0.5 * _ebbTree->maxSize();
-      else
+      if ( radius < 0 )
         radius = _ebbTree->maxSize() / pow( 2., _ebbTree->getHeight()) / 2;
       while ( suspectElems.empty() )
       {
@@ -773,8 +789,7 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
     TIDSortedElemSet::iterator elem = suspectElems.begin();
     for ( ; elem != suspectElems.end(); ++elem )
     {
-      double dist = SMESH_MeshAlgos::GetDistance( dynamic_cast<const SMDS_MeshFace*>(*elem),
-                                                   point );
+      double dist = SMESH_MeshAlgos::GetDistance( *elem, point );
       if ( dist < minDist + 1e-10)
       {
         minDist = dist;
@@ -1078,32 +1093,22 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
 
   // get ordered nodes
 
-  vector< gp_XYZ > xyz;
-  vector<const SMDS_MeshNode*> nodeList;
+  vector< SMESH_TNodeXYZ > xyz;
 
-  SMDS_ElemIteratorPtr nodeIt = element->nodesIterator();
-  if ( element->IsQuadratic() ) {
-    nodeIt = element->interlacedNodesElemIterator();
-    // if (const SMDS_VtkFace* f=dynamic_cast<const SMDS_VtkFace*>(element))
-    //   nodeIt = f->interlacedNodesElemIterator();
-    // else if (const SMDS_VtkEdge*  e =dynamic_cast<const SMDS_VtkEdge*>(element))
-    //   nodeIt = e->interlacedNodesElemIterator();
-  }
+  SMDS_ElemIteratorPtr nodeIt = element->interlacedNodesElemIterator();
   while ( nodeIt->more() )
   {
     SMESH_TNodeXYZ node = nodeIt->next();
     xyz.push_back( node );
-    nodeList.push_back(node._node);
   }
 
-  int i, nbNodes = (int) nodeList.size(); // central node of biquadratic is missing
+  int i, nbNodes = (int) xyz.size(); // central node of biquadratic is missing
 
   if ( element->GetType() == SMDSAbs_Face ) // --------------------------------------------------
   {
     // compute face normal
     gp_Vec faceNorm(0,0,0);
     xyz.push_back( xyz.front() );
-    nodeList.push_back( nodeList.front() );
     for ( i = 0; i < nbNodes; ++i )
     {
       gp_Vec edge1( xyz[i+1], xyz[i]);
@@ -1116,7 +1121,7 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
       // degenerated face: point is out if it is out of all face edges
       for ( i = 0; i < nbNodes; ++i )
       {
-        SMDS_LinearEdge edge( nodeList[i], nodeList[i+1] );
+        SMDS_LinearEdge edge( xyz[i]._node, xyz[i+1]._node );
         if ( !IsOut( &edge, point, tol ))
           return false;
       }
@@ -1203,13 +1208,22 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
     // (we consider quadratic edge as being composed of two straight parts)
     for ( i = 1; i < nbNodes; ++i )
     {
-      gp_Vec edge( xyz[i-1], xyz[i]);
-      gp_Vec n1p ( xyz[i-1], point);
-      double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude();
-      if ( dist > tol )
+      gp_Vec edge( xyz[i-1], xyz[i] );
+      gp_Vec n1p ( xyz[i-1], point  );
+      double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
+      if ( u <= 0. ) {
+        if ( n1p.SquareMagnitude() < tol * tol )
+          return false;
         continue;
-      gp_Vec n2p( xyz[i], point );
-      if ( fabs( edge.Magnitude() - n1p.Magnitude() - n2p.Magnitude()) > tol )
+      }
+      if ( u >= 1. ) {
+        if ( point.SquareDistance( xyz[i] ) < tol * tol )
+          return false;
+        continue;
+      }
+      gp_XYZ proj = ( 1. - u ) * xyz[i-1] + u * xyz[i]; // projection of the point on the edge
+      double dist2 = point.SquareDistance( proj );
+      if ( dist2 > tol * tol )
         continue;
       return false; // point is ON this part
     }
@@ -1218,7 +1232,7 @@ bool SMESH_MeshAlgos::IsOut( const SMDS_MeshElement* element, const gp_Pnt& poin
   // Node or 0D element -------------------------------------------------------------------------
   {
     gp_Vec n2p ( xyz[0], point );
-    return n2p.Magnitude() <= tol;
+    return n2p.SquareMagnitude() <= tol * tol;
   }
   return true;
 }
@@ -1290,6 +1304,31 @@ namespace
   }
 }
 
+//=======================================================================
+/*!
+ * \brief Return minimal distance from a point to an element
+ *
+ * Currently we ignore non-planarity and 2nd order of face
+ */
+//=======================================================================
+
+double SMESH_MeshAlgos::GetDistance( const SMDS_MeshElement* elem,
+                                     const gp_Pnt&           point )
+{
+  switch ( elem->GetType() )
+  {
+  case SMDSAbs_Volume:
+    return GetDistance( dynamic_cast<const SMDS_MeshVolume*>( elem ), point);
+  case SMDSAbs_Face:
+    return GetDistance( dynamic_cast<const SMDS_MeshFace*>( elem ), point);
+  case SMDSAbs_Edge:
+    return GetDistance( dynamic_cast<const SMDS_MeshEdge*>( elem ), point);
+  case SMDSAbs_Node:
+    return point.Distance( SMESH_TNodeXYZ( elem ));
+  }
+  return -1;
+}
+
 //=======================================================================
 /*!
  * \brief Return minimal distance from a point to a face
@@ -1386,6 +1425,68 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face,
   return badDistance;
 }
 
+//=======================================================================
+/*!
+ * \brief Return minimal distance from a point to an edge
+ */
+//=======================================================================
+
+double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* edge, const gp_Pnt& point )
+{
+  throw SALOME_Exception(LOCALIZED("not implemented so far"));
+}
+
+//=======================================================================
+/*!
+ * \brief Return minimal distance from a point to a volume
+ *
+ * Currently we ignore non-planarity and 2nd order
+ */
+//=======================================================================
+
+double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt& point )
+{
+  SMDS_VolumeTool vTool( volume );
+  vTool.SetExternalNormal();
+  const int iQ = volume->IsQuadratic() ? 2 : 1;
+
+  double n[3], bc[3];
+  double minDist = 1e100, dist;
+  for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+  {
+    // skip a facet with normal not "looking at" the point
+    if ( !vTool.GetFaceNormal( iF, n[0], n[1], n[2] ) ||
+         !vTool.GetFaceBaryCenter( iF, bc[0], bc[1], bc[2] ))
+      continue;
+    gp_XYZ bcp = point.XYZ() - gp_XYZ( bc[0], bc[1], bc[2] );
+    if ( gp_XYZ( n[0], n[1], n[2] ) * bcp < 1e-6 )
+      continue;
+
+    // find distance to a facet
+    const SMDS_MeshNode** nodes = vTool.GetFaceNodes( iF );
+    switch ( vTool.NbFaceNodes( iF ) / iQ ) {
+    case 3:
+    {
+      SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ] );
+      dist = GetDistance( &tmpFace, point );
+      break;
+    }
+    case 4:
+    {
+      SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ], nodes[ 3*iQ ]);
+      dist = GetDistance( &tmpFace, point );
+      break;
+    }
+    default:
+      vector<const SMDS_MeshNode *> nvec( nodes, nodes + vTool.NbFaceNodes( iF ));
+      SMDS_PolygonalFaceOfNodes tmpFace( nvec );
+      dist = GetDistance( &tmpFace, point );
+    }
+    minDist = Min( minDist, dist );
+  }
+  return minDist;
+}
+
 //================================================================================
 /*!
  * \brief Returns barycentric coordinates of a point within a triangle.
@@ -1440,10 +1541,8 @@ SMESH_MeshAlgos::FindFaceInSet(const SMDS_MeshNode*    n1,
   const SMDS_MeshElement* face = 0;
 
   SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
-  //MESSAGE("n1->GetInverseElementIterator(SMDSAbs_Face) " << invElemIt);
   while ( invElemIt->more() && !face ) // loop on inverse faces of n1
   {
-    //MESSAGE("in while ( invElemIt->more() && !face )");
     const SMDS_MeshElement* elem = invElemIt->next();
     if (avoidSet.count( elem ))
       continue;
@@ -1462,9 +1561,6 @@ SMESH_MeshAlgos::FindFaceInSet(const SMDS_MeshNode*    n1,
     if ( !face && elem->IsQuadratic())
     {
       // analysis for quadratic elements using all nodes
-      // const SMDS_VtkFace* F = dynamic_cast<const SMDS_VtkFace*>(elem);
-      // if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
-      // use special nodes iterator
       SMDS_ElemIteratorPtr anIter = elem->interlacedNodesElemIterator();
       const SMDS_MeshNode* prevN = static_cast<const SMDS_MeshNode*>( anIter->next() );
       for ( i1 = -1, i2 = 0; anIter->more() && !face; i1++, i2++ )
@@ -1499,7 +1595,7 @@ bool SMESH_MeshAlgos::FaceNormal(const SMDS_MeshElement* F, gp_XYZ& normal, bool
     return false;
 
   normal.SetCoord(0,0,0);
-  int nbNodes = F->IsQuadratic() ? F->NbNodes()/2 : F->NbNodes();
+  int nbNodes = F->NbCornerNodes();
   for ( int i = 0; i < nbNodes-2; ++i )
   {
     gp_XYZ p[3];
@@ -1550,9 +1646,10 @@ SMESH_NodeSearcher* SMESH_MeshAlgos::GetNodeSearcher(SMDS_Mesh& mesh)
  */
 //=======================================================================
 
-SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh& mesh)
+SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh& mesh,
+                                                           double     tolerance)
 {
-  return new SMESH_ElementSearcherImpl( mesh );
+  return new SMESH_ElementSearcherImpl( mesh, tolerance );
 }
 
 //=======================================================================
@@ -1562,7 +1659,8 @@ SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh& mesh)
 //=======================================================================
 
 SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh&           mesh,
-                                                           SMDS_ElemIteratorPtr elemIt)
+                                                           SMDS_ElemIteratorPtr elemIt,
+                                                           double               tolerance)
 {
-  return new SMESH_ElementSearcherImpl( mesh, elemIt );
+  return new SMESH_ElementSearcherImpl( mesh, tolerance, elemIt );
 }