X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FSMESHUtils%2FSMESH_MeshAlgos.cxx;h=e8b1187986c2b9d48745e3cd29fe0fa896521c6c;hb=274fd4f2db8d3a7fa23701764280ea1d175b194b;hp=9f0c88e20e9c849f9e1efe4dc7c1c24ea94853f7;hpb=b09372829929f8f561495d6c16527134971a1909;p=modules%2Fsmesh.git diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 9f0c88e20..e8b118798 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2021 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 @@ -253,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; @@ -286,10 +287,8 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() 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() ) @@ -342,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; + // } } } @@ -2252,6 +2260,56 @@ std::vector< const SMDS_MeshNode*> SMESH_MeshAlgos::GetCommonNodes(const SMDS_Me 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 @@ -2516,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; + } + } +}