Salome HOME
22316: EDF 2719 SMESH: Split hexas into prisms
[modules/smesh.git] / src / SMESHUtils / SMESH_MeshAlgos.cxx
index d0d6008d373f426e51dbfc49d9891fef1ebfb276..24481ca468f26ff6d6af5f5a947ef7482a6686d1 100644 (file)
 
 #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>
@@ -210,7 +212,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     void getElementsInSphere ( const gp_XYZ& center,
                                const double  radius, TIDSortedElemSet& foundElems);
     size_t getSize() { return std::max( _size, _elements.size() ); }
-    ~ElementBndBoxTree();
+    virtual ~ElementBndBoxTree();
 
   protected:
     ElementBndBoxTree():_size(0) {}
@@ -410,6 +412,10 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
  */
 //=======================================================================
 
+SMESH_ElementSearcher::~SMESH_ElementSearcher()
+{
+}
+
 struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
 {
   SMDS_Mesh*                   _mesh;
@@ -423,7 +429,7 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
 
   SMESH_ElementSearcherImpl( SMDS_Mesh& mesh, SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr())
     : _mesh(&mesh),_meshPartIt(elemIt),_ebbTree(0),_nodeSearcher(0),_tolerance(-1),_outerFacesFound(false) {}
-  ~SMESH_ElementSearcherImpl()
+  virtual ~SMESH_ElementSearcherImpl()
   {
     if ( _ebbTree )      delete _ebbTree;      _ebbTree      = 0;
     if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0;
@@ -739,7 +745,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 )
     {
@@ -769,8 +775,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;
@@ -1286,6 +1291,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
@@ -1382,6 +1412,101 @@ 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.
+ *        A not returned bc2 = 1. - bc0 - bc1.
+ *        The point lies within the triangle if ( bc0 >= 0 && bc1 >= 0 && bc0+bc1 <= 1 )
+ */
+//================================================================================
+
+void SMESH_MeshAlgos::GetBarycentricCoords( const gp_XY& p,
+                                            const gp_XY& t0,
+                                            const gp_XY& t1,
+                                            const gp_XY& t2,
+                                            double &     bc0,
+                                            double &     bc1)
+{
+  const double // matrix 2x2
+    T11 = t0.X()-t2.X(), T12 = t1.X()-t2.X(),
+    T21 = t0.Y()-t2.Y(), T22 = t1.Y()-t2.Y();
+  const double Tdet = T11*T22 - T12*T21; // matrix determinant
+  if ( Abs( Tdet ) < std::numeric_limits<double>::min() )
+  {
+    bc0 = bc1 = 2.;
+    return;
+  }
+  // matrix inverse
+  const double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
+  // vector
+  const double r11 = p.X()-t2.X(), r12 = p.Y()-t2.Y();
+  // barycentric coordinates: mutiply matrix by vector
+  bc0 = (t11 * r11 + t12 * r12)/Tdet;
+  bc1 = (t21 * r11 + t22 * r12)/Tdet;
+}
+
 //=======================================================================
 //function : FindFaceInSet
 //purpose  : Return a face having linked nodes n1 and n2 and which is