X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FSMESHUtils%2FSMESH_MeshAlgos.cxx;h=24958cc81fd77b243c71cd3142d8caf096d5b2d7;hb=HEAD;hp=55c5f50253bed26288668235f25414c9e26e4549;hpb=6ca4db2d7c49474a93e30c02bef83f05b354e540;p=modules%2Fsmesh.git diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 55c5f5025..fac231ec8 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2019 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2024 CEA, EDF, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -238,6 +238,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() void getElementsInBox ( const Bnd_B3d& box, TElemSeq& foundElems ); void getElementsInSphere ( const gp_XYZ& center, const double radius, TElemSeq& foundElems ); ElementBndBoxTree* getLeafAtPoint( const gp_XYZ& point ); + int getNbElements(); protected: ElementBndBoxTree() {} @@ -252,6 +253,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() void init(const SMDS_MeshElement* elem, double tolerance); }; std::vector< ElementBox* > _elements; + //std::string _name; typedef ObjectPool< ElementBox > TElementBoxPool; @@ -280,15 +282,13 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() double tolerance) :SMESH_Octree( new LimitAndPool() ) { - int nbElems = mesh.GetMeshInfo().NbElements( elemType ); + smIdType nbElems = mesh.GetMeshInfo().NbElements( elemType ); _elements.reserve( nbElems ); TElementBoxPool& elBoPool = getLimitAndPool()->_elBoPool; -#ifdef _DEBUG_ - if ( theElemIt && !theElemIt->more() ) + if (SALOME::VerbosityActivated() && theElemIt && !theElemIt->more() ) std::cout << "WARNING: ElementBndBoxTree constructed on empty iterator!" << std::endl; -#endif SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType ); while ( elemIt->more() ) @@ -316,7 +316,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() //================================================================================ /*! - * \brief Redistrubute element boxes among children + * \brief Redistribute element boxes among children */ //================================================================================ @@ -341,6 +341,15 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() if ( child->isLeaf() && child->_elements.capacity() > child->_elements.size() ) SMESHUtils::CompactVector( child->_elements ); + + // child->_name = _name + char('0' + j); + // if ( child->isLeaf() && child->_elements.size() ) + // { + // cout << child->_name << " "; + // for ( size_t i = 0; i < child->_elements.size(); ++i ) + // cout << child->_elements[i]->_element->GetID() << " "; + // cout << endl; + // } } } @@ -466,6 +475,27 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() return 0; } + //================================================================================ + /*! + * \brief Return number of elements + */ + //================================================================================ + + int ElementBndBoxTree::getNbElements() + { + int nb = 0; + if ( isLeaf() ) + { + nb = _elements.size(); + } + else + { + for (int i = 0; i < 8; i++) + nb += ((ElementBndBoxTree*) myChildren[i])->getNbElements(); + } + return nb; + } + //================================================================================ /*! * \brief Construct the element box @@ -524,7 +554,8 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher { delete _ebbTree[i]; _ebbTree[i] = NULL; } - if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0; + if ( _nodeSearcher ) delete _nodeSearcher; + _nodeSearcher = 0; } virtual int FindElementsByPoint(const gp_Pnt& point, SMDSAbs_ElementType type, @@ -665,6 +696,9 @@ bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin& lin GC_MakeSegment edge( SMESH_TNodeXYZ( face->GetNode( i )), SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) )); anExtCC.Init( lineCurve, edge.Value() ); + if ( !anExtCC.Extrema().IsDone() || + anExtCC.Extrema().IsParallel() ) + continue; if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol) { Standard_Real pl, pe; @@ -848,7 +882,7 @@ FindElementsByPoint(const gp_Pnt& point, /*! * \brief Find an element of given type most close to the given point * - * WARNING: Only face search is implemeneted so far + * WARNING: Only edge, face and volume search is implemented so far */ //======================================================================= @@ -1076,7 +1110,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) const SMDS_MeshElement* prevFace = u_int1->second._face; while ( ok && u_int2->second._coincides ) { - if ( SMESH_MeshAlgos::GetCommonNodes(prevFace , u_int2->second._face).empty() ) + if ( SMESH_MeshAlgos::NbCommonNodes(prevFace , u_int2->second._face) == 0 ) ok = false; else { @@ -1302,6 +1336,23 @@ gp_XYZ SMESH_ElementSearcherImpl::Project(const gp_Pnt& point, minDist = d; } } + if ( minDist > radius ) + { + ElementBndBoxTree::TElemSeq elems2; + ebbTree->getElementsInSphere( p, minDist, elems2 ); + for ( e = elems2.begin(); e != elems2.end(); ++e ) + { + if ( elems.count( *e )) + continue; + double d = SMESH_MeshAlgos::GetDistance( *e, point, &proj ); + if ( d < minDist ) + { + bestProj = proj; + elem = *e; + minDist = d; + } + } + } if ( closestElem ) *closestElem = elem; return bestProj; @@ -1634,7 +1685,7 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face, try { tgtCS = gp_Ax3( xyz[0], OZ, OX ); } - catch ( Standard_Failure ) { + catch ( Standard_Failure& ) { return badDistance; } trsf.SetTransformation( tgtCS ); @@ -1938,10 +1989,14 @@ SMESH_MeshAlgos::FindSharpEdges( SMDS_Mesh* theMesh, std::vector< Edge > resultEdges; if ( !theMesh ) return resultEdges; - typedef std::pair< bool, const SMDS_MeshNode* > TIsSharpAndMedium; - typedef NCollection_DataMap< SMESH_TLink, TIsSharpAndMedium, SMESH_TLink > TLinkSharpMap; + typedef std::pair< bool, const SMDS_MeshNode* > TIsSharpAndMedium; + typedef NCollection_DataMap< SMESH_TLink, TIsSharpAndMedium, SMESH_TLinkHasher > TLinkSharpMap; + + TLinkSharpMap linkIsSharp; + Standard_Integer nbBuckets = FromSmIdType( theMesh->NbFaces() ); + if ( nbBuckets > 0 ) + linkIsSharp.ReSize( nbBuckets ); - TLinkSharpMap linkIsSharp( theMesh->NbFaces() ); TIsSharpAndMedium sharpMedium( true, 0 ); bool & isSharp = sharpMedium.first; const SMDS_MeshNode* & nMedium = sharpMedium.second; @@ -2046,9 +2101,12 @@ SMESH_MeshAlgos::SeparateFacesByEdges( SMDS_Mesh* theMesh, const std::vector< Ed // build map of face edges (SMESH_TLink) and their faces - typedef std::vector< const SMDS_MeshElement* > TFaceVec; - typedef NCollection_DataMap< SMESH_TLink, TFaceVec, SMESH_TLink > TFacesByLinks; - TFacesByLinks facesByLink( theMesh->NbFaces() ); + typedef std::vector< const SMDS_MeshElement* > TFaceVec; + typedef NCollection_DataMap< SMESH_TLink, TFaceVec, SMESH_TLinkHasher > TFacesByLinks; + TFacesByLinks facesByLink; + Standard_Integer nbBuckets = FromSmIdType( theMesh->NbFaces() ); + if ( nbBuckets > 0 ) + facesByLink.ReSize( nbBuckets ); std::vector< const SMDS_MeshNode* > faceNodes; for ( SMDS_FaceIteratorPtr faceIt = theMesh->facesIterator(); faceIt->more(); ) @@ -2171,10 +2229,26 @@ bool SMESH_MeshAlgos::FaceNormal(const SMDS_MeshElement* F, gp_XYZ& normal, bool return ok; } -//======================================================================= -//function : GetCommonNodes -//purpose : Return nodes common to two elements -//======================================================================= +//================================================================================ +/*! + * \brief Return nodes common to two elements + */ +//================================================================================ + +int SMESH_MeshAlgos::NbCommonNodes(const SMDS_MeshElement* e1, + const SMDS_MeshElement* e2) +{ + int nb = 0; + for ( int i = 0 ; i < e1->NbNodes(); ++i ) + nb += ( e2->GetNodeIndex( e1->GetNode( i )) >= 0 ); + return nb; +} + +//================================================================================ +/*! + * \brief Return nodes common to two elements + */ +//================================================================================ std::vector< const SMDS_MeshNode*> SMESH_MeshAlgos::GetCommonNodes(const SMDS_MeshElement* e1, const SMDS_MeshElement* e2) @@ -2185,6 +2259,57 @@ std::vector< const SMDS_MeshNode*> SMESH_MeshAlgos::GetCommonNodes(const SMDS_Me common.push_back( e1->GetNode( i )); return common; } + +//================================================================================ +/*! + * \brief Return true if a node is on a boundary of 2D mesh. + * Optionally returns two neighboring boundary nodes (or more in non-manifold mesh) + */ +//================================================================================ + +bool SMESH_MeshAlgos::IsOn2DBoundary( const SMDS_MeshNode* theNode, + std::vector< const SMDS_MeshNode*> * theNeibors ) +{ + typedef NCollection_DataMap< SMESH_TLink, int, SMESH_TLinkHasher > TLinkCountMap; + TLinkCountMap linkCountMap( 10 ); + + int nbFreeLinks = 0; + for ( SMDS_ElemIteratorPtr fIt = theNode->GetInverseElementIterator(SMDSAbs_Face); fIt->more(); ) + { + const SMDS_MeshElement* face = fIt->next(); + const int nbCorners = face->NbCornerNodes(); + + int iN = face->GetNodeIndex( theNode ); + int iPrev = ( iN - 1 + nbCorners ) % nbCorners; + int iNext = ( iN + 1 ) % nbCorners; + + for ( int i : { iPrev, iNext } ) + { + SMESH_TLink link( theNode, face->GetNode( i )); + int* count = linkCountMap.ChangeSeek( link ); + if ( count ) ++( *count ); + else linkCountMap.Bind( link, 1 ); + + if ( !count ) ++nbFreeLinks; + else --nbFreeLinks; + } + } + + if ( theNeibors ) + { + theNeibors->clear(); + theNeibors->reserve( nbFreeLinks ); + for ( TLinkCountMap::Iterator linkIt( linkCountMap ); linkIt.More(); linkIt.Next() ) + if ( linkIt.Value() == 1 ) + { + theNeibors->push_back( linkIt.Key().node1() ); + if ( theNeibors->back() == theNode ) + theNeibors->back() = linkIt.Key().node2(); + } + } + return nbFreeLinks > 0; +} + //================================================================================ /*! * \brief Return true if node1 encounters first in the face and node2, after @@ -2449,3 +2574,111 @@ SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh& { return new SMESH_ElementSearcherImpl( mesh, tolerance, elemIt ); } + + +//================================================================================ +/*! + * \brief Intersect a ray with a convex volume + * \param [in] ray - the ray + * \param [in] rayLen - ray length + * \param [in] vol - the volume + * \param [out] tMin - return a ray parameter where the ray enters the volume + * \param [out] tMax - return a ray parameter where the ray exit the volume + * \param [out] iFacetMin - facet index where the ray enters the volume + * \param [out] iFacetMax - facet index where the ray exit the volume + * \return bool - true if the ray intersects the volume + */ +//================================================================================ + +bool SMESH_MeshAlgos::IntersectRayVolume( const gp_Ax1& ray, + const double rayLen, + const SMDS_MeshElement* vol, + double & tMin, + double & tMax, + int & iFacetMin, + int & iFacetMax) +{ + /* Ray-Convex Polyhedron Intersection Test by Eric Haines, erich@eye.com + * + * This test checks the ray against each face of a polyhedron, checking whether + * the set of intersection points found for each ray-plane intersection + * overlaps the previous intersection results. If there is no overlap (i.e. + * no line segment along the ray that is inside the polyhedron), then the + * ray misses and returns false; else true. + */ + SMDS_VolumeTool vTool; + if ( !vTool.Set( vol )) + return false; + + tMin = -Precision::Infinite() ; + tMax = rayLen ; + + /* Test each plane in polyhedron */ + for ( int iF = 0; iF < vTool.NbFaces(); ++iF ) + { + const SMDS_MeshNode** fNodes = vTool.GetFaceNodes( iF ); + gp_XYZ normal; + vTool.GetFaceNormal( iF, + normal.ChangeCoord(1), + normal.ChangeCoord(2), + normal.ChangeCoord(3) ); + double D = - ( normal * SMESH_NodeXYZ( fNodes[0] )); + + /* Compute intersection point T and sidedness */ + double vd = ray.Direction().XYZ() * normal; + double vn = ray.Location().XYZ() * normal + D; + if ( vd == 0.0 ) { + /* ray is parallel to plane - check if ray origin is inside plane's + half-space */ + if ( vn > 0.0 ) + /* ray origin is outside half-space */ + return false; + } + else + { + /* ray not parallel - get distance to plane */ + double t = -vn / vd ; + if ( vd < 0.0 ) + { + /* front face - T is a near point */ + if ( t > tMax ) return false; + if ( t > tMin ) { + /* hit near face */ + tMin = t ; + iFacetMin = iF; + } + } + else + { + /* back face - T is a far point */ + if ( t < tMin ) return false; + if ( t < tMax ) { + /* hit far face */ + tMax = t ; + iFacetMax = iF; + } + } + } + } + + /* survived all tests */ + /* Note: if ray originates on polyhedron, may want to change 0.0 to some + * epsilon to avoid intersecting the originating face. + */ + if ( tMin >= 0.0 ) { + /* outside, hitting front face */ + return true; + } + else + { + if ( tMax < rayLen ) { + /* inside, hitting back face */ + return true; + } + else + { + /* inside, but back face beyond tmax */ + return false; + } + } +}