src/Controls/SMESH_Controls.cxx

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