-/*//================================================================================
-/*!
- * \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<TopoDS_Edge>& theWire,
- std::vector<TopoDS_Vertex>& theVertices,
- int & theNbDegenEdges,
- const bool theConsiderMesh)
-{
- 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 ( !theConsiderMesh || SMESH_Algo::VertexNode( v, helper.GetMeshDS() ))
- {
- double angle = SMESH_MesherHelper::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 ( 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 ( 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 ( 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 = ( 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;
- int 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 is 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() );
- int 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 ( 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 coresponding 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 ]]);
- 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;
-}
-