src/Controls/SMESH_Controls.cxx

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