src/Controls/SMESH_Controls.cxx

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