Salome HOME
Fix regression of smesh/2D_mesh_Polygons_00/A2
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
index f9540e0046d526846c6e05291de07438572c5509..f86b5638719dfb201d6e15ab0a6e4487a1b0f174 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  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
@@ -6,7 +6,7 @@
 // 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.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -23,6 +23,7 @@
 #include "SMESH_ControlsDef.hxx"
 
 #include "SMDS_BallElement.hxx"
+#include "SMDS_FacePosition.hxx"
 #include "SMDS_Iterator.hxx"
 #include "SMDS_Mesh.hxx"
 #include "SMDS_MeshElement.hxx"
 #include "SMDS_QuadraticFaceOfNodes.hxx"
 #include "SMDS_VolumeTool.hxx"
 #include "SMESHDS_GroupBase.hxx"
+#include "SMESHDS_GroupOnFilter.hxx"
 #include "SMESHDS_Mesh.hxx"
+#include "SMESH_MeshAlgos.hxx"
 #include "SMESH_OctreeNode.hxx"
 
+#include <Basics_Utils.hxx>
+
 #include <BRepAdaptor_Surface.hxx>
+#include <BRepBndLib.hxx>
+#include <BRepBuilderAPI_Copy.hxx>
+#include <BRepClass3d_SolidClassifier.hxx>
 #include <BRepClass_FaceClassifier.hxx>
 #include <BRep_Tool.hxx>
+#include <GeomLib_IsPlanarSurface.hxx>
 #include <Geom_CylindricalSurface.hxx>
 #include <Geom_Plane.hxx>
 #include <Geom_Surface.hxx>
+#include <NCollection_Map.hxx>
 #include <Precision.hxx>
+#include <ShapeAnalysis_Surface.hxx>
 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
 #include <TColStd_MapOfInteger.hxx>
 #include <TColStd_SequenceOfAsciiString.hxx>
 #include <TColgp_Array1OfXYZ.hxx>
 #include <TopAbs.hxx>
+#include <TopExp.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Edge.hxx>
 #include <TopoDS_Face.hxx>
 #include <set>
 #include <limits>
 
-
 /*
                             AUXILIARY METHODS
 */
 
-namespace{
+namespace {
+
+  const double theEps = 1e-100;
+  const double theInf = 1e+100;
 
   inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
   {
@@ -85,6 +99,15 @@ namespace{
       v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
   }
 
+  inline double getCos2( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
+  {
+    gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
+    double dot = v1 * v2, len1 = v1.SquareMagnitude(), len2 = v2.SquareMagnitude();
+
+    return ( dot < 0 || len1 < gp::Resolution() || len2 < gp::Resolution() ? -1 :
+             dot * dot / len1 / len2 );
+  }
+
   inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
   {
     gp_Vec aVec1( P2 - P1 );
@@ -124,13 +147,13 @@ namespace{
     //  +-----+------+  +-----+------+ 
     //  |            |  |            |
     //  |            |  |            |
-    // result sould be 2 in both cases
+    // result should be 2 in both cases
     //
     int aResult0 = 0, aResult1 = 0;
      // last node, it is a medium one in a quadratic edge
     const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
-    const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 );
-    const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 );
+    const SMDS_MeshNode*    aNode0 = anEdge->GetNode( 0 );
+    const SMDS_MeshNode*    aNode1 = anEdge->GetNode( 1 );
     if ( aNode1 == aLastNode ) aNode1 = 0;
 
     SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator();
@@ -152,29 +175,6 @@ namespace{
     }
     int aResult = std::max ( aResult0, aResult1 );
 
-//     TColStd_MapOfInteger aMap;
-
-//     SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
-//     if ( anIter != 0 ) {
-//       while( anIter->more() ) {
-//      const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
-//      if ( aNode == 0 )
-//        return 0;
-//      SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
-//      while( anElemIter->more() ) {
-//        const SMDS_MeshElement* anElem = anElemIter->next();
-//        if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
-//          int anId = anElem->GetID();
-
-//          if ( anIter->more() )              // i.e. first node
-//            aMap.Add( anId );
-//          else if ( aMap.Contains( anId ) )
-//            aResult++;
-//        }
-//      }
-//       }
-//     }
-
     return aResult;
   }
 
@@ -190,7 +190,7 @@ namespace{
       n += q2 ^ q3;
     }
     double len = n.Modulus();
-    bool zeroLen = ( len <= numeric_limits<double>::min());
+    bool zeroLen = ( len <= std::numeric_limits<double>::min());
     if ( !zeroLen )
       n /= len;
 
@@ -208,10 +208,13 @@ using namespace SMESH::Controls;
  *                               FUNCTORS
  */
 
+//================================================================================
 /*
   Class       : NumericalFunctor
   Description : Base class for numerical functors
 */
+//================================================================================
+
 NumericalFunctor::NumericalFunctor():
   myMesh(NULL)
 {
@@ -223,7 +226,7 @@ void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh )
   myMesh = theMesh;
 }
 
-bool NumericalFunctor::GetPoints(const int theId,
+bool NumericalFunctor::GetPoints(const int       theId,
                                  TSequenceOfXYZ& theRes ) const
 {
   theRes.clear();
@@ -231,7 +234,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,10 +246,11 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
 {
   theRes.clear();
 
-  if ( anElem == 0)
+  if ( anElem == 0 )
     return false;
 
   theRes.reserve( anElem->NbNodes() );
+  theRes.setElement( anElem );
 
   // Get nodes of the element
   SMDS_ElemIteratorPtr anIter;
@@ -259,7 +267,6 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
       break;
     default:
       anIter = anElem->nodesIterator();
-      //return false;
     }
   }
   else {
@@ -267,9 +274,10 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
   }
 
   if ( anIter ) {
+    SMESH_NodeXYZ p;
     while( anIter->more() ) {
-      if ( const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( anIter->next() ))
-        theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
+      if ( p.Set( anIter->next() ))
+        theRes.push_back( p );
     }
   }
 
@@ -294,7 +302,7 @@ double NumericalFunctor::GetValue( long theId )
   myCurrElement = myMesh->FindElement( theId );
 
   TSequenceOfXYZ P;
-  if ( GetPoints( theId, P ))
+  if ( GetPoints( theId, P )) // elem type is checked here
     aVal = Round( GetValue( P ));
 
   return aVal;
@@ -316,11 +324,12 @@ double NumericalFunctor::Round( const double & aVal )
  */
 //================================================================================
 
-void NumericalFunctor::GetHistogram(int                  nbIntervals,
-                                    std::vector<int>&    nbEvents,
-                                    std::vector<double>& funValues,
-                                    const vector<int>&   elements,
-                                    const double*        minmax)
+void NumericalFunctor::GetHistogram(int                     nbIntervals,
+                                    std::vector<int>&       nbEvents,
+                                    std::vector<double>&    funValues,
+                                    const std::vector<int>& elements,
+                                    const double*           minmax,
+                                    const bool              isLogarithmic)
 {
   if ( nbIntervals < 1 ||
        !myMesh ||
@@ -333,13 +342,13 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
   std::multiset< double > values;
   if ( elements.empty() )
   {
-    SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType());
+    SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator( GetType() );
     while ( elemIt->more() )
       values.insert( GetValue( elemIt->next()->GetID() ));
   }
   else
   {
-    vector<int>::const_iterator id = elements.begin();
+    std::vector<int>::const_iterator id = elements.begin();
     for ( ; id != elements.end(); ++id )
       values.insert( GetValue( *id ));
   }
@@ -373,8 +382,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] )
@@ -390,9 +406,11 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
 }
 
 //=======================================================================
-//function : GetValue
-//purpose  : 
-//=======================================================================
+/*
+  Class       : Volume
+  Description : Functor calculating volume of a 3D element
+*/
+//================================================================================
 
 double Volume::GetValue( long theElementId )
 {
@@ -404,86 +422,93 @@ 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));
+  }
+  // Diagonals are undefined for concave polygons
+  // else if ( P.getElementEntity() == SMDSEntity_Quad_Polygon && P.size() > 2 ) // quad polygon
+  // {
+  //   // sides
+  //   aVal = getDistance( P( 1 ), P( P.size() )) + getDistance( P( P.size() ), P( P.size()-1 ));
+  //   for ( size_t i = 1; i < P.size()-1; i += 2 )
+  //   {
+  //     double L = getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 ));
+  //     aVal = Max( aVal, L );
+  //   }
+  //   // diagonals
+  //   for ( int i = P.size()-5; i > 0; i -= 2 )
+  //     for ( int j = i + 4; j < P.size() + i - 2; i += 2 )
+  //     {
+  //       double D = getDistance( P( i ), P( j ));
+  //       aVal = Max( aVal, D );
+  //     }
+  // }
+  // { // polygons
+    
+  // }
+
+  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 || 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));
-        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
@@ -496,10 +521,12 @@ SMDSAbs_ElementType MaxElementLength2D::GetType() const
   return SMDSAbs_Face;
 }
 
+//=======================================================================
 /*
   Class       : MaxElementLength3D
   Description : Functor calculating maximum length of 3D element
 */
+//================================================================================
 
 double MaxElementLength3D::GetValue( long theElementId )
 {
@@ -507,148 +534,165 @@ double MaxElementLength3D::GetValue( long theElementId )
   if( GetPoints( theElementId, P ) ) {
     double aVal = 0;
     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
-    SMDSAbs_ElementType aType = aElem->GetType();
+    SMDSAbs_EntityType      aType = aElem->GetEntityType();
     int len = P.size();
-    switch( aType ) {
-    case SMDSAbs_Volume:
-      if( len == 4 ) { // tetras
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 1 ));
-        double L4 = getDistance(P( 1 ),P( 4 ));
-        double L5 = getDistance(P( 2 ),P( 4 ));
-        double L6 = getDistance(P( 3 ),P( 4 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        break;
-      }
-      else if( len == 5 ) { // pyramids
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 4 ));
-        double L4 = getDistance(P( 4 ),P( 1 ));
-        double L5 = getDistance(P( 1 ),P( 5 ));
-        double L6 = getDistance(P( 2 ),P( 5 ));
-        double L7 = getDistance(P( 3 ),P( 5 ));
-        double L8 = getDistance(P( 4 ),P( 5 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(L7,L8));
-        break;
-      }
-      else if( len == 6 ) { // pentas
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 1 ));
-        double L4 = getDistance(P( 4 ),P( 5 ));
-        double L5 = getDistance(P( 5 ),P( 6 ));
-        double L6 = getDistance(P( 6 ),P( 4 ));
-        double L7 = getDistance(P( 1 ),P( 4 ));
-        double L8 = getDistance(P( 2 ),P( 5 ));
-        double L9 = getDistance(P( 3 ),P( 6 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(Max(L7,L8),L9));
-        break;
-      }
-      else if( len == 8 ) { // hexas
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 4 ));
-        double L4 = getDistance(P( 4 ),P( 1 ));
-        double L5 = getDistance(P( 5 ),P( 6 ));
-        double L6 = getDistance(P( 6 ),P( 7 ));
-        double L7 = getDistance(P( 7 ),P( 8 ));
-        double L8 = getDistance(P( 8 ),P( 5 ));
-        double L9 = getDistance(P( 1 ),P( 5 ));
-        double L10= getDistance(P( 2 ),P( 6 ));
-        double L11= getDistance(P( 3 ),P( 7 ));
-        double L12= getDistance(P( 4 ),P( 8 ));
-        double D1 = getDistance(P( 1 ),P( 7 ));
-        double D2 = getDistance(P( 2 ),P( 8 ));
-        double D3 = getDistance(P( 3 ),P( 5 ));
-        double D4 = getDistance(P( 4 ),P( 6 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
-        aVal = Max(aVal,Max(L11,L12));
-        aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
-        break;
-      }
-      else if( len == 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 ));
-        double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
-        double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
-        double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
-        double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        break;
-      }
-      else if( len == 13 ) { // quadratic pyramids
-        double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
-        double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
-        double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
-        double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
-        double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
-        double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(L7,L8));
-        break;
-      }
-      else if( len == 15 ) { // quadratic pentas
-        double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
-        double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
-        double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
-        double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
-        double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
-        double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
-        double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(Max(L7,L8),L9));
-        break;
-      }
-      else if( len == 20 || 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 ));
-        double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
-        double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
-        double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
-        double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
-        double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
-        double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
-        double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
-        double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
-        double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
-        double D1 = getDistance(P( 1 ),P( 7 ));
-        double D2 = getDistance(P( 2 ),P( 8 ));
-        double D3 = getDistance(P( 3 ),P( 5 ));
-        double D4 = getDistance(P( 4 ),P( 6 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
-        aVal = Max(aVal,Max(L11,L12));
-        aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
-        break;
-      }
-      else if( len > 1 && aElem->IsPoly() ) { // polys
-        // get the maximum distance between all pairs of nodes
-        for( int i = 1; i <= len; i++ ) {
-          for( int j = 1; j <= len; j++ ) {
-            if( j > i ) { // optimization of the loop
-              double D = getDistance( P(i), P(j) );
-              aVal = Max( aVal, D );
-            }
+    switch ( aType ) {
+    case SMDSEntity_Tetra: { // tetras
+      double L1 = getDistance(P( 1 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 1 ));
+      double L4 = getDistance(P( 1 ),P( 4 ));
+      double L5 = getDistance(P( 2 ),P( 4 ));
+      double L6 = getDistance(P( 3 ),P( 4 ));
+      aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+      break;
+    }
+    case SMDSEntity_Pyramid: { // pyramids
+      double L1 = getDistance(P( 1 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 4 ));
+      double L4 = getDistance(P( 4 ),P( 1 ));
+      double L5 = getDistance(P( 1 ),P( 5 ));
+      double L6 = getDistance(P( 2 ),P( 5 ));
+      double L7 = getDistance(P( 3 ),P( 5 ));
+      double L8 = getDistance(P( 4 ),P( 5 ));
+      aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+      aVal = Max(aVal,Max(L7,L8));
+      break;
+    }
+    case SMDSEntity_Penta: { // pentas
+      double L1 = getDistance(P( 1 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 1 ));
+      double L4 = getDistance(P( 4 ),P( 5 ));
+      double L5 = getDistance(P( 5 ),P( 6 ));
+      double L6 = getDistance(P( 6 ),P( 4 ));
+      double L7 = getDistance(P( 1 ),P( 4 ));
+      double L8 = getDistance(P( 2 ),P( 5 ));
+      double L9 = getDistance(P( 3 ),P( 6 ));
+      aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+      aVal = Max(aVal,Max(Max(L7,L8),L9));
+      break;
+    }
+    case SMDSEntity_Hexa: { // hexas
+      double L1 = getDistance(P( 1 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 4 ));
+      double L4 = getDistance(P( 4 ),P( 1 ));
+      double L5 = getDistance(P( 5 ),P( 6 ));
+      double L6 = getDistance(P( 6 ),P( 7 ));
+      double L7 = getDistance(P( 7 ),P( 8 ));
+      double L8 = getDistance(P( 8 ),P( 5 ));
+      double L9 = getDistance(P( 1 ),P( 5 ));
+      double L10= getDistance(P( 2 ),P( 6 ));
+      double L11= getDistance(P( 3 ),P( 7 ));
+      double L12= getDistance(P( 4 ),P( 8 ));
+      double D1 = getDistance(P( 1 ),P( 7 ));
+      double D2 = getDistance(P( 2 ),P( 8 ));
+      double D3 = getDistance(P( 3 ),P( 5 ));
+      double D4 = getDistance(P( 4 ),P( 6 ));
+      aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+      aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+      aVal = Max(aVal,Max(L11,L12));
+      aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
+      break;
+    }
+    case SMDSEntity_Hexagonal_Prism: { // 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;
+    }
+    case SMDSEntity_Quad_Tetra: { // quadratic tetras
+      double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
+      double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+      double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
+      double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
+      aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+      break;
+    }
+    case SMDSEntity_Quad_Pyramid: { // quadratic pyramids
+      double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+      double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+      double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+      double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
+      double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
+      double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
+      aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+      aVal = Max(aVal,Max(L7,L8));
+      break;
+    }
+    case SMDSEntity_Quad_Penta:
+    case SMDSEntity_BiQuad_Penta: { // quadratic pentas
+      double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+      double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+      double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
+      double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
+      double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
+      double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
+      double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
+      aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+      aVal = Max(aVal,Max(Max(L7,L8),L9));
+      break;
+    }
+    case SMDSEntity_Quad_Hexa:
+    case SMDSEntity_TriQuad_Hexa: { // quadratic hexas
+      double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
+      double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
+      double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
+      double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
+      double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
+      double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
+      double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
+      double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
+      double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
+      double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
+      double D1 = getDistance(P( 1 ),P( 7 ));
+      double D2 = getDistance(P( 2 ),P( 8 ));
+      double D3 = getDistance(P( 3 ),P( 5 ));
+      double D4 = getDistance(P( 4 ),P( 6 ));
+      aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+      aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+      aVal = Max(aVal,Max(L11,L12));
+      aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
+      break;
+    }
+    case SMDSEntity_Quad_Polyhedra:
+    case SMDSEntity_Polyhedra: { // polys
+      // get the maximum distance between all pairs of nodes
+      for( int i = 1; i <= len; i++ ) {
+        for( int j = 1; j <= len; j++ ) {
+          if( j > i ) { // optimization of the loop
+            double D = getDistance( P(i), P(j) );
+            aVal = Max( aVal, D );
           }
         }
       }
+      break;
     }
+    case SMDSEntity_Node:
+    case SMDSEntity_0D:
+    case SMDSEntity_Edge:
+    case SMDSEntity_Quad_Edge:
+    case SMDSEntity_Triangle:
+    case SMDSEntity_Quad_Triangle:
+    case SMDSEntity_BiQuad_Triangle:
+    case SMDSEntity_Quadrangle:
+    case SMDSEntity_Quad_Quadrangle:
+    case SMDSEntity_BiQuad_Quadrangle:
+    case SMDSEntity_Polygon:
+    case SMDSEntity_Quad_Polygon:
+    case SMDSEntity_Ball:
+    case SMDSEntity_Last: return 0;
+    } // switch ( aType )
 
     if( myPrecision >= 0 )
     {
@@ -670,28 +714,34 @@ SMDSAbs_ElementType MaxElementLength3D::GetType() const
   return SMDSAbs_Volume;
 }
 
-
+//=======================================================================
 /*
   Class       : MinimumAngle
   Description : Functor for calculation of minimum angle
 */
+//================================================================================
 
 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
 {
-  double aMin;
-
-  if (P.size() <3)
+  if ( P.size() < 3 )
     return 0.;
 
-  aMin = getAngle(P( P.size() ), P( 1 ), P( 2 ));
-  aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 )));
+  double aMaxCos2;
 
-  for (int i=2; i<P.size();i++){
-      double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) );
-    aMin = Min(aMin,A0);
+  aMaxCos2 = getCos2( P( P.size() ), P( 1 ), P( 2 ));
+  aMaxCos2 = Max( aMaxCos2, getCos2( P( P.size()-1 ), P( P.size() ), P( 1 )));
+
+  for ( size_t i = 2; i < P.size(); i++ )
+  {
+    double A0 = getCos2( P( i-1 ), P( i ), P( i+1 ) );
+    aMaxCos2 = Max( aMaxCos2, A0 );
   }
+  if ( aMaxCos2 < 0 )
+    return 0; // all nodes coincide
 
-  return aMin * 180.0 / M_PI;
+  double cos = sqrt( aMaxCos2 );
+  if ( cos >=  1 ) return 0;
+  return acos( cos ) * 180.0 / M_PI;
 }
 
 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
@@ -707,10 +757,13 @@ SMDSAbs_ElementType MinimumAngle::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : AspectRatio
   Description : Functor for calculating aspect ratio
 */
+//================================================================================
+
 double AspectRatio::GetValue( long theId )
 {
   double aVal = 0;
@@ -747,58 +800,51 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
 
   if ( nbNodes == 3 ) {
     // Compute lengths of the sides
-    std::vector< double > aLen (nbNodes);
-    for ( int i = 0; i < nbNodes - 1; i++ )
-      aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
-    aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
+    double aLen1 = getDistance( P( 1 ), P( 2 ));
+    double aLen2 = getDistance( P( 2 ), P( 3 ));
+    double aLen3 = getDistance( P( 3 ), P( 1 ));
     // Q = alfa * h * p / S, where
     //
     // alfa = sqrt( 3 ) / 6
     // h - length of the longest edge
     // p - half perimeter
     // S - triangle surface
-    const double alfa = sqrt( 3. ) / 6.;
-    double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
-    double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
-    double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
-    if ( anArea <= Precision::Confusion() )
-      return 0.;
+    const double     alfa = sqrt( 3. ) / 6.;
+    double         maxLen = Max( aLen1, Max( aLen2, aLen3 ));
+    double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.;
+    double         anArea = getArea( P( 1 ), P( 2 ), P( 3 ));
+    if ( anArea <= theEps  )
+      return theInf;
     return alfa * maxLen * half_perimeter / anArea;
   }
   else if ( nbNodes == 6 ) { // quadratic triangles
     // Compute lengths of the sides
-    std::vector< double > aLen (3);
-    aLen[0] = getDistance( P(1), P(3) );
-    aLen[1] = getDistance( P(3), P(5) );
-    aLen[2] = getDistance( P(5), P(1) );
-    // Q = alfa * h * p / S, where
-    //
-    // alfa = sqrt( 3 ) / 6
-    // h - length of the longest edge
-    // p - half perimeter
-    // S - triangle surface
-    const double alfa = sqrt( 3. ) / 6.;
-    double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
-    double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
-    double anArea = getArea( P(1), P(3), P(5) );
-    if ( anArea <= Precision::Confusion() )
-      return 0.;
+    double aLen1 = getDistance( P( 1 ), P( 3 ));
+    double aLen2 = getDistance( P( 3 ), P( 5 ));
+    double aLen3 = getDistance( P( 5 ), P( 1 ));
+    // algo same as for the linear triangle
+    const double     alfa = sqrt( 3. ) / 6.;
+    double         maxLen = Max( aLen1, Max( aLen2, aLen3 ));
+    double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.;
+    double         anArea = getArea( P( 1 ), P( 3 ), P( 5 ));
+    if ( anArea <= theEps )
+      return theInf;
     return alfa * maxLen * half_perimeter / anArea;
   }
   else if( nbNodes == 4 ) { // quadrangle
     // Compute lengths of the sides
-    std::vector< double > aLen (4);
+    double aLen[4];
     aLen[0] = getDistance( P(1), P(2) );
     aLen[1] = getDistance( P(2), P(3) );
     aLen[2] = getDistance( P(3), P(4) );
     aLen[3] = getDistance( P(4), P(1) );
     // Compute lengths of the diagonals
-    std::vector< double > aDia (2);
+    double aDia[2];
     aDia[0] = getDistance( P(1), P(3) );
     aDia[1] = getDistance( P(2), P(4) );
     // Compute areas of all triangles which can be built
     // taking three nodes of the quadrangle
-    std::vector< double > anArea (4);
+    double anArea[4];
     anArea[0] = getArea( P(1), P(2), P(3) );
     anArea[1] = getArea( P(1), P(2), P(4) );
     anArea[2] = getArea( P(1), P(3), P(4) );
@@ -814,35 +860,35 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
     // Si - areas of the triangles
     const double alpha = sqrt( 1 / 32. );
     double L = Max( aLen[ 0 ],
-                 Max( aLen[ 1 ],
-                   Max( aLen[ 2 ],
-                     Max( aLen[ 3 ],
-                       Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
+                    Max( aLen[ 1 ],
+                         Max( aLen[ 2 ],
+                              Max( aLen[ 3 ],
+                                   Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
     double C1 = sqrt( ( aLen[0] * aLen[0] +
                         aLen[1] * aLen[1] +
                         aLen[2] * aLen[2] +
                         aLen[3] * aLen[3] ) / 4. );
     double C2 = Min( anArea[ 0 ],
-                  Min( anArea[ 1 ],
-                    Min( anArea[ 2 ], anArea[ 3 ] ) ) );
-    if ( C2 <= Precision::Confusion() )
-      return 0.;
+                     Min( anArea[ 1 ],
+                          Min( anArea[ 2 ], anArea[ 3 ] ) ) );
+    if ( C2 <= theEps )
+      return theInf;
     return alpha * L * C1 / C2;
   }
   else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle
     // Compute lengths of the sides
-    std::vector< double > aLen (4);
+    double aLen[4];
     aLen[0] = getDistance( P(1), P(3) );
     aLen[1] = getDistance( P(3), P(5) );
     aLen[2] = getDistance( P(5), P(7) );
     aLen[3] = getDistance( P(7), P(1) );
     // Compute lengths of the diagonals
-    std::vector< double > aDia (2);
+    double aDia[2];
     aDia[0] = getDistance( P(1), P(5) );
     aDia[1] = getDistance( P(3), P(7) );
     // Compute areas of all triangles which can be built
     // taking three nodes of the quadrangle
-    std::vector< double > anArea (4);
+    double anArea[4];
     anArea[0] = getArea( P(1), P(3), P(5) );
     anArea[1] = getArea( P(1), P(3), P(7) );
     anArea[2] = getArea( P(1), P(5), P(7) );
@@ -869,8 +915,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;
@@ -891,10 +937,13 @@ SMDSAbs_ElementType AspectRatio::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : AspectRatio3D
   Description : Functor for calculating aspect ratio
 */
+//================================================================================
+
 namespace{
 
   inline double getHalfPerimeter(double theTria[3]){
@@ -1261,10 +1310,13 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Warping
   Description : Functor for calculating warping
 */
+//================================================================================
+
 double Warping::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
@@ -1277,7 +1329,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,
@@ -1288,8 +1344,8 @@ 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;
@@ -1318,37 +1374,44 @@ SMDSAbs_ElementType Warping::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Taper
   Description : Functor for calculating taper
 */
+//================================================================================
+
 double Taper::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
     return 0.;
 
   // Compute taper
-  double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2.;
-  double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2.;
-  double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2.;
-  double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2.;
+  double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) );
+  double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) );
+  double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) );
+  double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) );
 
   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
 {
   // the taper is in the range [0.0,1.0]
-  // 0.0  = good (no taper)
+  // 0.0 = good (no taper)
   // 1.0 = bad  (les cotes opposes sont allignes)
   return Value;
 }
@@ -1358,11 +1421,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.;
@@ -1380,7 +1445,7 @@ double Skew::GetValue( const TSequenceOfXYZ& P )
     return 0.;
 
   // Compute skew
-  static double PI2 = M_PI / 2.;
+  const double PI2 = M_PI / 2.;
   if ( P.size() == 3 )
   {
     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
@@ -1400,11 +1465,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;
 
-    return A * 180. / M_PI;
+    const double eps = 0.1; // val is in degrees
+
+    return val < eps ? 0. : val;
   }
 }
 
@@ -1422,20 +1487,26 @@ SMDSAbs_ElementType Skew::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Area
   Description : Functor for calculating area
 */
+//================================================================================
+
 double Area::GetValue( const TSequenceOfXYZ& P )
 {
   double val = 0.0;
-  if ( P.size() > 2 ) {
+  if ( P.size() > 2 )
+  {
     gp_Vec aVec1( P(2) - P(1) );
     gp_Vec aVec2( P(3) - P(1) );
     gp_Vec SumVec = aVec1 ^ aVec2;
-    for (int i=4; i<=P.size(); i++) {
+
+    for (size_t i=4; i<=P.size(); i++)
+    {
       gp_Vec aVec1( P(i-1) - P(1) );
-      gp_Vec aVec2( P(i) - P(1) );
+      gp_Vec aVec2( P(i  ) - P(1) );
       gp_Vec tmp = aVec1 ^ aVec2;
       SumVec.Add(tmp);
     }
@@ -1455,11 +1526,13 @@ SMDSAbs_ElementType Area::GetType() const
   return SMDSAbs_Face;
 }
 
-
+//================================================================================
 /*
   Class       : Length
   Description : Functor for calculating length of edge
 */
+//================================================================================
+
 double Length::GetValue( const TSequenceOfXYZ& P )
 {
   switch ( P.size() ) {
@@ -1480,210 +1553,252 @@ SMDSAbs_ElementType Length::GetType() const
   return SMDSAbs_Edge;
 }
 
+//================================================================================
 /*
   Class       : Length2D
-  Description : Functor for calculating length of edge
+  Description : Functor for calculating minimal length of edge
 */
+//================================================================================
 
-double Length2D::GetValue( long theElementId)
+double Length2D::GetValue( const TSequenceOfXYZ& P )
 {
-  TSequenceOfXYZ P;
-
-  //cout<<"Length2D::GetValue"<<endl;
-  if (GetPoints(theElementId,P)){
-    //for(int jj=1; jj<=P.size(); jj++)
-    //  cout<<"jj="<<jj<<" P("<<P(jj).X()<<","<<P(jj).Y()<<","<<P(jj).Z()<<")"<<endl;
-
-    double  aVal;// = GetValue( P );
-    const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
-    SMDSAbs_ElementType aType = aElem->GetType();
-
-    int len = P.size();
-
-    switch (aType){
-    case SMDSAbs_All:
-    case SMDSAbs_Node:
-    case SMDSAbs_Edge:
-      if (len == 2){
-        aVal = getDistance( P( 1 ), P( 2 ) );
-        break;
-      }
-      else if (len == 3){ // quadratic edge
-        aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
-        break;
-      }
-    case SMDSAbs_Face:
-      if (len == 3){ // triangles
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 1 ));
-        aVal = Max(L1,Max(L2,L3));
-        break;
-      }
-      else if (len == 4){ // quadrangles
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 4 ));
-        double L4 = getDistance(P( 4 ),P( 1 ));
-        aVal = Max(Max(L1,L2),Max(L3,L4));
-        break;
-      }
-      if (len == 6){ // quadratic triangles
-        double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
-        double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
-        double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
-        aVal = Max(L1,Max(L2,L3));
-        //cout<<"L1="<<L1<<" L2="<<L2<<"L3="<<L3<<" aVal="<<aVal<<endl;
-        break;
-      }
-      else if (len == 8){ // quadratic quadrangles
-        double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
-        double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
-        double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
-        double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
-        aVal = Max(Max(L1,L2),Max(L3,L4));
-        break;
-      }
-    case SMDSAbs_Volume:
-      if (len == 4){ // tetraidrs
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 1 ));
-        double L4 = getDistance(P( 1 ),P( 4 ));
-        double L5 = getDistance(P( 2 ),P( 4 ));
-        double L6 = getDistance(P( 3 ),P( 4 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        break;
-      }
-      else if (len == 5){ // piramids
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 4 ));
-        double L4 = getDistance(P( 4 ),P( 1 ));
-        double L5 = getDistance(P( 1 ),P( 5 ));
-        double L6 = getDistance(P( 2 ),P( 5 ));
-        double L7 = getDistance(P( 3 ),P( 5 ));
-        double L8 = getDistance(P( 4 ),P( 5 ));
-
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(L7,L8));
-        break;
-      }
-      else if (len == 6){ // pentaidres
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 1 ));
-        double L4 = getDistance(P( 4 ),P( 5 ));
-        double L5 = getDistance(P( 5 ),P( 6 ));
-        double L6 = getDistance(P( 6 ),P( 4 ));
-        double L7 = getDistance(P( 1 ),P( 4 ));
-        double L8 = getDistance(P( 2 ),P( 5 ));
-        double L9 = getDistance(P( 3 ),P( 6 ));
-
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(Max(L7,L8),L9));
-        break;
-      }
-      else if (len == 8){ // hexaider
-        double L1 = getDistance(P( 1 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 4 ));
-        double L4 = getDistance(P( 4 ),P( 1 ));
-        double L5 = getDistance(P( 5 ),P( 6 ));
-        double L6 = getDistance(P( 6 ),P( 7 ));
-        double L7 = getDistance(P( 7 ),P( 8 ));
-        double L8 = getDistance(P( 8 ),P( 5 ));
-        double L9 = getDistance(P( 1 ),P( 5 ));
-        double L10= getDistance(P( 2 ),P( 6 ));
-        double L11= getDistance(P( 3 ),P( 7 ));
-        double L12= getDistance(P( 4 ),P( 8 ));
-
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
-        aVal = Max(aVal,Max(L11,L12));
-        break;
-
-      }
-
-      if (len == 10){ // quadratic tetraidrs
-        double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
-        double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
-        double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
-        double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        break;
-      }
-      else if (len == 13){ // quadratic piramids
-        double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
-        double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
-        double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
-        double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
-        double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
-        double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(L7,L8));
-        break;
-      }
-      else if (len == 15){ // quadratic pentaidres
-        double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
-        double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
-        double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
-        double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
-        double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
-        double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
-        double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(Max(L7,L8),L9));
-        break;
-      }
-      else if (len == 20){ // quadratic hexaider
-        double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
-        double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
-        double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
-        double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
-        double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
-        double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
-        double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
-        double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
-        double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
-        double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
-        double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
-        double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
-        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-        aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
-        aVal = Max(aVal,Max(L11,L12));
-        break;
-
-      }
+  double aVal = 0;
+  int len = P.size();
+  SMDSAbs_EntityType aType = P.getElementEntity();
 
-    default: aVal=-1;
+  switch (aType) {
+  case SMDSEntity_Edge:
+    if (len == 2)
+      aVal = getDistance( P( 1 ), P( 2 ) );
+    break;
+  case SMDSEntity_Quad_Edge:
+    if (len == 3) // quadratic edge
+      aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
+    break;
+  case SMDSEntity_Triangle:
+    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 = Min(L1,Min(L2,L3));
     }
-
-    if (aVal <0){
-      return 0.;
+    break;
+  case SMDSEntity_Quadrangle:
+    if (len == 4){ // quadrangles
+      double L1 = getDistance(P( 1 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 4 ));
+      double L4 = getDistance(P( 4 ),P( 1 ));
+      aVal = Min(Min(L1,L2),Min(L3,L4));
     }
-
-    if ( myPrecision >= 0 )
-    {
-      double prec = pow( 10., (double)( myPrecision ) );
-      aVal = floor( aVal * prec + 0.5 ) / prec;
+    break;
+  case SMDSEntity_Quad_Triangle:
+  case SMDSEntity_BiQuad_Triangle:
+    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 = Min(L1,Min(L2,L3));
+    }
+    break;
+  case SMDSEntity_Quad_Quadrangle:
+  case SMDSEntity_BiQuad_Quadrangle:
+    if (len >= 8){ // quadratic quadrangles
+      double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+      double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+      double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
+      double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
+      aVal = Min(Min(L1,L2),Min(L3,L4));
+    }
+    break;
+  case SMDSEntity_Tetra:
+    if (len == 4){ // tetrahedra
+      double L1 = getDistance(P( 1 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 1 ));
+      double L4 = getDistance(P( 1 ),P( 4 ));
+      double L5 = getDistance(P( 2 ),P( 4 ));
+      double L6 = getDistance(P( 3 ),P( 4 ));
+      aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+    }
+    break;
+  case SMDSEntity_Pyramid:
+    if (len == 5){ // pyramid
+      double L1 = getDistance(P( 1 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 4 ));
+      double L4 = getDistance(P( 4 ),P( 1 ));
+      double L5 = getDistance(P( 1 ),P( 5 ));
+      double L6 = getDistance(P( 2 ),P( 5 ));
+      double L7 = getDistance(P( 3 ),P( 5 ));
+      double L8 = getDistance(P( 4 ),P( 5 ));
+
+      aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+      aVal = Min(aVal,Min(L7,L8));
+    }
+    break;
+  case SMDSEntity_Penta:
+    if (len == 6) { // pentahedron
+      double L1 = getDistance(P( 1 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 1 ));
+      double L4 = getDistance(P( 4 ),P( 5 ));
+      double L5 = getDistance(P( 5 ),P( 6 ));
+      double L6 = getDistance(P( 6 ),P( 4 ));
+      double L7 = getDistance(P( 1 ),P( 4 ));
+      double L8 = getDistance(P( 2 ),P( 5 ));
+      double L9 = getDistance(P( 3 ),P( 6 ));
+
+      aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+      aVal = Min(aVal,Min(Min(L7,L8),L9));
+    }
+    break;
+  case SMDSEntity_Hexa:
+    if (len == 8){ // hexahedron
+      double L1 = getDistance(P( 1 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 4 ));
+      double L4 = getDistance(P( 4 ),P( 1 ));
+      double L5 = getDistance(P( 5 ),P( 6 ));
+      double L6 = getDistance(P( 6 ),P( 7 ));
+      double L7 = getDistance(P( 7 ),P( 8 ));
+      double L8 = getDistance(P( 8 ),P( 5 ));
+      double L9 = getDistance(P( 1 ),P( 5 ));
+      double L10= getDistance(P( 2 ),P( 6 ));
+      double L11= getDistance(P( 3 ),P( 7 ));
+      double L12= getDistance(P( 4 ),P( 8 ));
+
+      aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+      aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
+      aVal = Min(aVal,Min(L11,L12));
     }
+    break;
+  case SMDSEntity_Quad_Tetra:
+    if (len == 10){ // quadratic tetrahedron
+      double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
+      double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+      double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
+      double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
+      aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+    }
+    break;
+  case SMDSEntity_Quad_Pyramid:
+    if (len == 13){ // quadratic pyramid
+      double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+      double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+      double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+      double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
+      double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
+      double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
+      aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+      aVal = Min(aVal,Min(L7,L8));
+    }
+    break;
+  case SMDSEntity_Quad_Penta:
+  case SMDSEntity_BiQuad_Penta:
+    if (len >= 15){ // quadratic pentahedron
+      double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+      double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+      double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
+      double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
+      double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
+      double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
+      double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
+      aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+      aVal = Min(aVal,Min(Min(L7,L8),L9));
+    }
+    break;
+  case SMDSEntity_Quad_Hexa:
+  case SMDSEntity_TriQuad_Hexa:
+    if (len >= 20) { // quadratic hexahedron
+      double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
+      double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
+      double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
+      double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
+      double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
+      double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
+      double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
+      double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
+      double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
+      double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
+      double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
+      double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
+      aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+      aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
+      aVal = Min(aVal,Min(L11,L12));
+    }
+    break;
+  case SMDSEntity_Polygon:
+    if ( len > 1 ) {
+      aVal = getDistance( P(1), P( P.size() ));
+      for ( size_t i = 1; i < P.size(); ++i )
+        aVal = Min( aVal, getDistance( P( i ), P( i+1 )));
+    }
+    break;
+  case SMDSEntity_Quad_Polygon:
+    if ( len > 2 ) {
+      aVal = getDistance( P(1), P( P.size() )) + getDistance( P(P.size()), P( P.size()-1 ));
+      for ( size_t i = 1; i < P.size()-1; i += 2 )
+        aVal = Min( aVal, getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 )));
+    }
+    break;
+  case SMDSEntity_Hexagonal_Prism:
+    if (len == 12) { // hexagonal prism
+      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( 5 ));
+      double L5 = getDistance(P( 5 ),P( 6 ));
+      double L6 = getDistance(P( 6 ),P( 1 ));
+
+      double L7 = getDistance(P( 7 ), P( 8 ));
+      double L8 = getDistance(P( 8 ), P( 9 ));
+      double L9 = getDistance(P( 9 ), P( 10 ));
+      double L10= getDistance(P( 10 ),P( 11 ));
+      double L11= getDistance(P( 11 ),P( 12 ));
+      double L12= getDistance(P( 12 ),P( 7 ));
+
+      double L13 = getDistance(P( 1 ),P( 7 ));
+      double L14 = getDistance(P( 2 ),P( 8 ));
+      double L15 = getDistance(P( 3 ),P( 9 ));
+      double L16 = getDistance(P( 4 ),P( 10 ));
+      double L17 = getDistance(P( 5 ),P( 11 ));
+      double L18 = getDistance(P( 6 ),P( 12 ));
+      aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+      aVal = Min(aVal, Min(Min(Min(L7,L8),Min(L9,L10)),Min(L11,L12)));
+      aVal = Min(aVal, Min(Min(Min(L13,L14),Min(L15,L16)),Min(L17,L18)));
+    }
+    break;
+  case SMDSEntity_Polyhedra:
+  {
+  }
+  break;
+  default:
+    return 0;
+  }
 
-    return aVal;
+  if (aVal < 0 ) {
+    return 0.;
+  }
 
+  if ( myPrecision >= 0 )
+  {
+    double prec = pow( 10., (double)( myPrecision ) );
+    aVal = floor( aVal * prec + 0.5 ) / prec;
   }
-  return 0.;
+
+  return aVal;
 }
 
 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
 {
-  // meaningless as it is not quality control functor
+  // meaningless as it is not quality control functor
   return Value;
 }
 
@@ -1701,14 +1816,16 @@ Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
   }
 }
 
-bool Length2D::Value::operator<(const Length2D::Value& x) const{
+bool Length2D::Value::operator<(const Length2D::Value& x) const
+{
   if(myPntId[0] < x.myPntId[0]) return true;
   if(myPntId[0] == x.myPntId[0])
     if(myPntId[1] < x.myPntId[1]) return true;
   return false;
 }
 
-void Length2D::GetValues(TValues& theValues){
+void Length2D::GetValues(TValues& theValues)
+{
   TValues aValues;
   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
   for(; anIter->more(); ){
@@ -1719,10 +1836,10 @@ void Length2D::GetValues(TValues& theValues){
         dynamic_cast<const SMDS_VtkFace*>(anElem);
       // use special nodes iterator
       SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
-      long aNodeId[4];
+      long aNodeId[4] = { 0,0,0,0 };
       gp_Pnt P[4];
 
-      double aLength;
+      double aLength = 0;
       const SMDS_MeshElement* aNode;
       if(anIter->more()){
         aNode = anIter->next();
@@ -1756,7 +1873,7 @@ void Length2D::GetValues(TValues& theValues){
     }
     else {
       SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
-      long aNodeId[2];
+      long aNodeId[2] = {0,0};
       gp_Pnt P[3];
 
       double aLength;
@@ -1791,10 +1908,104 @@ void Length2D::GetValues(TValues& theValues){
   }
 }
 
+//================================================================================
+/*
+  Class       : Deflection2D
+  Description : Functor for calculating number of faces conneted to the edge
+*/
+//================================================================================
+
+double Deflection2D::GetValue( const TSequenceOfXYZ& P )
+{
+  if ( myMesh && P.getElement() )
+  {
+    // get underlying surface
+    if ( myShapeIndex != P.getElement()->getshapeId() )
+    {
+      mySurface.Nullify();
+      myShapeIndex = P.getElement()->getshapeId();
+      const TopoDS_Shape& S =
+        static_cast< const SMESHDS_Mesh* >( myMesh )->IndexToShape( myShapeIndex );
+      if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE )
+      {
+        mySurface = new ShapeAnalysis_Surface( BRep_Tool::Surface( TopoDS::Face( S )));
+
+        GeomLib_IsPlanarSurface isPlaneCheck( mySurface->Surface() );
+        if ( isPlaneCheck.IsPlanar() )
+          myPlane.reset( new gp_Pln( isPlaneCheck.Plan() ));
+        else
+          myPlane.reset();
+      }
+    }
+    // project gravity center to the surface
+    if ( !mySurface.IsNull() )
+    {
+      gp_XYZ gc(0,0,0);
+      gp_XY  uv(0,0);
+      int nbUV = 0;
+      for ( size_t i = 0; i < P.size(); ++i )
+      {
+        gc += P(i+1);
+
+        if ( const SMDS_FacePosition* fPos = dynamic_cast<const SMDS_FacePosition*>
+             ( P.getElement()->GetNode( i )->GetPosition() ))
+        {
+          uv.ChangeCoord(1) += fPos->GetUParameter();
+          uv.ChangeCoord(2) += fPos->GetVParameter();
+          ++nbUV;
+        }
+      }
+      gc /= P.size();
+      if ( nbUV ) uv /= nbUV;
+
+      double maxLen = MaxElementLength2D().GetValue( P );
+      double    tol = 1e-3 * maxLen;
+      double dist;
+      if ( myPlane )
+      {
+        dist = myPlane->Distance( gc );
+        if ( dist < tol )
+          dist = 0;
+      }
+      else
+      {
+        if ( uv.X() != 0 && uv.Y() != 0 ) // faster way
+          mySurface->NextValueOfUV( uv, gc, tol, 0.5 * maxLen );
+        else
+          mySurface->ValueOfUV( gc, tol );
+        dist = mySurface->Gap();
+      }
+      return Round( dist );
+    }
+  }
+  return 0;
+}
+
+void Deflection2D::SetMesh( const SMDS_Mesh* theMesh )
+{
+  NumericalFunctor::SetMesh( dynamic_cast<const SMESHDS_Mesh* >( theMesh ));
+  myShapeIndex = -100;
+  myPlane.reset();
+}
+
+SMDSAbs_ElementType Deflection2D::GetType() const
+{
+  return SMDSAbs_Face;
+}
+
+double Deflection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+  // meaningless as it is not quality control functor
+  return Value;
+}
+
+//================================================================================
 /*
   Class       : MultiConnection
   Description : Functor for calculating number of faces conneted to the edge
 */
+//================================================================================
+
 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
 {
   return 0;
@@ -1815,10 +2026,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;
@@ -1838,7 +2052,7 @@ double MultiConnection2D::GetValue( long theElementId )
       SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
       if (!anIter) break;
 
-      const SMDS_MeshNode *aNode, *aNode0;
+      const SMDS_MeshNode *aNode, *aNode0 = 0;
       TColStd_MapOfInteger aMap, aMapPrev;
 
       for (i = 0; i <= len; i++) {
@@ -1899,14 +2113,17 @@ MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
   }
 }
 
-bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{
+bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const
+{
   if(myPntId[0] < x.myPntId[0]) return true;
   if(myPntId[0] == x.myPntId[0])
     if(myPntId[1] < x.myPntId[1]) return true;
   return false;
 }
 
-void MultiConnection2D::GetValues(MValues& theValues){
+void MultiConnection2D::GetValues(MValues& theValues)
+{
+  if ( !myMesh ) return;
   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
   for(; anIter->more(); ){
     const SMDS_MeshFace* anElem = anIter->next();
@@ -1916,7 +2133,7 @@ void MultiConnection2D::GetValues(MValues& theValues){
         (anElem)->interlacedNodesElemIterator();
     else
       aNodesIter = anElem->nodesIterator();
-    long aNodeId[3];
+    long aNodeId[3] = {0,0,0};
 
     //int aNbConnects=0;
     const SMDS_MeshNode* aNode0;
@@ -1962,10 +2179,13 @@ void MultiConnection2D::GetValues(MValues& theValues){
 
 }
 
+//================================================================================
 /*
   Class       : BallDiameter
   Description : Functor returning diameter of a ball element
 */
+//================================================================================
+
 double BallDiameter::GetValue( long theId )
 {
   double diameter = 0;
@@ -1989,15 +2209,53 @@ SMDSAbs_ElementType BallDiameter::GetType() const
   return SMDSAbs_Ball;
 }
 
+//================================================================================
+/*
+  Class       : NodeConnectivityNumber
+  Description : Functor returning number of elements connected to a node
+*/
+//================================================================================
+
+double NodeConnectivityNumber::GetValue( long theId )
+{
+  double nb = 0;
+
+  if ( const SMDS_MeshNode* node = myMesh->FindNode( theId ))
+  {
+    SMDSAbs_ElementType type;
+    if ( myMesh->NbVolumes() > 0 )
+      type = SMDSAbs_Volume;
+    else if ( myMesh->NbFaces() > 0 )
+      type = SMDSAbs_Face;
+    else if ( myMesh->NbEdges() > 0 )
+      type = SMDSAbs_Edge;
+    else
+      return 0;
+    nb = node->NbInverseElements( type );
+  }
+  return nb;
+}
+
+double NodeConnectivityNumber::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+  return Value;
+}
+
+SMDSAbs_ElementType NodeConnectivityNumber::GetType() const
+{
+  return SMDSAbs_Node;
+}
 
 /*
                             PREDICATES
 */
 
+//================================================================================
 /*
   Class       : BadOrientedVolume
   Description : Predicate bad oriented volumes
 */
+//================================================================================
 
 BadOrientedVolume::BadOrientedVolume()
 {
@@ -2036,7 +2294,7 @@ bool BareBorderVolume::IsSatisfy(long theElementId )
       if ( myTool.IsFreeFace( iF ))
       {
         const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
-        vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
+        std::vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
         if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
           return true;
       }
@@ -2044,9 +2302,11 @@ bool BareBorderVolume::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : BareBorderFace
 */
+//================================================================================
 
 bool BareBorderFace::IsSatisfy(long theElementId )
 {
@@ -2084,9 +2344,11 @@ bool BareBorderFace::IsSatisfy(long theElementId )
   return ok;
 }
 
+//================================================================================
 /*
   Class       : OverConstrainedVolume
 */
+//================================================================================
 
 bool OverConstrainedVolume::IsSatisfy(long theElementId )
 {
@@ -2104,9 +2366,11 @@ bool OverConstrainedVolume::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : OverConstrainedFace
 */
+//================================================================================
 
 bool OverConstrainedFace::IsSatisfy(long theElementId )
 {
@@ -2137,10 +2401,12 @@ bool OverConstrainedFace::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : CoincidentNodes
   Description : Predicate of Coincident nodes
 */
+//================================================================================
 
 CoincidentNodes::CoincidentNodes()
 {
@@ -2167,26 +2433,28 @@ void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
     while ( nIt->more() )
       nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
 
-    list< list< const SMDS_MeshNode*> > nodeGroups;
+    std::list< std::list< const SMDS_MeshNode*> > nodeGroups;
     SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
 
     myCoincidentIDs.Clear();
-    list< list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
+    std::list< std::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();
+      std::list< const SMDS_MeshNode*>& coincNodes = *groupIt;
+      std::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()
 {
@@ -2205,7 +2473,7 @@ bool CoincidentElements::IsSatisfy( long theElementId )
   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() );
+    std::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() )
@@ -2237,10 +2505,12 @@ SMDSAbs_ElementType CoincidentElements3D::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeBorders
   Description : Predicate for free borders
 */
+//================================================================================
 
 FreeBorders::FreeBorders()
 {
@@ -2263,10 +2533,13 @@ SMDSAbs_ElementType FreeBorders::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeEdges
   Description : Predicate for free Edges
 */
+//================================================================================
+
 FreeEdges::FreeEdges()
 {
   myMesh = 0;
@@ -2305,26 +2578,15 @@ bool FreeEdges::IsSatisfy( long theId )
   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
     return false;
 
-  SMDS_ElemIteratorPtr anIter;
-  if ( aFace->IsQuadratic() ) {
-    anIter = dynamic_cast<const SMDS_VtkFace*>
-      (aFace)->interlacedNodesElemIterator();
-  }
-  else {
-    anIter = aFace->nodesIterator();
-  }
+  SMDS_NodeIteratorPtr anIter = aFace->interlacedNodesIterator();
   if ( !anIter )
     return false;
 
   int i = 0, nbNodes = aFace->NbNodes();
   std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
   while( anIter->more() )
-  {
-    const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
-    if ( aNode == 0 )
+    if ( ! ( aNodes[ i++ ] = anIter->next() ))
       return false;
-    aNodes[ i++ ] = aNode;
-  }
   aNodes[ nbNodes ] = aNodes[ 0 ];
 
   for ( i = 0; i < nbNodes; i++ )
@@ -2380,7 +2642,7 @@ void FreeEdges::GetBoreders(TBorders& theBorders)
         interlacedNodesElemIterator();
     else
       aNodesIter = anElem->nodesIterator();
-    long aNodeId[2];
+    long aNodeId[2] = {0,0};
     const SMDS_MeshElement* aNode;
     if(aNodesIter->more()){
       aNode = aNodesIter->next();
@@ -2398,11 +2660,12 @@ void FreeEdges::GetBoreders(TBorders& theBorders)
   }
 }
 
-
+//================================================================================
 /*
   Class       : FreeNodes
   Description : Predicate for free nodes
 */
+//================================================================================
 
 FreeNodes::FreeNodes()
 {
@@ -2429,10 +2692,12 @@ SMDSAbs_ElementType FreeNodes::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeFaces
   Description : Predicate for free faces
 */
+//================================================================================
 
 FreeFaces::FreeFaces()
 {
@@ -2454,19 +2719,21 @@ bool FreeFaces::IsSatisfy( long theId )
 
   int nbNode = aFace->NbNodes();
 
-  // collect volumes check that number of volumss with count equal nbNode not less than 2
-  typedef map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
-  typedef map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
+  // collect volumes to check that number of volumes with count equal nbNode not less than 2
+  typedef std::map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
+  typedef std::map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
   TMapOfVolume mapOfVol;
 
   SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
-  while ( nodeItr->more() ) {
+  while ( nodeItr->more() )
+  {
     const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
     if ( !aNode ) continue;
     SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
-    while ( volItr->more() ) {
+    while ( volItr->more() )
+    {
       SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
-      TItrMapOfVolume itr = mapOfVol.insert(make_pair(aVol, 0)).first;
+      TItrMapOfVolume    itr = mapOfVol.insert( std::make_pair( aVol, 0 )).first;
       (*itr).second++;
     } 
   }
@@ -2476,7 +2743,7 @@ bool FreeFaces::IsSatisfy( long theId )
   for ( ; volItr != volEnd; ++volItr )
     if ( (*volItr).second >= nbNode )
        nbVol++;
-  // face is not free if number of volumes constructed on thier nodes more than one
+  // face is not free if number of volumes constructed on their nodes more than one
   return (nbVol < 2);
 }
 
@@ -2485,10 +2752,12 @@ SMDSAbs_ElementType FreeFaces::GetType() const
   return SMDSAbs_Face;
 }
 
+//================================================================================
 /*
   Class       : LinearOrQuadratic
   Description : Predicate to verify whether a mesh element is linear
 */
+//================================================================================
 
 LinearOrQuadratic::LinearOrQuadratic()
 {
@@ -2519,10 +2788,12 @@ SMDSAbs_ElementType LinearOrQuadratic::GetType() const
   return myType;
 }
 
+//================================================================================
 /*
   Class       : GroupColor
-  Description : Functor for check color of group to whic mesh element belongs to
+  Description : Functor for check color of group to which mesh element belongs to
 */
+//================================================================================
 
 GroupColor::GroupColor()
 {
@@ -2530,7 +2801,7 @@ GroupColor::GroupColor()
 
 bool GroupColor::IsSatisfy( long theId )
 {
-  return (myIDs.find( theId ) != myIDs.end());
+  return myIDs.count( theId );
 }
 
 void GroupColor::SetType( SMDSAbs_ElementType theType )
@@ -2548,16 +2819,15 @@ static bool isEqual( const Quantity_Color& theColor1,
 {
   // tolerance to compare colors
   const double tol = 5*1e-3;
-  return ( fabs( theColor1.Red() - theColor2.Red() ) < tol &&
+  return ( fabs( theColor1.Red()   - theColor2.Red() )   < tol &&
            fabs( theColor1.Green() - theColor2.Green() ) < tol &&
-           fabs( theColor1.Blue() - theColor2.Blue() ) < tol );
+           fabs( theColor1.Blue()  - theColor2.Blue() )  < tol );
 }
 
-
 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
 {
   myIDs.clear();
-  
+
   const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
   if ( !aMesh )
     return;
@@ -2565,20 +2835,24 @@ void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
   int nbGrp = aMesh->GetNbGroups();
   if ( !nbGrp )
     return;
-  
+
   // iterates on groups and find necessary elements ids
-  const std::set<SMESHDS_GroupBase*>& aGroups = aMesh->GetGroups();
-  set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
-  for (; GrIt != aGroups.end(); GrIt++) {
+  const std::set<SMESHDS_GroupBase*>&       aGroups = aMesh->GetGroups();
+  std::set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
+  for (; GrIt != aGroups.end(); GrIt++)
+  {
     SMESHDS_GroupBase* aGrp = (*GrIt);
     if ( !aGrp )
       continue;
     // check type and color of group
-    if ( !isEqual( myColor, aGrp->GetColor() ) )
-      continue;
-    if ( myType != SMDSAbs_All && myType != (SMDSAbs_ElementType)aGrp->GetType() )
+    if ( !isEqual( myColor, aGrp->GetColor() ))
       continue;
 
+    // IPAL52867 (prevent infinite recursion via GroupOnFilter)
+    if ( SMESHDS_GroupOnFilter * gof = dynamic_cast< SMESHDS_GroupOnFilter* >( aGrp ))
+      if ( gof->GetPredicate().get() == this )
+        continue;
+
     SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
     if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
       // add elements IDS into control
@@ -2591,6 +2865,7 @@ void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
 
 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
 {
+  Kernel_Utils::Localizer loc;
   TCollection_AsciiString aStr = theStr;
   aStr.RemoveAll( ' ' );
   aStr.RemoveAll( '\t' );
@@ -2611,6 +2886,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();
@@ -2619,10 +2895,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()
 {
@@ -2645,43 +2923,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 == 9 ) : 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 == 27 ): aNbNode == 8));
-    else if ( myGeomType == SMDSGeom_HEXAGONAL_PRISM )
-      isOk = (anElem->GetEntityType() == SMDSEntity_Hexagonal_Prism );
-     else if ( myGeomType == SMDSGeom_POLYHEDRA )
-      isOk = anElem->IsPoly();
-    break;
-    default: break;
-  }
+  bool isOk = ( anElem->GetGeomType() == myGeomType );
   return isOk;
 }
 
@@ -2705,12 +2947,214 @@ 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;
+  if ( myType == SMDSAbs_Node )
+    return myMesh->FindNode( theId );
+  const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
+  return ( anElem &&
+           myEntityType == anElem->GetEntityType() );
+}
+
+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 ConnectedElements
+ */
+//================================================================================
+
+ConnectedElements::ConnectedElements():
+  myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {}
+
+SMDSAbs_ElementType ConnectedElements::GetType() const
+{ return myType; }
+
+int ConnectedElements::GetNode() const
+{ return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
+
+std::vector<double> ConnectedElements::GetPoint() const
+{ return myXYZ; }
+
+void ConnectedElements::clearOkIDs()
+{ myOkIDsReady = false; myOkIDs.clear(); }
+
+void ConnectedElements::SetType( SMDSAbs_ElementType theType )
+{
+  if ( myType != theType || myMeshModifTracer.IsMeshModified() )
+    clearOkIDs();
+  myType = theType;
+}
+
+void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh )
+{
+  myMeshModifTracer.SetMesh( theMesh );
+  if ( myMeshModifTracer.IsMeshModified() )
+  {
+    clearOkIDs();
+    if ( !myXYZ.empty() )
+      SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh
+  }
+}
+
+void ConnectedElements::SetNode( int nodeID )
+{
+  myNodeID = nodeID;
+  myXYZ.clear();
+
+  bool isSameDomain = false;
+  if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() )
+    if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID ))
+    {
+      SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType );
+      while ( !isSameDomain && eIt->more() )
+        isSameDomain = IsSatisfy( eIt->next()->GetID() );
+    }
+  if ( !isSameDomain )
+    clearOkIDs();
+}
+
+void ConnectedElements::SetPoint( double x, double y, double z )
+{
+  myXYZ.resize(3);
+  myXYZ[0] = x;
+  myXYZ[1] = y;
+  myXYZ[2] = z;
+  myNodeID = 0;
+
+  bool isSameDomain = false;
+
+  // find myNodeID by myXYZ if possible
+  if ( myMeshModifTracer.GetMesh() )
+  {
+    SMESHUtils::Deleter<SMESH_ElementSearcher> searcher
+      ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
+
+    std::vector< const SMDS_MeshElement* > foundElems;
+    searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
+
+    if ( !foundElems.empty() )
+    {
+      myNodeID = foundElems[0]->GetNode(0)->GetID();
+      if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() )
+        isSameDomain = IsSatisfy( foundElems[0]->GetID() );
+    }
+  }
+  if ( !isSameDomain )
+    clearOkIDs();
+}
+
+bool ConnectedElements::IsSatisfy( long theElementId )
+{
+  // Here we do NOT check if the mesh has changed, we do it in Set...() only!!!
+
+  if ( !myOkIDsReady )
+  {
+    if ( !myMeshModifTracer.GetMesh() )
+      return false;
+    const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID );
+    if ( !node0 )
+      return false;
+
+    std::list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
+    std::set< int > checkedNodeIDs;
+    // algo:
+    // foreach node in nodeQueue:
+    //   foreach element sharing a node:
+    //     add ID of an element of myType to myOkIDs;
+    //     push all element nodes absent from checkedNodeIDs to nodeQueue;
+    while ( !nodeQueue.empty() )
+    {
+      const SMDS_MeshNode* node = nodeQueue.front();
+      nodeQueue.pop_front();
+
+      // loop on elements sharing the node
+      SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
+      while ( eIt->more() )
+      {
+        // keep elements of myType
+        const SMDS_MeshElement* element = eIt->next();
+        if ( element->GetType() == myType )
+          myOkIDs.insert( myOkIDs.end(), element->GetID() );
+
+        // enqueue nodes of the element
+        SMDS_ElemIteratorPtr nIt = element->nodesIterator();
+        while ( nIt->more() )
+        {
+          const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() );
+          if ( checkedNodeIDs.insert( n->GetID() ).second )
+            nodeQueue.push_back( n );
+        }
+      }
+    }
+    if ( myType == SMDSAbs_Node )
+      std::swap( myOkIDs, checkedNodeIDs );
+
+    size_t totalNbElems = myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType );
+    if ( myOkIDs.size() == totalNbElems )
+      myOkIDs.clear();
+
+    myOkIDsReady = true;
+  }
+
+  return myOkIDs.empty() ? true : myOkIDs.count( theElementId );
+}
+
 //================================================================================
 /*!
  * \brief Class CoplanarFaces
  */
 //================================================================================
 
+namespace
+{
+  inline bool isLessAngle( const gp_Vec& v1, const gp_Vec& v2, const double cos )
+  {
+    double dot = v1 * v2; // cos * |v1| * |v2|
+    double l1  = v1.SquareMagnitude();
+    double l2  = v2.SquareMagnitude();
+    return (( dot * cos >= 0 ) && 
+            ( dot * dot ) / l1 / l2 >= ( cos * cos ));
+  }
+}
 CoplanarFaces::CoplanarFaces()
   : myFaceID(0), myToler(0)
 {
@@ -2722,7 +3166,7 @@ void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
   {
     // Build a set of coplanar face ids
 
-    myCoplanarIDs.clear();
+    myCoplanarIDs.Clear();
 
     if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler )
       return;
@@ -2736,41 +3180,48 @@ void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
     if (!normOK)
       return;
 
-    const double radianTol = myToler * M_PI / 180.;
-    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 cosTol = Cos( myToler * M_PI / 180. );
+    NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks;
+
+    std::list< std::pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
+    faceQueue.push_back( std::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.Add( SMESH_TLink( n1, n2 )))
+          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 || isLessAngle( myNorm, norm, cosTol))
+            {
+              myCoplanarIDs.Add( f->GetID() );
+              faceQueue.push_back( std::make_pair( f, norm ));
+            }
           }
-          faceQueue.insert( faceQueue.end(), neighborFaces.begin(), neighborFaces.end() );
         }
       }
-      faceQueue.pop_front();
     }
   }
 }
 bool CoplanarFaces::IsSatisfy( long theElementId )
 {
-  return myCoplanarIDs.count( theElementId );
+  return myCoplanarIDs.Contains( theElementId );
 }
 
 /*
 *Class       : RangeOfIds
+ *Class       : RangeOfIds
   *Description : Predicate for Range of Ids.
   *              Range may be specified with two ways.
   *              1. Using AddToRange method
@@ -2894,11 +3345,13 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
   myIds.Clear();
 
   TCollection_AsciiString aStr = theStr;
+  for ( int i = 1; i <= aStr.Length(); ++i )
+  {
+    char c = aStr.Value( i );
+    if ( !isdigit( c ) && c != ',' && c != '-' )
+      aStr.SetValue( i, ',');
+  }
   aStr.RemoveAll( ' ' );
-  aStr.RemoveAll( '\t' );
-
-  for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
-    aStr.Remove( aPos, 2 );
 
   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
   int i = 1;
@@ -3169,6 +3622,31 @@ bool LogicalOR::IsSatisfy( long theId )
                               FILTER
 */
 
+// #ifdef WITH_TBB
+// #include <tbb/parallel_for.h>
+// #include <tbb/enumerable_thread_specific.h>
+
+// namespace Parallel
+// {
+//   typedef tbb::enumerable_thread_specific< TIdSequence > TIdSeq;
+
+//   struct Predicate
+//   {
+//     const SMDS_Mesh* myMesh;
+//     PredicatePtr     myPredicate;
+//     TIdSeq &         myOKIds;
+//     Predicate( const SMDS_Mesh* m, PredicatePtr p, TIdSeq & ids ):
+//       myMesh(m), myPredicate(p->Duplicate()), myOKIds(ids) {}
+//     void operator() ( const tbb::blocked_range<size_t>& r ) const
+//     {
+//       for ( size_t i = r.begin(); i != r.end(); ++i )
+//         if ( myPredicate->IsSatisfy( i ))
+//           myOKIds.local().push_back();
+//     }
+//   }
+// }
+// #endif
+
 Filter::Filter()
 {}
 
@@ -3349,7 +3827,7 @@ bool ManifoldPart::process()
       myMapIds.Add( aFaceId );
     }
 
-    if ( fi == ( myAllFacePtr.size() - 1 ) )
+    if ( fi == int( myAllFacePtr.size() - 1 ))
       fi = 0;
   } // end run on vector of faces
   return !myMapIds.IsEmpty();
@@ -3551,35 +4029,82 @@ void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
   }
 }
 
-
 /*
-   ElementsOnSurface
+  Class       : BelongToMeshGroup
+  Description : Verify whether a mesh element is included into a mesh group
 */
-
-ElementsOnSurface::ElementsOnSurface()
+BelongToMeshGroup::BelongToMeshGroup(): myGroup( 0 )
 {
-  myMesh = 0;
-  myIds.Clear();
-  myType = SMDSAbs_All;
-  mySurf.Nullify();
-  myToler = Precision::Confusion();
-  myUseBoundaries = false;
 }
 
-ElementsOnSurface::~ElementsOnSurface()
+void BelongToMeshGroup::SetGroup( SMESHDS_GroupBase* g )
 {
-  myMesh = 0;
+  myGroup = g;
 }
 
-void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
+void BelongToMeshGroup::SetStoreName( const std::string& sn )
 {
-  if ( myMesh == theMesh )
-    return;
-  myMesh = theMesh;
-  process();
+  myStoreName = sn;
 }
 
-bool ElementsOnSurface::IsSatisfy( long theElementId )
+void BelongToMeshGroup::SetMesh( const SMDS_Mesh* theMesh )
+{
+  if ( myGroup && myGroup->GetMesh() != theMesh )
+  {
+    myGroup = 0;
+  }
+  if ( !myGroup && !myStoreName.empty() )
+  {
+    if ( const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh))
+    {
+      const std::set<SMESHDS_GroupBase*>& grps = aMesh->GetGroups();
+      std::set<SMESHDS_GroupBase*>::const_iterator g = grps.begin();
+      for ( ; g != grps.end() && !myGroup; ++g )
+        if ( *g && myStoreName == (*g)->GetStoreName() )
+          myGroup = *g;
+    }
+  }
+  if ( myGroup )
+  {
+    myGroup->IsEmpty(); // make GroupOnFilter update its predicate
+  }
+}
+
+bool BelongToMeshGroup::IsSatisfy( long theElementId )
+{
+  return myGroup ? myGroup->Contains( theElementId ) : false;
+}
+
+SMDSAbs_ElementType BelongToMeshGroup::GetType() const
+{
+  return myGroup ? myGroup->GetType() : SMDSAbs_All;
+}
+
+//================================================================================
+//  ElementsOnSurface
+//================================================================================
+
+ElementsOnSurface::ElementsOnSurface()
+{
+  myIds.Clear();
+  myType = SMDSAbs_All;
+  mySurf.Nullify();
+  myToler = Precision::Confusion();
+  myUseBoundaries = false;
+}
+
+ElementsOnSurface::~ElementsOnSurface()
+{
+}
+
+void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
+{
+  myMeshModifTracer.SetMesh( theMesh );
+  if ( myMeshModifTracer.IsMeshModified())
+    process();
+}
+
+bool ElementsOnSurface::IsSatisfy( long theElementId )
 {
   return myIds.Contains( theElementId );
 }
@@ -3631,32 +4156,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 )
@@ -3709,33 +4216,96 @@ bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
 }
 
 
-/*
-  ElementsOnShape
-*/
+//================================================================================
+//  ElementsOnShape
+//================================================================================
 
-ElementsOnShape::ElementsOnShape()
-  : //myMesh(0),
-    myType(SMDSAbs_All),
-    myToler(Precision::Confusion()),
-    myAllNodesFlag(false)
+namespace {
+  const int theIsCheckedFlag = 0x0000100;
+}
+
+struct ElementsOnShape::Classifier
+{
+  Classifier() { mySolidClfr = 0; myFlags = 0; }
+  ~Classifier();
+  void Init(const TopoDS_Shape& s, double tol, const Bnd_B3d* box = 0 );
+  bool IsOut(const gp_Pnt& p)        { return SetChecked( true ), (this->*myIsOutFun)( p ); }
+  TopAbs_ShapeEnum ShapeType() const { return myShape.ShapeType(); }
+  const TopoDS_Shape& Shape() const  { return myShape; }
+  const Bnd_B3d* GetBndBox() const   { return & myBox; }
+  bool IsChecked()                   { return myFlags & theIsCheckedFlag; }
+  bool IsSetFlag( int flag ) const   { return myFlags & flag; }
+  void SetChecked( bool is ) { is ? SetFlag( theIsCheckedFlag ) : UnsetFlag( theIsCheckedFlag ); }
+  void SetFlag  ( int flag ) { myFlags |= flag; }
+  void UnsetFlag( int flag ) { myFlags &= ~flag; }
+
+private:
+  bool isOutOfSolid (const gp_Pnt& p);
+  bool isOutOfBox   (const gp_Pnt& p);
+  bool isOutOfFace  (const gp_Pnt& p);
+  bool isOutOfEdge  (const gp_Pnt& p);
+  bool isOutOfVertex(const gp_Pnt& p);
+  bool isBox        (const TopoDS_Shape& s);
+
+  bool (Classifier::*          myIsOutFun)(const gp_Pnt& p);
+  BRepClass3d_SolidClassifier* mySolidClfr; // ptr because of a run-time forbidden copy-constructor
+  Bnd_B3d                      myBox;
+  GeomAPI_ProjectPointOnSurf   myProjFace;
+  GeomAPI_ProjectPointOnCurve  myProjEdge;
+  gp_Pnt                       myVertexXYZ;
+  TopoDS_Shape                 myShape;
+  double                       myTol;
+  int                          myFlags;
+};
+
+struct ElementsOnShape::OctreeClassifier : public SMESH_Octree
+{
+  OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers );
+  OctreeClassifier( const OctreeClassifier*                           otherTree,
+                    const std::vector< ElementsOnShape::Classifier >& clsOther,
+                    std::vector< ElementsOnShape::Classifier >&       cls );
+  void GetClassifiersAtPoint( const gp_XYZ& p,
+                              std::vector< ElementsOnShape::Classifier* >& classifiers );
+protected:
+  OctreeClassifier() {}
+  SMESH_Octree* newChild() const { return new OctreeClassifier; }
+  void          buildChildrenData();
+  Bnd_B3d*      buildRootBox();
+
+  std::vector< ElementsOnShape::Classifier* > myClassifiers;
+};
+
+
+ElementsOnShape::ElementsOnShape():
+  myOctree(0),
+  myType(SMDSAbs_All),
+  myToler(Precision::Confusion()),
+  myAllNodesFlag(false)
 {
-  myCurShapeType = TopAbs_SHAPE;
 }
 
 ElementsOnShape::~ElementsOnShape()
 {
+  clearClassifiers();
 }
 
-void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
-{
-  myMeshModifTracer.SetMesh( theMesh );
-  if ( myMeshModifTracer.IsMeshModified())
-    SetShape(myShape, myType);
-}
-
-bool ElementsOnShape::IsSatisfy (long theElementId)
+Predicate* ElementsOnShape::clone() const
 {
-  return myIds.Contains(theElementId);
+  ElementsOnShape* cln = new ElementsOnShape();
+  cln->SetAllNodes ( myAllNodesFlag );
+  cln->SetTolerance( myToler );
+  cln->SetMesh     ( myMeshModifTracer.GetMesh() );
+  cln->myShape = myShape; // avoid creation of myClassifiers
+  cln->SetShape    ( myShape, myType );
+  cln->myClassifiers.resize( myClassifiers.size() );
+  for ( size_t i = 0; i < myClassifiers.size(); ++i )
+    cln->myClassifiers[ i ].Init( BRepBuilderAPI_Copy( myClassifiers[ i ].Shape()),
+                                  myToler, myClassifiers[ i ].GetBndBox() );
+  if ( myOctree ) // copy myOctree
+  {
+    cln->myOctree = new OctreeClassifier( myOctree, myClassifiers, cln->myClassifiers );
+  }
+  return cln;
 }
 
 SMDSAbs_ElementType ElementsOnShape::GetType() const
@@ -3758,202 +4328,814 @@ 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)
+{
+  myMeshModifTracer.SetMesh( theMesh );
+  if ( myMeshModifTracer.IsMeshModified())
+  {
+    size_t nbNodes = theMesh ? theMesh->NbNodes() : 0;
+    if ( myNodeIsChecked.size() == nbNodes )
+    {
+      std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
+    }
+    else
+    {
+      SMESHUtils::FreeVector( myNodeIsChecked );
+      SMESHUtils::FreeVector( myNodeIsOut );
+      myNodeIsChecked.resize( nbNodes, false );
+      myNodeIsOut.resize( nbNodes );
+    }
+  }
+}
+
+bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut )
+{
+  if ( n->GetID() >= (int) myNodeIsChecked.size() ||
+       !myNodeIsChecked[ n->GetID() ])
+    return false;
+
+  isOut = myNodeIsOut[ n->GetID() ];
+  return true;
+}
+
+void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool  isOut )
+{
+  if ( n->GetID() < (int) myNodeIsChecked.size() )
+  {
+    myNodeIsChecked[ n->GetID() ] = true;
+    myNodeIsOut    [ n->GetID() ] = isOut;
   }
 }
 
 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
                                 const SMDSAbs_ElementType theType)
 {
-  myType = theType;
+  bool shapeChanges = ( myShape != theShape );
+  myType  = theType;
   myShape = theShape;
-  myIds.Clear();
+  if ( myShape.IsNull() ) return;
 
-  const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh();
-  
-  if ( !myMesh ) return;
+  if ( shapeChanges )
+  {
+    // find most complex shapes
+    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 )
+    {
+      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() );
+    }
 
-  switch (myType)
+    clearClassifiers();
+    myClassifiers.resize( shapesMap.Extent() );
+    for ( int i = 0; i < shapesMap.Extent(); ++i )
+      myClassifiers[ i ].Init( shapesMap( i+1 ), myToler );
+  }
+
+  if ( theType == SMDSAbs_Node )
   {
-  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;
+    SMESHUtils::FreeVector( myNodeIsChecked );
+    SMESHUtils::FreeVector( myNodeIsOut );
   }
+  else
+  {
+    std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
+  }
+}
+
+void ElementsOnShape::clearClassifiers()
+{
+  // for ( size_t i = 0; i < myClassifiers.size(); ++i )
+  //   delete myClassifiers[ i ];
+  myClassifiers.clear();
 
-  myShapesMap.Clear();
-  addShape(myShape);
+  delete myOctree;
+  myOctree = 0;
 }
 
-void ElementsOnShape::addShape (const TopoDS_Shape& theShape)
+bool ElementsOnShape::IsSatisfy( long elemId )
 {
-  if (theShape.IsNull() || myMeshModifTracer.GetMesh() == 0)
-    return;
+  if ( myClassifiers.empty() )
+    return false;
+
+  const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh();
+  if ( myType == SMDSAbs_Node )
+    return IsSatisfy( mesh->FindNode( elemId ));
+  return IsSatisfy( mesh->FindElement( elemId ));
+}
+
+bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem)
+{
+  if ( !elem )
+    return false;
+
+  bool isSatisfy = myAllNodesFlag, isNodeOut;
 
-  if (!myShapesMap.Add(theShape)) return;
+  gp_XYZ centerXYZ (0, 0, 0);
 
-  myCurShapeType = theShape.ShapeType();
-  switch (myCurShapeType)
+  if ( !myOctree && myClassifiers.size() > 5 )
+  {
+    myWorkClassifiers.resize( myClassifiers.size() );
+    for ( size_t i = 0; i < myClassifiers.size(); ++i )
+      myWorkClassifiers[ i ] = & myClassifiers[ i ];
+    myOctree = new OctreeClassifier( myWorkClassifiers );
+  }
+
+  SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
+  while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
   {
-  case TopAbs_COMPOUND:
-  case TopAbs_COMPSOLID:
-  case TopAbs_SHELL:
-  case TopAbs_WIRE:
+    SMESH_TNodeXYZ aPnt( aNodeItr->next() );
+    centerXYZ += aPnt;
+
+    isNodeOut = true;
+    if ( !getNodeIsOut( aPnt._node, isNodeOut ))
     {
-      TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
-      for (; anIt.More(); anIt.Next()) addShape(anIt.Value());
+      if ( myOctree )
+      {
+        myWorkClassifiers.clear();
+        myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers );
+
+        for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
+          myWorkClassifiers[i]->SetChecked( false );
+
+        for ( size_t i = 0; i < myWorkClassifiers.size() && isNodeOut; ++i )
+          if ( !myWorkClassifiers[i]->IsChecked() )
+            isNodeOut = myWorkClassifiers[i]->IsOut( aPnt );
+      }
+      else
+      {
+        for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i )
+          isNodeOut = myClassifiers[i].IsOut( aPnt );
+      }
+      setNodeIsOut( aPnt._node, isNodeOut );
     }
-    break;
+    isSatisfy = !isNodeOut;
+  }
+
+  // Check the center point for volumes MantisBug 0020168
+  if ( isSatisfy &&
+       myAllNodesFlag &&
+       myClassifiers[0].ShapeType() == TopAbs_SOLID )
+  {
+    centerXYZ /= elem->NbNodes();
+    isSatisfy = false;
+    if ( myOctree )
+      for ( size_t i = 0; i < myWorkClassifiers.size() && !isSatisfy; ++i )
+        isSatisfy = ! myWorkClassifiers[i]->IsOut( centerXYZ );
+    else
+      for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i )
+        isSatisfy = ! myClassifiers[i].IsOut( centerXYZ );
+  }
+
+  return isSatisfy;
+}
+
+bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node,
+                                 TopoDS_Shape*        okShape)
+{
+  if ( !node )
+    return false;
+
+  if ( !myOctree && myClassifiers.size() > 5 )
+  {
+    myWorkClassifiers.resize( myClassifiers.size() );
+    for ( size_t i = 0; i < myClassifiers.size(); ++i )
+      myWorkClassifiers[ i ] = & myClassifiers[ i ];
+    myOctree = new OctreeClassifier( myWorkClassifiers );
+  }
+
+  bool isNodeOut = true;
+
+  if ( okShape || !getNodeIsOut( node, isNodeOut ))
+  {
+    SMESH_NodeXYZ aPnt = node;
+    if ( myOctree )
+    {
+      myWorkClassifiers.clear();
+      myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers );
+
+      for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
+        myWorkClassifiers[i]->SetChecked( false );
+
+      for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
+        if ( !myWorkClassifiers[i]->IsChecked() &&
+             !myWorkClassifiers[i]->IsOut( aPnt ))
+        {
+          isNodeOut = false;
+          if ( okShape )
+            *okShape = myWorkClassifiers[i]->Shape();
+          break;
+        }
+    }
+    else
+    {
+      for ( size_t i = 0; i < myClassifiers.size(); ++i )
+        if ( !myClassifiers[i].IsOut( aPnt ))
+        {
+          isNodeOut = false;
+          if ( okShape )
+            *okShape = myWorkClassifiers[i]->Shape();
+          break;
+        }
+    }
+    setNodeIsOut( node, isNodeOut );
+  }
+
+  return !isNodeOut;
+}
+
+void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape,
+                                        double              theTol,
+                                        const Bnd_B3d*      theBox )
+{
+  myShape = theShape;
+  myTol   = theTol;
+  myFlags = 0;
+
+  bool isShapeBox = false;
+  switch ( myShape.ShapeType() )
+  {
   case TopAbs_SOLID:
+  {
+    if (( isShapeBox = isBox( theShape )))
     {
-      myCurSC.Load(theShape);
-      process();
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfBox;
     }
-    break;
-  case TopAbs_FACE:
+    else
     {
-      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();
+      mySolidClfr = new BRepClass3d_SolidClassifier(theShape);
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid;
     }
     break;
+  }
+  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::Classifier::isOutOfFace;
+    break;
+  }
   case TopAbs_EDGE:
-    {
-      TopoDS_Edge anEdge = TopoDS::Edge(theShape);
-      Standard_Real u1, u2;
-      Handle(Geom_Curve) curve = BRep_Tool::Curve(anEdge, u1, u2);
-      myCurProjEdge.Init(curve, u1, u2);
-      process();
-    }
+  {
+    Standard_Real u1, u2;
+    Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2);
+    myProjEdge.Init(curve, u1, u2);
+    myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
     break;
+  }
   case TopAbs_VERTEX:
-    {
-      TopoDS_Vertex aV = TopoDS::Vertex(theShape);
-      myCurPnt = BRep_Tool::Pnt(aV);
-      process();
-    }
+  {
+    myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
+    myIsOutFun = & ElementsOnShape::Classifier::isOutOfVertex;
     break;
+  }
   default:
-    break;
+    throw SALOME_Exception("Programmer error in usage of ElementsOnShape::Classifier");
+  }
+
+  if ( !isShapeBox )
+  {
+    if ( theBox )
+    {
+      myBox = *theBox;
+    }
+    else
+    {
+      Bnd_Box box;
+      BRepBndLib::Add( myShape, box );
+      myBox.Clear();
+      myBox.Add( box.CornerMin() );
+      myBox.Add( box.CornerMax() );
+      gp_XYZ halfSize = 0.5 * ( box.CornerMax().XYZ() - box.CornerMin().XYZ() );
+      for ( int iDim = 1; iDim <= 3; ++iDim )
+      {
+        double x = halfSize.Coord( iDim );
+        halfSize.SetCoord( iDim, x + Max( myTol, 1e-2 * x ));
+      }
+      myBox.SetHSize( halfSize );
+    }
   }
 }
 
-void ElementsOnShape::process()
+ElementsOnShape::Classifier::~Classifier()
 {
-  const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh();
-  if (myShape.IsNull() || myMesh == 0)
-    return;
+  delete mySolidClfr; mySolidClfr = 0;
+}
 
-  SMDS_ElemIteratorPtr anIter = myMesh->elementsIterator(myType);
-  while (anIter->more())
-    process(anIter->next());
+bool ElementsOnShape::Classifier::isOutOfSolid (const gp_Pnt& p)
+{
+  mySolidClfr->Perform( p, myTol );
+  return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON );
 }
 
-void ElementsOnShape::process (const SMDS_MeshElement* theElemPtr)
+bool ElementsOnShape::Classifier::isOutOfBox (const gp_Pnt& p)
 {
-  if (myShape.IsNull())
-    return;
+  return myBox.IsOut( p.XYZ() );
+}
 
-  SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
-  bool isSatisfy = myAllNodesFlag;
+bool ElementsOnShape::Classifier::isOutOfFace  (const gp_Pnt& p)
+{
+  myProjFace.Perform( p );
+  if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
+  {
+    // check relatively to the face
+    Standard_Real 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;
+}
 
-  gp_XYZ centerXYZ (0, 0, 0);
+bool ElementsOnShape::Classifier::isOutOfEdge  (const gp_Pnt& p)
+{
+  myProjEdge.Perform( p );
+  return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
+}
 
-  while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
+bool ElementsOnShape::Classifier::isOutOfVertex(const gp_Pnt& p)
+{
+  return ( myVertexXYZ.Distance( p ) > myTol );
+}
+
+bool ElementsOnShape::Classifier::isBox (const TopoDS_Shape& theShape)
+{
+  TopTools_IndexedMapOfShape vMap;
+  TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
+  if ( vMap.Extent() != 8 )
+    return false;
+
+  myBox.Clear();
+  for ( int i = 1; i <= 8; ++i )
+    myBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))).XYZ() );
+
+  gp_XYZ pMin = myBox.CornerMin(), pMax = myBox.CornerMax();
+  for ( int i = 1; i <= 8; ++i )
   {
-    SMESH_TNodeXYZ aPnt ( aNodeItr->next() );
-    centerXYZ += aPnt;
+    gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i )));
+    for ( int iC = 1; iC <= 3; ++ iC )
+    {
+      double d1 = Abs( pMin.Coord( iC ) - p.Coord( iC ));
+      double d2 = Abs( pMax.Coord( iC ) - p.Coord( iC ));
+      if ( Min( d1, d2 ) > myTol )
+        return false;
+    }
+  }
+  myBox.Enlarge( myTol );
+  return true;
+}
+
+ElementsOnShape::
+OctreeClassifier::OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers )
+  :SMESH_Octree( new SMESH_TreeLimit )
+{
+  myClassifiers = classifiers;
+  compute();
+}
+
+ElementsOnShape::
+OctreeClassifier::OctreeClassifier( const OctreeClassifier*                           otherTree,
+                                    const std::vector< ElementsOnShape::Classifier >& clsOther,
+                                    std::vector< ElementsOnShape::Classifier >&       cls )
+  :SMESH_Octree( new SMESH_TreeLimit )
+{
+  myBox = new Bnd_B3d( *otherTree->getBox() );
+
+  if (( myIsLeaf = otherTree->isLeaf() ))
+  {
+    myClassifiers.resize( otherTree->myClassifiers.size() );
+    for ( size_t i = 0; i < otherTree->myClassifiers.size(); ++i )
+    {
+      int ind = otherTree->myClassifiers[i] - & clsOther[0];
+      myClassifiers[ i ] = & cls[ ind ];
+    }
+  }
+  else if ( otherTree->myChildren )
+  {
+    myChildren = new SMESH_Tree< Bnd_B3d, 8 > * [ 8 ];
+    for ( int i = 0; i < nbChildren(); i++ )
+      myChildren[i] =
+        new OctreeClassifier( static_cast<const OctreeClassifier*>( otherTree->myChildren[i]),
+                              clsOther, cls );
+  }
+}
+
+void ElementsOnShape::
+OctreeClassifier::GetClassifiersAtPoint( const gp_XYZ& point,
+                                         std::vector< ElementsOnShape::Classifier* >& result )
+{
+  if ( getBox()->IsOut( point ))
+    return;
+
+  if ( isLeaf() )
+  {
+    for ( size_t i = 0; i < myClassifiers.size(); ++i )
+      if ( !myClassifiers[i]->GetBndBox()->IsOut( point ))
+        result.push_back( myClassifiers[i] );
+  }
+  else
+  {
+    for (int i = 0; i < nbChildren(); i++)
+      ((OctreeClassifier*) myChildren[i])->GetClassifiersAtPoint( point, result );
+  }
+}
+
+void ElementsOnShape::OctreeClassifier::buildChildrenData()
+{
+  // distribute myClassifiers among myChildren
+
+  const int childFlag[8] = { 0x0000001,
+                             0x0000002,
+                             0x0000004,
+                             0x0000008,
+                             0x0000010,
+                             0x0000020,
+                             0x0000040,
+                             0x0000080 };
+  int nbInChild[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
 
-    switch (myCurShapeType)
+  for ( size_t i = 0; i < myClassifiers.size(); ++i )
+  {
+    for ( int j = 0; j < nbChildren(); j++ )
     {
-    case TopAbs_SOLID:
+      if ( !myClassifiers[i]->GetBndBox()->IsOut( *myChildren[j]->getBox() ))
       {
-        myCurSC.Perform(aPnt, myToler);
-        isSatisfy = (myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON);
+        myClassifiers[i]->SetFlag( childFlag[ j ]);
+        ++nbInChild[ j ];
       }
-      break;
-    case TopAbs_FACE:
+    }
+  }
+
+  for ( int j = 0; j < nbChildren(); j++ )
+  {
+    OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ j ]);
+    child->myClassifiers.resize( nbInChild[ j ]);
+    for ( size_t i = 0; nbInChild[ j ] && i < myClassifiers.size(); ++i )
+    {
+      if ( myClassifiers[ i ]->IsSetFlag( childFlag[ j ]))
       {
-        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);
-        }
+        --nbInChild[ j ];
+        child->myClassifiers[ nbInChild[ j ]] = myClassifiers[ i ];
+        myClassifiers[ i ]->UnsetFlag( childFlag[ j ]);
       }
-      break;
-    case TopAbs_EDGE:
-      {
-        myCurProjEdge.Perform(aPnt);
-        isSatisfy = (myCurProjEdge.NbPoints() > 0 && myCurProjEdge.LowerDistance() <= myToler);
+    }
+  }
+  SMESHUtils::FreeVector( myClassifiers );
+
+  // define if a child isLeaf()
+  for ( int i = 0; i < nbChildren(); i++ )
+  {
+    OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ i ]);
+    child->myIsLeaf = ( child->myClassifiers.size() <= 5  );
+  }
+}
+
+Bnd_B3d* ElementsOnShape::OctreeClassifier::buildRootBox()
+{
+  Bnd_B3d* box = new Bnd_B3d;
+  for ( size_t i = 0; i < myClassifiers.size(); ++i )
+    box->Add( *myClassifiers[i]->GetBndBox() );
+  return box;
+}
+
+/*
+  Class       : BelongToGeom
+  Description : Predicate for verifying whether entity belongs to
+                specified geometrical support
+*/
+
+BelongToGeom::BelongToGeom()
+  : myMeshDS(NULL),
+    myType(SMDSAbs_NbElementTypes),
+    myIsSubshape(false),
+    myTolerance(Precision::Confusion())
+{}
+
+Predicate* BelongToGeom::clone() const
+{
+  BelongToGeom* cln = new BelongToGeom( *this );
+  cln->myElementsOnShapePtr.reset( static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ));
+  return cln;
+}
+
+void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
+{
+  if ( myMeshDS != theMesh )
+  {
+    myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
+    init();
+  }
+}
+
+void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
+{
+  if ( myShape != theShape )
+  {
+    myShape = theShape;
+    init();
+  }
+}
+
+static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
+                        const TopoDS_Shape&               theShape)
+{
+  if (theMap.Contains(theShape)) return true;
+
+  if (theShape.ShapeType() == TopAbs_COMPOUND ||
+      theShape.ShapeType() == TopAbs_COMPSOLID)
+  {
+    TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
+    for (; anIt.More(); anIt.Next())
+    {
+      if (!IsSubShape(theMap, anIt.Value())) {
+        return false;
       }
-      break;
-    case TopAbs_VERTEX:
+    }
+    return true;
+  }
+
+  return false;
+}
+
+void BelongToGeom::init()
+{
+  if ( !myMeshDS || myShape.IsNull() ) return;
+
+  // is sub-shape of main shape?
+  TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
+  if (aMainShape.IsNull()) {
+    myIsSubshape = false;
+  }
+  else {
+    TopTools_IndexedMapOfShape aMap;
+    TopExp::MapShapes( aMainShape, aMap );
+    myIsSubshape = IsSubShape( aMap, myShape );
+    if ( myIsSubshape )
+    {
+      aMap.Clear();
+      TopExp::MapShapes( myShape, aMap );
+      mySubShapesIDs.Clear();
+      for ( int i = 1; i <= aMap.Extent(); ++i )
       {
-        isSatisfy = (myCurPnt.Distance(aPnt) <= myToler);
+        int subID = myMeshDS->ShapeToIndex( aMap( i ));
+        if ( subID > 0 )
+          mySubShapesIDs.Add( subID );
       }
-      break;
-    default:
+    }
+  }
+
+  //if (!myIsSubshape) // to be always ready to check an element not bound to geometry
+  {
+    if ( !myElementsOnShapePtr )
+      myElementsOnShapePtr.reset( new ElementsOnShape() );
+    myElementsOnShapePtr->SetTolerance( myTolerance );
+    myElementsOnShapePtr->SetAllNodes( true ); // "belong", while false means "lays on"
+    myElementsOnShapePtr->SetMesh( myMeshDS );
+    myElementsOnShapePtr->SetShape( myShape, myType );
+  }
+}
+
+bool BelongToGeom::IsSatisfy (long theId)
+{
+  if (myMeshDS == 0 || myShape.IsNull())
+    return false;
+
+  if (!myIsSubshape)
+  {
+    return myElementsOnShapePtr->IsSatisfy(theId);
+  }
+
+  // Case of sub-mesh
+
+  if (myType == SMDSAbs_Node)
+  {
+    if ( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ))
+    {
+      if ( aNode->getshapeId() < 1 )
+        return myElementsOnShapePtr->IsSatisfy(theId);
+      else
+        return mySubShapesIDs.Contains( aNode->getshapeId() );
+    }
+  }
+  else
+  {
+    if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
+    {
+      if ( anElem->GetType() == myType )
       {
-        isSatisfy = false;
+        if ( anElem->getshapeId() < 1 )
+          return myElementsOnShapePtr->IsSatisfy(theId);
+        else
+          return mySubShapesIDs.Contains( anElem->getshapeId() );
       }
     }
   }
 
-  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;
+  return false;
+}
+
+void BelongToGeom::SetType (SMDSAbs_ElementType theType)
+{
+  if ( myType != theType )
+  {
+    myType = theType;
+    init();
   }
+}
+
+SMDSAbs_ElementType BelongToGeom::GetType() const
+{
+  return myType;
+}
+
+TopoDS_Shape BelongToGeom::GetShape()
+{
+  return myShape;
+}
+
+const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const
+{
+  return myMeshDS;
+}
 
-  if (isSatisfy)
-    myIds.Add(theElemPtr->GetID());
+void BelongToGeom::SetTolerance (double theTolerance)
+{
+  myTolerance = theTolerance;
+  init();
 }
 
-TSequenceOfXYZ::TSequenceOfXYZ()
+double BelongToGeom::GetTolerance()
+{
+  return myTolerance;
+}
+
+/*
+  Class       : LyingOnGeom
+  Description : Predicate for verifying whether entiy lying or partially lying on
+  specified geometrical support
+*/
+
+LyingOnGeom::LyingOnGeom()
+  : myMeshDS(NULL),
+    myType(SMDSAbs_NbElementTypes),
+    myIsSubshape(false),
+    myTolerance(Precision::Confusion())
 {}
 
-TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n)
+Predicate* LyingOnGeom::clone() const
+{
+  LyingOnGeom* cln = new LyingOnGeom( *this );
+  cln->myElementsOnShapePtr.reset( static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ));
+  return cln;
+}
+
+void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
+{
+  if ( myMeshDS != theMesh )
+  {
+    myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
+    init();
+  }
+}
+
+void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
+{
+  if ( myShape != theShape )
+  {
+    myShape = theShape;
+    init();
+  }
+}
+
+void LyingOnGeom::init()
+{
+  if (!myMeshDS || myShape.IsNull()) return;
+
+  // is sub-shape of main shape?
+  TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
+  if (aMainShape.IsNull()) {
+    myIsSubshape = false;
+  }
+  else {
+    myIsSubshape = myMeshDS->IsGroupOfSubShapes( myShape );
+  }
+
+  if (myIsSubshape)
+  {
+    TopTools_IndexedMapOfShape shapes;
+    TopExp::MapShapes( myShape, shapes );
+    mySubShapesIDs.Clear();
+    for ( int i = 1; i <= shapes.Extent(); ++i )
+    {
+      int subID = myMeshDS->ShapeToIndex( shapes( i ));
+      if ( subID > 0 )
+        mySubShapesIDs.Add( subID );
+    }
+  }
+  // else // to be always ready to check an element not bound to geometry
+  {
+    if ( !myElementsOnShapePtr )
+      myElementsOnShapePtr.reset( new ElementsOnShape() );
+    myElementsOnShapePtr->SetTolerance( myTolerance );
+    myElementsOnShapePtr->SetAllNodes( false ); // lays on, while true means "belong"
+    myElementsOnShapePtr->SetMesh( myMeshDS );
+    myElementsOnShapePtr->SetShape( myShape, myType );
+  }
+}
+
+bool LyingOnGeom::IsSatisfy( long theId )
+{
+  if ( myMeshDS == 0 || myShape.IsNull() )
+    return false;
+
+  if (!myIsSubshape)
+  {
+    return myElementsOnShapePtr->IsSatisfy(theId);
+  }
+
+  // Case of sub-mesh
+
+  const SMDS_MeshElement* elem =
+    ( myType == SMDSAbs_Node ) ? myMeshDS->FindNode( theId ) : myMeshDS->FindElement( theId );
+
+  if ( mySubShapesIDs.Contains( elem->getshapeId() ))
+    return true;
+
+  if ( elem->GetType() != SMDSAbs_Node && elem->GetType() == myType )
+  {
+    SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator();
+    while ( nodeItr->more() )
+    {
+      const SMDS_MeshElement* aNode = nodeItr->next();
+      if ( mySubShapesIDs.Contains( aNode->getshapeId() ))
+        return true;
+    }
+  }
+
+  return false;
+}
+
+void LyingOnGeom::SetType( SMDSAbs_ElementType theType )
+{
+  if ( myType != theType )
+  {
+    myType = theType;
+    init();
+  }
+}
+
+SMDSAbs_ElementType LyingOnGeom::GetType() const
+{
+  return myType;
+}
+
+TopoDS_Shape LyingOnGeom::GetShape()
+{
+  return myShape;
+}
+
+const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const
+{
+  return myMeshDS;
+}
+
+void LyingOnGeom::SetTolerance (double theTolerance)
+{
+  myTolerance = theTolerance;
+  init();
+}
+
+double LyingOnGeom::GetTolerance()
+{
+  return myTolerance;
+}
+
+TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0)
+{}
+
+TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n), myElem(0)
 {}
 
-TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t)
+TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t), myElem(0)
 {}
 
-TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray)
+TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray), myElem(theSequenceOfXYZ.myElem)
 {}
 
 template <class InputIterator>
-TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd)
+TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd), myElem(0)
 {}
 
 TSequenceOfXYZ::~TSequenceOfXYZ()
@@ -3962,6 +5144,7 @@ TSequenceOfXYZ::~TSequenceOfXYZ()
 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
 {
   myArray = theSequenceOfXYZ.myArray;
+  myElem  = theSequenceOfXYZ.myElem;
   return *this;
 }
 
@@ -3995,6 +5178,11 @@ TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
   return myArray.size();
 }
 
+SMDSAbs_EntityType TSequenceOfXYZ::getElementEntity() const
+{
+  return myElem ? myElem->GetEntityType() : SMDSEntity_Last;
+}
+
 TMeshModifTracer::TMeshModifTracer():
   myMeshModifTime(0), myMesh(0)
 {