src/Controls/SMESH_Controls.cxx

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