Salome HOME
Join modifications from branch OCC_development_for_3_2_0a2
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
index 2bf885a7a29a0d1500bc03e04f0818c4e2eda230..c00cbae43554bcf0c5ea8f40f2d46537bad8eed7 100644 (file)
@@ -1,20 +1,20 @@
 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-//  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS 
-// 
-//  This library is free software; you can redistribute it and/or 
-//  modify it under the terms of the GNU Lesser General Public 
+//  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
 //  License as published by the Free Software Foundation; either
-//  version 2.1 of the License. 
-// 
-//  This library is distributed in the hope that it will be useful, 
+//  version 2.1 of the License.
+//
+//  This library is distributed in the hope that it will be useful,
 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 //  Lesser General Public License for more details.
-// 
-//  You should have received a copy of the GNU Lesser General Public 
-//  License along with this library; if not, write to the Free Software 
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA 
-// 
+//
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
 
 #include "SMESH_ControlsDef.hxx"
 #include "SMDS_Iterator.hxx"
 #include "SMDS_MeshElement.hxx"
 #include "SMDS_MeshNode.hxx"
+#include "SMDS_VolumeTool.hxx"
+#include "SMDS_QuadraticFaceOfNodes.hxx"
+#include "SMDS_QuadraticEdge.hxx"
 
 
 /*
-                            AUXILIARY METHODS 
+                            AUXILIARY METHODS
 */
 
 namespace{
   inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
   {
     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
-    
+
     return v1.Magnitude() < gp::Resolution() ||
       v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
   }
@@ -80,39 +83,76 @@ namespace{
     return aDist;
   }
 
-  int getNbMultiConnection( SMDS_Mesh* theMesh, const int theId )
+  int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId )
   {
     if ( theMesh == 0 )
       return 0;
-    
+
     const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
-    if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 )
+    if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */)
       return 0;
-    
-    TColStd_MapOfInteger aMap;
-    
-    int aResult = 0;
-    SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
-    if ( anIter != 0 ) {
-      while( anIter->more() ) {
-       const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
-       if ( aNode == 0 )
-         return 0;
-       SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
-       while( anElemIter->more() ) {
-         const SMDS_MeshElement* anElem = anElemIter->next();
-         if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
-           int anId = anElem->GetID();
-           
-           if ( anIter->more() )              // i.e. first node
-             aMap.Add( anId );
-           else if ( aMap.Contains( anId ) )
-             aResult++;
-         }
-       }
+
+    // for each pair of nodes in anEdge (there are 2 pairs in a quadratic edge)
+    // count elements containing both nodes of the pair.
+    // Note that there may be such cases for a quadratic edge (a horizontal line):
+    //
+    //  Case 1          Case 2
+    //  |     |      |        |      |
+    //  |     |      |        |      |
+    //  +-----+------+  +-----+------+ 
+    //  |            |  |            |
+    //  |            |  |            |
+    // result sould be 2 in both cases
+    //
+    int aResult0 = 0, aResult1 = 0;
+     // last node, it is a medium one in a quadratic edge
+    const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
+    const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 );
+    const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 );
+    if ( aNode1 == aLastNode ) aNode1 = 0;
+
+    SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator();
+    while( anElemIter->more() ) {
+      const SMDS_MeshElement* anElem = anElemIter->next();
+      if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
+        SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
+        while ( anIter->more() ) {
+          if ( const SMDS_MeshElement* anElemNode = anIter->next() ) {
+            if ( anElemNode == aNode0 ) {
+              aResult0++;
+              if ( !aNode1 ) break; // not a quadratic edge
+            }
+            else if ( anElemNode == aNode1 )
+              aResult1++;
+          }
+        }
       }
     }
-    
+    int aResult = max ( aResult0, aResult1 );
+
+//     TColStd_MapOfInteger aMap;
+
+//     SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
+//     if ( anIter != 0 ) {
+//       while( anIter->more() ) {
+//     const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
+//     if ( aNode == 0 )
+//       return 0;
+//     SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
+//     while( anElemIter->more() ) {
+//       const SMDS_MeshElement* anElem = anElemIter->next();
+//       if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
+//         int anId = anElem->GetID();
+
+//         if ( anIter->more() )              // i.e. first node
+//           aMap.Add( anId );
+//         else if ( aMap.Contains( anId ) )
+//           aResult++;
+//       }
+//     }
+//       }
+//     }
+
     return aResult;
   }
 
@@ -136,7 +176,7 @@ NumericalFunctor::NumericalFunctor():
   myPrecision = -1;
 }
 
-void NumericalFunctor::SetMesh( SMDS_Mesh* theMesh )
+void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh )
 {
   myMesh = theMesh;
 }
@@ -153,23 +193,41 @@ bool NumericalFunctor::GetPoints(const int theId,
 }
 
 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
-                                 TSequenceOfXYZ& theRes )
+                                 TSequenceOfXYZ&         theRes )
 {
   theRes.clear();
 
   if ( anElem == 0)
     return false;
 
+  theRes.reserve( anElem->NbNodes() );
+
   // Get nodes of the element
-  SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
-  if ( anIter != 0 )
-  {
-    while( anIter->more() )
-    {
-      const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
-      if ( aNode != 0 ){
+  SMDS_ElemIteratorPtr anIter;
+
+  if ( anElem->IsQuadratic() ) {
+    switch ( anElem->GetType() ) {
+    case SMDSAbs_Edge:
+      anIter = static_cast<const SMDS_QuadraticEdge*>
+        (anElem)->interlacedNodesElemIterator();
+      break;
+    case SMDSAbs_Face:
+      anIter = static_cast<const SMDS_QuadraticFaceOfNodes*>
+        (anElem)->interlacedNodesElemIterator();
+      break;
+    default:
+      anIter = anElem->nodesIterator();
+      //return false;
+    }
+  }
+  else {
+    anIter = anElem->nodesIterator();
+  }
+
+  if ( anIter ) {
+    while( anIter->more() ) {
+      if ( const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( anIter->next() ))
         theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
-      }
     }
   }
 
@@ -203,6 +261,42 @@ double NumericalFunctor::GetValue( long theId )
   return 0.;
 }
 
+//=======================================================================
+//function : GetValue
+//purpose  : 
+//=======================================================================
+
+double Volume::GetValue( long theElementId )
+{
+  if ( theElementId && myMesh ) {
+    SMDS_VolumeTool aVolumeTool;
+    if ( aVolumeTool.Set( myMesh->FindElement( theElementId )))
+      return aVolumeTool.GetSize();
+  }
+  return 0;
+}
+
+//=======================================================================
+//function : GetBadRate
+//purpose  : meaningless as it is not quality control functor
+//=======================================================================
+
+double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+  return Value;
+}
+
+//=======================================================================
+//function : GetType
+//purpose  : 
+//=======================================================================
+
+SMDSAbs_ElementType Volume::GetType() const
+{
+  return SMDSAbs_Volume;
+}
+
+
 /*
   Class       : MinimumAngle
   Description : Functor for calculation of minimum angle
@@ -212,32 +306,24 @@ double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
 {
   double aMin;
 
-  if ( P.size() == 3 )
-  {
-    double A0 = getAngle( P( 3 ), P( 1 ), P( 2 ) );
-    double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
-    double A2 = getAngle( P( 2 ), P( 3 ), P( 1 ) );
+  if (P.size() <3)
+    return 0.;
 
-    aMin = Min( A0, Min( A1, A2 ) );
-  }
-  else if ( P.size() == 4 )
-  {
-    double A0 = getAngle( P( 4 ), P( 1 ), P( 2 ) );
-    double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
-    double A2 = getAngle( P( 2 ), P( 3 ), P( 4 ) );
-    double A3 = getAngle( P( 3 ), P( 4 ), P( 1 ) );
-    
-    aMin = Min( Min( A0, A1 ), Min( A2, A3 ) );
+  aMin = getAngle(P( P.size() ), P( 1 ), P( 2 ));
+  aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 )));
+
+  for (int i=2; i<P.size();i++){
+      double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) );
+    aMin = Min(aMin,A0);
   }
-  else
-    return 0.;
-  
-  return aMin * 180 / PI;
+
+  return aMin * 180.0 / PI;
 }
 
 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
 {
-  const double aBestAngle = PI / nbNodes;
+  //const double aBestAngle = PI / nbNodes;
+  const double aBestAngle = 180.0 - ( 360.0 / double(nbNodes) );
   return ( fabs( aBestAngle - Value ));
 }
 
@@ -253,38 +339,65 @@ SMDSAbs_ElementType MinimumAngle::GetType() const
 */
 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
 {
+  // According to "Mesh quality control" by Nadir Bouhamau referring to
+  // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis.
+  // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4
+  // PAL10872
+
   int nbNodes = P.size();
 
-  if ( nbNodes != 3 && nbNodes != 4 )
+  if ( nbNodes < 3 )
     return 0;
 
   // Compute lengths of the sides
 
-  double aLen[ nbNodes ];
+  vector< double > aLen (nbNodes);
+
   for ( int i = 0; i < nbNodes - 1; i++ )
     aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
   aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
 
   // Compute aspect ratio
 
-  if ( nbNodes == 3 ) 
+  if ( nbNodes == 3 )
   {
+    // Q = alfa * h * p / S, where
+    //
+    // alfa = sqrt( 3 ) / 6
+    // h - length of the longest edge
+    // p - half perimeter
+    // S - triangle surface
+
+    const double alfa = sqrt( 3. ) / 6.;
+    double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
+    double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
     double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
     if ( anArea <= Precision::Confusion() )
       return 0.;
-    double aMaxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
-    static double aCoef = sqrt( 3. ) / 4;
 
-    return aCoef * aMaxLen * aMaxLen / anArea;
+    return alfa * maxLen * half_perimeter / anArea;
   }
   else
   {
-    double aMinLen = Min( Min( aLen[ 0 ], aLen[ 1 ] ), Min( aLen[ 2 ], aLen[ 3 ] ) );
-    if ( aMinLen <= Precision::Confusion() )
-      return 0.;
-    double aMaxLen = Max( Max( aLen[ 0 ], aLen[ 1 ] ), Max( aLen[ 2 ], aLen[ 3 ] ) );
-    
-    return aMaxLen / aMinLen;
+    // return aspect ratio of the worst triange which can be built
+    // taking three nodes of the quadrangle
+    TSequenceOfXYZ triaPnts(3);
+    // triangle on nodes 1 3 2
+    triaPnts(1) = P(1);
+    triaPnts(2) = P(3);
+    triaPnts(3) = P(2);
+    double ar = GetValue( triaPnts );
+    // triangle on nodes 1 3 4
+    triaPnts(3) = P(4);
+    ar = Max ( ar, GetValue( triaPnts ));
+    // triangle on nodes 1 2 4
+    triaPnts(2) = P(2);
+    ar = Max ( ar, GetValue( triaPnts ));
+    // triangle on nodes 3 2 4
+    triaPnts(1) = P(3);
+    ar = Max ( ar, GetValue( triaPnts ));
+
+    return ar;
   }
 }
 
@@ -344,7 +457,7 @@ namespace{
     double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
     double R = c2*d2*(a2+b2+e2+f2-c2-d2);
     double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
-    
+
     return sqrt(P+Q+R-S)/12.0;
   }
 
@@ -353,7 +466,7 @@ namespace{
     gp_Vec aVec2( P( 3 ) - P( 1 ) );
     gp_Vec aVec3( P( 4 ) - P( 1 ) );
     gp_Vec anAreaVec( aVec1 ^ aVec2 );
-    return abs(aVec3 * anAreaVec) / 6.0;
+    return fabs(aVec3 * anAreaVec) / 6.0;
   }
 
   inline double getMaxHeight(double theLen[6])
@@ -404,7 +517,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
     double aVolume = getVolume(P);
     //double aVolume = getVolume(aLen);
     double aHeight = getMaxHeight(aLen);
-    static double aCoeff = sqrt(6.0)/36.0;
+    static double aCoeff = sqrt(2.0)/12.0;
     aQuality = aCoeff*aHeight*aSumArea/aVolume;
     break;
   }
@@ -590,6 +703,21 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
     break;
   }
   }
+  if ( nbNodes > 4 ) {
+    // avaluate aspect ratio of quadranle faces
+    AspectRatio aspect2D;
+    SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
+    int nbFaces = SMDS_VolumeTool::NbFaces( type );
+    TSequenceOfXYZ points(4);
+    for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume
+      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
+        points( p + 1 ) = P( pInd[ p ] + 1 );
+      aQuality = max( aQuality, aspect2D.GetValue( points ));
+    }
+  }
   return aQuality;
 }
 
@@ -637,8 +765,8 @@ double Warping::ComputeA( const gp_XYZ& thePnt1,
   if ( L < Precision::Confusion())
     return 0.;
 
-  gp_XYZ GI = ( thePnt2 - thePnt1 ) / 2. - theG;
-  gp_XYZ GJ = ( thePnt3 - thePnt2 ) / 2. - theG;
+  gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
+  gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
   gp_XYZ N  = GI.Crossed( GJ );
 
   if ( N.Modulus() < gp::Resolution() )
@@ -735,7 +863,7 @@ double Skew::GetValue( const TSequenceOfXYZ& P )
 
     return Max( A0, Max( A1, A2 ) ) * 180 / PI;
   }
-  else 
+  else
   {
     gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2;
     gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2;
@@ -770,16 +898,21 @@ SMDSAbs_ElementType Skew::GetType() const
 */
 double Area::GetValue( const TSequenceOfXYZ& P )
 {
-  if ( P.size() == 3 )
-    return getArea( P( 1 ), P( 2 ), P( 3 ) );
-  else if ( P.size() == 4 )
-    return getArea( P( 1 ), P( 2 ), P( 3 ) ) + getArea( P( 1 ), P( 3 ), P( 4 ) );
-  else
-    return 0;
+  gp_Vec aVec1( P(2) - P(1) );
+  gp_Vec aVec2( P(3) - P(1) );
+  gp_Vec SumVec = aVec1 ^ aVec2;
+  for (int i=4; i<=P.size(); i++) {
+    gp_Vec aVec1( P(i-1) - P(1) );
+    gp_Vec aVec2( P(i) - P(1) );
+    gp_Vec tmp = aVec1 ^ aVec2;
+    SumVec.Add(tmp);
+  }
+  return SumVec.Magnitude() * 0.5;
 }
 
 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
 {
+  // meaningless as it is not a quality control functor
   return Value;
 }
 
@@ -795,11 +928,16 @@ SMDSAbs_ElementType Area::GetType() const
 */
 double Length::GetValue( const TSequenceOfXYZ& P )
 {
-  return ( P.size() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0 );
+  switch ( P.size() ) {
+  case 2:  return getDistance( P( 1 ), P( 2 ) );
+  case 3:  return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) );
+  default: return 0.;
+  }
 }
 
 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
 {
+  // meaningless as it is not quality control functor
   return Value;
 }
 
@@ -817,21 +955,28 @@ double Length2D::GetValue( long theElementId)
 {
   TSequenceOfXYZ P;
 
+  //cout<<"Length2D::GetValue"<<endl;
   if (GetPoints(theElementId,P)){
-    
+    //for(int jj=1; jj<=P.size(); jj++)
+    //  cout<<"jj="<<jj<<" P("<<P(jj).X()<<","<<P(jj).Y()<<","<<P(jj).Z()<<")"<<endl;
+
     double  aVal;// = GetValue( P );
     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
     SMDSAbs_ElementType aType = aElem->GetType();
-    
+
     int len = P.size();
-    
+
     switch (aType){
     case SMDSAbs_All:
-    case SMDSAbs_Node: 
+    case SMDSAbs_Node:
     case SMDSAbs_Edge:
       if (len == 2){
        aVal = getDistance( P( 1 ), P( 2 ) );
-       break;
+        break;
+      }
+      else if (len == 3){ // quadratic edge
+       aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
+        break;
       }
     case SMDSAbs_Face:
       if (len == 3){ // triangles
@@ -849,6 +994,22 @@ double Length2D::GetValue( long theElementId)
        aVal = Max(Max(L1,L2),Max(L3,L4));
        break;
       }
+      if (len == 6){ // quadratic triangles
+       double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+       double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+       double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
+       aVal = Max(L1,Max(L2,L3));
+        //cout<<"L1="<<L1<<" L2="<<L2<<"L3="<<L3<<" aVal="<<aVal<<endl;
+       break;
+      }
+      else if (len == 8){ // quadratic quadrangles
+       double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+       double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+       double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
+       double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
+       aVal = Max(Max(L1,L2),Max(L3,L4));
+       break;
+      }
     case SMDSAbs_Volume:
       if (len == 4){ // tetraidrs
        double L1 = getDistance(P( 1 ),P( 2 ));
@@ -859,7 +1020,7 @@ double Length2D::GetValue( long theElementId)
        double L6 = getDistance(P( 3 ),P( 4 ));
        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
        break;
-      } 
+      }
       else if (len == 5){ // piramids
        double L1 = getDistance(P( 1 ),P( 2 ));
        double L2 = getDistance(P( 2 ),P( 3 ));
@@ -869,7 +1030,7 @@ double Length2D::GetValue( long theElementId)
        double L6 = getDistance(P( 2 ),P( 5 ));
        double L7 = getDistance(P( 3 ),P( 5 ));
        double L8 = getDistance(P( 4 ),P( 5 ));
-      
+
        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
        aVal = Max(aVal,Max(L7,L8));
        break;
@@ -884,7 +1045,7 @@ double Length2D::GetValue( long theElementId)
        double L7 = getDistance(P( 1 ),P( 4 ));
        double L8 = getDistance(P( 2 ),P( 5 ));
        double L9 = getDistance(P( 3 ),P( 6 ));
-      
+
        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
        aVal = Max(aVal,Max(Max(L7,L8),L9));
        break;
@@ -902,17 +1063,74 @@ double Length2D::GetValue( long theElementId)
        double L10= getDistance(P( 2 ),P( 6 ));
        double L11= getDistance(P( 3 ),P( 7 ));
        double L12= getDistance(P( 4 ),P( 8 ));
-      
+
        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
        aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
        aVal = Max(aVal,Max(L11,L12));
        break;
-       
+
       }
-      
-    default: aVal=-1; 
+
+      if (len == 10){ // quadratic tetraidrs
+       double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
+       double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
+       double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
+       double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+       double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
+       double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
+       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+       break;
+      }
+      else if (len == 13){ // quadratic piramids
+       double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
+       double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
+       double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
+       double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+       double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+       double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
+       double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
+       double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
+       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+       aVal = Max(aVal,Max(L7,L8));
+       break;
+      }
+      else if (len == 15){ // quadratic pentaidres
+       double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
+       double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
+       double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+       double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+       double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
+       double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
+       double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
+       double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
+       double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
+       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+       aVal = Max(aVal,Max(Max(L7,L8),L9));
+       break;
+      }
+      else if (len == 20){ // quadratic hexaider
+       double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
+       double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
+       double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
+       double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
+       double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
+       double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
+       double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
+       double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
+       double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
+       double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
+       double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
+       double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
+       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+       aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+       aVal = Max(aVal,Max(L11,L12));
+       break;
+
+      }
+
+    default: aVal=-1;
     }
-    
+
     if (aVal <0){
       return 0.;
     }
@@ -922,7 +1140,7 @@ double Length2D::GetValue( long theElementId)
       double prec = pow( 10., (double)( myPrecision ) );
       aVal = floor( aVal * prec + 0.5 ) / prec;
     }
-    
+
     return aVal;
 
   }
@@ -931,6 +1149,7 @@ double Length2D::GetValue( long theElementId)
 
 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
 {
+  // meaningless as it is not quality control functor
   return Value;
 }
 
@@ -960,38 +1179,81 @@ void Length2D::GetValues(TValues& theValues){
   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
   for(; anIter->more(); ){
     const SMDS_MeshFace* anElem = anIter->next();
-    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];
+    if(anElem->IsQuadratic()) {
+      const SMDS_QuadraticFaceOfNodes* F =
+        static_cast<const SMDS_QuadraticFaceOfNodes*>(anElem);
+      // use special nodes iterator
+      SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
+      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);
     }
-    
-    aLength = P[0].Distance(P[1]);
-    
-    Value aValue(aLength,aNodeId[0],aNodeId[1]);
-    theValues.insert(aValue);
   }
 }
 
@@ -1010,6 +1272,7 @@ double MultiConnection::GetValue( long theId )
 
 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
 {
+  // meaningless as it is not quality control functor
   return Value;
 }
 
@@ -1031,20 +1294,19 @@ double MultiConnection2D::GetValue( long theElementId )
 {
   TSequenceOfXYZ P;
   int aResult = 0;
-  
+
   if (GetPoints(theElementId,P)){
-    double  aVal;
     const SMDS_MeshElement* anFaceElem = myMesh->FindElement( theElementId );
     SMDSAbs_ElementType aType = anFaceElem->GetType();
-    
+
     int len = P.size();
-    
+
     TColStd_MapOfInteger aMap;
     int aResult = 0;
-    
+
     switch (aType){
     case SMDSAbs_All:
-    case SMDSAbs_Node: 
+    case SMDSAbs_Node:
     case SMDSAbs_Edge:
     case SMDSAbs_Face:
       if (len == 3){ // triangles
@@ -1074,20 +1336,21 @@ double MultiConnection2D::GetValue( long theElementId )
            }
          }
        }
-       
+
        aResult = Max(Max(Nb[0],Nb[1]),Nb[2]);
       }
       break;
     case SMDSAbs_Volume:
     default: aResult=0;
     }
-    
+
   }
   return aResult;//getNbMultiConnection( myMesh, theId );
 }
 
 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
 {
+  // meaningless as it is not quality control functor
   return Value;
 }
 
@@ -1115,8 +1378,12 @@ void MultiConnection2D::GetValues(MValues& theValues){
   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
   for(; anIter->more(); ){
     const SMDS_MeshFace* anElem = anIter->next();
-    long anElemId = anElem->GetID();
-    SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
+    SMDS_ElemIteratorPtr aNodesIter;
+    if ( anElem->IsQuadratic() )
+      aNodesIter = static_cast<const SMDS_QuadraticFaceOfNodes*>
+        (anElem)->interlacedNodesElemIterator();
+    else
+      aNodesIter = anElem->nodesIterator();
     long aNodeId[3];
 
     //int aNbConnects=0;
@@ -1129,17 +1396,18 @@ void MultiConnection2D::GetValues(MValues& theValues){
       const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
       aNodeId[0] = aNodeId[1] = aNodes->GetID();
     }
-    for(; aNodesIter->more(); ){
+    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 {
+      }
+      else {
        theValues[aValue] = 1;
        //aNbConnects = 1;
       }
@@ -1149,10 +1417,11 @@ void MultiConnection2D::GetValues(MValues& theValues){
     }
     Value aValue(aNodeId[0],aNodeId[2]);
     MValues::iterator aItr = theValues.find(aValue);
-    if (aItr != theValues.end()){
+    if (aItr != theValues.end()) {
       aItr->second += 1;
       //aNbConnects = nb;
-    } else {
+    }
+    else {
       theValues[aValue] = 1;
       //aNbConnects = 1;
     }
@@ -1165,6 +1434,37 @@ void MultiConnection2D::GetValues(MValues& theValues){
                             PREDICATES
 */
 
+/*
+  Class       : BadOrientedVolume
+  Description : Predicate bad oriented volumes
+*/
+
+BadOrientedVolume::BadOrientedVolume()
+{
+  myMesh = 0;
+}
+
+void BadOrientedVolume::SetMesh( const SMDS_Mesh* theMesh )
+{
+  myMesh = theMesh;
+}
+
+bool BadOrientedVolume::IsSatisfy( long theId )
+{
+  if ( myMesh == 0 )
+    return false;
+
+  SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
+  return !vTool.IsForward();
+}
+
+SMDSAbs_ElementType BadOrientedVolume::GetType() const
+{
+  return SMDSAbs_Volume;
+}
+
+
+
 /*
   Class       : FreeBorders
   Description : Predicate for free borders
@@ -1175,7 +1475,7 @@ FreeBorders::FreeBorders()
   myMesh = 0;
 }
 
-void FreeBorders::SetMesh( SMDS_Mesh* theMesh )
+void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
 {
   myMesh = theMesh;
 }
@@ -1200,7 +1500,7 @@ FreeEdges::FreeEdges()
   myMesh = 0;
 }
 
-void FreeEdges::SetMesh( SMDS_Mesh* theMesh )
+void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
 {
   myMesh = theMesh;
 }
@@ -1218,7 +1518,7 @@ bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId
       {
         int anId = anElem->GetID();
 
-        if ( i == 0 ) 
+        if ( i == 0 )
           aMap.Add( anId );
         else if ( aMap.Contains( anId ) && anId != theFaceId )
           return false;
@@ -1237,29 +1537,33 @@ bool FreeEdges::IsSatisfy( long theId )
   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
     return false;
 
-  int nbNodes = aFace->NbNodes();
-  const SMDS_MeshNode* aNodes[ nbNodes ];
-  int i = 0;
-  SMDS_ElemIteratorPtr anIter = aFace->nodesIterator();
+  SMDS_ElemIteratorPtr anIter;
+  if ( aFace->IsQuadratic() ) {
+    anIter = static_cast<const SMDS_QuadraticFaceOfNodes*>
+      (aFace)->interlacedNodesElemIterator();
+  }
+  else {
+    anIter = aFace->nodesIterator();
+  }
   if ( anIter != 0 )
+    return false;
+
+  int i = 0, nbNodes = aFace->NbNodes();
+  vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
+  while( anIter->more() )
   {
-    while( anIter->more() )
-    {
-      const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
-      if ( aNode == 0 )
-        return false;
-      aNodes[ i++ ] = aNode;
-    }
+    const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
+    if ( aNode == 0 )
+      return false;
+    aNodes[ i++ ] = aNode;
   }
+  aNodes[ nbNodes ] = aNodes[ 0 ];
 
-  for ( int i = 0; i < nbNodes - 1; i++ )
+  for ( i = 0; i < nbNodes; i++ )
     if ( IsFreeEdge( &aNodes[ i ], theId ) )
       return true;
 
-  aNodes[ 1 ] = aNodes[ nbNodes - 1 ];
-  
-  return IsFreeEdge( &aNodes[ 0 ], theId );
-
+  return false;
 }
 
 SMDSAbs_ElementType FreeEdges::GetType() const
@@ -1284,7 +1588,7 @@ bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
 }
 
 inline void UpdateBorders(const FreeEdges::Border& theBorder,
-                         FreeEdges::TBorders& theRegistry, 
+                         FreeEdges::TBorders& theRegistry,
                          FreeEdges::TBorders& theContainer)
 {
   if(theRegistry.find(theBorder) == theRegistry.end()){
@@ -1302,13 +1606,18 @@ void FreeEdges::GetBoreders(TBorders& theBorders)
   for(; anIter->more(); ){
     const SMDS_MeshFace* anElem = anIter->next();
     long anElemId = anElem->GetID();
-    SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
+    SMDS_ElemIteratorPtr aNodesIter;
+    if ( anElem->IsQuadratic() )
+      aNodesIter = static_cast<const SMDS_QuadraticFaceOfNodes*>(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();
@@ -1345,9 +1654,9 @@ RangeOfIds::RangeOfIds()
 
 //=======================================================================
 // name    : SetMesh
-// Purpose : Set mesh 
+// Purpose : Set mesh
 //=======================================================================
-void RangeOfIds::SetMesh( SMDS_Mesh* theMesh )
+void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh )
 {
   myMesh = theMesh;
 }
@@ -1391,9 +1700,9 @@ void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
     TCollection_AsciiString aStr;
     if ( aMinId != IntegerFirst() )
       aStr += aMinId;
-      
+
     aStr += "-";
-      
+
     if ( aMaxId != IntegerLast() )
       aStr += aMaxId;
 
@@ -1461,7 +1770,7 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
   {
     tmpStr = aStr.Token( ",", i++ );
     int aPos = tmpStr.Search( '-' );
-    
+
     if ( aPos == -1 )
     {
       if ( tmpStr.IsIntegerValue() )
@@ -1473,14 +1782,14 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
     {
       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
       TCollection_AsciiString aMinStr = tmpStr;
-      
+
       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
 
       if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() ||
            !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() )
         return false;
-           
+
       myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
       myMax.Append( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
     }
@@ -1527,7 +1836,7 @@ bool RangeOfIds::IsSatisfy( long theId )
     if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All )
       return false;
   }
-    
+
   if ( myIds.Contains( theId ) )
     return true;
 
@@ -1549,7 +1858,7 @@ Comparator::Comparator():
 Comparator::~Comparator()
 {}
 
-void Comparator::SetMesh( SMDS_Mesh* theMesh )
+void Comparator::SetMesh( const SMDS_Mesh* theMesh )
 {
   if ( myFunctor )
     myFunctor->SetMesh( theMesh );
@@ -1634,7 +1943,7 @@ bool LogicalNOT::IsSatisfy( long theId )
   return myPredicate && !myPredicate->IsSatisfy( theId );
 }
 
-void LogicalNOT::SetMesh( SMDS_Mesh* theMesh )
+void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
 {
   if ( myPredicate )
     myPredicate->SetMesh( theMesh );
@@ -1661,7 +1970,7 @@ LogicalBinary::LogicalBinary()
 LogicalBinary::~LogicalBinary()
 {}
 
-void LogicalBinary::SetMesh( SMDS_Mesh* theMesh )
+void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
 {
   if ( myPredicate1 )
     myPredicate1->SetMesh( theMesh );
@@ -1698,10 +2007,10 @@ SMDSAbs_ElementType LogicalBinary::GetType() const
 */
 bool LogicalAND::IsSatisfy( long theId )
 {
-  return 
-    myPredicate1 && 
-    myPredicate2 && 
-    myPredicate1->IsSatisfy( theId ) && 
+  return
+    myPredicate1 &&
+    myPredicate2 &&
+    myPredicate1->IsSatisfy( theId ) &&
     myPredicate2->IsSatisfy( theId );
 }
 
@@ -1712,10 +2021,10 @@ bool LogicalAND::IsSatisfy( long theId )
 */
 bool LogicalOR::IsSatisfy( long theId )
 {
-  return 
-    myPredicate1 && 
-    myPredicate2 && 
-    myPredicate1->IsSatisfy( theId ) || 
+  return
+    myPredicate1 &&
+    myPredicate2 &&
+    myPredicate1->IsSatisfy( theId ) ||
     myPredicate2->IsSatisfy( theId );
 }
 
@@ -1735,11 +2044,10 @@ void Filter::SetPredicate( PredicatePtr thePredicate )
   myPredicate = thePredicate;
 }
 
-
-template<class TElement, class TIterator, class TPredicate> 
-void FillSequence(const TIterator& theIterator,
-                 TPredicate& thePredicate,
-                 Filter::TIdSequence& theSequence)
+template<class TElement, class TIterator, class TPredicate>
+inline void FillSequence(const TIterator& theIterator,
+                        TPredicate& thePredicate,
+                        Filter::TIdSequence& theSequence)
 {
   if ( theIterator ) {
     while( theIterator->more() ) {
@@ -1751,40 +2059,46 @@ void FillSequence(const TIterator& theIterator,
   }
 }
 
-Filter::TIdSequence
-Filter::GetElementsId( SMDS_Mesh* theMesh )
+void
+Filter::
+GetElementsId( const SMDS_Mesh* theMesh,
+              PredicatePtr thePredicate,
+              TIdSequence& theSequence )
 {
-  TIdSequence aSequence;
-  if ( !theMesh || !myPredicate ) return aSequence;
+  theSequence.clear();
+
+  if ( !theMesh || !thePredicate )
+    return;
 
-  myPredicate->SetMesh( theMesh );
+  thePredicate->SetMesh( theMesh );
 
-  SMDSAbs_ElementType aType = myPredicate->GetType();
+  SMDSAbs_ElementType aType = thePredicate->GetType();
   switch(aType){
-  case SMDSAbs_Node:{
-    FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),myPredicate,aSequence);
+  case SMDSAbs_Node:
+    FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),thePredicate,theSequence);
     break;
-  }
-  case SMDSAbs_Edge:{
-    FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
+  case SMDSAbs_Edge:
+    FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
     break;
-  }
-  case SMDSAbs_Face:{
-    FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
+  case SMDSAbs_Face:
+    FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
     break;
-  }
-  case SMDSAbs_Volume:{
-    FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
+  case SMDSAbs_Volume:
+    FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),thePredicate,theSequence);
     break;
-  }
-  case SMDSAbs_All:{
-    FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
-    FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
-    FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
+  case SMDSAbs_All:
+    FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
+    FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
+    FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),thePredicate,theSequence);
     break;
   }
-  }
-  return aSequence;
+}
+
+void
+Filter::GetElementsId( const SMDS_Mesh* theMesh,
+                      Filter::TIdSequence& theSequence )
+{
+  GetElementsId(theMesh,myPredicate,theSequence);
 }
 
 /*
@@ -1793,9 +2107,9 @@ Filter::GetElementsId( SMDS_Mesh* theMesh )
 
 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
 
-/*  
+/*
    Internal class Link
-*/ 
+*/
 
 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
                           SMDS_MeshNode* theNode2 )
@@ -1832,7 +2146,7 @@ bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
 
 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
                             const ManifoldPart::Link& theLink2 )
-{ 
+{
   return theLink1.IsEqual( theLink2 );
 }
 
@@ -1848,7 +2162,7 @@ ManifoldPart::~ManifoldPart()
   myMesh = 0;
 }
 
-void ManifoldPart::SetMesh( SMDS_Mesh* theMesh )
+void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
 {
   myMesh = theMesh;
   process();
@@ -1878,7 +2192,7 @@ bool ManifoldPart::process()
 {
   myMapIds.Clear();
   myMapBadGeomIds.Clear();
-  
+
   myAllFacePtr.clear();
   myAllFacePtrIntDMap.clear();
   if ( !myMesh )
@@ -1910,7 +2224,7 @@ bool ManifoldPart::process()
     if ( fi == aStartIndx )
       isStartTreat = true;
     // as result next time when fi will be equal to aStartIndx
-    
+
     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
     if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
       continue;
@@ -1943,7 +2257,7 @@ static void getLinks( const SMDS_MeshFace* theFace,
   SMDS_MeshNode* aNode = 0;
   for ( ; aNodeItr->more() && i <= aNbNode; )
   {
-    
+
     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
     if ( i == 1 )
       aNode = aN1;
@@ -1962,17 +2276,15 @@ static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
   TColgp_Array1OfXYZ anArrOfXYZ(1,4);
   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
   int i = 1;
-  for ( ; aNodeItr->more() && i <= 4; i++ )
-  {
+  for ( ; aNodeItr->more() && i <= 4; i++ ) {
     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
     anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
   }
-  
+
   gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1);
   gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1);
   n  = q1 ^ q2;
-  if ( aNbNode > 3 )
-  {
+  if ( aNbNode > 3 ) {
     gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
     n += q2 ^ q3;
   }
@@ -1992,7 +2304,7 @@ bool ManifoldPart::findConnected
   theResFaces.Clear();
   if ( !theAllFacePtrInt.size() )
     return false;
-  
+
   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
   {
     myMapBadGeomIds.Add( theStartFace->GetID() );
@@ -2004,7 +2316,7 @@ bool ManifoldPart::findConnected
   theResFaces.Add( theStartFace->GetID() );
   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
 
-  expandBoundary( aMapOfBoundary, aSeqOfBoundary, 
+  expandBoundary( aMapOfBoundary, aSeqOfBoundary,
                  aDMapLinkFace, theNonManifold, theStartFace );
 
   bool isDone = false;
@@ -2022,7 +2334,7 @@ bool ManifoldPart::findConnected
 
       ManifoldPart::TVectorOfFacePtr aFaces;
       // find next
-      if ( myIsOnlyManifold && 
+      if ( myIsOnlyManifold &&
            (theNonManifold.find( aLink ) != theNonManifold.end()) )
         continue;
       else
@@ -2046,7 +2358,7 @@ bool ManifoldPart::findConnected
           continue;
         }
       }
-      
+
       // compare normal with normals of neighbor element
       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
@@ -2065,7 +2377,7 @@ bool ManifoldPart::findConnected
           continue;
         // add new element to connected and extend the boundaries.
         theResFaces.Add( anNextFaceID );
-        expandBoundary( aMapOfBoundary, aSeqOfBoundary, 
+        expandBoundary( aMapOfBoundary, aSeqOfBoundary,
                         aDMapLinkFace, theNonManifold, aNextFace );
         isToReset = true;
       }
@@ -2101,7 +2413,7 @@ void ManifoldPart::expandBoundary
 {
   ManifoldPart::TVectorOfLink aLinks;
   getLinks( theNextFace, aLinks );
-  int aNbLink = aLinks.size();
+  int aNbLink = (int)aLinks.size();
   for ( int i = 0; i < aNbLink; i++ )
   {
     ManifoldPart::Link aLink = aLinks[ i ];
@@ -2177,8 +2489,8 @@ ElementsOnSurface::~ElementsOnSurface()
   myMesh = 0;
 }
 
-void ElementsOnSurface::SetMesh( SMDS_Mesh* theMesh )
-{ 
+void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
+{
   if ( myMesh == theMesh )
     return;
   myMesh = theMesh;