Salome HOME
Merge from V5_1_main branch 24/11/2010
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
index e6946a84aa003b27ed9ebea83326e7f271102131..5238b982f3c55ee99f83bf4048d2096b2f866e3c 100644 (file)
@@ -1,4 +1,4 @@
-//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
+//  Copyright (C) 2007-2010  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
 //
 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+
 #include "SMESH_ControlsDef.hxx"
 
 #include <set>
+#include <limits>
 
 #include <BRepAdaptor_Surface.hxx>
 #include <BRepClass_FaceClassifier.hxx>
 */
 
 namespace{
+
+  inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
+  {
+    return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() );
+  }
+
   inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
   {
     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
@@ -149,27 +157,47 @@ namespace{
 //     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++;
-//       }
-//     }
+//      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;
   }
 
+  gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 )
+  {
+    int aNbNode = theFace->NbNodes();
+
+    gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0));
+    gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0));
+    gp_XYZ n  = q1 ^ q2;
+    if ( aNbNode > 3 ) {
+      gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0));
+      n += q2 ^ q3;
+    }
+    double len = n.Modulus();
+    bool zeroLen = ( len <= numeric_limits<double>::min());
+    if ( !zeroLen )
+      n /= len;
+
+    if (ok) *ok = !zeroLen;
+
+    return n;
+  }
 }
 
 
@@ -177,8 +205,8 @@ namespace{
 using namespace SMESH::Controls;
 
 /*
-                                FUNCTORS
-*/
*                               FUNCTORS
+ */
 
 /*
   Class       : NumericalFunctor
@@ -276,6 +304,63 @@ double NumericalFunctor::GetValue( long theId )
   return 0.;
 }
 
+//================================================================================
+/*!
+ * \brief Return histogram of functor values
+ *  \param nbIntervals - number of intervals
+ *  \param nbEvents - number of mesh elements having values within i-th interval
+ *  \param funValues - boundaries of intervals
+ */
+//================================================================================
+
+void NumericalFunctor::GetHistogram(int                  nbIntervals,
+                                    std::vector<int>&    nbEvents,
+                                    std::vector<double>& funValues)
+{
+  if ( nbIntervals < 1 ||
+       !myMesh ||
+       !myMesh->GetMeshInfo().NbElements( GetType() ))
+    return;
+  nbEvents.resize( nbIntervals, 0 );
+  funValues.resize( nbIntervals+1 );
+
+  // get all values sorted
+  std::multiset< double > values;
+  SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType());
+  while ( elemIt->more() )
+    values.insert( GetValue( elemIt->next()->GetID() ));
+
+  // case nbIntervals == 1
+  funValues[0] = *values.begin();
+  funValues[nbIntervals] = *values.rbegin();
+  if ( nbIntervals == 1 )
+  {
+    nbEvents[0] = values.size();
+    return;
+  }
+  // case of 1 value
+  if (funValues.front() == funValues.back())
+  {
+    nbEvents.resize( 1 );
+    nbEvents[0] = values.size();
+    funValues[1] = funValues.back();
+    funValues.resize( 2 );
+  }
+  // generic case
+  std::multiset< double >::iterator min = values.begin(), max;
+  for ( int i = 0; i < nbIntervals; ++i )
+  {
+    double r = (i+1) / double( nbIntervals );
+    funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
+    if ( min != values.end() && *min <= funValues[i+1] )
+    {
+      max = values.upper_bound( funValues[i+1] ); // greater than funValues[i+1], or end()
+      nbEvents[i] = std::distance( min, max );
+      min = max;
+    }
+  }
+}
+
 //=======================================================================
 //function : GetValue
 //purpose  : 
@@ -312,6 +397,246 @@ SMDSAbs_ElementType Volume::GetType() const
 }
 
 
+/*
+  Class       : MaxElementLength2D
+  Description : Functor calculating maximum length of 2D element
+*/
+
+double MaxElementLength2D::GetValue( long theElementId )
+{
+  TSequenceOfXYZ P;
+  if( GetPoints( theElementId, P ) ) {
+    double aVal = 0;
+    const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
+    SMDSAbs_ElementType aType = aElem->GetType();
+    int len = P.size();
+    switch( aType ) {
+    case SMDSAbs_Face:
+      if( len == 3 ) { // triangles
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 1 ));
+        aVal = Max(L1,Max(L2,L3));
+        break;
+      }
+      else if( len == 4 ) { // quadrangles
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 4 ));
+        double L4 = getDistance(P( 4 ),P( 1 ));
+        double D1 = getDistance(P( 1 ),P( 3 ));
+        double D2 = getDistance(P( 2 ),P( 4 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
+        break;
+      }
+      else if( len == 6 ) { // quadratic triangles
+        double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+        double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+        double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
+        aVal = Max(L1,Max(L2,L3));
+        break;
+      }
+      else if( len == 8 ) { // quadratic quadrangles
+        double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+        double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+        double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
+        double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
+        double D1 = getDistance(P( 1 ),P( 5 ));
+        double D2 = getDistance(P( 3 ),P( 7 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
+        break;
+      }
+    }
+
+    if( myPrecision >= 0 )
+    {
+      double prec = pow( 10., (double)myPrecision );
+      aVal = floor( aVal * prec + 0.5 ) / prec;
+    }
+    return aVal;
+  }
+  return 0.;
+}
+
+double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+  return Value;
+}
+
+SMDSAbs_ElementType MaxElementLength2D::GetType() const
+{
+  return SMDSAbs_Face;
+}
+
+/*
+  Class       : MaxElementLength3D
+  Description : Functor calculating maximum length of 3D element
+*/
+
+double MaxElementLength3D::GetValue( long theElementId )
+{
+  TSequenceOfXYZ P;
+  if( GetPoints( theElementId, P ) ) {
+    double aVal = 0;
+    const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
+    SMDSAbs_ElementType aType = aElem->GetType();
+    int len = P.size();
+    switch( aType ) {
+    case SMDSAbs_Volume:
+      if( len == 4 ) { // tetras
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 1 ));
+        double L4 = getDistance(P( 1 ),P( 4 ));
+        double L5 = getDistance(P( 2 ),P( 4 ));
+        double L6 = getDistance(P( 3 ),P( 4 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        break;
+      }
+      else if( len == 5 ) { // pyramids
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 4 ));
+        double L4 = getDistance(P( 4 ),P( 1 ));
+        double L5 = getDistance(P( 1 ),P( 5 ));
+        double L6 = getDistance(P( 2 ),P( 5 ));
+        double L7 = getDistance(P( 3 ),P( 5 ));
+        double L8 = getDistance(P( 4 ),P( 5 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(L7,L8));
+        break;
+      }
+      else if( len == 6 ) { // pentas
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 1 ));
+        double L4 = getDistance(P( 4 ),P( 5 ));
+        double L5 = getDistance(P( 5 ),P( 6 ));
+        double L6 = getDistance(P( 6 ),P( 4 ));
+        double L7 = getDistance(P( 1 ),P( 4 ));
+        double L8 = getDistance(P( 2 ),P( 5 ));
+        double L9 = getDistance(P( 3 ),P( 6 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(Max(L7,L8),L9));
+        break;
+      }
+      else if( len == 8 ) { // hexas
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 4 ));
+        double L4 = getDistance(P( 4 ),P( 1 ));
+        double L5 = getDistance(P( 5 ),P( 6 ));
+        double L6 = getDistance(P( 6 ),P( 7 ));
+        double L7 = getDistance(P( 7 ),P( 8 ));
+        double L8 = getDistance(P( 8 ),P( 5 ));
+        double L9 = getDistance(P( 1 ),P( 5 ));
+        double L10= getDistance(P( 2 ),P( 6 ));
+        double L11= getDistance(P( 3 ),P( 7 ));
+        double L12= getDistance(P( 4 ),P( 8 ));
+        double D1 = getDistance(P( 1 ),P( 7 ));
+        double D2 = getDistance(P( 2 ),P( 8 ));
+        double D3 = getDistance(P( 3 ),P( 5 ));
+        double D4 = getDistance(P( 4 ),P( 6 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+        aVal = Max(aVal,Max(L11,L12));
+        aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
+        break;
+      }
+      else if( len == 10 ) { // quadratic tetras
+        double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
+        double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+        double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
+        double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        break;
+      }
+      else if( len == 13 ) { // quadratic pyramids
+        double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+        double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+        double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+        double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
+        double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
+        double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(L7,L8));
+        break;
+      }
+      else if( len == 15 ) { // quadratic pentas
+        double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+        double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+        double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
+        double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
+        double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
+        double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
+        double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(Max(L7,L8),L9));
+        break;
+      }
+      else if( len == 20 ) { // quadratic hexas
+        double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
+        double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
+        double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
+        double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
+        double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
+        double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
+        double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
+        double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
+        double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
+        double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
+        double D1 = getDistance(P( 1 ),P( 7 ));
+        double D2 = getDistance(P( 2 ),P( 8 ));
+        double D3 = getDistance(P( 3 ),P( 5 ));
+        double D4 = getDistance(P( 4 ),P( 6 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+        aVal = Max(aVal,Max(L11,L12));
+        aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
+        break;
+      }
+      else if( len > 1 && aElem->IsPoly() ) { // polys
+        // get the maximum distance between all pairs of nodes
+        for( int i = 1; i <= len; i++ ) {
+          for( int j = 1; j <= len; j++ ) {
+            if( j > i ) { // optimization of the loop
+              double D = getDistance( P(i), P(j) );
+              aVal = Max( aVal, D );
+            }
+          }
+        }
+      }
+    }
+
+    if( myPrecision >= 0 )
+    {
+      double prec = pow( 10., (double)myPrecision );
+      aVal = floor( aVal * prec + 0.5 ) / prec;
+    }
+    return aVal;
+  }
+  return 0.;
+}
+
+double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+  return Value;
+}
+
+SMDSAbs_ElementType MaxElementLength3D::GetType() const
+{
+  return SMDSAbs_Volume;
+}
+
+
 /*
   Class       : MinimumAngle
   Description : Functor for calculation of minimum angle
@@ -407,47 +732,94 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
     return alfa * maxLen * half_perimeter / anArea;
   }
   else if( nbNodes == 4 ) { // quadrangle
-    // return aspect ratio of the worst triange which can be built
+    // Compute lengths of the sides
+    std::vector< 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);
+    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
-    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;
-  }
-  else { // nbNodes==8 - quadratic quadrangle
-    // return aspect ratio of the worst triange which can be built
+    std::vector< 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) );
+    anArea[3] = getArea( P(2), P(3), P(4) );
+    // Q = alpha * L * C1 / C2, where
+    //
+    // alpha = sqrt( 1/32 )
+    // L = max( L1, L2, L3, L4, D1, D2 )
+    // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+    // 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. );
+    double C2 = Min( anArea[ 0 ],
+                  Min( anArea[ 1 ],
+                    Min( anArea[ 2 ], anArea[ 3 ] ) ) );
+    if ( C2 <= Precision::Confusion() )
+      return 0.;
+    return alpha * L * C1 / C2;
+  }
+  else if( nbNodes == 8 ){ // nbNodes==8 - quadratic quadrangle
+    // Compute lengths of the sides
+    std::vector< 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);
+    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
-    TSequenceOfXYZ triaPnts(3);
-    // triangle on nodes 1 3 2
-    triaPnts(1) = P(1);
-    triaPnts(2) = P(5);
-    triaPnts(3) = P(3);
-    double ar = GetValue( triaPnts );
-    // triangle on nodes 1 3 4
-    triaPnts(3) = P(7);
-    ar = Max ( ar, GetValue( triaPnts ));
-    // triangle on nodes 1 2 4
-    triaPnts(2) = P(3);
-    ar = Max ( ar, GetValue( triaPnts ));
-    // triangle on nodes 3 2 4
-    triaPnts(1) = P(5);
-    ar = Max ( ar, GetValue( triaPnts ));
-
-    return ar;
+    std::vector< 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) );
+    anArea[3] = getArea( P(3), P(5), P(7) );
+    // Q = alpha * L * C1 / C2, where
+    //
+    // alpha = sqrt( 1/32 )
+    // L = max( L1, L2, L3, L4, D1, D2 )
+    // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+    // 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. );
+    double C2 = Min( anArea[ 0 ],
+                  Min( anArea[ 1 ],
+                    Min( anArea[ 2 ], anArea[ 3 ] ) ) );
+    if ( C2 <= Precision::Confusion() )
+      return 0.;
+    return alpha * L * C1 / C2;
   }
+  return 0;
 }
 
 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
@@ -476,9 +848,9 @@ namespace{
 
   inline double getArea(double theHalfPerim, double theTria[3]){
     return sqrt(theHalfPerim*
-               (theHalfPerim-theTria[0])*
-               (theHalfPerim-theTria[1])*
-               (theHalfPerim-theTria[2]));
+                (theHalfPerim-theTria[0])*
+                (theHalfPerim-theTria[1])*
+                (theHalfPerim-theTria[2]));
   }
 
   inline double getVolume(double theLen[6]){
@@ -963,16 +1335,20 @@ SMDSAbs_ElementType Skew::GetType() const
 */
 double Area::GetValue( const TSequenceOfXYZ& P )
 {
-  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);
+  double val = 0.0;
+  if ( P.size() > 2 ) {
+    gp_Vec aVec1( P(2) - P(1) );
+    gp_Vec aVec2( P(3) - P(1) );
+    gp_Vec SumVec = aVec1 ^ aVec2;
+    for (int i=4; i<=P.size(); i++) {
+      gp_Vec aVec1( P(i-1) - P(1) );
+      gp_Vec aVec2( P(i) - P(1) );
+      gp_Vec tmp = aVec1 ^ aVec2;
+      SumVec.Add(tmp);
+    }
+    val = SumVec.Magnitude() * 0.5;
   }
-  return SumVec.Magnitude() * 0.5;
+  return val;
 }
 
 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
@@ -1036,160 +1412,160 @@ double Length2D::GetValue( long theElementId)
     case SMDSAbs_Node:
     case SMDSAbs_Edge:
       if (len == 2){
-       aVal = getDistance( P( 1 ), P( 2 ) );
+        aVal = getDistance( P( 1 ), P( 2 ) );
         break;
       }
       else if (len == 3){ // quadratic edge
-       aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
+        aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
         break;
       }
     case SMDSAbs_Face:
       if (len == 3){ // triangles
-       double L1 = getDistance(P( 1 ),P( 2 ));
-       double L2 = getDistance(P( 2 ),P( 3 ));
-       double L3 = getDistance(P( 3 ),P( 1 ));
-       aVal = Max(L1,Max(L2,L3));
-       break;
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 1 ));
+        aVal = Max(L1,Max(L2,L3));
+        break;
       }
       else if (len == 4){ // quadrangles
-       double L1 = getDistance(P( 1 ),P( 2 ));
-       double L2 = getDistance(P( 2 ),P( 3 ));
-       double L3 = getDistance(P( 3 ),P( 4 ));
-       double L4 = getDistance(P( 4 ),P( 1 ));
-       aVal = Max(Max(L1,L2),Max(L3,L4));
-       break;
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 4 ));
+        double L4 = getDistance(P( 4 ),P( 1 ));
+        aVal = Max(Max(L1,L2),Max(L3,L4));
+        break;
       }
       if (len == 6){ // quadratic triangles
-       double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
-       double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
-       double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
-       aVal = Max(L1,Max(L2,L3));
+        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;
+        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;
+        double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+        double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+        double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
+        double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
+        aVal = Max(Max(L1,L2),Max(L3,L4));
+        break;
       }
     case SMDSAbs_Volume:
       if (len == 4){ // tetraidrs
-       double L1 = getDistance(P( 1 ),P( 2 ));
-       double L2 = getDistance(P( 2 ),P( 3 ));
-       double L3 = getDistance(P( 3 ),P( 1 ));
-       double L4 = getDistance(P( 1 ),P( 4 ));
-       double L5 = getDistance(P( 2 ),P( 4 ));
-       double L6 = getDistance(P( 3 ),P( 4 ));
-       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-       break;
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 1 ));
+        double L4 = getDistance(P( 1 ),P( 4 ));
+        double L5 = getDistance(P( 2 ),P( 4 ));
+        double L6 = getDistance(P( 3 ),P( 4 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        break;
       }
       else if (len == 5){ // piramids
-       double L1 = getDistance(P( 1 ),P( 2 ));
-       double L2 = getDistance(P( 2 ),P( 3 ));
-       double L3 = getDistance(P( 3 ),P( 1 ));
-       double L4 = getDistance(P( 4 ),P( 1 ));
-       double L5 = getDistance(P( 1 ),P( 5 ));
-       double L6 = getDistance(P( 2 ),P( 5 ));
-       double L7 = getDistance(P( 3 ),P( 5 ));
-       double L8 = getDistance(P( 4 ),P( 5 ));
-
-       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-       aVal = Max(aVal,Max(L7,L8));
-       break;
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 4 ));
+        double L4 = getDistance(P( 4 ),P( 1 ));
+        double L5 = getDistance(P( 1 ),P( 5 ));
+        double L6 = getDistance(P( 2 ),P( 5 ));
+        double L7 = getDistance(P( 3 ),P( 5 ));
+        double L8 = getDistance(P( 4 ),P( 5 ));
+
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(L7,L8));
+        break;
       }
       else if (len == 6){ // pentaidres
-       double L1 = getDistance(P( 1 ),P( 2 ));
-       double L2 = getDistance(P( 2 ),P( 3 ));
-       double L3 = getDistance(P( 3 ),P( 1 ));
-       double L4 = getDistance(P( 4 ),P( 5 ));
-       double L5 = getDistance(P( 5 ),P( 6 ));
-       double L6 = getDistance(P( 6 ),P( 4 ));
-       double L7 = getDistance(P( 1 ),P( 4 ));
-       double L8 = getDistance(P( 2 ),P( 5 ));
-       double L9 = getDistance(P( 3 ),P( 6 ));
-
-       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-       aVal = Max(aVal,Max(Max(L7,L8),L9));
-       break;
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 1 ));
+        double L4 = getDistance(P( 4 ),P( 5 ));
+        double L5 = getDistance(P( 5 ),P( 6 ));
+        double L6 = getDistance(P( 6 ),P( 4 ));
+        double L7 = getDistance(P( 1 ),P( 4 ));
+        double L8 = getDistance(P( 2 ),P( 5 ));
+        double L9 = getDistance(P( 3 ),P( 6 ));
+
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(Max(L7,L8),L9));
+        break;
       }
       else if (len == 8){ // hexaider
-       double L1 = getDistance(P( 1 ),P( 2 ));
-       double L2 = getDistance(P( 2 ),P( 3 ));
-       double L3 = getDistance(P( 3 ),P( 4 ));
-       double L4 = getDistance(P( 4 ),P( 1 ));
-       double L5 = getDistance(P( 5 ),P( 6 ));
-       double L6 = getDistance(P( 6 ),P( 7 ));
-       double L7 = getDistance(P( 7 ),P( 8 ));
-       double L8 = getDistance(P( 8 ),P( 5 ));
-       double L9 = getDistance(P( 1 ),P( 5 ));
-       double L10= getDistance(P( 2 ),P( 6 ));
-       double L11= getDistance(P( 3 ),P( 7 ));
-       double L12= getDistance(P( 4 ),P( 8 ));
-
-       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-       aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
-       aVal = Max(aVal,Max(L11,L12));
-       break;
+        double L1 = getDistance(P( 1 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 4 ));
+        double L4 = getDistance(P( 4 ),P( 1 ));
+        double L5 = getDistance(P( 5 ),P( 6 ));
+        double L6 = getDistance(P( 6 ),P( 7 ));
+        double L7 = getDistance(P( 7 ),P( 8 ));
+        double L8 = getDistance(P( 8 ),P( 5 ));
+        double L9 = getDistance(P( 1 ),P( 5 ));
+        double L10= getDistance(P( 2 ),P( 6 ));
+        double L11= getDistance(P( 3 ),P( 7 ));
+        double L12= getDistance(P( 4 ),P( 8 ));
+
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+        aVal = Max(aVal,Max(L11,L12));
+        break;
 
       }
 
       if (len == 10){ // quadratic tetraidrs
-       double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
-       double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
-       double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
-       double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
-       double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
-       double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
-       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-       break;
+        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;
+        double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+        double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+        double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+        double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
+        double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
+        double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(L7,L8));
+        break;
       }
       else if (len == 15){ // quadratic pentaidres
-       double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
-       double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
-       double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
-       double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
-       double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
-       double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
-       double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
-       double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
-       double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
-       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-       aVal = Max(aVal,Max(Max(L7,L8),L9));
-       break;
+        double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
+        double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
+        double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+        double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+        double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
+        double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
+        double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
+        double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
+        double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
+        aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+        aVal = Max(aVal,Max(Max(L7,L8),L9));
+        break;
       }
       else if (len == 20){ // quadratic hexaider
-       double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
-       double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
-       double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
-       double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
-       double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
-       double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
-       double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
-       double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
-       double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
-       double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
-       double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
-       double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
-       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
-       aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
-       aVal = Max(aVal,Max(L11,L12));
-       break;
+        double 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;
 
       }
 
@@ -1467,12 +1843,12 @@ void MultiConnection2D::GetValues(MValues& theValues){
       Value aValue(aNodeId[1],aNodeId[2]);
       MValues::iterator aItr = theValues.find(aValue);
       if (aItr != theValues.end()){
-       aItr->second += 1;
-       //aNbConnects = nb;
+        aItr->second += 1;
+        //aNbConnects = nb;
       }
       else {
-       theValues[aValue] = 1;
-       //aNbConnects = 1;
+        theValues[aValue] = 1;
+        //aNbConnects = 1;
       }
       //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
       aNodeId[1] = aNodeId[2];
@@ -1651,8 +2027,8 @@ bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
 }
 
 inline void UpdateBorders(const FreeEdges::Border& theBorder,
-                         FreeEdges::TBorders& theRegistry,
-                         FreeEdges::TBorders& theContainer)
+                          FreeEdges::TBorders& theRegistry,
+                          FreeEdges::TBorders& theContainer)
 {
   if(theRegistry.find(theBorder) == theRegistry.end()){
     theRegistry.insert(theBorder);
@@ -1938,8 +2314,10 @@ bool ElemGeomType::IsSatisfy( long theId )
 {
   if (!myMesh) return false;
   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
+  if ( !anElem )
+    return false;
   const SMDSAbs_ElementType anElemType = anElem->GetType();
-  if ( !anElem || (myType != SMDSAbs_All && anElemType != myType) )
+  if ( myType != SMDSAbs_All && anElemType != myType )
     return false;
   const int aNbNode = anElem->NbNodes();
   bool isOk = false;
@@ -1955,22 +2333,22 @@ bool ElemGeomType::IsSatisfy( long theId )
 
   case SMDSAbs_Face:
     if ( myGeomType == SMDSGeom_TRIANGLE )
-      isOk = (!anElem->IsPoly() && aNbNode == 3);
+      isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 6 : aNbNode == 3));
     else if ( myGeomType == SMDSGeom_QUADRANGLE )
-      isOk = (!anElem->IsPoly() && aNbNode == 4);
+      isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 8 : aNbNode == 4));
     else if ( myGeomType == SMDSGeom_POLYGON )
       isOk = anElem->IsPoly();
     break;
 
   case SMDSAbs_Volume:
     if ( myGeomType == SMDSGeom_TETRA )
-      isOk = (!anElem->IsPoly() && aNbNode == 4);
+      isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 10 : aNbNode == 4));
     else if ( myGeomType == SMDSGeom_PYRAMID )
-      isOk = (!anElem->IsPoly() && aNbNode == 5);
+      isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 13 : aNbNode == 5));
     else if ( myGeomType == SMDSGeom_PENTA )
-      isOk = (!anElem->IsPoly() && aNbNode == 6);
+      isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 15 : aNbNode == 6));
     else if ( myGeomType == SMDSGeom_HEXA )
-      isOk = (!anElem->IsPoly() && aNbNode == 8);
+      isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 20 : aNbNode == 8));
      else if ( myGeomType == SMDSGeom_POLYHEDRA )
       isOk = anElem->IsPoly();
     break;
@@ -1999,13 +2377,71 @@ SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
   return myGeomType;
 }
 
+//================================================================================
+/*!
+ * \brief Class CoplanarFaces
+ */
+//================================================================================
+
+CoplanarFaces::CoplanarFaces()
+  : myMesh(0), myFaceID(0), myToler(0)
+{
+}
+bool CoplanarFaces::IsSatisfy( long theElementId )
+{
+  if ( myCoplanarIDs.empty() )
+  {
+    // Build a set of coplanar face ids
+
+    if ( !myMesh || !myFaceID || !myToler )
+      return false;
+
+    const SMDS_MeshElement* face = myMesh->FindElement( myFaceID );
+    if ( !face || face->GetType() != SMDSAbs_Face )
+      return false;
+
+    bool normOK;
+    gp_Vec myNorm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
+    if (!normOK)
+      return false;
+
+    const double radianTol = myToler * PI180;
+    typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TFaceIt;
+    std::set<const SMDS_MeshElement*> checkedFaces, checkedNodes;
+    std::list<const SMDS_MeshElement*> faceQueue( 1, face );
+    while ( !faceQueue.empty() )
+    {
+      face = faceQueue.front();
+      if ( checkedFaces.insert( face ).second )
+      {
+        gp_Vec norm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
+        if (!normOK || myNorm.Angle( norm ) <= radianTol)
+        {
+          myCoplanarIDs.insert( face->GetID() );
+          std::set<const SMDS_MeshElement*> neighborFaces;
+          for ( int i = 0; i < face->NbCornerNodes(); ++i )
+          {
+            const SMDS_MeshNode* n = face->GetNode( i );
+            if ( checkedNodes.insert( n ).second )
+              neighborFaces.insert( TFaceIt( n->GetInverseElementIterator(SMDSAbs_Face)),
+                                    TFaceIt());
+          }
+          faceQueue.insert( faceQueue.end(), neighborFaces.begin(), neighborFaces.end() );
+        }
+      }
+      faceQueue.pop_front();
+    }
+  }
+  return myCoplanarIDs.count( theElementId );
+}
+
 /*
-  Class       : RangeOfIds
-  Description : Predicate for Range of Ids.
-                Range may be specified with two ways.
-                1. Using AddToRange method
-                2. With SetRangeStr method. Parameter of this method is a string
-                   like as "1,2,3,50-60,63,67,70-"
+  *Class       : RangeOfIds
+  *Description : Predicate for Range of Ids.
+  *              Range may be specified with two ways.
+  *              1. Using AddToRange method
+  *              2. With SetRangeStr method. Parameter of this method is a string
+  *                 like as "1,2,3,50-60,63,67,70-"
 */
 
 //=======================================================================
@@ -2412,15 +2848,15 @@ void Filter::SetPredicate( PredicatePtr thePredicate )
 
 template<class TElement, class TIterator, class TPredicate>
 inline void FillSequence(const TIterator& theIterator,
-                        TPredicate& thePredicate,
-                        Filter::TIdSequence& theSequence)
+                         TPredicate& thePredicate,
+                         Filter::TIdSequence& theSequence)
 {
   if ( theIterator ) {
     while( theIterator->more() ) {
       TElement anElem = theIterator->next();
       long anId = anElem->GetID();
       if ( thePredicate->IsSatisfy( anId ) )
-       theSequence.push_back( anId );
+        theSequence.push_back( anId );
     }
   }
 }
@@ -2428,8 +2864,8 @@ inline void FillSequence(const TIterator& theIterator,
 void
 Filter::
 GetElementsId( const SMDS_Mesh* theMesh,
-              PredicatePtr thePredicate,
-              TIdSequence& theSequence )
+               PredicatePtr thePredicate,
+               TIdSequence& theSequence )
 {
   theSequence.clear();
 
@@ -2462,7 +2898,7 @@ GetElementsId( const SMDS_Mesh* theMesh,
 
 void
 Filter::GetElementsId( const SMDS_Mesh* theMesh,
-                      Filter::TIdSequence& theSequence )
+                       Filter::TIdSequence& theSequence )
 {
   GetElementsId(theMesh,myPredicate,theSequence);
 }
@@ -2635,32 +3071,6 @@ static void getLinks( const SMDS_MeshFace* theFace,
   }
 }
 
-static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
-{
-  gp_XYZ n;
-  int aNbNode = theFace->NbNodes();
-  TColgp_Array1OfXYZ anArrOfXYZ(1,4);
-  SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
-  int i = 1;
-  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 ) {
-    gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
-    n += q2 ^ q3;
-  }
-  double len = n.Modulus();
-  if ( len > 0 )
-    n /= len;
-
-  return n;
-}
-
 bool ManifoldPart::findConnected
                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
                   SMDS_MeshFace*                           theStartFace,
@@ -3254,31 +3664,54 @@ void ElementsOnShape::process (const SMDS_MeshElement* theElemPtr)
 TSequenceOfXYZ::TSequenceOfXYZ()
 {}
 
-TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : std::vector<gp_XYZ>(n)
+TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n)
 {}
 
-TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const value_type& t) : std::vector<gp_XYZ>(n,t)
+TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t)
 {}
 
-TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : std::vector<gp_XYZ>(theSequenceOfXYZ)
+TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray)
 {}
 
 template <class InputIterator>
-TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): std::vector<gp_XYZ>(theBegin,theEnd)
+TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd)
+{}
+
+TSequenceOfXYZ::~TSequenceOfXYZ()
 {}
 
 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
 {
-  std::vector<gp_XYZ>::operator=(theSequenceOfXYZ);
+  myArray = theSequenceOfXYZ.myArray;
   return *this;
 }
 
-std::vector<gp_XYZ>::reference TSequenceOfXYZ::operator()(size_type n)
+gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
+{
+  return myArray[n-1];
+}
+
+const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
+{
+  return myArray[n-1];
+}
+
+void TSequenceOfXYZ::clear()
+{
+  myArray.clear();
+}
+
+void TSequenceOfXYZ::reserve(size_type n)
+{
+  myArray.reserve(n);
+}
+
+void TSequenceOfXYZ::push_back(const gp_XYZ& v)
 {
-  return std::vector<gp_XYZ>::operator[](n-1);
+  myArray.push_back(v);
 }
 
-std::vector<gp_XYZ>::const_reference TSequenceOfXYZ::operator()(size_type n) const
+TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
 {
-  return std::vector<gp_XYZ>::operator[](n-1);
+  return myArray.size();
 }