src/Controls/SMESH_Controls.cxx

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