X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=5238b982f3c55ee99f83bf4048d2096b2f866e3c;hb=d8f644ca3d4ce62f2ef41d4aacb52f5bb1221df3;hp=7edaa0e2169f315eea09cf9f44b7583cb1026821;hpb=57b43b4d010e2d0a1529d3c131bbb9d416e63258;p=modules%2Fsmesh.git diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 7edaa0e21..5238b982f 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -1,4 +1,6 @@ -// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS // // This library is free software; you can redistribute it and/or @@ -15,31 +17,43 @@ // License along with this library; if not, write to the Free Software // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// #include "SMESH_ControlsDef.hxx" #include +#include +#include +#include #include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + #include #include #include -#include #include +#include #include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include #include "SMDS_Mesh.hxx" #include "SMDS_Iterator.hxx" @@ -49,12 +63,20 @@ #include "SMDS_QuadraticFaceOfNodes.hxx" #include "SMDS_QuadraticEdge.hxx" +#include "SMESHDS_Mesh.hxx" +#include "SMESHDS_GroupBase.hxx" /* AUXILIARY METHODS */ namespace{ + + inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode ) + { + return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() ); + } + inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) { gp_Vec v1( P1 - P2 ), v2( P3 - P2 ); @@ -128,34 +150,54 @@ namespace{ } } } - int aResult = max ( aResult0, aResult1 ); + int aResult = std::max ( aResult0, aResult1 ); // TColStd_MapOfInteger aMap; // SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator(); // if ( anIter != 0 ) { // while( anIter->more() ) { -// const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); -// if ( aNode == 0 ) -// return 0; -// SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); -// while( anElemIter->more() ) { -// const SMDS_MeshElement* anElem = anElemIter->next(); -// if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) { -// int anId = anElem->GetID(); - -// if ( anIter->more() ) // i.e. first node -// aMap.Add( anId ); -// else if ( aMap.Contains( anId ) ) -// aResult++; -// } -// } +// const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); +// if ( aNode == 0 ) +// return 0; +// SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); +// while( anElemIter->more() ) { +// const SMDS_MeshElement* anElem = anElemIter->next(); +// if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) { +// int anId = anElem->GetID(); + +// if ( anIter->more() ) // i.e. first node +// aMap.Add( anId ); +// else if ( aMap.Contains( anId ) ) +// aResult++; +// } +// } // } // } return aResult; } + gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 ) + { + int aNbNode = theFace->NbNodes(); + + gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0)); + gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0)); + gp_XYZ n = q1 ^ q2; + if ( aNbNode > 3 ) { + gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0)); + n += q2 ^ q3; + } + double len = n.Modulus(); + bool zeroLen = ( len <= numeric_limits::min()); + if ( !zeroLen ) + n /= len; + + if (ok) *ok = !zeroLen; + + return n; + } } @@ -163,8 +205,8 @@ namespace{ using namespace SMESH::Controls; /* - FUNCTORS -*/ + * FUNCTORS + */ /* Class : NumericalFunctor @@ -262,6 +304,63 @@ double NumericalFunctor::GetValue( long theId ) return 0.; } +//================================================================================ +/*! + * \brief Return histogram of functor values + * \param nbIntervals - number of intervals + * \param nbEvents - number of mesh elements having values within i-th interval + * \param funValues - boundaries of intervals + */ +//================================================================================ + +void NumericalFunctor::GetHistogram(int nbIntervals, + std::vector& nbEvents, + std::vector& funValues) +{ + if ( nbIntervals < 1 || + !myMesh || + !myMesh->GetMeshInfo().NbElements( GetType() )) + return; + nbEvents.resize( nbIntervals, 0 ); + funValues.resize( nbIntervals+1 ); + + // get all values sorted + std::multiset< double > values; + SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType()); + while ( elemIt->more() ) + values.insert( GetValue( elemIt->next()->GetID() )); + + // case nbIntervals == 1 + funValues[0] = *values.begin(); + funValues[nbIntervals] = *values.rbegin(); + if ( nbIntervals == 1 ) + { + nbEvents[0] = values.size(); + return; + } + // case of 1 value + if (funValues.front() == funValues.back()) + { + nbEvents.resize( 1 ); + nbEvents[0] = values.size(); + funValues[1] = funValues.back(); + funValues.resize( 2 ); + } + // generic case + std::multiset< double >::iterator min = values.begin(), max; + for ( int i = 0; i < nbIntervals; ++i ) + { + double r = (i+1) / double( nbIntervals ); + funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r; + if ( min != values.end() && *min <= funValues[i+1] ) + { + max = values.upper_bound( funValues[i+1] ); // greater than funValues[i+1], or end() + nbEvents[i] = std::distance( min, max ); + min = max; + } + } +} + //======================================================================= //function : GetValue //purpose : @@ -298,6 +397,246 @@ SMDSAbs_ElementType Volume::GetType() const } +/* + Class : MaxElementLength2D + Description : Functor calculating maximum length of 2D element +*/ + +double MaxElementLength2D::GetValue( long theElementId ) +{ + TSequenceOfXYZ P; + if( GetPoints( theElementId, P ) ) { + double aVal = 0; + const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId ); + SMDSAbs_ElementType aType = aElem->GetType(); + int len = P.size(); + switch( aType ) { + 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 )); + double D1 = getDistance(P( 1 ),P( 3 )); + double D2 = getDistance(P( 2 ),P( 4 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2)); + break; + } + else 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)); + 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 )); + double D1 = getDistance(P( 1 ),P( 5 )); + double D2 = getDistance(P( 3 ),P( 7 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2)); + break; + } + } + + if( myPrecision >= 0 ) + { + double prec = pow( 10., (double)myPrecision ); + aVal = floor( aVal * prec + 0.5 ) / prec; + } + return aVal; + } + return 0.; +} + +double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; +} + +SMDSAbs_ElementType MaxElementLength2D::GetType() const +{ + return SMDSAbs_Face; +} + +/* + Class : MaxElementLength3D + Description : Functor calculating maximum length of 3D element +*/ + +double MaxElementLength3D::GetValue( long theElementId ) +{ + TSequenceOfXYZ P; + if( GetPoints( theElementId, P ) ) { + double aVal = 0; + const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId ); + SMDSAbs_ElementType aType = aElem->GetType(); + int len = P.size(); + switch( aType ) { + case SMDSAbs_Volume: + if( len == 4 ) { // tetras + 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 ) { // pyramids + 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 ) { // pentas + 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 ) { // hexas + 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 )); + double D1 = getDistance(P( 1 ),P( 7 )); + double D2 = getDistance(P( 2 ),P( 8 )); + double D3 = getDistance(P( 3 ),P( 5 )); + double D4 = getDistance(P( 4 ),P( 6 )); + 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)); + aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); + break; + } + else if( len == 10 ) { // quadratic tetras + 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 pyramids + 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 pentas + 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 hexas + 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 )); + double D1 = getDistance(P( 1 ),P( 7 )); + double D2 = getDistance(P( 2 ),P( 8 )); + double D3 = getDistance(P( 3 ),P( 5 )); + double D4 = getDistance(P( 4 ),P( 6 )); + 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)); + aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); + break; + } + else if( len > 1 && aElem->IsPoly() ) { // polys + // get the maximum distance between all pairs of nodes + for( int i = 1; i <= len; i++ ) { + for( int j = 1; j <= len; j++ ) { + if( j > i ) { // optimization of the loop + double D = getDistance( P(i), P(j) ); + aVal = Max( aVal, D ); + } + } + } + } + } + + if( myPrecision >= 0 ) + { + double prec = pow( 10., (double)myPrecision ); + aVal = floor( aVal * prec + 0.5 ) / prec; + } + return aVal; + } + return 0.; +} + +double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; +} + +SMDSAbs_ElementType MaxElementLength3D::GetType() const +{ + return SMDSAbs_Volume; +} + + /* Class : MinimumAngle Description : Functor for calculation of minimum angle @@ -350,56 +689,137 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) if ( nbNodes < 3 ) return 0; - // Compute lengths of the sides - - vector< double > aLen (nbNodes); - - for ( int i = 0; i < nbNodes - 1; i++ ) - aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) ); - aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) ); - // Compute aspect ratio - if ( nbNodes == 3 ) - { + if ( nbNodes == 3 ) { + // Compute lengths of the sides + std::vector< double > aLen (nbNodes); + for ( int i = 0; i < nbNodes - 1; i++ ) + aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) ); + aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) ); // Q = alfa * h * p / S, where // // alfa = sqrt( 3 ) / 6 // h - length of the longest edge // p - half perimeter // S - triangle surface - const double alfa = sqrt( 3. ) / 6.; double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) ); double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.; double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) ); if ( anArea <= Precision::Confusion() ) return 0.; - return alfa * maxLen * half_perimeter / anArea; } - else - { - // return aspect ratio of the worst triange which can be built + else if ( nbNodes == 6 ) { // quadratic triangles + // Compute lengths of the sides + std::vector< double > aLen (3); + aLen[0] = getDistance( P(1), P(3) ); + aLen[1] = getDistance( P(3), P(5) ); + aLen[2] = getDistance( P(5), P(1) ); + // Q = alfa * h * p / S, where + // + // alfa = sqrt( 3 ) / 6 + // h - length of the longest edge + // p - half perimeter + // S - triangle surface + const double alfa = sqrt( 3. ) / 6.; + double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) ); + double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.; + double anArea = getArea( P(1), P(3), P(5) ); + if ( anArea <= Precision::Confusion() ) + return 0.; + return alfa * maxLen * half_perimeter / anArea; + } + else if( nbNodes == 4 ) { // quadrangle + // Compute lengths of the sides + std::vector< double > aLen (4); + aLen[0] = getDistance( P(1), P(2) ); + aLen[1] = getDistance( P(2), P(3) ); + aLen[2] = getDistance( P(3), P(4) ); + aLen[3] = getDistance( P(4), P(1) ); + // Compute lengths of the diagonals + std::vector< double > aDia (2); + aDia[0] = getDistance( P(1), P(3) ); + aDia[1] = getDistance( P(2), P(4) ); + // Compute areas of all triangles which can be built + // taking three nodes of the quadrangle + std::vector< double > anArea (4); + anArea[0] = getArea( P(1), P(2), P(3) ); + anArea[1] = getArea( P(1), P(2), P(4) ); + anArea[2] = getArea( P(1), P(3), P(4) ); + anArea[3] = getArea( P(2), P(3), P(4) ); + // Q = alpha * L * C1 / C2, where + // + // alpha = sqrt( 1/32 ) + // L = max( L1, L2, L3, L4, D1, D2 ) + // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 ) + // C2 = min( S1, S2, S3, S4 ) + // Li - lengths of the edges + // Di - lengths of the diagonals + // Si - areas of the triangles + const double alpha = sqrt( 1 / 32. ); + double L = Max( aLen[ 0 ], + Max( aLen[ 1 ], + Max( aLen[ 2 ], + Max( aLen[ 3 ], + Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) ); + double C1 = sqrt( ( aLen[0] * aLen[0] + + aLen[1] * aLen[1] + + aLen[2] * aLen[2] + + aLen[3] * aLen[3] ) / 4. ); + double C2 = Min( anArea[ 0 ], + Min( anArea[ 1 ], + Min( anArea[ 2 ], anArea[ 3 ] ) ) ); + if ( C2 <= Precision::Confusion() ) + return 0.; + return alpha * L * C1 / C2; + } + else if( nbNodes == 8 ){ // nbNodes==8 - quadratic quadrangle + // Compute lengths of the sides + std::vector< double > aLen (4); + aLen[0] = getDistance( P(1), P(3) ); + aLen[1] = getDistance( P(3), P(5) ); + aLen[2] = getDistance( P(5), P(7) ); + aLen[3] = getDistance( P(7), P(1) ); + // Compute lengths of the diagonals + std::vector< double > aDia (2); + aDia[0] = getDistance( P(1), P(5) ); + aDia[1] = getDistance( P(3), P(7) ); + // Compute areas of all triangles which can be built // taking three nodes of the quadrangle - TSequenceOfXYZ triaPnts(3); - // triangle on nodes 1 3 2 - triaPnts(1) = P(1); - triaPnts(2) = P(3); - triaPnts(3) = P(2); - double ar = GetValue( triaPnts ); - // triangle on nodes 1 3 4 - triaPnts(3) = P(4); - ar = Max ( ar, GetValue( triaPnts )); - // triangle on nodes 1 2 4 - triaPnts(2) = P(2); - ar = Max ( ar, GetValue( triaPnts )); - // triangle on nodes 3 2 4 - triaPnts(1) = P(3); - ar = Max ( ar, GetValue( triaPnts )); - - return ar; + std::vector< double > anArea (4); + anArea[0] = getArea( P(1), P(3), P(5) ); + anArea[1] = getArea( P(1), P(3), P(7) ); + anArea[2] = getArea( P(1), P(5), P(7) ); + anArea[3] = getArea( P(3), P(5), P(7) ); + // Q = alpha * L * C1 / C2, where + // + // alpha = sqrt( 1/32 ) + // L = max( L1, L2, L3, L4, D1, D2 ) + // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 ) + // C2 = min( S1, S2, S3, S4 ) + // Li - lengths of the edges + // Di - lengths of the diagonals + // Si - areas of the triangles + const double alpha = sqrt( 1 / 32. ); + double L = Max( aLen[ 0 ], + Max( aLen[ 1 ], + Max( aLen[ 2 ], + Max( aLen[ 3 ], + Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) ); + double C1 = sqrt( ( aLen[0] * aLen[0] + + aLen[1] * aLen[1] + + aLen[2] * aLen[2] + + aLen[3] * aLen[3] ) / 4. ); + double C2 = Min( anArea[ 0 ], + Min( anArea[ 1 ], + Min( anArea[ 2 ], anArea[ 3 ] ) ) ); + if ( C2 <= Precision::Confusion() ) + return 0.; + return alpha * L * C1 / C2; } + return 0; } double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -428,9 +848,9 @@ namespace{ inline double getArea(double theHalfPerim, double theTria[3]){ return sqrt(theHalfPerim* - (theHalfPerim-theTria[0])* - (theHalfPerim-theTria[1])* - (theHalfPerim-theTria[2])); + (theHalfPerim-theTria[0])* + (theHalfPerim-theTria[1])* + (theHalfPerim-theTria[2])); } inline double getVolume(double theLen[6]){ @@ -472,11 +892,11 @@ namespace{ inline double getMaxHeight(double theLen[6]) { - double aHeight = max(theLen[0],theLen[1]); - aHeight = max(aHeight,theLen[2]); - aHeight = max(aHeight,theLen[3]); - aHeight = max(aHeight,theLen[4]); - aHeight = max(aHeight,theLen[5]); + 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; } @@ -486,7 +906,17 @@ 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 return aQuality; + } + switch(nbNodes){ case 4:{ double aLen[6] = { @@ -520,187 +950,188 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) //double aVolume = getVolume(aLen); double aHeight = getMaxHeight(aLen); static double aCoeff = sqrt(2.0)/12.0; - aQuality = aCoeff*aHeight*aSumArea/aVolume; + 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 = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + 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 = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + 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 = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } break; } @@ -717,7 +1148,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) 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 = max( aQuality, aspect2D.GetValue( points )); + aQuality = std::max( aQuality, aspect2D.GetValue( points )); } } return aQuality; @@ -746,7 +1177,7 @@ double Warping::GetValue( const TSequenceOfXYZ& P ) if ( P.size() != 4 ) return 0; - gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4; + 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 ); @@ -777,7 +1208,7 @@ double Warping::ComputeA( const gp_XYZ& thePnt1, N.Normalize(); double H = ( thePnt2 - theG ).Dot( N ); - return asin( fabs( H / L ) ) * 180 / PI; + return asin( fabs( H / L ) ) * 180. / PI; } double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -801,13 +1232,13 @@ SMDSAbs_ElementType Warping::GetType() const double Taper::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 4 ) - return 0; + 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 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() ) @@ -841,42 +1272,46 @@ SMDSAbs_ElementType Taper::GetType() const */ 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_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 ); + 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; + return 0.; // Compute skew - static double PI2 = PI / 2; + static double PI2 = 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 / PI; + return Max( A0, Max( A1, A2 ) ) * 180. / 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_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 ) ); + ? 0. : fabs( PI2 - v1.Angle( v2 ) ); - return A * 180 / PI; + //BUG SWP12743 + if ( A < Precision::Angular() ) + return 0.; + + return A * 180. / PI; } } @@ -900,16 +1335,20 @@ SMDSAbs_ElementType Skew::GetType() const */ double Area::GetValue( const TSequenceOfXYZ& P ) { - 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); + 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 SumVec.Magnitude() * 0.5; + return val; } double Area::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -973,160 +1412,160 @@ double Length2D::GetValue( long theElementId) case SMDSAbs_Node: case SMDSAbs_Edge: if (len == 2){ - aVal = getDistance( P( 1 ), P( 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 )); + 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; + 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; + 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)); + 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="<FindNode( theNodeId ); + if (!aNode) + return false; + + return (aNode->NbInverseElements() < 1); } -//======================================================================= -// name : GetRangeStr -// Purpose : Get range as a string. -// Example: "1,2,3,50-60,63,67,70-" -//======================================================================= -void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr ) +SMDSAbs_ElementType FreeNodes::GetType() const { - theResStr.Clear(); + return SMDSAbs_Node; +} - TColStd_SequenceOfInteger anIntSeq; - TColStd_SequenceOfAsciiString aStrSeq; - TColStd_MapIteratorOfMapOfInteger anIter( myIds ); - for ( ; anIter.More(); anIter.Next() ) - { - int anId = anIter.Key(); - TCollection_AsciiString aStr( anId ); - anIntSeq.Append( anId ); - aStrSeq.Append( aStr ); - } +/* + Class : FreeFaces + Description : Predicate for free faces +*/ - for ( int i = 1, n = myMin.Length(); i <= n; i++ ) - { - int aMinId = myMin( i ); - int aMaxId = myMax( i ); +FreeFaces::FreeFaces() +{ + myMesh = 0; +} - TCollection_AsciiString aStr; - if ( aMinId != IntegerFirst() ) - aStr += aMinId; +void FreeFaces::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; +} - aStr += "-"; +bool FreeFaces::IsSatisfy( long theId ) +{ + if (!myMesh) return false; + // check that faces nodes refers to less than two common volumes + const SMDS_MeshElement* aFace = myMesh->FindElement( theId ); + if ( !aFace || aFace->GetType() != SMDSAbs_Face ) + return false; - if ( aMaxId != IntegerLast() ) - aStr += aMaxId; + int nbNode = aFace->NbNodes(); - // find position of the string in result sequence and insert string in it - if ( anIntSeq.Length() == 0 ) - { - anIntSeq.Append( aMinId ); - aStrSeq.Append( aStr ); - } - else - { - if ( aMinId < anIntSeq.First() ) + // collect volumes check that number of volumss with count equal nbNode not less than 2 + typedef map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters + typedef map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator + TMapOfVolume mapOfVol; + + SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator(); + while ( nodeItr->more() ) { + const SMDS_MeshNode* aNode = static_cast(nodeItr->next()); + if ( !aNode ) continue; + SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume); + while ( volItr->more() ) { + SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next(); + TItrMapOfVolume itr = mapOfVol.insert(make_pair(aVol, 0)).first; + (*itr).second++; + } + } + int nbVol = 0; + TItrMapOfVolume volItr = mapOfVol.begin(); + TItrMapOfVolume volEnd = mapOfVol.end(); + for ( ; volItr != volEnd; ++volItr ) + if ( (*volItr).second >= nbNode ) + nbVol++; + // face is not free if number of volumes constructed on thier nodes more than one + return (nbVol < 2); +} + +SMDSAbs_ElementType FreeFaces::GetType() const +{ + return SMDSAbs_Face; +} + +/* + Class : LinearOrQuadratic + Description : Predicate to verify whether a mesh element is linear +*/ + +LinearOrQuadratic::LinearOrQuadratic() +{ + myMesh = 0; +} + +void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; +} + +bool LinearOrQuadratic::IsSatisfy( long theId ) +{ + if (!myMesh) return false; + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) ) + return false; + return (!anElem->IsQuadratic()); +} + +void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType ) +{ + myType = theType; +} + +SMDSAbs_ElementType LinearOrQuadratic::GetType() const +{ + return myType; +} + +/* + Class : GroupColor + Description : Functor for check color of group to whic mesh element belongs to +*/ + +GroupColor::GroupColor() +{ +} + +bool GroupColor::IsSatisfy( long theId ) +{ + return (myIDs.find( theId ) != myIDs.end()); +} + +void GroupColor::SetType( SMDSAbs_ElementType theType ) +{ + myType = theType; +} + +SMDSAbs_ElementType GroupColor::GetType() const +{ + return myType; +} + +static bool isEqual( const Quantity_Color& theColor1, + const Quantity_Color& theColor2 ) +{ + // tolerance to compare colors + const double tol = 5*1e-3; + return ( fabs( theColor1.Red() - theColor2.Red() ) < tol && + fabs( theColor1.Green() - theColor2.Green() ) < tol && + fabs( theColor1.Blue() - theColor2.Blue() ) < tol ); +} + + +void GroupColor::SetMesh( const SMDS_Mesh* theMesh ) +{ + myIDs.clear(); + + const SMESHDS_Mesh* aMesh = dynamic_cast(theMesh); + if ( !aMesh ) + return; + + int nbGrp = aMesh->GetNbGroups(); + if ( !nbGrp ) + return; + + // iterates on groups and find necessary elements ids + const std::set& aGroups = aMesh->GetGroups(); + set::const_iterator GrIt = aGroups.begin(); + for (; GrIt != aGroups.end(); GrIt++) { + SMESHDS_GroupBase* aGrp = (*GrIt); + if ( !aGrp ) + continue; + // check type and color of group + if ( !isEqual( myColor, aGrp->GetColor() ) ) + continue; + if ( myType != SMDSAbs_All && myType != (SMDSAbs_ElementType)aGrp->GetType() ) + continue; + + SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType(); + if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) { + // add elements IDS into control + int aSize = aGrp->Extent(); + for (int i = 0; i < aSize; i++) + myIDs.insert( aGrp->GetID(i+1) ); + } + } +} + +void GroupColor::SetColorStr( const TCollection_AsciiString& theStr ) +{ + TCollection_AsciiString aStr = theStr; + aStr.RemoveAll( ' ' ); + aStr.RemoveAll( '\t' ); + for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) ) + aStr.Remove( aPos, 2 ); + Standard_Real clr[3]; + clr[0] = clr[1] = clr[2] = 0.; + for ( int i = 0; i < 3; i++ ) { + TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 ); + if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() ) + clr[i] = tmpStr.RealValue(); + } + myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB ); +} + +//======================================================================= +// name : GetRangeStr +// Purpose : Get range as a string. +// Example: "1,2,3,50-60,63,67,70-" +//======================================================================= +void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const +{ + theResStr.Clear(); + theResStr += TCollection_AsciiString( myColor.Red() ); + theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() ); + theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() ); +} + +/* + Class : ElemGeomType + Description : Predicate to check element geometry type +*/ + +ElemGeomType::ElemGeomType() +{ + myMesh = 0; + myType = SMDSAbs_All; + myGeomType = SMDSGeom_TRIANGLE; +} + +void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; +} + +bool ElemGeomType::IsSatisfy( long theId ) +{ + if (!myMesh) return false; + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + if ( !anElem ) + return false; + const SMDSAbs_ElementType anElemType = anElem->GetType(); + if ( myType != SMDSAbs_All && anElemType != myType ) + return false; + const int aNbNode = anElem->NbNodes(); + bool isOk = false; + switch( anElemType ) + { + case SMDSAbs_Node: + isOk = (myGeomType == SMDSGeom_POINT); + break; + + case SMDSAbs_Edge: + isOk = (myGeomType == SMDSGeom_EDGE); + break; + + case SMDSAbs_Face: + if ( myGeomType == SMDSGeom_TRIANGLE ) + isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 6 : aNbNode == 3)); + else if ( myGeomType == SMDSGeom_QUADRANGLE ) + isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 8 : aNbNode == 4)); + else if ( myGeomType == SMDSGeom_POLYGON ) + isOk = anElem->IsPoly(); + break; + + case SMDSAbs_Volume: + if ( myGeomType == SMDSGeom_TETRA ) + isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 10 : aNbNode == 4)); + else if ( myGeomType == SMDSGeom_PYRAMID ) + isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 13 : aNbNode == 5)); + else if ( myGeomType == SMDSGeom_PENTA ) + isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 15 : aNbNode == 6)); + else if ( myGeomType == SMDSGeom_HEXA ) + isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 20 : aNbNode == 8)); + else if ( myGeomType == SMDSGeom_POLYHEDRA ) + isOk = anElem->IsPoly(); + break; + default: break; + } + return isOk; +} + +void ElemGeomType::SetType( SMDSAbs_ElementType theType ) +{ + myType = theType; +} + +SMDSAbs_ElementType ElemGeomType::GetType() const +{ + return myType; +} + +void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType ) +{ + myGeomType = theType; +} + +SMDSAbs_GeometryType ElemGeomType::GetGeomType() const +{ + return myGeomType; +} + +//================================================================================ +/*! + * \brief Class CoplanarFaces + */ +//================================================================================ + +CoplanarFaces::CoplanarFaces() + : myMesh(0), myFaceID(0), myToler(0) +{ +} +bool CoplanarFaces::IsSatisfy( long theElementId ) +{ + if ( myCoplanarIDs.empty() ) + { + // Build a set of coplanar face ids + + if ( !myMesh || !myFaceID || !myToler ) + return false; + + const SMDS_MeshElement* face = myMesh->FindElement( myFaceID ); + if ( !face || face->GetType() != SMDSAbs_Face ) + return false; + + bool normOK; + gp_Vec myNorm = getNormale( static_cast(face), &normOK ); + if (!normOK) + return false; + + const double radianTol = myToler * PI180; + typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TFaceIt; + std::set checkedFaces, checkedNodes; + std::list faceQueue( 1, face ); + while ( !faceQueue.empty() ) + { + face = faceQueue.front(); + if ( checkedFaces.insert( face ).second ) + { + gp_Vec norm = getNormale( static_cast(face), &normOK ); + if (!normOK || myNorm.Angle( norm ) <= radianTol) + { + myCoplanarIDs.insert( face->GetID() ); + std::set neighborFaces; + for ( int i = 0; i < face->NbCornerNodes(); ++i ) + { + const SMDS_MeshNode* n = face->GetNode( i ); + if ( checkedNodes.insert( n ).second ) + neighborFaces.insert( TFaceIt( n->GetInverseElementIterator(SMDSAbs_Face)), + TFaceIt()); + } + faceQueue.insert( faceQueue.end(), neighborFaces.begin(), neighborFaces.end() ); + } + } + faceQueue.pop_front(); + } + } + return myCoplanarIDs.count( theElementId ); +} + +/* + *Class : RangeOfIds + *Description : Predicate for Range of Ids. + * Range may be specified with two ways. + * 1. Using AddToRange method + * 2. With SetRangeStr method. Parameter of this method is a string + * like as "1,2,3,50-60,63,67,70-" +*/ + +//======================================================================= +// name : RangeOfIds +// Purpose : Constructor +//======================================================================= +RangeOfIds::RangeOfIds() +{ + myMesh = 0; + myType = SMDSAbs_All; +} + +//======================================================================= +// name : SetMesh +// Purpose : Set mesh +//======================================================================= +void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; +} + +//======================================================================= +// name : AddToRange +// Purpose : Add ID to the range +//======================================================================= +bool RangeOfIds::AddToRange( long theEntityId ) +{ + myIds.Add( theEntityId ); + return true; +} + +//======================================================================= +// name : GetRangeStr +// Purpose : Get range as a string. +// Example: "1,2,3,50-60,63,67,70-" +//======================================================================= +void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr ) +{ + theResStr.Clear(); + + TColStd_SequenceOfInteger anIntSeq; + TColStd_SequenceOfAsciiString aStrSeq; + + TColStd_MapIteratorOfMapOfInteger anIter( myIds ); + for ( ; anIter.More(); anIter.Next() ) + { + int anId = anIter.Key(); + TCollection_AsciiString aStr( anId ); + anIntSeq.Append( anId ); + aStrSeq.Append( aStr ); + } + + for ( int i = 1, n = myMin.Length(); i <= n; i++ ) + { + int aMinId = myMin( i ); + int aMaxId = myMax( i ); + + TCollection_AsciiString aStr; + if ( aMinId != IntegerFirst() ) + aStr += aMinId; + + aStr += "-"; + + if ( aMaxId != IntegerLast() ) + aStr += aMaxId; + + // find position of the string in result sequence and insert string in it + if ( anIntSeq.Length() == 0 ) + { + anIntSeq.Append( aMinId ); + aStrSeq.Append( aStr ); + } + else + { + if ( aMinId < anIntSeq.First() ) { anIntSeq.Prepend( aMinId ); aStrSeq.Prepend( aStr ); @@ -2046,15 +2848,15 @@ void Filter::SetPredicate( PredicatePtr thePredicate ) template inline void FillSequence(const TIterator& theIterator, - TPredicate& thePredicate, - Filter::TIdSequence& theSequence) + TPredicate& thePredicate, + Filter::TIdSequence& theSequence) { if ( theIterator ) { while( theIterator->more() ) { TElement anElem = theIterator->next(); long anId = anElem->GetID(); if ( thePredicate->IsSatisfy( anId ) ) - theSequence.push_back( anId ); + theSequence.push_back( anId ); } } } @@ -2062,8 +2864,8 @@ inline void FillSequence(const TIterator& theIterator, void Filter:: GetElementsId( const SMDS_Mesh* theMesh, - PredicatePtr thePredicate, - TIdSequence& theSequence ) + PredicatePtr thePredicate, + TIdSequence& theSequence ) { theSequence.clear(); @@ -2096,7 +2898,7 @@ GetElementsId( const SMDS_Mesh* theMesh, void Filter::GetElementsId( const SMDS_Mesh* theMesh, - Filter::TIdSequence& theSequence ) + Filter::TIdSequence& theSequence ) { GetElementsId(theMesh,myPredicate,theSequence); } @@ -2269,32 +3071,6 @@ static void getLinks( const SMDS_MeshFace* theFace, } } -static gp_XYZ getNormale( const SMDS_MeshFace* theFace ) -{ - gp_XYZ n; - int aNbNode = theFace->NbNodes(); - TColgp_Array1OfXYZ anArrOfXYZ(1,4); - SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator(); - int i = 1; - for ( ; aNodeItr->more() && i <= 4; i++ ) { - SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); - anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); - } - - gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1); - gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1); - n = q1 ^ q2; - if ( aNbNode > 3 ) { - gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1); - n += q2 ^ q3; - } - double len = n.Modulus(); - if ( len > 0 ) - n /= len; - - return n; -} - bool ManifoldPart::findConnected ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt, SMDS_MeshFace* theStartFace, @@ -2482,6 +3258,7 @@ ElementsOnSurface::ElementsOnSurface() myType = SMDSAbs_All; mySurf.Nullify(); myToler = Precision::Confusion(); + myUseBoundaries = false; } ElementsOnSurface::~ElementsOnSurface() @@ -2494,7 +3271,6 @@ void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh ) if ( myMesh == theMesh ) return; myMesh = theMesh; - myIds.Clear(); process(); } @@ -2507,25 +3283,41 @@ SMDSAbs_ElementType ElementsOnSurface::GetType() const { return myType; } void ElementsOnSurface::SetTolerance( const double theToler ) -{ myToler = theToler; } +{ + if ( myToler != theToler ) + myIds.Clear(); + myToler = theToler; +} double ElementsOnSurface::GetTolerance() const +{ return myToler; } + +void ElementsOnSurface::SetUseBoundaries( bool theUse ) { - return myToler; + if ( myUseBoundaries != theUse ) { + myUseBoundaries = theUse; + SetSurface( mySurf, myType ); + } } void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape, const SMDSAbs_ElementType theType ) { + myIds.Clear(); myType = theType; mySurf.Nullify(); if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE ) - { - mySurf.Nullify(); return; - } - TopoDS_Face aFace = TopoDS::Face( theShape ); - mySurf = BRep_Tool::Surface( aFace ); + mySurf = TopoDS::Face( theShape ); + BRepAdaptor_Surface SA( mySurf, myUseBoundaries ); + Standard_Real + u1 = SA.FirstUParameter(), + u2 = SA.LastUParameter(), + v1 = SA.FirstVParameter(), + v2 = SA.LastVParameter(); + Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf ); + myProjector.Init( surf, u1,u2, v1,v2 ); + process(); } void ElementsOnSurface::process() @@ -2539,6 +3331,7 @@ void ElementsOnSurface::process() if ( myType == SMDSAbs_Face || myType == SMDSAbs_All ) { + myIds.ReSize( myMesh->NbFaces() ); SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); for(; anIter->more(); ) process( anIter->next() ); @@ -2546,6 +3339,7 @@ void ElementsOnSurface::process() if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All ) { + myIds.ReSize( myIds.Extent() + myMesh->NbEdges() ); SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator(); for(; anIter->more(); ) process( anIter->next() ); @@ -2553,6 +3347,7 @@ void ElementsOnSurface::process() if ( myType == SMDSAbs_Node ) { + myIds.ReSize( myMesh->NbNodes() ); SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator(); for(; anIter->more(); ) process( anIter->next() ); @@ -2576,32 +3371,347 @@ void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr ) myIds.Add( theElemPtr->GetID() ); } -bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) const +bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) { if ( mySurf.IsNull() ) return false; gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() ); - double aToler2 = myToler * myToler; - if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane))) + // double aToler2 = myToler * myToler; +// if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane))) +// { +// gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln(); +// if ( aPln.SquareDistance( aPnt ) > aToler2 ) +// return false; +// } +// else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface))) +// { +// gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder(); +// double aRad = aCyl.Radius(); +// gp_Ax3 anAxis = aCyl.Position(); +// gp_XYZ aLoc = aCyl.Location().XYZ(); +// double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc ); +// double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc ); +// if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 ) +// return false; +// } +// else +// return false; + myProjector.Perform( aPnt ); + bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler ); + + return isOn; +} + + +/* + ElementsOnShape +*/ + +ElementsOnShape::ElementsOnShape() + : myMesh(0), + myType(SMDSAbs_All), + myToler(Precision::Confusion()), + myAllNodesFlag(false) +{ + myCurShapeType = TopAbs_SHAPE; +} + +ElementsOnShape::~ElementsOnShape() +{ +} + +void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh) +{ + if (myMesh != theMesh) { + myMesh = theMesh; + SetShape(myShape, myType); + } +} + +bool ElementsOnShape::IsSatisfy (long theElementId) +{ + return myIds.Contains(theElementId); +} + +SMDSAbs_ElementType ElementsOnShape::GetType() const +{ + return myType; +} + +void ElementsOnShape::SetTolerance (const double theToler) +{ + if (myToler != theToler) { + myToler = theToler; + SetShape(myShape, myType); + } +} + +double ElementsOnShape::GetTolerance() const +{ + return myToler; +} + +void ElementsOnShape::SetAllNodes (bool theAllNodes) +{ + if (myAllNodesFlag != theAllNodes) { + myAllNodesFlag = theAllNodes; + SetShape(myShape, myType); + } +} + +void ElementsOnShape::SetShape (const TopoDS_Shape& theShape, + const SMDSAbs_ElementType theType) +{ + myType = theType; + myShape = theShape; + myIds.Clear(); + + if (myMesh == 0) return; + + switch (myType) { - gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln(); - if ( aPln.SquareDistance( aPnt ) > aToler2 ) - return false; + case SMDSAbs_All: + myIds.ReSize(myMesh->NbEdges() + myMesh->NbFaces() + myMesh->NbVolumes()); + break; + case SMDSAbs_Node: + myIds.ReSize(myMesh->NbNodes()); + break; + case SMDSAbs_Edge: + myIds.ReSize(myMesh->NbEdges()); + break; + case SMDSAbs_Face: + myIds.ReSize(myMesh->NbFaces()); + break; + case SMDSAbs_Volume: + myIds.ReSize(myMesh->NbVolumes()); + break; + default: + break; } - else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface))) + + myShapesMap.Clear(); + addShape(myShape); +} + +void ElementsOnShape::addShape (const TopoDS_Shape& theShape) +{ + if (theShape.IsNull() || myMesh == 0) + return; + + if (!myShapesMap.Add(theShape)) return; + + myCurShapeType = theShape.ShapeType(); + switch (myCurShapeType) { - gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder(); - double aRad = aCyl.Radius(); - gp_Ax3 anAxis = aCyl.Position(); - gp_XYZ aLoc = aCyl.Location().XYZ(); - double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc ); - double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc ); - if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 ) - return false; + case TopAbs_COMPOUND: + case TopAbs_COMPSOLID: + case TopAbs_SHELL: + case TopAbs_WIRE: + { + TopoDS_Iterator anIt (theShape, Standard_True, Standard_True); + for (; anIt.More(); anIt.Next()) addShape(anIt.Value()); + } + break; + case TopAbs_SOLID: + { + myCurSC.Load(theShape); + process(); + } + break; + case TopAbs_FACE: + { + TopoDS_Face aFace = TopoDS::Face(theShape); + BRepAdaptor_Surface SA (aFace, true); + Standard_Real + u1 = SA.FirstUParameter(), + u2 = SA.LastUParameter(), + v1 = SA.FirstVParameter(), + v2 = SA.LastVParameter(); + Handle(Geom_Surface) surf = BRep_Tool::Surface(aFace); + myCurProjFace.Init(surf, u1,u2, v1,v2); + myCurFace = aFace; + process(); + } + break; + case TopAbs_EDGE: + { + TopoDS_Edge anEdge = TopoDS::Edge(theShape); + Standard_Real u1, u2; + Handle(Geom_Curve) curve = BRep_Tool::Curve(anEdge, u1, u2); + myCurProjEdge.Init(curve, u1, u2); + process(); + } + break; + case TopAbs_VERTEX: + { + TopoDS_Vertex aV = TopoDS::Vertex(theShape); + myCurPnt = BRep_Tool::Pnt(aV); + process(); + } + break; + default: + break; + } +} + +void ElementsOnShape::process() +{ + if (myShape.IsNull() || myMesh == 0) + return; + + if (myType == SMDSAbs_Node) + { + SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator(); + while (anIter->more()) + process(anIter->next()); } else - return false; + { + if (myType == SMDSAbs_Edge || myType == SMDSAbs_All) + { + SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator(); + while (anIter->more()) + process(anIter->next()); + } - return true; + if (myType == SMDSAbs_Face || myType == SMDSAbs_All) + { + SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); + while (anIter->more()) { + process(anIter->next()); + } + } + + if (myType == SMDSAbs_Volume || myType == SMDSAbs_All) + { + SMDS_VolumeIteratorPtr anIter = myMesh->volumesIterator(); + while (anIter->more()) + process(anIter->next()); + } + } +} + +void ElementsOnShape::process (const SMDS_MeshElement* theElemPtr) +{ + if (myShape.IsNull()) + return; + + SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator(); + bool isSatisfy = myAllNodesFlag; + + gp_XYZ centerXYZ (0, 0, 0); + + while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) + { + SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); + gp_Pnt aPnt (aNode->X(), aNode->Y(), aNode->Z()); + centerXYZ += aPnt.XYZ(); + + switch (myCurShapeType) + { + case TopAbs_SOLID: + { + myCurSC.Perform(aPnt, myToler); + isSatisfy = (myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON); + } + break; + case TopAbs_FACE: + { + myCurProjFace.Perform(aPnt); + isSatisfy = (myCurProjFace.IsDone() && myCurProjFace.LowerDistance() <= myToler); + if (isSatisfy) + { + // check relatively the face + Quantity_Parameter u, v; + myCurProjFace.LowerDistanceParameters(u, v); + gp_Pnt2d aProjPnt (u, v); + BRepClass_FaceClassifier aClsf (myCurFace, aProjPnt, myToler); + isSatisfy = (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON); + } + } + break; + case TopAbs_EDGE: + { + myCurProjEdge.Perform(aPnt); + isSatisfy = (myCurProjEdge.NbPoints() > 0 && myCurProjEdge.LowerDistance() <= myToler); + } + break; + case TopAbs_VERTEX: + { + isSatisfy = (aPnt.Distance(myCurPnt) <= myToler); + } + break; + default: + { + isSatisfy = false; + } + } + } + + if (isSatisfy && myCurShapeType == TopAbs_SOLID) { // Check the center point for volumes MantisBug 0020168 + centerXYZ /= theElemPtr->NbNodes(); + gp_Pnt aCenterPnt (centerXYZ); + myCurSC.Perform(aCenterPnt, myToler); + if ( !(myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON)) + isSatisfy = false; + } + + if (isSatisfy) + myIds.Add(theElemPtr->GetID()); +} + +TSequenceOfXYZ::TSequenceOfXYZ() +{} + +TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n) +{} + +TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t) +{} + +TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray) +{} + +template +TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd) +{} + +TSequenceOfXYZ::~TSequenceOfXYZ() +{} + +TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ) +{ + myArray = theSequenceOfXYZ.myArray; + return *this; +} + +gp_XYZ& TSequenceOfXYZ::operator()(size_type n) +{ + return myArray[n-1]; +} + +const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const +{ + return myArray[n-1]; +} + +void TSequenceOfXYZ::clear() +{ + myArray.clear(); +} + +void TSequenceOfXYZ::reserve(size_type n) +{ + myArray.reserve(n); +} + +void TSequenceOfXYZ::push_back(const gp_XYZ& v) +{ + myArray.push_back(v); +} + +TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const +{ + return myArray.size(); }