X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESHUtils%2FSMESH_MeshAlgos.cxx;h=a5171703b8f692055813d3ef2a70c29d2d364510;hp=2ecd29a9044f5b483c2645a8b13e7be1a4b700cf;hb=674c0d8b9d98776136d216ec8f1bad56acac5bf5;hpb=f7aba4830d53719b963fdb7fccc98b760fdef2d1 diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 2ecd29a90..a5171703b 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-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 @@ -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 @@ -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 _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(*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 nodeList; + vector< SMESH_TNodeXYZ > xyz; - SMDS_ElemIteratorPtr nodeIt = element->nodesIterator(); - if ( element->IsQuadratic() ) { - nodeIt = element->interlacedNodesElemIterator(); - // if (const SMDS_VtkFace* f=dynamic_cast(element)) - // nodeIt = f->interlacedNodesElemIterator(); - // else if (const SMDS_VtkEdge* e =dynamic_cast(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( 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 +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 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(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( 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 ); }