X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FSMESHUtils%2FSMESH_MeshAlgos.cxx;h=e8b1187986c2b9d48745e3cd29fe0fa896521c6c;hb=274fd4f2db8d3a7fa23701764280ea1d175b194b;hp=26b3974dc52784ffe2cd98078def12191c27245f;hpb=6d32f944a0a115b6419184c50b57bf7c4eef5786;p=modules%2Fsmesh.git diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 26b3974dc..e8b118798 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-2023 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 { @@ -1254,18 +1288,43 @@ gp_XYZ SMESH_ElementSearcherImpl::Project(const gp_Pnt& point, gp_XYZ p = point.XYZ(); ElementBndBoxTree* ebbLeaf = ebbTree->getLeafAtPoint( p ); const Bnd_B3d* box = ebbLeaf ? ebbLeaf->getBox() : ebbTree->getBox(); - double radius = ( box->CornerMax() - box->CornerMin() ).Modulus(); + gp_XYZ pMin = box->CornerMin(), pMax = box->CornerMax(); + double radius = Precision::Infinite(); + if ( ebbLeaf || !box->IsOut( p )) + { + for ( int i = 1; i <= 3; ++i ) + { + double d = 0.5 * ( pMax.Coord(i) - pMin.Coord(i) ); + if ( d > Precision::Confusion() ) + radius = Min( d, radius ); + } + if ( !ebbLeaf ) + radius /= ebbTree->getHeight( /*full=*/true ); + } + else // p outside of box + { + for ( int i = 1; i <= 3; ++i ) + { + double d = 0; + if ( point.Coord(i) < pMin.Coord(i) ) + d = pMin.Coord(i) - point.Coord(i); + else if ( point.Coord(i) > pMax.Coord(i) ) + d = point.Coord(i) - pMax.Coord(i); + if ( d > Precision::Confusion() ) + radius = Min( d, radius ); + } + } ElementBndBoxTree::TElemSeq elems; ebbTree->getElementsInSphere( p, radius, elems ); while ( elems.empty() && radius < 1e100 ) { - radius *= 1.5; + radius *= 1.1; ebbTree->getElementsInSphere( p, radius, elems ); } gp_XYZ proj, bestProj; const SMDS_MeshElement* elem = 0; - double minDist = 2 * radius; + double minDist = Precision::Infinite(); ElementBndBoxTree::TElemSeq::iterator e = elems.begin(); for ( ; e != elems.end(); ++e ) { @@ -1277,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; @@ -1609,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 ); @@ -1916,7 +1992,11 @@ SMESH_MeshAlgos::FindSharpEdges( SMDS_Mesh* theMesh, typedef std::pair< bool, const SMDS_MeshNode* > TIsSharpAndMedium; typedef NCollection_DataMap< SMESH_TLink, TIsSharpAndMedium, SMESH_TLink > TLinkSharpMap; - TLinkSharpMap linkIsSharp( theMesh->NbFaces() ); + TLinkSharpMap linkIsSharp; + Standard_Integer nbBuckets = FromSmIdType( theMesh->NbFaces() ); + if ( nbBuckets > 0 ) + linkIsSharp.ReSize( nbBuckets ); + TIsSharpAndMedium sharpMedium( true, 0 ); bool & isSharp = sharpMedium.first; const SMDS_MeshNode* & nMedium = sharpMedium.second; @@ -2023,7 +2103,10 @@ SMESH_MeshAlgos::SeparateFacesByEdges( SMDS_Mesh* theMesh, const std::vector< Ed typedef std::vector< const SMDS_MeshElement* > TFaceVec; typedef NCollection_DataMap< SMESH_TLink, TFaceVec, SMESH_TLink > TFacesByLinks; - TFacesByLinks facesByLink( theMesh->NbFaces() ); + 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(); ) @@ -2146,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) @@ -2160,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_TLink > 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 @@ -2347,28 +2497,16 @@ void SMESH_MeshAlgos::Get1DBranches( SMDS_ElemIteratorPtr theEdgeIt, if ( nbBranches == 2 && !startIsBranchEnd ) // join two branches starting at the same node { - if ( nodeBranches[0].back() == nodeBranches[1].back() ) - { - // it is a closed branch, keep theStartNode first - nodeBranches[0].pop_back(); - nodeBranches[0].reserve( nodeBranches[0].size() + nodeBranches[1].size() ); - nodeBranches[0].insert( nodeBranches[0].end(), - nodeBranches[1].rbegin(), nodeBranches[1].rend() ); - branches[0].reserve( branches[0].size() + branches[1].size() ); - branches[0].insert( branches[0].end(), branches[1].rbegin(), branches[1].rend() ); - } - else - { - std::reverse( nodeBranches[0].begin(), nodeBranches[0].end() ); - nodeBranches[0].pop_back(); - nodeBranches[0].reserve( nodeBranches[0].size() + nodeBranches[1].size() ); - nodeBranches[0].insert( nodeBranches[0].end(), - nodeBranches[1].begin(), nodeBranches[1].end() ); - - std::reverse( branches[0].begin(), branches[0].end() ); - branches[0].reserve( branches[0].size() + branches[1].size() ); - branches[0].insert( branches[0].end(), branches[1].begin(), branches[1].end() ); - } + std::reverse( nodeBranches[0].begin(), nodeBranches[0].end() ); + nodeBranches[0].pop_back(); + nodeBranches[0].reserve( nodeBranches[0].size() + nodeBranches[1].size() ); + nodeBranches[0].insert( nodeBranches[0].end(), + nodeBranches[1].begin(), nodeBranches[1].end() ); + + std::reverse( branches[0].begin(), branches[0].end() ); + branches[0].reserve( branches[0].size() + branches[1].size() ); + branches[0].insert( branches[0].end(), branches[1].begin(), branches[1].end() ); + nodeBranches[1].clear(); branches[1].clear(); } @@ -2436,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; + } + } +}