From: eap Date: Mon, 26 Feb 2018 16:06:46 +0000 (+0300) Subject: 23521: EDF 16246 - problems with quadrangles in use case X-Git-Tag: V8_5_0a1~3 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=f7712f9c03c5553fc299b9bd9d28d81ce1c83666;p=modules%2Fsmesh.git 23521: EDF 16246 - problems with quadrangles in use case --- diff --git a/src/StdMeshers/StdMeshers_FaceSide.hxx b/src/StdMeshers/StdMeshers_FaceSide.hxx index 1964638fe..86fe12c62 100644 --- a/src/StdMeshers/StdMeshers_FaceSide.hxx +++ b/src/StdMeshers/StdMeshers_FaceSide.hxx @@ -207,7 +207,7 @@ public: double constValue = 0) const; /*! * \brief Return nodes in the order they encounter while walking along - * the while side or a specified EDGE. For a closed side, the 1st point repeats at end. + * the whole side or a specified EDGE. For a closed side, the 1st point repeats at end. * \param iE - index of the EDGE. Default is "all EDGEs". */ std::vector GetOrderedNodes( int iE = ALL_EDGES ) const; diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 679769a8f..d3d071863 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -65,10 +65,11 @@ #include "Utils_ExceptHandlers.hxx" #include +#include typedef NCollection_Array2 StdMeshers_Array2OfNode; -typedef gp_XY gp_UV; +typedef gp_XY gp_UV; typedef SMESH_Comment TComm; using namespace std; @@ -1020,30 +1021,429 @@ bool StdMeshers_Quadrangle_2D::IsApplicable( const TopoDS_Shape & aShape, bool t return ( toCheckAll && nbFoundFaces != 0 ); } +namespace +{ + //================================================================================ + /*! + * \brief Return true if only two given edges meat at their common vertex + */ + //================================================================================ + + bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1, + const TopoDS_Edge& e2, + SMESH_Mesh & mesh) + { + TopoDS_Vertex v; + if (!TopExp::CommonVertex(e1, e2, v)) + return false; + TopTools_ListIteratorOfListOfShape ancestIt(mesh.GetAncestors(v)); + for (; ancestIt.More() ; ancestIt.Next()) + if (ancestIt.Value().ShapeType() == TopAbs_EDGE) + if (!e1.IsSame(ancestIt.Value()) && !e2.IsSame(ancestIt.Value())) + return false; + return true; + } + + //-------------------------------------------------------------------------------- + /*! + * \brief EDGE of a FACE + */ + struct Edge + { + TopoDS_Edge myEdge; + TopoDS_Vertex my1stVertex; + int myIndex; + double myAngle; // angle at my1stVertex + int myNbSegments; // discretization + Edge* myPrev; // preceding EDGE + Edge* myNext; // next EDGE + + // traits used by boost::intrusive::circular_list_algorithms + typedef Edge node; + typedef Edge * node_ptr; + typedef const Edge * const_node_ptr; + static node_ptr get_next(const_node_ptr n) { return n->myNext; } + static void set_next(node_ptr n, node_ptr next) { n->myNext = next; } + static node_ptr get_previous(const_node_ptr n) { return n->myPrev; } + static void set_previous(node_ptr n, node_ptr prev){ n->myPrev = prev; } + }; + + //-------------------------------------------------------------------------------- + /*! + * \brief Four sides of a quadrangle evaluating its quality + */ + struct QuadQuality + { + typedef std::set< QuadQuality, QuadQuality > set; + + Edge* myCornerE[4]; + int myNbSeg [4]; + + // quality criteria to minimize + int myOppDiff; + double myQuartDiff; + double mySumAngle; + + // Compute quality criateria and add self to the set of variants + // + void AddSelf( QuadQuality::set& theVariants ) + { + if ( myCornerE[2] == myCornerE[1] || // exclude invalid variants + myCornerE[2] == myCornerE[3] ) + return; + + // count nb segments between corners + mySumAngle = 0; + double totNbSeg = 0; + for ( int i1 = 3, i2 = 0; i2 < 4; i1 = i2++ ) + { + myNbSeg[ i1 ] = 0; + for ( Edge* e = myCornerE[ i1 ]; e != myCornerE[ i2 ]; e = e->myNext ) + myNbSeg[ i1 ] += e->myNbSegments; + mySumAngle -= myCornerE[ i1 ]->myAngle / M_PI; // [-1,1] + totNbSeg += myNbSeg[ i1 ]; + } + + myOppDiff = ( Abs( myNbSeg[0] - myNbSeg[2] ) + + Abs( myNbSeg[1] - myNbSeg[3] )); + + double nbSideIdeal = totNbSeg / 4.; + myQuartDiff = -( Min( Min( myNbSeg[0], myNbSeg[1] ), + Min( myNbSeg[1], myNbSeg[2] )) / nbSideIdeal ); + + theVariants.insert( *this ); + +#ifndef _DEBUG_ + if ( theVariants.size() > 1 ) // erase a worse variant + theVariants.erase( ++theVariants.begin() ); +#endif + }; + + // first criterion - equality of nbSeg of opposite sides + int crit1() const { return myOppDiff; } + + // second criterion - equality of nbSeg of adjacent sides and sharpness of angles + double crit2() const { return myQuartDiff + mySumAngle; } + + bool operator () ( const QuadQuality& q1, const QuadQuality& q2) const + { + if ( q1.crit1() < q2.crit1() ) + return true; + if ( q1.crit1() > q2.crit1() ) + return false; + return q1.crit2() < q2.crit2(); + } + }; + + //================================================================================ + /*! + * \brief Unite EDGEs to get a required number of sides + * \param [in] theNbCorners - the required number of sides + * \param [in] theConsiderMesh - to considered only meshed VERTEXes + * \param [in] theFaceSide - the FACE EDGEs + * \param [out] theVertices - the found corner vertices + */ + //================================================================================ + + void uniteEdges( const int theNbCorners, + const bool theConsiderMesh, + const StdMeshers_FaceSide& theFaceSide, + const TopoDS_Shape& theBaseVertex, + std::vector& theVertices ) + { + theVertices.clear(); + + // form a circular list of EDGEs + std::vector< Edge > edges( theFaceSide.NbEdges() ); + boost::intrusive::circular_list_algorithms< Edge > circularList; + circularList.init_header( &edges[0] ); + edges[0].myEdge = theFaceSide.Edge( 0 ); + edges[0].myIndex = 0; + edges[0].myNbSegments = 0; + for ( int i = 1; i < theFaceSide.NbEdges(); ++i ) + { + edges[ i ].myEdge = theFaceSide.Edge( i ); + edges[ i ].myIndex = i; + edges[ i ].myNbSegments = 0; + circularList.link_after( &edges[ i-1 ], &edges[ i ] ); + } + // remove degenerated edges + int nbEdges = edges.size(); + Edge* edge0 = &edges[0]; + for ( size_t i = 0; i < edges.size(); ++i ) + if ( SMESH_Algo::isDegenerated( edges[i].myEdge )) + { + edge0 = circularList.unlink( &edges[i] ); + --nbEdges; + } + + // sort edges by angle + std::multimap< double, Edge* > edgeByAngle; + int i, iBase = -1, nbConvexAngles = 0; + Edge* e = edge0; + for ( i = 0; i < nbEdges; ++i, e = e->myNext ) + { + e->my1stVertex = SMESH_MesherHelper::IthVertex( 0, e->myEdge ); + if ( e->my1stVertex.IsSame( theBaseVertex )) + iBase = e->myIndex; + + e->myAngle = -2 * M_PI; + if ( !theConsiderMesh || theFaceSide.VertexNode( e->myIndex )) + { + e->myAngle = SMESH_MesherHelper::GetAngle( e->myPrev->myEdge, e->myEdge, + theFaceSide.Face(), e->my1stVertex ); + if ( e->myAngle > 2 * M_PI ) // GetAngle() failed + e->myAngle *= -1.; + } + edgeByAngle.insert( std::make_pair( e->myAngle, e )); + nbConvexAngles += ( e->myAngle > 0 ); + } + + if ( !theConsiderMesh || theNbCorners < 4 || nbConvexAngles <= theNbCorners ) + { + // return corners with maximal angles + + std::set< int > cornerIndices; + if ( iBase != -1 ) + cornerIndices.insert( iBase ); + + std::multimap< double, Edge* >::reverse_iterator a2e = edgeByAngle.rbegin(); + for (; (int) cornerIndices.size() < theNbCorners; ++a2e ) + cornerIndices.insert( a2e->second->myIndex ); + + std::set< int >::iterator i = cornerIndices.begin(); + for ( ; i != cornerIndices.end(); ++i ) + theVertices.push_back( edges[ *i ].my1stVertex ); + + return; + } + + // get nb of segments + int totNbSeg = 0; // tatal nb segments + std::vector nodes; + for ( i = 0, e = edge0; i < nbEdges; ++i, e = e->myNext ) + { + nodes.clear(); + theFaceSide.GetEdgeNodes( e->myIndex, nodes, /*addVertex=*/false, false ); + e->myNbSegments += nodes.size() + 1; + totNbSeg += nodes.size() + 1; + + // join with the previous edge those edges with concave angles + if ( e->myAngle <= 0 ) + { + e->myPrev->myNbSegments += e->myNbSegments; + e = circularList.unlink( e )->myPrev; + --nbEdges; + --i; + } + } + + if ( edge0->myNext->myPrev != edge0 ) // edge0 removed, find another edge0 + for ( size_t i = 0; i < edges.size(); ++i ) + if ( edges[i].myNext->myPrev == & edges[i] ) + { + edge0 = &edges[i]; + break; + } + + + // sort different variants by quality + + QuadQuality::set quadVariants; + + // find index of a corner most opposite to corner of edge0 + int iOpposite0, nbHalf = 0; + for ( e = edge0; nbHalf <= totNbSeg / 2; e = e->myNext ) + nbHalf += e->myNbSegments; + iOpposite0 = e->myIndex; + + // compose different variants of quadrangles + QuadQuality quad; + for ( ; edge0->myIndex != iOpposite0; edge0 = edge0->myNext ) + { + quad.myCornerE[ 0 ] = edge0; + + // find opposite corner 2 + for ( nbHalf = 0, e = edge0; nbHalf < totNbSeg / 2; e = e->myNext ) + nbHalf += e->myNbSegments; + if ( e == edge0->myNext ) // no space for corner 1 + e = e->myNext; + quad.myCornerE[ 2 ] = e; + + bool moreVariants2 = ( totNbSeg % 2 || nbHalf != totNbSeg / 2 ); + + // enumerate different variants of corners 1 and 3 + for ( Edge* e1 = edge0->myNext; e1 != quad.myCornerE[ 2 ]; e1 = e1->myNext ) + { + quad.myCornerE[ 1 ] = e1; + + // find opposite corner 3 + for ( nbHalf = 0, e = e1; nbHalf < totNbSeg / 2; e = e->myNext ) + nbHalf += e->myNbSegments; + if ( e == quad.myCornerE[ 2 ] ) + e = e->myNext; + quad.myCornerE[ 3 ] = e; + + bool moreVariants3 = ( totNbSeg % 2 || nbHalf != totNbSeg / 2 ); + + quad.AddSelf( quadVariants ); + + // another variants + if ( moreVariants2 ) + { + quad.myCornerE[ 2 ] = quad.myCornerE[ 2 ]->myPrev; + quad.AddSelf( quadVariants ); + quad.myCornerE[ 2 ] = quad.myCornerE[ 2 ]->myNext; + } + if ( moreVariants3 ) + { + quad.myCornerE[ 3 ] = quad.myCornerE[ 3 ]->myPrev; + quad.AddSelf( quadVariants ); + + if ( moreVariants2 ) + { + quad.myCornerE[ 2 ] = quad.myCornerE[ 2 ]->myPrev; + quad.AddSelf( quadVariants ); + quad.myCornerE[ 2 ] = quad.myCornerE[ 2 ]->myNext; + } + } + } + } + + const QuadQuality& bestQuad = *quadVariants.begin(); + theVertices.resize( 4 ); + theVertices[ 0 ] = bestQuad.myCornerE[ 0 ]->my1stVertex; + theVertices[ 1 ] = bestQuad.myCornerE[ 1 ]->my1stVertex; + theVertices[ 2 ] = bestQuad.myCornerE[ 2 ]->my1stVertex; + theVertices[ 3 ] = bestQuad.myCornerE[ 3 ]->my1stVertex; + + return; + } + +} // namespace + //================================================================================ /*! - * \brief Return true if only two given edges meat at their common vertex + * \brief Finds vertices at the most sharp face corners + * \param [in] theFace - the FACE + * \param [in,out] theWire - the ordered edges of the face. It can be modified to + * have the first VERTEX of the first EDGE in \a vertices + * \param [out] theVertices - the found corner vertices in the order corresponding to + * the order of EDGEs in \a theWire + * \param [out] theNbDegenEdges - nb of degenerated EDGEs in theFace + * \param [in] theConsiderMesh - if \c true, only meshed VERTEXes are considered + * as possible corners + * \return int - number of quad sides found: 0, 3 or 4 */ //================================================================================ -static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1, - const TopoDS_Edge& e2, - SMESH_Mesh & mesh) +int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, + SMESH_Mesh & theMesh, + std::list& theWire, + std::vector& theVertices, + int & theNbDegenEdges, + const bool theConsiderMesh) { - TopoDS_Vertex v; - if (!TopExp::CommonVertex(e1, e2, v)) - return false; - TopTools_ListIteratorOfListOfShape ancestIt(mesh.GetAncestors(v)); - for (; ancestIt.More() ; ancestIt.Next()) - if (ancestIt.Value().ShapeType() == TopAbs_EDGE) - if (!e1.IsSame(ancestIt.Value()) && !e2.IsSame(ancestIt.Value())) - return false; - return true; + theNbDegenEdges = 0; + + SMESH_MesherHelper helper( theMesh ); + if ( myHelper ) + helper.CopySubShapeInfo( *myHelper ); + + StdMeshers_FaceSide faceSide( theFace, theWire, &theMesh, + /*isFwd=*/true, /*skipMedium=*/true, &helper ); + + // count degenerated EDGEs and possible corner VERTEXes + for ( int iE = 0; iE < faceSide.NbEdges(); ++iE ) + { + if ( SMESH_Algo::isDegenerated( faceSide.Edge( iE ))) + ++theNbDegenEdges; + else if ( !theConsiderMesh || faceSide.VertexNode( iE )) + theVertices.push_back( faceSide.FirstVertex( iE )); + } + + // find out required nb of corners (3 or 4) + int nbCorners = 4; + TopoDS_Shape triaVertex = helper.GetMeshDS()->IndexToShape( myTriaVertexID ); + if ( !triaVertex.IsNull() && + triaVertex.ShapeType() == TopAbs_VERTEX && + helper.IsSubShape( triaVertex, theFace ) && + theVertices.size() != 4 ) + nbCorners = 3; + else + triaVertex.Nullify(); + + // check nb of available EDGEs + if ( faceSide.NbEdges() < nbCorners ) + return error(COMPERR_BAD_SHAPE, + TComm("Face must have 4 sides and not ") << faceSide.NbEdges() ); + + if ( theConsiderMesh ) + { + const int nbSegments = Max( faceSide.NbPoints()-1, faceSide.NbSegments() ); + if ( nbSegments < nbCorners ) + return error(COMPERR_BAD_INPUT_MESH, TComm("Too few boundary nodes: ") << nbSegments); + } + + if ( nbCorners == 3 ) + { + if ( theVertices.size() < 3 ) + return error(COMPERR_BAD_SHAPE, + TComm("Face must have 3 meshed sides and not ") << theVertices.size() ); + } + else // triaVertex not defined or invalid + { + if ( theVertices.size() == 3 && theNbDegenEdges == 0 ) + { + if ( myTriaVertexID < 1 ) + return error(COMPERR_BAD_PARMETERS, + "No Base vertex provided for a trilateral geometrical face"); + + TComm comment("Invalid Base vertex: "); + comment << myTriaVertexID << ", which is not in [ "; + comment << helper.GetMeshDS()->ShapeToIndex( faceSide.FirstVertex(0) ) << ", "; + comment << helper.GetMeshDS()->ShapeToIndex( faceSide.FirstVertex(1) ) << ", "; + comment << helper.GetMeshDS()->ShapeToIndex( faceSide.FirstVertex(2) ) << " ]"; + return error(COMPERR_BAD_PARMETERS, comment ); + } + if ( theVertices.size() + theNbDegenEdges < 4 ) + return error(COMPERR_BAD_SHAPE, + TComm("Face must have 4 meshed sides and not ") << theVertices.size() ); + } + + if ((int) theVertices.size() > nbCorners ) + { + // there are more EDGEs than required nb of sides; + // unite some EDGEs to fix the nb of sides + uniteEdges( nbCorners, theConsiderMesh, faceSide, triaVertex, theVertices ); + } + + if ( nbCorners == 3 && !triaVertex.IsSame( theVertices[0] )) + { + // make theVertices begin from triaVertex + for ( size_t i = 0; i < theVertices.size(); ++i ) + if ( triaVertex.IsSame( theVertices[i] )) + { + theVertices.erase( theVertices.begin(), theVertices.begin() + i ); + break; + } + else + { + theVertices.push_back( theVertices[i] ); + } + } + + // make theWire begin from the 1st corner vertex + while ( !theVertices[0].IsSame( helper.IthVertex( 0, theWire.front() )) || + SMESH_Algo::isDegenerated( theWire.front() )) + theWire.splice( theWire.end(), theWire, theWire.begin() ); + + return nbCorners; } //============================================================================= /*! - * + * */ //============================================================================= @@ -4253,388 +4653,6 @@ bool StdMeshers_Quadrangle_2D::check() return isOK; } -//================================================================================ -/*! - * \brief Finds vertices at the most sharp face corners - * \param [in] theFace - the FACE - * \param [in,out] theWire - the ordered edges of the face. It can be modified to - * have the first VERTEX of the first EDGE in \a vertices - * \param [out] theVertices - the found corner vertices in the order corresponding to - * the order of EDGEs in \a theWire - * \param [out] theNbDegenEdges - nb of degenerated EDGEs in theFace - * \param [in] theConsiderMesh - if \c true, only meshed VERTEXes are considered - * as possible corners - * \return int - number of quad sides found: 0, 3 or 4 - */ -//================================================================================ - -int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, - SMESH_Mesh & theMesh, - std::list& theWire, - std::vector& theVertices, - int & theNbDegenEdges, - const bool theConsiderMesh) -{ - theNbDegenEdges = 0; - - SMESH_MesherHelper helper( theMesh ); - if ( myHelper ) - helper.CopySubShapeInfo( *myHelper ); - StdMeshers_FaceSide faceSide( theFace, theWire, &theMesh, - /*isFwd=*/true, /*skipMedium=*/true, &helper ); - - // sort theVertices by angle - multimap vertexByAngle; - TopTools_DataMapOfShapeReal angleByVertex; - TopoDS_Edge prevE = theWire.back(); - if ( SMESH_Algo::isDegenerated( prevE )) - { - list::reverse_iterator edge = ++theWire.rbegin(); - while ( SMESH_Algo::isDegenerated( *edge ) /*|| helper.IsRealSeam( *edge )*/) - ++edge; - if ( edge == theWire.rend() ) - return false; - prevE = *edge; - } - list::iterator edge = theWire.begin(); - for ( int iE = 0; edge != theWire.end(); ++edge, ++iE ) - { - if ( SMESH_Algo::isDegenerated( *edge ) /*|| helper.IsRealSeam( *edge )*/) - { - ++theNbDegenEdges; - continue; - } - if ( !theConsiderMesh || faceSide.VertexNode( iE )) - { - TopoDS_Vertex v = helper.IthVertex( 0, *edge ); - double angle = helper.GetAngle( prevE, *edge, theFace, v ); - vertexByAngle.insert( make_pair( angle, v )); - angleByVertex.Bind( v, angle ); - } - prevE = *edge; - } - - // find out required nb of corners (3 or 4) - int nbCorners = 4; - TopoDS_Shape triaVertex = helper.GetMeshDS()->IndexToShape( myTriaVertexID ); - if ( !triaVertex.IsNull() && - triaVertex.ShapeType() == TopAbs_VERTEX && - helper.IsSubShape( triaVertex, theFace ) && - ( vertexByAngle.size() != 4 || vertexByAngle.begin()->first < 5 * M_PI/180. )) - nbCorners = 3; - else - triaVertex.Nullify(); - - // check nb of available corners - if ( faceSide.NbEdges() < nbCorners ) - return error(COMPERR_BAD_SHAPE, - TComm("Face must have 4 sides but not ") << faceSide.NbEdges() ); - - if ( theConsiderMesh ) - { - const int nbSegments = Max( faceSide.NbPoints()-1, faceSide.NbSegments() ); - if ( nbSegments < nbCorners ) - return error(COMPERR_BAD_INPUT_MESH, TComm("Too few boundary nodes: ") << nbSegments); - } - - if ( nbCorners == 3 ) - { - if ( vertexByAngle.size() < 3 ) - return error(COMPERR_BAD_SHAPE, - TComm("Face must have 3 sides but not ") << vertexByAngle.size() ); - } - else - { - if ( vertexByAngle.size() == 3 && theNbDegenEdges == 0 ) - { - if ( myTriaVertexID < 1 ) - return error(COMPERR_BAD_PARMETERS, - "No Base vertex provided for a trilateral geometrical face"); - - TComm comment("Invalid Base vertex: "); - comment << myTriaVertexID << " its ID is not among [ "; - multimap::iterator a2v = vertexByAngle.begin(); - comment << helper.GetMeshDS()->ShapeToIndex( a2v->second ) << ", "; a2v++; - comment << helper.GetMeshDS()->ShapeToIndex( a2v->second ) << ", "; a2v++; - comment << helper.GetMeshDS()->ShapeToIndex( a2v->second ) << " ]"; - return error(COMPERR_BAD_PARMETERS, comment ); - } - if ( vertexByAngle.size() + ( theNbDegenEdges > 0 ) < 4 && - vertexByAngle.size() + theNbDegenEdges != 4 ) - return error(COMPERR_BAD_SHAPE, - TComm("Face must have 4 sides but not ") << vertexByAngle.size() ); - } - - // put all corner vertices in a map - TopTools_MapOfShape vMap; - if ( nbCorners == 3 ) - vMap.Add( triaVertex ); - multimap::reverse_iterator a2v = vertexByAngle.rbegin(); - for ( int iC = 0; a2v != vertexByAngle.rend() && iC < nbCorners; ++a2v, ++iC ) - vMap.Add( (*a2v).second ); - - // check if there are possible variations in choosing corners - bool haveVariants = false; - if ((int) vertexByAngle.size() > nbCorners ) - { - double lostAngle = a2v->first; - double lastAngle = ( --a2v, a2v->first ); - haveVariants = ( lostAngle * 1.1 >= lastAngle ); - } - - const double angleTol = 5.* M_PI/180; - myCheckOri = ( (int)vertexByAngle.size() > nbCorners || - vertexByAngle.begin()->first < angleTol ); - - // make theWire begin from a corner vertex or triaVertex - if ( nbCorners == 3 ) - while ( !triaVertex.IsSame( ( helper.IthVertex( 0, theWire.front() ))) || - SMESH_Algo::isDegenerated( theWire.front() )) - theWire.splice( theWire.end(), theWire, theWire.begin() ); - else - while ( !vMap.Contains( helper.IthVertex( 0, theWire.front() )) || - SMESH_Algo::isDegenerated( theWire.front() )) - theWire.splice( theWire.end(), theWire, theWire.begin() ); - - // fill the result vector and prepare for its refinement - theVertices.clear(); - vector< double > angles; - vector< TopoDS_Edge > edgeVec; - vector< int > cornerInd, nbSeg; - int nbSegTot = 0; - angles .reserve( vertexByAngle.size() ); - edgeVec.reserve( vertexByAngle.size() ); - nbSeg .reserve( vertexByAngle.size() ); - cornerInd.reserve( nbCorners ); - for ( edge = theWire.begin(); edge != theWire.end(); ++edge ) - { - if ( SMESH_Algo::isDegenerated( *edge )) - continue; - TopoDS_Vertex v = helper.IthVertex( 0, *edge ); - bool isCorner = vMap.Contains( v ); - if ( isCorner ) - { - theVertices.push_back( v ); - cornerInd.push_back( angles.size() ); - } - angles .push_back( angleByVertex.IsBound( v ) ? angleByVertex( v ) : -M_PI ); - edgeVec.push_back( *edge ); - if ( theConsiderMesh && haveVariants ) - { - if ( SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( *edge )) - nbSeg.push_back( sm->NbNodes() + 1 ); - else - nbSeg.push_back( 0 ); - nbSegTot += nbSeg.back(); - } - } - - // refine the result vector - make sides equal by length if - // there are several equal angles - if ( haveVariants ) - { - if ( nbCorners == 3 ) - angles[0] = 2 * M_PI; // not to move the base triangle VERTEX - - // here we refer to VERTEX'es and EDGEs by indices in angles and edgeVec vectors - typedef int TGeoIndex; - - // for each vertex find a vertex till which there are nbSegHalf segments - const int nbSegHalf = ( nbSegTot % 2 || nbCorners == 3 ) ? 0 : nbSegTot / 2; - vector< TGeoIndex > halfDivider( angles.size(), -1 ); - int nbHalfDividers = 0; - if ( nbSegHalf ) - { - // get min angle of corners - double minAngle = 10.; - for ( size_t iC = 0; iC < cornerInd.size(); ++iC ) - minAngle = Min( minAngle, angles[ cornerInd[ iC ]]); - - // find halfDivider's - for ( TGeoIndex iV1 = 0; iV1 < TGeoIndex( angles.size() ); ++iV1 ) - { - int nbSegs = 0; - TGeoIndex iV2 = iV1; - do { - nbSegs += nbSeg[ iV2 ]; - iV2 = helper.WrapIndex( iV2 + 1, nbSeg.size() ); - } while ( nbSegs < nbSegHalf ); - - if ( nbSegs == nbSegHalf && - angles[ iV1 ] + angleTol >= minAngle && - angles[ iV2 ] + angleTol >= minAngle ) - { - halfDivider[ iV1 ] = iV2; - ++nbHalfDividers; - } - } - } - - set< TGeoIndex > refinedCorners, treatedCorners; - for ( size_t iC = 0; iC < cornerInd.size(); ++iC ) - { - TGeoIndex iV = cornerInd[iC]; - if ( !treatedCorners.insert( iV ).second ) - continue; - list< TGeoIndex > equVerts; // inds of vertices that can become corners - equVerts.push_back( iV ); - int nbC[2] = { 0, 0 }; - // find equal angles backward and forward from the iV-th corner vertex - for ( int isFwd = 0; isFwd < 2; ++isFwd ) - { - int dV = isFwd ? +1 : -1; - int iCNext = helper.WrapIndex( iC + dV, cornerInd.size() ); - TGeoIndex iVNext = helper.WrapIndex( iV + dV, angles.size() ); - while ( iVNext != iV ) - { - bool equal = Abs( angles[iV] - angles[iVNext] ) < angleTol; - if ( equal ) - equVerts.insert( isFwd ? equVerts.end() : equVerts.begin(), iVNext ); - if ( iVNext == cornerInd[ iCNext ]) - { - if ( !equal ) - { - if ( angles[iV] < angles[iVNext] ) - refinedCorners.insert( iVNext ); - break; - } - nbC[ isFwd ]++; - treatedCorners.insert( cornerInd[ iCNext ] ); - iCNext = helper.WrapIndex( iCNext + dV, cornerInd.size() ); - } - iVNext = helper.WrapIndex( iVNext + dV, angles.size() ); - } - if ( iVNext == iV ) - break; // all angles equal - } - - const bool allCornersSame = ( nbC[0] == 3 ); - if ( allCornersSame && nbHalfDividers > 0 ) - { - // select two halfDivider's as corners - TGeoIndex hd1, hd2 = -1; - size_t iC2; - for ( iC2 = 0; iC2 < cornerInd.size() && hd2 < 0; ++iC2 ) - { - hd1 = cornerInd[ iC2 ]; - hd2 = halfDivider[ hd1 ]; - if ( std::find( equVerts.begin(), equVerts.end(), hd2 ) == equVerts.end() ) - hd2 = -1; // hd2-th vertex can't become a corner - else - break; - } - if ( hd2 >= 0 ) - { - angles[ hd1 ] = 2 * M_PI; // make hd1-th vertex no more "equal" - angles[ hd2 ] = 2 * M_PI; - refinedCorners.insert( hd1 ); - refinedCorners.insert( hd2 ); - treatedCorners = refinedCorners; - // update cornerInd - equVerts.push_front( equVerts.back() ); - equVerts.push_back( equVerts.front() ); - list< TGeoIndex >::iterator hdPos = - std::find( equVerts.begin(), equVerts.end(), hd2 ); - if ( hdPos == equVerts.end() ) break; - cornerInd[ helper.WrapIndex( iC2 + 0, cornerInd.size()) ] = hd1; - cornerInd[ helper.WrapIndex( iC2 + 1, cornerInd.size()) ] = *( --hdPos ); - cornerInd[ helper.WrapIndex( iC2 + 2, cornerInd.size()) ] = hd2; - cornerInd[ helper.WrapIndex( iC2 + 3, cornerInd.size()) ] = *( ++hdPos, ++hdPos ); - - theVertices[ 0 ] = helper.IthVertex( 0, edgeVec[ cornerInd[0] ]); - theVertices[ 1 ] = helper.IthVertex( 0, edgeVec[ cornerInd[1] ]); - theVertices[ 2 ] = helper.IthVertex( 0, edgeVec[ cornerInd[2] ]); - theVertices[ 3 ] = helper.IthVertex( 0, edgeVec[ cornerInd[3] ]); - iC = -1; - continue; - } - } - - // move corners to make sides equal by length - int nbEqualV = equVerts.size(); - int nbExcessV = nbEqualV - ( 1 + nbC[0] + nbC[1] ); - if ( nbExcessV > 0 ) // there are nbExcessV vertices that can become corners - { - // calculate normalized length of each "side" enclosed between neighbor equVerts - vector< double > accuLength; - double totalLen = 0; - vector< TGeoIndex > evVec( equVerts.begin(), equVerts.end() ); - size_t iEV = 0; - TGeoIndex iE = cornerInd[ helper.WrapIndex( iC - nbC[0] - 1, cornerInd.size() )]; - TGeoIndex iEEnd = cornerInd[ helper.WrapIndex( iC + nbC[1] + 1, cornerInd.size() )]; - while ((int) accuLength.size() < nbEqualV + int( !allCornersSame ) ) - { - // accumulate length of edges before iEV-th equal vertex - accuLength.push_back( totalLen ); - do { - accuLength.back() += SMESH_Algo::EdgeLength( edgeVec[ iE ]); - iE = helper.WrapIndex( iE + 1, edgeVec.size()); - if ( iEV < evVec.size() && iE == evVec[ iEV ] ) { - iEV++; - break; // equal vertex reached - } - } - while( iE != iEEnd ); - totalLen = accuLength.back(); - } - accuLength.resize( equVerts.size() ); - for ( size_t iS = 0; iS < accuLength.size(); ++iS ) - accuLength[ iS ] /= totalLen; - - // find equVerts most close to the ideal sub-division of all sides - int iBestEV = 0; - int iCorner = helper.WrapIndex( iC - nbC[0], cornerInd.size() ); - int nbSides = Min( nbCorners, 2 + nbC[0] + nbC[1] ); - for ( int iS = 1; iS < nbSides; ++iS, ++iBestEV ) - { - double idealLen = iS / double( nbSides ); - double d, bestDist = 2.; - for ( iEV = iBestEV; iEV < accuLength.size(); ++iEV ) - { - d = Abs( idealLen - accuLength[ iEV ]); - - // take into account presence of a corresponding halfDivider - const double cornerWgt = 0.5 / nbSides; - const double vertexWgt = 0.25 / nbSides; - TGeoIndex hd = halfDivider[ evVec[ iEV ]]; - if ( hd < 0 ) - d += vertexWgt; - else if( refinedCorners.count( hd )) - d -= cornerWgt; - else - d -= vertexWgt; - - // choose vertex with the best d - if ( d < bestDist ) - { - bestDist = d; - iBestEV = iEV; - } - } - if ( iBestEV > iS-1 + nbExcessV ) - iBestEV = iS-1 + nbExcessV; - theVertices[ iCorner ] = helper.IthVertex( 0, edgeVec[ evVec[ iBestEV ]]); - cornerInd [ iCorner ] = evVec[ iBestEV ]; - refinedCorners.insert( evVec[ iBestEV ]); - iCorner = helper.WrapIndex( iCorner + 1, cornerInd.size() ); - } - - } // if ( nbExcessV > 0 ) - else - { - refinedCorners.insert( cornerInd[ iC ]); - } - } // loop on cornerInd - - // make theWire begin from the cornerInd[0]-th EDGE - while ( !theWire.front().IsSame( edgeVec[ cornerInd[0] ])) - theWire.splice( theWire.begin(), theWire, --theWire.end() ); - - } // if ( haveVariants ) - - return nbCorners; -} - //================================================================================ /*! * \brief Constructor of a side of quad