From 6ca4db2d7c49474a93e30c02bef83f05b354e540 Mon Sep 17 00:00:00 2001 From: eap Date: Fri, 22 Feb 2019 21:17:20 +0300 Subject: [PATCH] Fix SALOME_TESTS/Grids/smesh/imps_09/K0 Fix Triangulator for the case of self-intersecting but valid polygon + some fixes for #16469 1) ElementsOnShape: fix too high octree of classifiers in case of large tolerance 2) SMESH_MeshEditor::SewFreeBorder() SIGSEGV on over-constrained elements 3) Project(): adjust radius to avoid checking too many elements 4) Protect SMESH_Gen_i::Compute() from CORBA error in case of a removed mesh 5) smeshBuilder.Mesh.FindCoincidentNodesOnPart() - fix for a case of ID list --- src/Controls/SMESH_Controls.cxx | 9 +- src/SMESH/SMESH_MeshEditor.cxx | 45 ++-- src/SMESHUtils/SMESH_MeshAlgos.cxx | 63 +++-- src/SMESHUtils/SMESH_MeshAlgos.hxx | 24 +- src/SMESHUtils/SMESH_Offset.cxx | 74 +++++- src/SMESHUtils/SMESH_PolyLine.cxx | 274 ++++++++++++++++---- src/SMESHUtils/SMESH_Triangulate.cxx | 85 +++++- src/SMESH_I/SMESH_Gen_i.cxx | 18 +- src/SMESH_SWIG/smeshBuilder.py | 2 +- src/StdMeshers/StdMeshers_Projection_2D.cxx | 2 +- 10 files changed, 443 insertions(+), 153 deletions(-) diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index f18cadd34..ecd12a0b5 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -4242,6 +4242,7 @@ struct ElementsOnShape::Classifier TopAbs_ShapeEnum ShapeType() const { return myShape.ShapeType(); } const TopoDS_Shape& Shape() const { return myShape; } const Bnd_B3d* GetBndBox() const { return & myBox; } + double Tolerance() const { return myTol; } bool IsChecked() { return myFlags & theIsCheckedFlag; } bool IsSetFlag( int flag ) const { return myFlags & flag; } void SetChecked( bool is ) { is ? SetFlag( theIsCheckedFlag ) : UnsetFlag( theIsCheckedFlag ); } @@ -4458,10 +4459,9 @@ bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem) myOctree = new OctreeClassifier( myWorkClassifiers ); } - SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator(); - while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) + for ( int i = 0, nb = elem->NbNodes(); i < nb && (isSatisfy == myAllNodesFlag); ++i ) { - SMESH_TNodeXYZ aPnt( aNodeItr->next() ); + SMESH_TNodeXYZ aPnt( elem->GetNode( i )); centerXYZ += aPnt; isNodeOut = true; @@ -4816,7 +4816,8 @@ void ElementsOnShape::OctreeClassifier::buildChildrenData() for ( int i = 0; i < nbChildren(); i++ ) { OctreeClassifier* child = static_cast( myChildren[ i ]); - child->myIsLeaf = ( child->myClassifiers.size() <= 5 ); + child->myIsLeaf = ( child->myClassifiers.size() <= 5 || + child->maxSize() < child->myClassifiers[0]->Tolerance() ); } } diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index feacce422..284e83277 100644 --- a/src/SMESH/SMESH_MeshEditor.cxx +++ b/src/SMESH/SMESH_MeshEditor.cxx @@ -7588,35 +7588,37 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirst theNodes.push_back( theSecondNode ); const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode; - TIDSortedElemSet foundElems; + //TIDSortedElemSet foundElems; bool needTheLast = ( theLastNode != 0 ); + vector nodes; + while ( nStart != theLastNode ) { if ( nStart == theFirstNode ) return !needTheLast; - // find all free border faces sharing form nStart + // find all free border faces sharing nStart list< const SMDS_MeshElement* > curElemList; list< const SMDS_MeshNode* > nStartList; SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face); while ( invElemIt->more() ) { const SMDS_MeshElement* e = invElemIt->next(); - if ( e == curElem || foundElems.insert( e ).second ) { + //if ( e == curElem || foundElems.insert( e ).second ) // e can encounter twice in border + { // get nodes - int iNode = 0, nbNodes = e->NbNodes(); - vector nodes( nbNodes+1 ); nodes.assign( SMDS_MeshElement::iterator( e->interlacedNodesIterator() ), SMDS_MeshElement::iterator() ); nodes.push_back( nodes[ 0 ]); // check 2 links + int iNode = 0, nbNodes = nodes.size() - 1; for ( iNode = 0; iNode < nbNodes; iNode++ ) - if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) || - (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) && - ControlFreeBorder( &nodes[ iNode ], e->GetID() )) + if ((( nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) || + ( nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) && + ( ControlFreeBorder( &nodes[ iNode ], e->GetID() ))) { - nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]); + nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart )]); curElemList.push_back( e ); } } @@ -7668,7 +7670,7 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirst else if ( !contNodes[0].empty() && !contNodes[1].empty() ) { // choice: clear a worse one int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 ); - int iWorse = ( needTheLast ? 1 - iLongest : iLongest ); + int iWorse = ( needTheLast ? 1 - iLongest : iLongest ); contNodes[ iWorse ].clear(); contFaces[ iWorse ].clear(); } @@ -7682,10 +7684,8 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirst theNodes.pop_back(); // remove nIgnore theNodes.pop_back(); // remove nStart theFaces.pop_back(); // remove curElem - list< const SMDS_MeshNode* >::iterator nIt = cNL->begin(); - list< const SMDS_MeshElement* >::iterator fIt = cFL->begin(); - for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt ); - for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt ); + theNodes.splice( theNodes.end(), *cNL ); + theFaces.splice( theFaces.end(), *cFL ); return true; } // several continuations found @@ -8026,6 +8026,10 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode, nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin(); nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin(); + // element can be split while iterating on border if it has two edges in the border + std::map< const SMDS_MeshElement* , const SMDS_MeshElement* > elemReplaceMap; + std::map< const SMDS_MeshElement* , const SMDS_MeshElement* >::iterator elemReplaceMapIt; + TElemOfNodeListMap insertMap; TElemOfNodeListMap::iterator insertMapIt; // insertMap is @@ -8073,12 +8077,15 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode, const SMDS_MeshNode* nIns = *nIt [ 1 - intoBord ]; if ( intoBord == 1 ) { // move node of the border to be on a link of elem of the side - gp_XYZ p1 (n1->X(), n1->Y(), n1->Z()); - gp_XYZ p2 (n2->X(), n2->Y(), n2->Z()); + SMESH_NodeXYZ p1( n1 ), p2( n2 ); double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]); gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio; GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() ); } + elemReplaceMapIt = elemReplaceMap.find( elem ); + if ( elemReplaceMapIt != elemReplaceMap.end() ) + elem = elemReplaceMapIt->second; + insertMapIt = insertMap.find( elem ); bool notFound = ( insertMapIt == insertMap.end() ); bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 ); @@ -8101,8 +8108,10 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode, UpdateVolumes(n12, n22, nodeList); } // 3. find an element appeared on n1 and n2 after the insertion - insertMap.erase( elem ); - elem = findAdjacentFace( n1, n2, 0 ); + insertMap.erase( insertMapIt ); + const SMDS_MeshElement* elem2 = findAdjacentFace( n1, n2, 0 ); + elemReplaceMap.insert( std::make_pair( elem, elem2 )); + elem = elem2; } if ( notFound || otherLink ) { // add element and nodes of the side into the insertMap diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 26b3974dc..55c5f5025 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -1254,18 +1254,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 ) { @@ -2347,28 +2372,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(); } diff --git a/src/SMESHUtils/SMESH_MeshAlgos.hxx b/src/SMESHUtils/SMESH_MeshAlgos.hxx index 03e695a1d..beebbe796 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.hxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.hxx @@ -443,28 +443,12 @@ namespace SMESH_MeshAlgos bool triangulate( std::vector< const SMDS_MeshNode*>& nodes, const size_t nbNodes ); - /*! - * \brief Vertex of a polygon. Together with 2 neighbor Vertices represents a triangle - */ - struct PolyVertex - { - SMESH_NodeXYZ _nxyz; - size_t _index; - gp_XY _xy; - PolyVertex* _prev; - PolyVertex* _next; - - void SetNodeAndNext( const SMDS_MeshNode* n, PolyVertex& v, size_t index ); - void GetTriaNodes( const SMDS_MeshNode** nodes, size_t* nodeIndices) const; - double TriaArea() const; - bool IsInsideTria( const PolyVertex* v ); - PolyVertex* Delete(); - }; + struct PolyVertex; struct Optimizer; + struct Data; - std::vector< PolyVertex > _pv; - std::vector< size_t > _nodeIndex; - Optimizer* _optimizer; + Data* _data; + Optimizer* _optimizer; }; // structure used in MakePolyLine() to define a cutting plane diff --git a/src/SMESHUtils/SMESH_Offset.cxx b/src/SMESHUtils/SMESH_Offset.cxx index 44ea8dc85..b1fb119e6 100644 --- a/src/SMESHUtils/SMESH_Offset.cxx +++ b/src/SMESHUtils/SMESH_Offset.cxx @@ -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 )]; } }; @@ -2558,7 +2576,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. */ //================================================================================ @@ -2607,14 +2626,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; @@ -2689,7 +2743,7 @@ namespace while ( !theLoops.AllEdgesUsed() ) { - theLoops.AddNewLoop(); + EdgeLoop& loop = theLoops.AddNewLoop(); // add 1st edge to a new loop size_t i1; @@ -2720,7 +2774,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 ) @@ -2751,6 +2805,10 @@ namespace } while ( lastNode != firstNode ); + + if ( twinEdge == & myLinks[ i1 ]) + loop.myHasPending = true; + } // while ( !theLoops.AllEdgesUsed() ) return; diff --git a/src/SMESHUtils/SMESH_PolyLine.cxx b/src/SMESHUtils/SMESH_PolyLine.cxx index 29f6c81c0..c252afb2c 100644 --- a/src/SMESHUtils/SMESH_PolyLine.cxx +++ b/src/SMESHUtils/SMESH_PolyLine.cxx @@ -52,7 +52,17 @@ namespace int mySrcPntInd; //!< start point index TIDSortedElemSet myElemSet, myAvoidSet; - Path(): myLength(0.0), myFace(0) {} + Path(const SMDS_MeshElement* face=0, int srcInd=-1): + myLength(0.0), myFace(face), mySrcPntInd( srcInd ) {} + + void CopyNoPoints( const Path& other ); + + bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig, std::vector< Path > * paths = 0 ); + + bool SetCutAtCorner( const SMESH_NodeXYZ& cornerNode, + const gp_XYZ& plnNorm, + const gp_XYZ& plnOrig, + std::vector< Path >* paths); bool SetCutAtCorner( const SMESH_NodeXYZ& cornerNode, const SMDS_MeshElement* face, @@ -61,8 +71,6 @@ namespace void AddPoint( const gp_XYZ& p ); - bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig ); - bool ReachSamePoint( const Path& other ); static void Remove( std::vector< Path > & paths, size_t& i ); @@ -80,6 +88,25 @@ namespace myFace == other.myFace ); } + //================================================================================ + /*! + * \brief Copy data except points + */ + //================================================================================ + + void Path::CopyNoPoints( const Path& other ) + { + myLength = other.myLength; + mySrcPntInd = other.mySrcPntInd; + myFace = other.myFace; + myNode1 = other.myNode1; + myNode2 = other.myNode2; + myNodeInd1 = other.myNodeInd1; + myNodeInd2 = other.myNodeInd2; + myDot1 = other.myDot1; + myDot2 = other.myDot2; + } + //================================================================================ /*! * \brief Remove a path from a vector @@ -93,16 +120,8 @@ namespace size_t j = paths.size() - 1; // last item to be removed if ( i < j ) { + paths[ i ].CopyNoPoints ( paths[ j ]); paths[ i ].myPoints.swap( paths[ j ].myPoints ); - paths[ i ].myLength = paths[ j ].myLength; - paths[ i ].mySrcPntInd = paths[ j ].mySrcPntInd; - paths[ i ].myFace = paths[ j ].myFace; - paths[ i ].myNode1 = paths[ j ].myNode1; - paths[ i ].myNode2 = paths[ j ].myNode2; - paths[ i ].myNodeInd1 = paths[ j ].myNodeInd1; - paths[ i ].myNodeInd2 = paths[ j ].myNodeInd2; - paths[ i ].myDot1 = paths[ j ].myDot1; - paths[ i ].myDot2 = paths[ j ].myDot2; } } paths.pop_back(); @@ -110,6 +129,62 @@ namespace --i; } + //================================================================================ + /*! + * \brief Try to extend self by a point located at a node. + * Return a success flag. + */ + //================================================================================ + + bool Path::SetCutAtCorner( const SMESH_NodeXYZ& cornerNode, + const gp_XYZ& plnNorm, + const gp_XYZ& plnOrig, + std::vector< Path > * paths ) + { + bool ok = false; + const bool isContinuation = myFace; // extend this path or find all possible paths? + const SMDS_MeshElement* lastFace = myFace; + + myAvoidSet.clear(); + + SMDS_ElemIteratorPtr fIt = cornerNode->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() ) + { + Path path( lastFace, mySrcPntInd ); + if ( !path.SetCutAtCorner( cornerNode, fIt->next(), plnNorm, plnOrig )) + continue; + + if ( !myAvoidSet.insert( path.myNode1.Node() ).second || + !myAvoidSet.insert( path.myNode2.Node() ).second ) + continue; + + if ( isContinuation ) + { + if ( ok ) // non-manifold continuation + { + path.myPoints = myPoints; + path.myLength = myLength; + path.AddPoint( cornerNode ); + paths->push_back( path ); + } + else + { + double len = myLength; + this->CopyNoPoints( path ); + this->myLength = len; + this->AddPoint( path.myPoints.back() ); + } + } + else + { + paths->push_back( path ); + } + ok = true; + } + + return ok; + } + //================================================================================ /*! * \brief Store a point that is at a node of a face if the face is intersected by plane. @@ -169,11 +244,14 @@ namespace * \brief Try to find the next point * \param [in] plnNorm - cutting plane normal * \param [in] plnOrig - cutting plane origin + * \param [in] paths - all paths */ //================================================================================ - bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig ) + bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig, std::vector< Path > * paths ) { + bool ok = false; + int nodeInd3 = ( myNodeInd1 + 1 ) % myFace->NbCornerNodes(); if ( myNodeInd2 == nodeInd3 ) nodeInd3 = ( myNodeInd1 + 2 ) % myFace->NbCornerNodes(); @@ -195,20 +273,13 @@ namespace } else if ( dot3 == 0. ) { - SMDS_ElemIteratorPtr fIt = node3._node->GetInverseElementIterator(SMDSAbs_Face); - while ( fIt->more() ) - if ( SetCutAtCorner( node3, fIt->next(), plnNorm, plnOrig )) - return true; - return false; + ok = SetCutAtCorner( node3, plnNorm, plnOrig, paths ); + return ok; } else if ( myDot2 == 0. ) { - SMESH_NodeXYZ node2 = myNode2; // copy as myNode2 changes in SetCutAtCorner() - SMDS_ElemIteratorPtr fIt = node2._node->GetInverseElementIterator(SMDSAbs_Face); - while ( fIt->more() ) - if ( SetCutAtCorner( node2, fIt->next(), plnNorm, plnOrig )) - return true; - return false; + ok = SetCutAtCorner( myNode2, plnNorm, plnOrig, paths ); + return ok; } double r = Abs( myDot1 / ( myDot2 - myDot1 )); @@ -216,10 +287,32 @@ namespace myAvoidSet.clear(); myAvoidSet.insert( myFace ); - myFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node, - myElemSet, myAvoidSet, - &myNodeInd1, &myNodeInd2 ); - return myFace; + const SMDS_MeshElement* nextFace; + int ind1, ind2; + while (( nextFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node, + myElemSet, myAvoidSet, + &ind1, &ind2 ))) + { + if ( ok ) // non-manifold continuation + { + paths->push_back( *this ); + paths->back().myFace = nextFace; + paths->back().myNodeInd1 = ind1; + paths->back().myNodeInd2 = ind2; + } + else + { + myFace = nextFace; + myNodeInd1 = ind1; + myNodeInd2 = ind2; + } + ok = true; + if ( !paths ) + break; + myAvoidSet.insert( nextFace ); + } + + return ok; } //================================================================================ @@ -263,6 +356,13 @@ namespace { SMESH_MeshAlgos::PolySegment& polySeg = mySegments[ iSeg ]; + if ( ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ).SquareModulus() == 0 ) + { + myPaths[ iSeg ].AddPoint( polySeg.myXYZ[0] ); + myPaths[ iSeg ].AddPoint( polySeg.myXYZ[1] ); + return; + } + // the cutting plane gp_XYZ plnNorm = ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ) ^ polySeg.myVector.XYZ(); gp_XYZ plnOrig = polySeg.myXYZ[1]; @@ -275,14 +375,13 @@ namespace for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points { - Path path; - path.mySrcPntInd = iP; + Path path( 0, iP ); size_t nbPaths = paths.size(); if ( polySeg.myFace[ iP ]) // the end point lies on polySeg.myFace[ iP ] { // check coincidence of polySeg.myXYZ[ iP ] with nodes - const double tol = 1e-20; + const double tol = 1e-17; SMESH_NodeXYZ nodes[4]; for ( int i = 0; i < 3 && !polySeg.myNode1[ iP ]; ++i ) { @@ -292,6 +391,11 @@ namespace } nodes[ 3 ] = nodes[ 0 ]; + double dot[ 4 ]; + for ( int i = 0; i < 3; ++i ) + dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig ); + dot[ 3 ] = dot[ 0 ]; + // check coincidence of polySeg.myXYZ[ iP ] with edges for ( int i = 0; i < 3 && !polySeg.myNode1[ iP ]; ++i ) { @@ -300,16 +404,35 @@ namespace { polySeg.myNode1[ iP ] = nodes[ i ].Node(); polySeg.myNode2[ iP ] = nodes[ i + 1 ].Node(); + + int i3 = ( i + 2 ) % 3; + if ( dot[ i ] * dot [ i3 ] > 0 && + dot[ i+1 ] * dot [ i3 ] > 0 ) // point iP is inside a neighbor triangle + { + path.myAvoidSet.insert( polySeg.myFace[ iP ]); + const SMDS_MeshElement* face2 = + SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ], + polySeg.myNode2[ iP ], + path.myElemSet, + path.myAvoidSet ); + if ( face2 ) + polySeg.myFace[ iP ] = face2; + else + ;// ?? + for ( int i = 0; i < 3; ++i ) + { + nodes[ i ] = polySeg.myFace[ iP ]->GetNode( i ); + dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig ); + } + dot[ 3 ] = dot[ 0 ]; + polySeg.myNode1[ iP ] = polySeg.myNode2[ iP ] = 0; + break; + } } } if ( !polySeg.myNode1[ iP ] ) // polySeg.myXYZ[ iP ] is within polySeg.myFace[ iP ] { - double dot[ 4 ]; - for ( int i = 0; i < 3; ++i ) - dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig ); - dot[ 3 ] = dot[ 0 ]; - int iCut = 0; // index of a cut edge if ( dot[ 1 ] * dot[ 2 ] < 0. ) iCut = 1; else if ( dot[ 2 ] * dot[ 3 ] < 0. ) iCut = 2; @@ -364,27 +487,45 @@ namespace path.AddPoint( polySeg.myXYZ[ iP ]); path.myAvoidSet.insert( path.myFace ); paths.push_back( path ); + std::swap( polySeg.myNode1[ iP ], polySeg.myNode2[ iP ]); } if ( nbPaths == paths.size() ) throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1 << " in a PolySegment " << iSeg ); - } - else if ( polySeg.myNode1[ iP ] ) // the end point is at a node - { - std::set nodes; - SMDS_ElemIteratorPtr fIt = polySeg.myNode1[ iP ]->GetInverseElementIterator(SMDSAbs_Face); - while ( fIt->more() ) + + if ( path.myDot1 == 0. && + path.myDot2 == 0. && + paths.size() - nbPaths >= 2 ) // use a face non-parallel to the plane { - path.myPoints.clear(); - if ( path.SetCutAtCorner( polySeg.myNode1[ iP ], fIt->next(), plnNorm, plnOrig )) + const SMDS_MeshElement* goodFace = 0; + for ( size_t j = nbPaths; j < paths.size(); ++j ) { - if (( path.myDot1 * path.myDot2 != 0 ) || - ( nodes.insert( path.myDot1 == 0 ? path.myNode1._node : path.myNode2._node ).second )) - paths.push_back( path ); + path = paths[j]; + if ( path.Extend( plnNorm, plnOrig )) + goodFace = paths[j].myFace; + else + paths[j].myFace = 0; } + if ( !goodFace ) + throw SALOME_Exception ( SMESH_Comment("Cant move from point ") << iP+1 + << " of a PolySegment " << iSeg ); + for ( size_t j = nbPaths; j < paths.size(); ++j ) + if ( !paths[j].myFace ) + { + paths[j].myFace = goodFace; + paths[j].myNodeInd1 = goodFace->GetNodeIndex( paths[j].myNode1.Node() ); + paths[j].myNodeInd2 = goodFace->GetNodeIndex( paths[j].myNode2.Node() ); + } } } + else if ( polySeg.myNode1[ iP ] ) // the end point is at a node + { + path.myFace = 0; + path.SetCutAtCorner( polySeg.myNode1[ iP ], plnNorm, plnOrig, &paths ); + } + + // look for a one-segment path for ( size_t i = 0; i < nbPaths; ++i ) for ( size_t j = nbPaths; j < paths.size(); ++j ) @@ -394,7 +535,9 @@ namespace myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] ); paths.clear(); } - } + + } // loop on the polySeg end points to initialize all possible paths + // 2) extend paths and compose the shortest one connecting the two points @@ -405,8 +548,8 @@ namespace for ( size_t i = 0; i < paths.size(); ++i ) { Path& path = paths[ i ]; - if ( !path.Extend( plnNorm, plnOrig ) || // path reached a mesh boundary - path.myLength > myPaths[ iSeg ].myLength ) // path is longer than others + if ( !path.Extend( plnNorm, plnOrig, &paths ) || // path reached a mesh boundary + path.myLength > myPaths[ iSeg ].myLength ) // path is longer than others { Path::Remove( paths, i ); continue; @@ -428,8 +571,10 @@ namespace paths[j].myPoints.rbegin(), paths[j].myPoints.rend() ); } + if ( i < j ) std::swap( i, j ); Path::Remove( paths, i ); Path::Remove( paths, j ); + break; } } } @@ -504,7 +649,7 @@ void SMESH_MeshAlgos::MakePolyLine( SMDS_Mesh* theMes gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ(); isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits::min() ); - if ( !isVectorOK[ iSeg ]) + if ( !isVectorOK[ iSeg ] && ( p1 - p2 ).SquareModulus() > 0. ) { gp_XYZ pMid = 0.5 * ( p1 + p2 ); const SMDS_MeshElement* face; @@ -512,14 +657,35 @@ void SMESH_MeshAlgos::MakePolyLine( SMDS_Mesh* theMes polySeg.myVector = polySeg.myMidProjPoint.XYZ() - pMid; gp_XYZ faceNorm; - SMESH_MeshAlgos::FaceNormal( face, faceNorm ); + SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false ); - if ( polySeg.myVector.Magnitude() < Precision::Confusion() || - polySeg.myVector * faceNorm < Precision::Confusion() ) + const double tol = Precision::Confusion(); + if ( polySeg.myVector.Magnitude() < tol || polySeg.myVector * faceNorm < tol ) { polySeg.myVector = faceNorm; polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef; } + plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ(); + if ( plnNorm.SquareModulus() == 0 ) // p1-p2 perpendicular to mesh + { + double radius = faceNorm.Modulus(); + std::vector< const SMDS_MeshElement* > foundElems; + while ( plnNorm.SquareModulus() == 0 && radius < 1e200 ) + { + foundElems.clear(); + searcher->GetElementsInSphere( p1, radius, SMDSAbs_Face, foundElems ); + searcher->GetElementsInSphere( p2, radius, SMDSAbs_Face, foundElems ); + radius *= 2; + polySeg.myVector.SetCoord( 0,0,0 ); + for ( size_t i = 0; i < foundElems.size(); ++i ) + { + SMESH_MeshAlgos::FaceNormal( foundElems[i], faceNorm ); + polySeg.myVector += faceNorm / foundElems.size(); + } + plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ(); + } + } + } else { diff --git a/src/SMESHUtils/SMESH_Triangulate.cxx b/src/SMESHUtils/SMESH_Triangulate.cxx index 990e376d3..873eb37d7 100644 --- a/src/SMESHUtils/SMESH_Triangulate.cxx +++ b/src/SMESHUtils/SMESH_Triangulate.cxx @@ -60,13 +60,47 @@ namespace } bool operator<(const Node& other) const { return _triaIndex < other._triaIndex; } }; - typedef boost::container::flat_set< Node > TNodeSet; + typedef boost::container::flat_set< Node > TriaNodeSet; } +/*! + * \brief Vertex of a polygon. Together with 2 neighbor Vertices represents a triangle + */ +struct Triangulate::PolyVertex +{ + SMESH_NodeXYZ _nxyz; + size_t _index; + gp_XY _xy; + PolyVertex* _prev; + PolyVertex* _next; + + void SetNodeAndNext( const SMDS_MeshNode* n, PolyVertex& v, size_t index ); + void GetTriaNodes( const SMDS_MeshNode** nodes, size_t* nodeIndices) const; + double TriaArea() const; + bool IsInsideTria( const PolyVertex* v ); + PolyVertex* Delete(); + + struct Compare // compare PolyVertex'es by node + { + bool operator()(const PolyVertex* a, const PolyVertex* b) const + { + return ( a->_nxyz.Node() < b->_nxyz.Node() ); + } + }; + // set of PolyVertex sorted by mesh node + typedef boost::container::flat_set< PolyVertex*, Compare > PVSet; +}; + +struct Triangulate::Data +{ + std::vector< PolyVertex > _pv; + std::vector< size_t > _nodeIndex; + PolyVertex::PVSet _uniqueNodePV; +}; struct Triangulate::Optimizer { - std::vector< TNodeSet > _nodeUsage; // inclusions of a node in triangles + std::vector< TriaNodeSet > _nodeUsage; // inclusions of a node in triangles //================================================================================ /*! @@ -107,19 +141,19 @@ struct Triangulate::Optimizer size_t i2 = iTria + ( i + 1 ) % 3; size_t ind1 = nodeIndices[ i1 ]; // node index in points size_t ind2 = nodeIndices[ i2 ]; - TNodeSet & usage1 = _nodeUsage[ ind1 ]; // triangles using a node - TNodeSet & usage2 = _nodeUsage[ ind2 ]; + TriaNodeSet & usage1 = _nodeUsage[ ind1 ]; // triangles using a node + TriaNodeSet & usage2 = _nodeUsage[ ind2 ]; if ( usage1.size() < 2 || usage2.size() < 2 ) continue; // look for another triangle using two nodes - TNodeSet::iterator usIt1 = usage1.begin(); + TriaNodeSet::iterator usIt1 = usage1.begin(); for ( ; usIt1 != usage1.end(); ++usIt1 ) { if ( usIt1->_triaIndex == iTria ) continue; // current triangle - TNodeSet::iterator usIt2 = usage2.find( *usIt1 ); + TriaNodeSet::iterator usIt2 = usage2.find( *usIt1 ); if ( usIt2 == usage2.end() ) continue; // no common _triaIndex in two usages @@ -138,13 +172,13 @@ struct Triangulate::Optimizer // swap edge by modifying nodeIndices nodeIndices[ i2 ] = ind4; - _nodeUsage[ ind2 ].erase ({ iTria, i2 - iTria }); _nodeUsage[ ind4 ].insert({ iTria, i2 - iTria }); + _nodeUsage[ ind2 ].erase ({ iTria, i2 - iTria }); i1 = usIt1->Index(); nodeIndices[ i1 ] = ind3; - _nodeUsage[ ind1 ].erase ( *usIt1 ); _nodeUsage[ ind3 ].insert( *usIt1 ); + _nodeUsage[ ind1 ].erase ( *usIt1 ); --i; // to re-check a current edge badness1 = badness3; @@ -170,16 +204,16 @@ struct Triangulate::Optimizer std::vector< PolyVertex > & points, bool checkArea = false ) { - //if ( checkArea ) + if ( checkArea ) { points[ i2 ]._prev = & points[ i1 ]; points[ i2 ]._next = & points[ i3 ]; double a = points[ i2 ].TriaArea(); - if ( a < 0 ) - return std::numeric_limits::max(); - return 1. / a; + // if ( a < 0 ) + // return std::numeric_limits::max(); + // return 1. / a; - if ( points[ i2 ].TriaArea() < 0 ) + if ( a < 0 ) return 2; } const gp_XY & p1 = points[ i1 ]._xy; @@ -311,12 +345,28 @@ bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v ) bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes, const size_t nbNodes) { + std::vector< PolyVertex >& _pv = _data->_pv; + std::vector< size_t >& _nodeIndex = _data->_nodeIndex; + PolyVertex::PVSet& _uniqueNodePV = _data->_uniqueNodePV; + // connect nodes into a ring _pv.resize( nbNodes ); for ( size_t i = 1; i < nbNodes; ++i ) _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i], i-1 ); _pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], _pv[0], nbNodes-1 ); + // assure correctness of PolyVertex::_index as a node can encounter more than once + // within a polygon boundary + if ( _optimizer && nbNodes > 4 ) + { + _uniqueNodePV.clear(); + for ( size_t i = 0; i < nbNodes; ++i ) + { + PolyVertex::PVSet::iterator pv = _uniqueNodePV.insert( &_pv[i] ).first; + _pv[i]._index = (*pv)->_index; + } + } + // get a polygon normal gp_XYZ normal(0,0,0), p0,v01,v02; p0 = _pv[0]._nxyz; @@ -342,11 +392,16 @@ bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes, _pv[i]._xy.SetY( axes.YDirection().XYZ() * p ); } + // compute minimal triangle area + double sumArea = 0; + for ( size_t i = 0; i < nbNodes; ++i ) + sumArea += _pv[i].TriaArea(); + const double minArea = 1e-6 * sumArea / ( nbNodes - 2 ); + // in a loop, find triangles with positive area and having no vertices inside int iN = 0, nbTria = nbNodes - 2; nodes.resize( nbTria * 3 ); _nodeIndex.resize( nbTria * 3 ); - const double minArea = 1e-6; PolyVertex* v = &_pv[0], *vi; int nbVertices = nbNodes, nbBadTria = 0, isGoodTria; while ( nbBadTria < nbVertices ) @@ -430,6 +485,7 @@ bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes, Triangulate::Triangulate( bool optimize ): _optimizer(0) { + _data = new Data; if ( optimize ) _optimizer = new Optimizer; } @@ -442,6 +498,7 @@ Triangulate::Triangulate( bool optimize ): _optimizer(0) Triangulate::~Triangulate() { + delete _data; delete _optimizer; _optimizer = 0; } diff --git a/src/SMESH_I/SMESH_Gen_i.cxx b/src/SMESH_I/SMESH_Gen_i.cxx index dafa206c1..ba4898e97 100644 --- a/src/SMESH_I/SMESH_Gen_i.cxx +++ b/src/SMESH_I/SMESH_Gen_i.cxx @@ -1993,14 +1993,16 @@ CORBA::Boolean SMESH_Gen_i::Compute( SMESH::SMESH_Mesh_ptr theMesh, void SMESH_Gen_i::CancelCompute( SMESH::SMESH_Mesh_ptr theMesh, GEOM::GEOM_Object_ptr theShapeObject ) { - SMESH_Mesh_i* meshServant = dynamic_cast( GetServant( theMesh ).in() ); - ::SMESH_Mesh& myLocMesh = meshServant->GetImpl(); - TopoDS_Shape myLocShape; - if(theMesh->HasShapeToMesh()) - myLocShape = GeomObjectToShape( theShapeObject ); - else - myLocShape = SMESH_Mesh::PseudoShape(); - myGen.CancelCompute( myLocMesh, myLocShape); + if ( SMESH_Mesh_i* meshServant = dynamic_cast( GetServant( theMesh ).in() )) + { + ::SMESH_Mesh& myLocMesh = meshServant->GetImpl(); + TopoDS_Shape myLocShape; + if(theMesh->HasShapeToMesh()) + myLocShape = GeomObjectToShape( theShapeObject ); + else + myLocShape = SMESH_Mesh::PseudoShape(); + myGen.CancelCompute( myLocMesh, myLocShape); + } } //============================================================================= diff --git a/src/SMESH_SWIG/smeshBuilder.py b/src/SMESH_SWIG/smeshBuilder.py index 132ca3d43..381bb20ca 100644 --- a/src/SMESH_SWIG/smeshBuilder.py +++ b/src/SMESH_SWIG/smeshBuilder.py @@ -6300,7 +6300,7 @@ class Mesh(metaclass = MeshMeta): SubMeshOrGroup = [ obj.GetMesh() ] break if isinstance( obj, int ): - SubMeshOrGroup = self.GetIDSource( SubMeshOrGroup, SMESH.NODE ) + SubMeshOrGroup = [ self.GetIDSource( SubMeshOrGroup, SMESH.NODE )] unRegister.set( SubMeshOrGroup ) break diff --git a/src/StdMeshers/StdMeshers_Projection_2D.cxx b/src/StdMeshers/StdMeshers_Projection_2D.cxx index 9dd87361f..b88130cb2 100644 --- a/src/StdMeshers/StdMeshers_Projection_2D.cxx +++ b/src/StdMeshers/StdMeshers_Projection_2D.cxx @@ -1694,7 +1694,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& if ( !projDone || is1DComputed ) // ---------------------------------------------------------------- // The mapper can create distorted faces by placing nodes out of the FACE - // boundary, also bad face can be created if EDGEs already discretized + // boundary, also bad faces can be created if EDGEs already discretized // --> fix bad faces by smoothing // ---------------------------------------------------------------- if ( helper.IsDistorted2D( tgtSubMesh, /*checkUV=*/false, &helper )) -- 2.39.2