Salome HOME
[bos #23982] EDF 22984 - aspect ratio of hexa
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
index 0623399e9dba8246b5958f7f32117d1937ec5324..1ef3a5f2f4018efa8b7c14ca2454928338ecaa43 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2021  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
 #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_MeshNode.hxx"
-#include "SMDS_QuadraticEdge.hxx"
-#include "SMDS_QuadraticFaceOfNodes.hxx"
 #include "SMDS_VolumeTool.hxx"
 #include "SMESHDS_GroupBase.hxx"
 #include "SMESHDS_GroupOnFilter.hxx"
 #include "SMESH_MeshAlgos.hxx"
 #include "SMESH_OctreeNode.hxx"
 
+#include <GEOMUtils.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>
@@ -93,6 +98,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 );
@@ -129,13 +143,13 @@ namespace {
     //  Case 1          Case 2
     //  |     |      |        |      |
     //  |     |      |        |      |
-    //  +-----+------+  +-----+------+ 
+    //  +-----+------+  +-----+------+
     //  |            |  |            |
     //  |            |  |            |
-    // 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
+    // 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 );
@@ -175,7 +189,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;
 
@@ -220,7 +234,7 @@ bool NumericalFunctor::GetPoints(const int       theId,
     return false;
 
   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
-  if ( !anElem || anElem->GetType() != this->GetType() )
+  if ( !IsApplicable( anElem ))
     return false;
 
   return GetPoints( anElem, theRes );
@@ -238,34 +252,12 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
   theRes.setElement( anElem );
 
   // Get nodes of the element
-  SMDS_ElemIteratorPtr anIter;
-
-  if ( anElem->IsQuadratic() ) {
-    switch ( anElem->GetType() ) {
-    case SMDSAbs_Edge:
-      anIter = dynamic_cast<const SMDS_VtkEdge*>
-        (anElem)->interlacedNodesElemIterator();
-      break;
-    case SMDSAbs_Face:
-      anIter = dynamic_cast<const SMDS_VtkFace*>
-        (anElem)->interlacedNodesElemIterator();
-      break;
-    default:
-      anIter = anElem->nodesIterator();
-    }
-  }
-  else {
-    anIter = anElem->nodesIterator();
-  }
-
+  SMDS_NodeIteratorPtr anIter= anElem->interlacedNodesIterator();
   if ( anIter ) {
-    double xyz[3];
+    SMESH_NodeXYZ p;
     while( anIter->more() ) {
-      if ( const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( anIter->next() ))
-      {
-        aNode->GetXYZ( xyz );
-        theRes.push_back( gp_XYZ( xyz[0], xyz[1], xyz[2] ));
-      }
+      if ( p.Set( anIter->next() ))
+        theRes.push_back( p );
     }
   }
 
@@ -301,6 +293,24 @@ double NumericalFunctor::Round( const double & aVal )
   return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal;
 }
 
+//================================================================================
+/*!
+ * \brief Return true if a value can be computed for a given element.
+ *        Some NumericalFunctor's are meaningful for elements of a certain
+ *        geometry only.
+ */
+//================================================================================
+
+bool NumericalFunctor::IsApplicable( const SMDS_MeshElement* element ) const
+{
+  return element && element->GetType() == this->GetType();
+}
+
+bool NumericalFunctor::IsApplicable( long theElementId ) const
+{
+  return IsApplicable( myMesh->FindElement( theElementId ));
+}
+
 //================================================================================
 /*!
  * \brief Return histogram of functor values
@@ -312,12 +322,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,
-                                    const bool           isLogarithmic)
+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 ||
@@ -336,7 +346,7 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
   }
   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 ));
   }
@@ -482,7 +492,7 @@ double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P )
   //     }
   // }
   // { // polygons
-    
+
   // }
 
   if( myPrecision >= 0 )
@@ -614,7 +624,8 @@ double MaxElementLength3D::GetValue( long theElementId )
       aVal = Max(aVal,Max(L7,L8));
       break;
     }
-    case SMDSEntity_Quad_Penta: { // quadratic pentas
+    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 ));
@@ -710,21 +721,25 @@ SMDSAbs_ElementType MaxElementLength3D::GetType() const
 
 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;
+
+  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 = getAngle( P( i-1 ), P( i ), P( i+1 ) );
-    aMin = Min(aMin,A0);
+    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
@@ -751,19 +766,9 @@ double AspectRatio::GetValue( long theId )
 {
   double aVal = 0;
   myCurrElement = myMesh->FindElement( theId );
-  if ( myCurrElement && myCurrElement->GetVtkType() == VTK_QUAD )
-  {
-    // issue 21723
-    vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
-    if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
-      aVal = Round( vtkMeshQuality::QuadAspectRatio( avtkCell ));
-  }
-  else
-  {
-    TSequenceOfXYZ P;
-    if ( GetPoints( myCurrElement, P ))
-      aVal = Round( GetValue( P ));
-  }
+  TSequenceOfXYZ P;
+  if ( GetPoints( myCurrElement, P ))
+    aVal = Round( GetValue( P ));
   return aVal;
 }
 
@@ -783,58 +788,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 ) );
+    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) );
+    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) );
@@ -843,42 +841,42 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
     //
     // alpha = sqrt( 1/32 )
     // L = max( L1, L2, L3, L4, D1, D2 )
-    // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+    // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
     // C2 = min( S1, S2, S3, S4 )
     // Li - lengths of the edges
     // Di - lengths of the diagonals
     // Si - areas of the triangles
     const double alpha = sqrt( 1 / 32. );
     double L = Max( aLen[ 0 ],
-                 Max( aLen[ 1 ],
-                   Max( aLen[ 2 ],
-                     Max( aLen[ 3 ],
-                       Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
-    double C1 = sqrt( aLen[0] * aLen[0] +
-                        aLen[1] * aLen[1] +
-                        aLen[2] * aLen[2] +
-                        aLen[3] * aLen[3] ) / 4. );
+                    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] );
     double C2 = Min( anArea[ 0 ],
-                  Min( anArea[ 1 ],
-                    Min( anArea[ 2 ], anArea[ 3 ] ) ) );
+                     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) );
@@ -887,24 +885,24 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
     //
     // alpha = sqrt( 1/32 )
     // L = max( L1, L2, L3, L4, D1, D2 )
-    // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+    // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
     // C2 = min( S1, S2, S3, S4 )
     // Li - lengths of the edges
     // Di - lengths of the diagonals
     // Si - areas of the triangles
     const double alpha = sqrt( 1 / 32. );
     double L = Max( aLen[ 0 ],
-                 Max( aLen[ 1 ],
-                   Max( aLen[ 2 ],
-                     Max( aLen[ 3 ],
-                       Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
-    double C1 = sqrt( aLen[0] * aLen[0] +
-                        aLen[1] * aLen[1] +
-                        aLen[2] * aLen[2] +
-                        aLen[3] * aLen[3] ) / 4. );
+                    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] );
     double C2 = Min( anArea[ 0 ],
-                  Min( anArea[ 1 ],
-                    Min( anArea[ 2 ], anArea[ 3 ] ) ) );
+                     Min( anArea[ 1 ],
+                          Min( anArea[ 2 ], anArea[ 3 ] ) ) );
     if ( C2 <= theEps )
       return theInf;
     return alpha * L * C1 / C2;
@@ -912,6 +910,11 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
   return 0;
 }
 
+bool AspectRatio::IsApplicable( const SMDS_MeshElement* element ) const
+{
+  return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() );
+}
+
 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
 {
   // the aspect ratio is in the range [1.0,infinity]
@@ -994,6 +997,102 @@ namespace{
     return aHeight;
   }
 
+  //================================================================================
+  /*!
+   * \brief Standard quality of a tetrahedron but not normalized
+   */
+  //================================================================================
+
+  double tetQualityByHomardMethod( const gp_XYZ & p1,
+                                   const gp_XYZ & p2,
+                                   const gp_XYZ & p3,
+                                   const gp_XYZ & p4 )
+  {
+    gp_XYZ edgeVec[6];
+    edgeVec[0] = ( p1 - p2 );
+    edgeVec[1] = ( p2 - p3 );
+    edgeVec[2] = ( p3 - p1 );
+    edgeVec[3] = ( p4 - p1 );
+    edgeVec[4] = ( p4 - p2 );
+    edgeVec[5] = ( p4 - p3 );
+
+    double maxEdgeLen2            = edgeVec[0].SquareModulus();
+    maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[1].SquareModulus() );
+    maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[2].SquareModulus() );
+    maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[3].SquareModulus() );
+    maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[4].SquareModulus() );
+    maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[5].SquareModulus() );
+    double maxEdgeLen = Sqrt( maxEdgeLen2 );
+
+    gp_XYZ cross01 = edgeVec[0] ^ edgeVec[1];
+    double sumArea = ( cross01                 ).Modulus(); // actually double area
+    sumArea       += ( edgeVec[0] ^ edgeVec[3] ).Modulus();
+    sumArea       += ( edgeVec[1] ^ edgeVec[4] ).Modulus();
+    sumArea       += ( edgeVec[2] ^ edgeVec[5] ).Modulus();
+
+    double sixVolume = Abs( cross01 * edgeVec[4] ); // 6 * volume
+    double quality   = maxEdgeLen * sumArea / sixVolume; // not normalized!!!
+    return quality;
+  }
+
+  //================================================================================
+  /*!
+   * \brief HOMARD method of hexahedron quality
+   * 1. Decompose the hexa into 24 tetra: each face is splitted into 4 triangles by
+   *    adding the diagonals and every triangle is connected to the center of the hexa.
+   * 2. Compute the quality of every tetra with the same formula as for the standard quality,
+   *    except that the factor for the normalization is not the same because the final goal
+   *    is to have a quality equal to 1 for a perfect cube. So the formula is:
+   *    qual = max(lengthes of 6 edges) * (sum of surfaces of 4 faces) / (7.6569*6*volume)
+   * 3. The quality of the hexa is the highest value of the qualities of the 24 tetra
+   */
+  //================================================================================
+
+  double hexQualityByHomardMethod( const TSequenceOfXYZ& P )
+  {
+    gp_XYZ quadCenter[6];
+    quadCenter[0] = ( P(1) + P(2) + P(3) + P(4) ) / 4.;
+    quadCenter[1] = ( P(5) + P(6) + P(7) + P(8) ) / 4.;
+    quadCenter[2] = ( P(1) + P(2) + P(6) + P(5) ) / 4.;
+    quadCenter[3] = ( P(2) + P(3) + P(7) + P(6) ) / 4.;
+    quadCenter[4] = ( P(3) + P(4) + P(8) + P(7) ) / 4.;
+    quadCenter[5] = ( P(1) + P(4) + P(8) + P(5) ) / 4.;
+
+    gp_XYZ hexCenter = ( P(1) + P(2) + P(3) + P(4) + P(5) + P(6) + P(7) + P(8) ) / 8.;
+
+    // quad 1 ( 1 2 3 4 )
+    double quality =        tetQualityByHomardMethod( P(1), P(2), quadCenter[0], hexCenter );
+    quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[0], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[0], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(4), P(1), quadCenter[0], hexCenter ));
+    // quad 2 ( 5 6 7 8 )
+    quality = Max( quality, tetQualityByHomardMethod( P(5), P(6), quadCenter[1], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(6), P(7), quadCenter[1], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(7), P(8), quadCenter[1], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[1], hexCenter ));
+    // quad 3 ( 1 2 6 5 )
+    quality = Max( quality, tetQualityByHomardMethod( P(1), P(2), quadCenter[2], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(2), P(6), quadCenter[2], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(6), P(5), quadCenter[2], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[2], hexCenter ));
+    // quad 4 ( 2 3 7 6 )
+    quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[3], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(3), P(7), quadCenter[3], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(7), P(6), quadCenter[3], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(6), P(2), quadCenter[3], hexCenter ));
+    // quad 5 ( 3 4 8 7 )
+    quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[4], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[4], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(8), P(7), quadCenter[4], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(7), P(3), quadCenter[4], hexCenter ));
+    // quad 6 ( 1 4 8 5 )
+    quality = Max( quality, tetQualityByHomardMethod( P(1), P(4), quadCenter[5], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[5], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[5], hexCenter ));
+    quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[5], hexCenter ));
+
+    return quality / 7.65685424949;
+  }
 }
 
 double AspectRatio3D::GetValue( long theId )
@@ -1005,8 +1104,8 @@ double AspectRatio3D::GetValue( long theId )
     // Action from CoTech | ACTION 31.3:
     // EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with
     // those of ParaView. The library used by ParaView for those calculations can be reused in SMESH.
-    vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
-    if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
+    vtkUnstructuredGrid* grid = const_cast<SMDS_Mesh*>( myMesh )->GetGrid();
+    if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->GetVtkID() ))
       aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
   }
   else
@@ -1018,6 +1117,11 @@ double AspectRatio3D::GetValue( long theId )
   return aVal;
 }
 
+bool AspectRatio3D::IsApplicable( const SMDS_MeshElement* element ) const
+{
+  return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() );
+}
+
 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
 {
   double aQuality = 0.0;
@@ -1025,12 +1129,12 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
 
   int nbNodes = P.size();
 
-  if(myCurrElement->IsQuadratic()) {
-    if(nbNodes==10) nbNodes=4; // quadratic tetrahedron
+  if( myCurrElement->IsQuadratic() ) {
+    if     (nbNodes==10) nbNodes=4; // quadratic tetrahedron
     else if(nbNodes==13) nbNodes=5; // quadratic pyramid
     else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
     else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
-    else if(nbNodes==27) nbNodes=8; // quadratic hexahedron
+    else if(nbNodes==27) nbNodes=8; // tri-quadratic hexahedron
     else return aQuality;
   }
 
@@ -1074,7 +1178,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
   case 5:{
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
-      aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+      aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
     }
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
@@ -1093,7 +1197,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
   case 6:{
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
-      aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+      aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
     }
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
@@ -1118,9 +1222,13 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
     break;
   }
   case 8:{
+
+    return hexQualityByHomardMethod( P ); // bos #23982
+
+
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
-      aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+      aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
     }
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
@@ -1255,7 +1363,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
   case 12:
     {
       gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
-      aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
+      aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8]));
     }
     {
       gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
@@ -1269,7 +1377,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
   } // switch(nbNodes)
 
   if ( nbNodes > 4 ) {
-    // avaluate aspect ratio of quadranle faces
+    // evaluate aspect ratio of quadrangle faces
     AspectRatio aspect2D;
     SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
     int nbFaces = SMDS_VolumeTool::NbFaces( type );
@@ -1278,7 +1386,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
       if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
         continue;
       const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
-      for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face
+      for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadrangle face
         points( p + 1 ) = P( pInd[ p ] + 1 );
       aQuality = std::max( aQuality, aspect2D.GetValue( points ));
     }
@@ -1307,6 +1415,11 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const
 */
 //================================================================================
 
+bool Warping::IsApplicable( const SMDS_MeshElement* element ) const
+{
+  return NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4;
+}
+
 double Warping::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
@@ -1371,6 +1484,11 @@ SMDSAbs_ElementType Warping::GetType() const
 */
 //================================================================================
 
+bool Taper::IsApplicable( const SMDS_MeshElement* element ) const
+{
+  return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4 );
+}
+
 double Taper::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
@@ -1429,6 +1547,11 @@ static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ
   return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
 }
 
+bool Skew::IsApplicable( const SMDS_MeshElement* element ) const
+{
+  return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() <= 4 );
+}
+
 double Skew::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 3 && P.size() != 4 )
@@ -1545,444 +1668,555 @@ SMDSAbs_ElementType Length::GetType() const
 
 //================================================================================
 /*
-  Class       : Length2D
-  Description : Functor for calculating minimal length of edge
+  Class       : Length3D
+  Description : Functor for calculating minimal length of element edge
 */
 //================================================================================
 
-double Length2D::GetValue( long theElementId )
-{
-  TSequenceOfXYZ P;
-
-  if ( GetPoints( theElementId, P ))
-  {
-    double aVal = 0;
-    int len = P.size();
-    SMDSAbs_EntityType aType = P.getElementEntity();
-
-    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));
-      }
-      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));
-      }
-      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){ // 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
-        aVal = Min(aVal,Min(L7,L8));
-      }
-      break;
-    case SMDSEntity_Penta:
-      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 = 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 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
-      }
-      break;
-    case SMDSEntity_Quad_Pyramid:
-      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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
-        aVal = Min(aVal,Min(L7,L8));
-      }
-      break;
-    case SMDSEntity_Quad_Penta:
-      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 = 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 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 = 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;
-    }
-
-    if (aVal < 0 ) {
-      return 0.;
-    }
-
-    if ( myPrecision >= 0 )
-    {
-      double prec = pow( 10., (double)( myPrecision ) );
-      aVal = floor( aVal * prec + 0.5 ) / prec;
-    }
-
-    return aVal;
-
-  }
-  return 0.;
-}
-
-double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
-{
-  // meaningless as it is not a quality control functor
-  return Value;
-}
-
-SMDSAbs_ElementType Length2D::GetType() const
-{
-  return SMDSAbs_Face;
-}
-
-Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
-  myLength(theLength)
-{
-  myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
-  if(thePntId1 > thePntId2){
-    myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
-  }
-}
-
-bool Length2D::Value::operator<(const Length2D::Value& x) const
-{
-  if(myPntId[0] < x.myPntId[0]) return true;
-  if(myPntId[0] == x.myPntId[0])
-    if(myPntId[1] < x.myPntId[1]) return true;
-  return false;
-}
-
-void Length2D::GetValues(TValues& theValues)
+Length3D::Length3D():
+  Length2D ( SMDSAbs_Volume )
 {
-  TValues aValues;
-  SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
-  for(; anIter->more(); ){
-    const SMDS_MeshFace* anElem = anIter->next();
-
-    if(anElem->IsQuadratic()) {
-      const SMDS_VtkFace* F =
-        dynamic_cast<const SMDS_VtkFace*>(anElem);
-      // use special nodes iterator
-      SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
-      long aNodeId[4];
-      gp_Pnt P[4];
-
-      double aLength;
-      const SMDS_MeshElement* aNode;
-      if(anIter->more()){
-        aNode = anIter->next();
-        const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
-        P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
-        aNodeId[0] = aNodeId[1] = aNode->GetID();
-        aLength = 0;
-      }
-      for(; anIter->more(); ){
-        const SMDS_MeshNode* N1 = static_cast<const SMDS_MeshNode*> (anIter->next());
-        P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z());
-        aNodeId[2] = N1->GetID();
-        aLength = P[1].Distance(P[2]);
-        if(!anIter->more()) break;
-        const SMDS_MeshNode* N2 = static_cast<const SMDS_MeshNode*> (anIter->next());
-        P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z());
-        aNodeId[3] = N2->GetID();
-        aLength += P[2].Distance(P[3]);
-        Value aValue1(aLength,aNodeId[1],aNodeId[2]);
-        Value aValue2(aLength,aNodeId[2],aNodeId[3]);
-        P[1] = P[3];
-        aNodeId[1] = aNodeId[3];
-        theValues.insert(aValue1);
-        theValues.insert(aValue2);
-      }
-      aLength += P[2].Distance(P[0]);
-      Value aValue1(aLength,aNodeId[1],aNodeId[2]);
-      Value aValue2(aLength,aNodeId[2],aNodeId[0]);
-      theValues.insert(aValue1);
-      theValues.insert(aValue2);
-    }
-    else {
-      SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
-      long aNodeId[2];
-      gp_Pnt P[3];
-
-      double aLength;
-      const SMDS_MeshElement* aNode;
-      if(aNodesIter->more()){
-        aNode = aNodesIter->next();
-        const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
-        P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
-        aNodeId[0] = aNodeId[1] = aNode->GetID();
-        aLength = 0;
-      }
-      for(; aNodesIter->more(); ){
-        aNode = aNodesIter->next();
-        const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
-        long anId = aNode->GetID();
-        
-        P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
-        
-        aLength = P[1].Distance(P[2]);
-        
-        Value aValue(aLength,aNodeId[1],anId);
-        aNodeId[1] = anId;
-        P[1] = P[2];
-        theValues.insert(aValue);
-      }
-
-      aLength = P[0].Distance(P[1]);
-
-      Value aValue(aLength,aNodeId[0],aNodeId[1]);
-      theValues.insert(aValue);
-    }
-  }
 }
 
 //================================================================================
 /*
-  Class       : MultiConnection
-  Description : Functor for calculating number of faces conneted to the edge
+  Class       : Length2D
+  Description : Functor for calculating minimal length of element edge
 */
 //================================================================================
 
-double MultiConnection::GetValue( const TSequenceOfXYZ& P )
-{
-  return 0;
-}
-double MultiConnection::GetValue( long theId )
-{
-  return getNbMultiConnection( myMesh, theId );
-}
-
-double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
-{
-  // meaningless as it is not quality control functor
-  return Value;
-}
-
-SMDSAbs_ElementType MultiConnection::GetType() const
+Length2D::Length2D( SMDSAbs_ElementType type ):
+  myType ( type )
 {
-  return SMDSAbs_Edge;
 }
 
-//================================================================================
-/*
-  Class       : MultiConnection2D
-  Description : Functor for calculating number of faces conneted to the edge
-*/
-//================================================================================
-
-double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
+bool Length2D::IsApplicable( const SMDS_MeshElement* element ) const
 {
-  return 0;
+  return ( NumericalFunctor::IsApplicable( element ) &&
+           element->GetEntityType() != SMDSEntity_Polyhedra );
 }
 
-double MultiConnection2D::GetValue( long theElementId )
+double Length2D::GetValue( const TSequenceOfXYZ& P )
 {
-  int aResult = 0;
-
-  const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
-  SMDSAbs_ElementType aType = aFaceElem->GetType();
+  double aVal = 0;
+  int len = P.size();
+  SMDSAbs_EntityType aType = P.getElementEntity();
 
   switch (aType) {
-  case SMDSAbs_Face:
-    {
-      int i = 0, len = aFaceElem->NbNodes();
-      SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
-      if (!anIter) break;
-
-      const SMDS_MeshNode *aNode, *aNode0;
-      TColStd_MapOfInteger aMap, aMapPrev;
-
-      for (i = 0; i <= len; i++) {
-        aMapPrev = aMap;
-        aMap.Clear();
-
-        int aNb = 0;
-        if (anIter->more()) {
-          aNode = (SMDS_MeshNode*)anIter->next();
-        } else {
-          if (i == len)
-            aNode = aNode0;
-          else
-            break;
-        }
-        if (!aNode) break;
-        if (i == 0) aNode0 = aNode;
-
-        SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
-        while (anElemIter->more()) {
-          const SMDS_MeshElement* anElem = anElemIter->next();
-          if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
-            int anId = anElem->GetID();
-
-            aMap.Add(anId);
+  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));
+    }
+    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));
+    }
+    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;
+  }
+
+  if (aVal < 0 ) {
+    return 0.;
+  }
+
+  if ( myPrecision >= 0 )
+  {
+    double prec = pow( 10., (double)( myPrecision ) );
+    aVal = floor( aVal * prec + 0.5 ) / prec;
+  }
+
+  return aVal;
+}
+
+double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+  // meaningless as it is not a quality control functor
+  return Value;
+}
+
+SMDSAbs_ElementType Length2D::GetType() const
+{
+  return myType;
+}
+
+Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
+  myLength(theLength)
+{
+  myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
+  if(thePntId1 > thePntId2){
+    myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
+  }
+}
+
+bool Length2D::Value::operator<(const Length2D::Value& x) const
+{
+  if(myPntId[0] < x.myPntId[0]) return true;
+  if(myPntId[0] == x.myPntId[0])
+    if(myPntId[1] < x.myPntId[1]) return true;
+  return false;
+}
+
+void Length2D::GetValues(TValues& theValues)
+{
+  if ( myType == SMDSAbs_Face )
+  {
+    for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
+    {
+      const SMDS_MeshFace* anElem = anIter->next();
+      if ( anElem->IsQuadratic() )
+      {
+        // use special nodes iterator
+        SMDS_NodeIteratorPtr anIter = anElem->interlacedNodesIterator();
+        long aNodeId[4] = { 0,0,0,0 };
+        gp_Pnt P[4];
+
+        double aLength = 0;
+        if ( anIter->more() )
+        {
+          const SMDS_MeshNode* aNode = anIter->next();
+          P[0] = P[1] = SMESH_NodeXYZ( aNode );
+          aNodeId[0] = aNodeId[1] = aNode->GetID();
+          aLength = 0;
+        }
+        for ( ; anIter->more(); )
+        {
+          const SMDS_MeshNode* N1 = anIter->next();
+          P[2] = SMESH_NodeXYZ( N1 );
+          aNodeId[2] = N1->GetID();
+          aLength = P[1].Distance(P[2]);
+          if(!anIter->more()) break;
+          const SMDS_MeshNode* N2 = anIter->next();
+          P[3] = SMESH_NodeXYZ( N2 );
+          aNodeId[3] = N2->GetID();
+          aLength += P[2].Distance(P[3]);
+          Value aValue1(aLength,aNodeId[1],aNodeId[2]);
+          Value aValue2(aLength,aNodeId[2],aNodeId[3]);
+          P[1] = P[3];
+          aNodeId[1] = aNodeId[3];
+          theValues.insert(aValue1);
+          theValues.insert(aValue2);
+        }
+        aLength += P[2].Distance(P[0]);
+        Value aValue1(aLength,aNodeId[1],aNodeId[2]);
+        Value aValue2(aLength,aNodeId[2],aNodeId[0]);
+        theValues.insert(aValue1);
+        theValues.insert(aValue2);
+      }
+      else {
+        SMDS_NodeIteratorPtr aNodesIter = anElem->nodeIterator();
+        long aNodeId[2] = {0,0};
+        gp_Pnt P[3];
+
+        double aLength;
+        const SMDS_MeshElement* aNode;
+        if ( aNodesIter->more())
+        {
+          aNode = aNodesIter->next();
+          P[0] = P[1] = SMESH_NodeXYZ( aNode );
+          aNodeId[0] = aNodeId[1] = aNode->GetID();
+          aLength = 0;
+        }
+        for( ; aNodesIter->more(); )
+        {
+          aNode = aNodesIter->next();
+          long anId = aNode->GetID();
+
+          P[2] = SMESH_NodeXYZ( aNode );
+
+          aLength = P[1].Distance(P[2]);
+
+          Value aValue(aLength,aNodeId[1],anId);
+          aNodeId[1] = anId;
+          P[1] = P[2];
+          theValues.insert(aValue);
+        }
+
+        aLength = P[0].Distance(P[1]);
+
+        Value aValue(aLength,aNodeId[0],aNodeId[1]);
+        theValues.insert(aValue);
+      }
+    }
+  }
+  else
+  {
+    // not implemented
+  }
+}
+
+//================================================================================
+/*
+  Class       : Deflection2D
+  Description : computes distance between a face center and an underlying surface
+*/
+//================================================================================
+
+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 ( SMDS_FacePositionPtr fPos = 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;
+}
+double MultiConnection::GetValue( long theId )
+{
+  return getNbMultiConnection( myMesh, theId );
+}
+
+double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+  // meaningless as it is not quality control functor
+  return Value;
+}
+
+SMDSAbs_ElementType MultiConnection::GetType() const
+{
+  return SMDSAbs_Edge;
+}
+
+//================================================================================
+/*
+  Class       : MultiConnection2D
+  Description : Functor for calculating number of faces conneted to the edge
+*/
+//================================================================================
+
+double MultiConnection2D::GetValue( const TSequenceOfXYZ& /*P*/ )
+{
+  return 0;
+}
+
+double MultiConnection2D::GetValue( long theElementId )
+{
+  int aResult = 0;
+
+  const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
+  SMDSAbs_ElementType aType = aFaceElem->GetType();
+
+  switch (aType) {
+  case SMDSAbs_Face:
+    {
+      int i = 0, len = aFaceElem->NbNodes();
+      SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
+      if (!anIter) break;
+
+      const SMDS_MeshNode *aNode, *aNode0 = 0;
+      TColStd_MapOfInteger aMap, aMapPrev;
+
+      for (i = 0; i <= len; i++) {
+        aMapPrev = aMap;
+        aMap.Clear();
+
+        int aNb = 0;
+        if (anIter->more()) {
+          aNode = (SMDS_MeshNode*)anIter->next();
+        } else {
+          if (i == len)
+            aNode = aNode0;
+          else
+            break;
+        }
+        if (!aNode) break;
+        if (i == 0) aNode0 = aNode;
+
+        SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
+        while (anElemIter->more()) {
+          const SMDS_MeshElement* anElem = anElemIter->next();
+          if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
+            int anId = anElem->GetID();
+
+            aMap.Add(anId);
             if (aMapPrev.Contains(anId)) {
               aNb++;
             }
@@ -2029,59 +2263,24 @@ bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) cons
 void MultiConnection2D::GetValues(MValues& theValues)
 {
   if ( !myMesh ) return;
-  SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
-  for(; anIter->more(); ){
-    const SMDS_MeshFace* anElem = anIter->next();
-    SMDS_ElemIteratorPtr aNodesIter;
-    if ( anElem->IsQuadratic() )
-      aNodesIter = dynamic_cast<const SMDS_VtkFace*>
-        (anElem)->interlacedNodesElemIterator();
-    else
-      aNodesIter = anElem->nodesIterator();
-    long aNodeId[3];
+  for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
+  {
+    const SMDS_MeshFace*     anElem = anIter->next();
+    SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
 
-    //int aNbConnects=0;
-    const SMDS_MeshNode* aNode0;
-    const SMDS_MeshNode* aNode1;
+    const SMDS_MeshNode* aNode1 = anElem->GetNode( anElem->NbNodes() - 1 );
     const SMDS_MeshNode* aNode2;
-    if(aNodesIter->more()){
-      aNode0 = (SMDS_MeshNode*) aNodesIter->next();
-      aNode1 = aNode0;
-      const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
-      aNodeId[0] = aNodeId[1] = aNodes->GetID();
-    }
-    for(; aNodesIter->more(); ) {
-      aNode2 = (SMDS_MeshNode*) aNodesIter->next();
-      long anId = aNode2->GetID();
-      aNodeId[2] = anId;
-
-      Value aValue(aNodeId[1],aNodeId[2]);
-      MValues::iterator aItr = theValues.find(aValue);
-      if (aItr != theValues.end()){
-        aItr->second += 1;
-        //aNbConnects = nb;
-      }
-      else {
-        theValues[aValue] = 1;
-        //aNbConnects = 1;
-      }
-      //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
-      aNodeId[1] = aNodeId[2];
+    for ( ; aNodesIter->more(); )
+    {
+      aNode2 = aNodesIter->next();
+
+      Value aValue ( aNode1->GetID(), aNode2->GetID() );
+      MValues::iterator aItr = theValues.insert( std::make_pair( aValue, 0 )).first;
+      aItr->second++;
       aNode1 = aNode2;
     }
-    Value aValue(aNodeId[0],aNodeId[2]);
-    MValues::iterator aItr = theValues.find(aValue);
-    if (aItr != theValues.end()) {
-      aItr->second += 1;
-      //aNbConnects = nb;
-    }
-    else {
-      theValues[aValue] = 1;
-      //aNbConnects = 1;
-    }
-    //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
   }
-
+  return;
 }
 
 //================================================================================
@@ -2096,7 +2295,7 @@ double BallDiameter::GetValue( long theId )
   double diameter = 0;
 
   if ( const SMDS_BallElement* ball =
-       dynamic_cast<const SMDS_BallElement*>( myMesh->FindElement( theId )))
+       myMesh->DownCast< SMDS_BallElement >( myMesh->FindElement( theId )))
   {
     diameter = ball->GetDiameter();
   }
@@ -2114,6 +2313,42 @@ 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
@@ -2142,7 +2377,19 @@ bool BadOrientedVolume::IsSatisfy( long theId )
     return false;
 
   SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
-  return !vTool.IsForward();
+
+  bool isOk = true;
+  if ( vTool.IsPoly() )
+  {
+    isOk = true;
+    for ( int i = 0; i < vTool.NbFaces() && isOk; ++i )
+      isOk = vTool.IsFaceExternal( i );
+  }
+  else
+  {
+    isOk = vTool.IsForward();
+  }
+  return !isOk;
 }
 
 SMDSAbs_ElementType BadOrientedVolume::GetType() const
@@ -2163,7 +2410,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;
       }
@@ -2292,25 +2539,34 @@ SMDSAbs_ElementType CoincidentNodes::GetType() const
   return SMDSAbs_Node;
 }
 
+void CoincidentNodes::SetTolerance( const double theToler )
+{
+  if ( myToler != theToler )
+  {
+    SetMesh(0);
+    myToler = theToler;
+  }
+}
+
 void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
 {
   myMeshModifTracer.SetMesh( theMesh );
   if ( myMeshModifTracer.IsMeshModified() )
   {
     TIDSortedNodeSet nodesToCheck;
-    SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
+    SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
     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() );
     }
@@ -2342,7 +2598,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() )
@@ -2421,18 +2677,14 @@ void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
 
 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId  )
 {
-  TColStd_MapOfInteger aMap;
-  for ( int i = 0; i < 2; i++ )
+  SMDS_ElemIteratorPtr anElemIter = theNodes[ 0 ]->GetInverseElementIterator(SMDSAbs_Face);
+  while( anElemIter->more() )
   {
-    SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face);
-    while( anElemIter->more() )
+    if ( const SMDS_MeshElement* anElem = anElemIter->next())
     {
-      if ( const SMDS_MeshElement* anElem = anElemIter->next())
-      {
-        const int anId = anElem->GetID();
-        if ( anId != theFaceId && !aMap.Add( anId ))
-          return false;
-      }
+      const int anId = anElem->GetID();
+      if ( anId != theFaceId && anElem->GetNodeIndex( theNodes[1] ) >= 0 )
+        return false;
     }
   }
   return true;
@@ -2501,31 +2753,21 @@ inline void UpdateBorders(const FreeEdges::Border& theBorder,
 void FreeEdges::GetBoreders(TBorders& theBorders)
 {
   TBorders aRegistry;
-  SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
-  for(; anIter->more(); ){
+  for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
+  {
     const SMDS_MeshFace* anElem = anIter->next();
     long anElemId = anElem->GetID();
-    SMDS_ElemIteratorPtr aNodesIter;
-    if ( anElem->IsQuadratic() )
-      aNodesIter = static_cast<const SMDS_VtkFace*>(anElem)->
-        interlacedNodesElemIterator();
-    else
-      aNodesIter = anElem->nodesIterator();
-    long aNodeId[2];
-    const SMDS_MeshElement* aNode;
-    if(aNodesIter->more()){
-      aNode = aNodesIter->next();
-      aNodeId[0] = aNodeId[1] = aNode->GetID();
-    }
-    for(; aNodesIter->more(); ){
-      aNode = aNodesIter->next();
-      long anId = aNode->GetID();
-      Border aBorder(anElemId,aNodeId[1],anId);
-      aNodeId[1] = anId;
-      UpdateBorders(aBorder,aRegistry,theBorders);
+    SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
+    if ( !aNodesIter->more() ) continue;
+    long aNodeId[2] = {0,0};
+    aNodeId[0] = anElem->GetNode( anElem->NbNodes()-1 )->GetID();
+    for ( ; aNodesIter->more(); )
+    {
+      aNodeId[1] = aNodesIter->next()->GetID();
+      Border aBorder( anElemId, aNodeId[0], aNodeId[1] );
+      UpdateBorders( aBorder, aRegistry, theBorders );
+      aNodeId[0] = aNodeId[1];
     }
-    Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
-    UpdateBorders(aBorder,aRegistry,theBorders);
   }
 }
 
@@ -2589,28 +2831,30 @@ bool FreeFaces::IsSatisfy( long theId )
   int nbNode = aFace->NbNodes();
 
   // collect volumes to check that number of volumes 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
+  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++;
-    } 
+    }
   }
   int nbVol = 0;
   TItrMapOfVolume volItr = mapOfVol.begin();
   TItrMapOfVolume volEnd = mapOfVol.end();
   for ( ; volItr != volEnd; ++volItr )
     if ( (*volItr).second >= nbNode )
-       nbVol++;
-  // face is not free if number of volumes constructed on thier nodes more than one
+      nbVol++;
+  // face is not free if number of volumes constructed on their nodes more than one
   return (nbVol < 2);
 }
 
@@ -2658,7 +2902,7 @@ SMDSAbs_ElementType LinearOrQuadratic::GetType() const
 //================================================================================
 /*
   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
 */
 //================================================================================
 
@@ -2704,8 +2948,8 @@ void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
     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();
+  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);
@@ -2932,10 +3176,10 @@ void ConnectedElements::SetPoint( double x, double y, double z )
   // find myNodeID by myXYZ if possible
   if ( myMeshModifTracer.GetMesh() )
   {
-    auto_ptr<SMESH_ElementSearcher> searcher
+    SMESHUtils::Deleter<SMESH_ElementSearcher> searcher
       ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
 
-    vector< const SMDS_MeshElement* > foundElems;
+    std::vector< const SMDS_MeshElement* > foundElems;
     searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
 
     if ( !foundElems.empty() )
@@ -2961,7 +3205,7 @@ bool ConnectedElements::IsSatisfy( long theElementId )
     if ( !node0 )
       return false;
 
-    list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
+    std::list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
     std::set< int > checkedNodeIDs;
     // algo:
     // foreach node in nodeQueue:
@@ -2979,7 +3223,7 @@ bool ConnectedElements::IsSatisfy( long theElementId )
       {
         // keep elements of myType
         const SMDS_MeshElement* element = eIt->next();
-        if ( element->GetType() == myType )
+        if ( myType == SMDSAbs_All || element->GetType() == myType )
           myOkIDs.insert( myOkIDs.end(), element->GetID() );
 
         // enqueue nodes of the element
@@ -3018,7 +3262,7 @@ namespace
     double dot = v1 * v2; // cos * |v1| * |v2|
     double l1  = v1.SquareMagnitude();
     double l2  = v2.SquareMagnitude();
-    return (( dot * cos >= 0 ) && 
+    return (( dot * cos >= 0 ) &&
             ( dot * dot ) / l1 / l2 >= ( cos * cos ));
   }
 }
@@ -3050,8 +3294,8 @@ void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
     const double cosTol = Cos( myToler * M_PI / 180. );
     NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks;
 
-    std::list< pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
-    faceQueue.push_back( make_pair( face, myNorm ));
+    std::list< std::pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
+    faceQueue.push_back( std::make_pair( face, myNorm ));
     while ( !faceQueue.empty() )
     {
       face   = faceQueue.front().first;
@@ -3074,7 +3318,7 @@ void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
             if (!normOK || isLessAngle( myNorm, norm, cosTol))
             {
               myCoplanarIDs.Add( f->GetID() );
-              faceQueue.push_back( make_pair( f, norm ));
+              faceQueue.push_back( std::make_pair( f, norm ));
             }
           }
         }
@@ -3216,7 +3460,7 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
   {
     char c = aStr.Value( i );
     if ( !isdigit( c ) && c != ',' && c != '-' )
-      aStr.SetValue( i, ' ');
+      aStr.SetValue( i, ',');
   }
   aStr.RemoveAll( ' ' );
 
@@ -3525,9 +3769,10 @@ void Filter::SetPredicate( PredicatePtr thePredicate )
   myPredicate = thePredicate;
 }
 
-void Filter::GetElementsId( const SMDS_Mesh* theMesh,
-                            PredicatePtr     thePredicate,
-                            TIdSequence&     theSequence )
+void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
+                            PredicatePtr         thePredicate,
+                            TIdSequence&         theSequence,
+                            SMDS_ElemIteratorPtr theElements )
 {
   theSequence.clear();
 
@@ -3536,21 +3781,28 @@ void Filter::GetElementsId( const SMDS_Mesh* theMesh,
 
   thePredicate->SetMesh( theMesh );
 
-  SMDS_ElemIteratorPtr elemIt = theMesh->elementsIterator( thePredicate->GetType() );
-  if ( elemIt ) {
-    while ( elemIt->more() ) {
-      const SMDS_MeshElement* anElem = elemIt->next();
-      long anId = anElem->GetID();
-      if ( thePredicate->IsSatisfy( anId ) )
-        theSequence.push_back( anId );
+  if ( !theElements )
+    theElements = theMesh->elementsIterator( thePredicate->GetType() );
+
+  if ( theElements ) {
+    while ( theElements->more() ) {
+      const SMDS_MeshElement* anElem = theElements->next();
+      if ( thePredicate->GetType() == SMDSAbs_All ||
+           thePredicate->GetType() == anElem->GetType() )
+      {
+        long anId = anElem->GetID();
+        if ( thePredicate->IsSatisfy( anId ) )
+          theSequence.push_back( anId );
+      }
     }
   }
 }
 
 void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
-                            Filter::TIdSequence& theSequence )
+                            Filter::TIdSequence& theSequence,
+                            SMDS_ElemIteratorPtr theElements )
 {
-  GetElementsId(theMesh,myPredicate,theSequence);
+  GetElementsId(theMesh,myPredicate,theSequence,theElements);
 }
 
 /*
@@ -3875,24 +4127,20 @@ void ManifoldPart::expandBoundary
 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
 {
-  std::set<SMDS_MeshCell *> aSetOfFaces;
+
   // take all faces that shared first node
-  SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
-  for ( ; anItr->more(); )
-  {
-    SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
-    if ( !aFace )
-      continue;
-    aSetOfFaces.insert( aFace );
-  }
+  SMDS_ElemIteratorPtr anItr = theLink.myNode1->GetInverseElementIterator( SMDSAbs_Face );
+  SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > faces( anItr ), facesEnd;
+  std::set<const SMDS_MeshElement *> aSetOfFaces( faces, facesEnd );
+
   // take all faces that shared second node
-  anItr = theLink.myNode2->facesIterator();
+  anItr = theLink.myNode2->GetInverseElementIterator( SMDSAbs_Face );
   // find the common part of two sets
   for ( ; anItr->more(); )
   {
-    SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
-    if ( aSetOfFaces.count( aFace ) )
-      theFaces.push_back( aFace );
+    const SMDS_MeshElement* aFace = anItr->next();
+    if ( aSetOfFaces.count( aFace ))
+      theFaces.push_back( (SMDS_MeshFace*) aFace );
   }
 }
 
@@ -3947,9 +4195,9 @@ SMDSAbs_ElementType BelongToMeshGroup::GetType() const
   return myGroup ? myGroup->GetType() : SMDSAbs_All;
 }
 
-/*
-  ElementsOnSurface
-*/
+//================================================================================
+//  ElementsOnSurface
+//================================================================================
 
 ElementsOnSurface::ElementsOnSurface()
 {
@@ -3982,8 +4230,10 @@ SMDSAbs_ElementType ElementsOnSurface::GetType() const
 void ElementsOnSurface::SetTolerance( const double theToler )
 {
   if ( myToler != theToler )
-    myIds.Clear();
-  myToler = theToler;
+  {
+    myToler = theToler;
+    process();
+  }
 }
 
 double ElementsOnSurface::GetTolerance() const
@@ -4083,15 +4333,77 @@ 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; }
+  double Tolerance() const           { return myTol; }
+  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 isOutOfNone  (const gp_Pnt& /*p*/) { return true; }
+  bool isBox        (const TopoDS_Shape& s);
+
+  TopoDS_Shape prepareSolid( const TopoDS_Shape& theSolid );
+
+  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 );
+  size_t GetSize();
+
+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)
 {
 }
 
@@ -4100,6 +4412,40 @@ ElementsOnShape::~ElementsOnShape()
   clearClassifiers();
 }
 
+Predicate* ElementsOnShape::clone() const
+{
+  size_t size = sizeof( *this );
+  if ( myOctree )
+    size += myOctree->GetSize();
+  if ( !myClassifiers.empty() )
+    size += sizeof( myClassifiers[0] ) * myClassifiers.size();
+  if ( !myWorkClassifiers.empty() )
+    size += sizeof( myWorkClassifiers[0] ) * myWorkClassifiers.size();
+  if ( size > 1e+9 ) // 1G
+  {
+#ifdef _DEBUG_
+    std::cout << "Avoid ElementsOnShape::clone(), too large: " << size << " bytes " << std::endl;
+#endif
+    return 0;
+  }
+
+  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
 {
   return myType;
@@ -4165,27 +4511,32 @@ void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool  isOut )
 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
                                 const SMDSAbs_ElementType theType)
 {
+  bool shapeChanges = ( myShape != theShape );
   myType  = theType;
   myShape = theShape;
   if ( myShape.IsNull() ) return;
 
-  TopTools_IndexedMapOfShape shapesMap;
-  TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
-  TopExp_Explorer sub;
-  for ( int i = 0; i < 4; ++i )
+  if ( shapeChanges )
   {
-    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() );
-  }
+    // 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() );
+    }
 
-  clearClassifiers();
-  myClassifiers.resize( shapesMap.Extent() );
-  for ( int i = 0; i < shapesMap.Extent(); ++i )
-    myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler );
+    clearClassifiers();
+    myClassifiers.resize( shapesMap.Extent() );
+    for ( int i = 0; i < shapesMap.Extent(); ++i )
+      myClassifiers[ i ].Init( shapesMap( i+1 ), myToler );
+  }
 
   if ( theType == SMDSAbs_Node )
   {
@@ -4200,125 +4551,295 @@ void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
 
 void ElementsOnShape::clearClassifiers()
 {
-  for ( size_t i = 0; i < myClassifiers.size(); ++i )
-    delete myClassifiers[ i ];
+  // for ( size_t i = 0; i < myClassifiers.size(); ++i )
+  //   delete myClassifiers[ i ];
   myClassifiers.clear();
+
+  delete myOctree;
+  myOctree = 0;
 }
 
-bool ElementsOnShape::IsSatisfy (long elemId)
+bool ElementsOnShape::IsSatisfy( long elemId )
 {
-  const SMDS_Mesh*        mesh = myMeshModifTracer.GetMesh();
-  const SMDS_MeshElement* elem =
-    ( myType == SMDSAbs_Node ? mesh->FindNode( elemId ) : mesh->FindElement( elemId ));
-  if ( !elem || myClassifiers.empty() )
+  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;
 
   gp_XYZ centerXYZ (0, 0, 0);
 
-  SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
-  while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
+  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 );
+
+    SMESHUtils::FreeVector( myWorkClassifiers );
+  }
+
+  for ( int i = 0, nb = elem->NbNodes(); i < nb  && (isSatisfy == myAllNodesFlag); ++i )
   {
-    SMESH_TNodeXYZ aPnt( aNodeItr->next() );
+    SMESH_TNodeXYZ aPnt( elem->GetNode( i ));
     centerXYZ += aPnt;
 
     isNodeOut = true;
     if ( !getNodeIsOut( aPnt._node, isNodeOut ))
     {
-      for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i )
-        isNodeOut = myClassifiers[i]->IsOut( aPnt );
+      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 );
     }
     isSatisfy = !isNodeOut;
   }
 
   // Check the center point for volumes MantisBug 0020168
-  if (isSatisfy &&
-      myAllNodesFlag &&
-      myClassifiers[0]->ShapeType() == TopAbs_SOLID)
+  if ( isSatisfy &&
+       myAllNodesFlag &&
+       myClassifiers[0].ShapeType() == TopAbs_SOLID )
   {
     centerXYZ /= elem->NbNodes();
     isSatisfy = false;
-    for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i )
-      isSatisfy = ! myClassifiers[i]->IsOut( centerXYZ );
+    if ( myOctree )
+    {
+      myWorkClassifiers.clear();
+      myOctree->GetClassifiersAtPoint( centerXYZ, myWorkClassifiers );
+      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;
 }
 
-TopAbs_ShapeEnum ElementsOnShape::TClassifier::ShapeType() const
-{
-  return myShape.ShapeType();
-}
+//================================================================================
+/*!
+ * \brief Check and optionally return a satisfying shape
+ */
+//================================================================================
 
-bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p)
+bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node,
+                                 TopoDS_Shape*        okShape)
 {
-  return (this->*myIsOutFun)( p );
+  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 = myClassifiers[i].Shape();
+          break;
+        }
+    }
+    setNodeIsOut( node, isNodeOut );
+  }
+
+  return !isNodeOut;
 }
 
-void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double theTol)
+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 ( isBox( theShape ))
+  case TopAbs_SOLID:
+  {
+    if (( isShapeBox = isBox( theShape )))
     {
-      myIsOutFun = & ElementsOnShape::TClassifier::isOutOfBox;
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfBox;
     }
     else
     {
-      mySolidClfr.Load(theShape);
-      myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid;
+      mySolidClfr = new BRepClass3d_SolidClassifier( prepareSolid( theShape ));
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid;
     }
     break;
   }
-  case TopAbs_FACE:  {
+  case TopAbs_FACE:
+  {
     Standard_Real u1,u2,v1,v2;
     Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
-    surf->Bounds( u1,u2,v1,v2 );
-    myProjFace.Init(surf, u1,u2, v1,v2, myTol );
-    myIsOutFun = & ElementsOnShape::TClassifier::isOutOfFace;
+    if ( surf.IsNull() )
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
+    else
+    {
+      surf->Bounds( u1,u2,v1,v2 );
+      myProjFace.Init(surf, u1,u2, v1,v2, myTol );
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace;
+    }
     break;
   }
-  case TopAbs_EDGE:  {
+  case TopAbs_EDGE:
+  {
     Standard_Real u1, u2;
-    Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge(theShape), u1, u2);
-    myProjEdge.Init(curve, u1, u2);
-    myIsOutFun = & ElementsOnShape::TClassifier::isOutOfEdge;
+    Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2);
+    if ( curve.IsNull() )
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
+    else
+    {
+      myProjEdge.Init(curve, u1, u2);
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
+    }
     break;
   }
-  case TopAbs_VERTEX:{
+  case TopAbs_VERTEX:
+  {
     myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
-    myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex;
+    myIsOutFun = & ElementsOnShape::Classifier::isOutOfVertex;
     break;
   }
   default:
-    throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier");
+    throw SALOME_Exception("Programmer error in usage of ElementsOnShape::Classifier");
+  }
+
+  if ( !isShapeBox )
+  {
+    if ( theBox )
+    {
+      myBox = *theBox;
+    }
+    else
+    {
+      Bnd_Box box;
+      if ( myShape.ShapeType() == TopAbs_FACE )
+      {
+        BRepAdaptor_Surface SA( TopoDS::Face( myShape ), /*useBoundaries=*/false );
+        if ( SA.GetType() == GeomAbs_BSplineSurface )
+          BRepBndLib::AddOptimal( myShape, box,
+                                  /*useTriangulation=*/true, /*useShapeTolerance=*/true );
+      }
+      if ( box.IsVoid() )
+        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 );
+    }
   }
 }
 
-bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p)
+ElementsOnShape::Classifier::~Classifier()
+{
+  delete mySolidClfr; mySolidClfr = 0;
+}
+
+TopoDS_Shape ElementsOnShape::Classifier::prepareSolid( const TopoDS_Shape& theSolid )
 {
-  mySolidClfr.Perform( p, myTol );
-  return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON );
+  // try to limit tolerance of theSolid down to myTol (issue #19026)
+
+  // check if tolerance of theSolid is more than myTol
+  bool tolIsOk = true; // max tolerance is at VERTEXes
+  for ( TopExp_Explorer exp( theSolid, TopAbs_VERTEX ); exp.More() &&  tolIsOk; exp.Next() )
+    tolIsOk = ( myTol >= BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current() )));
+  if ( tolIsOk )
+    return theSolid;
+
+  // make a copy to prevent the original shape from changes
+  TopoDS_Shape resultShape = BRepBuilderAPI_Copy( theSolid );
+
+  if ( !GEOMUtils::FixShapeTolerance( resultShape, TopAbs_SHAPE, myTol ))
+    return theSolid;
+  return resultShape;
+}
+
+bool ElementsOnShape::Classifier::isOutOfSolid( const gp_Pnt& p )
+{
+  if ( isOutOfBox( p )) return true;
+  mySolidClfr->Perform( p, myTol );
+  return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON );
 }
 
-bool ElementsOnShape::TClassifier::isOutOfBox (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfBox( const gp_Pnt& p )
 {
   return myBox.IsOut( p.XYZ() );
 }
 
-bool ElementsOnShape::TClassifier::isOutOfFace  (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfFace( const gp_Pnt& p )
 {
+  if ( isOutOfBox( p )) return true;
   myProjFace.Perform( p );
   if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
   {
     // check relatively to the face
-    Quantity_Parameter u, v;
+    Standard_Real u, v;
     myProjFace.LowerDistanceParameters(u, v);
     gp_Pnt2d aProjPnt (u, v);
     BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
@@ -4328,18 +4849,19 @@ bool ElementsOnShape::TClassifier::isOutOfFace  (const gp_Pnt& p)
   return true;
 }
 
-bool ElementsOnShape::TClassifier::isOutOfEdge  (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfEdge( const gp_Pnt& p )
 {
+  if ( isOutOfBox( p )) return true;
   myProjEdge.Perform( p );
   return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
 }
 
-bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfVertex( const gp_Pnt& p )
 {
   return ( myVertexXYZ.Distance( p ) > myTol );
 }
 
-bool ElementsOnShape::TClassifier::isBox (const TopoDS_Shape& theShape)
+bool ElementsOnShape::Classifier::isBox(const TopoDS_Shape& theShape )
 {
   TopTools_IndexedMapOfShape vMap;
   TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
@@ -4366,6 +4888,132 @@ bool ElementsOnShape::TClassifier::isBox (const TopoDS_Shape& theShape)
   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 );
+  }
+}
+
+size_t ElementsOnShape::OctreeClassifier::GetSize()
+{
+  size_t res = sizeof( *this );
+  if ( !myClassifiers.empty() )
+    res += sizeof( myClassifiers[0] ) * myClassifiers.size();
+
+  if ( !isLeaf() )
+    for (int i = 0; i < nbChildren(); i++)
+      res += ((OctreeClassifier*) myChildren[i])->GetSize();
+
+  return res;
+}
+
+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 };
+
+  for ( size_t i = 0; i < myClassifiers.size(); ++i )
+  {
+    for ( int j = 0; j < nbChildren(); j++ )
+    {
+      if ( !myClassifiers[i]->GetBndBox()->IsOut( *myChildren[j]->getBox() ))
+      {
+        myClassifiers[i]->SetFlag( childFlag[ j ]);
+        ++nbInChild[ j ];
+      }
+    }
+  }
+
+  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 ]))
+      {
+        --nbInChild[ j ];
+        child->myClassifiers[ nbInChild[ j ]] = myClassifiers[ i ];
+        myClassifiers[ i ]->UnsetFlag( childFlag[ j ]);
+      }
+    }
+  }
+  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  ||
+                        child->maxSize() < child->myClassifiers[0]->Tolerance() );
+  }
+}
+
+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
@@ -4375,25 +5023,45 @@ bool ElementsOnShape::TClassifier::isBox (const TopoDS_Shape& theShape)
 
 BelongToGeom::BelongToGeom()
   : myMeshDS(NULL),
-    myType(SMDSAbs_All),
+    myType(SMDSAbs_NbElementTypes),
     myIsSubshape(false),
     myTolerance(Precision::Confusion())
 {}
 
+Predicate* BelongToGeom::clone() const
+{
+  BelongToGeom* cln = 0;
+  if ( myElementsOnShapePtr )
+    if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
+    {
+      cln = new BelongToGeom( *this );
+      cln->myElementsOnShapePtr.reset( eos );
+    }
+  return cln;
+}
+
 void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
 {
-  myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
-  init();
+  if ( myMeshDS != theMesh )
+  {
+    myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
+    init();
+  }
+  if ( myElementsOnShapePtr )
+    myElementsOnShapePtr->SetMesh( myMeshDS );
 }
 
 void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
 {
-  myShape = theShape;
-  init();
+  if ( myShape != theShape )
+  {
+    myShape = theShape;
+    init();
+  }
 }
 
 static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
-                        const TopoDS_Shape& theShape)
+                        const TopoDS_Shape&               theShape)
 {
   if (theMap.Contains(theShape)) return true;
 
@@ -4415,7 +5083,7 @@ static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
 
 void BelongToGeom::init()
 {
-  if (!myMeshDS || myShape.IsNull()) return;
+  if ( !myMeshDS || myShape.IsNull() ) return;
 
   // is sub-shape of main shape?
   TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
@@ -4424,38 +5092,31 @@ void BelongToGeom::init()
   }
   else {
     TopTools_IndexedMapOfShape aMap;
-    TopExp::MapShapes(aMainShape, aMap);
-    myIsSubshape = IsSubShape(aMap, myShape);
+    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 )
+      {
+        int subID = myMeshDS->ShapeToIndex( aMap( i ));
+        if ( subID > 0 )
+          mySubShapesIDs.Add( subID );
+      }
+    }
   }
 
   //if (!myIsSubshape) // to be always ready to check an element not bound to geometry
   {
-    myElementsOnShapePtr.reset(new ElementsOnShape());
-    myElementsOnShapePtr->SetTolerance(myTolerance);
-    myElementsOnShapePtr->SetAllNodes(true); // "belong", while false means "lays on"
-    myElementsOnShapePtr->SetMesh(myMeshDS);
-    myElementsOnShapePtr->SetShape(myShape, myType);
-  }
-}
-
-static bool IsContains( const SMESHDS_Mesh*     theMeshDS,
-                        const TopoDS_Shape&     theShape,
-                        const SMDS_MeshElement* theElem,
-                        TopAbs_ShapeEnum        theFindShapeEnum,
-                        TopAbs_ShapeEnum        theAvoidShapeEnum = TopAbs_SHAPE )
-{
-  TopExp_Explorer anExp( theShape,theFindShapeEnum,theAvoidShapeEnum );
-
-  while( anExp.More() )
-  {
-    const TopoDS_Shape& aShape = anExp.Current();
-    if( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ) ){
-      if( aSubMesh->Contains( theElem ) )
-        return true;
-    }
-    anExp.Next();
+    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 );
   }
-  return false;
 }
 
 bool BelongToGeom::IsSatisfy (long theId)
@@ -4468,51 +5129,28 @@ bool BelongToGeom::IsSatisfy (long theId)
     return myElementsOnShapePtr->IsSatisfy(theId);
   }
 
-  // Case of submesh
+  // Case of sub-mesh
+
   if (myType == SMDSAbs_Node)
   {
-    if( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ) )
+    if ( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ))
     {
       if ( aNode->getshapeId() < 1 )
         return myElementsOnShapePtr->IsSatisfy(theId);
-
-      const SMDS_PositionPtr& aPosition = aNode->GetPosition();
-      SMDS_TypeOfPosition aTypeOfPosition = aPosition->GetTypeOfPosition();
-      switch( aTypeOfPosition )
-      {
-      case SMDS_TOP_VERTEX : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_VERTEX ));
-      case SMDS_TOP_EDGE   : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_EDGE ));
-      case SMDS_TOP_FACE   : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_FACE ));
-      case SMDS_TOP_3DSPACE: return ( IsContains( myMeshDS,myShape,aNode,TopAbs_SOLID ) ||
-                                      IsContains( myMeshDS,myShape,aNode,TopAbs_SHELL ));
-      default:;
-      }
+      else
+        return mySubShapesIDs.Contains( aNode->getshapeId() );
     }
   }
   else
   {
     if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
     {
-      if ( anElem->getshapeId() < 1 )
-        return myElementsOnShapePtr->IsSatisfy(theId);
-
-      if( myType == SMDSAbs_All )
-      {
-        return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ) ||
-                 IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ) ||
-                 IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
-                 IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
-      }
-      else if( myType == anElem->GetType() )
+      if ( myType == SMDSAbs_All || anElem->GetType() == myType )
       {
-        switch( myType )
-        {
-        case SMDSAbs_Edge  : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ));
-        case SMDSAbs_Face  : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ));
-        case SMDSAbs_Volume: return ( IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
-                                      IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
-        default:;
-        }
+        if ( anElem->getshapeId() < 1 )
+          return myElementsOnShapePtr->IsSatisfy(theId);
+        else
+          return mySubShapesIDs.Contains( anElem->getshapeId() );
       }
     }
   }
@@ -4522,8 +5160,11 @@ bool BelongToGeom::IsSatisfy (long theId)
 
 void BelongToGeom::SetType (SMDSAbs_ElementType theType)
 {
-  myType = theType;
-  init();
+  if ( myType != theType )
+  {
+    myType = theType;
+    init();
+  }
 }
 
 SMDSAbs_ElementType BelongToGeom::GetType() const
@@ -4544,8 +5185,7 @@ const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const
 void BelongToGeom::SetTolerance (double theTolerance)
 {
   myTolerance = theTolerance;
-  if (!myIsSubshape)
-    init();
+  init();
 }
 
 double BelongToGeom::GetTolerance()
@@ -4556,26 +5196,46 @@ double BelongToGeom::GetTolerance()
 /*
   Class       : LyingOnGeom
   Description : Predicate for verifying whether entiy lying or partially lying on
-                specified geometrical support
+  specified geometrical support
 */
 
 LyingOnGeom::LyingOnGeom()
   : myMeshDS(NULL),
-    myType(SMDSAbs_All),
+    myType(SMDSAbs_NbElementTypes),
     myIsSubshape(false),
     myTolerance(Precision::Confusion())
 {}
 
+Predicate* LyingOnGeom::clone() const
+{
+  LyingOnGeom* cln = 0;
+  if ( myElementsOnShapePtr )
+    if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
+    {
+      cln = new LyingOnGeom( *this );
+      cln->myElementsOnShapePtr.reset( eos );
+    }
+  return cln;
+}
+
 void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
 {
-  myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
-  init();
+  if ( myMeshDS != theMesh )
+  {
+    myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
+    init();
+  }
+  if ( myElementsOnShapePtr )
+    myElementsOnShapePtr->SetMesh( myMeshDS );
 }
 
 void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
 {
-  myShape = theShape;
-  init();
+  if ( myShape != theShape )
+  {
+    myShape = theShape;
+    init();
+  }
 }
 
 void LyingOnGeom::init()
@@ -4603,13 +5263,14 @@ void LyingOnGeom::init()
         mySubShapesIDs.Add( subID );
     }
   }
-  else
+  // else // to be always ready to check an element not bound to geometry
   {
-    myElementsOnShapePtr.reset(new ElementsOnShape());
-    myElementsOnShapePtr->SetTolerance(myTolerance);
-    myElementsOnShapePtr->SetAllNodes(false); // lays on, while true means "belong"
-    myElementsOnShapePtr->SetMesh(myMeshDS);
-    myElementsOnShapePtr->SetShape(myShape, myType);
+    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 );
   }
 }
 
@@ -4631,7 +5292,8 @@ bool LyingOnGeom::IsSatisfy( long theId )
   if ( mySubShapesIDs.Contains( elem->getshapeId() ))
     return true;
 
-  if ( elem->GetType() != SMDSAbs_Node )
+  if (( elem->GetType() != SMDSAbs_Node ) &&
+      ( myType == SMDSAbs_All || elem->GetType() == myType ))
   {
     SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator();
     while ( nodeItr->more() )
@@ -4647,8 +5309,11 @@ bool LyingOnGeom::IsSatisfy( long theId )
 
 void LyingOnGeom::SetType( SMDSAbs_ElementType theType )
 {
-  myType = theType;
-  init();
+  if ( myType != theType )
+  {
+    myType = theType;
+    init();
+  }
 }
 
 SMDSAbs_ElementType LyingOnGeom::GetType() const
@@ -4669,8 +5334,7 @@ const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const
 void LyingOnGeom::SetTolerance (double theTolerance)
 {
   myTolerance = theTolerance;
-  if (!myIsSubshape)
-    init();
+  init();
 }
 
 double LyingOnGeom::GetTolerance()
@@ -4678,39 +5342,6 @@ double LyingOnGeom::GetTolerance()
   return myTolerance;
 }
 
-bool LyingOnGeom::Contains( const SMESHDS_Mesh*     theMeshDS,
-                            const TopoDS_Shape&     theShape,
-                            const SMDS_MeshElement* theElem,
-                            TopAbs_ShapeEnum        theFindShapeEnum,
-                            TopAbs_ShapeEnum        theAvoidShapeEnum )
-{
-  // if (IsContains(theMeshDS, theShape, theElem, theFindShapeEnum, theAvoidShapeEnum))
-  //   return true;
-
-  // TopTools_MapOfShape aSubShapes;
-  // TopExp_Explorer exp( theShape, theFindShapeEnum, theAvoidShapeEnum );
-  // for ( ; exp.More(); exp.Next() )
-  // {
-  //   const TopoDS_Shape& aShape = exp.Current();
-  //   if ( !aSubShapes.Add( aShape )) continue;
-
-  //   if ( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ))
-  //   {
-  //     if ( aSubMesh->Contains( theElem ))
-  //       return true;
-
-  //     SMDS_ElemIteratorPtr nodeItr = theElem->nodesIterator();
-  //     while ( nodeItr->more() )
-  //     {
-  //       const SMDS_MeshElement* aNode = nodeItr->next();
-  //       if ( aSubMesh->Contains( aNode ))
-  //         return true;
-  //     }
-  //   }
-  // }
-  return false;
-}
-
 TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0)
 {}