X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=fa70861643ae1c33a859f1e342589d6212af79e5;hb=5601792b855f6e3f3f6d6cb22a37fb69026b7345;hp=198ea8da423e4742a3197bc86f5c27fedac1b6e6;hpb=8fa039a796957b302d86d90b22afc0a998573f83;p=modules%2Fsmesh.git diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 198ea8da4..fa7086164 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -1,52 +1,57 @@ -// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2013 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 +// 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 -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 2.1 of the License. +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. // -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// You should have received a copy of the GNU Lesser General Public +// 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.salome-platform.org/ or email : webmaster.salome@opencascade.com +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // #include "SMESH_ControlsDef.hxx" -#include -#include +#include "SMDS_BallElement.hxx" +#include "SMDS_Iterator.hxx" +#include "SMDS_Mesh.hxx" +#include "SMDS_MeshElement.hxx" +#include "SMDS_MeshNode.hxx" +#include "SMDS_QuadraticEdge.hxx" +#include "SMDS_QuadraticFaceOfNodes.hxx" +#include "SMDS_VolumeTool.hxx" +#include "SMESHDS_GroupBase.hxx" +#include "SMESHDS_Mesh.hxx" +#include "SMESH_OctreeNode.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 @@ -55,22 +60,21 @@ #include #include -#include "SMDS_Mesh.hxx" -#include "SMDS_Iterator.hxx" -#include "SMDS_MeshElement.hxx" -#include "SMDS_MeshNode.hxx" -#include "SMDS_VolumeTool.hxx" -#include "SMDS_QuadraticFaceOfNodes.hxx" -#include "SMDS_QuadraticEdge.hxx" +#include -#include "SMESHDS_Mesh.hxx" -#include "SMESHDS_GroupBase.hxx" +#include +#include + +#include /* AUXILIARY METHODS */ -namespace{ +namespace { + + const double theEps = 1e-100; + const double theInf = 1e+100; inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode ) { @@ -208,10 +212,13 @@ using namespace SMESH::Controls; * FUNCTORS */ +//================================================================================ /* Class : NumericalFunctor Description : Base class for numerical functors */ +//================================================================================ + NumericalFunctor::NumericalFunctor(): myMesh(NULL) { @@ -231,7 +238,11 @@ bool NumericalFunctor::GetPoints(const int theId, if ( myMesh == 0 ) return false; - return GetPoints( myMesh->FindElement( theId ), theRes ); + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + if ( !anElem || anElem->GetType() != this->GetType() ) + return false; + + return GetPoints( anElem, theRes ); } bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, @@ -239,7 +250,7 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, { theRes.clear(); - if ( anElem == 0) + if ( anElem == 0 ) return false; theRes.reserve( anElem->NbNodes() ); @@ -284,24 +295,25 @@ long NumericalFunctor::GetPrecision() const void NumericalFunctor::SetPrecision( const long thePrecision ) { myPrecision = thePrecision; + myPrecisionValue = pow( 10., (double)( myPrecision ) ); } double NumericalFunctor::GetValue( long theId ) { + double aVal = 0; + myCurrElement = myMesh->FindElement( theId ); + TSequenceOfXYZ P; if ( GetPoints( theId, P )) - { - double aVal = GetValue( P ); - if ( myPrecision >= 0 ) - { - double prec = pow( 10., (double)( myPrecision ) ); - aVal = floor( aVal * prec + 0.5 ) / prec; - } - return aVal; - } + aVal = Round( GetValue( P )); - return 0.; + return aVal; +} + +double NumericalFunctor::Round( const double & aVal ) +{ + return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal; } //================================================================================ @@ -310,12 +322,17 @@ double NumericalFunctor::GetValue( long theId ) * \param nbIntervals - number of intervals * \param nbEvents - number of mesh elements having values within i-th interval * \param funValues - boundaries of intervals + * \param elements - elements to check vulue of; empty list means "of all" + * \param minmax - boundaries of diapason of values to divide into intervals */ //================================================================================ void NumericalFunctor::GetHistogram(int nbIntervals, std::vector& nbEvents, - std::vector& funValues) + std::vector& funValues, + const vector& elements, + const double* minmax, + const bool isLogarithmic) { if ( nbIntervals < 1 || !myMesh || @@ -326,13 +343,30 @@ void NumericalFunctor::GetHistogram(int nbIntervals, // get all values sorted std::multiset< double > values; - SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType()); - while ( elemIt->more() ) - values.insert( GetValue( elemIt->next()->GetID() )); + if ( elements.empty() ) + { + SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType()); + while ( elemIt->more() ) + values.insert( GetValue( elemIt->next()->GetID() )); + } + else + { + vector::const_iterator id = elements.begin(); + for ( ; id != elements.end(); ++id ) + values.insert( GetValue( *id )); + } + if ( minmax ) + { + funValues[0] = minmax[0]; + funValues[nbIntervals] = minmax[1]; + } + else + { + funValues[0] = *values.begin(); + funValues[nbIntervals] = *values.rbegin(); + } // case nbIntervals == 1 - funValues[0] = *values.begin(); - funValues[nbIntervals] = *values.rbegin(); if ( nbIntervals == 1 ) { nbEvents[0] = values.size(); @@ -350,21 +384,36 @@ void NumericalFunctor::GetHistogram(int nbIntervals, 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; + // find end value of i-th interval + double r = (i+1) / double(nbIntervals); + if (isLogarithmic && funValues.front() > 1e-07 && funValues.back() > 1e-07) { + double logmin = log10(funValues.front()); + double lval = logmin + r * (log10(funValues.back()) - logmin); + funValues[i+1] = pow(10.0, lval); + } + else { + funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r; + } + + // count values in the i-th interval if there are any if ( min != values.end() && *min <= funValues[i+1] ) { - max = values.upper_bound( funValues[i+1] ); // greater than funValues[i+1], or end() + // find the first value out of the interval + max = values.upper_bound( funValues[i+1] ); // max is greater than funValues[i+1], or end() nbEvents[i] = std::distance( min, max ); min = max; } } + // add values larger than minmax[1] + nbEvents.back() += std::distance( min, values.end() ); } //======================================================================= -//function : GetValue -//purpose : -//======================================================================= +/* + Class : Volume + Description : Functor calculating volume of a 3D element +*/ +//================================================================================ double Volume::GetValue( long theElementId ) { @@ -376,86 +425,72 @@ double Volume::GetValue( long theElementId ) return 0; } -//======================================================================= -//function : GetBadRate -//purpose : meaningless as it is not quality control functor -//======================================================================= - double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const { return Value; } -//======================================================================= -//function : GetType -//purpose : -//======================================================================= - SMDSAbs_ElementType Volume::GetType() const { return SMDSAbs_Volume; } - +//======================================================================= /* Class : MaxElementLength2D Description : Functor calculating maximum length of 2D element */ +//================================================================================ + +double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P ) +{ + if(P.size() == 0) + return 0.; + double aVal = 0; + int len = P.size(); + 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)); + } + 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)); + } + 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)); + } + else if( len == 8 || len == 9 ) { // 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)); + } + + if( myPrecision >= 0 ) + { + double prec = pow( 10., (double)myPrecision ); + aVal = floor( aVal * prec + 0.5 ) / prec; + } + return aVal; +} 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.; + return GetPoints( theElementId, P ) ? GetValue(P) : 0.0; } double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -468,10 +503,12 @@ SMDSAbs_ElementType MaxElementLength2D::GetType() const return SMDSAbs_Face; } +//======================================================================= /* Class : MaxElementLength3D Description : Functor calculating maximum length of 3D element */ +//================================================================================ double MaxElementLength3D::GetValue( long theElementId ) { @@ -543,6 +580,12 @@ double MaxElementLength3D::GetValue( long theElementId ) aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); break; } + else if( len == 12 ) { // hexagonal prism + for ( int i1 = 1; i1 < 12; ++i1 ) + for ( int i2 = i1+1; i1 <= 12; ++i1 ) + aVal = Max( aVal, getDistance(P( i1 ),P( i2 ))); + 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 )); @@ -580,7 +623,7 @@ double MaxElementLength3D::GetValue( long theElementId ) aVal = Max(aVal,Max(Max(L7,L8),L9)); break; } - else if( len == 20 ) { // quadratic hexas + else if( len == 20 || len == 27 ) { // 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 )); @@ -636,11 +679,12 @@ SMDSAbs_ElementType MaxElementLength3D::GetType() const return SMDSAbs_Volume; } - +//======================================================================= /* Class : MinimumAngle Description : Functor for calculation of minimum angle */ +//================================================================================ double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) { @@ -657,7 +701,7 @@ double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) aMin = Min(aMin,A0); } - return aMin * 180.0 / PI; + return aMin * 180.0 / M_PI; } double MinimumAngle::GetBadRate( double Value, int nbNodes ) const @@ -673,10 +717,33 @@ SMDSAbs_ElementType MinimumAngle::GetType() const } +//================================================================================ /* Class : AspectRatio Description : Functor for calculating aspect ratio */ +//================================================================================ + +double AspectRatio::GetValue( long theId ) +{ + double aVal = 0; + myCurrElement = myMesh->FindElement( theId ); + if ( myCurrElement && myCurrElement->GetVtkType() == VTK_QUAD ) + { + // issue 21723 + vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid(); + if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() )) + aVal = Round( vtkMeshQuality::QuadAspectRatio( avtkCell )); + } + else + { + TSequenceOfXYZ P; + if ( GetPoints( myCurrElement, P )) + aVal = Round( GetValue( P )); + } + return aVal; +} + double AspectRatio::GetValue( const TSequenceOfXYZ& P ) { // According to "Mesh quality control" by Nadir Bouhamau referring to @@ -707,8 +774,8 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) 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.; + if ( anArea <= theEps ) + return theInf; return alfa * maxLen * half_perimeter / anArea; } else if ( nbNodes == 6 ) { // quadratic triangles @@ -727,8 +794,8 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) 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.; + if ( anArea <= theEps ) + return theInf; return alfa * maxLen * half_perimeter / anArea; } else if( nbNodes == 4 ) { // quadrangle @@ -771,11 +838,11 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) double C2 = Min( anArea[ 0 ], Min( anArea[ 1 ], Min( anArea[ 2 ], anArea[ 3 ] ) ) ); - if ( C2 <= Precision::Confusion() ) - return 0.; + if ( C2 <= theEps ) + return theInf; return alpha * L * C1 / C2; } - else if( nbNodes == 8 ){ // nbNodes==8 - quadratic quadrangle + else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle // Compute lengths of the sides std::vector< double > aLen (4); aLen[0] = getDistance( P(1), P(3) ); @@ -815,8 +882,8 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) double C2 = Min( anArea[ 0 ], Min( anArea[ 1 ], Min( anArea[ 2 ], anArea[ 3 ] ) ) ); - if ( C2 <= Precision::Confusion() ) - return 0.; + if ( C2 <= theEps ) + return theInf; return alpha * L * C1 / C2; } return 0; @@ -825,9 +892,10 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const { // the aspect ratio is in the range [1.0,infinity] + // < 1.0 = very bad, zero area // 1.0 = good // infinity = bad - return Value / 1000.; + return ( Value < 0.9 ) ? 1000 : Value / 1000.; } SMDSAbs_ElementType AspectRatio::GetType() const @@ -836,10 +904,13 @@ SMDSAbs_ElementType AspectRatio::GetType() const } +//================================================================================ /* Class : AspectRatio3D Description : Functor for calculating aspect ratio */ +//================================================================================ + namespace{ inline double getHalfPerimeter(double theTria[3]){ @@ -902,6 +973,28 @@ namespace{ } +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; @@ -914,10 +1007,11 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) 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){ + switch(nbNodes) { case 4:{ double aLen[6] = { getDistance(P( 1 ),P( 2 )), // a @@ -1135,7 +1229,22 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) } 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; @@ -1168,10 +1277,13 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const } +//================================================================================ /* Class : Warping Description : Functor for calculating warping */ +//================================================================================ + double Warping::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 4 ) @@ -1184,7 +1296,11 @@ double Warping::GetValue( const TSequenceOfXYZ& P ) 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 val = Max( Max( A1, A2 ), Max( A3, A4 ) ); + + const double eps = 0.1; // val is in degrees + + return val < eps ? 0. : val; } double Warping::ComputeA( const gp_XYZ& thePnt1, @@ -1195,20 +1311,20 @@ double Warping::ComputeA( const gp_XYZ& thePnt1, 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.; + if ( L < theEps ) + return theInf; 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 PI / 2; + return M_PI / 2; N.Normalize(); double H = ( thePnt2 - theG ).Dot( N ); - return asin( fabs( H / L ) ) * 180. / PI; + return asin( fabs( H / L ) ) * 180. / M_PI; } double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -1225,10 +1341,13 @@ SMDSAbs_ElementType Warping::GetType() const } +//================================================================================ /* Class : Taper Description : Functor for calculating taper */ +//================================================================================ + double Taper::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 4 ) @@ -1241,15 +1360,19 @@ double Taper::GetValue( const TSequenceOfXYZ& P ) double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2.; double JA = 0.25 * ( J1 + J2 + J3 + J4 ); - if ( JA <= Precision::Confusion() ) - return 0.; + if ( JA <= theEps ) + return theInf; 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 val = Max( Max( T1, T2 ), Max( T3, T4 ) ); + + const double eps = 0.01; + + return val < eps ? 0. : val; } double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -1265,11 +1388,13 @@ 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.; @@ -1287,14 +1412,14 @@ double Skew::GetValue( const TSequenceOfXYZ& P ) return 0.; // Compute skew - static double PI2 = PI / 2.; + const 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. / PI; + return Max( A0, Max( A1, A2 ) ) * 180. / M_PI; } else { @@ -1307,11 +1432,11 @@ double Skew::GetValue( const TSequenceOfXYZ& P ) double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution() ? 0. : fabs( PI2 - v1.Angle( v2 ) ); - //BUG SWP12743 - if ( A < Precision::Angular() ) - return 0.; + double val = A * 180. / M_PI; + + const double eps = 0.1; // val is in degrees - return A * 180. / PI; + return val < eps ? 0. : val; } } @@ -1329,10 +1454,13 @@ SMDSAbs_ElementType Skew::GetType() const } +//================================================================================ /* Class : Area Description : Functor for calculating area */ +//================================================================================ + double Area::GetValue( const TSequenceOfXYZ& P ) { double val = 0.0; @@ -1362,11 +1490,13 @@ SMDSAbs_ElementType Area::GetType() const return SMDSAbs_Face; } - +//================================================================================ /* Class : Length - Description : Functor for calculating length off edge + Description : Functor for calculating length of edge */ +//================================================================================ + double Length::GetValue( const TSequenceOfXYZ& P ) { switch ( P.size() ) { @@ -1387,10 +1517,12 @@ SMDSAbs_ElementType Length::GetType() const return SMDSAbs_Edge; } +//================================================================================ /* Class : Length2D Description : Functor for calculating length of edge */ +//================================================================================ double Length2D::GetValue( long theElementId) { @@ -1698,10 +1830,13 @@ void Length2D::GetValues(TValues& theValues){ } } +//================================================================================ /* Class : MultiConnection Description : Functor for calculating number of faces conneted to the edge */ +//================================================================================ + double MultiConnection::GetValue( const TSequenceOfXYZ& P ) { return 0; @@ -1722,10 +1857,13 @@ 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; @@ -1869,14 +2007,47 @@ void MultiConnection2D::GetValues(MValues& theValues){ } +//================================================================================ +/* + 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( 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() { @@ -1902,21 +2073,245 @@ SMDSAbs_ElementType BadOrientedVolume::GetType() const return SMDSAbs_Volume; } - - /* - Class : FreeBorders - Description : Predicate for free borders + Class : BareBorderVolume */ -FreeBorders::FreeBorders() +bool BareBorderVolume::IsSatisfy(long theElementId ) { - myMesh = 0; + 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; } -void FreeBorders::SetMesh( const SMDS_Mesh* theMesh ) +//================================================================================ +/* + Class : BareBorderFace +*/ +//================================================================================ + +bool BareBorderFace::IsSatisfy(long theElementId ) { - myMesh = theMesh; + 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 ) @@ -1930,10 +2325,13 @@ SMDSAbs_ElementType FreeBorders::GetType() const } +//================================================================================ /* Class : FreeEdges Description : Predicate for free Edges */ +//================================================================================ + FreeEdges::FreeEdges() { myMesh = 0; @@ -1949,17 +2347,13 @@ bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId TColStd_MapOfInteger aMap; for ( int i = 0; i < 2; i++ ) { - SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(); + SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face); while( anElemIter->more() ) { - const SMDS_MeshElement* anElem = anElemIter->next(); - if ( anElem != 0 && anElem->GetType() == SMDSAbs_Face ) + if ( const SMDS_MeshElement* anElem = anElemIter->next()) { - int anId = anElem->GetID(); - - if ( i == 0 ) - aMap.Add( anId ); - else if ( aMap.Contains( anId ) && anId != theFaceId ) + const int anId = anElem->GetID(); + if ( anId != theFaceId && !aMap.Add( anId )) return false; } } @@ -1984,7 +2378,7 @@ bool FreeEdges::IsSatisfy( long theId ) else { anIter = aFace->nodesIterator(); } - if ( anIter == 0 ) + if ( !anIter ) return false; int i = 0, nbNodes = aFace->NbNodes(); @@ -2062,21 +2456,19 @@ void FreeEdges::GetBoreders(TBorders& theBorders) long anId = aNode->GetID(); Border aBorder(anElemId,aNodeId[1],anId); aNodeId[1] = anId; - //std::cout<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; - } + bool isOk = ( anElem->GetGeomType() == myGeomType ); return isOk; } @@ -2377,6 +2745,54 @@ SMDSAbs_GeometryType ElemGeomType::GetGeomType() const return myGeomType; } +//================================================================================ +/* + Class : ElemEntityType + Description : Predicate to check element entity type +*/ +//================================================================================ + +ElemEntityType::ElemEntityType(): + myMesh( 0 ), + myType( SMDSAbs_All ), + myEntityType( SMDSEntity_0D ) +{ +} + +void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; +} + +bool ElemEntityType::IsSatisfy( long theId ) +{ + if ( !myMesh ) return false; + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + return ( anElem && + myEntityType == anElem->GetEntityType() && + ( myType == SMDSAbs_Edge || myType == SMDSAbs_Face || myType == SMDSAbs_Volume )); +} + +void ElemEntityType::SetType( SMDSAbs_ElementType theType ) +{ + myType = theType; +} + +SMDSAbs_ElementType ElemEntityType::GetType() const +{ + return myType; +} + +void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType ) +{ + myEntityType = theEntityType; +} + +SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const +{ + return myEntityType; +} + //================================================================================ /*! * \brief Class CoplanarFaces @@ -2384,59 +2800,72 @@ SMDSAbs_GeometryType ElemGeomType::GetGeomType() const //================================================================================ CoplanarFaces::CoplanarFaces() - : myMesh(0), myFaceID(0), myToler(0) + : myFaceID(0), myToler(0) { } -bool CoplanarFaces::IsSatisfy( long theElementId ) +void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh ) { - if ( myCoplanarIDs.empty() ) + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified() ) { // Build a set of coplanar face ids - if ( !myMesh || !myFaceID || !myToler ) - return false; + myCoplanarIDs.clear(); + + if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler ) + return; - const SMDS_MeshElement* face = myMesh->FindElement( myFaceID ); + const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID ); if ( !face || face->GetType() != SMDSAbs_Face ) - return false; + return; bool normOK; gp_Vec myNorm = getNormale( static_cast(face), &normOK ); if (!normOK) - return false; + return; + + const double radianTol = myToler * M_PI / 180.; + std::set< SMESH_TLink > checkedLinks; - const double radianTol = myToler * PI180; - typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TFaceIt; - std::set checkedFaces, checkedNodes; - std::list faceQueue( 1, face ); + std::list< pair< const SMDS_MeshElement*, gp_Vec > > faceQueue; + faceQueue.push_back( make_pair( face, myNorm )); while ( !faceQueue.empty() ) { - face = faceQueue.front(); - if ( checkedFaces.insert( face ).second ) + face = faceQueue.front().first; + myNorm = faceQueue.front().second; + faceQueue.pop_front(); + + for ( int i = 0, nbN = face->NbCornerNodes(); i < nbN; ++i ) { - gp_Vec norm = getNormale( static_cast(face), &normOK ); - if (!normOK || myNorm.Angle( norm ) <= radianTol) + const SMDS_MeshNode* n1 = face->GetNode( i ); + const SMDS_MeshNode* n2 = face->GetNode(( i+1 )%nbN); + if ( !checkedLinks.insert( SMESH_TLink( n1, n2 )).second ) + continue; + SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() ) { - myCoplanarIDs.insert( face->GetID() ); - std::set neighborFaces; - for ( int i = 0; i < face->NbCornerNodes(); ++i ) + const SMDS_MeshElement* f = fIt->next(); + if ( f->GetNodeIndex( n2 ) > -1 ) { - const SMDS_MeshNode* n = face->GetNode( i ); - if ( checkedNodes.insert( n ).second ) - neighborFaces.insert( TFaceIt( n->GetInverseElementIterator(SMDSAbs_Face)), - TFaceIt()); + gp_Vec norm = getNormale( static_cast(f), &normOK ); + if (!normOK || myNorm.Angle( norm ) <= radianTol) + { + myCoplanarIDs.insert( f->GetID() ); + faceQueue.push_back( make_pair( f, norm )); + } } - faceQueue.insert( faceQueue.end(), neighborFaces.begin(), neighborFaces.end() ); } } - faceQueue.pop_front(); } } +} +bool CoplanarFaces::IsSatisfy( long theElementId ) +{ return myCoplanarIDs.count( theElementId ); } /* - *Class : RangeOfIds + *Class : RangeOfIds *Description : Predicate for Range of Ids. * Range may be specified with two ways. * 1. Using AddToRange method @@ -2588,8 +3017,8 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr ) while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' ); while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' ); - if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() || - !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() ) + if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) || + (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) ) return false; myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() ); @@ -2635,7 +3064,7 @@ bool RangeOfIds::IsSatisfy( long theId ) else { const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); - if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All ) + if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All )) return false; } @@ -2826,8 +3255,8 @@ bool LogicalOR::IsSatisfy( long theId ) return myPredicate1 && myPredicate2 && - myPredicate1->IsSatisfy( theId ) || - myPredicate2->IsSatisfy( theId ); + (myPredicate1->IsSatisfy( theId ) || + myPredicate2->IsSatisfy( theId )); } @@ -2846,26 +3275,9 @@ void Filter::SetPredicate( PredicatePtr thePredicate ) myPredicate = thePredicate; } -template -inline void FillSequence(const TIterator& theIterator, - 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 ); - } - } -} - -void -Filter:: -GetElementsId( const SMDS_Mesh* theMesh, - PredicatePtr thePredicate, - TIdSequence& theSequence ) +void Filter::GetElementsId( const SMDS_Mesh* theMesh, + PredicatePtr thePredicate, + TIdSequence& theSequence ) { theSequence.clear(); @@ -2874,31 +3286,19 @@ GetElementsId( const SMDS_Mesh* theMesh, thePredicate->SetMesh( theMesh ); - SMDSAbs_ElementType aType = thePredicate->GetType(); - switch(aType){ - case SMDSAbs_Node: - FillSequence(theMesh->nodesIterator(),thePredicate,theSequence); - break; - case SMDSAbs_Edge: - FillSequence(theMesh->edgesIterator(),thePredicate,theSequence); - break; - case SMDSAbs_Face: - FillSequence(theMesh->facesIterator(),thePredicate,theSequence); - break; - case SMDSAbs_Volume: - FillSequence(theMesh->volumesIterator(),thePredicate,theSequence); - break; - case SMDSAbs_All: - FillSequence(theMesh->edgesIterator(),thePredicate,theSequence); - FillSequence(theMesh->facesIterator(),thePredicate,theSequence); - FillSequence(theMesh->volumesIterator(),thePredicate,theSequence); - break; + SMDS_ElemIteratorPtr elemIt = theMesh->elementsIterator( thePredicate->GetType() ); + if ( elemIt ) { + while ( elemIt->more() ) { + const SMDS_MeshElement* anElem = elemIt->next(); + long anId = anElem->GetID(); + if ( thePredicate->IsSatisfy( anId ) ) + theSequence.push_back( anId ); + } } } -void -Filter::GetElementsId( const SMDS_Mesh* theMesh, - Filter::TIdSequence& theSequence ) +void Filter::GetElementsId( const SMDS_Mesh* theMesh, + Filter::TIdSequence& theSequence ) { GetElementsId(theMesh,myPredicate,theSequence); } @@ -3253,7 +3653,6 @@ void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink, ElementsOnSurface::ElementsOnSurface() { - myMesh = 0; myIds.Clear(); myType = SMDSAbs_All; mySurf.Nullify(); @@ -3263,15 +3662,13 @@ ElementsOnSurface::ElementsOnSurface() ElementsOnSurface::~ElementsOnSurface() { - myMesh = 0; } void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh ) { - if ( myMesh == theMesh ) - return; - myMesh = theMesh; - process(); + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified()) + process(); } bool ElementsOnSurface::IsSatisfy( long theElementId ) @@ -3326,32 +3723,14 @@ void ElementsOnSurface::process() if ( mySurf.IsNull() ) return; - if ( myMesh == 0 ) + if ( !myMeshModifTracer.GetMesh() ) return; - if ( myType == SMDSAbs_Face || myType == SMDSAbs_All ) - { - myIds.ReSize( myMesh->NbFaces() ); - SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); - for(; anIter->more(); ) - process( anIter->next() ); - } + myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType )); - if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All ) - { - myIds.ReSize( myIds.Extent() + myMesh->NbEdges() ); - SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator(); - for(; anIter->more(); ) - process( anIter->next() ); - } - - if ( myType == SMDSAbs_Node ) - { - myIds.ReSize( myMesh->NbNodes() ); - SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator(); - for(; anIter->more(); ) - process( anIter->next() ); - } + SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType ); + for(; anIter->more(); ) + process( anIter->next() ); } void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr ) @@ -3409,29 +3788,16 @@ bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) */ ElementsOnShape::ElementsOnShape() - : myMesh(0), + : //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); + clearClassifiers(); } SMDSAbs_ElementType ElementsOnShape::GetType() const @@ -3454,213 +3820,163 @@ double ElementsOnShape::GetTolerance() const void ElementsOnShape::SetAllNodes (bool theAllNodes) { - if (myAllNodesFlag != theAllNodes) { - myAllNodesFlag = theAllNodes; - SetShape(myShape, myType); - } + myAllNodesFlag = theAllNodes; +} + +void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh) +{ + myMesh = theMesh; } void ElementsOnShape::SetShape (const TopoDS_Shape& theShape, const SMDSAbs_ElementType theType) { - myType = theType; + myType = theType; myShape = theShape; - myIds.Clear(); - - if (myMesh == 0) return; - - switch (myType) + if ( myShape.IsNull() ) return; + + TopTools_IndexedMapOfShape shapesMap; + TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }; + TopExp_Explorer sub; + for ( int i = 0; i < 4; ++i ) { - 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; + if ( shapesMap.IsEmpty() ) + for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() ) + shapesMap.Add( sub.Current() ); + if ( i > 0 ) + for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() ) + shapesMap.Add( sub.Current() ); } - myShapesMap.Clear(); - addShape(myShape); + clearClassifiers(); + myClassifiers.resize( shapesMap.Extent() ); + for ( int i = 0; i < shapesMap.Extent(); ++i ) + myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler ); } -void ElementsOnShape::addShape (const TopoDS_Shape& theShape) +void ElementsOnShape::clearClassifiers() { - if (theShape.IsNull() || myMesh == 0) - return; + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + delete myClassifiers[ i ]; + myClassifiers.clear(); +} - if (!myShapesMap.Add(theShape)) return; +bool ElementsOnShape::IsSatisfy (long elemId) +{ + const SMDS_MeshElement* elem = + ( myType == SMDSAbs_Node ? myMesh->FindNode( elemId ) : myMesh->FindElement( elemId )); + if ( !elem || myClassifiers.empty() ) + return false; - myCurShapeType = theShape.ShapeType(); - switch (myCurShapeType) + for ( size_t i = 0; i < myClassifiers.size(); ++i ) { - 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: + SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator(); + bool isSatisfy = myAllNodesFlag; + + gp_XYZ centerXYZ (0, 0, 0); + + while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) { - 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(); + SMESH_TNodeXYZ aPnt ( aNodeItr->next() ); + centerXYZ += aPnt; + isSatisfy = ! myClassifiers[i]->IsOut( aPnt ); } - break; - case TopAbs_VERTEX: + + // Check the center point for volumes MantisBug 0020168 + if (isSatisfy && + myAllNodesFlag && + myClassifiers[i]->ShapeType() == TopAbs_SOLID) { - TopoDS_Vertex aV = TopoDS::Vertex(theShape); - myCurPnt = BRep_Tool::Pnt(aV); - process(); + centerXYZ /= elem->NbNodes(); + isSatisfy = ! myClassifiers[i]->IsOut( centerXYZ ); } - break; - default: - break; + if ( isSatisfy ) + return true; } + + return false; } -void ElementsOnShape::process() +TopAbs_ShapeEnum ElementsOnShape::TClassifier::ShapeType() const { - if (myShape.IsNull() || myMesh == 0) - return; + return myShape.ShapeType(); +} + +bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p) +{ + return (this->*myIsOutFun)( p ); +} - if (myType == SMDSAbs_Node) +void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double theTol) +{ + myShape = theShape; + myTol = theTol; + switch ( myShape.ShapeType() ) { - SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator(); - while (anIter->more()) - process(anIter->next()); + case TopAbs_SOLID: { + mySolidClfr.Load(theShape); + myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid; + break; } - else - { - if (myType == SMDSAbs_Edge || myType == SMDSAbs_All) - { - SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator(); - while (anIter->more()) - process(anIter->next()); - } - - 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()); - } + case TopAbs_FACE: { + Standard_Real u1,u2,v1,v2; + Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape )); + surf->Bounds( u1,u2,v1,v2 ); + myProjFace.Init(surf, u1,u2, v1,v2, myTol ); + myIsOutFun = & ElementsOnShape::TClassifier::isOutOfFace; + break; + } + case TopAbs_EDGE: { + Standard_Real u1, u2; + Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge(theShape), u1, u2); + myProjEdge.Init(curve, u1, u2); + myIsOutFun = & ElementsOnShape::TClassifier::isOutOfEdge; + break; + } + case TopAbs_VERTEX:{ + myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) ); + myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex; + break; + } + default: + throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier"); } } -void ElementsOnShape::process (const SMDS_MeshElement* theElemPtr) +bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p) { - if (myShape.IsNull()) - return; - - SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator(); - bool isSatisfy = myAllNodesFlag; - - gp_XYZ centerXYZ (0, 0, 0); + mySolidClfr.Perform( p, myTol ); + return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON ); +} - while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) +bool ElementsOnShape::TClassifier::isOutOfFace (const gp_Pnt& p) +{ + myProjFace.Perform( p ); + if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol ) { - 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; - } - } + // check relatively to the face + Quantity_Parameter u, v; + myProjFace.LowerDistanceParameters(u, v); + gp_Pnt2d aProjPnt (u, v); + BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol ); + if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON ) + return false; } + return true; +} - 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; - } +bool ElementsOnShape::TClassifier::isOutOfEdge (const gp_Pnt& p) +{ + myProjEdge.Perform( p ); + return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol ); +} - if (isSatisfy) - myIds.Add(theElemPtr->GetID()); +bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p) +{ + return ( myVertexXYZ.Distance( p ) > myTol ); } + TSequenceOfXYZ::TSequenceOfXYZ() {} @@ -3715,3 +4031,24 @@ TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const { return myArray.size(); } + +TMeshModifTracer::TMeshModifTracer(): + myMeshModifTime(0), myMesh(0) +{ +} +void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh ) +{ + if ( theMesh != myMesh ) + myMeshModifTime = 0; + myMesh = theMesh; +} +bool TMeshModifTracer::IsMeshModified() +{ + bool modified = false; + if ( myMesh ) + { + modified = ( myMeshModifTime != myMesh->GetMTime() ); + myMeshModifTime = myMesh->GetMTime(); + } + return modified; +}