src/Controls/SMESH_Controls.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21
  22 #include "SMESH_ControlsDef.hxx"
  23
  24 #include "SMDS_BallElement.hxx"
  25 #include "SMDS_Iterator.hxx"
  26 #include "SMDS_Mesh.hxx"
  27 #include "SMDS_MeshElement.hxx"
  28 #include "SMDS_MeshNode.hxx"
  29 #include "SMDS_QuadraticEdge.hxx"
  30 #include "SMDS_QuadraticFaceOfNodes.hxx"
  31 #include "SMDS_VolumeTool.hxx"
  32 #include "SMESHDS_GroupBase.hxx"
  33 #include "SMESHDS_Mesh.hxx"
  34 #include "SMESH_OctreeNode.hxx"
  35
  36 #include <BRepAdaptor_Surface.hxx>
  37 #include <BRepClass_FaceClassifier.hxx>
  38 #include <BRep_Tool.hxx>
  39 #include <Geom_CylindricalSurface.hxx>
  40 #include <Geom_Plane.hxx>
  41 #include <Geom_Surface.hxx>
  42 #include <Precision.hxx>
  43 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
  44 #include <TColStd_MapOfInteger.hxx>
  45 #include <TColStd_SequenceOfAsciiString.hxx>
  46 #include <TColgp_Array1OfXYZ.hxx>
  47 #include <TopAbs.hxx>
  48 #include <TopoDS.hxx>
  49 #include <TopoDS_Edge.hxx>
  50 #include <TopoDS_Face.hxx>
  51 #include <TopoDS_Iterator.hxx>
  52 #include <TopoDS_Shape.hxx>
  53 #include <TopoDS_Vertex.hxx>
  54 #include <gp_Ax3.hxx>
  55 #include <gp_Cylinder.hxx>
  56 #include <gp_Dir.hxx>
  57 #include <gp_Pln.hxx>
  58 #include <gp_Pnt.hxx>
  59 #include <gp_Vec.hxx>
  60 #include <gp_XYZ.hxx>
  61
  62 #include <vtkMeshQuality.h>
  63
  64 #include <set>
  65 #include <limits>
  66
  67
  68 /*
  69                             AUXILIARY METHODS
  70 */
  71
  72 namespace {
  73
  74   const double theEps = 1e-100;
  75   const double theInf = 1e+100;
  76
  77   inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
  78   {
  79     return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() );
  80   }
  81
  82   inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
  83   {
  84     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
  85
  86     return v1.Magnitude() < gp::Resolution() ||
  87       v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
  88   }
  89
  90   inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
  91   {
  92     gp_Vec aVec1( P2 - P1 );
  93     gp_Vec aVec2( P3 - P1 );
  94     return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
  95   }
  96
  97   inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
  98   {
  99     return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
 100   }
 101
 102
 103
 104   inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
 105   {
 106     double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
 107     return aDist;
 108   }
 109
 110   int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId )
 111   {
 112     if ( theMesh == 0 )
 113       return 0;
 114
 115     const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
 116     if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */)
 117       return 0;
 118
 119     // for each pair of nodes in anEdge (there are 2 pairs in a quadratic edge)
 120     // count elements containing both nodes of the pair.
 121     // Note that there may be such cases for a quadratic edge (a horizontal line):
 122     //
 123     //  Case 1          Case 2
 124     //  |     |      |        |      |
 125     //  |     |      |        |      |
 126     //  +-----+------+  +-----+------+
 127     //  |            |  |            |
 128     //  |            |  |            |
 129     // result sould be 2 in both cases
 130     //
 131     int aResult0 = 0, aResult1 = 0;
 132      // last node, it is a medium one in a quadratic edge
 133     const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
 134     const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 );
 135     const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 );
 136     if ( aNode1 == aLastNode ) aNode1 = 0;
 137
 138     SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator();
 139     while( anElemIter->more() ) {
 140       const SMDS_MeshElement* anElem = anElemIter->next();
 141       if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
 142         SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
 143         while ( anIter->more() ) {
 144           if ( const SMDS_MeshElement* anElemNode = anIter->next() ) {
 145             if ( anElemNode == aNode0 ) {
 146               aResult0++;
 147               if ( !aNode1 ) break; // not a quadratic edge
 148             }
 149             else if ( anElemNode == aNode1 )
 150               aResult1++;
 151           }
 152         }
 153       }
 154     }
 155     int aResult = std::max ( aResult0, aResult1 );
 156
 157 //     TColStd_MapOfInteger aMap;
 158
 159 //     SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
 160 //     if ( anIter != 0 ) {
 161 //       while( anIter->more() ) {
 162 //      const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
 163 //      if ( aNode == 0 )
 164 //        return 0;
 165 //      SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
 166 //      while( anElemIter->more() ) {
 167 //        const SMDS_MeshElement* anElem = anElemIter->next();
 168 //        if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
 169 //          int anId = anElem->GetID();
 170
 171 //          if ( anIter->more() )              // i.e. first node
 172 //            aMap.Add( anId );
 173 //          else if ( aMap.Contains( anId ) )
 174 //            aResult++;
 175 //        }
 176 //      }
 177 //       }
 178 //     }
 179
 180     return aResult;
 181   }
 182
 183   gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 )
 184   {
 185     int aNbNode = theFace->NbNodes();
 186
 187     gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0));
 188     gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0));
 189     gp_XYZ n  = q1 ^ q2;
 190     if ( aNbNode > 3 ) {
 191       gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0));
 192       n += q2 ^ q3;
 193     }
 194     double len = n.Modulus();
 195     bool zeroLen = ( len <= numeric_limits<double>::min());
 196     if ( !zeroLen )
 197       n /= len;
 198
 199     if (ok) *ok = !zeroLen;
 200
 201     return n;
 202   }
 203 }
 204
 205
 206
 207 using namespace SMESH::Controls;
 208
 209 /*
 210  *                               FUNCTORS
 211  */
 212
 213 //================================================================================
 214 /*
 215   Class       : NumericalFunctor
 216   Description : Base class for numerical functors
 217 */
 218 //================================================================================
 219
 220 NumericalFunctor::NumericalFunctor():
 221   myMesh(NULL)
 222 {
 223   myPrecision = -1;
 224 }
 225
 226 void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh )
 227 {
 228   myMesh = theMesh;
 229 }
 230
 231 bool NumericalFunctor::GetPoints(const int theId,
 232                                  TSequenceOfXYZ& theRes ) const
 233 {
 234   theRes.clear();
 235
 236   if ( myMesh == 0 )
 237     return false;
 238
 239   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
 240   if ( !anElem || anElem->GetType() != this->GetType() )
 241     return false;
 242
 243   return GetPoints( anElem, theRes );
 244 }
 245
 246 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
 247                                  TSequenceOfXYZ&         theRes )
 248 {
 249   theRes.clear();
 250
 251   if ( anElem == 0 )
 252     return false;
 253
 254   theRes.reserve( anElem->NbNodes() );
 255
 256   // Get nodes of the element
 257   SMDS_ElemIteratorPtr anIter;
 258
 259   if ( anElem->IsQuadratic() ) {
 260     switch ( anElem->GetType() ) {
 261     case SMDSAbs_Edge:
 262       anIter = dynamic_cast<const SMDS_VtkEdge*>
 263         (anElem)->interlacedNodesElemIterator();
 264       break;
 265     case SMDSAbs_Face:
 266       anIter = dynamic_cast<const SMDS_VtkFace*>
 267         (anElem)->interlacedNodesElemIterator();
 268       break;
 269     default:
 270       anIter = anElem->nodesIterator();
 271       //return false;
 272     }
 273   }
 274   else {
 275     anIter = anElem->nodesIterator();
 276   }
 277
 278   if ( anIter ) {
 279     while( anIter->more() ) {
 280       if ( const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( anIter->next() ))
 281         theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
 282     }
 283   }
 284
 285   return true;
 286 }
 287
 288 long  NumericalFunctor::GetPrecision() const
 289 {
 290   return myPrecision;
 291 }
 292
 293 void  NumericalFunctor::SetPrecision( const long thePrecision )
 294 {
 295   myPrecision = thePrecision;
 296   myPrecisionValue = pow( 10., (double)( myPrecision ) );
 297 }
 298
 299 double NumericalFunctor::GetValue( long theId )
 300 {
 301   double aVal = 0;
 302
 303   myCurrElement = myMesh->FindElement( theId );
 304
 305   TSequenceOfXYZ P;
 306   if ( GetPoints( theId, P ))
 307     aVal = Round( GetValue( P ));
 308
 309   return aVal;
 310 }
 311
 312 double NumericalFunctor::Round( const double & aVal )
 313 {
 314   return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal;
 315 }
 316
 317 //================================================================================
 318 /*!
 319  * \brief Return histogram of functor values
 320  *  \param nbIntervals - number of intervals
 321  *  \param nbEvents - number of mesh elements having values within i-th interval
 322  *  \param funValues - boundaries of intervals
 323  *  \param elements - elements to check vulue of; empty list means "of all"
 324  *  \param minmax - boundaries of diapason of values to divide into intervals
 325  */
 326 //================================================================================
 327
 328 void NumericalFunctor::GetHistogram(int                  nbIntervals,
 329                                     std::vector<int>&    nbEvents,
 330                                     std::vector<double>& funValues,
 331                                     const vector<int>&   elements,
 332                                     const double*        minmax,
 333                                     const bool           isLogarithmic)
 334 {
 335   if ( nbIntervals < 1 ||
 336        !myMesh ||
 337        !myMesh->GetMeshInfo().NbElements( GetType() ))
 338     return;
 339   nbEvents.resize( nbIntervals, 0 );
 340   funValues.resize( nbIntervals+1 );
 341
 342   // get all values sorted
 343   std::multiset< double > values;
 344   if ( elements.empty() )
 345   {
 346     SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType());
 347     while ( elemIt->more() )
 348       values.insert( GetValue( elemIt->next()->GetID() ));
 349   }
 350   else
 351   {
 352     vector<int>::const_iterator id = elements.begin();
 353     for ( ; id != elements.end(); ++id )
 354       values.insert( GetValue( *id ));
 355   }
 356
 357   if ( minmax )
 358   {
 359     funValues[0] = minmax[0];
 360     funValues[nbIntervals] = minmax[1];
 361   }
 362   else
 363   {
 364     funValues[0] = *values.begin();
 365     funValues[nbIntervals] = *values.rbegin();
 366   }
 367   // case nbIntervals == 1
 368   if ( nbIntervals == 1 )
 369   {
 370     nbEvents[0] = values.size();
 371     return;
 372   }
 373   // case of 1 value
 374   if (funValues.front() == funValues.back())
 375   {
 376     nbEvents.resize( 1 );
 377     nbEvents[0] = values.size();
 378     funValues[1] = funValues.back();
 379     funValues.resize( 2 );
 380   }
 381   // generic case
 382   std::multiset< double >::iterator min = values.begin(), max;
 383   for ( int i = 0; i < nbIntervals; ++i )
 384   {
 385     // find end value of i-th interval
 386     double r = (i+1) / double(nbIntervals);
 387     if (isLogarithmic && funValues.front() > 1e-07 && funValues.back() > 1e-07) {
 388       double logmin = log10(funValues.front());
 389       double lval = logmin + r * (log10(funValues.back()) - logmin);
 390       funValues[i+1] = pow(10.0, lval);
 391     }
 392     else {
 393       funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
 394     }
 395
 396     // count values in the i-th interval if there are any
 397     if ( min != values.end() && *min <= funValues[i+1] )
 398     {
 399       // find the first value out of the interval
 400       max = values.upper_bound( funValues[i+1] ); // max is greater than funValues[i+1], or end()
 401       nbEvents[i] = std::distance( min, max );
 402       min = max;
 403     }
 404   }
 405   // add values larger than minmax[1]
 406   nbEvents.back() += std::distance( min, values.end() );
 407 }
 408
 409 //=======================================================================
 410 /*
 411   Class       : Volume
 412   Description : Functor calculating volume of a 3D element
 413 */
 414 //================================================================================
 415
 416 double Volume::GetValue( long theElementId )
 417 {
 418   if ( theElementId && myMesh ) {
 419     SMDS_VolumeTool aVolumeTool;
 420     if ( aVolumeTool.Set( myMesh->FindElement( theElementId )))
 421       return aVolumeTool.GetSize();
 422   }
 423   return 0;
 424 }
 425
 426 double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
 427 {
 428   return Value;
 429 }
 430
 431 SMDSAbs_ElementType Volume::GetType() const
 432 {
 433   return SMDSAbs_Volume;
 434 }
 435
 436 //=======================================================================
 437 /*
 438   Class       : MaxElementLength2D
 439   Description : Functor calculating maximum length of 2D element
 440 */
 441 //================================================================================
 442
 443 double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P )
 444 {
 445   if(P.size() == 0)
 446     return 0.;
 447   double aVal = 0;
 448   int len = P.size();
 449   if( len == 3 ) { // triangles
 450     double L1 = getDistance(P( 1 ),P( 2 ));
 451     double L2 = getDistance(P( 2 ),P( 3 ));
 452     double L3 = getDistance(P( 3 ),P( 1 ));
 453     aVal = Max(L1,Max(L2,L3));
 454   }
 455   else if( len == 4 ) { // quadrangles
 456     double L1 = getDistance(P( 1 ),P( 2 ));
 457     double L2 = getDistance(P( 2 ),P( 3 ));
 458     double L3 = getDistance(P( 3 ),P( 4 ));
 459     double L4 = getDistance(P( 4 ),P( 1 ));
 460     double D1 = getDistance(P( 1 ),P( 3 ));
 461     double D2 = getDistance(P( 2 ),P( 4 ));
 462     aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
 463   }
 464   else if( len == 6 ) { // quadratic triangles
 465     double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
 466     double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
 467     double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
 468     aVal = Max(L1,Max(L2,L3));
 469   }
 470   else if( len == 8 || len == 9 ) { // quadratic quadrangles
 471     double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
 472     double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
 473     double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
 474     double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
 475     double D1 = getDistance(P( 1 ),P( 5 ));
 476     double D2 = getDistance(P( 3 ),P( 7 ));
 477     aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
 478   }
 479
 480   if( myPrecision >= 0 )
 481   {
 482     double prec = pow( 10., (double)myPrecision );
 483     aVal = floor( aVal * prec + 0.5 ) / prec;
 484   }
 485   return aVal;
 486 }
 487
 488 double MaxElementLength2D::GetValue( long theElementId )
 489 {
 490   TSequenceOfXYZ P;
 491   return GetPoints( theElementId, P ) ? GetValue(P) : 0.0;
 492 }
 493
 494 double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const
 495 {
 496   return Value;
 497 }
 498
 499 SMDSAbs_ElementType MaxElementLength2D::GetType() const
 500 {
 501   return SMDSAbs_Face;
 502 }
 503
 504 //=======================================================================
 505 /*
 506   Class       : MaxElementLength3D
 507   Description : Functor calculating maximum length of 3D element
 508 */
 509 //================================================================================
 510
 511 double MaxElementLength3D::GetValue( long theElementId )
 512 {
 513   TSequenceOfXYZ P;
 514   if( GetPoints( theElementId, P ) ) {
 515     double aVal = 0;
 516     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
 517     SMDSAbs_ElementType aType = aElem->GetType();
 518     int len = P.size();
 519     switch( aType ) {
 520     case SMDSAbs_Volume:
 521       if( len == 4 ) { // tetras
 522         double L1 = getDistance(P( 1 ),P( 2 ));
 523         double L2 = getDistance(P( 2 ),P( 3 ));
 524         double L3 = getDistance(P( 3 ),P( 1 ));
 525         double L4 = getDistance(P( 1 ),P( 4 ));
 526         double L5 = getDistance(P( 2 ),P( 4 ));
 527         double L6 = getDistance(P( 3 ),P( 4 ));
 528         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 529         break;
 530       }
 531       else if( len == 5 ) { // pyramids
 532         double L1 = getDistance(P( 1 ),P( 2 ));
 533         double L2 = getDistance(P( 2 ),P( 3 ));
 534         double L3 = getDistance(P( 3 ),P( 4 ));
 535         double L4 = getDistance(P( 4 ),P( 1 ));
 536         double L5 = getDistance(P( 1 ),P( 5 ));
 537         double L6 = getDistance(P( 2 ),P( 5 ));
 538         double L7 = getDistance(P( 3 ),P( 5 ));
 539         double L8 = getDistance(P( 4 ),P( 5 ));
 540         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 541         aVal = Max(aVal,Max(L7,L8));
 542         break;
 543       }
 544       else if( len == 6 ) { // pentas
 545         double L1 = getDistance(P( 1 ),P( 2 ));
 546         double L2 = getDistance(P( 2 ),P( 3 ));
 547         double L3 = getDistance(P( 3 ),P( 1 ));
 548         double L4 = getDistance(P( 4 ),P( 5 ));
 549         double L5 = getDistance(P( 5 ),P( 6 ));
 550         double L6 = getDistance(P( 6 ),P( 4 ));
 551         double L7 = getDistance(P( 1 ),P( 4 ));
 552         double L8 = getDistance(P( 2 ),P( 5 ));
 553         double L9 = getDistance(P( 3 ),P( 6 ));
 554         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 555         aVal = Max(aVal,Max(Max(L7,L8),L9));
 556         break;
 557       }
 558       else if( len == 8 ) { // hexas
 559         double L1 = getDistance(P( 1 ),P( 2 ));
 560         double L2 = getDistance(P( 2 ),P( 3 ));
 561         double L3 = getDistance(P( 3 ),P( 4 ));
 562         double L4 = getDistance(P( 4 ),P( 1 ));
 563         double L5 = getDistance(P( 5 ),P( 6 ));
 564         double L6 = getDistance(P( 6 ),P( 7 ));
 565         double L7 = getDistance(P( 7 ),P( 8 ));
 566         double L8 = getDistance(P( 8 ),P( 5 ));
 567         double L9 = getDistance(P( 1 ),P( 5 ));
 568         double L10= getDistance(P( 2 ),P( 6 ));
 569         double L11= getDistance(P( 3 ),P( 7 ));
 570         double L12= getDistance(P( 4 ),P( 8 ));
 571         double D1 = getDistance(P( 1 ),P( 7 ));
 572         double D2 = getDistance(P( 2 ),P( 8 ));
 573         double D3 = getDistance(P( 3 ),P( 5 ));
 574         double D4 = getDistance(P( 4 ),P( 6 ));
 575         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 576         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
 577         aVal = Max(aVal,Max(L11,L12));
 578         aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
 579         break;
 580       }
 581       else if( len == 12 ) { // hexagonal prism
 582         for ( int i1 = 1; i1 < 12; ++i1 )
 583           for ( int i2 = i1+1; i1 <= 12; ++i1 )
 584             aVal = Max( aVal, getDistance(P( i1 ),P( i2 )));
 585         break;
 586       }
 587       else if( len == 10 ) { // quadratic tetras
 588         double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
 589         double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
 590         double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
 591         double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
 592         double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
 593         double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
 594         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 595         break;
 596       }
 597       else if( len == 13 ) { // quadratic pyramids
 598         double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
 599         double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
 600         double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
 601         double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
 602         double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
 603         double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
 604         double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
 605         double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
 606         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 607         aVal = Max(aVal,Max(L7,L8));
 608         break;
 609       }
 610       else if( len == 15 ) { // quadratic pentas
 611         double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
 612         double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
 613         double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
 614         double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
 615         double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
 616         double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
 617         double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
 618         double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
 619         double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
 620         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 621         aVal = Max(aVal,Max(Max(L7,L8),L9));
 622         break;
 623       }
 624       else if( len == 20 || len == 27 ) { // quadratic hexas
 625         double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
 626         double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
 627         double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
 628         double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
 629         double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
 630         double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
 631         double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
 632         double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
 633         double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
 634         double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
 635         double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
 636         double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
 637         double D1 = getDistance(P( 1 ),P( 7 ));
 638         double D2 = getDistance(P( 2 ),P( 8 ));
 639         double D3 = getDistance(P( 3 ),P( 5 ));
 640         double D4 = getDistance(P( 4 ),P( 6 ));
 641         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 642         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
 643         aVal = Max(aVal,Max(L11,L12));
 644         aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
 645         break;
 646       }
 647       else if( len > 1 && aElem->IsPoly() ) { // polys
 648         // get the maximum distance between all pairs of nodes
 649         for( int i = 1; i <= len; i++ ) {
 650           for( int j = 1; j <= len; j++ ) {
 651             if( j > i ) { // optimization of the loop
 652               double D = getDistance( P(i), P(j) );
 653               aVal = Max( aVal, D );
 654             }
 655           }
 656         }
 657       }
 658     }
 659
 660     if( myPrecision >= 0 )
 661     {
 662       double prec = pow( 10., (double)myPrecision );
 663       aVal = floor( aVal * prec + 0.5 ) / prec;
 664     }
 665     return aVal;
 666   }
 667   return 0.;
 668 }
 669
 670 double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const
 671 {
 672   return Value;
 673 }
 674
 675 SMDSAbs_ElementType MaxElementLength3D::GetType() const
 676 {
 677   return SMDSAbs_Volume;
 678 }
 679
 680 //=======================================================================
 681 /*
 682   Class       : MinimumAngle
 683   Description : Functor for calculation of minimum angle
 684 */
 685 //================================================================================
 686
 687 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
 688 {
 689   double aMin;
 690
 691   if (P.size() <3)
 692     return 0.;
 693
 694   aMin = getAngle(P( P.size() ), P( 1 ), P( 2 ));
 695   aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 )));
 696
 697   for (int i=2; i<P.size();i++){
 698       double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) );
 699     aMin = Min(aMin,A0);
 700   }
 701
 702   return aMin * 180.0 / M_PI;
 703 }
 704
 705 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
 706 {
 707   //const double aBestAngle = PI / nbNodes;
 708   const double aBestAngle = 180.0 - ( 360.0 / double(nbNodes) );
 709   return ( fabs( aBestAngle - Value ));
 710 }
 711
 712 SMDSAbs_ElementType MinimumAngle::GetType() const
 713 {
 714   return SMDSAbs_Face;
 715 }
 716
 717
 718 //================================================================================
 719 /*
 720   Class       : AspectRatio
 721   Description : Functor for calculating aspect ratio
 722 */
 723 //================================================================================
 724
 725 double AspectRatio::GetValue( long theId )
 726 {
 727   double aVal = 0;
 728   myCurrElement = myMesh->FindElement( theId );
 729   if ( myCurrElement && myCurrElement->GetVtkType() == VTK_QUAD )
 730   {
 731     // issue 21723
 732     vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
 733     if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
 734       aVal = Round( vtkMeshQuality::QuadAspectRatio( avtkCell ));
 735   }
 736   else
 737   {
 738     TSequenceOfXYZ P;
 739     if ( GetPoints( myCurrElement, P ))
 740       aVal = Round( GetValue( P ));
 741   }
 742   return aVal;
 743 }
 744
 745 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
 746 {
 747   // According to "Mesh quality control" by Nadir Bouhamau referring to
 748   // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis.
 749   // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4
 750   // PAL10872
 751
 752   int nbNodes = P.size();
 753
 754   if ( nbNodes < 3 )
 755     return 0;
 756
 757   // Compute aspect ratio
 758
 759   if ( nbNodes == 3 ) {
 760     // Compute lengths of the sides
 761     std::vector< double > aLen (nbNodes);
 762     for ( int i = 0; i < nbNodes - 1; i++ )
 763       aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
 764     aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
 765     // Q = alfa * h * p / S, where
 766     //
 767     // alfa = sqrt( 3 ) / 6
 768     // h - length of the longest edge
 769     // p - half perimeter
 770     // S - triangle surface
 771     const double alfa = sqrt( 3. ) / 6.;
 772     double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
 773     double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
 774     double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
 775     if ( anArea <= theEps  )
 776       return theInf;
 777     return alfa * maxLen * half_perimeter / anArea;
 778   }
 779   else if ( nbNodes == 6 ) { // quadratic triangles
 780     // Compute lengths of the sides
 781     std::vector< double > aLen (3);
 782     aLen[0] = getDistance( P(1), P(3) );
 783     aLen[1] = getDistance( P(3), P(5) );
 784     aLen[2] = getDistance( P(5), P(1) );
 785     // Q = alfa * h * p / S, where
 786     //
 787     // alfa = sqrt( 3 ) / 6
 788     // h - length of the longest edge
 789     // p - half perimeter
 790     // S - triangle surface
 791     const double alfa = sqrt( 3. ) / 6.;
 792     double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
 793     double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
 794     double anArea = getArea( P(1), P(3), P(5) );
 795     if ( anArea <= theEps )
 796       return theInf;
 797     return alfa * maxLen * half_perimeter / anArea;
 798   }
 799   else if( nbNodes == 4 ) { // quadrangle
 800     // Compute lengths of the sides
 801     std::vector< double > aLen (4);
 802     aLen[0] = getDistance( P(1), P(2) );
 803     aLen[1] = getDistance( P(2), P(3) );
 804     aLen[2] = getDistance( P(3), P(4) );
 805     aLen[3] = getDistance( P(4), P(1) );
 806     // Compute lengths of the diagonals
 807     std::vector< double > aDia (2);
 808     aDia[0] = getDistance( P(1), P(3) );
 809     aDia[1] = getDistance( P(2), P(4) );
 810     // Compute areas of all triangles which can be built
 811     // taking three nodes of the quadrangle
 812     std::vector< double > anArea (4);
 813     anArea[0] = getArea( P(1), P(2), P(3) );
 814     anArea[1] = getArea( P(1), P(2), P(4) );
 815     anArea[2] = getArea( P(1), P(3), P(4) );
 816     anArea[3] = getArea( P(2), P(3), P(4) );
 817     // Q = alpha * L * C1 / C2, where
 818     //
 819     // alpha = sqrt( 1/32 )
 820     // L = max( L1, L2, L3, L4, D1, D2 )
 821     // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
 822     // C2 = min( S1, S2, S3, S4 )
 823     // Li - lengths of the edges
 824     // Di - lengths of the diagonals
 825     // Si - areas of the triangles
 826     const double alpha = sqrt( 1 / 32. );
 827     double L = Max( aLen[ 0 ],
 828                  Max( aLen[ 1 ],
 829                    Max( aLen[ 2 ],
 830                      Max( aLen[ 3 ],
 831                        Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
 832     double C1 = sqrt( ( aLen[0] * aLen[0] +
 833                         aLen[1] * aLen[1] +
 834                         aLen[2] * aLen[2] +
 835                         aLen[3] * aLen[3] ) / 4. );
 836     double C2 = Min( anArea[ 0 ],
 837                   Min( anArea[ 1 ],
 838                     Min( anArea[ 2 ], anArea[ 3 ] ) ) );
 839     if ( C2 <= theEps )
 840       return theInf;
 841     return alpha * L * C1 / C2;
 842   }
 843   else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle
 844     // Compute lengths of the sides
 845     std::vector< double > aLen (4);
 846     aLen[0] = getDistance( P(1), P(3) );
 847     aLen[1] = getDistance( P(3), P(5) );
 848     aLen[2] = getDistance( P(5), P(7) );
 849     aLen[3] = getDistance( P(7), P(1) );
 850     // Compute lengths of the diagonals
 851     std::vector< double > aDia (2);
 852     aDia[0] = getDistance( P(1), P(5) );
 853     aDia[1] = getDistance( P(3), P(7) );
 854     // Compute areas of all triangles which can be built
 855     // taking three nodes of the quadrangle
 856     std::vector< double > anArea (4);
 857     anArea[0] = getArea( P(1), P(3), P(5) );
 858     anArea[1] = getArea( P(1), P(3), P(7) );
 859     anArea[2] = getArea( P(1), P(5), P(7) );
 860     anArea[3] = getArea( P(3), P(5), P(7) );
 861     // Q = alpha * L * C1 / C2, where
 862     //
 863     // alpha = sqrt( 1/32 )
 864     // L = max( L1, L2, L3, L4, D1, D2 )
 865     // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
 866     // C2 = min( S1, S2, S3, S4 )
 867     // Li - lengths of the edges
 868     // Di - lengths of the diagonals
 869     // Si - areas of the triangles
 870     const double alpha = sqrt( 1 / 32. );
 871     double L = Max( aLen[ 0 ],
 872                  Max( aLen[ 1 ],
 873                    Max( aLen[ 2 ],
 874                      Max( aLen[ 3 ],
 875                        Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
 876     double C1 = sqrt( ( aLen[0] * aLen[0] +
 877                         aLen[1] * aLen[1] +
 878                         aLen[2] * aLen[2] +
 879                         aLen[3] * aLen[3] ) / 4. );
 880     double C2 = Min( anArea[ 0 ],
 881                   Min( anArea[ 1 ],
 882                     Min( anArea[ 2 ], anArea[ 3 ] ) ) );
 883     if ( C2 <= theEps )
 884       return theInf;
 885     return alpha * L * C1 / C2;
 886   }
 887   return 0;
 888 }
 889
 890 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
 891 {
 892   // the aspect ratio is in the range [1.0,infinity]
 893   // < 1.0 = very bad, zero area
 894   // 1.0 = good
 895   // infinity = bad
 896   return ( Value < 0.9 ) ? 1000 : Value / 1000.;
 897 }
 898
 899 SMDSAbs_ElementType AspectRatio::GetType() const
 900 {
 901   return SMDSAbs_Face;
 902 }
 903
 904
 905 //================================================================================
 906 /*
 907   Class       : AspectRatio3D
 908   Description : Functor for calculating aspect ratio
 909 */
 910 //================================================================================
 911
 912 namespace{
 913
 914   inline double getHalfPerimeter(double theTria[3]){
 915     return (theTria[0] + theTria[1] + theTria[2])/2.0;
 916   }
 917
 918   inline double getArea(double theHalfPerim, double theTria[3]){
 919     return sqrt(theHalfPerim*
 920                 (theHalfPerim-theTria[0])*
 921                 (theHalfPerim-theTria[1])*
 922                 (theHalfPerim-theTria[2]));
 923   }
 924
 925   inline double getVolume(double theLen[6]){
 926     double a2 = theLen[0]*theLen[0];
 927     double b2 = theLen[1]*theLen[1];
 928     double c2 = theLen[2]*theLen[2];
 929     double d2 = theLen[3]*theLen[3];
 930     double e2 = theLen[4]*theLen[4];
 931     double f2 = theLen[5]*theLen[5];
 932     double P = 4.0*a2*b2*d2;
 933     double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
 934     double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
 935     return sqrt(P-Q+R)/12.0;
 936   }
 937
 938   inline double getVolume2(double theLen[6]){
 939     double a2 = theLen[0]*theLen[0];
 940     double b2 = theLen[1]*theLen[1];
 941     double c2 = theLen[2]*theLen[2];
 942     double d2 = theLen[3]*theLen[3];
 943     double e2 = theLen[4]*theLen[4];
 944     double f2 = theLen[5]*theLen[5];
 945
 946     double P = a2*e2*(b2+c2+d2+f2-a2-e2);
 947     double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
 948     double R = c2*d2*(a2+b2+e2+f2-c2-d2);
 949     double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
 950
 951     return sqrt(P+Q+R-S)/12.0;
 952   }
 953
 954   inline double getVolume(const TSequenceOfXYZ& P){
 955     gp_Vec aVec1( P( 2 ) - P( 1 ) );
 956     gp_Vec aVec2( P( 3 ) - P( 1 ) );
 957     gp_Vec aVec3( P( 4 ) - P( 1 ) );
 958     gp_Vec anAreaVec( aVec1 ^ aVec2 );
 959     return fabs(aVec3 * anAreaVec) / 6.0;
 960   }
 961
 962   inline double getMaxHeight(double theLen[6])
 963   {
 964     double aHeight = std::max(theLen[0],theLen[1]);
 965     aHeight = std::max(aHeight,theLen[2]);
 966     aHeight = std::max(aHeight,theLen[3]);
 967     aHeight = std::max(aHeight,theLen[4]);
 968     aHeight = std::max(aHeight,theLen[5]);
 969     return aHeight;
 970   }
 971
 972 }
 973
 974 double AspectRatio3D::GetValue( long theId )
 975 {
 976   double aVal = 0;
 977   myCurrElement = myMesh->FindElement( theId );
 978   if ( myCurrElement && myCurrElement->GetVtkType() == VTK_TETRA )
 979   {
 980     // Action from CoTech | ACTION 31.3:
 981     // EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with
 982     // those of ParaView. The library used by ParaView for those calculations can be reused in SMESH.
 983     vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
 984     if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
 985       aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
 986   }
 987   else
 988   {
 989     TSequenceOfXYZ P;
 990     if ( GetPoints( myCurrElement, P ))
 991       aVal = Round( GetValue( P ));
 992   }
 993   return aVal;
 994 }
 995
 996 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
 997 {
 998   double aQuality = 0.0;
 999   if(myCurrElement->IsPoly()) return aQuality;
1000
1001   int nbNodes = P.size();
1002
1003   if(myCurrElement->IsQuadratic()) {
1004     if(nbNodes==10) nbNodes=4; // quadratic tetrahedron
1005     else if(nbNodes==13) nbNodes=5; // quadratic pyramid
1006     else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
1007     else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
1008     else if(nbNodes==27) nbNodes=8; // quadratic hexahedron
1009     else return aQuality;
1010   }
1011
1012   switch(nbNodes) {
1013   case 4:{
1014     double aLen[6] = {
1015       getDistance(P( 1 ),P( 2 )), // a
1016       getDistance(P( 2 ),P( 3 )), // b
1017       getDistance(P( 3 ),P( 1 )), // c
1018       getDistance(P( 2 ),P( 4 )), // d
1019       getDistance(P( 3 ),P( 4 )), // e
1020       getDistance(P( 1 ),P( 4 ))  // f
1021     };
1022     double aTria[4][3] = {
1023       {aLen[0],aLen[1],aLen[2]}, // abc
1024       {aLen[0],aLen[3],aLen[5]}, // adf
1025       {aLen[1],aLen[3],aLen[4]}, // bde
1026       {aLen[2],aLen[4],aLen[5]}  // cef
1027     };
1028     double aSumArea = 0.0;
1029     double aHalfPerimeter = getHalfPerimeter(aTria[0]);
1030     double anArea = getArea(aHalfPerimeter,aTria[0]);
1031     aSumArea += anArea;
1032     aHalfPerimeter = getHalfPerimeter(aTria[1]);
1033     anArea = getArea(aHalfPerimeter,aTria[1]);
1034     aSumArea += anArea;
1035     aHalfPerimeter = getHalfPerimeter(aTria[2]);
1036     anArea = getArea(aHalfPerimeter,aTria[2]);
1037     aSumArea += anArea;
1038     aHalfPerimeter = getHalfPerimeter(aTria[3]);
1039     anArea = getArea(aHalfPerimeter,aTria[3]);
1040     aSumArea += anArea;
1041     double aVolume = getVolume(P);
1042     //double aVolume = getVolume(aLen);
1043     double aHeight = getMaxHeight(aLen);
1044     static double aCoeff = sqrt(2.0)/12.0;
1045     if ( aVolume > DBL_MIN )
1046       aQuality = aCoeff*aHeight*aSumArea/aVolume;
1047     break;
1048   }
1049   case 5:{
1050     {
1051       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
1052       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1053     }
1054     {
1055       gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
1056       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1057     }
1058     {
1059       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
1060       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1061     }
1062     {
1063       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
1064       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1065     }
1066     break;
1067   }
1068   case 6:{
1069     {
1070       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
1071       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1072     }
1073     {
1074       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
1075       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1076     }
1077     {
1078       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
1079       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1080     }
1081     {
1082       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1083       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1084     }
1085     {
1086       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
1087       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1088     }
1089     {
1090       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
1091       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1092     }
1093     break;
1094   }
1095   case 8:{
1096     {
1097       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1098       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1099     }
1100     {
1101       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
1102       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1103     }
1104     {
1105       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
1106       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1107     }
1108     {
1109       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
1110       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1111     }
1112     {
1113       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
1114       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1115     }
1116     {
1117       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
1118       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1119     }
1120     {
1121       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
1122       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1123     }
1124     {
1125       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
1126       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1127     }
1128     {
1129       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
1130       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1131     }
1132     {
1133       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
1134       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1135     }
1136     {
1137       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
1138       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1139     }
1140     {
1141       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
1142       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1143     }
1144     {
1145       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
1146       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1147     }
1148     {
1149       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
1150       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1151     }
1152     {
1153       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
1154       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1155     }
1156     {
1157       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
1158       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1159     }
1160     {
1161       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
1162       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1163     }
1164     {
1165       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
1166       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1167     }
1168     {
1169       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
1170       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1171     }
1172     {
1173       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
1174       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1175     }
1176     {
1177       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
1178       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1179     }
1180     {
1181       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1182       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1183     }
1184     {
1185       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
1186       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1187     }
1188     {
1189       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
1190       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1191     }
1192     {
1193       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1194       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1195     }
1196     {
1197       gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
1198       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1199     }
1200     {
1201       gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
1202       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1203     }
1204     {
1205       gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
1206       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1207     }
1208     {
1209       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
1210       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1211     }
1212     {
1213       gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
1214       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1215     }
1216     {
1217       gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
1218       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1219     }
1220     {
1221       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
1222       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1223     }
1224     {
1225       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
1226       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1227     }
1228     break;
1229   }
1230   case 12:
1231     {
1232       gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
1233       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1234     }
1235     {
1236       gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
1237       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1238     }
1239     {
1240       gp_XYZ aXYZ[8] = {P( 3 ),P( 4 ),P( 6 ),P( 1 ),P( 9 ),P( 10 ),P( 12 ),P( 7 )};
1241       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1242     }
1243     break;
1244   } // switch(nbNodes)
1245
1246   if ( nbNodes > 4 ) {
1247     // avaluate aspect ratio of quadranle faces
1248     AspectRatio aspect2D;
1249     SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
1250     int nbFaces = SMDS_VolumeTool::NbFaces( type );
1251     TSequenceOfXYZ points(4);
1252     for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume
1253       if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
1254         continue;
1255       const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
1256       for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face
1257         points( p + 1 ) = P( pInd[ p ] + 1 );
1258       aQuality = std::max( aQuality, aspect2D.GetValue( points ));
1259     }
1260   }
1261   return aQuality;
1262 }
1263
1264 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
1265 {
1266   // the aspect ratio is in the range [1.0,infinity]
1267   // 1.0 = good
1268   // infinity = bad
1269   return Value / 1000.;
1270 }
1271
1272 SMDSAbs_ElementType AspectRatio3D::GetType() const
1273 {
1274   return SMDSAbs_Volume;
1275 }
1276
1277
1278 //================================================================================
1279 /*
1280   Class       : Warping
1281   Description : Functor for calculating warping
1282 */
1283 //================================================================================
1284
1285 double Warping::GetValue( const TSequenceOfXYZ& P )
1286 {
1287   if ( P.size() != 4 )
1288     return 0;
1289
1290   gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.;
1291
1292   double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
1293   double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
1294   double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
1295   double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
1296
1297   double val = Max( Max( A1, A2 ), Max( A3, A4 ) );
1298
1299   const double eps = 0.1; // val is in degrees
1300
1301   return val < eps ? 0. : val;
1302 }
1303
1304 double Warping::ComputeA( const gp_XYZ& thePnt1,
1305                           const gp_XYZ& thePnt2,
1306                           const gp_XYZ& thePnt3,
1307                           const gp_XYZ& theG ) const
1308 {
1309   double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
1310   double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
1311   double L = Min( aLen1, aLen2 ) * 0.5;
1312   if ( L < theEps )
1313     return theInf;
1314
1315   gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
1316   gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
1317   gp_XYZ N  = GI.Crossed( GJ );
1318
1319   if ( N.Modulus() < gp::Resolution() )
1320     return M_PI / 2;
1321
1322   N.Normalize();
1323
1324   double H = ( thePnt2 - theG ).Dot( N );
1325   return asin( fabs( H / L ) ) * 180. / M_PI;
1326 }
1327
1328 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
1329 {
1330   // the warp is in the range [0.0,PI/2]
1331   // 0.0 = good (no warp)
1332   // PI/2 = bad  (face pliee)
1333   return Value;
1334 }
1335
1336 SMDSAbs_ElementType Warping::GetType() const
1337 {
1338   return SMDSAbs_Face;
1339 }
1340
1341
1342 //================================================================================
1343 /*
1344   Class       : Taper
1345   Description : Functor for calculating taper
1346 */
1347 //================================================================================
1348
1349 double Taper::GetValue( const TSequenceOfXYZ& P )
1350 {
1351   if ( P.size() != 4 )
1352     return 0.;
1353
1354   // Compute taper
1355   double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2.;
1356   double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2.;
1357   double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2.;
1358   double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2.;
1359
1360   double JA = 0.25 * ( J1 + J2 + J3 + J4 );
1361   if ( JA <= theEps )
1362     return theInf;
1363
1364   double T1 = fabs( ( J1 - JA ) / JA );
1365   double T2 = fabs( ( J2 - JA ) / JA );
1366   double T3 = fabs( ( J3 - JA ) / JA );
1367   double T4 = fabs( ( J4 - JA ) / JA );
1368
1369   double val = Max( Max( T1, T2 ), Max( T3, T4 ) );
1370
1371   const double eps = 0.01;
1372
1373   return val < eps ? 0. : val;
1374 }
1375
1376 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
1377 {
1378   // the taper is in the range [0.0,1.0]
1379   // 0.0  = good (no taper)
1380   // 1.0 = bad  (les cotes opposes sont allignes)
1381   return Value;
1382 }
1383
1384 SMDSAbs_ElementType Taper::GetType() const
1385 {
1386   return SMDSAbs_Face;
1387 }
1388
1389 //================================================================================
1390 /*
1391   Class       : Skew
1392   Description : Functor for calculating skew in degrees
1393 */
1394 //================================================================================
1395
1396 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
1397 {
1398   gp_XYZ p12 = ( p2 + p1 ) / 2.;
1399   gp_XYZ p23 = ( p3 + p2 ) / 2.;
1400   gp_XYZ p31 = ( p3 + p1 ) / 2.;
1401
1402   gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
1403
1404   return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
1405 }
1406
1407 double Skew::GetValue( const TSequenceOfXYZ& P )
1408 {
1409   if ( P.size() != 3 && P.size() != 4 )
1410     return 0.;
1411
1412   // Compute skew
1413   const double PI2 = M_PI / 2.;
1414   if ( P.size() == 3 )
1415   {
1416     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
1417     double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
1418     double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
1419
1420     return Max( A0, Max( A1, A2 ) ) * 180. / M_PI;
1421   }
1422   else
1423   {
1424     gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.;
1425     gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.;
1426     gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.;
1427     gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.;
1428
1429     gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
1430     double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
1431       ? 0. : fabs( PI2 - v1.Angle( v2 ) );
1432
1433     double val = A * 180. / M_PI;
1434
1435     const double eps = 0.1; // val is in degrees
1436
1437     return val < eps ? 0. : val;
1438   }
1439 }
1440
1441 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
1442 {
1443   // the skew is in the range [0.0,PI/2].
1444   // 0.0 = good
1445   // PI/2 = bad
1446   return Value;
1447 }
1448
1449 SMDSAbs_ElementType Skew::GetType() const
1450 {
1451   return SMDSAbs_Face;
1452 }
1453
1454
1455 //================================================================================
1456 /*
1457   Class       : Area
1458   Description : Functor for calculating area
1459 */
1460 //================================================================================
1461
1462 double Area::GetValue( const TSequenceOfXYZ& P )
1463 {
1464   double val = 0.0;
1465   if ( P.size() > 2 ) {
1466     gp_Vec aVec1( P(2) - P(1) );
1467     gp_Vec aVec2( P(3) - P(1) );
1468     gp_Vec SumVec = aVec1 ^ aVec2;
1469     for (int i=4; i<=P.size(); i++) {
1470       gp_Vec aVec1( P(i-1) - P(1) );
1471       gp_Vec aVec2( P(i) - P(1) );
1472       gp_Vec tmp = aVec1 ^ aVec2;
1473       SumVec.Add(tmp);
1474     }
1475     val = SumVec.Magnitude() * 0.5;
1476   }
1477   return val;
1478 }
1479
1480 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
1481 {
1482   // meaningless as it is not a quality control functor
1483   return Value;
1484 }
1485
1486 SMDSAbs_ElementType Area::GetType() const
1487 {
1488   return SMDSAbs_Face;
1489 }
1490
1491 //================================================================================
1492 /*
1493   Class       : Length
1494   Description : Functor for calculating length of edge
1495 */
1496 //================================================================================
1497
1498 double Length::GetValue( const TSequenceOfXYZ& P )
1499 {
1500   switch ( P.size() ) {
1501   case 2:  return getDistance( P( 1 ), P( 2 ) );
1502   case 3:  return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) );
1503   default: return 0.;
1504   }
1505 }
1506
1507 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
1508 {
1509   // meaningless as it is not quality control functor
1510   return Value;
1511 }
1512
1513 SMDSAbs_ElementType Length::GetType() const
1514 {
1515   return SMDSAbs_Edge;
1516 }
1517
1518 //================================================================================
1519 /*
1520   Class       : Length2D
1521   Description : Functor for calculating length of edge
1522 */
1523 //================================================================================
1524
1525 double Length2D::GetValue( long theElementId)
1526 {
1527   TSequenceOfXYZ P;
1528
1529   //cout<<"Length2D::GetValue"<<endl;
1530   if (GetPoints(theElementId,P)){
1531     //for(int jj=1; jj<=P.size(); jj++)
1532     //  cout<<"jj="<<jj<<" P("<<P(jj).X()<<","<<P(jj).Y()<<","<<P(jj).Z()<<")"<<endl;
1533
1534     double  aVal;// = GetValue( P );
1535     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
1536     SMDSAbs_ElementType aType = aElem->GetType();
1537
1538     int len = P.size();
1539
1540     switch (aType){
1541     case SMDSAbs_All:
1542     case SMDSAbs_Node:
1543     case SMDSAbs_Edge:
1544       if (len == 2){
1545         aVal = getDistance( P( 1 ), P( 2 ) );
1546         break;
1547       }
1548       else if (len == 3){ // quadratic edge
1549         aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
1550         break;
1551       }
1552     case SMDSAbs_Face:
1553       if (len == 3){ // triangles
1554         double L1 = getDistance(P( 1 ),P( 2 ));
1555         double L2 = getDistance(P( 2 ),P( 3 ));
1556         double L3 = getDistance(P( 3 ),P( 1 ));
1557         aVal = Max(L1,Max(L2,L3));
1558         break;
1559       }
1560       else if (len == 4){ // quadrangles
1561         double L1 = getDistance(P( 1 ),P( 2 ));
1562         double L2 = getDistance(P( 2 ),P( 3 ));
1563         double L3 = getDistance(P( 3 ),P( 4 ));
1564         double L4 = getDistance(P( 4 ),P( 1 ));
1565         aVal = Max(Max(L1,L2),Max(L3,L4));
1566         break;
1567       }
1568       if (len == 6){ // quadratic triangles
1569         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1570         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1571         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
1572         aVal = Max(L1,Max(L2,L3));
1573         //cout<<"L1="<<L1<<" L2="<<L2<<"L3="<<L3<<" aVal="<<aVal<<endl;
1574         break;
1575       }
1576       else if (len == 8){ // quadratic quadrangles
1577         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1578         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1579         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
1580         double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
1581         aVal = Max(Max(L1,L2),Max(L3,L4));
1582         break;
1583       }
1584     case SMDSAbs_Volume:
1585       if (len == 4){ // tetraidrs
1586         double L1 = getDistance(P( 1 ),P( 2 ));
1587         double L2 = getDistance(P( 2 ),P( 3 ));
1588         double L3 = getDistance(P( 3 ),P( 1 ));
1589         double L4 = getDistance(P( 1 ),P( 4 ));
1590         double L5 = getDistance(P( 2 ),P( 4 ));
1591         double L6 = getDistance(P( 3 ),P( 4 ));
1592         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1593         break;
1594       }
1595       else if (len == 5){ // piramids
1596         double L1 = getDistance(P( 1 ),P( 2 ));
1597         double L2 = getDistance(P( 2 ),P( 3 ));
1598         double L3 = getDistance(P( 3 ),P( 4 ));
1599         double L4 = getDistance(P( 4 ),P( 1 ));
1600         double L5 = getDistance(P( 1 ),P( 5 ));
1601         double L6 = getDistance(P( 2 ),P( 5 ));
1602         double L7 = getDistance(P( 3 ),P( 5 ));
1603         double L8 = getDistance(P( 4 ),P( 5 ));
1604
1605         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1606         aVal = Max(aVal,Max(L7,L8));
1607         break;
1608       }
1609       else if (len == 6){ // pentaidres
1610         double L1 = getDistance(P( 1 ),P( 2 ));
1611         double L2 = getDistance(P( 2 ),P( 3 ));
1612         double L3 = getDistance(P( 3 ),P( 1 ));
1613         double L4 = getDistance(P( 4 ),P( 5 ));
1614         double L5 = getDistance(P( 5 ),P( 6 ));
1615         double L6 = getDistance(P( 6 ),P( 4 ));
1616         double L7 = getDistance(P( 1 ),P( 4 ));
1617         double L8 = getDistance(P( 2 ),P( 5 ));
1618         double L9 = getDistance(P( 3 ),P( 6 ));
1619
1620         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1621         aVal = Max(aVal,Max(Max(L7,L8),L9));
1622         break;
1623       }
1624       else if (len == 8){ // hexaider
1625         double L1 = getDistance(P( 1 ),P( 2 ));
1626         double L2 = getDistance(P( 2 ),P( 3 ));
1627         double L3 = getDistance(P( 3 ),P( 4 ));
1628         double L4 = getDistance(P( 4 ),P( 1 ));
1629         double L5 = getDistance(P( 5 ),P( 6 ));
1630         double L6 = getDistance(P( 6 ),P( 7 ));
1631         double L7 = getDistance(P( 7 ),P( 8 ));
1632         double L8 = getDistance(P( 8 ),P( 5 ));
1633         double L9 = getDistance(P( 1 ),P( 5 ));
1634         double L10= getDistance(P( 2 ),P( 6 ));
1635         double L11= getDistance(P( 3 ),P( 7 ));
1636         double L12= getDistance(P( 4 ),P( 8 ));
1637
1638         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1639         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
1640         aVal = Max(aVal,Max(L11,L12));
1641         break;
1642
1643       }
1644
1645       if (len == 10){ // quadratic tetraidrs
1646         double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
1647         double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
1648         double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
1649         double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1650         double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
1651         double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
1652         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1653         break;
1654       }
1655       else if (len == 13){ // quadratic piramids
1656         double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
1657         double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
1658         double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1659         double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1660         double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1661         double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
1662         double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
1663         double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
1664         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1665         aVal = Max(aVal,Max(L7,L8));
1666         break;
1667       }
1668       else if (len == 15){ // quadratic pentaidres
1669         double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
1670         double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
1671         double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1672         double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1673         double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
1674         double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
1675         double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
1676         double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
1677         double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
1678         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1679         aVal = Max(aVal,Max(Max(L7,L8),L9));
1680         break;
1681       }
1682       else if (len == 20){ // quadratic hexaider
1683         double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
1684         double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
1685         double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
1686         double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
1687         double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
1688         double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
1689         double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
1690         double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
1691         double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
1692         double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
1693         double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
1694         double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
1695         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1696         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
1697         aVal = Max(aVal,Max(L11,L12));
1698         break;
1699
1700       }
1701
1702     default: aVal=-1;
1703     }
1704
1705     if (aVal <0){
1706       return 0.;
1707     }
1708
1709     if ( myPrecision >= 0 )
1710     {
1711       double prec = pow( 10., (double)( myPrecision ) );
1712       aVal = floor( aVal * prec + 0.5 ) / prec;
1713     }
1714
1715     return aVal;
1716
1717   }
1718   return 0.;
1719 }
1720
1721 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1722 {
1723   // meaningless as it is not quality control functor
1724   return Value;
1725 }
1726
1727 SMDSAbs_ElementType Length2D::GetType() const
1728 {
1729   return SMDSAbs_Face;
1730 }
1731
1732 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
1733   myLength(theLength)
1734 {
1735   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
1736   if(thePntId1 > thePntId2){
1737     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
1738   }
1739 }
1740
1741 bool Length2D::Value::operator<(const Length2D::Value& x) const{
1742   if(myPntId[0] < x.myPntId[0]) return true;
1743   if(myPntId[0] == x.myPntId[0])
1744     if(myPntId[1] < x.myPntId[1]) return true;
1745   return false;
1746 }
1747
1748 void Length2D::GetValues(TValues& theValues){
1749   TValues aValues;
1750   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1751   for(; anIter->more(); ){
1752     const SMDS_MeshFace* anElem = anIter->next();
1753
1754     if(anElem->IsQuadratic()) {
1755       const SMDS_VtkFace* F =
1756         dynamic_cast<const SMDS_VtkFace*>(anElem);
1757       // use special nodes iterator
1758       SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
1759       long aNodeId[4];
1760       gp_Pnt P[4];
1761
1762       double aLength;
1763       const SMDS_MeshElement* aNode;
1764       if(anIter->more()){
1765         aNode = anIter->next();
1766         const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1767         P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1768         aNodeId[0] = aNodeId[1] = aNode->GetID();
1769         aLength = 0;
1770       }
1771       for(; anIter->more(); ){
1772         const SMDS_MeshNode* N1 = static_cast<const SMDS_MeshNode*> (anIter->next());
1773         P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z());
1774         aNodeId[2] = N1->GetID();
1775         aLength = P[1].Distance(P[2]);
1776         if(!anIter->more()) break;
1777         const SMDS_MeshNode* N2 = static_cast<const SMDS_MeshNode*> (anIter->next());
1778         P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z());
1779         aNodeId[3] = N2->GetID();
1780         aLength += P[2].Distance(P[3]);
1781         Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1782         Value aValue2(aLength,aNodeId[2],aNodeId[3]);
1783         P[1] = P[3];
1784         aNodeId[1] = aNodeId[3];
1785         theValues.insert(aValue1);
1786         theValues.insert(aValue2);
1787       }
1788       aLength += P[2].Distance(P[0]);
1789       Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1790       Value aValue2(aLength,aNodeId[2],aNodeId[0]);
1791       theValues.insert(aValue1);
1792       theValues.insert(aValue2);
1793     }
1794     else {
1795       SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1796       long aNodeId[2];
1797       gp_Pnt P[3];
1798
1799       double aLength;
1800       const SMDS_MeshElement* aNode;
1801       if(aNodesIter->more()){
1802         aNode = aNodesIter->next();
1803         const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1804         P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1805         aNodeId[0] = aNodeId[1] = aNode->GetID();
1806         aLength = 0;
1807       }
1808       for(; aNodesIter->more(); ){
1809         aNode = aNodesIter->next();
1810         const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1811         long anId = aNode->GetID();
1812
1813         P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1814
1815         aLength = P[1].Distance(P[2]);
1816
1817         Value aValue(aLength,aNodeId[1],anId);
1818         aNodeId[1] = anId;
1819         P[1] = P[2];
1820         theValues.insert(aValue);
1821       }
1822
1823       aLength = P[0].Distance(P[1]);
1824
1825       Value aValue(aLength,aNodeId[0],aNodeId[1]);
1826       theValues.insert(aValue);
1827     }
1828   }
1829 }
1830
1831 //================================================================================
1832 /*
1833   Class       : MultiConnection
1834   Description : Functor for calculating number of faces conneted to the edge
1835 */
1836 //================================================================================
1837
1838 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
1839 {
1840   return 0;
1841 }
1842 double MultiConnection::GetValue( long theId )
1843 {
1844   return getNbMultiConnection( myMesh, theId );
1845 }
1846
1847 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
1848 {
1849   // meaningless as it is not quality control functor
1850   return Value;
1851 }
1852
1853 SMDSAbs_ElementType MultiConnection::GetType() const
1854 {
1855   return SMDSAbs_Edge;
1856 }
1857
1858 //================================================================================
1859 /*
1860   Class       : MultiConnection2D
1861   Description : Functor for calculating number of faces conneted to the edge
1862 */
1863 //================================================================================
1864
1865 double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
1866 {
1867   return 0;
1868 }
1869
1870 double MultiConnection2D::GetValue( long theElementId )
1871 {
1872   int aResult = 0;
1873
1874   const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
1875   SMDSAbs_ElementType aType = aFaceElem->GetType();
1876
1877   switch (aType) {
1878   case SMDSAbs_Face:
1879     {
1880       int i = 0, len = aFaceElem->NbNodes();
1881       SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
1882       if (!anIter) break;
1883
1884       const SMDS_MeshNode *aNode, *aNode0;
1885       TColStd_MapOfInteger aMap, aMapPrev;
1886
1887       for (i = 0; i <= len; i++) {
1888         aMapPrev = aMap;
1889         aMap.Clear();
1890
1891         int aNb = 0;
1892         if (anIter->more()) {
1893           aNode = (SMDS_MeshNode*)anIter->next();
1894         } else {
1895           if (i == len)
1896             aNode = aNode0;
1897           else
1898             break;
1899         }
1900         if (!aNode) break;
1901         if (i == 0) aNode0 = aNode;
1902
1903         SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
1904         while (anElemIter->more()) {
1905           const SMDS_MeshElement* anElem = anElemIter->next();
1906           if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
1907             int anId = anElem->GetID();
1908
1909             aMap.Add(anId);
1910             if (aMapPrev.Contains(anId)) {
1911               aNb++;
1912             }
1913           }
1914         }
1915         aResult = Max(aResult, aNb);
1916       }
1917     }
1918     break;
1919   default:
1920     aResult = 0;
1921   }
1922
1923   return aResult;
1924 }
1925
1926 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1927 {
1928   // meaningless as it is not quality control functor
1929   return Value;
1930 }
1931
1932 SMDSAbs_ElementType MultiConnection2D::GetType() const
1933 {
1934   return SMDSAbs_Face;
1935 }
1936
1937 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
1938 {
1939   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
1940   if(thePntId1 > thePntId2){
1941     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
1942   }
1943 }
1944
1945 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{
1946   if(myPntId[0] < x.myPntId[0]) return true;
1947   if(myPntId[0] == x.myPntId[0])
1948     if(myPntId[1] < x.myPntId[1]) return true;
1949   return false;
1950 }
1951
1952 void MultiConnection2D::GetValues(MValues& theValues){
1953   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1954   for(; anIter->more(); ){
1955     const SMDS_MeshFace* anElem = anIter->next();
1956     SMDS_ElemIteratorPtr aNodesIter;
1957     if ( anElem->IsQuadratic() )
1958       aNodesIter = dynamic_cast<const SMDS_VtkFace*>
1959         (anElem)->interlacedNodesElemIterator();
1960     else
1961       aNodesIter = anElem->nodesIterator();
1962     long aNodeId[3];
1963
1964     //int aNbConnects=0;
1965     const SMDS_MeshNode* aNode0;
1966     const SMDS_MeshNode* aNode1;
1967     const SMDS_MeshNode* aNode2;
1968     if(aNodesIter->more()){
1969       aNode0 = (SMDS_MeshNode*) aNodesIter->next();
1970       aNode1 = aNode0;
1971       const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
1972       aNodeId[0] = aNodeId[1] = aNodes->GetID();
1973     }
1974     for(; aNodesIter->more(); ) {
1975       aNode2 = (SMDS_MeshNode*) aNodesIter->next();
1976       long anId = aNode2->GetID();
1977       aNodeId[2] = anId;
1978
1979       Value aValue(aNodeId[1],aNodeId[2]);
1980       MValues::iterator aItr = theValues.find(aValue);
1981       if (aItr != theValues.end()){
1982         aItr->second += 1;
1983         //aNbConnects = nb;
1984       }
1985       else {
1986         theValues[aValue] = 1;
1987         //aNbConnects = 1;
1988       }
1989       //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1990       aNodeId[1] = aNodeId[2];
1991       aNode1 = aNode2;
1992     }
1993     Value aValue(aNodeId[0],aNodeId[2]);
1994     MValues::iterator aItr = theValues.find(aValue);
1995     if (aItr != theValues.end()) {
1996       aItr->second += 1;
1997       //aNbConnects = nb;
1998     }
1999     else {
2000       theValues[aValue] = 1;
2001       //aNbConnects = 1;
2002     }
2003     //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
2004   }
2005
2006 }
2007
2008 //================================================================================
2009 /*
2010   Class       : BallDiameter
2011   Description : Functor returning diameter of a ball element
2012 */
2013 //================================================================================
2014
2015 double BallDiameter::GetValue( long theId )
2016 {
2017   double diameter = 0;
2018
2019   if ( const SMDS_BallElement* ball =
2020        dynamic_cast<const SMDS_BallElement*>( myMesh->FindElement( theId )))
2021   {
2022     diameter = ball->GetDiameter();
2023   }
2024   return diameter;
2025 }
2026
2027 double BallDiameter::GetBadRate( double Value, int /*nbNodes*/ ) const
2028 {
2029   // meaningless as it is not a quality control functor
2030   return Value;
2031 }
2032
2033 SMDSAbs_ElementType BallDiameter::GetType() const
2034 {
2035   return SMDSAbs_Ball;
2036 }
2037
2038
2039 /*
2040                             PREDICATES
2041 */
2042
2043 //================================================================================
2044 /*
2045   Class       : BadOrientedVolume
2046   Description : Predicate bad oriented volumes
2047 */
2048 //================================================================================
2049
2050 BadOrientedVolume::BadOrientedVolume()
2051 {
2052   myMesh = 0;
2053 }
2054
2055 void BadOrientedVolume::SetMesh( const SMDS_Mesh* theMesh )
2056 {
2057   myMesh = theMesh;
2058 }
2059
2060 bool BadOrientedVolume::IsSatisfy( long theId )
2061 {
2062   if ( myMesh == 0 )
2063     return false;
2064
2065   SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
2066   return !vTool.IsForward();
2067 }
2068
2069 SMDSAbs_ElementType BadOrientedVolume::GetType() const
2070 {
2071   return SMDSAbs_Volume;
2072 }
2073
2074 /*
2075   Class       : BareBorderVolume
2076 */
2077
2078 bool BareBorderVolume::IsSatisfy(long theElementId )
2079 {
2080   SMDS_VolumeTool  myTool;
2081   if ( myTool.Set( myMesh->FindElement(theElementId)))
2082   {
2083     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2084       if ( myTool.IsFreeFace( iF ))
2085       {
2086         const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
2087         vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
2088         if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
2089           return true;
2090       }
2091   }
2092   return false;
2093 }
2094
2095 //================================================================================
2096 /*
2097   Class       : BareBorderFace
2098 */
2099 //================================================================================
2100
2101 bool BareBorderFace::IsSatisfy(long theElementId )
2102 {
2103   bool ok = false;
2104   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2105   {
2106     if ( face->GetType() == SMDSAbs_Face )
2107     {
2108       int nbN = face->NbCornerNodes();
2109       for ( int i = 0; i < nbN && !ok; ++i )
2110       {
2111         // check if a link is shared by another face
2112         const SMDS_MeshNode* n1 = face->GetNode( i );
2113         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2114         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2115         bool isShared = false;
2116         while ( !isShared && fIt->more() )
2117         {
2118           const SMDS_MeshElement* f = fIt->next();
2119           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2120         }
2121         if ( !isShared )
2122         {
2123           const int iQuad = face->IsQuadratic();
2124           myLinkNodes.resize( 2 + iQuad);
2125           myLinkNodes[0] = n1;
2126           myLinkNodes[1] = n2;
2127           if ( iQuad )
2128             myLinkNodes[2] = face->GetNode( i+nbN );
2129           ok = !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
2130         }
2131       }
2132     }
2133   }
2134   return ok;
2135 }
2136
2137 //================================================================================
2138 /*
2139   Class       : OverConstrainedVolume
2140 */
2141 //================================================================================
2142
2143 bool OverConstrainedVolume::IsSatisfy(long theElementId )
2144 {
2145   // An element is over-constrained if it has N-1 free borders where
2146   // N is the number of edges/faces for a 2D/3D element.
2147   SMDS_VolumeTool  myTool;
2148   if ( myTool.Set( myMesh->FindElement(theElementId)))
2149   {
2150     int nbSharedFaces = 0;
2151     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2152       if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
2153         break;
2154     return ( nbSharedFaces == 1 );
2155   }
2156   return false;
2157 }
2158
2159 //================================================================================
2160 /*
2161   Class       : OverConstrainedFace
2162 */
2163 //================================================================================
2164
2165 bool OverConstrainedFace::IsSatisfy(long theElementId )
2166 {
2167   // An element is over-constrained if it has N-1 free borders where
2168   // N is the number of edges/faces for a 2D/3D element.
2169   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2170     if ( face->GetType() == SMDSAbs_Face )
2171     {
2172       int nbSharedBorders = 0;
2173       int nbN = face->NbCornerNodes();
2174       for ( int i = 0; i < nbN; ++i )
2175       {
2176         // check if a link is shared by another face
2177         const SMDS_MeshNode* n1 = face->GetNode( i );
2178         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2179         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2180         bool isShared = false;
2181         while ( !isShared && fIt->more() )
2182         {
2183           const SMDS_MeshElement* f = fIt->next();
2184           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2185         }
2186         if ( isShared && ++nbSharedBorders > 1 )
2187           break;
2188       }
2189       return ( nbSharedBorders == 1 );
2190     }
2191   return false;
2192 }
2193
2194 //================================================================================
2195 /*
2196   Class       : CoincidentNodes
2197   Description : Predicate of Coincident nodes
2198 */
2199 //================================================================================
2200
2201 CoincidentNodes::CoincidentNodes()
2202 {
2203   myToler = 1e-5;
2204 }
2205
2206 bool CoincidentNodes::IsSatisfy( long theElementId )
2207 {
2208   return myCoincidentIDs.Contains( theElementId );
2209 }
2210
2211 SMDSAbs_ElementType CoincidentNodes::GetType() const
2212 {
2213   return SMDSAbs_Node;
2214 }
2215
2216 void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
2217 {
2218   myMeshModifTracer.SetMesh( theMesh );
2219   if ( myMeshModifTracer.IsMeshModified() )
2220   {
2221     TIDSortedNodeSet nodesToCheck;
2222     SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
2223     while ( nIt->more() )
2224       nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
2225
2226     list< list< const SMDS_MeshNode*> > nodeGroups;
2227     SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
2228
2229     myCoincidentIDs.Clear();
2230     list< list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
2231     for ( ; groupIt != nodeGroups.end(); ++groupIt )
2232     {
2233       list< const SMDS_MeshNode*>& coincNodes = *groupIt;
2234       list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
2235       for ( ; n != coincNodes.end(); ++n )
2236         myCoincidentIDs.Add( (*n)->GetID() );
2237     }
2238   }
2239 }
2240
2241 //================================================================================
2242 /*
2243   Class       : CoincidentElements
2244   Description : Predicate of Coincident Elements
2245   Note        : This class is suitable only for visualization of Coincident Elements
2246 */
2247 //================================================================================
2248
2249 CoincidentElements::CoincidentElements()
2250 {
2251   myMesh = 0;
2252 }
2253
2254 void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
2255 {
2256   myMesh = theMesh;
2257 }
2258
2259 bool CoincidentElements::IsSatisfy( long theElementId )
2260 {
2261   if ( !myMesh ) return false;
2262
2263   if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
2264   {
2265     if ( e->GetType() != GetType() ) return false;
2266     set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
2267     const int nbNodes = e->NbNodes();
2268     SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
2269     while ( invIt->more() )
2270     {
2271       const SMDS_MeshElement* e2 = invIt->next();
2272       if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
2273
2274       bool sameNodes = true;
2275       for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
2276         sameNodes = ( elemNodes.count( e2->GetNode( i )));
2277       if ( sameNodes )
2278         return true;
2279     }
2280   }
2281   return false;
2282 }
2283
2284 SMDSAbs_ElementType CoincidentElements1D::GetType() const
2285 {
2286   return SMDSAbs_Edge;
2287 }
2288 SMDSAbs_ElementType CoincidentElements2D::GetType() const
2289 {
2290   return SMDSAbs_Face;
2291 }
2292 SMDSAbs_ElementType CoincidentElements3D::GetType() const
2293 {
2294   return SMDSAbs_Volume;
2295 }
2296
2297
2298 //================================================================================
2299 /*
2300   Class       : FreeBorders
2301   Description : Predicate for free borders
2302 */
2303 //================================================================================
2304
2305 FreeBorders::FreeBorders()
2306 {
2307   myMesh = 0;
2308 }
2309
2310 void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
2311 {
2312   myMesh = theMesh;
2313 }
2314
2315 bool FreeBorders::IsSatisfy( long theId )
2316 {
2317   return getNbMultiConnection( myMesh, theId ) == 1;
2318 }
2319
2320 SMDSAbs_ElementType FreeBorders::GetType() const
2321 {
2322   return SMDSAbs_Edge;
2323 }
2324
2325
2326 //================================================================================
2327 /*
2328   Class       : FreeEdges
2329   Description : Predicate for free Edges
2330 */
2331 //================================================================================
2332
2333 FreeEdges::FreeEdges()
2334 {
2335   myMesh = 0;
2336 }
2337
2338 void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
2339 {
2340   myMesh = theMesh;
2341 }
2342
2343 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId  )
2344 {
2345   TColStd_MapOfInteger aMap;
2346   for ( int i = 0; i < 2; i++ )
2347   {
2348     SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face);
2349     while( anElemIter->more() )
2350     {
2351       if ( const SMDS_MeshElement* anElem = anElemIter->next())
2352       {
2353         const int anId = anElem->GetID();
2354         if ( anId != theFaceId && !aMap.Add( anId ))
2355           return false;
2356       }
2357     }
2358   }
2359   return true;
2360 }
2361
2362 bool FreeEdges::IsSatisfy( long theId )
2363 {
2364   if ( myMesh == 0 )
2365     return false;
2366
2367   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2368   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
2369     return false;
2370
2371   SMDS_ElemIteratorPtr anIter;
2372   if ( aFace->IsQuadratic() ) {
2373     anIter = dynamic_cast<const SMDS_VtkFace*>
2374       (aFace)->interlacedNodesElemIterator();
2375   }
2376   else {
2377     anIter = aFace->nodesIterator();
2378   }
2379   if ( !anIter )
2380     return false;
2381
2382   int i = 0, nbNodes = aFace->NbNodes();
2383   std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
2384   while( anIter->more() )
2385   {
2386     const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
2387     if ( aNode == 0 )
2388       return false;
2389     aNodes[ i++ ] = aNode;
2390   }
2391   aNodes[ nbNodes ] = aNodes[ 0 ];
2392
2393   for ( i = 0; i < nbNodes; i++ )
2394     if ( IsFreeEdge( &aNodes[ i ], theId ) )
2395       return true;
2396
2397   return false;
2398 }
2399
2400 SMDSAbs_ElementType FreeEdges::GetType() const
2401 {
2402   return SMDSAbs_Face;
2403 }
2404
2405 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
2406   myElemId(theElemId)
2407 {
2408   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
2409   if(thePntId1 > thePntId2){
2410     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
2411   }
2412 }
2413
2414 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
2415   if(myPntId[0] < x.myPntId[0]) return true;
2416   if(myPntId[0] == x.myPntId[0])
2417     if(myPntId[1] < x.myPntId[1]) return true;
2418   return false;
2419 }
2420
2421 inline void UpdateBorders(const FreeEdges::Border& theBorder,
2422                           FreeEdges::TBorders& theRegistry,
2423                           FreeEdges::TBorders& theContainer)
2424 {
2425   if(theRegistry.find(theBorder) == theRegistry.end()){
2426     theRegistry.insert(theBorder);
2427     theContainer.insert(theBorder);
2428   }else{
2429     theContainer.erase(theBorder);
2430   }
2431 }
2432
2433 void FreeEdges::GetBoreders(TBorders& theBorders)
2434 {
2435   TBorders aRegistry;
2436   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
2437   for(; anIter->more(); ){
2438     const SMDS_MeshFace* anElem = anIter->next();
2439     long anElemId = anElem->GetID();
2440     SMDS_ElemIteratorPtr aNodesIter;
2441     if ( anElem->IsQuadratic() )
2442       aNodesIter = static_cast<const SMDS_VtkFace*>(anElem)->
2443         interlacedNodesElemIterator();
2444     else
2445       aNodesIter = anElem->nodesIterator();
2446     long aNodeId[2];
2447     const SMDS_MeshElement* aNode;
2448     if(aNodesIter->more()){
2449       aNode = aNodesIter->next();
2450       aNodeId[0] = aNodeId[1] = aNode->GetID();
2451     }
2452     for(; aNodesIter->more(); ){
2453       aNode = aNodesIter->next();
2454       long anId = aNode->GetID();
2455       Border aBorder(anElemId,aNodeId[1],anId);
2456       aNodeId[1] = anId;
2457       UpdateBorders(aBorder,aRegistry,theBorders);
2458     }
2459     Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
2460     UpdateBorders(aBorder,aRegistry,theBorders);
2461   }
2462 }
2463
2464 //================================================================================
2465 /*
2466   Class       : FreeNodes
2467   Description : Predicate for free nodes
2468 */
2469 //================================================================================
2470
2471 FreeNodes::FreeNodes()
2472 {
2473   myMesh = 0;
2474 }
2475
2476 void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
2477 {
2478   myMesh = theMesh;
2479 }
2480
2481 bool FreeNodes::IsSatisfy( long theNodeId )
2482 {
2483   const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
2484   if (!aNode)
2485     return false;
2486
2487   return (aNode->NbInverseElements() < 1);
2488 }
2489
2490 SMDSAbs_ElementType FreeNodes::GetType() const
2491 {
2492   return SMDSAbs_Node;
2493 }
2494
2495
2496 //================================================================================
2497 /*
2498   Class       : FreeFaces
2499   Description : Predicate for free faces
2500 */
2501 //================================================================================
2502
2503 FreeFaces::FreeFaces()
2504 {
2505   myMesh = 0;
2506 }
2507
2508 void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
2509 {
2510   myMesh = theMesh;
2511 }
2512
2513 bool FreeFaces::IsSatisfy( long theId )
2514 {
2515   if (!myMesh) return false;
2516   // check that faces nodes refers to less than two common volumes
2517   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2518   if ( !aFace || aFace->GetType() != SMDSAbs_Face )
2519     return false;
2520
2521   int nbNode = aFace->NbNodes();
2522
2523   // collect volumes check that number of volumss with count equal nbNode not less than 2
2524   typedef map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
2525   typedef map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
2526   TMapOfVolume mapOfVol;
2527
2528   SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
2529   while ( nodeItr->more() ) {
2530     const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
2531     if ( !aNode ) continue;
2532     SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
2533     while ( volItr->more() ) {
2534       SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
2535       TItrMapOfVolume itr = mapOfVol.insert(make_pair(aVol, 0)).first;
2536       (*itr).second++;
2537     }
2538   }
2539   int nbVol = 0;
2540   TItrMapOfVolume volItr = mapOfVol.begin();
2541   TItrMapOfVolume volEnd = mapOfVol.end();
2542   for ( ; volItr != volEnd; ++volItr )
2543     if ( (*volItr).second >= nbNode )
2544        nbVol++;
2545   // face is not free if number of volumes constructed on thier nodes more than one
2546   return (nbVol < 2);
2547 }
2548
2549 SMDSAbs_ElementType FreeFaces::GetType() const
2550 {
2551   return SMDSAbs_Face;
2552 }
2553
2554 //================================================================================
2555 /*
2556   Class       : LinearOrQuadratic
2557   Description : Predicate to verify whether a mesh element is linear
2558 */
2559 //================================================================================
2560
2561 LinearOrQuadratic::LinearOrQuadratic()
2562 {
2563   myMesh = 0;
2564 }
2565
2566 void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
2567 {
2568   myMesh = theMesh;
2569 }
2570
2571 bool LinearOrQuadratic::IsSatisfy( long theId )
2572 {
2573   if (!myMesh) return false;
2574   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2575   if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
2576     return false;
2577   return (!anElem->IsQuadratic());
2578 }
2579
2580 void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
2581 {
2582   myType = theType;
2583 }
2584
2585 SMDSAbs_ElementType LinearOrQuadratic::GetType() const
2586 {
2587   return myType;
2588 }
2589
2590 //================================================================================
2591 /*
2592   Class       : GroupColor
2593   Description : Functor for check color of group to whic mesh element belongs to
2594 */
2595 //================================================================================
2596
2597 GroupColor::GroupColor()
2598 {
2599 }
2600
2601 bool GroupColor::IsSatisfy( long theId )
2602 {
2603   return (myIDs.find( theId ) != myIDs.end());
2604 }
2605
2606 void GroupColor::SetType( SMDSAbs_ElementType theType )
2607 {
2608   myType = theType;
2609 }
2610
2611 SMDSAbs_ElementType GroupColor::GetType() const
2612 {
2613   return myType;
2614 }
2615
2616 static bool isEqual( const Quantity_Color& theColor1,
2617                      const Quantity_Color& theColor2 )
2618 {
2619   // tolerance to compare colors
2620   const double tol = 5*1e-3;
2621   return ( fabs( theColor1.Red() - theColor2.Red() ) < tol &&
2622            fabs( theColor1.Green() - theColor2.Green() ) < tol &&
2623            fabs( theColor1.Blue() - theColor2.Blue() ) < tol );
2624 }
2625
2626
2627 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
2628 {
2629   myIDs.clear();
2630
2631   const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
2632   if ( !aMesh )
2633     return;
2634
2635   int nbGrp = aMesh->GetNbGroups();
2636   if ( !nbGrp )
2637     return;
2638
2639   // iterates on groups and find necessary elements ids
2640   const std::set<SMESHDS_GroupBase*>& aGroups = aMesh->GetGroups();
2641   set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
2642   for (; GrIt != aGroups.end(); GrIt++) {
2643     SMESHDS_GroupBase* aGrp = (*GrIt);
2644     if ( !aGrp )
2645       continue;
2646     // check type and color of group
2647     if ( !isEqual( myColor, aGrp->GetColor() ) )
2648       continue;
2649     if ( myType != SMDSAbs_All && myType != (SMDSAbs_ElementType)aGrp->GetType() )
2650       continue;
2651
2652     SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
2653     if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
2654       // add elements IDS into control
2655       int aSize = aGrp->Extent();
2656       for (int i = 0; i < aSize; i++)
2657         myIDs.insert( aGrp->GetID(i+1) );
2658     }
2659   }
2660 }
2661
2662 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
2663 {
2664   TCollection_AsciiString aStr = theStr;
2665   aStr.RemoveAll( ' ' );
2666   aStr.RemoveAll( '\t' );
2667   for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
2668     aStr.Remove( aPos, 2 );
2669   Standard_Real clr[3];
2670   clr[0] = clr[1] = clr[2] = 0.;
2671   for ( int i = 0; i < 3; i++ ) {
2672     TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
2673     if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
2674       clr[i] = tmpStr.RealValue();
2675   }
2676   myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
2677 }
2678
2679 //=======================================================================
2680 // name    : GetRangeStr
2681 // Purpose : Get range as a string.
2682 //           Example: "1,2,3,50-60,63,67,70-"
2683 //=======================================================================
2684
2685 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
2686 {
2687   theResStr.Clear();
2688   theResStr += TCollection_AsciiString( myColor.Red() );
2689   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
2690   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
2691 }
2692
2693 //================================================================================
2694 /*
2695   Class       : ElemGeomType
2696   Description : Predicate to check element geometry type
2697 */
2698 //================================================================================
2699
2700 ElemGeomType::ElemGeomType()
2701 {
2702   myMesh = 0;
2703   myType = SMDSAbs_All;
2704   myGeomType = SMDSGeom_TRIANGLE;
2705 }
2706
2707 void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
2708 {
2709   myMesh = theMesh;
2710 }
2711
2712 bool ElemGeomType::IsSatisfy( long theId )
2713 {
2714   if (!myMesh) return false;
2715   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2716   if ( !anElem )
2717     return false;
2718   const SMDSAbs_ElementType anElemType = anElem->GetType();
2719   if ( myType != SMDSAbs_All && anElemType != myType )
2720     return false;
2721   bool isOk = ( anElem->GetGeomType() == myGeomType );
2722   return isOk;
2723 }
2724
2725 void ElemGeomType::SetType( SMDSAbs_ElementType theType )
2726 {
2727   myType = theType;
2728 }
2729
2730 SMDSAbs_ElementType ElemGeomType::GetType() const
2731 {
2732   return myType;
2733 }
2734
2735 void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
2736 {
2737   myGeomType = theType;
2738 }
2739
2740 SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
2741 {
2742   return myGeomType;
2743 }
2744
2745 //================================================================================
2746 /*
2747   Class       : ElemEntityType
2748   Description : Predicate to check element entity type
2749 */
2750 //================================================================================
2751
2752 ElemEntityType::ElemEntityType():
2753   myMesh( 0 ),
2754   myType( SMDSAbs_All ),
2755   myEntityType( SMDSEntity_0D )
2756 {
2757 }
2758
2759 void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh )
2760 {
2761   myMesh = theMesh;
2762 }
2763
2764 bool ElemEntityType::IsSatisfy( long theId )
2765 {
2766   if ( !myMesh ) return false;
2767   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2768   return ( anElem &&
2769            myEntityType == anElem->GetEntityType() &&
2770            ( myType == SMDSAbs_Edge || myType == SMDSAbs_Face || myType ==  SMDSAbs_Volume ));
2771 }
2772
2773 void ElemEntityType::SetType( SMDSAbs_ElementType theType )
2774 {
2775   myType = theType;
2776 }
2777
2778 SMDSAbs_ElementType ElemEntityType::GetType() const
2779 {
2780   return myType;
2781 }
2782
2783 void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType )
2784 {
2785   myEntityType = theEntityType;
2786 }
2787
2788 SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const
2789 {
2790   return myEntityType;
2791 }
2792
2793 //================================================================================
2794 /*!
2795  * \brief Class CoplanarFaces
2796  */
2797 //================================================================================
2798
2799 CoplanarFaces::CoplanarFaces()
2800   : myFaceID(0), myToler(0)
2801 {
2802 }
2803 void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
2804 {
2805   myMeshModifTracer.SetMesh( theMesh );
2806   if ( myMeshModifTracer.IsMeshModified() )
2807   {
2808     // Build a set of coplanar face ids
2809
2810     myCoplanarIDs.clear();
2811
2812     if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler )
2813       return;
2814
2815     const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID );
2816     if ( !face || face->GetType() != SMDSAbs_Face )
2817       return;
2818
2819     bool normOK;
2820     gp_Vec myNorm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
2821     if (!normOK)
2822       return;
2823
2824     const double radianTol = myToler * M_PI / 180.;
2825     std::set< SMESH_TLink > checkedLinks;
2826
2827     std::list< pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
2828     faceQueue.push_back( make_pair( face, myNorm ));
2829     while ( !faceQueue.empty() )
2830     {
2831       face   = faceQueue.front().first;
2832       myNorm = faceQueue.front().second;
2833       faceQueue.pop_front();
2834
2835       for ( int i = 0, nbN = face->NbCornerNodes(); i < nbN; ++i )
2836       {
2837         const SMDS_MeshNode*  n1 = face->GetNode( i );
2838         const SMDS_MeshNode*  n2 = face->GetNode(( i+1 )%nbN);
2839         if ( !checkedLinks.insert( SMESH_TLink( n1, n2 )).second )
2840           continue;
2841         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face);
2842         while ( fIt->more() )
2843         {
2844           const SMDS_MeshElement* f = fIt->next();
2845           if ( f->GetNodeIndex( n2 ) > -1 )
2846           {
2847             gp_Vec norm = getNormale( static_cast<const SMDS_MeshFace*>(f), &normOK );
2848             if (!normOK || myNorm.Angle( norm ) <= radianTol)
2849             {
2850               myCoplanarIDs.insert( f->GetID() );
2851               faceQueue.push_back( make_pair( f, norm ));
2852             }
2853           }
2854         }
2855       }
2856     }
2857   }
2858 }
2859 bool CoplanarFaces::IsSatisfy( long theElementId )
2860 {
2861   return myCoplanarIDs.count( theElementId );
2862 }
2863
2864 /*
2865  *Class       : RangeOfIds
2866   *Description : Predicate for Range of Ids.
2867   *              Range may be specified with two ways.
2868   *              1. Using AddToRange method
2869   *              2. With SetRangeStr method. Parameter of this method is a string
2870   *                 like as "1,2,3,50-60,63,67,70-"
2871 */
2872
2873 //=======================================================================
2874 // name    : RangeOfIds
2875 // Purpose : Constructor
2876 //=======================================================================
2877 RangeOfIds::RangeOfIds()
2878 {
2879   myMesh = 0;
2880   myType = SMDSAbs_All;
2881 }
2882
2883 //=======================================================================
2884 // name    : SetMesh
2885 // Purpose : Set mesh
2886 //=======================================================================
2887 void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh )
2888 {
2889   myMesh = theMesh;
2890 }
2891
2892 //=======================================================================
2893 // name    : AddToRange
2894 // Purpose : Add ID to the range
2895 //=======================================================================
2896 bool RangeOfIds::AddToRange( long theEntityId )
2897 {
2898   myIds.Add( theEntityId );
2899   return true;
2900 }
2901
2902 //=======================================================================
2903 // name    : GetRangeStr
2904 // Purpose : Get range as a string.
2905 //           Example: "1,2,3,50-60,63,67,70-"
2906 //=======================================================================
2907 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
2908 {
2909   theResStr.Clear();
2910
2911   TColStd_SequenceOfInteger     anIntSeq;
2912   TColStd_SequenceOfAsciiString aStrSeq;
2913
2914   TColStd_MapIteratorOfMapOfInteger anIter( myIds );
2915   for ( ; anIter.More(); anIter.Next() )
2916   {
2917     int anId = anIter.Key();
2918     TCollection_AsciiString aStr( anId );
2919     anIntSeq.Append( anId );
2920     aStrSeq.Append( aStr );
2921   }
2922
2923   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
2924   {
2925     int aMinId = myMin( i );
2926     int aMaxId = myMax( i );
2927
2928     TCollection_AsciiString aStr;
2929     if ( aMinId != IntegerFirst() )
2930       aStr += aMinId;
2931
2932     aStr += "-";
2933
2934     if ( aMaxId != IntegerLast() )
2935       aStr += aMaxId;
2936
2937     // find position of the string in result sequence and insert string in it
2938     if ( anIntSeq.Length() == 0 )
2939     {
2940       anIntSeq.Append( aMinId );
2941       aStrSeq.Append( aStr );
2942     }
2943     else
2944     {
2945       if ( aMinId < anIntSeq.First() )
2946       {
2947         anIntSeq.Prepend( aMinId );
2948         aStrSeq.Prepend( aStr );
2949       }
2950       else if ( aMinId > anIntSeq.Last() )
2951       {
2952         anIntSeq.Append( aMinId );
2953         aStrSeq.Append( aStr );
2954       }
2955       else
2956         for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
2957           if ( aMinId < anIntSeq( j ) )
2958           {
2959             anIntSeq.InsertBefore( j, aMinId );
2960             aStrSeq.InsertBefore( j, aStr );
2961             break;
2962           }
2963     }
2964   }
2965
2966   if ( aStrSeq.Length() == 0 )
2967     return;
2968
2969   theResStr = aStrSeq( 1 );
2970   for ( int j = 2, k = aStrSeq.Length(); j <= k; j++  )
2971   {
2972     theResStr += ",";
2973     theResStr += aStrSeq( j );
2974   }
2975 }
2976
2977 //=======================================================================
2978 // name    : SetRangeStr
2979 // Purpose : Define range with string
2980 //           Example of entry string: "1,2,3,50-60,63,67,70-"
2981 //=======================================================================
2982 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
2983 {
2984   myMin.Clear();
2985   myMax.Clear();
2986   myIds.Clear();
2987
2988   TCollection_AsciiString aStr = theStr;
2989   aStr.RemoveAll( ' ' );
2990   aStr.RemoveAll( '\t' );
2991
2992   for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
2993     aStr.Remove( aPos, 2 );
2994
2995   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
2996   int i = 1;
2997   while ( tmpStr != "" )
2998   {
2999     tmpStr = aStr.Token( ",", i++ );
3000     int aPos = tmpStr.Search( '-' );
3001
3002     if ( aPos == -1 )
3003     {
3004       if ( tmpStr.IsIntegerValue() )
3005         myIds.Add( tmpStr.IntegerValue() );
3006       else
3007         return false;
3008     }
3009     else
3010     {
3011       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
3012       TCollection_AsciiString aMinStr = tmpStr;
3013
3014       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
3015       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
3016
3017       if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) ||
3018            (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
3019         return false;
3020
3021       myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
3022       myMax.Append( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
3023     }
3024   }
3025
3026   return true;
3027 }
3028
3029 //=======================================================================
3030 // name    : GetType
3031 // Purpose : Get type of supported entities
3032 //=======================================================================
3033 SMDSAbs_ElementType RangeOfIds::GetType() const
3034 {
3035   return myType;
3036 }
3037
3038 //=======================================================================
3039 // name    : SetType
3040 // Purpose : Set type of supported entities
3041 //=======================================================================
3042 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
3043 {
3044   myType = theType;
3045 }
3046
3047 //=======================================================================
3048 // name    : IsSatisfy
3049 // Purpose : Verify whether entity satisfies to this rpedicate
3050 //=======================================================================
3051 bool RangeOfIds::IsSatisfy( long theId )
3052 {
3053   if ( !myMesh )
3054     return false;
3055
3056   if ( myType == SMDSAbs_Node )
3057   {
3058     if ( myMesh->FindNode( theId ) == 0 )
3059       return false;
3060   }
3061   else
3062   {
3063     const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3064     if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All ))
3065       return false;
3066   }
3067
3068   if ( myIds.Contains( theId ) )
3069     return true;
3070
3071   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
3072     if ( theId >= myMin( i ) && theId <= myMax( i ) )
3073       return true;
3074
3075   return false;
3076 }
3077
3078 /*
3079   Class       : Comparator
3080   Description : Base class for comparators
3081 */
3082 Comparator::Comparator():
3083   myMargin(0)
3084 {}
3085
3086 Comparator::~Comparator()
3087 {}
3088
3089 void Comparator::SetMesh( const SMDS_Mesh* theMesh )
3090 {
3091   if ( myFunctor )
3092     myFunctor->SetMesh( theMesh );
3093 }
3094
3095 void Comparator::SetMargin( double theValue )
3096 {
3097   myMargin = theValue;
3098 }
3099
3100 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
3101 {
3102   myFunctor = theFunct;
3103 }
3104
3105 SMDSAbs_ElementType Comparator::GetType() const
3106 {
3107   return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
3108 }
3109
3110 double Comparator::GetMargin()
3111 {
3112   return myMargin;
3113 }
3114
3115
3116 /*
3117   Class       : LessThan
3118   Description : Comparator "<"
3119 */
3120 bool LessThan::IsSatisfy( long theId )
3121 {
3122   return myFunctor && myFunctor->GetValue( theId ) < myMargin;
3123 }
3124
3125
3126 /*
3127   Class       : MoreThan
3128   Description : Comparator ">"
3129 */
3130 bool MoreThan::IsSatisfy( long theId )
3131 {
3132   return myFunctor && myFunctor->GetValue( theId ) > myMargin;
3133 }
3134
3135
3136 /*
3137   Class       : EqualTo
3138   Description : Comparator "="
3139 */
3140 EqualTo::EqualTo():
3141   myToler(Precision::Confusion())
3142 {}
3143
3144 bool EqualTo::IsSatisfy( long theId )
3145 {
3146   return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
3147 }
3148
3149 void EqualTo::SetTolerance( double theToler )
3150 {
3151   myToler = theToler;
3152 }
3153
3154 double EqualTo::GetTolerance()
3155 {
3156   return myToler;
3157 }
3158
3159 /*
3160   Class       : LogicalNOT
3161   Description : Logical NOT predicate
3162 */
3163 LogicalNOT::LogicalNOT()
3164 {}
3165
3166 LogicalNOT::~LogicalNOT()
3167 {}
3168
3169 bool LogicalNOT::IsSatisfy( long theId )
3170 {
3171   return myPredicate && !myPredicate->IsSatisfy( theId );
3172 }
3173
3174 void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
3175 {
3176   if ( myPredicate )
3177     myPredicate->SetMesh( theMesh );
3178 }
3179
3180 void LogicalNOT::SetPredicate( PredicatePtr thePred )
3181 {
3182   myPredicate = thePred;
3183 }
3184
3185 SMDSAbs_ElementType LogicalNOT::GetType() const
3186 {
3187   return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
3188 }
3189
3190
3191 /*
3192   Class       : LogicalBinary
3193   Description : Base class for binary logical predicate
3194 */
3195 LogicalBinary::LogicalBinary()
3196 {}
3197
3198 LogicalBinary::~LogicalBinary()
3199 {}
3200
3201 void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
3202 {
3203   if ( myPredicate1 )
3204     myPredicate1->SetMesh( theMesh );
3205
3206   if ( myPredicate2 )
3207     myPredicate2->SetMesh( theMesh );
3208 }
3209
3210 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
3211 {
3212   myPredicate1 = thePredicate;
3213 }
3214
3215 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
3216 {
3217   myPredicate2 = thePredicate;
3218 }
3219
3220 SMDSAbs_ElementType LogicalBinary::GetType() const
3221 {
3222   if ( !myPredicate1 || !myPredicate2 )
3223     return SMDSAbs_All;
3224
3225   SMDSAbs_ElementType aType1 = myPredicate1->GetType();
3226   SMDSAbs_ElementType aType2 = myPredicate2->GetType();
3227
3228   return aType1 == aType2 ? aType1 : SMDSAbs_All;
3229 }
3230
3231
3232 /*
3233   Class       : LogicalAND
3234   Description : Logical AND
3235 */
3236 bool LogicalAND::IsSatisfy( long theId )
3237 {
3238   return
3239     myPredicate1 &&
3240     myPredicate2 &&
3241     myPredicate1->IsSatisfy( theId ) &&
3242     myPredicate2->IsSatisfy( theId );
3243 }
3244
3245
3246 /*
3247   Class       : LogicalOR
3248   Description : Logical OR
3249 */
3250 bool LogicalOR::IsSatisfy( long theId )
3251 {
3252   return
3253     myPredicate1 &&
3254     myPredicate2 &&
3255     (myPredicate1->IsSatisfy( theId ) ||
3256     myPredicate2->IsSatisfy( theId ));
3257 }
3258
3259
3260 /*
3261                               FILTER
3262 */
3263
3264 Filter::Filter()
3265 {}
3266
3267 Filter::~Filter()
3268 {}
3269
3270 void Filter::SetPredicate( PredicatePtr thePredicate )
3271 {
3272   myPredicate = thePredicate;
3273 }
3274
3275 void Filter::GetElementsId( const SMDS_Mesh* theMesh,
3276                             PredicatePtr     thePredicate,
3277                             TIdSequence&     theSequence )
3278 {
3279   theSequence.clear();
3280
3281   if ( !theMesh || !thePredicate )
3282     return;
3283
3284   thePredicate->SetMesh( theMesh );
3285
3286   SMDS_ElemIteratorPtr elemIt = theMesh->elementsIterator( thePredicate->GetType() );
3287   if ( elemIt ) {
3288     while ( elemIt->more() ) {
3289       const SMDS_MeshElement* anElem = elemIt->next();
3290       long anId = anElem->GetID();
3291       if ( thePredicate->IsSatisfy( anId ) )
3292         theSequence.push_back( anId );
3293     }
3294   }
3295 }
3296
3297 void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
3298                             Filter::TIdSequence& theSequence )
3299 {
3300   GetElementsId(theMesh,myPredicate,theSequence);
3301 }
3302
3303 /*
3304                               ManifoldPart
3305 */
3306
3307 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
3308
3309 /*
3310    Internal class Link
3311 */
3312
3313 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
3314                           SMDS_MeshNode* theNode2 )
3315 {
3316   myNode1 = theNode1;
3317   myNode2 = theNode2;
3318 }
3319
3320 ManifoldPart::Link::~Link()
3321 {
3322   myNode1 = 0;
3323   myNode2 = 0;
3324 }
3325
3326 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
3327 {
3328   if ( myNode1 == theLink.myNode1 &&
3329        myNode2 == theLink.myNode2 )
3330     return true;
3331   else if ( myNode1 == theLink.myNode2 &&
3332             myNode2 == theLink.myNode1 )
3333     return true;
3334   else
3335     return false;
3336 }
3337
3338 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
3339 {
3340   if(myNode1 < x.myNode1) return true;
3341   if(myNode1 == x.myNode1)
3342     if(myNode2 < x.myNode2) return true;
3343   return false;
3344 }
3345
3346 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
3347                             const ManifoldPart::Link& theLink2 )
3348 {
3349   return theLink1.IsEqual( theLink2 );
3350 }
3351
3352 ManifoldPart::ManifoldPart()
3353 {
3354   myMesh = 0;
3355   myAngToler = Precision::Angular();
3356   myIsOnlyManifold = true;
3357 }
3358
3359 ManifoldPart::~ManifoldPart()
3360 {
3361   myMesh = 0;
3362 }
3363
3364 void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
3365 {
3366   myMesh = theMesh;
3367   process();
3368 }
3369
3370 SMDSAbs_ElementType ManifoldPart::GetType() const
3371 { return SMDSAbs_Face; }
3372
3373 bool ManifoldPart::IsSatisfy( long theElementId )
3374 {
3375   return myMapIds.Contains( theElementId );
3376 }
3377
3378 void ManifoldPart::SetAngleTolerance( const double theAngToler )
3379 { myAngToler = theAngToler; }
3380
3381 double ManifoldPart::GetAngleTolerance() const
3382 { return myAngToler; }
3383
3384 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
3385 { myIsOnlyManifold = theIsOnly; }
3386
3387 void ManifoldPart::SetStartElem( const long  theStartId )
3388 { myStartElemId = theStartId; }
3389
3390 bool ManifoldPart::process()
3391 {
3392   myMapIds.Clear();
3393   myMapBadGeomIds.Clear();
3394
3395   myAllFacePtr.clear();
3396   myAllFacePtrIntDMap.clear();
3397   if ( !myMesh )
3398     return false;
3399
3400   // collect all faces into own map
3401   SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
3402   for (; anFaceItr->more(); )
3403   {
3404     SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
3405     myAllFacePtr.push_back( aFacePtr );
3406     myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
3407   }
3408
3409   SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
3410   if ( !aStartFace )
3411     return false;
3412
3413   // the map of non manifold links and bad geometry
3414   TMapOfLink aMapOfNonManifold;
3415   TColStd_MapOfInteger aMapOfTreated;
3416
3417   // begin cycle on faces from start index and run on vector till the end
3418   //  and from begin to start index to cover whole vector
3419   const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
3420   bool isStartTreat = false;
3421   for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
3422   {
3423     if ( fi == aStartIndx )
3424       isStartTreat = true;
3425     // as result next time when fi will be equal to aStartIndx
3426
3427     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
3428     if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
3429       continue;
3430
3431     aMapOfTreated.Add( aFacePtr->GetID() );
3432     TColStd_MapOfInteger aResFaces;
3433     if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
3434                          aMapOfNonManifold, aResFaces ) )
3435       continue;
3436     TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
3437     for ( ; anItr.More(); anItr.Next() )
3438     {
3439       int aFaceId = anItr.Key();
3440       aMapOfTreated.Add( aFaceId );
3441       myMapIds.Add( aFaceId );
3442     }
3443
3444     if ( fi == ( myAllFacePtr.size() - 1 ) )
3445       fi = 0;
3446   } // end run on vector of faces
3447   return !myMapIds.IsEmpty();
3448 }
3449
3450 static void getLinks( const SMDS_MeshFace* theFace,
3451                       ManifoldPart::TVectorOfLink& theLinks )
3452 {
3453   int aNbNode = theFace->NbNodes();
3454   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
3455   int i = 1;
3456   SMDS_MeshNode* aNode = 0;
3457   for ( ; aNodeItr->more() && i <= aNbNode; )
3458   {
3459
3460     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
3461     if ( i == 1 )
3462       aNode = aN1;
3463     i++;
3464     SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
3465     i++;
3466     ManifoldPart::Link aLink( aN1, aN2 );
3467     theLinks.push_back( aLink );
3468   }
3469 }
3470
3471 bool ManifoldPart::findConnected
3472                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
3473                   SMDS_MeshFace*                           theStartFace,
3474                   ManifoldPart::TMapOfLink&                theNonManifold,
3475                   TColStd_MapOfInteger&                    theResFaces )
3476 {
3477   theResFaces.Clear();
3478   if ( !theAllFacePtrInt.size() )
3479     return false;
3480
3481   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
3482   {
3483     myMapBadGeomIds.Add( theStartFace->GetID() );
3484     return false;
3485   }
3486
3487   ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
3488   ManifoldPart::TVectorOfLink aSeqOfBoundary;
3489   theResFaces.Add( theStartFace->GetID() );
3490   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
3491
3492   expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3493                  aDMapLinkFace, theNonManifold, theStartFace );
3494
3495   bool isDone = false;
3496   while ( !isDone && aMapOfBoundary.size() != 0 )
3497   {
3498     bool isToReset = false;
3499     ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
3500     for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
3501     {
3502       ManifoldPart::Link aLink = *pLink;
3503       if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
3504         continue;
3505       // each link could be treated only once
3506       aMapToSkip.insert( aLink );
3507
3508       ManifoldPart::TVectorOfFacePtr aFaces;
3509       // find next
3510       if ( myIsOnlyManifold &&
3511            (theNonManifold.find( aLink ) != theNonManifold.end()) )
3512         continue;
3513       else
3514       {
3515         getFacesByLink( aLink, aFaces );
3516         // filter the element to keep only indicated elements
3517         ManifoldPart::TVectorOfFacePtr aFiltered;
3518         ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3519         for ( ; pFace != aFaces.end(); ++pFace )
3520         {
3521           SMDS_MeshFace* aFace = *pFace;
3522           if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
3523             aFiltered.push_back( aFace );
3524         }
3525         aFaces = aFiltered;
3526         if ( aFaces.size() < 2 )  // no neihgbour faces
3527           continue;
3528         else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
3529         {
3530           theNonManifold.insert( aLink );
3531           continue;
3532         }
3533       }
3534
3535       // compare normal with normals of neighbor element
3536       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
3537       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3538       for ( ; pFace != aFaces.end(); ++pFace )
3539       {
3540         SMDS_MeshFace* aNextFace = *pFace;
3541         if ( aPrevFace == aNextFace )
3542           continue;
3543         int anNextFaceID = aNextFace->GetID();
3544         if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
3545          // should not be with non manifold restriction. probably bad topology
3546           continue;
3547         // check if face was treated and skipped
3548         if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
3549              !isInPlane( aPrevFace, aNextFace ) )
3550           continue;
3551         // add new element to connected and extend the boundaries.
3552         theResFaces.Add( anNextFaceID );
3553         expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3554                         aDMapLinkFace, theNonManifold, aNextFace );
3555         isToReset = true;
3556       }
3557     }
3558     isDone = !isToReset;
3559   }
3560
3561   return !theResFaces.IsEmpty();
3562 }
3563
3564 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
3565                               const SMDS_MeshFace* theFace2 )
3566 {
3567   gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
3568   gp_XYZ aNorm2XYZ = getNormale( theFace2 );
3569   if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
3570   {
3571     myMapBadGeomIds.Add( theFace2->GetID() );
3572     return false;
3573   }
3574   if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
3575     return true;
3576
3577   return false;
3578 }
3579
3580 void ManifoldPart::expandBoundary
3581                    ( ManifoldPart::TMapOfLink&            theMapOfBoundary,
3582                      ManifoldPart::TVectorOfLink&         theSeqOfBoundary,
3583                      ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
3584                      ManifoldPart::TMapOfLink&            theNonManifold,
3585                      SMDS_MeshFace*                       theNextFace ) const
3586 {
3587   ManifoldPart::TVectorOfLink aLinks;
3588   getLinks( theNextFace, aLinks );
3589   int aNbLink = (int)aLinks.size();
3590   for ( int i = 0; i < aNbLink; i++ )
3591   {
3592     ManifoldPart::Link aLink = aLinks[ i ];
3593     if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
3594       continue;
3595     if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
3596     {
3597       if ( myIsOnlyManifold )
3598       {
3599         // remove from boundary
3600         theMapOfBoundary.erase( aLink );
3601         ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
3602         for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
3603         {
3604           ManifoldPart::Link aBoundLink = *pLink;
3605           if ( aBoundLink.IsEqual( aLink ) )
3606           {
3607             theSeqOfBoundary.erase( pLink );
3608             break;
3609           }
3610         }
3611       }
3612     }
3613     else
3614     {
3615       theMapOfBoundary.insert( aLink );
3616       theSeqOfBoundary.push_back( aLink );
3617       theDMapLinkFacePtr[ aLink ] = theNextFace;
3618     }
3619   }
3620 }
3621
3622 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
3623                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
3624 {
3625   std::set<SMDS_MeshCell *> aSetOfFaces;
3626   // take all faces that shared first node
3627   SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
3628   for ( ; anItr->more(); )
3629   {
3630     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3631     if ( !aFace )
3632       continue;
3633     aSetOfFaces.insert( aFace );
3634   }
3635   // take all faces that shared second node
3636   anItr = theLink.myNode2->facesIterator();
3637   // find the common part of two sets
3638   for ( ; anItr->more(); )
3639   {
3640     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3641     if ( aSetOfFaces.count( aFace ) )
3642       theFaces.push_back( aFace );
3643   }
3644 }
3645
3646
3647 /*
3648    ElementsOnSurface
3649 */
3650
3651 ElementsOnSurface::ElementsOnSurface()
3652 {
3653   myIds.Clear();
3654   myType = SMDSAbs_All;
3655   mySurf.Nullify();
3656   myToler = Precision::Confusion();
3657   myUseBoundaries = false;
3658 }
3659
3660 ElementsOnSurface::~ElementsOnSurface()
3661 {
3662 }
3663
3664 void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
3665 {
3666   myMeshModifTracer.SetMesh( theMesh );
3667   if ( myMeshModifTracer.IsMeshModified())
3668     process();
3669 }
3670
3671 bool ElementsOnSurface::IsSatisfy( long theElementId )
3672 {
3673   return myIds.Contains( theElementId );
3674 }
3675
3676 SMDSAbs_ElementType ElementsOnSurface::GetType() const
3677 { return myType; }
3678
3679 void ElementsOnSurface::SetTolerance( const double theToler )
3680 {
3681   if ( myToler != theToler )
3682     myIds.Clear();
3683   myToler = theToler;
3684 }
3685
3686 double ElementsOnSurface::GetTolerance() const
3687 { return myToler; }
3688
3689 void ElementsOnSurface::SetUseBoundaries( bool theUse )
3690 {
3691   if ( myUseBoundaries != theUse ) {
3692     myUseBoundaries = theUse;
3693     SetSurface( mySurf, myType );
3694   }
3695 }
3696
3697 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
3698                                     const SMDSAbs_ElementType theType )
3699 {
3700   myIds.Clear();
3701   myType = theType;
3702   mySurf.Nullify();
3703   if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
3704     return;
3705   mySurf = TopoDS::Face( theShape );
3706   BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
3707   Standard_Real
3708     u1 = SA.FirstUParameter(),
3709     u2 = SA.LastUParameter(),
3710     v1 = SA.FirstVParameter(),
3711     v2 = SA.LastVParameter();
3712   Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
3713   myProjector.Init( surf, u1,u2, v1,v2 );
3714   process();
3715 }
3716
3717 void ElementsOnSurface::process()
3718 {
3719   myIds.Clear();
3720   if ( mySurf.IsNull() )
3721     return;
3722
3723   if ( !myMeshModifTracer.GetMesh() )
3724     return;
3725
3726   myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
3727
3728   SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType );
3729   for(; anIter->more(); )
3730     process( anIter->next() );
3731 }
3732
3733 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
3734 {
3735   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
3736   bool isSatisfy = true;
3737   for ( ; aNodeItr->more(); )
3738   {
3739     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
3740     if ( !isOnSurface( aNode ) )
3741     {
3742       isSatisfy = false;
3743       break;
3744     }
3745   }
3746   if ( isSatisfy )
3747     myIds.Add( theElemPtr->GetID() );
3748 }
3749
3750 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
3751 {
3752   if ( mySurf.IsNull() )
3753     return false;
3754
3755   gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
3756   //  double aToler2 = myToler * myToler;
3757 //   if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
3758 //   {
3759 //     gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
3760 //     if ( aPln.SquareDistance( aPnt ) > aToler2 )
3761 //       return false;
3762 //   }
3763 //   else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
3764 //   {
3765 //     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
3766 //     double aRad = aCyl.Radius();
3767 //     gp_Ax3 anAxis = aCyl.Position();
3768 //     gp_XYZ aLoc = aCyl.Location().XYZ();
3769 //     double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
3770 //     double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
3771 //     if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
3772 //       return false;
3773 //   }
3774 //   else
3775 //     return false;
3776   myProjector.Perform( aPnt );
3777   bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
3778
3779   return isOn;
3780 }
3781
3782
3783 /*
3784   ElementsOnShape
3785 */
3786
3787 ElementsOnShape::ElementsOnShape()
3788   : //myMesh(0),
3789     myType(SMDSAbs_All),
3790     myToler(Precision::Confusion()),
3791     myAllNodesFlag(false)
3792 {
3793 }
3794
3795 ElementsOnShape::~ElementsOnShape()
3796 {
3797   clearClassifiers();
3798 }
3799
3800 SMDSAbs_ElementType ElementsOnShape::GetType() const
3801 {
3802   return myType;
3803 }
3804
3805 void ElementsOnShape::SetTolerance (const double theToler)
3806 {
3807   if (myToler != theToler) {
3808     myToler = theToler;
3809     SetShape(myShape, myType);
3810   }
3811 }
3812
3813 double ElementsOnShape::GetTolerance() const
3814 {
3815   return myToler;
3816 }
3817
3818 void ElementsOnShape::SetAllNodes (bool theAllNodes)
3819 {
3820   myAllNodesFlag = theAllNodes;
3821 }
3822
3823 void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
3824 {
3825   myMesh = theMesh;
3826 }
3827
3828 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
3829                                 const SMDSAbs_ElementType theType)
3830 {
3831   myType  = theType;
3832   myShape = theShape;
3833   if ( myShape.IsNull() ) return;
3834
3835   TopTools_IndexedMapOfShape shapesMap;
3836   TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
3837   TopExp_Explorer sub;
3838   for ( int i = 0; i < 4; ++i )
3839   {
3840     if ( shapesMap.IsEmpty() )
3841       for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() )
3842         shapesMap.Add( sub.Current() );
3843     if ( i > 0 )
3844       for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() )
3845         shapesMap.Add( sub.Current() );
3846   }
3847
3848   clearClassifiers();
3849   myClassifiers.resize( shapesMap.Extent() );
3850   for ( int i = 0; i < shapesMap.Extent(); ++i )
3851     myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler );
3852 }
3853
3854 void ElementsOnShape::clearClassifiers()
3855 {
3856   for ( size_t i = 0; i < myClassifiers.size(); ++i )
3857     delete myClassifiers[ i ];
3858   myClassifiers.clear();
3859 }
3860
3861 bool ElementsOnShape::IsSatisfy (long elemId)
3862 {
3863   const SMDS_MeshElement* elem =
3864     ( myType == SMDSAbs_Node ? myMesh->FindNode( elemId ) : myMesh->FindElement( elemId ));
3865   if ( !elem || myClassifiers.empty() )
3866     return false;
3867
3868   for ( size_t i = 0; i < myClassifiers.size(); ++i )
3869   {
3870     SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
3871     bool isSatisfy = myAllNodesFlag;
3872
3873     gp_XYZ centerXYZ (0, 0, 0);
3874
3875     while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
3876     {
3877       SMESH_TNodeXYZ aPnt ( aNodeItr->next() );
3878       centerXYZ += aPnt;
3879       isSatisfy = ! myClassifiers[i]->IsOut( aPnt );
3880     }
3881
3882     // Check the center point for volumes MantisBug 0020168
3883     if (isSatisfy &&
3884         myAllNodesFlag &&
3885         myClassifiers[i]->ShapeType() == TopAbs_SOLID)
3886     {
3887       centerXYZ /= elem->NbNodes();
3888       isSatisfy = ! myClassifiers[i]->IsOut( centerXYZ );
3889     }
3890     if ( isSatisfy )
3891       return true;
3892   }
3893
3894   return false;
3895 }
3896
3897 TopAbs_ShapeEnum ElementsOnShape::TClassifier::ShapeType() const
3898 {
3899   return myShape.ShapeType();
3900 }
3901
3902 bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p)
3903 {
3904   return (this->*myIsOutFun)( p );
3905 }
3906
3907 void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double theTol)
3908 {
3909   myShape = theShape;
3910   myTol   = theTol;
3911   switch ( myShape.ShapeType() )
3912   {
3913   case TopAbs_SOLID: {
3914     mySolidClfr.Load(theShape);
3915     myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid;
3916     break;
3917   }
3918   case TopAbs_FACE:  {
3919     Standard_Real u1,u2,v1,v2;
3920     Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
3921     surf->Bounds( u1,u2,v1,v2 );
3922     myProjFace.Init(surf, u1,u2, v1,v2, myTol );
3923     myIsOutFun = & ElementsOnShape::TClassifier::isOutOfFace;
3924     break;
3925   }
3926   case TopAbs_EDGE:  {
3927     Standard_Real u1, u2;
3928     Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge(theShape), u1, u2);
3929     myProjEdge.Init(curve, u1, u2);
3930     myIsOutFun = & ElementsOnShape::TClassifier::isOutOfEdge;
3931     break;
3932   }
3933   case TopAbs_VERTEX:{
3934     myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
3935     myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex;
3936     break;
3937   }
3938   default:
3939     throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier");
3940   }
3941 }
3942
3943 bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p)
3944 {
3945   mySolidClfr.Perform( p, myTol );
3946   return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON );
3947 }
3948
3949 bool ElementsOnShape::TClassifier::isOutOfFace  (const gp_Pnt& p)
3950 {
3951   myProjFace.Perform( p );
3952   if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
3953   {
3954     // check relatively to the face
3955     Quantity_Parameter u, v;
3956     myProjFace.LowerDistanceParameters(u, v);
3957     gp_Pnt2d aProjPnt (u, v);
3958     BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
3959     if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
3960       return false;
3961   }
3962   return true;
3963 }
3964
3965 bool ElementsOnShape::TClassifier::isOutOfEdge  (const gp_Pnt& p)
3966 {
3967   myProjEdge.Perform( p );
3968   return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
3969 }
3970
3971 bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p)
3972 {
3973   return ( myVertexXYZ.Distance( p ) > myTol );
3974 }
3975
3976
3977 TSequenceOfXYZ::TSequenceOfXYZ()
3978 {}
3979
3980 TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n)
3981 {}
3982
3983 TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t)
3984 {}
3985
3986 TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray)
3987 {}
3988
3989 template <class InputIterator>
3990 TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd)
3991 {}
3992
3993 TSequenceOfXYZ::~TSequenceOfXYZ()
3994 {}
3995
3996 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
3997 {
3998   myArray = theSequenceOfXYZ.myArray;
3999   return *this;
4000 }
4001
4002 gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
4003 {
4004   return myArray[n-1];
4005 }
4006
4007 const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
4008 {
4009   return myArray[n-1];
4010 }
4011
4012 void TSequenceOfXYZ::clear()
4013 {
4014   myArray.clear();
4015 }
4016
4017 void TSequenceOfXYZ::reserve(size_type n)
4018 {
4019   myArray.reserve(n);
4020 }
4021
4022 void TSequenceOfXYZ::push_back(const gp_XYZ& v)
4023 {
4024   myArray.push_back(v);
4025 }
4026
4027 TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
4028 {
4029   return myArray.size();
4030 }
4031
4032 TMeshModifTracer::TMeshModifTracer():
4033   myMeshModifTime(0), myMesh(0)
4034 {
4035 }
4036 void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh )
4037 {
4038   if ( theMesh != myMesh )
4039     myMeshModifTime = 0;
4040   myMesh = theMesh;
4041 }
4042 bool TMeshModifTracer::IsMeshModified()
4043 {
4044   bool modified = false;
4045   if ( myMesh )
4046   {
4047     modified = ( myMeshModifTime != myMesh->GetMTime() );
4048     myMeshModifTime = myMesh->GetMTime();
4049   }
4050   return modified;
4051 }