+ return SMDSAbs_Face;
+}
+
+
+/*
+ Class : AspectRatio3D
+ Description : Functor for calculating aspect ratio
+*/
+namespace{
+
+ inline double getHalfPerimeter(double theTria[3]){
+ return (theTria[0] + theTria[1] + theTria[2])/2.0;
+ }
+
+ inline double getArea(double theHalfPerim, double theTria[3]){
+ return sqrt(theHalfPerim*
+ (theHalfPerim-theTria[0])*
+ (theHalfPerim-theTria[1])*
+ (theHalfPerim-theTria[2]));
+ }
+
+ inline double getVolume(double theLen[6]){
+ double a2 = theLen[0]*theLen[0];
+ double b2 = theLen[1]*theLen[1];
+ double c2 = theLen[2]*theLen[2];
+ double d2 = theLen[3]*theLen[3];
+ double e2 = theLen[4]*theLen[4];
+ double f2 = theLen[5]*theLen[5];
+ double P = 4.0*a2*b2*d2;
+ double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
+ double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
+ return sqrt(P-Q+R)/12.0;
+ }
+
+ inline double getVolume2(double theLen[6]){
+ double a2 = theLen[0]*theLen[0];
+ double b2 = theLen[1]*theLen[1];
+ double c2 = theLen[2]*theLen[2];
+ double d2 = theLen[3]*theLen[3];
+ double e2 = theLen[4]*theLen[4];
+ double f2 = theLen[5]*theLen[5];
+
+ double P = a2*e2*(b2+c2+d2+f2-a2-e2);
+ double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
+ double R = c2*d2*(a2+b2+e2+f2-c2-d2);
+ double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
+
+ return sqrt(P+Q+R-S)/12.0;
+ }
+
+ inline double getVolume(const TSequenceOfXYZ& P){
+ gp_Vec aVec1( P( 2 ) - P( 1 ) );
+ gp_Vec aVec2( P( 3 ) - P( 1 ) );
+ gp_Vec aVec3( P( 4 ) - P( 1 ) );
+ gp_Vec anAreaVec( aVec1 ^ aVec2 );
+ return fabs(aVec3 * anAreaVec) / 6.0;
+ }
+
+ inline double getMaxHeight(double theLen[6])
+ {
+ double aHeight = std::max(theLen[0],theLen[1]);
+ aHeight = std::max(aHeight,theLen[2]);
+ aHeight = std::max(aHeight,theLen[3]);
+ aHeight = std::max(aHeight,theLen[4]);
+ aHeight = std::max(aHeight,theLen[5]);
+ return aHeight;
+ }
+
+}
+
+double AspectRatio3D::GetValue( long theId )
+{
+ double aVal = 0;
+ myCurrElement = myMesh->FindElement( theId );
+ if ( myCurrElement && myCurrElement->GetVtkType() == VTK_TETRA )
+ {
+ // Action from CoTech | ACTION 31.3:
+ // EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with
+ // those of ParaView. The library used by ParaView for those calculations can be reused in SMESH.
+ vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
+ if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
+ aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
+ }
+ else
+ {
+ TSequenceOfXYZ P;
+ if ( GetPoints( myCurrElement, P ))
+ aVal = Round( GetValue( P ));
+ }
+ return aVal;
+}
+
+double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
+{
+ double aQuality = 0.0;
+ if(myCurrElement->IsPoly()) return aQuality;
+
+ int nbNodes = P.size();
+
+ if(myCurrElement->IsQuadratic()) {
+ if(nbNodes==10) nbNodes=4; // quadratic tetrahedron
+ else if(nbNodes==13) nbNodes=5; // quadratic pyramid
+ else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
+ else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
+ else if(nbNodes==27) nbNodes=8; // quadratic hexahedron
+ else return aQuality;
+ }
+
+ switch(nbNodes) {
+ case 4:{
+ double aLen[6] = {
+ getDistance(P( 1 ),P( 2 )), // a
+ getDistance(P( 2 ),P( 3 )), // b
+ getDistance(P( 3 ),P( 1 )), // c
+ getDistance(P( 2 ),P( 4 )), // d
+ getDistance(P( 3 ),P( 4 )), // e
+ getDistance(P( 1 ),P( 4 )) // f
+ };
+ double aTria[4][3] = {
+ {aLen[0],aLen[1],aLen[2]}, // abc
+ {aLen[0],aLen[3],aLen[5]}, // adf
+ {aLen[1],aLen[3],aLen[4]}, // bde
+ {aLen[2],aLen[4],aLen[5]} // cef
+ };
+ double aSumArea = 0.0;
+ double aHalfPerimeter = getHalfPerimeter(aTria[0]);
+ double anArea = getArea(aHalfPerimeter,aTria[0]);
+ aSumArea += anArea;
+ aHalfPerimeter = getHalfPerimeter(aTria[1]);
+ anArea = getArea(aHalfPerimeter,aTria[1]);
+ aSumArea += anArea;
+ aHalfPerimeter = getHalfPerimeter(aTria[2]);
+ anArea = getArea(aHalfPerimeter,aTria[2]);
+ aSumArea += anArea;
+ aHalfPerimeter = getHalfPerimeter(aTria[3]);
+ anArea = getArea(aHalfPerimeter,aTria[3]);
+ aSumArea += anArea;
+ double aVolume = getVolume(P);
+ //double aVolume = getVolume(aLen);
+ double aHeight = getMaxHeight(aLen);
+ static double aCoeff = sqrt(2.0)/12.0;
+ if ( aVolume > DBL_MIN )
+ aQuality = aCoeff*aHeight*aSumArea/aVolume;
+ break;
+ }
+ case 5:{
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ break;
+ }
+ case 6:{
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ break;
+ }
+ case 8:{
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ }
+ break;
+ }
+ case 12:
+ {
+ gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[8] = {P( 3 ),P( 4 ),P( 6 ),P( 1 ),P( 9 ),P( 10 ),P( 12 ),P( 7 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
+ }
+ break;
+ } // switch(nbNodes)
+
+ if ( nbNodes > 4 ) {
+ // avaluate aspect ratio of quadranle faces
+ AspectRatio aspect2D;
+ SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
+ int nbFaces = SMDS_VolumeTool::NbFaces( type );
+ TSequenceOfXYZ points(4);
+ for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume
+ if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
+ continue;
+ const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
+ for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face
+ points( p + 1 ) = P( pInd[ p ] + 1 );
+ aQuality = std::max( aQuality, aspect2D.GetValue( points ));
+ }
+ }
+ return aQuality;
+}
+
+double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // the aspect ratio is in the range [1.0,infinity]
+ // 1.0 = good
+ // infinity = bad
+ return Value / 1000.;
+}
+
+SMDSAbs_ElementType AspectRatio3D::GetType() const
+{
+ return SMDSAbs_Volume;
+}
+
+
+/*
+ Class : Warping
+ Description : Functor for calculating warping
+*/
+double Warping::GetValue( const TSequenceOfXYZ& P )
+{
+ if ( P.size() != 4 )
+ return 0;
+
+ gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.;
+
+ double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
+ double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
+ double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
+ double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
+
+ return Max( Max( A1, A2 ), Max( A3, A4 ) );
+}
+
+double Warping::ComputeA( const gp_XYZ& thePnt1,
+ const gp_XYZ& thePnt2,
+ const gp_XYZ& thePnt3,
+ const gp_XYZ& theG ) const
+{
+ double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
+ double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
+ double L = Min( aLen1, aLen2 ) * 0.5;
+ if ( L < Precision::Confusion())
+ return 0.;
+
+ gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
+ gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
+ gp_XYZ N = GI.Crossed( GJ );
+
+ if ( N.Modulus() < gp::Resolution() )
+ return M_PI / 2;
+
+ N.Normalize();
+
+ double H = ( thePnt2 - theG ).Dot( N );
+ return asin( fabs( H / L ) ) * 180. / M_PI;
+}
+
+double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // the warp is in the range [0.0,PI/2]
+ // 0.0 = good (no warp)
+ // PI/2 = bad (face pliee)
+ return Value;
+}
+
+SMDSAbs_ElementType Warping::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
+
+/*
+ Class : Taper
+ Description : Functor for calculating taper
+*/
+double Taper::GetValue( const TSequenceOfXYZ& P )
+{
+ if ( P.size() != 4 )
+ return 0.;
+
+ // Compute taper
+ double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2.;
+ double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2.;
+ double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2.;
+ double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2.;
+
+ double JA = 0.25 * ( J1 + J2 + J3 + J4 );
+ if ( JA <= Precision::Confusion() )
+ return 0.;
+
+ double T1 = fabs( ( J1 - JA ) / JA );
+ double T2 = fabs( ( J2 - JA ) / JA );
+ double T3 = fabs( ( J3 - JA ) / JA );
+ double T4 = fabs( ( J4 - JA ) / JA );
+
+ return Max( Max( T1, T2 ), Max( T3, T4 ) );
+}
+
+double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // the taper is in the range [0.0,1.0]
+ // 0.0 = good (no taper)
+ // 1.0 = bad (les cotes opposes sont allignes)
+ return Value;
+}
+
+SMDSAbs_ElementType Taper::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
+
+/*
+ Class : Skew
+ Description : Functor for calculating skew in degrees
+*/
+static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
+{
+ gp_XYZ p12 = ( p2 + p1 ) / 2.;
+ gp_XYZ p23 = ( p3 + p2 ) / 2.;
+ gp_XYZ p31 = ( p3 + p1 ) / 2.;
+
+ gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
+
+ return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
+}
+
+double Skew::GetValue( const TSequenceOfXYZ& P )
+{
+ if ( P.size() != 3 && P.size() != 4 )
+ return 0.;
+
+ // Compute skew
+ static double PI2 = M_PI / 2.;
+ if ( P.size() == 3 )
+ {
+ double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
+ double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
+ double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
+
+ return Max( A0, Max( A1, A2 ) ) * 180. / M_PI;
+ }
+ else
+ {
+ gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.;
+ gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.;
+ gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.;
+ gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.;
+
+ gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
+ double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
+ ? 0. : fabs( PI2 - v1.Angle( v2 ) );
+
+ //BUG SWP12743
+ if ( A < Precision::Angular() )
+ return 0.;
+
+ return A * 180. / M_PI;
+ }
+}
+
+double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // the skew is in the range [0.0,PI/2].
+ // 0.0 = good
+ // PI/2 = bad
+ return Value;
+}
+
+SMDSAbs_ElementType Skew::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
+
+/*
+ Class : Area
+ Description : Functor for calculating area
+*/
+double Area::GetValue( const TSequenceOfXYZ& P )
+{
+ double val = 0.0;
+ if ( P.size() > 2 ) {
+ gp_Vec aVec1( P(2) - P(1) );
+ gp_Vec aVec2( P(3) - P(1) );
+ gp_Vec SumVec = aVec1 ^ aVec2;
+ for (int i=4; i<=P.size(); i++) {
+ gp_Vec aVec1( P(i-1) - P(1) );
+ gp_Vec aVec2( P(i) - P(1) );
+ gp_Vec tmp = aVec1 ^ aVec2;
+ SumVec.Add(tmp);
+ }
+ val = SumVec.Magnitude() * 0.5;
+ }
+ return val;
+}
+
+double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // meaningless as it is not a quality control functor
+ return Value;
+}
+
+SMDSAbs_ElementType Area::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
+
+/*
+ Class : Length
+ Description : Functor for calculating length of edge
+*/
+double Length::GetValue( const TSequenceOfXYZ& P )
+{
+ switch ( P.size() ) {
+ case 2: return getDistance( P( 1 ), P( 2 ) );
+ case 3: return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) );
+ default: return 0.;
+ }
+}
+
+double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // meaningless as it is not quality control functor
+ return Value;
+}
+
+SMDSAbs_ElementType Length::GetType() const
+{
+ return SMDSAbs_Edge;
+}
+
+/*
+ Class : Length2D
+ Description : Functor for calculating length of edge
+*/
+
+double Length2D::GetValue( long theElementId)
+{
+ TSequenceOfXYZ P;
+
+ //cout<<"Length2D::GetValue"<<endl;
+ if (GetPoints(theElementId,P)){
+ //for(int jj=1; jj<=P.size(); jj++)
+ // cout<<"jj="<<jj<<" P("<<P(jj).X()<<","<<P(jj).Y()<<","<<P(jj).Z()<<")"<<endl;
+
+ double aVal;// = GetValue( P );
+ const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
+ SMDSAbs_ElementType aType = aElem->GetType();
+
+ int len = P.size();
+
+ switch (aType){
+ case SMDSAbs_All:
+ case SMDSAbs_Node:
+ case SMDSAbs_Edge:
+ if (len == 2){
+ aVal = getDistance( P( 1 ), P( 2 ) );
+ break;
+ }
+ else if (len == 3){ // quadratic edge
+ aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
+ break;
+ }
+ case SMDSAbs_Face:
+ if (len == 3){ // triangles
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 1 ));
+ aVal = Max(L1,Max(L2,L3));
+ break;
+ }
+ else if (len == 4){ // quadrangles
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 1 ));
+ aVal = Max(Max(L1,L2),Max(L3,L4));
+ break;
+ }
+ if (len == 6){ // quadratic triangles
+ double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+ double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+ double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
+ aVal = Max(L1,Max(L2,L3));
+ //cout<<"L1="<<L1<<" L2="<<L2<<"L3="<<L3<<" aVal="<<aVal<<endl;
+ break;
+ }
+ else if (len == 8){ // quadratic quadrangles
+ double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+ double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+ double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
+ double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
+ aVal = Max(Max(L1,L2),Max(L3,L4));
+ break;
+ }
+ case SMDSAbs_Volume:
+ if (len == 4){ // tetraidrs
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 1 ));
+ double L4 = getDistance(P( 1 ),P( 4 ));
+ double L5 = getDistance(P( 2 ),P( 4 ));
+ double L6 = getDistance(P( 3 ),P( 4 ));
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ break;
+ }
+ else if (len == 5){ // piramids
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 1 ));
+ double L5 = getDistance(P( 1 ),P( 5 ));
+ double L6 = getDistance(P( 2 ),P( 5 ));
+ double L7 = getDistance(P( 3 ),P( 5 ));
+ double L8 = getDistance(P( 4 ),P( 5 ));
+
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(L7,L8));
+ break;
+ }
+ else if (len == 6){ // pentaidres
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 1 ));
+ double L4 = getDistance(P( 4 ),P( 5 ));
+ double L5 = getDistance(P( 5 ),P( 6 ));
+ double L6 = getDistance(P( 6 ),P( 4 ));
+ double L7 = getDistance(P( 1 ),P( 4 ));
+ double L8 = getDistance(P( 2 ),P( 5 ));
+ double L9 = getDistance(P( 3 ),P( 6 ));
+
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),L9));
+ break;
+ }
+ else if (len == 8){ // hexaider
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 1 ));
+ double L5 = getDistance(P( 5 ),P( 6 ));
+ double L6 = getDistance(P( 6 ),P( 7 ));
+ double L7 = getDistance(P( 7 ),P( 8 ));
+ double L8 = getDistance(P( 8 ),P( 5 ));
+ double L9 = getDistance(P( 1 ),P( 5 ));
+ double L10= getDistance(P( 2 ),P( 6 ));
+ double L11= getDistance(P( 3 ),P( 7 ));
+ double L12= getDistance(P( 4 ),P( 8 ));
+
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+ aVal = Max(aVal,Max(L11,L12));
+ break;
+
+ }
+
+ if (len == 10){ // quadratic tetraidrs
+ double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
+ double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+ double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
+ double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ break;
+ }
+ else if (len == 13){ // quadratic piramids
+ double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+ double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+ double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
+ double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
+ double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(L7,L8));
+ break;
+ }
+ else if (len == 15){ // quadratic pentaidres
+ double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+ double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+ double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
+ double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
+ double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
+ double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
+ double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),L9));
+ break;
+ }
+ else if (len == 20){ // quadratic hexaider
+ double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
+ double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
+ double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
+ double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
+ double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
+ double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
+ double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
+ double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
+ double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+ aVal = Max(aVal,Max(L11,L12));
+ break;
+
+ }
+
+ default: aVal=-1;
+ }
+
+ if (aVal <0){
+ return 0.;
+ }
+
+ if ( myPrecision >= 0 )
+ {
+ double prec = pow( 10., (double)( myPrecision ) );
+ aVal = floor( aVal * prec + 0.5 ) / prec;
+ }
+
+ return aVal;
+
+ }
+ return 0.;
+}
+
+double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // meaningless as it is not quality control functor
+ return Value;
+}
+
+SMDSAbs_ElementType Length2D::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
+Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
+ myLength(theLength)
+{
+ myPntId[0] = thePntId1; myPntId[1] = thePntId2;
+ if(thePntId1 > thePntId2){
+ myPntId[1] = thePntId1; myPntId[0] = thePntId2;
+ }
+}
+
+bool Length2D::Value::operator<(const Length2D::Value& x) const{
+ if(myPntId[0] < x.myPntId[0]) return true;
+ if(myPntId[0] == x.myPntId[0])
+ if(myPntId[1] < x.myPntId[1]) return true;
+ return false;
+}
+
+void Length2D::GetValues(TValues& theValues){
+ TValues aValues;
+ SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
+ for(; anIter->more(); ){
+ const SMDS_MeshFace* anElem = anIter->next();
+
+ if(anElem->IsQuadratic()) {
+ const SMDS_VtkFace* F =
+ dynamic_cast<const SMDS_VtkFace*>(anElem);
+ // use special nodes iterator
+ SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
+ long aNodeId[4];
+ gp_Pnt P[4];
+
+ double aLength;
+ const SMDS_MeshElement* aNode;
+ if(anIter->more()){
+ aNode = anIter->next();
+ const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
+ P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
+ aNodeId[0] = aNodeId[1] = aNode->GetID();
+ aLength = 0;
+ }
+ for(; anIter->more(); ){
+ const SMDS_MeshNode* N1 = static_cast<const SMDS_MeshNode*> (anIter->next());
+ P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z());
+ aNodeId[2] = N1->GetID();
+ aLength = P[1].Distance(P[2]);
+ if(!anIter->more()) break;
+ const SMDS_MeshNode* N2 = static_cast<const SMDS_MeshNode*> (anIter->next());
+ P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z());
+ aNodeId[3] = N2->GetID();
+ aLength += P[2].Distance(P[3]);
+ Value aValue1(aLength,aNodeId[1],aNodeId[2]);
+ Value aValue2(aLength,aNodeId[2],aNodeId[3]);
+ P[1] = P[3];
+ aNodeId[1] = aNodeId[3];
+ theValues.insert(aValue1);
+ theValues.insert(aValue2);
+ }
+ aLength += P[2].Distance(P[0]);
+ Value aValue1(aLength,aNodeId[1],aNodeId[2]);
+ Value aValue2(aLength,aNodeId[2],aNodeId[0]);
+ theValues.insert(aValue1);
+ theValues.insert(aValue2);
+ }
+ else {
+ SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
+ long aNodeId[2];
+ gp_Pnt P[3];
+
+ double aLength;
+ const SMDS_MeshElement* aNode;
+ if(aNodesIter->more()){
+ aNode = aNodesIter->next();
+ const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
+ P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
+ aNodeId[0] = aNodeId[1] = aNode->GetID();
+ aLength = 0;
+ }
+ for(; aNodesIter->more(); ){
+ aNode = aNodesIter->next();
+ const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
+ long anId = aNode->GetID();
+
+ P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
+
+ aLength = P[1].Distance(P[2]);
+
+ Value aValue(aLength,aNodeId[1],anId);
+ aNodeId[1] = anId;
+ P[1] = P[2];
+ theValues.insert(aValue);
+ }
+
+ aLength = P[0].Distance(P[1]);
+
+ Value aValue(aLength,aNodeId[0],aNodeId[1]);
+ theValues.insert(aValue);
+ }
+ }
+}
+
+/*
+ Class : MultiConnection
+ Description : Functor for calculating number of faces conneted to the edge
+*/
+double MultiConnection::GetValue( const TSequenceOfXYZ& P )
+{
+ return 0;
+}
+double MultiConnection::GetValue( long theId )
+{
+ return getNbMultiConnection( myMesh, theId );
+}
+
+double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // meaningless as it is not quality control functor
+ return Value;
+}
+
+SMDSAbs_ElementType MultiConnection::GetType() const
+{
+ return SMDSAbs_Edge;
+}
+
+/*
+ Class : MultiConnection2D
+ Description : Functor for calculating number of faces conneted to the edge
+*/
+double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
+{
+ return 0;
+}
+
+double MultiConnection2D::GetValue( long theElementId )
+{
+ int aResult = 0;
+
+ const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
+ SMDSAbs_ElementType aType = aFaceElem->GetType();
+
+ switch (aType) {
+ case SMDSAbs_Face:
+ {
+ int i = 0, len = aFaceElem->NbNodes();
+ SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
+ if (!anIter) break;
+
+ const SMDS_MeshNode *aNode, *aNode0;
+ TColStd_MapOfInteger aMap, aMapPrev;
+
+ for (i = 0; i <= len; i++) {
+ aMapPrev = aMap;
+ aMap.Clear();
+
+ int aNb = 0;
+ if (anIter->more()) {
+ aNode = (SMDS_MeshNode*)anIter->next();
+ } else {
+ if (i == len)
+ aNode = aNode0;
+ else
+ break;
+ }
+ if (!aNode) break;
+ if (i == 0) aNode0 = aNode;
+
+ SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
+ while (anElemIter->more()) {
+ const SMDS_MeshElement* anElem = anElemIter->next();
+ if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
+ int anId = anElem->GetID();
+
+ aMap.Add(anId);
+ if (aMapPrev.Contains(anId)) {
+ aNb++;
+ }
+ }
+ }
+ aResult = Max(aResult, aNb);
+ }
+ }
+ break;
+ default:
+ aResult = 0;
+ }
+
+ return aResult;
+}
+
+double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // meaningless as it is not quality control functor
+ return Value;
+}
+
+SMDSAbs_ElementType MultiConnection2D::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
+MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
+{
+ myPntId[0] = thePntId1; myPntId[1] = thePntId2;
+ if(thePntId1 > thePntId2){
+ myPntId[1] = thePntId1; myPntId[0] = thePntId2;
+ }
+}
+
+bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{
+ if(myPntId[0] < x.myPntId[0]) return true;
+ if(myPntId[0] == x.myPntId[0])
+ if(myPntId[1] < x.myPntId[1]) return true;
+ return false;
+}
+
+void MultiConnection2D::GetValues(MValues& theValues){
+ SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
+ for(; anIter->more(); ){
+ const SMDS_MeshFace* anElem = anIter->next();
+ SMDS_ElemIteratorPtr aNodesIter;
+ if ( anElem->IsQuadratic() )
+ aNodesIter = dynamic_cast<const SMDS_VtkFace*>
+ (anElem)->interlacedNodesElemIterator();
+ else
+ aNodesIter = anElem->nodesIterator();
+ long aNodeId[3];
+
+ //int aNbConnects=0;
+ const SMDS_MeshNode* aNode0;
+ const SMDS_MeshNode* aNode1;
+ const SMDS_MeshNode* aNode2;
+ if(aNodesIter->more()){
+ aNode0 = (SMDS_MeshNode*) aNodesIter->next();
+ aNode1 = aNode0;
+ const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
+ aNodeId[0] = aNodeId[1] = aNodes->GetID();
+ }
+ for(; aNodesIter->more(); ) {
+ aNode2 = (SMDS_MeshNode*) aNodesIter->next();
+ long anId = aNode2->GetID();
+ aNodeId[2] = anId;
+
+ Value aValue(aNodeId[1],aNodeId[2]);
+ MValues::iterator aItr = theValues.find(aValue);
+ if (aItr != theValues.end()){
+ aItr->second += 1;
+ //aNbConnects = nb;
+ }
+ else {
+ theValues[aValue] = 1;
+ //aNbConnects = 1;
+ }
+ //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
+ aNodeId[1] = aNodeId[2];
+ aNode1 = aNode2;
+ }
+ Value aValue(aNodeId[0],aNodeId[2]);
+ MValues::iterator aItr = theValues.find(aValue);
+ if (aItr != theValues.end()) {
+ aItr->second += 1;
+ //aNbConnects = nb;
+ }
+ else {
+ theValues[aValue] = 1;
+ //aNbConnects = 1;
+ }
+ //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
+ }
+
+}
+
+/*
+ Class : BallDiameter
+ Description : Functor returning diameter of a ball element
+*/
+double BallDiameter::GetValue( long theId )
+{
+ double diameter = 0;
+
+ if ( const SMDS_BallElement* ball =
+ dynamic_cast<const SMDS_BallElement*>( myMesh->FindElement( theId )))
+ {
+ diameter = ball->GetDiameter();
+ }
+ return diameter;
+}
+
+double BallDiameter::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // meaningless as it is not a quality control functor
+ return Value;
+}
+
+SMDSAbs_ElementType BallDiameter::GetType() const
+{
+ return SMDSAbs_Ball;
+}
+
+
+/*
+ PREDICATES
+*/
+
+/*
+ Class : BadOrientedVolume
+ Description : Predicate bad oriented volumes
+*/
+
+BadOrientedVolume::BadOrientedVolume()
+{
+ myMesh = 0;
+}
+
+void BadOrientedVolume::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMesh = theMesh;
+}
+
+bool BadOrientedVolume::IsSatisfy( long theId )
+{
+ if ( myMesh == 0 )
+ return false;
+
+ SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
+ return !vTool.IsForward();
+}
+
+SMDSAbs_ElementType BadOrientedVolume::GetType() const
+{
+ return SMDSAbs_Volume;
+}
+
+/*
+ Class : BareBorderVolume
+*/
+
+bool BareBorderVolume::IsSatisfy(long theElementId )
+{
+ SMDS_VolumeTool myTool;
+ if ( myTool.Set( myMesh->FindElement(theElementId)))
+ {
+ for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
+ if ( myTool.IsFreeFace( iF ))
+ {
+ const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
+ vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
+ if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
+ return true;
+ }
+ }
+ return false;
+}
+
+/*
+ Class : BareBorderFace
+*/
+
+bool BareBorderFace::IsSatisfy(long theElementId )
+{
+ bool ok = false;
+ if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
+ {
+ if ( face->GetType() == SMDSAbs_Face )
+ {
+ int nbN = face->NbCornerNodes();
+ for ( int i = 0; i < nbN && !ok; ++i )
+ {
+ // check if a link is shared by another face
+ const SMDS_MeshNode* n1 = face->GetNode( i );
+ const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
+ SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
+ bool isShared = false;
+ while ( !isShared && fIt->more() )
+ {
+ const SMDS_MeshElement* f = fIt->next();
+ isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
+ }
+ if ( !isShared )
+ {
+ const int iQuad = face->IsQuadratic();
+ myLinkNodes.resize( 2 + iQuad);
+ myLinkNodes[0] = n1;
+ myLinkNodes[1] = n2;
+ if ( iQuad )
+ myLinkNodes[2] = face->GetNode( i+nbN );
+ ok = !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
+ }
+ }
+ }
+ }
+ return ok;
+}
+
+/*
+ Class : OverConstrainedVolume
+*/
+
+bool OverConstrainedVolume::IsSatisfy(long theElementId )
+{
+ // An element is over-constrained if it has N-1 free borders where
+ // N is the number of edges/faces for a 2D/3D element.
+ SMDS_VolumeTool myTool;
+ if ( myTool.Set( myMesh->FindElement(theElementId)))
+ {
+ int nbSharedFaces = 0;
+ for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
+ if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
+ break;
+ return ( nbSharedFaces == 1 );
+ }
+ return false;
+}
+
+/*
+ Class : OverConstrainedFace
+*/
+
+bool OverConstrainedFace::IsSatisfy(long theElementId )
+{
+ // An element is over-constrained if it has N-1 free borders where
+ // N is the number of edges/faces for a 2D/3D element.
+ if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
+ if ( face->GetType() == SMDSAbs_Face )
+ {
+ int nbSharedBorders = 0;
+ int nbN = face->NbCornerNodes();
+ for ( int i = 0; i < nbN; ++i )
+ {
+ // check if a link is shared by another face
+ const SMDS_MeshNode* n1 = face->GetNode( i );
+ const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
+ SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
+ bool isShared = false;
+ while ( !isShared && fIt->more() )
+ {
+ const SMDS_MeshElement* f = fIt->next();
+ isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
+ }
+ if ( isShared && ++nbSharedBorders > 1 )
+ break;
+ }
+ return ( nbSharedBorders == 1 );
+ }
+ return false;
+}
+
+/*
+ Class : CoincidentNodes
+ Description : Predicate of Coincident nodes
+*/
+
+CoincidentNodes::CoincidentNodes()
+{
+ myToler = 1e-5;
+}
+
+bool CoincidentNodes::IsSatisfy( long theElementId )
+{
+ return myCoincidentIDs.Contains( theElementId );
+}
+
+SMDSAbs_ElementType CoincidentNodes::GetType() const
+{
+ return SMDSAbs_Node;
+}
+
+void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMeshModifTracer.SetMesh( theMesh );
+ if ( myMeshModifTracer.IsMeshModified() )
+ {
+ TIDSortedNodeSet nodesToCheck;
+ SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
+ while ( nIt->more() )
+ nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
+
+ list< list< const SMDS_MeshNode*> > nodeGroups;
+ SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
+
+ myCoincidentIDs.Clear();
+ list< list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
+ for ( ; groupIt != nodeGroups.end(); ++groupIt )
+ {
+ list< const SMDS_MeshNode*>& coincNodes = *groupIt;
+ list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
+ for ( ; n != coincNodes.end(); ++n )
+ myCoincidentIDs.Add( (*n)->GetID() );
+ }
+ }
+}
+
+/*
+ Class : CoincidentElements
+ Description : Predicate of Coincident Elements
+ Note : This class is suitable only for visualization of Coincident Elements
+*/
+
+CoincidentElements::CoincidentElements()
+{
+ myMesh = 0;
+}
+
+void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMesh = theMesh;
+}
+
+bool CoincidentElements::IsSatisfy( long theElementId )
+{
+ if ( !myMesh ) return false;
+
+ if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
+ {
+ if ( e->GetType() != GetType() ) return false;
+ set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
+ const int nbNodes = e->NbNodes();
+ SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e2 = invIt->next();
+ if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
+
+ bool sameNodes = true;
+ for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
+ sameNodes = ( elemNodes.count( e2->GetNode( i )));
+ if ( sameNodes )
+ return true;
+ }
+ }
+ return false;
+}
+
+SMDSAbs_ElementType CoincidentElements1D::GetType() const
+{
+ return SMDSAbs_Edge;
+}
+SMDSAbs_ElementType CoincidentElements2D::GetType() const
+{
+ return SMDSAbs_Face;
+}
+SMDSAbs_ElementType CoincidentElements3D::GetType() const
+{
+ return SMDSAbs_Volume;
+}
+
+
+/*
+ Class : FreeBorders
+ Description : Predicate for free borders
+*/
+
+FreeBorders::FreeBorders()
+{
+ myMesh = 0;
+}
+
+void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMesh = theMesh;
+}
+
+bool FreeBorders::IsSatisfy( long theId )
+{
+ return getNbMultiConnection( myMesh, theId ) == 1;