+//=======================================================================
+/*!
+ * \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;
+}
+