src/Controls/SMESH_Controls.cxx

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