+
+/*//================================================================================
+/*!
+ * \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
+ * \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<TopoDS_Edge>& theWire,
+ std::vector<TopoDS_Vertex>& theVertices,
+ int & theNbDegenEdges)
+{
+ theNbDegenEdges = 0;
+
+ SMESH_MesherHelper helper( theMesh );
+
+ // sort theVertices by angle
+ multimap<double, TopoDS_Vertex> vertexByAngle;
+ TopTools_DataMapOfShapeReal angleByVertex;
+ TopoDS_Edge prevE = theWire.back();
+ if ( SMESH_Algo::isDegenerated( prevE ))
+ {
+ list<TopoDS_Edge>::reverse_iterator edge = ++theWire.rbegin();
+ while ( SMESH_Algo::isDegenerated( *edge ))
+ ++edge;
+ if ( edge == theWire.rend() )
+ return false;
+ prevE = *edge;
+ }
+ list<TopoDS_Edge>::iterator edge = theWire.begin();
+ for ( ; edge != theWire.end(); ++edge )
+ {
+ if ( SMESH_Algo::isDegenerated( *edge ))
+ {
+ ++theNbDegenEdges;
+ continue;
+ }
+ TopoDS_Vertex v = helper.IthVertex( 0, *edge );
+ if ( SMESH_Algo::VertexNode( v, helper.GetMeshDS() ))
+ {
+ double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace );
+ 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 ))
+ nbCorners = 3;
+ else
+ triaVertex.Nullify();
+
+ // check nb of available corners
+ 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<double, TopoDS_Vertex>::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<double, TopoDS_Vertex>::reverse_iterator a2v = vertexByAngle.rbegin();
+ for ( ; a2v != vertexByAngle.rend() && vMap.Extent() < nbCorners; ++a2v )
+ vMap.Add( (*a2v).second );
+
+ // check if there are possible variations in choosing corners
+ bool isThereVariants = false;
+ if ( vertexByAngle.size() > nbCorners )
+ {
+ double lostAngle = a2v->first;
+ double lastAngle = ( --a2v, a2v->first );
+ isThereVariants = ( lostAngle * 1.1 >= lastAngle );
+ }
+
+ // 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;
+ angles.reserve( vertexByAngle.size() );
+ edgeVec.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( v ));
+ edgeVec.push_back( *edge );
+ }
+
+ // refine the result vector - make sides elual by length if
+ // there are several equal angles
+ if ( isThereVariants )
+ {
+ if ( nbCorners == 3 )
+ angles[0] = 2 * M_PI; // not to move the base triangle VERTEX
+
+ set< int > refinedCorners;
+ for ( size_t iC = 0; iC < cornerInd.size(); ++iC )
+ {
+ int iV = cornerInd[iC];
+ if ( !refinedCorners.insert( iV ).second )
+ continue;
+ list< int > equalVertices;
+ equalVertices.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() );
+ int iVNext = helper.WrapIndex( iV + dV, angles.size() );
+ while ( iVNext != iV )
+ {
+ bool equal = Abs( angles[iV] - angles[iVNext] ) < 0.1 * angles[iV];
+ if ( equal )
+ equalVertices.insert( isFwd ? equalVertices.end() : equalVertices.begin(), iVNext );
+ if ( iVNext == cornerInd[ iCNext ])
+ {
+ if ( !equal )
+ break;
+ nbC[ isFwd ]++;
+ refinedCorners.insert( cornerInd[ iCNext ] );
+ iCNext = helper.WrapIndex( iCNext + dV, cornerInd.size() );
+ }
+ iVNext = helper.WrapIndex( iVNext + dV, angles.size() );
+ }
+ }
+ // move corners to make sides equal by length
+ int nbEqualV = equalVertices.size();
+ int nbExcessV = nbEqualV - ( 1 + nbC[0] + nbC[1] );
+ if ( nbExcessV > 0 )
+ {
+ // calculate normalized length of each side enclosed between neighbor equalVertices
+ vector< double > curLengths;
+ double totalLen = 0;
+ vector< int > evVec( equalVertices.begin(), equalVertices.end() );
+ int iEV = 0;
+ int iE = cornerInd[ helper.WrapIndex( iC - nbC[0] - 1, cornerInd.size() )];
+ int iEEnd = cornerInd[ helper.WrapIndex( iC + nbC[1] + 1, cornerInd.size() )];
+ while ( curLengths.size() < nbEqualV + 1 )
+ {
+ curLengths.push_back( totalLen );
+ do {
+ curLengths.back() += SMESH_Algo::EdgeLength( edgeVec[ iE ]);
+ iE = helper.WrapIndex( iE + 1, edgeVec.size());
+ if ( iEV < evVec.size() && iE == evVec[ iEV++ ] )
+ break;
+ }
+ while( iE != iEEnd );
+ totalLen = curLengths.back();
+ }
+ curLengths.resize( equalVertices.size() );
+ for ( size_t iS = 0; iS < curLengths.size(); ++iS )
+ curLengths[ iS ] /= totalLen;
+
+ // find equalVertices most close to the ideal sub-division of all sides
+ int iBestEV = 0;
+ int iCorner = helper.WrapIndex( iC - nbC[0], cornerInd.size() );
+ int nbSides = 2 + nbC[0] + nbC[1];
+ for ( int iS = 1; iS < nbSides; ++iS, ++iBestEV )
+ {
+ double idealLen = iS / double( nbSides );
+ double d, bestDist = 1.;
+ for ( iEV = iBestEV; iEV < curLengths.size(); ++iEV )
+ if (( d = Abs( idealLen - curLengths[ iEV ])) < bestDist )
+ {
+ bestDist = d;
+ iBestEV = iEV;
+ }
+ if ( iBestEV > iS-1 + nbExcessV )
+ iBestEV = iS-1 + nbExcessV;
+ theVertices[ iCorner ] = helper.IthVertex( 0, edgeVec[ evVec[ iBestEV ]]);
+ iCorner = helper.WrapIndex( iCorner + 1, cornerInd.size() );
+ }
+ }
+ }
+ }
+
+ return nbCorners;
+}
