X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FSMESHUtils%2FSMESH_MeshAlgos.cxx;h=957828474e306b0d3bb07286179f454aaf297c1c;hb=69f064e1222d71a85b8ac1e11d82e0480b6d6cea;hp=45acf33b115cb4067aa23516ee4a774f59d5fdec;hpb=ef59152514df17c64ad01c8abab01331618a1fc4;p=modules%2Fsmesh.git diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 45acf33b1..957828474 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2014 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 @@ -28,9 +28,11 @@ #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 @@ -743,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 ) { @@ -757,10 +759,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 +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(*elem), - point ); + double dist = SMESH_MeshAlgos::GetDistance( *elem, point ); if ( dist < minDist + 1e-10) { minDist = dist; @@ -1290,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( elem ), point); + case SMDSAbs_Face: + return GetDistance( dynamic_cast( elem ), point); + case SMDSAbs_Edge: + return GetDistance( dynamic_cast( 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 +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 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::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