X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FSMESHUtils%2FSMESH_MeshAlgos.cxx;h=e8b1187986c2b9d48745e3cd29fe0fa896521c6c;hb=274fd4f2db8d3a7fa23701764280ea1d175b194b;hp=019f57762e2202d03910ca782171f1cb672dad20;hpb=05bdaa6d2e34f8faf44db72a387557bfdfc52486;p=modules%2Fsmesh.git diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 019f57762..e8b118798 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 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 @@ -23,7 +23,7 @@ // Created : Tue Apr 30 18:00:36 2013 // Author : Edward AGAPOV (eap) -// This file holds some low level algorithms extracted from SMESH_MeshEditor +// Initially this file held some low level algorithms extracted from SMESH_MeshEditor // to make them accessible from Controls package #include "SMESH_MeshAlgos.hxx" @@ -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,11 +282,14 @@ 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; + if (SALOME::VerbosityActivated() && theElemIt && !theElemIt->more() ) + std::cout << "WARNING: ElementBndBoxTree constructed on empty iterator!" << std::endl; + SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType ); while ( elemIt->more() ) { @@ -311,7 +316,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() //================================================================================ /*! - * \brief Redistrubute element boxes among children + * \brief Redistribute element boxes among children */ //================================================================================ @@ -336,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; + // } } } @@ -461,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 @@ -519,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, @@ -660,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; @@ -843,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 */ //======================================================================= @@ -874,7 +913,7 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point, radius = point.Distance( boxCenter ) - 0.5 * ebbTree->maxSize(); if ( radius < 0 ) radius = ebbTree->maxSize() / pow( 2., getTreeHeight()) / 2; - while ( suspectElems.empty() ) + while ( suspectElems.empty() && radius < 1e100 ) { ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems ); radius *= 1.1; @@ -1071,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 { @@ -1249,22 +1288,47 @@ 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() ) + 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 ) { - double d = SMESH_MeshAlgos::GetDistance( *e, p, &proj ); + double d = SMESH_MeshAlgos::GetDistance( *e, point, &proj ); if ( d < minDist ) { bestProj = proj; @@ -1272,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; @@ -1460,7 +1541,8 @@ namespace // . RIGHT . // . . enum PositionName { POS_LEFT = 1, POS_VERTEX = 2, POS_RIGHT = 4, //POS_ON = 8, - POS_ALL = POS_LEFT | POS_RIGHT | POS_VERTEX }; + POS_ALL = POS_LEFT | POS_RIGHT | POS_VERTEX, + POS_MAX = POS_RIGHT }; struct PointPos { PositionName _name; @@ -1558,10 +1640,35 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face, const double badDistance = -1; if ( !face ) return badDistance; + int nbCorners = face->NbCornerNodes(); + if ( nbCorners > 3 ) + { + std::vector< const SMDS_MeshNode* > nodes; + int nbTria = SMESH_MeshAlgos::Triangulate().GetTriangles( face, nodes ); + + double minDist = Precision::Infinite(); + gp_XYZ cp; + for ( int i = 0; i < 3 * nbTria; i += 3 ) + { + SMDS_FaceOfNodes triangle( nodes[i], nodes[i+1], nodes[i+2] ); + double dist = GetDistance( &triangle, point, closestPnt ); + if ( dist < minDist ) + { + minDist = dist; + if ( closestPnt ) + cp = *closestPnt; + } + } + + if ( closestPnt ) + *closestPnt = cp; + return minDist; + } + // coordinates of nodes (medium nodes, if any, ignored) typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator; std::vector xyz( TXyzIterator( face->nodesIterator()), TXyzIterator() ); - xyz.resize( face->NbCornerNodes()+1 ); + xyz.resize( 4 ); // transformation to get xyz[0] lies on the origin, xyz[1] lies on the Z axis, // and xyz[2] lies in the XZ plane. This is to pass to 2D space on XZ plane. @@ -1578,14 +1685,14 @@ 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 ); // move all the nodes to 2D std::vector xy( xyz.size() ); - for ( size_t i = 0;i < xyz.size()-1; ++i ) + for ( size_t i = 0; i < 3; ++i ) { gp_XYZ p3d = xyz[i]; trsf.Transforms( p3d ); @@ -1600,71 +1707,63 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face, gp_XY point2D( tmpPnt.X(), tmpPnt.Z() ); // loop on edges of the face to analyze point position ralative to the face - std::set< PointPos > pntPosSet; + std::vector< PointPos > pntPosByType[ POS_MAX + 1 ]; for ( size_t i = 1; i < xy.size(); ++i ) { PointPos pos = getPointPosition( point2D, &xy[0], i-1 ); - pntPosSet.insert( pos ); + pntPosByType[ pos._name ].push_back( pos ); } // compute distance - double minDist2 = Precision::Infinite(); - for ( std::set< PointPos >::iterator posIt = pntPosSet.begin(); posIt != pntPosSet.end(); ++posIt) + double dist = badDistance; + + if ( pntPosByType[ POS_LEFT ].size() > 0 ) // point is most close to an edge { - PointPos pos = *posIt; - if ( pos._name != pntPosSet.begin()->_name ) - break; - switch ( pos._name ) - { - case POS_LEFT: // point is most close to an edge - { - gp_Vec edge( xyz[ pos._index ], xyz[ pos._index+1 ]); - gp_Vec n1p ( xyz[ pos._index ], point ); - double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge - // projection of the point on the edge - gp_XYZ proj = xyz[ pos._index ] + u * edge.XYZ(); - double dist2 = point.SquareDistance( proj ); - if ( dist2 < minDist2 ) - { - if ( closestPnt ) *closestPnt = proj; - minDist2 = dist2; - } - break; - } + PointPos& pos = pntPosByType[ POS_LEFT ][0]; + + gp_Vec edge( xyz[ pos._index ], xyz[ pos._index+1 ]); + gp_Vec n1p ( xyz[ pos._index ], point ); + double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge + gp_XYZ proj = xyz[ pos._index ] + u * edge.XYZ(); // projection on the edge + dist = point.Distance( proj ); + if ( closestPnt ) *closestPnt = proj; + } - case POS_RIGHT: // point is inside the face + else if ( pntPosByType[ POS_RIGHT ].size() >= 2 ) // point is inside the face + { + dist = Abs( tmpPnt.Y() ); + if ( closestPnt ) { - double distToFacePlane = Abs( tmpPnt.Y() ); - if ( closestPnt ) - { - if ( distToFacePlane < std::numeric_limits::min() ) { - *closestPnt = point.XYZ(); - } - else { - tmpPnt.SetY( 0 ); - trsf.Inverted().Transforms( tmpPnt ); - *closestPnt = tmpPnt; - } + if ( dist < std::numeric_limits::min() ) { + *closestPnt = point.XYZ(); + } + else { + tmpPnt.SetY( 0 ); + trsf.Inverted().Transforms( tmpPnt ); + *closestPnt = tmpPnt; } - return distToFacePlane; } + } - case POS_VERTEX: // point is most close to a node + else if ( pntPosByType[ POS_VERTEX ].size() > 0 ) // point is most close to a node + { + double minDist2 = Precision::Infinite(); + for ( size_t i = 0; i < pntPosByType[ POS_VERTEX ].size(); ++i ) { - double dist2 = point.SquareDistance( xyz[ pos._index ]); - if ( dist2 < minDist2 ) + PointPos& pos = pntPosByType[ POS_VERTEX ][i]; + + double d2 = point.SquareDistance( xyz[ pos._index ]); + if ( minDist2 > d2 ) { if ( closestPnt ) *closestPnt = xyz[ pos._index ]; - minDist2 = dist2; + minDist2 = d2; } - break; - } - default:; - return badDistance; } + dist = Sqrt( minDist2 ); } - return Sqrt( minDist2 ); + + return dist; } //======================================================================= @@ -1741,7 +1840,7 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, !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 ) + if ( gp_XYZ( n[0], n[1], n[2] ) * bcp < -1e-12 ) continue; // find distance to a facet @@ -1893,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; @@ -2000,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(); ) @@ -2123,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) @@ -2137,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 @@ -2166,6 +2339,195 @@ bool SMESH_MeshAlgos::IsRightOrder( const SMDS_MeshElement* face, return ( diff == 1 ) || ( diff == -face->NbNodes()+1 ); } +//======================================================================= +/*! + * \brief Partition given 1D elements into groups of contiguous edges. + * A node where number of meeting edges != 2 is a group end. + * An optional startNode is used to orient groups it belongs to. + * \return a list of edge groups and a list of corresponding node groups. + * If a group is closed, the first and last nodes of the group are same. + */ +//======================================================================= + +void SMESH_MeshAlgos::Get1DBranches( SMDS_ElemIteratorPtr theEdgeIt, + TElemGroupVector& theEdgeGroups, + TNodeGroupVector& theNodeGroups, + const SMDS_MeshNode* theStartNode ) +{ + if ( !theEdgeIt ) + return; + + // build map of nodes and their adjacent edges + + typedef std::vector< const SMDS_MeshNode* > TNodeVec; + typedef std::vector< const SMDS_MeshElement* > TEdgeVec; + typedef NCollection_DataMap< const SMDS_MeshNode*, TEdgeVec, SMESH_Hasher > TEdgesByNodeMap; + TEdgesByNodeMap edgesByNode; + + while ( theEdgeIt->more() ) + { + const SMDS_MeshElement* edge = theEdgeIt->next(); + if ( edge->GetType() != SMDSAbs_Edge ) + continue; + + const SMDS_MeshNode* nodes[2] = { edge->GetNode(0), edge->GetNode(1) }; + for ( int i = 0; i < 2; ++i ) + { + TEdgeVec* nodeEdges = edgesByNode.ChangeSeek( nodes[i] ); + if ( !nodeEdges ) + { + nodeEdges = edgesByNode.Bound( nodes[i], TEdgeVec() ); + nodeEdges->reserve(2); + } + nodeEdges->push_back( edge ); + } + } + + if ( edgesByNode.IsEmpty() ) + return; + + + // build edge branches + + TElemGroupVector branches(2); + TNodeGroupVector nodeBranches(2); + + while ( !edgesByNode.IsEmpty() ) + { + if ( !theStartNode || !edgesByNode.IsBound( theStartNode )) + { + theStartNode = TEdgesByNodeMap::Iterator( edgesByNode ).Key(); + } + + size_t nbBranches = 0; + bool startIsBranchEnd = false; + + while ( edgesByNode.IsBound( theStartNode )) + { + // initialize a new branch + + ++nbBranches; + if ( branches.size() < nbBranches ) + { + branches.push_back ( TEdgeVec() ); + nodeBranches.push_back( TNodeVec() ); + } + TEdgeVec & branch = branches [ nbBranches - 1 ]; + TNodeVec & nodeBranch = nodeBranches[ nbBranches - 1 ]; + branch.clear(); + nodeBranch.clear(); + { + TEdgeVec& edges = edgesByNode( theStartNode ); + startIsBranchEnd = ( edges.size() != 2 ); + + int nbEdges = 0; + const SMDS_MeshElement* startEdge = 0; + for ( size_t i = 0; i < edges.size(); ++i ) + { + if ( !startEdge && edges[i] ) + { + startEdge = edges[i]; + edges[i] = 0; + } + nbEdges += bool( edges[i] ); + } + if ( nbEdges == 0 ) + edgesByNode.UnBind( theStartNode ); + if ( !startEdge ) + continue; + + branch.push_back( startEdge ); + + nodeBranch.push_back( theStartNode ); + nodeBranch.push_back( branch.back()->GetNode(0) ); + if ( nodeBranch.back() == theStartNode ) + nodeBranch.back() = branch.back()->GetNode(1); + } + + // fill the branch + + bool isBranchEnd = false; + TEdgeVec* edgesPtr; + + while (( !isBranchEnd ) && ( edgesPtr = edgesByNode.ChangeSeek( nodeBranch.back() ))) + { + TEdgeVec& edges = *edgesPtr; + + isBranchEnd = ( edges.size() != 2 ); + + const SMDS_MeshNode* lastNode = nodeBranch.back(); + + switch ( edges.size() ) + { + case 1: + edgesByNode.UnBind( lastNode ); + break; + + case 2: + { + if ( const SMDS_MeshElement* nextEdge = edges[ edges[0] == branch.back() ]) + { + branch.push_back( nextEdge ); + + const SMDS_MeshNode* nextNode = nextEdge->GetNode(0); + if ( nodeBranch.back() == nextNode ) + nextNode = nextEdge->GetNode(1); + nodeBranch.push_back( nextNode ); + } + edgesByNode.UnBind( lastNode ); + break; + } + + default: + int nbEdges = 0; + for ( size_t i = 0; i < edges.size(); ++i ) + { + if ( edges[i] == branch.back() ) + edges[i] = 0; + nbEdges += bool( edges[i] ); + } + if ( nbEdges == 0 ) + edgesByNode.UnBind( lastNode ); + } + } + } // while ( edgesByNode.IsBound( theStartNode )) + + + // put the found branches to the result + + if ( nbBranches == 2 && !startIsBranchEnd ) // join two branches starting at the same node + { + 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(); + } + + for ( size_t i = 0; i < nbBranches; ++i ) + { + if ( branches[i].empty() ) + continue; + + theEdgeGroups.push_back( TEdgeVec() ); + theEdgeGroups.back().swap( branches[i] ); + + theNodeGroups.push_back( TNodeVec() ); + theNodeGroups.back().swap( nodeBranches[i] ); + } + + } // while ( !edgesByNode.IsEmpty() ) + + return; +} + //======================================================================= /*! * \brief Return SMESH_NodeSearcher @@ -2212,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; + } + } +}