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   for ( int i = 1; i <= aStr.Length(); ++i )
3216   {
3217     char c = aStr.Value( i );
3218     if ( !isdigit( c ) && c != ',' && c != '-' )
3219       aStr.SetValue( i, ' ');
3220   }
3221   aStr.RemoveAll( ' ' );
3222
3223   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
3224   int i = 1;
3225   while ( tmpStr != "" )
3226   {
3227     tmpStr = aStr.Token( ",", i++ );
3228     int aPos = tmpStr.Search( '-' );
3229
3230     if ( aPos == -1 )
3231     {
3232       if ( tmpStr.IsIntegerValue() )
3233         myIds.Add( tmpStr.IntegerValue() );
3234       else
3235         return false;
3236     }
3237     else
3238     {
3239       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
3240       TCollection_AsciiString aMinStr = tmpStr;
3241
3242       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
3243       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
3244
3245       if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) ||
3246            (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
3247         return false;
3248
3249       myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
3250       myMax.Append( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
3251     }
3252   }
3253
3254   return true;
3255 }
3256
3257 //=======================================================================
3258 // name    : GetType
3259 // Purpose : Get type of supported entities
3260 //=======================================================================
3261 SMDSAbs_ElementType RangeOfIds::GetType() const
3262 {
3263   return myType;
3264 }
3265
3266 //=======================================================================
3267 // name    : SetType
3268 // Purpose : Set type of supported entities
3269 //=======================================================================
3270 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
3271 {
3272   myType = theType;
3273 }
3274
3275 //=======================================================================
3276 // name    : IsSatisfy
3277 // Purpose : Verify whether entity satisfies to this rpedicate
3278 //=======================================================================
3279 bool RangeOfIds::IsSatisfy( long theId )
3280 {
3281   if ( !myMesh )
3282     return false;
3283
3284   if ( myType == SMDSAbs_Node )
3285   {
3286     if ( myMesh->FindNode( theId ) == 0 )
3287       return false;
3288   }
3289   else
3290   {
3291     const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3292     if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All ))
3293       return false;
3294   }
3295
3296   if ( myIds.Contains( theId ) )
3297     return true;
3298
3299   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
3300     if ( theId >= myMin( i ) && theId <= myMax( i ) )
3301       return true;
3302
3303   return false;
3304 }
3305
3306 /*
3307   Class       : Comparator
3308   Description : Base class for comparators
3309 */
3310 Comparator::Comparator():
3311   myMargin(0)
3312 {}
3313
3314 Comparator::~Comparator()
3315 {}
3316
3317 void Comparator::SetMesh( const SMDS_Mesh* theMesh )
3318 {
3319   if ( myFunctor )
3320     myFunctor->SetMesh( theMesh );
3321 }
3322
3323 void Comparator::SetMargin( double theValue )
3324 {
3325   myMargin = theValue;
3326 }
3327
3328 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
3329 {
3330   myFunctor = theFunct;
3331 }
3332
3333 SMDSAbs_ElementType Comparator::GetType() const
3334 {
3335   return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
3336 }
3337
3338 double Comparator::GetMargin()
3339 {
3340   return myMargin;
3341 }
3342
3343
3344 /*
3345   Class       : LessThan
3346   Description : Comparator "<"
3347 */
3348 bool LessThan::IsSatisfy( long theId )
3349 {
3350   return myFunctor && myFunctor->GetValue( theId ) < myMargin;
3351 }
3352
3353
3354 /*
3355   Class       : MoreThan
3356   Description : Comparator ">"
3357 */
3358 bool MoreThan::IsSatisfy( long theId )
3359 {
3360   return myFunctor && myFunctor->GetValue( theId ) > myMargin;
3361 }
3362
3363
3364 /*
3365   Class       : EqualTo
3366   Description : Comparator "="
3367 */
3368 EqualTo::EqualTo():
3369   myToler(Precision::Confusion())
3370 {}
3371
3372 bool EqualTo::IsSatisfy( long theId )
3373 {
3374   return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
3375 }
3376
3377 void EqualTo::SetTolerance( double theToler )
3378 {
3379   myToler = theToler;
3380 }
3381
3382 double EqualTo::GetTolerance()
3383 {
3384   return myToler;
3385 }
3386
3387 /*
3388   Class       : LogicalNOT
3389   Description : Logical NOT predicate
3390 */
3391 LogicalNOT::LogicalNOT()
3392 {}
3393
3394 LogicalNOT::~LogicalNOT()
3395 {}
3396
3397 bool LogicalNOT::IsSatisfy( long theId )
3398 {
3399   return myPredicate && !myPredicate->IsSatisfy( theId );
3400 }
3401
3402 void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
3403 {
3404   if ( myPredicate )
3405     myPredicate->SetMesh( theMesh );
3406 }
3407
3408 void LogicalNOT::SetPredicate( PredicatePtr thePred )
3409 {
3410   myPredicate = thePred;
3411 }
3412
3413 SMDSAbs_ElementType LogicalNOT::GetType() const
3414 {
3415   return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
3416 }
3417
3418
3419 /*
3420   Class       : LogicalBinary
3421   Description : Base class for binary logical predicate
3422 */
3423 LogicalBinary::LogicalBinary()
3424 {}
3425
3426 LogicalBinary::~LogicalBinary()
3427 {}
3428
3429 void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
3430 {
3431   if ( myPredicate1 )
3432     myPredicate1->SetMesh( theMesh );
3433
3434   if ( myPredicate2 )
3435     myPredicate2->SetMesh( theMesh );
3436 }
3437
3438 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
3439 {
3440   myPredicate1 = thePredicate;
3441 }
3442
3443 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
3444 {
3445   myPredicate2 = thePredicate;
3446 }
3447
3448 SMDSAbs_ElementType LogicalBinary::GetType() const
3449 {
3450   if ( !myPredicate1 || !myPredicate2 )
3451     return SMDSAbs_All;
3452
3453   SMDSAbs_ElementType aType1 = myPredicate1->GetType();
3454   SMDSAbs_ElementType aType2 = myPredicate2->GetType();
3455
3456   return aType1 == aType2 ? aType1 : SMDSAbs_All;
3457 }
3458
3459
3460 /*
3461   Class       : LogicalAND
3462   Description : Logical AND
3463 */
3464 bool LogicalAND::IsSatisfy( long theId )
3465 {
3466   return
3467     myPredicate1 &&
3468     myPredicate2 &&
3469     myPredicate1->IsSatisfy( theId ) &&
3470     myPredicate2->IsSatisfy( theId );
3471 }
3472
3473
3474 /*
3475   Class       : LogicalOR
3476   Description : Logical OR
3477 */
3478 bool LogicalOR::IsSatisfy( long theId )
3479 {
3480   return
3481     myPredicate1 &&
3482     myPredicate2 &&
3483     (myPredicate1->IsSatisfy( theId ) ||
3484     myPredicate2->IsSatisfy( theId ));
3485 }
3486
3487
3488 /*
3489                               FILTER
3490 */
3491
3492 // #ifdef WITH_TBB
3493 // #include <tbb/parallel_for.h>
3494 // #include <tbb/enumerable_thread_specific.h>
3495
3496 // namespace Parallel
3497 // {
3498 //   typedef tbb::enumerable_thread_specific< TIdSequence > TIdSeq;
3499
3500 //   struct Predicate
3501 //   {
3502 //     const SMDS_Mesh* myMesh;
3503 //     PredicatePtr     myPredicate;
3504 //     TIdSeq &         myOKIds;
3505 //     Predicate( const SMDS_Mesh* m, PredicatePtr p, TIdSeq & ids ):
3506 //       myMesh(m), myPredicate(p->Duplicate()), myOKIds(ids) {}
3507 //     void operator() ( const tbb::blocked_range<size_t>& r ) const
3508 //     {
3509 //       for ( size_t i = r.begin(); i != r.end(); ++i )
3510 //         if ( myPredicate->IsSatisfy( i ))
3511 //           myOKIds.local().push_back();
3512 //     }
3513 //   }
3514 // }
3515 // #endif
3516
3517 Filter::Filter()
3518 {}
3519
3520 Filter::~Filter()
3521 {}
3522
3523 void Filter::SetPredicate( PredicatePtr thePredicate )
3524 {
3525   myPredicate = thePredicate;
3526 }
3527
3528 void Filter::GetElementsId( const SMDS_Mesh* theMesh,
3529                             PredicatePtr     thePredicate,
3530                             TIdSequence&     theSequence )
3531 {
3532   theSequence.clear();
3533
3534   if ( !theMesh || !thePredicate )
3535     return;
3536
3537   thePredicate->SetMesh( theMesh );
3538
3539   SMDS_ElemIteratorPtr elemIt = theMesh->elementsIterator( thePredicate->GetType() );
3540   if ( elemIt ) {
3541     while ( elemIt->more() ) {
3542       const SMDS_MeshElement* anElem = elemIt->next();
3543       long anId = anElem->GetID();
3544       if ( thePredicate->IsSatisfy( anId ) )
3545         theSequence.push_back( anId );
3546     }
3547   }
3548 }
3549
3550 void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
3551                             Filter::TIdSequence& theSequence )
3552 {
3553   GetElementsId(theMesh,myPredicate,theSequence);
3554 }
3555
3556 /*
3557                               ManifoldPart
3558 */
3559
3560 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
3561
3562 /*
3563    Internal class Link
3564 */
3565
3566 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
3567                           SMDS_MeshNode* theNode2 )
3568 {
3569   myNode1 = theNode1;
3570   myNode2 = theNode2;
3571 }
3572
3573 ManifoldPart::Link::~Link()
3574 {
3575   myNode1 = 0;
3576   myNode2 = 0;
3577 }
3578
3579 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
3580 {
3581   if ( myNode1 == theLink.myNode1 &&
3582        myNode2 == theLink.myNode2 )
3583     return true;
3584   else if ( myNode1 == theLink.myNode2 &&
3585             myNode2 == theLink.myNode1 )
3586     return true;
3587   else
3588     return false;
3589 }
3590
3591 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
3592 {
3593   if(myNode1 < x.myNode1) return true;
3594   if(myNode1 == x.myNode1)
3595     if(myNode2 < x.myNode2) return true;
3596   return false;
3597 }
3598
3599 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
3600                             const ManifoldPart::Link& theLink2 )
3601 {
3602   return theLink1.IsEqual( theLink2 );
3603 }
3604
3605 ManifoldPart::ManifoldPart()
3606 {
3607   myMesh = 0;
3608   myAngToler = Precision::Angular();
3609   myIsOnlyManifold = true;
3610 }
3611
3612 ManifoldPart::~ManifoldPart()
3613 {
3614   myMesh = 0;
3615 }
3616
3617 void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
3618 {
3619   myMesh = theMesh;
3620   process();
3621 }
3622
3623 SMDSAbs_ElementType ManifoldPart::GetType() const
3624 { return SMDSAbs_Face; }
3625
3626 bool ManifoldPart::IsSatisfy( long theElementId )
3627 {
3628   return myMapIds.Contains( theElementId );
3629 }
3630
3631 void ManifoldPart::SetAngleTolerance( const double theAngToler )
3632 { myAngToler = theAngToler; }
3633
3634 double ManifoldPart::GetAngleTolerance() const
3635 { return myAngToler; }
3636
3637 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
3638 { myIsOnlyManifold = theIsOnly; }
3639
3640 void ManifoldPart::SetStartElem( const long  theStartId )
3641 { myStartElemId = theStartId; }
3642
3643 bool ManifoldPart::process()
3644 {
3645   myMapIds.Clear();
3646   myMapBadGeomIds.Clear();
3647
3648   myAllFacePtr.clear();
3649   myAllFacePtrIntDMap.clear();
3650   if ( !myMesh )
3651     return false;
3652
3653   // collect all faces into own map
3654   SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
3655   for (; anFaceItr->more(); )
3656   {
3657     SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
3658     myAllFacePtr.push_back( aFacePtr );
3659     myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
3660   }
3661
3662   SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
3663   if ( !aStartFace )
3664     return false;
3665
3666   // the map of non manifold links and bad geometry
3667   TMapOfLink aMapOfNonManifold;
3668   TColStd_MapOfInteger aMapOfTreated;
3669
3670   // begin cycle on faces from start index and run on vector till the end
3671   //  and from begin to start index to cover whole vector
3672   const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
3673   bool isStartTreat = false;
3674   for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
3675   {
3676     if ( fi == aStartIndx )
3677       isStartTreat = true;
3678     // as result next time when fi will be equal to aStartIndx
3679
3680     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
3681     if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
3682       continue;
3683
3684     aMapOfTreated.Add( aFacePtr->GetID() );
3685     TColStd_MapOfInteger aResFaces;
3686     if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
3687                          aMapOfNonManifold, aResFaces ) )
3688       continue;
3689     TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
3690     for ( ; anItr.More(); anItr.Next() )
3691     {
3692       int aFaceId = anItr.Key();
3693       aMapOfTreated.Add( aFaceId );
3694       myMapIds.Add( aFaceId );
3695     }
3696
3697     if ( fi == int( myAllFacePtr.size() - 1 ))
3698       fi = 0;
3699   } // end run on vector of faces
3700   return !myMapIds.IsEmpty();
3701 }
3702
3703 static void getLinks( const SMDS_MeshFace* theFace,
3704                       ManifoldPart::TVectorOfLink& theLinks )
3705 {
3706   int aNbNode = theFace->NbNodes();
3707   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
3708   int i = 1;
3709   SMDS_MeshNode* aNode = 0;
3710   for ( ; aNodeItr->more() && i <= aNbNode; )
3711   {
3712
3713     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
3714     if ( i == 1 )
3715       aNode = aN1;
3716     i++;
3717     SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
3718     i++;
3719     ManifoldPart::Link aLink( aN1, aN2 );
3720     theLinks.push_back( aLink );
3721   }
3722 }
3723
3724 bool ManifoldPart::findConnected
3725                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
3726                   SMDS_MeshFace*                           theStartFace,
3727                   ManifoldPart::TMapOfLink&                theNonManifold,
3728                   TColStd_MapOfInteger&                    theResFaces )
3729 {
3730   theResFaces.Clear();
3731   if ( !theAllFacePtrInt.size() )
3732     return false;
3733
3734   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
3735   {
3736     myMapBadGeomIds.Add( theStartFace->GetID() );
3737     return false;
3738   }
3739
3740   ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
3741   ManifoldPart::TVectorOfLink aSeqOfBoundary;
3742   theResFaces.Add( theStartFace->GetID() );
3743   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
3744
3745   expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3746                  aDMapLinkFace, theNonManifold, theStartFace );
3747
3748   bool isDone = false;
3749   while ( !isDone && aMapOfBoundary.size() != 0 )
3750   {
3751     bool isToReset = false;
3752     ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
3753     for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
3754     {
3755       ManifoldPart::Link aLink = *pLink;
3756       if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
3757         continue;
3758       // each link could be treated only once
3759       aMapToSkip.insert( aLink );
3760
3761       ManifoldPart::TVectorOfFacePtr aFaces;
3762       // find next
3763       if ( myIsOnlyManifold &&
3764            (theNonManifold.find( aLink ) != theNonManifold.end()) )
3765         continue;
3766       else
3767       {
3768         getFacesByLink( aLink, aFaces );
3769         // filter the element to keep only indicated elements
3770         ManifoldPart::TVectorOfFacePtr aFiltered;
3771         ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3772         for ( ; pFace != aFaces.end(); ++pFace )
3773         {
3774           SMDS_MeshFace* aFace = *pFace;
3775           if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
3776             aFiltered.push_back( aFace );
3777         }
3778         aFaces = aFiltered;
3779         if ( aFaces.size() < 2 )  // no neihgbour faces
3780           continue;
3781         else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
3782         {
3783           theNonManifold.insert( aLink );
3784           continue;
3785         }
3786       }
3787
3788       // compare normal with normals of neighbor element
3789       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
3790       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3791       for ( ; pFace != aFaces.end(); ++pFace )
3792       {
3793         SMDS_MeshFace* aNextFace = *pFace;
3794         if ( aPrevFace == aNextFace )
3795           continue;
3796         int anNextFaceID = aNextFace->GetID();
3797         if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
3798          // should not be with non manifold restriction. probably bad topology
3799           continue;
3800         // check if face was treated and skipped
3801         if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
3802              !isInPlane( aPrevFace, aNextFace ) )
3803           continue;
3804         // add new element to connected and extend the boundaries.
3805         theResFaces.Add( anNextFaceID );
3806         expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3807                         aDMapLinkFace, theNonManifold, aNextFace );
3808         isToReset = true;
3809       }
3810     }
3811     isDone = !isToReset;
3812   }
3813
3814   return !theResFaces.IsEmpty();
3815 }
3816
3817 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
3818                               const SMDS_MeshFace* theFace2 )
3819 {
3820   gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
3821   gp_XYZ aNorm2XYZ = getNormale( theFace2 );
3822   if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
3823   {
3824     myMapBadGeomIds.Add( theFace2->GetID() );
3825     return false;
3826   }
3827   if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
3828     return true;
3829
3830   return false;
3831 }
3832
3833 void ManifoldPart::expandBoundary
3834                    ( ManifoldPart::TMapOfLink&            theMapOfBoundary,
3835                      ManifoldPart::TVectorOfLink&         theSeqOfBoundary,
3836                      ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
3837                      ManifoldPart::TMapOfLink&            theNonManifold,
3838                      SMDS_MeshFace*                       theNextFace ) const
3839 {
3840   ManifoldPart::TVectorOfLink aLinks;
3841   getLinks( theNextFace, aLinks );
3842   int aNbLink = (int)aLinks.size();
3843   for ( int i = 0; i < aNbLink; i++ )
3844   {
3845     ManifoldPart::Link aLink = aLinks[ i ];
3846     if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
3847       continue;
3848     if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
3849     {
3850       if ( myIsOnlyManifold )
3851       {
3852         // remove from boundary
3853         theMapOfBoundary.erase( aLink );
3854         ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
3855         for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
3856         {
3857           ManifoldPart::Link aBoundLink = *pLink;
3858           if ( aBoundLink.IsEqual( aLink ) )
3859           {
3860             theSeqOfBoundary.erase( pLink );
3861             break;
3862           }
3863         }
3864       }
3865     }
3866     else
3867     {
3868       theMapOfBoundary.insert( aLink );
3869       theSeqOfBoundary.push_back( aLink );
3870       theDMapLinkFacePtr[ aLink ] = theNextFace;
3871     }
3872   }
3873 }
3874
3875 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
3876                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
3877 {
3878   std::set<SMDS_MeshCell *> aSetOfFaces;
3879   // take all faces that shared first node
3880   SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
3881   for ( ; anItr->more(); )
3882   {
3883     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3884     if ( !aFace )
3885       continue;
3886     aSetOfFaces.insert( aFace );
3887   }
3888   // take all faces that shared second node
3889   anItr = theLink.myNode2->facesIterator();
3890   // find the common part of two sets
3891   for ( ; anItr->more(); )
3892   {
3893     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3894     if ( aSetOfFaces.count( aFace ) )
3895       theFaces.push_back( aFace );
3896   }
3897 }
3898
3899 /*
3900   Class       : BelongToMeshGroup
3901   Description : Verify whether a mesh element is included into a mesh group
3902 */
3903 BelongToMeshGroup::BelongToMeshGroup(): myGroup( 0 )
3904 {
3905 }
3906
3907 void BelongToMeshGroup::SetGroup( SMESHDS_GroupBase* g )
3908 {
3909   myGroup = g;
3910 }
3911
3912 void BelongToMeshGroup::SetStoreName( const std::string& sn )
3913 {
3914   myStoreName = sn;
3915 }
3916
3917 void BelongToMeshGroup::SetMesh( const SMDS_Mesh* theMesh )
3918 {
3919   if ( myGroup && myGroup->GetMesh() != theMesh )
3920   {
3921     myGroup = 0;
3922   }
3923   if ( !myGroup && !myStoreName.empty() )
3924   {
3925     if ( const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh))
3926     {
3927       const std::set<SMESHDS_GroupBase*>& grps = aMesh->GetGroups();
3928       std::set<SMESHDS_GroupBase*>::const_iterator g = grps.begin();
3929       for ( ; g != grps.end() && !myGroup; ++g )
3930         if ( *g && myStoreName == (*g)->GetStoreName() )
3931           myGroup = *g;
3932     }
3933   }
3934   if ( myGroup )
3935   {
3936     myGroup->IsEmpty(); // make GroupOnFilter update its predicate
3937   }
3938 }
3939
3940 bool BelongToMeshGroup::IsSatisfy( long theElementId )
3941 {
3942   return myGroup ? myGroup->Contains( theElementId ) : false;
3943 }
3944
3945 SMDSAbs_ElementType BelongToMeshGroup::GetType() const
3946 {
3947   return myGroup ? myGroup->GetType() : SMDSAbs_All;
3948 }
3949
3950 /*
3951   ElementsOnSurface
3952 */
3953
3954 ElementsOnSurface::ElementsOnSurface()
3955 {
3956   myIds.Clear();
3957   myType = SMDSAbs_All;
3958   mySurf.Nullify();
3959   myToler = Precision::Confusion();
3960   myUseBoundaries = false;
3961 }
3962
3963 ElementsOnSurface::~ElementsOnSurface()
3964 {
3965 }
3966
3967 void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
3968 {
3969   myMeshModifTracer.SetMesh( theMesh );
3970   if ( myMeshModifTracer.IsMeshModified())
3971     process();
3972 }
3973
3974 bool ElementsOnSurface::IsSatisfy( long theElementId )
3975 {
3976   return myIds.Contains( theElementId );
3977 }
3978
3979 SMDSAbs_ElementType ElementsOnSurface::GetType() const
3980 { return myType; }
3981
3982 void ElementsOnSurface::SetTolerance( const double theToler )
3983 {
3984   if ( myToler != theToler )
3985     myIds.Clear();
3986   myToler = theToler;
3987 }
3988
3989 double ElementsOnSurface::GetTolerance() const
3990 { return myToler; }
3991
3992 void ElementsOnSurface::SetUseBoundaries( bool theUse )
3993 {
3994   if ( myUseBoundaries != theUse ) {
3995     myUseBoundaries = theUse;
3996     SetSurface( mySurf, myType );
3997   }
3998 }
3999
4000 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
4001                                     const SMDSAbs_ElementType theType )
4002 {
4003   myIds.Clear();
4004   myType = theType;
4005   mySurf.Nullify();
4006   if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
4007     return;
4008   mySurf = TopoDS::Face( theShape );
4009   BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
4010   Standard_Real
4011     u1 = SA.FirstUParameter(),
4012     u2 = SA.LastUParameter(),
4013     v1 = SA.FirstVParameter(),
4014     v2 = SA.LastVParameter();
4015   Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
4016   myProjector.Init( surf, u1,u2, v1,v2 );
4017   process();
4018 }
4019
4020 void ElementsOnSurface::process()
4021 {
4022   myIds.Clear();
4023   if ( mySurf.IsNull() )
4024     return;
4025
4026   if ( !myMeshModifTracer.GetMesh() )
4027     return;
4028
4029   myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
4030
4031   SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType );
4032   for(; anIter->more(); )
4033     process( anIter->next() );
4034 }
4035
4036 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
4037 {
4038   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
4039   bool isSatisfy = true;
4040   for ( ; aNodeItr->more(); )
4041   {
4042     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
4043     if ( !isOnSurface( aNode ) )
4044     {
4045       isSatisfy = false;
4046       break;
4047     }
4048   }
4049   if ( isSatisfy )
4050     myIds.Add( theElemPtr->GetID() );
4051 }
4052
4053 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
4054 {
4055   if ( mySurf.IsNull() )
4056     return false;
4057
4058   gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
4059   //  double aToler2 = myToler * myToler;
4060 //   if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
4061 //   {
4062 //     gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
4063 //     if ( aPln.SquareDistance( aPnt ) > aToler2 )
4064 //       return false;
4065 //   }
4066 //   else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
4067 //   {
4068 //     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
4069 //     double aRad = aCyl.Radius();
4070 //     gp_Ax3 anAxis = aCyl.Position();
4071 //     gp_XYZ aLoc = aCyl.Location().XYZ();
4072 //     double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4073 //     double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4074 //     if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
4075 //       return false;
4076 //   }
4077 //   else
4078 //     return false;
4079   myProjector.Perform( aPnt );
4080   bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
4081
4082   return isOn;
4083 }
4084
4085
4086 /*
4087   ElementsOnShape
4088 */
4089
4090 ElementsOnShape::ElementsOnShape()
4091   : //myMesh(0),
4092     myType(SMDSAbs_All),
4093     myToler(Precision::Confusion()),
4094     myAllNodesFlag(false)
4095 {
4096 }
4097
4098 ElementsOnShape::~ElementsOnShape()
4099 {
4100   clearClassifiers();
4101 }
4102
4103 SMDSAbs_ElementType ElementsOnShape::GetType() const
4104 {
4105   return myType;
4106 }
4107
4108 void ElementsOnShape::SetTolerance (const double theToler)
4109 {
4110   if (myToler != theToler) {
4111     myToler = theToler;
4112     SetShape(myShape, myType);
4113   }
4114 }
4115
4116 double ElementsOnShape::GetTolerance() const
4117 {
4118   return myToler;
4119 }
4120
4121 void ElementsOnShape::SetAllNodes (bool theAllNodes)
4122 {
4123   myAllNodesFlag = theAllNodes;
4124 }
4125
4126 void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
4127 {
4128   myMeshModifTracer.SetMesh( theMesh );
4129   if ( myMeshModifTracer.IsMeshModified())
4130   {
4131     size_t nbNodes = theMesh ? theMesh->NbNodes() : 0;
4132     if ( myNodeIsChecked.size() == nbNodes )
4133     {
4134       std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4135     }
4136     else
4137     {
4138       SMESHUtils::FreeVector( myNodeIsChecked );
4139       SMESHUtils::FreeVector( myNodeIsOut );
4140       myNodeIsChecked.resize( nbNodes, false );
4141       myNodeIsOut.resize( nbNodes );
4142     }
4143   }
4144 }
4145
4146 bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut )
4147 {
4148   if ( n->GetID() >= (int) myNodeIsChecked.size() ||
4149        !myNodeIsChecked[ n->GetID() ])
4150     return false;
4151
4152   isOut = myNodeIsOut[ n->GetID() ];
4153   return true;
4154 }
4155
4156 void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool  isOut )
4157 {
4158   if ( n->GetID() < (int) myNodeIsChecked.size() )
4159   {
4160     myNodeIsChecked[ n->GetID() ] = true;
4161     myNodeIsOut    [ n->GetID() ] = isOut;
4162   }
4163 }
4164
4165 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
4166                                 const SMDSAbs_ElementType theType)
4167 {
4168   myType  = theType;
4169   myShape = theShape;
4170   if ( myShape.IsNull() ) return;
4171
4172   TopTools_IndexedMapOfShape shapesMap;
4173   TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
4174   TopExp_Explorer sub;
4175   for ( int i = 0; i < 4; ++i )
4176   {
4177     if ( shapesMap.IsEmpty() )
4178       for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() )
4179         shapesMap.Add( sub.Current() );
4180     if ( i > 0 )
4181       for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() )
4182         shapesMap.Add( sub.Current() );
4183   }
4184
4185   clearClassifiers();
4186   myClassifiers.resize( shapesMap.Extent() );
4187   for ( int i = 0; i < shapesMap.Extent(); ++i )
4188     myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler );
4189
4190   if ( theType == SMDSAbs_Node )
4191   {
4192     SMESHUtils::FreeVector( myNodeIsChecked );
4193     SMESHUtils::FreeVector( myNodeIsOut );
4194   }
4195   else
4196   {
4197     std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4198   }
4199 }
4200
4201 void ElementsOnShape::clearClassifiers()
4202 {
4203   for ( size_t i = 0; i < myClassifiers.size(); ++i )
4204     delete myClassifiers[ i ];
4205   myClassifiers.clear();
4206 }
4207
4208 bool ElementsOnShape::IsSatisfy (long elemId)
4209 {
4210   const SMDS_Mesh*        mesh = myMeshModifTracer.GetMesh();
4211   const SMDS_MeshElement* elem =
4212     ( myType == SMDSAbs_Node ? mesh->FindNode( elemId ) : mesh->FindElement( elemId ));
4213   if ( !elem || myClassifiers.empty() )
4214     return false;
4215
4216   bool isSatisfy = myAllNodesFlag, isNodeOut;
4217
4218   gp_XYZ centerXYZ (0, 0, 0);
4219
4220   SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
4221   while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
4222   {
4223     SMESH_TNodeXYZ aPnt( aNodeItr->next() );
4224     centerXYZ += aPnt;
4225
4226     isNodeOut = true;
4227     if ( !getNodeIsOut( aPnt._node, isNodeOut ))
4228     {
4229       for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i )
4230         isNodeOut = myClassifiers[i]->IsOut( aPnt );
4231
4232       setNodeIsOut( aPnt._node, isNodeOut );
4233     }
4234     isSatisfy = !isNodeOut;
4235   }
4236
4237   // Check the center point for volumes MantisBug 0020168
4238   if (isSatisfy &&
4239       myAllNodesFlag &&
4240       myClassifiers[0]->ShapeType() == TopAbs_SOLID)
4241   {
4242     centerXYZ /= elem->NbNodes();
4243     isSatisfy = false;
4244     for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i )
4245       isSatisfy = ! myClassifiers[i]->IsOut( centerXYZ );
4246   }
4247
4248   return isSatisfy;
4249 }
4250
4251 TopAbs_ShapeEnum ElementsOnShape::TClassifier::ShapeType() const
4252 {
4253   return myShape.ShapeType();
4254 }
4255
4256 bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p)
4257 {
4258   return (this->*myIsOutFun)( p );
4259 }
4260
4261 void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double theTol)
4262 {
4263   myShape = theShape;
4264   myTol   = theTol;
4265   switch ( myShape.ShapeType() )
4266   {
4267   case TopAbs_SOLID: {
4268     if ( isBox( theShape ))
4269     {
4270       myIsOutFun = & ElementsOnShape::TClassifier::isOutOfBox;
4271     }
4272     else
4273     {
4274       mySolidClfr.Load(theShape);
4275       myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid;
4276     }
4277     break;
4278   }
4279   case TopAbs_FACE:  {
4280     Standard_Real u1,u2,v1,v2;
4281     Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
4282     surf->Bounds( u1,u2,v1,v2 );
4283     myProjFace.Init(surf, u1,u2, v1,v2, myTol );
4284     myIsOutFun = & ElementsOnShape::TClassifier::isOutOfFace;
4285     break;
4286   }
4287   case TopAbs_EDGE:  {
4288     Standard_Real u1, u2;
4289     Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge(theShape), u1, u2);
4290     myProjEdge.Init(curve, u1, u2);
4291     myIsOutFun = & ElementsOnShape::TClassifier::isOutOfEdge;
4292     break;
4293   }
4294   case TopAbs_VERTEX:{
4295     myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
4296     myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex;
4297     break;
4298   }
4299   default:
4300     throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier");
4301   }
4302 }
4303
4304 bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p)
4305 {
4306   mySolidClfr.Perform( p, myTol );
4307   return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON );
4308 }
4309
4310 bool ElementsOnShape::TClassifier::isOutOfBox (const gp_Pnt& p)
4311 {
4312   return myBox.IsOut( p.XYZ() );
4313 }
4314
4315 bool ElementsOnShape::TClassifier::isOutOfFace  (const gp_Pnt& p)
4316 {
4317   myProjFace.Perform( p );
4318   if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
4319   {
4320     // check relatively to the face
4321     Quantity_Parameter u, v;
4322     myProjFace.LowerDistanceParameters(u, v);
4323     gp_Pnt2d aProjPnt (u, v);
4324     BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
4325     if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
4326       return false;
4327   }
4328   return true;
4329 }
4330
4331 bool ElementsOnShape::TClassifier::isOutOfEdge  (const gp_Pnt& p)
4332 {
4333   myProjEdge.Perform( p );
4334   return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
4335 }
4336
4337 bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p)
4338 {
4339   return ( myVertexXYZ.Distance( p ) > myTol );
4340 }
4341
4342 bool ElementsOnShape::TClassifier::isBox (const TopoDS_Shape& theShape)
4343 {
4344   TopTools_IndexedMapOfShape vMap;
4345   TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
4346   if ( vMap.Extent() != 8 )
4347     return false;
4348
4349   myBox.Clear();
4350   for ( int i = 1; i <= 8; ++i )
4351     myBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))).XYZ() );
4352
4353   gp_XYZ pMin = myBox.CornerMin(), pMax = myBox.CornerMax();
4354   for ( int i = 1; i <= 8; ++i )
4355   {
4356     gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i )));
4357     for ( int iC = 1; iC <= 3; ++ iC )
4358     {
4359       double d1 = Abs( pMin.Coord( iC ) - p.Coord( iC ));
4360       double d2 = Abs( pMax.Coord( iC ) - p.Coord( iC ));
4361       if ( Min( d1, d2 ) > myTol )
4362         return false;
4363     }
4364   }
4365   myBox.Enlarge( myTol );
4366   return true;
4367 }
4368
4369
4370 /*
4371   Class       : BelongToGeom
4372   Description : Predicate for verifying whether entity belongs to
4373                 specified geometrical support
4374 */
4375
4376 BelongToGeom::BelongToGeom()
4377   : myMeshDS(NULL),
4378     myType(SMDSAbs_All),
4379     myIsSubshape(false),
4380     myTolerance(Precision::Confusion())
4381 {}
4382
4383 void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
4384 {
4385   myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
4386   init();
4387 }
4388
4389 void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
4390 {
4391   myShape = theShape;
4392   init();
4393 }
4394
4395 static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
4396                         const TopoDS_Shape& theShape)
4397 {
4398   if (theMap.Contains(theShape)) return true;
4399
4400   if (theShape.ShapeType() == TopAbs_COMPOUND ||
4401       theShape.ShapeType() == TopAbs_COMPSOLID)
4402   {
4403     TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
4404     for (; anIt.More(); anIt.Next())
4405     {
4406       if (!IsSubShape(theMap, anIt.Value())) {
4407         return false;
4408       }
4409     }
4410     return true;
4411   }
4412
4413   return false;
4414 }
4415
4416 void BelongToGeom::init()
4417 {
4418   if (!myMeshDS || myShape.IsNull()) return;
4419
4420   // is sub-shape of main shape?
4421   TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
4422   if (aMainShape.IsNull()) {
4423     myIsSubshape = false;
4424   }
4425   else {
4426     TopTools_IndexedMapOfShape aMap;
4427     TopExp::MapShapes(aMainShape, aMap);
4428     myIsSubshape = IsSubShape(aMap, myShape);
4429   }
4430
4431   //if (!myIsSubshape) // to be always ready to check an element not bound to geometry
4432   {
4433     myElementsOnShapePtr.reset(new ElementsOnShape());
4434     myElementsOnShapePtr->SetTolerance(myTolerance);
4435     myElementsOnShapePtr->SetAllNodes(true); // "belong", while false means "lays on"
4436     myElementsOnShapePtr->SetMesh(myMeshDS);
4437     myElementsOnShapePtr->SetShape(myShape, myType);
4438   }
4439 }
4440
4441 static bool IsContains( const SMESHDS_Mesh*     theMeshDS,
4442                         const TopoDS_Shape&     theShape,
4443                         const SMDS_MeshElement* theElem,
4444                         TopAbs_ShapeEnum        theFindShapeEnum,
4445                         TopAbs_ShapeEnum        theAvoidShapeEnum = TopAbs_SHAPE )
4446 {
4447   TopExp_Explorer anExp( theShape,theFindShapeEnum,theAvoidShapeEnum );
4448
4449   while( anExp.More() )
4450   {
4451     const TopoDS_Shape& aShape = anExp.Current();
4452     if( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ) ){
4453       if( aSubMesh->Contains( theElem ) )
4454         return true;
4455     }
4456     anExp.Next();
4457   }
4458   return false;
4459 }
4460
4461 bool BelongToGeom::IsSatisfy (long theId)
4462 {
4463   if (myMeshDS == 0 || myShape.IsNull())
4464     return false;
4465
4466   if (!myIsSubshape)
4467   {
4468     return myElementsOnShapePtr->IsSatisfy(theId);
4469   }
4470
4471   // Case of submesh
4472   if (myType == SMDSAbs_Node)
4473   {
4474     if( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ) )
4475     {
4476       if ( aNode->getshapeId() < 1 )
4477         return myElementsOnShapePtr->IsSatisfy(theId);
4478
4479       const SMDS_PositionPtr& aPosition = aNode->GetPosition();
4480       SMDS_TypeOfPosition aTypeOfPosition = aPosition->GetTypeOfPosition();
4481       switch( aTypeOfPosition )
4482       {
4483       case SMDS_TOP_VERTEX : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_VERTEX ));
4484       case SMDS_TOP_EDGE   : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_EDGE ));
4485       case SMDS_TOP_FACE   : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_FACE ));
4486       case SMDS_TOP_3DSPACE: return ( IsContains( myMeshDS,myShape,aNode,TopAbs_SOLID ) ||
4487                                       IsContains( myMeshDS,myShape,aNode,TopAbs_SHELL ));
4488       default:;
4489       }
4490     }
4491   }
4492   else
4493   {
4494     if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
4495     {
4496       if ( anElem->getshapeId() < 1 )
4497         return myElementsOnShapePtr->IsSatisfy(theId);
4498
4499       if( myType == SMDSAbs_All )
4500       {
4501         return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ) ||
4502                  IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ) ||
4503                  IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
4504                  IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
4505       }
4506       else if( myType == anElem->GetType() )
4507       {
4508         switch( myType )
4509         {
4510         case SMDSAbs_Edge  : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ));
4511         case SMDSAbs_Face  : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ));
4512         case SMDSAbs_Volume: return ( IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
4513                                       IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
4514         default:;
4515         }
4516       }
4517     }
4518   }
4519
4520   return false;
4521 }
4522
4523 void BelongToGeom::SetType (SMDSAbs_ElementType theType)
4524 {
4525   myType = theType;
4526   init();
4527 }
4528
4529 SMDSAbs_ElementType BelongToGeom::GetType() const
4530 {
4531   return myType;
4532 }
4533
4534 TopoDS_Shape BelongToGeom::GetShape()
4535 {
4536   return myShape;
4537 }
4538
4539 const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const
4540 {
4541   return myMeshDS;
4542 }
4543
4544 void BelongToGeom::SetTolerance (double theTolerance)
4545 {
4546   myTolerance = theTolerance;
4547   if (!myIsSubshape)
4548     init();
4549 }
4550
4551 double BelongToGeom::GetTolerance()
4552 {
4553   return myTolerance;
4554 }
4555
4556 /*
4557   Class       : LyingOnGeom
4558   Description : Predicate for verifying whether entiy lying or partially lying on
4559                 specified geometrical support
4560 */
4561
4562 LyingOnGeom::LyingOnGeom()
4563   : myMeshDS(NULL),
4564     myType(SMDSAbs_All),
4565     myIsSubshape(false),
4566     myTolerance(Precision::Confusion())
4567 {}
4568
4569 void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
4570 {
4571   myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
4572   init();
4573 }
4574
4575 void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
4576 {
4577   myShape = theShape;
4578   init();
4579 }
4580
4581 void LyingOnGeom::init()
4582 {
4583   if (!myMeshDS || myShape.IsNull()) return;
4584
4585   // is sub-shape of main shape?
4586   TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
4587   if (aMainShape.IsNull()) {
4588     myIsSubshape = false;
4589   }
4590   else {
4591     myIsSubshape = myMeshDS->IsGroupOfSubShapes( myShape );
4592   }
4593
4594   if (myIsSubshape)
4595   {
4596     TopTools_IndexedMapOfShape shapes;
4597     TopExp::MapShapes( myShape, shapes );
4598     mySubShapesIDs.Clear();
4599     for ( int i = 1; i <= shapes.Extent(); ++i )
4600     {
4601       int subID = myMeshDS->ShapeToIndex( shapes( i ));
4602       if ( subID > 0 )
4603         mySubShapesIDs.Add( subID );
4604     }
4605   }
4606   else
4607   {
4608     myElementsOnShapePtr.reset(new ElementsOnShape());
4609     myElementsOnShapePtr->SetTolerance(myTolerance);
4610     myElementsOnShapePtr->SetAllNodes(false); // lays on, while true means "belong"
4611     myElementsOnShapePtr->SetMesh(myMeshDS);
4612     myElementsOnShapePtr->SetShape(myShape, myType);
4613   }
4614 }
4615
4616 bool LyingOnGeom::IsSatisfy( long theId )
4617 {
4618   if ( myMeshDS == 0 || myShape.IsNull() )
4619     return false;
4620
4621   if (!myIsSubshape)
4622   {
4623     return myElementsOnShapePtr->IsSatisfy(theId);
4624   }
4625
4626   // Case of sub-mesh
4627
4628   const SMDS_MeshElement* elem =
4629     ( myType == SMDSAbs_Node ) ? myMeshDS->FindNode( theId ) : myMeshDS->FindElement( theId );
4630
4631   if ( mySubShapesIDs.Contains( elem->getshapeId() ))
4632     return true;
4633
4634   if ( elem->GetType() != SMDSAbs_Node )
4635   {
4636     SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator();
4637     while ( nodeItr->more() )
4638     {
4639       const SMDS_MeshElement* aNode = nodeItr->next();
4640       if ( mySubShapesIDs.Contains( aNode->getshapeId() ))
4641         return true;
4642     }
4643   }
4644
4645   return false;
4646 }
4647
4648 void LyingOnGeom::SetType( SMDSAbs_ElementType theType )
4649 {
4650   myType = theType;
4651   init();
4652 }
4653
4654 SMDSAbs_ElementType LyingOnGeom::GetType() const
4655 {
4656   return myType;
4657 }
4658
4659 TopoDS_Shape LyingOnGeom::GetShape()
4660 {
4661   return myShape;
4662 }
4663
4664 const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const
4665 {
4666   return myMeshDS;
4667 }
4668
4669 void LyingOnGeom::SetTolerance (double theTolerance)
4670 {
4671   myTolerance = theTolerance;
4672   if (!myIsSubshape)
4673     init();
4674 }
4675
4676 double LyingOnGeom::GetTolerance()
4677 {
4678   return myTolerance;
4679 }
4680
4681 bool LyingOnGeom::Contains( const SMESHDS_Mesh*     theMeshDS,
4682                             const TopoDS_Shape&     theShape,
4683                             const SMDS_MeshElement* theElem,
4684                             TopAbs_ShapeEnum        theFindShapeEnum,
4685                             TopAbs_ShapeEnum        theAvoidShapeEnum )
4686 {
4687   // if (IsContains(theMeshDS, theShape, theElem, theFindShapeEnum, theAvoidShapeEnum))
4688   //   return true;
4689
4690   // TopTools_MapOfShape aSubShapes;
4691   // TopExp_Explorer exp( theShape, theFindShapeEnum, theAvoidShapeEnum );
4692   // for ( ; exp.More(); exp.Next() )
4693   // {
4694   //   const TopoDS_Shape& aShape = exp.Current();
4695   //   if ( !aSubShapes.Add( aShape )) continue;
4696
4697   //   if ( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ))
4698   //   {
4699   //     if ( aSubMesh->Contains( theElem ))
4700   //       return true;
4701
4702   //     SMDS_ElemIteratorPtr nodeItr = theElem->nodesIterator();
4703   //     while ( nodeItr->more() )
4704   //     {
4705   //       const SMDS_MeshElement* aNode = nodeItr->next();
4706   //       if ( aSubMesh->Contains( aNode ))
4707   //         return true;
4708   //     }
4709   //   }
4710   // }
4711   return false;
4712 }
4713
4714 TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0)
4715 {}
4716
4717 TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n), myElem(0)
4718 {}
4719
4720 TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t), myElem(0)
4721 {}
4722
4723 TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray), myElem(theSequenceOfXYZ.myElem)
4724 {}
4725
4726 template <class InputIterator>
4727 TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd), myElem(0)
4728 {}
4729
4730 TSequenceOfXYZ::~TSequenceOfXYZ()
4731 {}
4732
4733 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
4734 {
4735   myArray = theSequenceOfXYZ.myArray;
4736   myElem  = theSequenceOfXYZ.myElem;
4737   return *this;
4738 }
4739
4740 gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
4741 {
4742   return myArray[n-1];
4743 }
4744
4745 const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
4746 {
4747   return myArray[n-1];
4748 }
4749
4750 void TSequenceOfXYZ::clear()
4751 {
4752   myArray.clear();
4753 }
4754
4755 void TSequenceOfXYZ::reserve(size_type n)
4756 {
4757   myArray.reserve(n);
4758 }
4759
4760 void TSequenceOfXYZ::push_back(const gp_XYZ& v)
4761 {
4762   myArray.push_back(v);
4763 }
4764
4765 TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
4766 {
4767   return myArray.size();
4768 }
4769
4770 SMDSAbs_EntityType TSequenceOfXYZ::getElementEntity() const
4771 {
4772   return myElem ? myElem->GetEntityType() : SMDSEntity_Last;
4773 }
4774
4775 TMeshModifTracer::TMeshModifTracer():
4776   myMeshModifTime(0), myMesh(0)
4777 {
4778 }
4779 void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh )
4780 {
4781   if ( theMesh != myMesh )
4782     myMeshModifTime = 0;
4783   myMesh = theMesh;
4784 }
4785 bool TMeshModifTracer::IsMeshModified()
4786 {
4787   bool modified = false;
4788   if ( myMesh )
4789   {
4790     modified = ( myMeshModifTime != myMesh->GetMTime() );
4791     myMeshModifTime = myMesh->GetMTime();
4792   }
4793   return modified;
4794 }