Salome HOME
PR: synchro V7_main tag mergefrom_V6_main_06Mar13
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
index 04a587dba823b7d45d5c18dc90b83c1a568363a7..eb1e5bdb95e8f0d7879d1588ac4dc96b45f5977f 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2011  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012  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
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-//
 
 #include "SMESH_ControlsDef.hxx"
 
-#include <set>
-#include <limits>
+#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 <BRepAdaptor_Surface.hxx>
 #include <BRepClass_FaceClassifier.hxx>
 #include <BRep_Tool.hxx>
-
-#include <TopAbs.hxx>
-#include <TopoDS.hxx>
-#include <TopoDS_Edge.hxx>
-#include <TopoDS_Face.hxx>
-#include <TopoDS_Shape.hxx>
-#include <TopoDS_Vertex.hxx>
-#include <TopoDS_Iterator.hxx>
-
 #include <Geom_CylindricalSurface.hxx>
 #include <Geom_Plane.hxx>
 #include <Geom_Surface.hxx>
-
 #include <Precision.hxx>
 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
 #include <TColStd_MapOfInteger.hxx>
 #include <TColStd_SequenceOfAsciiString.hxx>
 #include <TColgp_Array1OfXYZ.hxx>
-
+#include <TopAbs.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Iterator.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS_Vertex.hxx>
 #include <gp_Ax3.hxx>
 #include <gp_Cylinder.hxx>
 #include <gp_Dir.hxx>
 #include <gp_Vec.hxx>
 #include <gp_XYZ.hxx>
 
-#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 <vtkMeshQuality.h>
 
-#include "SMESHDS_Mesh.hxx"
-#include "SMESHDS_GroupBase.hxx"
+#include <set>
+#include <limits>
 
-#include <vtkMeshQuality.h>
 
 /*
                             AUXILIARY METHODS
 */
 
-namespace{
+namespace {
+
+  const double theEps = 1e-100;
+  const double theInf = 1e+100;
 
   inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
   {
@@ -210,10 +210,13 @@ using namespace SMESH::Controls;
  *                               FUNCTORS
  */
 
+//================================================================================
 /*
   Class       : NumericalFunctor
   Description : Base class for numerical functors
 */
+//================================================================================
+
 NumericalFunctor::NumericalFunctor():
   myMesh(NULL)
 {
@@ -233,7 +236,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,
@@ -241,7 +248,7 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
 {
   theRes.clear();
 
-  if ( anElem == 0)
+  if ( anElem == 0 )
     return false;
 
   theRes.reserve( anElem->NbNodes() );
@@ -322,7 +329,8 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
                                     std::vector<int>&    nbEvents,
                                     std::vector<double>& funValues,
                                     const vector<int>&   elements,
-                                    const double*        minmax)
+                                    const double*        minmax,
+                                    const bool           isLogarithmic)
 {
   if ( nbIntervals < 1 ||
        !myMesh ||
@@ -375,8 +383,15 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
   for ( int i = 0; i < nbIntervals; ++i )
   {
     // find end value of i-th interval
-    double r = (i+1) / double( nbIntervals );
-    funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
+    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] )
@@ -392,9 +407,11 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
 }
 
 //=======================================================================
-//function : GetValue
-//purpose  : 
-//=======================================================================
+/*
+  Class       : Volume
+  Description : Functor calculating volume of a 3D element
+*/
+//================================================================================
 
 double Volume::GetValue( long theElementId )
 {
@@ -406,86 +423,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
@@ -498,10 +501,12 @@ SMDSAbs_ElementType MaxElementLength2D::GetType() const
   return SMDSAbs_Face;
 }
 
+//=======================================================================
 /*
   Class       : MaxElementLength3D
   Description : Functor calculating maximum length of 3D element
 */
+//================================================================================
 
 double MaxElementLength3D::GetValue( long theElementId )
 {
@@ -573,6 +578,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 ));
@@ -610,7 +621,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 ));
@@ -666,11 +677,12 @@ SMDSAbs_ElementType MaxElementLength3D::GetType() const
   return SMDSAbs_Volume;
 }
 
-
+//=======================================================================
 /*
   Class       : MinimumAngle
   Description : Functor for calculation of minimum angle
 */
+//================================================================================
 
 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
 {
@@ -687,7 +699,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
@@ -703,10 +715,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
@@ -737,8 +772,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
@@ -757,8 +792,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
@@ -801,11 +836,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) );
@@ -845,8 +880,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;
@@ -855,9 +890,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
@@ -866,10 +902,13 @@ SMDSAbs_ElementType AspectRatio::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : AspectRatio3D
   Description : Functor for calculating aspect ratio
 */
+//================================================================================
+
 namespace{
 
   inline double getHalfPerimeter(double theTria[3]){
@@ -966,10 +1005,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
@@ -1187,7 +1227,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;
@@ -1220,10 +1275,13 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Warping
   Description : Functor for calculating warping
 */
+//================================================================================
+
 double Warping::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
@@ -1247,20 +1305,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
@@ -1277,10 +1335,13 @@ SMDSAbs_ElementType Warping::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Taper
   Description : Functor for calculating taper
 */
+//================================================================================
+
 double Taper::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
@@ -1293,8 +1354,8 @@ 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 );
@@ -1317,11 +1378,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.;
@@ -1339,14 +1402,14 @@ double Skew::GetValue( const TSequenceOfXYZ& P )
     return 0.;
 
   // Compute skew
-  static double PI2 = PI / 2.;
+  static double PI2 = M_PI / 2.;
   if ( P.size() == 3 )
   {
     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
     double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
     double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
 
-    return Max( A0, Max( A1, A2 ) ) * 180. / PI;
+    return Max( A0, Max( A1, A2 ) ) * 180. / M_PI;
   }
   else
   {
@@ -1360,10 +1423,10 @@ double Skew::GetValue( const TSequenceOfXYZ& P )
       ? 0. : fabs( PI2 - v1.Angle( v2 ) );
 
     //BUG SWP12743
-    if ( A < Precision::Angular() )
-      return 0.;
+    if ( A < theEps )
+      return theInf;
 
-    return A * 180. / PI;
+    return A * 180. / M_PI;
   }
 }
 
@@ -1381,10 +1444,13 @@ SMDSAbs_ElementType Skew::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Area
   Description : Functor for calculating area
 */
+//================================================================================
+
 double Area::GetValue( const TSequenceOfXYZ& P )
 {
   double val = 0.0;
@@ -1414,11 +1480,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() ) {
@@ -1439,10 +1507,12 @@ SMDSAbs_ElementType Length::GetType() const
   return SMDSAbs_Edge;
 }
 
+//================================================================================
 /*
   Class       : Length2D
   Description : Functor for calculating length of edge
 */
+//================================================================================
 
 double Length2D::GetValue( long theElementId)
 {
@@ -1750,10 +1820,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;
@@ -1774,10 +1847,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;
@@ -1921,14 +1997,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<const SMDS_BallElement*>( myMesh->FindElement( theId )))
+  {
+    diameter = ball->GetDiameter();
+  }
+  return diameter;
+}
+
+double BallDiameter::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+  // meaningless as it is not a quality control functor
+  return Value;
+}
+
+SMDSAbs_ElementType BallDiameter::GetType() const
+{
+  return SMDSAbs_Ball;
+}
+
+
 /*
                             PREDICATES
 */
 
+//================================================================================
 /*
   Class       : BadOrientedVolume
   Description : Predicate bad oriented volumes
 */
+//================================================================================
 
 BadOrientedVolume::BadOrientedVolume()
 {
@@ -1975,9 +2084,11 @@ bool BareBorderVolume::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : BareBorderFace
 */
+//================================================================================
 
 bool BareBorderFace::IsSatisfy(long theElementId )
 {
@@ -2001,10 +2112,11 @@ bool BareBorderFace::IsSatisfy(long theElementId )
         }
         if ( !isShared )
         {
-          myLinkNodes.resize( 2 + face->IsQuadratic());
+          const int iQuad = face->IsQuadratic();
+          myLinkNodes.resize( 2 + iQuad);
           myLinkNodes[0] = n1;
           myLinkNodes[1] = n2;
-          if ( face->IsQuadratic() )
+          if ( iQuad )
             myLinkNodes[2] = face->GetNode( i+nbN );
           ok = !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
         }
@@ -2014,9 +2126,11 @@ bool BareBorderFace::IsSatisfy(long theElementId )
   return ok;
 }
 
+//================================================================================
 /*
   Class       : OverConstrainedVolume
 */
+//================================================================================
 
 bool OverConstrainedVolume::IsSatisfy(long theElementId )
 {
@@ -2034,9 +2148,11 @@ bool OverConstrainedVolume::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : OverConstrainedFace
 */
+//================================================================================
 
 bool OverConstrainedFace::IsSatisfy(long theElementId )
 {
@@ -2067,10 +2183,116 @@ bool OverConstrainedFace::IsSatisfy(long theElementId )
   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()
 {
@@ -2093,10 +2315,13 @@ SMDSAbs_ElementType FreeBorders::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeEdges
   Description : Predicate for free Edges
 */
+//================================================================================
+
 FreeEdges::FreeEdges()
 {
   myMesh = 0;
@@ -2112,17 +2337,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;
       }
     }
@@ -2147,7 +2368,7 @@ bool FreeEdges::IsSatisfy( long theId )
   else {
     anIter = aFace->nodesIterator();
   }
-  if ( anIter == 0 )
+  if ( !anIter )
     return false;
 
   int i = 0, nbNodes = aFace->NbNodes();
@@ -2225,21 +2446,19 @@ void FreeEdges::GetBoreders(TBorders& theBorders)
       long anId = aNode->GetID();
       Border aBorder(anElemId,aNodeId[1],anId);
       aNodeId[1] = anId;
-      //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
       UpdateBorders(aBorder,aRegistry,theBorders);
     }
     Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
-    //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
     UpdateBorders(aBorder,aRegistry,theBorders);
   }
-  //std::cout<<"theBorders.size() = "<<theBorders.size()<<endl;
 }
 
-
+//================================================================================
 /*
   Class       : FreeNodes
   Description : Predicate for free nodes
 */
+//================================================================================
 
 FreeNodes::FreeNodes()
 {
@@ -2266,10 +2485,12 @@ SMDSAbs_ElementType FreeNodes::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeFaces
   Description : Predicate for free faces
 */
+//================================================================================
 
 FreeFaces::FreeFaces()
 {
@@ -2322,10 +2543,12 @@ SMDSAbs_ElementType FreeFaces::GetType() const
   return SMDSAbs_Face;
 }
 
+//================================================================================
 /*
   Class       : LinearOrQuadratic
   Description : Predicate to verify whether a mesh element is linear
 */
+//================================================================================
 
 LinearOrQuadratic::LinearOrQuadratic()
 {
@@ -2356,10 +2579,12 @@ SMDSAbs_ElementType LinearOrQuadratic::GetType() const
   return myType;
 }
 
+//================================================================================
 /*
   Class       : GroupColor
   Description : Functor for check color of group to whic mesh element belongs to
 */
+//================================================================================
 
 GroupColor::GroupColor()
 {
@@ -2448,6 +2673,7 @@ void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
 // Purpose : Get range as a string.
 //           Example: "1,2,3,50-60,63,67,70-"
 //=======================================================================
+
 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
 {
   theResStr.Clear();
@@ -2456,10 +2682,12 @@ void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
 }
 
+//================================================================================
 /*
   Class       : ElemGeomType
   Description : Predicate to check element geometry type
 */
+//================================================================================
 
 ElemGeomType::ElemGeomType()
 {
@@ -2482,41 +2710,7 @@ bool ElemGeomType::IsSatisfy( long theId )
   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;
-  }
+  bool isOk = ( anElem->GetGeomType() == myGeomType );
   return isOk;
 }
 
@@ -2540,6 +2734,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
@@ -2547,59 +2789,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<const SMDS_MeshFace*>(face), &normOK );
     if (!normOK)
-      return false;
+      return;
 
-    const double radianTol = myToler * PI180;
-    typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TFaceIt;
-    std::set<const SMDS_MeshElement*> checkedFaces, checkedNodes;
-    std::list<const SMDS_MeshElement*> faceQueue( 1, face );
+    const double radianTol = myToler * M_PI / 180.;
+    std::set< SMESH_TLink > checkedLinks;
+
+    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<const SMDS_MeshFace*>(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<const SMDS_MeshElement*> 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<const SMDS_MeshFace*>(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
@@ -3009,26 +3264,9 @@ void Filter::SetPredicate( PredicatePtr thePredicate )
   myPredicate = thePredicate;
 }
 
-template<class TElement, class TIterator, class TPredicate>
-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();
 
@@ -3037,31 +3275,19 @@ GetElementsId( const SMDS_Mesh* theMesh,
 
   thePredicate->SetMesh( theMesh );
 
-  SMDSAbs_ElementType aType = thePredicate->GetType();
-  switch(aType){
-  case SMDSAbs_Node:
-    FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),thePredicate,theSequence);
-    break;
-  case SMDSAbs_Edge:
-    FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
-    break;
-  case SMDSAbs_Face:
-    FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
-    break;
-  case SMDSAbs_Volume:
-    FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),thePredicate,theSequence);
-    break;
-  case SMDSAbs_All:
-    FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
-    FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
-    FillSequence<const SMDS_MeshElement*>(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);
 }
@@ -3416,7 +3642,6 @@ void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
 
 ElementsOnSurface::ElementsOnSurface()
 {
-  myMesh = 0;
   myIds.Clear();
   myType = SMDSAbs_All;
   mySurf.Nullify();
@@ -3426,15 +3651,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 )
@@ -3489,32 +3712,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() );
-  }
-
-  if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All )
-  {
-    myIds.ReSize( myIds.Extent() + myMesh->NbEdges() );
-    SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
-    for(; anIter->more(); )
-      process( anIter->next() );
-  }
+  myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
 
-  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 )
@@ -3572,29 +3777,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
@@ -3617,213 +3809,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();
+}
 
-  if (myType == SMDSAbs_Node)
+bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p)
+{
+  return (this->*myIsOutFun)( p );
+}
+
+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()
 {}
 
@@ -3878,3 +4020,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;
+}