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