+ SaveFacet savedFacet( myCurFace );
+
+ // concave polyhedron
+
+ if ( myFwdLinks.empty() ) // get links of the least ambiguously oriented facet
+ {
+ for ( size_t i = 0; i < myPolyFacetOri.size() && !ori; ++i )
+ ori = myPolyFacetOri[ i ];
+
+ if ( !ori ) // none facet is oriented yet
+ {
+ // find the least ambiguously oriented facet
+ int faceMostConvex = -1;
+ std::map< double, int > convexity2face;
+ for ( size_t iF = 0; iF < myPolyFacetOri.size() && faceMostConvex < 0; ++iF )
+ {
+ if ( projectNodesToNormal( iF, minProj, maxProj ))
+ {
+ // all nodes are on the same side of the facet
+ me->myPolyFacetOri[ iF ] = ( minProj < 0 ? +1 : -1 );
+ faceMostConvex = iF;
+ }
+ else
+ {
+ ori = ( -minProj < maxProj ? -1 : +1 );
+ double convexity = std::min( -minProj, maxProj ) / std::max( -minProj, maxProj );
+ convexity2face.insert( make_pair( convexity, iF * ori ));
+ }
+ }
+ if ( faceMostConvex < 0 ) // none facet has nodes on the same side
+ {
+ // use the least ambiguous facet
+ faceMostConvex = convexity2face.begin()->second;
+ ori = ( faceMostConvex < 0 ? -1 : +1 );
+ faceMostConvex = std::abs( faceMostConvex );
+ me->myPolyFacetOri[ faceMostConvex ] = ori;
+ }
+ }
+ // collect links of the oriented facets in myFwdLinks
+ for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
+ {
+ ori = myPolyFacetOri[ iF ];
+ if ( !ori ) continue;
+ setFace( iF );
+ for ( int i = 0; i < myCurFace.myNbNodes; ++i )
+ {
+ NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
+ me->myFwdLinks.insert( make_pair( link, link.myOri ));
+ }
+ }
+ }
+
+ // compare orientation of links of the facet with myFwdLinks
+ ori = 0;
+ setFace( faceIndex );
+ vector< NLink > links( myCurFace.myNbNodes ), links2;
+ for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
+ {
+ NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
+ std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
+ if ( l2o != myFwdLinks.end() )
+ ori = link.myOri * l2o->second * -1;
+ links[ i ] = link;
+ }
+ while ( !ori ) // the facet has no common links with already oriented facets
+ {
+ // orient and collect links of other non-oriented facets
+ for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
+ {
+ if ( myPolyFacetOri[ iF ] ) continue; // already oriented
+ setFace( iF );
+ links2.clear();
+ ori = 0;
+ for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
+ {
+ NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
+ std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
+ if ( l2o != myFwdLinks.end() )
+ ori = link.myOri * l2o->second * -1;
+ links2.push_back( link );
+ }
+ if ( ori ) // one more facet oriented
+ {
+ me->myPolyFacetOri[ iF ] = ori;
+ for ( size_t i = 0; i < links2.size(); ++i )
+ me->myFwdLinks.insert( make_pair( links2[i], links2[i].myOri * ori ));
+ break;
+ }
+ }
+ if ( !ori )
+ return false; // error in algorithm: infinite loop
+
+ // try to orient the facet again
+ ori = 0;
+ for ( size_t i = 0; i < links.size() && !ori; ++i )
+ {
+ std::map<Link, int>::const_iterator l2o = myFwdLinks.find( links[i] );
+ if ( l2o != myFwdLinks.end() )
+ ori = links[i].myOri * l2o->second * -1;
+ }
+ me->myPolyFacetOri[ faceIndex ] = ori;
+ }
+
+ return ori > 0;
+}
+
+//=======================================================================
+//function : projectNodesToNormal
+//purpose : compute min and max projections of all nodes to normal of a facet.
+//=======================================================================
+
+bool SMDS_VolumeTool::projectNodesToNormal( int faceIndex,
+ double& minProj,
+ double& maxProj ) const
+{
+ minProj = std::numeric_limits<double>::max();
+ maxProj = std::numeric_limits<double>::min();
+
+ XYZ normal;
+ if ( !GetFaceNormal( faceIndex, normal.x, normal.y, normal.z ))
+ return false;
+ XYZ p0 ( myCurFace.myNodes[0] );
+ for ( size_t i = 0; i < myVolumeNodes.size(); ++i )
+ {
+ if ( std::find( myCurFace.myNodes.begin() + 1,
+ myCurFace.myNodes.end(),
+ myVolumeNodes[ i ] ) != myCurFace.myNodes.end() )
+ continue; // node of the faceIndex-th facet
+
+ double proj = normal.Dot( XYZ( myVolumeNodes[ i ]) - p0 );
+ if ( proj < minProj ) minProj = proj;
+ if ( proj > maxProj ) maxProj = proj;
+ }
+ const double tol = 1e-7;
+ minProj += tol;
+ maxProj -= tol;
+ bool diffSize = ( minProj * maxProj < 0 );
+ // if ( diffSize )
+ // {
+ // minProj = -minProj;