+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 Return angle between mesh segments of given EDGEs meeting at theVertexNode
+ */
+ //================================================================================
+
+ double getAngleByNodes( const int theE1Index,
+ const int theE2Index,
+ const SMDS_MeshNode* theVertexNode,
+ const StdMeshers_FaceSide& theFaceSide,
+ const gp_Vec& theFaceNormal)
+ {
+ int eID1 = theFaceSide.EdgeID( theE1Index );
+ int eID2 = theFaceSide.EdgeID( theE2Index );
+
+ const SMDS_MeshNode *n1 = 0, *n2 = 0;
+ bool is1st;
+ SMDS_ElemIteratorPtr segIt = theVertexNode->GetInverseElementIterator( SMDSAbs_Edge );
+ while ( segIt->more() )
+ {
+ const SMDS_MeshElement* seg = segIt->next();
+ int shapeID = seg->GetShapeID();
+ if ( shapeID == eID1 )
+ is1st = true;
+ else if ( shapeID == eID2 )
+ is1st = false;
+ else
+ continue;
+ ( is1st ? n1 : n2 ) = seg->GetNodeWrap( 1 + seg->GetNodeIndex( theVertexNode ));
+ }
+
+ if ( !n1 || !n2 )
+ {
+ std::vector<const SMDS_MeshNode*> nodes;
+ for ( int is2nd = 0; is2nd < 2; ++is2nd )
+ {
+ const SMDS_MeshNode* & n = is2nd ? n2 : n1;
+ if ( n ) continue;
+ nodes.clear();
+ if ( is2nd ) theFaceSide.GetEdgeNodes( theE2Index, nodes );
+ else theFaceSide.GetEdgeNodes( theE1Index, nodes );
+ if ( nodes.size() >= 2 )
+ {
+ if ( nodes[0] == theVertexNode )
+ n = nodes[1];
+ else
+ n = nodes[ nodes.size() - 2 ];
+ }
+ }
+ }
+ double angle = -2 * M_PI;
+ if ( n1 && n2 )
+ {
+ SMESH_NodeXYZ p1 = n1, p2 = theVertexNode, p3 = n2;
+ gp_Vec v1( p1, p2 ), v2( p2, p3 );
+ try
+ {
+ angle = v1.AngleWithRef( v2, theFaceNormal );
+ }
+ catch(...)
+ {
+ }
+ if ( std::isnan( angle ))
+ angle = -2 * M_PI;
+ }
+ return angle;
+ }
+
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief EDGE of a FACE
+ */
+ struct Edge
+ {
+ TopoDS_Edge myEdge;
+ TopoDS_Vertex my1stVertex;
+ int myIndex;
+ bool myIsCorner; // is fixed corner
+ 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;
+ int myIsFixedCorner;
+ 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] ||
+ myCornerE[0] == 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] ));
+
+ myIsFixedCorner = - totNbSeg * ( myCornerE[0]->myIsCorner +
+ myCornerE[1]->myIsCorner +
+ myCornerE[2]->myIsCorner +
+ myCornerE[3]->myIsCorner );
+
+ double nbSideIdeal = totNbSeg / 4.;
+ myQuartDiff = -( Min( Min( myNbSeg[0], myNbSeg[1] ),
+ Min( myNbSeg[2], myNbSeg[3] )) / 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 + myIsFixedCorner; }
+
+ // 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, 3 or 4
+ * \param [in] theConsiderMesh - to considered only meshed VERTEXes
+ * \param [in] theFaceSide - the FACE EDGEs
+ * \param [in] theFixedVertices - VERTEXes to be used as corners
+ * \param [out] theVertices - the found corner vertices
+ * \param [out] theHaveConcaveVertices - return if there are concave vertices
+ */
+ //================================================================================
+
+ void uniteEdges( const int theNbCorners,
+ const bool theConsiderMesh,
+ const StdMeshers_FaceSide& theFaceSide,
+ const TopTools_MapOfShape& theFixedVertices,
+ std::vector<TopoDS_Vertex>& theVertices,
+ bool& theHaveConcaveVertices)
+ {
+ // 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, nbConvexAngles = 0, nbSharpAngles = 0;
+ const SMDS_MeshNode* vertNode = 0;
+ gp_Vec faceNormal;
+ const double angTol = 5. / 180 * M_PI;
+ const double sharpAngle = 0.5 * M_PI - angTol;
+ Edge* e = edge0;
+ for ( i = 0; i < nbEdges; ++i, e = e->myNext )
+ {
+ e->my1stVertex = SMESH_MesherHelper::IthVertex( 0, e->myEdge );
+ e->myIsCorner = theFixedVertices.Contains( e->my1stVertex );
+
+ e->myAngle = -2 * M_PI;
+ if ( !theConsiderMesh || ( vertNode = theFaceSide.VertexNode( e->myIndex )))
+ {
+ e->myAngle = SMESH_MesherHelper::GetAngle( e->myPrev->myEdge, e->myEdge,
+ theFaceSide.Face(), e->my1stVertex,
+ &faceNormal );
+ if ( e->myAngle > 2 * M_PI ) // GetAngle() failed
+ e->myAngle *= -1.;
+ else if ( vertNode && ( 0. <= e->myAngle ) && ( e->myAngle <= angTol ))
+ e->myAngle = getAngleByNodes( e->myPrev->myIndex, e->myIndex,
+ vertNode, theFaceSide, faceNormal );
+ }
+ edgeByAngle.insert( std::make_pair( e->myAngle, e ));
+ nbConvexAngles += ( e->myAngle > angTol );
+ nbSharpAngles += ( e->myAngle > sharpAngle );
+ }
+
+ theHaveConcaveVertices = ( nbConvexAngles < nbEdges );
+
+ if ((int) theVertices.size() == theNbCorners )
+ return;
+
+ theVertices.clear();
+
+ if ( !theConsiderMesh || theNbCorners < 4 ||
+ nbConvexAngles <= theNbCorners ||
+ nbSharpAngles == theNbCorners )
+ {
+ if ( nbEdges == theNbCorners ) // return all vertices
+ {
+ for ( e = edge0; (int) theVertices.size() < theNbCorners; e = e->myNext )
+ theVertices.push_back( e->my1stVertex );
+ return;
+ }
+
+ // return corners with maximal angles
+
+ std::set< int > cornerIndices;
+ if ( !theFixedVertices.IsEmpty() )
+ for ( i = 0, e = edge0; i < nbEdges; ++i, e = e->myNext )
+ if ( e->myIsCorner )
+ cornerIndices.insert( e->myIndex );
+
+ 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<const SMDS_MeshNode*> nodes;
+ for ( i = 0, e = edge0; i < nbEdges; ++i, e = e->myNext )
+ {
+ nodes.clear();
+ theFaceSide.GetEdgeNodes( e->myIndex, nodes, /*addVertex=*/true, true );
+ if ( nodes.size() == 2 && nodes[0] == nodes[1] ) // all nodes merged
+ {
+ e->myAngle = -1; // to remove
+ }
+ else
+ {
+ 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;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Remove a seam and degenerated edge from a wire if the shape is
+ * a quadrangle with a seam inside.
+ */
+ //================================================================================
+
+ bool removeInternalSeam( std::list<TopoDS_Edge>& theWire,
+ SMESH_MesherHelper& theHelper)
+ {
+ if ( !theHelper.HasRealSeam() ||
+ theHelper.NbDegeneratedEdges() != 2 ) // 1 EDGE + 1 VERTEX
+ return false;
+
+ typedef std::list<TopoDS_Edge>::iterator TEdgeIter;
+ std::vector< TEdgeIter > edgesToRemove;
+ edgesToRemove.reserve( 5 );
+ for ( TEdgeIter eIt = theWire.begin(); eIt != theWire.end(); ++eIt )
+ {
+ int eID = theHelper.ShapeToIndex( *eIt );
+ if ( theHelper.IsRealSeam( eID ) || theHelper.IsDegenShape( eID ))
+ edgesToRemove.push_back( eIt );
+ }
+
+ if ( theWire.size() - edgesToRemove.size() < 4 )
+ return false; // cone e.g.
+
+ for ( size_t i = 0; i < edgesToRemove.size(); ++i )
+ theWire.erase( edgesToRemove[ i ]);
+
+ return true;
+ }
+
+} // namespace
+