X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESHUtils%2FSMESH_Offset.cxx;h=71f53213ded3a5c564c73528fa19379518f835bc;hp=f571535f2dcf73fd729f2a2e06d4850715fab4e5;hb=155c6427afd2f7517d00bb5e9cb56aa2da867ddf;hpb=ead27b045b448d98442ad5817fd47947178cd7e8 diff --git a/src/SMESHUtils/SMESH_Offset.cxx b/src/SMESHUtils/SMESH_Offset.cxx index f571535f2..71f53213d 100644 --- a/src/SMESHUtils/SMESH_Offset.cxx +++ b/src/SMESHUtils/SMESH_Offset.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2019 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 @@ -40,20 +40,20 @@ namespace { - const size_t theMaxNbFaces = 256; // max number of faces sharing a node + const int theMaxNbFaces = 256; // max number of faces sharing a node - typedef NCollection_DataMap< Standard_Address, const SMDS_MeshNode* > TNNMap; - typedef NCollection_Map< SMESH_Link, SMESH_Link > TLinkMap; + typedef NCollection_DataMap< const SMDS_MeshNode*, const SMDS_MeshNode*, SMESH_Hasher > TNNMap; + typedef NCollection_Map< SMESH_Link, SMESH_Link > TLinkMap; //-------------------------------------------------------------------------------- /*! * \brief Intersected face side storing a node created at this intersection - * and a intersected face + * and an intersected face */ struct CutLink { bool myReverse; - const SMDS_MeshNode* myNode[2]; // side nodes + const SMDS_MeshNode* myNode[2]; // side nodes. WARNING: don't set them directly, use Set() mutable SMESH_NodeXYZ myIntNode; // intersection node const SMDS_MeshElement* myFace; // cutter face int myIndex; // index of a node on the same link @@ -105,7 +105,7 @@ namespace int myIndex; // positive -> side index, negative -> State const SMDS_MeshElement* myFace; - enum State { _INTERNAL = -1, _COPLANAR = -2 }; + enum State { _INTERNAL = -1, _COPLANAR = -2, _PENDING = -3 }; void Set( const SMDS_MeshNode* Node1, const SMDS_MeshNode* Node2, @@ -141,11 +141,11 @@ namespace EdgeLoop() : SMDS_PolygonalFaceOfNodes( std::vector() ) {} void Clear() { myLinks.clear(); myIsBndConnected = false; myHasPending = false; } bool SetConnected() { bool was = myIsBndConnected; myIsBndConnected = true; return !was; } - bool Contains( const SMDS_MeshNode* n ) const + size_t Contains( const SMDS_MeshNode* n ) const { for ( size_t i = 0; i < myLinks.size(); ++i ) - if ( myLinks[i]->myNode1 == n ) return true; - return false; + if ( myLinks[i]->myNode1 == n ) return i + 1; + return 0; } virtual int NbNodes() const { return myLinks.size(); } virtual SMDS_ElemIteratorPtr nodesIterator() const @@ -224,6 +224,24 @@ namespace myLoopOfEdge[ Index( *loop->myLinks[ iE ] )] = 0; loop->Clear(); } + void Join( EdgeLoop& loop1, size_t iAfterConcact, + EdgeLoop& loop2, size_t iFromEdge2 ) + { + std::vector< const EdgePart* > linksAfterContact( loop1.myLinks.begin() + iAfterConcact, + loop1.myLinks.end() ); + loop1.myLinks.reserve( loop2.myLinks.size() + loop1.myLinks.size() ); + loop1.myLinks.resize( iAfterConcact ); + loop1.myLinks.insert( loop1.myLinks.end(), + loop2.myLinks.begin() + iFromEdge2, loop2.myLinks.end() ); + loop1.myLinks.insert( loop1.myLinks.end(), + loop2.myLinks.begin(), loop2.myLinks.begin() + iFromEdge2 ); + loop1.myLinks.insert( loop1.myLinks.end(), + linksAfterContact.begin(), linksAfterContact.end() ); + loop1.myIsBndConnected = loop2.myIsBndConnected; + loop2.Clear(); + for ( size_t iE = 0; iE < loop1.myLinks.size(); ++iE ) + myLoopOfEdge[ Index( *loop1.myLinks[ iE ] )] = & loop1; + } size_t Index( const EdgePart& edge ) const { return &edge - myEdge0; } EdgeLoop* GetLoopOf( const EdgePart* edge ) { return myLoopOfEdge[ Index( *edge )]; } }; @@ -575,7 +593,7 @@ namespace TIDSortedElemSet elemSet, avoidSet; int iFace = 0; const SMDS_MeshElement* f; - for ( ; faceIt->more(); faceIt->next() ) + for ( ; faceIt->more() && iFace < theMaxNbFaces; faceIt->next() ) { avoidSet.insert( faces[ iFace ].myFace ); f = SMESH_MeshAlgos::FindFaceInSet( theNewNode, faces[ iFace ].myNode2.Node(), @@ -597,7 +615,7 @@ namespace faces[ iFace ].SetNodes( i0, i1 ); faces[ iFace ].SetNormal( theFaceNormals ); } - int nbFaces = Min( iFace + 1, (int)theMaxNbFaces ); + int nbFaces = iFace + 1; theNewPos.SetCoord( 0, 0, 0 ); gp_XYZ oldXYZ = SMESH_NodeXYZ( theNewNode ); @@ -684,7 +702,7 @@ namespace { gp_XYZ p1 = oldXYZ + faces[ i ].Norm() * theOffset; gp_XYZ p2 = oldXYZ + faces[ iPrev ].Norm() * theOffset; - useOneNormal = ( p1 - p2 ).SquareModulus() > theTol * theTol; + useOneNormal = ( p1 - p2 ).SquareModulus() > 1e-12; } } if ( useOneNormal && theNewNode->isMarked() ) @@ -694,11 +712,15 @@ namespace return useOneNormal; } +} // namespace + +namespace SMESH_MeshAlgos +{ //-------------------------------------------------------------------------------- /*! * \brief Intersect faces of a mesh */ - struct Intersector + struct Intersector::Algo { SMDS_Mesh* myMesh; double myTol, myEps; @@ -716,7 +738,7 @@ namespace int myNbOnPlane1, myNbOnPlane2; TIntPointSet myIntPointSet; - Intersector( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals ) + Algo( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals ) : myMesh( mesh ), myTol( tol ), myEps( 1e-100 ), @@ -727,9 +749,17 @@ namespace void Cut( const SMDS_MeshElement* face1, const SMDS_MeshElement* face2, const int nbCommonNodes ); - void MakeNewFaces( SMESH_MeshAlgos::TEPairVec& theNew2OldFaces, - SMESH_MeshAlgos::TNPairVec& theNew2OldNodes, - const double theSign ); + void Cut( const SMDS_MeshElement* face, + SMESH_NodeXYZ& lineEnd1, + int edgeIndex1, + SMESH_NodeXYZ& lineEnd2, + int edgeIndex2 ); + void MakeNewFaces( TElemIntPairVec& theNew2OldFaces, + TNodeIntPairVec& theNew2OldNodes, + const double theSign, + const bool theOptimize ); + + void IntersectNewEdges( const CutFace& theCFace ); private: @@ -782,7 +812,6 @@ namespace bool & isCollinear ); bool intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint ); bool isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes ); - void intersectNewEdges( const CutFace& theCFace ); const SMDS_MeshNode* createNode( const gp_XYZ& p ); }; @@ -805,7 +834,7 @@ namespace */ //================================================================================ - const SMDS_MeshNode* Intersector::createNode( const gp_XYZ& p ) + const SMDS_MeshNode* Intersector::Algo::createNode( const gp_XYZ& p ) { const SMDS_MeshNode* n = myMesh->AddNode( p.X(), p.Y(), p.Z() ); n->setIsMarked( true ); // cut nodes are marked @@ -818,7 +847,7 @@ namespace */ //================================================================================ - void Intersector::addLink( CutLink& link ) + void Intersector::Algo::addLink( CutLink& link ) { link.myIndex = 0; const CutLink* added = & myCutLinks.Added( link ); @@ -844,7 +873,7 @@ namespace */ //================================================================================ - bool Intersector::findLink( CutLink& link ) + bool Intersector::Algo::findLink( CutLink& link ) { link.myIndex = 0; while ( myCutLinks.Contains( link )) @@ -872,12 +901,12 @@ namespace */ //================================================================================ - bool Intersector::isPlaneIntersected( const gp_XYZ& n2, - const double d2, - const std::vector< SMESH_NodeXYZ >& nodes1, - std::vector< double > & dist1, - int & nbOnPlane1, - int & iNotOnPlane1) + bool Intersector::Algo::isPlaneIntersected( const gp_XYZ& n2, + const double d2, + const std::vector< SMESH_NodeXYZ >& nodes1, + std::vector< double > & dist1, + int & nbOnPlane1, + int & iNotOnPlane1) { iNotOnPlane1 = nbOnPlane1 = 0; dist1.resize( nodes1.size() ); @@ -915,12 +944,12 @@ namespace */ //================================================================================ - void Intersector::computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes, - const std::vector< double >& dist, - const int nbOnPln, - const int iMaxCoo, - double * u, - int* iE) + void Intersector::Algo::computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes, + const std::vector< double >& dist, + const int nbOnPln, + const int iMaxCoo, + double * u, + int* iE) { if ( nbOnPln == 3 ) { @@ -961,9 +990,9 @@ namespace */ //================================================================================ - void Intersector::findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes, - const std::vector< double > & dist, - CutLink& link ) + void Intersector::Algo::findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes, + const std::vector< double > & dist, + CutLink& link ) { int i1 = ( dist[0] == 0 ? 0 : 1 ), i2 = ( dist[2] == 0 ? 2 : 1 ); CutLink link2 = link; @@ -978,11 +1007,11 @@ namespace */ //================================================================================ - void Intersector::intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1, - const std::vector< double > & dist1, - const int iEdge1, - const SMDS_MeshElement* face2, - CutLink& link1) + void Intersector::Algo::intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1, + const std::vector< double > & dist1, + const int iEdge1, + const SMDS_MeshElement* face2, + CutLink& link1) { const int iEdge2 = ( iEdge1 + 1 ) % nodes1.size(); const SMESH_NodeXYZ& p1 = nodes1[ iEdge1 ]; @@ -1028,15 +1057,15 @@ namespace */ //================================================================================ - void Intersector::replaceIntNode( const SMDS_MeshNode* nToKeep, - const SMDS_MeshNode* nToRemove ) + void Intersector::Algo::replaceIntNode( const SMDS_MeshNode* nToKeep, + const SMDS_MeshNode* nToRemove ) { if ( nToKeep == nToRemove ) return; if ( nToRemove->GetID() < nToKeep->GetID() ) // keep node with lower ID - myRemove2KeepNodes.Bind((void*) nToKeep, nToRemove ); + myRemove2KeepNodes.Bind( nToKeep, nToRemove ); else - myRemove2KeepNodes.Bind((void*) nToRemove, nToKeep ); + myRemove2KeepNodes.Bind( nToRemove, nToKeep ); } //================================================================================ @@ -1053,13 +1082,13 @@ namespace */ //================================================================================ - void Intersector::computeIntPoint( const double u1, - const double u2, - const int iE1, - const int iE2, - CutLink & link, - const SMDS_MeshNode* & node1, - const SMDS_MeshNode* & node2) + void Intersector::Algo::computeIntPoint( const double u1, + const double u2, + const int iE1, + const int iE2, + CutLink & link, + const SMDS_MeshNode* & node1, + const SMDS_MeshNode* & node2) { if ( u1 > u2 + myTol ) { @@ -1100,11 +1129,11 @@ namespace */ //================================================================================ - void Intersector::cutCollinearLink( const int iNotOnPlane1, - const std::vector< SMESH_NodeXYZ >& nodes1, - const SMDS_MeshElement* face2, - const CutLink& link1, - const CutLink& link2) + void Intersector::Algo::cutCollinearLink( const int iNotOnPlane1, + const std::vector< SMESH_NodeXYZ >& nodes1, + const SMDS_MeshElement* face2, + const CutLink& link1, + const CutLink& link2) { int iN1 = ( iNotOnPlane1 + 1 ) % 3; @@ -1128,7 +1157,7 @@ namespace */ //================================================================================ - void Intersector::setPlaneIndices( const gp_XYZ& planeNorm ) + void Intersector::Algo::setPlaneIndices( const gp_XYZ& planeNorm ) { switch ( MaxIndex( planeNorm )) { case 1: myInd1 = 2; myInd2 = 3; break; @@ -1143,9 +1172,9 @@ namespace */ //================================================================================ - void Intersector::Cut( const SMDS_MeshElement* face1, - const SMDS_MeshElement* face2, - const int nbCommonNodes) + void Intersector::Algo::Cut( const SMDS_MeshElement* face1, + const SMDS_MeshElement* face2, + const int nbCommonNodes) { myFace1 = face1; myFace2 = face2; @@ -1241,16 +1270,83 @@ namespace return; } + //================================================================================ + /*! + * \brief Store a face cut by a line given by its ends + * accompanied by indices of intersected face edges. + * Edge index is <0 if a line end is inside the face. + * \param [in] face - a face to cut + * \param [inout] lineEnd1 - line end coordinates + optional node existing at this point + * \param [in] edgeIndex1 - index of face edge cut by lineEnd1 + * \param [inout] lineEnd2 - line end coordinates + optional node existing at this point + * \param [in] edgeIndex2 - index of face edge cut by lineEnd2 + */ + //================================================================================ + + void Intersector::Algo::Cut( const SMDS_MeshElement* face, + SMESH_NodeXYZ& lineEnd1, + int edgeIndex1, + SMESH_NodeXYZ& lineEnd2, + int edgeIndex2 ) + { + if ( lineEnd1.Node() && lineEnd2.Node() && + face->GetNodeIndex( lineEnd1.Node() ) >= 0 && + face->GetNodeIndex( lineEnd2.Node() ) >= 0 ) + return; // intersection at a face node or edge + + if ((int) myNormals.size() <= face->GetID() ) + const_cast< std::vector< gp_XYZ >& >( myNormals ).resize( face->GetID() + 1 ); + + const CutFace& cf = myCutFaces.Added( CutFace( face )); + cf.InitLinks(); + + // look for intersection nodes coincident with line ends + CutLink links[2]; + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + SMESH_NodeXYZ& lineEnd = is2nd ? lineEnd2 : lineEnd1; + int edgeIndex = is2nd ? edgeIndex2 : edgeIndex1; + CutLink & link = links[ is2nd ]; + + link.myIntNode = lineEnd; + + for ( size_t i = ( edgeIndex < 0 ? 3 : 0 ); i < cf.myLinks.size(); ++i ) + if ( coincide( lineEnd, SMESH_NodeXYZ( cf.myLinks[i].myNode1 ), myTol )) + { + link.myIntNode = cf.myLinks[i].myNode1; + break; + } + + if ( edgeIndex >= 0 ) + { + link.Set( face->GetNode ( edgeIndex ), + face->GetNodeWrap( edgeIndex + 1 ), + /*cuttingFace=*/0); + findLink( link ); + } + + if ( !link.myIntNode ) + link.myIntNode.Set( createNode( lineEnd )); + + lineEnd._node = link.IntNode(); + + if ( link.myNode[0] ) + addLink( link ); + } + + cf.AddEdge( links[0], links[1], /*face=*/0, /*nbOnPlane=*/0, /*iNotOnPlane=*/-1 ); + } + //================================================================================ /*! * \brief Intersect two 2D line segments */ //================================================================================ - bool Intersector::intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1, - const gp_XY s2p0, const gp_XY s2p1, - double & t1, double & t2, - bool & isCollinear ) + bool Intersector::Algo::intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1, + const gp_XY s2p0, const gp_XY s2p1, + double & t1, double & t2, + bool & isCollinear ) { gp_XY u = s1p1 - s1p0; gp_XY v = s2p1 - s2p0; @@ -1306,7 +1402,7 @@ namespace */ //================================================================================ - bool Intersector::intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint ) + bool Intersector::Algo::intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint ) { int i01 = iE1, i11 = ( iE1 + 1 ) % 3; int i02 = iE2, i12 = ( iE2 + 1 ) % 3; @@ -1393,7 +1489,8 @@ namespace for ( int is2nd = 0; is2nd < 2; ++is2nd ) { const SMDS_MeshElement* f = is2nd ? myFace1 : myFace2; - const CutFace& cf = myCutFaces.Added( CutFace( is2nd ? myFace2 : myFace1 )); + if ( !f ) continue; + const CutFace& cf = myCutFaces.Added( CutFace( is2nd ? myFace2 : myFace1 )); for ( size_t i = 0; i < cf.myLinks.size(); ++i ) if ( cf.myLinks[i].myFace == f && //cf.myLinks[i].myIndex != EdgePart::_COPLANAR && @@ -1421,13 +1518,14 @@ namespace */ //================================================================================ - bool Intersector::isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes ) + bool Intersector::Algo::isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes ) { double bc1, bc2; SMESH_MeshAlgos::GetBarycentricCoords( p2D( p ), p2D( nodes[0] ), p2D( nodes[1] ), p2D( nodes[2] ), bc1, bc2 ); - return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. ); + //return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. ); + return ( myTol < bc1 && myTol < bc2 && bc1 + bc2 + myTol < 1. ); } //================================================================================ @@ -1436,7 +1534,7 @@ namespace */ //================================================================================ - void Intersector::cutCoplanar() + void Intersector::Algo::cutCoplanar() { // find intersections of edges @@ -1553,7 +1651,7 @@ namespace } return; - } // Intersector::cutCoplanar() + } // Intersector::Algo::cutCoplanar() //================================================================================ /*! @@ -1561,19 +1659,29 @@ namespace */ //================================================================================ - void Intersector::intersectNewEdges( const CutFace& cf ) + void Intersector::Algo::IntersectNewEdges( const CutFace& cf ) { IntPoint2D intPoint; if ( cf.NbInternalEdges() < 2 ) return; + if ( myNodes1.empty() ) + { + myNodes1.resize(2); + myNodes2.resize(2); + } + const gp_XYZ& faceNorm = myNormals[ cf.myInitFace->GetID() ]; setPlaneIndices( faceNorm ); // choose indices on an axis-aligned plane size_t limit = cf.myLinks.size() * cf.myLinks.size() * 2; - for ( size_t i1 = 3; i1 < cf.myLinks.size(); ++i1 ) + size_t i1 = 3; + while ( cf.myLinks[i1-1].IsInternal() && i1 > 0 ) + --i1; + + for ( ; i1 < cf.myLinks.size(); ++i1 ) { if ( !cf.myLinks[i1].IsInternal() ) continue; @@ -1752,7 +1860,7 @@ namespace } } if ( cf.myLinks.size() >= limit ) - throw SALOME_Exception( "Infinite loop in Intersector::intersectNewEdges()" ); + throw SALOME_Exception( "Infinite loop in Intersector::Algo::IntersectNewEdges()" ); } ++i1; // each internal edge encounters twice } @@ -1765,10 +1873,23 @@ namespace */ //================================================================================ - void Intersector::MakeNewFaces( SMESH_MeshAlgos::TEPairVec& theNew2OldFaces, - SMESH_MeshAlgos::TNPairVec& theNew2OldNodes, - const double theSign) + void Intersector::Algo::MakeNewFaces( SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces, + SMESH_MeshAlgos::TNodeIntPairVec& theNew2OldNodes, + const double theSign, + const bool theOptimize) { + // fill theNew2OldFaces if empty + TCutFaceMap::const_iterator cutFacesIt = myCutFaces.cbegin(); + if ( theNew2OldFaces.empty() ) + for ( ; cutFacesIt != myCutFaces.cend(); ++cutFacesIt ) + { + const CutFace& cf = *cutFacesIt; + int index = cf.myInitFace->GetID(); // index in theNew2OldFaces + if ((int) theNew2OldFaces.size() <= index ) + theNew2OldFaces.resize( index + 1 ); + theNew2OldFaces[ index ] = std::make_pair( cf.myInitFace, index ); + } + // unmark all nodes except intersection ones for ( SMDS_NodeIteratorPtr nIt = myMesh->nodesIterator(); nIt->more(); ) @@ -1789,22 +1910,21 @@ namespace // intersect edges added to myCutFaces - TCutFaceMap::const_iterator cutFacesIt = myCutFaces.cbegin(); - for ( ; cutFacesIt != myCutFaces.cend(); ++cutFacesIt ) + for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt ) { const CutFace& cf = *cutFacesIt; cf.ReplaceNodes( myRemove2KeepNodes ); - intersectNewEdges( cf ); + IntersectNewEdges( cf ); } // make new faces EdgeLoopSet loopSet; - SMESH_MeshAlgos::Triangulate triangulator; + SMESH_MeshAlgos::Triangulate triangulator( theOptimize ); std::vector< EdgePart > cutOffLinks; TLinkMap cutOffCoplanarLinks; std::vector< const CutFace* > touchedFaces; - SMESH_MeshAlgos::TEPairVec::value_type new2OldTria; + SMESH_MeshAlgos::TElemIntPairVec::value_type new2OldTria; CutFace cutFace(0); std::vector< const SMDS_MeshNode* > nodes; std::vector faces; @@ -1829,7 +1949,7 @@ namespace // avoid loops that are not connected to boundary edges of cf.myInitFace if ( cf.RemoveInternalLoops( loopSet )) { - intersectNewEdges( cf ); + IntersectNewEdges( cf ); cf.MakeLoops( loopSet, normal ); } // erase loops that are cut off by face intersections @@ -1876,13 +1996,36 @@ namespace // remove split faces for ( size_t id = 1; id < theNew2OldFaces.size(); ++id ) { - if ( theNew2OldFaces[id].first ) + if ( theNew2OldFaces[id].first || + theNew2OldFaces[id].second == 0 ) continue; if ( const SMDS_MeshElement* f = myMesh->FindElement( id )) myMesh->RemoveFreeElement( f ); } - // remove face connected to cut off parts of cf.myInitFace + // remove faces that are merged off + for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt ) + { + const CutFace& cf = *cutFacesIt; + if ( !cf.myLinks.empty() || cf.myInitFace->IsNull() ) + continue; + + nodes.assign( cf.myInitFace->begin_nodes(), cf.myInitFace->end_nodes() ); + for ( size_t i = 0; i < nodes.size(); ++i ) + { + const SMDS_MeshNode* n = nodes[ i ]; + while ( myRemove2KeepNodes.IsBound( n )) + n = myRemove2KeepNodes( n ); + if ( n != nodes[ i ] && cf.myInitFace->GetNodeIndex( n ) >= 0 ) + { + theNew2OldFaces[ cf.myInitFace->GetID() ].first = 0; + myMesh->RemoveFreeElement( cf.myInitFace ); + break; + } + } + } + + // remove faces connected to cut off parts of cf.myInitFace nodes.resize(2); for ( size_t i = 0; i < cutOffLinks.size(); ++i ) @@ -1894,9 +2037,9 @@ namespace if ( nodes[0] != nodes[1] && myMesh->GetElementsByNodes( nodes, faces )) { - if ( cutOffLinks[i].myFace && - cutOffLinks[i].myIndex != EdgePart::_COPLANAR && - faces.size() == 2 ) + if ( // cutOffLinks[i].myFace && + cutOffLinks[i].myIndex != EdgePart::_COPLANAR && + faces.size() != 1 ) continue; for ( size_t iF = 0; iF < faces.size(); ++iF ) { @@ -1929,13 +2072,22 @@ namespace for ( size_t i = 0; i < touchedFaces.size(); ++i ) { const CutFace& cf = *touchedFaces[i]; + if ( cf.myInitFace->IsNull() ) + continue; int index = cf.myInitFace->GetID(); // index in theNew2OldFaces if ( !theNew2OldFaces[ index ].first ) continue; // already cut off + cf.InitLinks(); if ( !cf.ReplaceNodes( myRemove2KeepNodes )) - continue; // just keep as is + { + if ( cf.myLinks.size() == 3 && + cf.myInitFace->GetNodeIndex( cf.myLinks[0].myNode1 ) >= 0 && + cf.myInitFace->GetNodeIndex( cf.myLinks[1].myNode1 ) >= 0 && + cf.myInitFace->GetNodeIndex( cf.myLinks[2].myNode1 ) >= 0 ) + continue; // just keep as is + } if ( cf.myLinks.size() == 3 ) { @@ -1953,8 +2105,8 @@ namespace // add used new nodes to theNew2OldNodes - SMESH_MeshAlgos::TNPairVec::value_type new2OldNode; - new2OldNode.second = NULL; + SMESH_MeshAlgos::TNodeIntPairVec::value_type new2OldNode; + new2OldNode.second = 0; for ( cutLinksIt = myCutLinks.cbegin(); cutLinksIt != myCutLinks.cend(); ++cutLinksIt ) { const CutLink& link = *cutLinksIt; @@ -1968,6 +2120,189 @@ namespace return; } + //================================================================================ + Intersector::Intersector( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals ) + { + myAlgo = new Algo( mesh, tol, normals ); + } + //================================================================================ + Intersector::~Intersector() + { + delete myAlgo; + } + //================================================================================ + //! compute cut of two faces of the mesh + void Intersector::Cut( const SMDS_MeshElement* face1, + const SMDS_MeshElement* face2, + const int nbCommonNodes ) + { + myAlgo->Cut( face1, face2, nbCommonNodes ); + } + //================================================================================ + //! store a face cut by a line given by its ends + // accompanied by indices of intersected face edges. + // Edge index is <0 if a line end is inside the face. + void Intersector::Cut( const SMDS_MeshElement* face, + SMESH_NodeXYZ& lineEnd1, + int edgeIndex1, + SMESH_NodeXYZ& lineEnd2, + int edgeIndex2 ) + { + myAlgo->Cut( face, lineEnd1, edgeIndex1, lineEnd2, edgeIndex2 ); + } + //================================================================================ + //! split all face intersected by Cut() methods + void Intersector::MakeNewFaces( SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces, + SMESH_MeshAlgos::TNodeIntPairVec& theNew2OldNodes, + const double theSign, + const bool theOptimize ) + { + myAlgo->MakeNewFaces( theNew2OldFaces, theNew2OldNodes, theSign, theOptimize ); + } + //================================================================================ + //! Cut a face by planes, whose normals point to parts to keep + bool Intersector::CutByPlanes(const SMDS_MeshElement* theFace, + const std::vector< gp_Ax1 > & thePlanes, + const double theTol, + std::vector< TFace > & theNewFaceConnectivity ) + { + theNewFaceConnectivity.clear(); + + // check if theFace is wholly cut off + std::vector< SMESH_NodeXYZ > facePoints( theFace->begin_nodes(), theFace->end_nodes() ); + facePoints.resize( theFace->NbCornerNodes() ); + for ( size_t iP = 0; iP < thePlanes.size(); ++iP ) + { + size_t nbOut = 0; + const gp_Pnt& O = thePlanes[iP].Location(); + for ( size_t i = 0; i < facePoints.size(); ++i ) + { + gp_Vec Op( O, facePoints[i] ); + nbOut += ( Op * thePlanes[iP].Direction() <= 0 ); + } + if ( nbOut == facePoints.size() ) + return true; + } + + // copy theFace into a temporary mesh + SMDS_Mesh mesh; + Bnd_B3d faceBox; + std::vector< const SMDS_MeshNode* > faceNodes; + faceNodes.resize( facePoints.size() ); + for ( size_t i = 0; i < facePoints.size(); ++i ) + { + const SMESH_NodeXYZ& n = facePoints[i]; + faceNodes[i] = mesh.AddNode( n.X(), n.Y(), n.Z() ); + faceBox.Add( n ); + } + const SMDS_MeshElement* faceToCut = 0; + switch ( theFace->NbCornerNodes() ) + { + case 3: + faceToCut = mesh.AddFace( faceNodes[0], faceNodes[1], faceNodes[2] ); + break; + case 4: + faceToCut = mesh.AddFace( faceNodes[0], faceNodes[1], faceNodes[2], faceNodes[3] ); + break; + default: + faceToCut = mesh.AddPolygonalFace( faceNodes ); + } + + std::vector< gp_XYZ > normals( 2 + thePlanes.size() ); + SMESH_MeshAlgos::FaceNormal( faceToCut, normals[ faceToCut->GetID() ]); + + // add faces corresponding to thePlanes + std::vector< const SMDS_MeshElement* > planeFaces; + double faceSize = Sqrt( faceBox.SquareExtent() ); + gp_XYZ center = 0.5 * ( faceBox.CornerMin() + faceBox.CornerMax() ); + for ( size_t i = 0; i < thePlanes.size(); ++i ) + { + gp_Ax2 plnAx( thePlanes[i].Location(), thePlanes[i].Direction() ); + gp_XYZ O = plnAx.Location().XYZ(); + gp_XYZ X = plnAx.XDirection().XYZ(); + gp_XYZ Y = plnAx.YDirection().XYZ(); + gp_XYZ Z = plnAx.Direction().XYZ(); + + double dot = ( O - center ) * Z; + gp_XYZ o = center + Z * dot; // center projected to a plane + + gp_XYZ p1 = o + X * faceSize * 2; + gp_XYZ p2 = o + Y * faceSize * 2; + gp_XYZ p3 = o - (X + Y ) * faceSize * 2; + + const SMDS_MeshNode* n1 = mesh.AddNode( p1.X(), p1.Y(), p1.Z() ); + const SMDS_MeshNode* n2 = mesh.AddNode( p2.X(), p2.Y(), p2.Z() ); + const SMDS_MeshNode* n3 = mesh.AddNode( p3.X(), p3.Y(), p3.Z() ); + planeFaces.push_back( mesh.AddFace( n1, n2, n3 )); + + normals[ planeFaces.back()->GetID() ] = thePlanes[i].Direction().XYZ(); + } + + // cut theFace + Algo algo ( &mesh, theTol, normals ); + for ( size_t i = 0; i < planeFaces.size(); ++i ) + { + algo.Cut( faceToCut, planeFaces[i], 0 ); + } + + // retrieve a result + SMESH_MeshAlgos::TElemIntPairVec new2OldFaces; + SMESH_MeshAlgos::TNodeIntPairVec new2OldNodes; + TCutFaceMap::const_iterator cutFacesIt= algo.myCutFaces.cbegin(); + for ( ; cutFacesIt != algo.myCutFaces.cend(); ++cutFacesIt ) + { + const CutFace& cf = *cutFacesIt; + if ( cf.myInitFace != faceToCut ) + continue; + + if ( !cf.IsCut() ) + { + theNewFaceConnectivity.push_back( facePoints ); + break; + } + + // intersect cut lines + algo.IntersectNewEdges( cf ); + + // form loops of new faces + EdgeLoopSet loopSet; + cf.MakeLoops( loopSet, normals[ faceToCut->GetID() ]); + + // erase loops that are cut off by thePlanes + const double sign = 1; + std::vector< EdgePart > cutOffLinks; + TLinkMap cutOffCoplanarLinks; + cf.CutOffLoops( loopSet, sign, normals, cutOffLinks, cutOffCoplanarLinks ); + + for ( size_t iL = 0; iL < loopSet.myNbLoops; ++iL ) + { + EdgeLoop& loop = loopSet.myLoops[ iL ]; + if ( loop.myLinks.size() > 0 ) + { + facePoints.clear(); + for ( SMDS_NodeIteratorPtr nIt = loop.nodeIterator(); nIt->more(); ) + { + const SMDS_MeshNode* n = nIt->next(); + facePoints.push_back( n ); + int iN = faceToCut->GetNodeIndex( n ); + if ( iN < 0 ) + facePoints.back()._node = 0; // an intersection point + else + facePoints.back()._node = theFace->GetNode( iN ); + } + theNewFaceConnectivity.push_back( facePoints ); + } + } + break; + } + + return theNewFaceConnectivity.empty(); + } + +} // namespace SMESH_MeshAlgos + +namespace +{ //================================================================================ /*! * \brief Debug @@ -2092,8 +2427,9 @@ namespace // if ( !myLinks[i].IsInternal() ) // myLinks[ i ].myFace = cutterFace; // else - myLinks[ i ].ReplaceCoplanar( newEdge ); - myLinks[ i+1 ].ReplaceCoplanar( newEdge ); + myLinks[ i ].ReplaceCoplanar( newEdge ); + if ( myLinks[i].IsInternal() && i+1 < myLinks.size() ) + myLinks[ i+1 ].ReplaceCoplanar( newEdge ); return; } i += myLinks[i].IsInternal(); @@ -2189,11 +2525,11 @@ namespace bool replaced = false; for ( size_t i = 0; i < myLinks.size(); ++i ) { - while ( theRm2KeepMap.IsBound((Standard_Address) myLinks[i].myNode1 )) - replaced = ( myLinks[i].myNode1 = theRm2KeepMap((Standard_Address) myLinks[i].myNode1 )); + while ( theRm2KeepMap.IsBound( myLinks[i].myNode1 )) + replaced = ( myLinks[i].myNode1 = theRm2KeepMap( myLinks[i].myNode1 )); - while ( theRm2KeepMap.IsBound((Standard_Address) myLinks[i].myNode2 )) - replaced = ( myLinks[i].myNode2 = theRm2KeepMap((Standard_Address) myLinks[i].myNode2 )); + while ( theRm2KeepMap.IsBound( myLinks[i].myNode2 )) + replaced = ( myLinks[i].myNode2 = theRm2KeepMap( myLinks[i].myNode2 )); } //if ( replaced ) // remove equal links @@ -2277,7 +2613,8 @@ namespace //================================================================================ /*! * \brief Remove loops that are not connected to boundary edges of myFace by - * adding edges connecting these loops to the boundary + * adding edges connecting these loops to the boundary. + * Such loops must be removed as they form polygons with ring topology. */ //================================================================================ @@ -2326,14 +2663,49 @@ namespace while ( prevNbReached < nbReachedLoops ); - // add links connecting internal loops with the boundary ones for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL ) { EdgeLoop& loop = theLoops.myLoops[ iL ]; - if ( loop.myIsBndConnected ) + if ( loop.myIsBndConnected || loop.myLinks.size() == 0 ) continue; + if ( loop.myHasPending ) + { + // try to join the loop to another one, with which it contacts at a node + + // look for a node where the loop reverses + const EdgePart* edgePrev = loop.myLinks.back(); + for ( size_t iE = 0; iE < loop.myLinks.size(); edgePrev = loop.myLinks[ iE++ ] ) + { + if ( !edgePrev->IsTwin( *loop.myLinks[ iE ])) + continue; + const SMDS_MeshNode* reverseNode = edgePrev->myNode2; + + // look for a loop including reverseNode + size_t iContactEdge2; // index(+1) of edge starting at reverseNode + for ( size_t iL2 = 0; iL2 < theLoops.myNbLoops; ++iL2 ) + { + if ( iL == iL2 ) + continue; + EdgeLoop& loop2 = theLoops.myLoops[ iL2 ]; + if ( ! ( iContactEdge2 = loop2.Contains( reverseNode ))) + continue; + + // insert loop2 into the loop + theLoops.Join( loop, iE, loop2, iContactEdge2 - 1 ); + break; + } + if ( loop.myIsBndConnected ) + break; + } + + if ( loop.myIsBndConnected ) + continue; + } + + // add links connecting internal loops with the boundary ones + // find a pair of closest nodes const SMDS_MeshNode *closestNode1, *closestNode2; double minDist = 1e100; @@ -2401,14 +2773,22 @@ namespace { theLoops.AddNewLoop(); theLoops.AddEdge( myLinks[0] ); - theLoops.AddEdge( myLinks[1] ); - theLoops.AddEdge( myLinks[2] ); + if ( myLinks[0].myNode2 == myLinks[1].myNode1 ) + { + theLoops.AddEdge( myLinks[1] ); + theLoops.AddEdge( myLinks[2] ); + } + else + { + theLoops.AddEdge( myLinks[2] ); + theLoops.AddEdge( myLinks[1] ); + } return; } while ( !theLoops.AllEdgesUsed() ) { - theLoops.AddNewLoop(); + EdgeLoop& loop = theLoops.AddNewLoop(); // add 1st edge to a new loop size_t i1; @@ -2439,7 +2819,7 @@ namespace // choose among candidates if ( theLoops.myCandidates.size() == 0 ) { - theLoops.GetLoopOf( lastEdge )->myHasPending = true; + loop.myHasPending = bool( twinEdge ); lastEdge = twinEdge; } else if ( theLoops.myCandidates.size() == 1 ) @@ -2470,6 +2850,10 @@ namespace } while ( lastNode != firstNode ); + + if ( twinEdge == & myLinks[ i1 ]) + loop.myHasPending = true; + } // while ( !theLoops.AllEdgesUsed() ) return; @@ -2488,6 +2872,8 @@ namespace TLinkMap& theCutOffCoplanarLinks) const { EdgePart sideEdge; + boost::container::flat_set< const SMDS_MeshElement* > checkedCoplanar; + for ( size_t i = 0; i < myLinks.size(); ++i ) { if ( !myLinks[i].myFace ) @@ -2535,6 +2921,21 @@ namespace loop = theLoops.GetLoopOf( twin ); toErase = ( loop && !loop->myLinks.empty() ); } + + if ( toErase ) // do not erase if cutFace is connected to a co-planar cutFace + { + checkedCoplanar.clear(); + for ( size_t iE = 0; iE < myLinks.size() && toErase; ++iE ) + { + if ( !myLinks[iE].myFace || myLinks[iE].myIndex != EdgePart::_COPLANAR ) + continue; + bool isAdded = checkedCoplanar.insert( myLinks[iE].myFace ).second; + if ( !isAdded ) + continue; + toErase = SMESH_MeshAlgos::GetCommonNodes( myLinks[i ].myFace, + myLinks[iE].myFace ).size() < 1; + } + } } if ( toErase ) @@ -2546,8 +2947,8 @@ namespace { if ( !loop->myLinks[ iE ]->myFace && !loop->myLinks[ iE ]->IsInternal() )// && - // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked - // !loop->myLinks[ iE ]->myNode2->isMarked() ) + // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked + // !loop->myLinks[ iE ]->myNode2->isMarked() ) { int i = loop->myLinks[ iE ]->myIndex; sideEdge.Set( myInitFace->GetNode ( i ), @@ -2614,7 +3015,7 @@ namespace //================================================================================ /*! - * \brief Replace _COPLANAR cut edge by _INTERNAL oe vice versa + * \brief Replace _COPLANAR cut edge by _INTERNAL or vice versa */ //================================================================================ @@ -2623,7 +3024,9 @@ namespace if ( myIndex + e.myIndex == _COPLANAR + _INTERNAL ) { //check if the faces are connected - int nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e.myFace, myFace ).size(); + int nbCommonNodes = 0; + if ( e.myFace && myFace ) + nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e.myFace, myFace ).size(); bool toReplace = (( myIndex == _INTERNAL && nbCommonNodes > 1 ) || ( myIndex == _COPLANAR && nbCommonNodes < 2 )); if ( toReplace ) @@ -2642,8 +3045,8 @@ namespace /*! * \brief Create an offsetMesh of given faces * \param [in] faceIt - the input faces - * \param [out] new2OldFaces - history of faces - * \param [out] new2OldNodes - history of nodes + * \param [out] new2OldFaces - history of faces (new face -> old face ID) + * \param [out] new2OldNodes - history of nodes (new node -> old node ID) * \return SMDS_Mesh* - the new offset mesh, a caller should delete */ //================================================================================ @@ -2652,21 +3055,20 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt, SMDS_Mesh& theSrcMesh, const double theOffset, const bool theFixIntersections, - TEPairVec& theNew2OldFaces, - TNPairVec& theNew2OldNodes) + TElemIntPairVec& theNew2OldFaces, + TNodeIntPairVec& theNew2OldNodes) { - SMDS_Mesh* newMesh = new SMDS_Mesh; - theNew2OldFaces.clear(); - theNew2OldNodes.clear(); - theNew2OldFaces.push_back - ( std::make_pair(( const SMDS_MeshElement*) 0, - ( const SMDS_MeshElement*) 0)); // to have index == face->GetID() - if ( theSrcMesh.GetMeshInfo().NbFaces( ORDER_QUADRATIC ) > 0 ) throw SALOME_Exception( "Offset of quadratic mesh not supported" ); if ( theSrcMesh.GetMeshInfo().NbFaces() > theSrcMesh.GetMeshInfo().NbTriangles() ) throw SALOME_Exception( "Offset of non-triangular mesh not supported" ); + SMDS_Mesh* newMesh = new SMDS_Mesh; + theNew2OldFaces.clear(); + theNew2OldNodes.clear(); + theNew2OldFaces.push_back + ( std::make_pair(( const SMDS_MeshElement*) 0, 0)); // to have index == face->GetID() + // copy input faces to the newMesh keeping IDs of nodes double minNodeDist = 1e100; @@ -2686,7 +3088,7 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt, { SMESH_NodeXYZ xyz( nodes[i] ); newNode = newMesh->AddNodeWithID( xyz.X(), xyz.Y(), xyz.Z(), nodes[i]->GetID() ); - theNew2OldNodes.push_back( std::make_pair( newNode, nodes[i] )); + theNew2OldNodes.push_back( std::make_pair( newNode, nodes[i]->GetID() )); nodes[i] = newNode; } } @@ -2724,14 +3126,14 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt, default: continue; } - theNew2OldFaces.push_back( std::make_pair( newFace, face )); + theNew2OldFaces.push_back( std::make_pair( newFace, face->GetID() )); SMESH_NodeXYZ pPrev = nodes.back(), p; for ( size_t i = 0; i < nodes.size(); ++i ) { p.Set( nodes[i] ); double dist = ( pPrev - p ).SquareModulus(); - if ( dist > std::numeric_limits::min() ) + if ( dist < minNodeDist && dist > std::numeric_limits::min() ) minNodeDist = dist; pPrev = p; } @@ -2742,12 +3144,13 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt, std::vector< gp_XYZ > normals( theNew2OldFaces.size() ); for ( size_t i = 1; i < normals.size(); ++i ) { - if ( !SMESH_MeshAlgos::FaceNormal( theNew2OldFaces[i].second, normals[i] )) + if ( !SMESH_MeshAlgos::FaceNormal( theNew2OldFaces[i].first, normals[i] )) normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors } - const double tol = 1e-3 * Sqrt( minNodeDist ); const double sign = ( theOffset < 0 ? -1 : +1 ); + const double tol = Min( 1e-3 * Sqrt( minNodeDist ), + 1e-2 * theOffset * sign ); // translate new nodes by normal to input faces gp_XYZ newXYZ; @@ -2790,7 +3193,7 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt, { newNode = newMesh->AddNode( newXYZ.X(), newXYZ.Y(), newXYZ.Z() ); newNode->setIsMarked( true ); - theNew2OldNodes.push_back( std::make_pair( newNode, theNew2OldNodes[i].second )); + theNew2OldNodes.push_back( std::make_pair( newNode, 0 )); multiPos.emplace_back( newNode ); } } @@ -2900,7 +3303,10 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt, break; if ( !isConcaveNode2 && nbCommonNodes > 0 ) - continue; + { + if ( normals[ newFace->GetID() ] * normals[ closeFace->GetID() ] < 1.0 ) + continue; // not co-planar + } } intersector.Cut( newFace, closeFace, nbCommonNodes );