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::SetMesh( const SMDS_Mesh* theMesh )
2440 {
2441   myMeshModifTracer.SetMesh( theMesh );
2442   if ( myMeshModifTracer.IsMeshModified() )
2443   {
2444     TIDSortedNodeSet nodesToCheck;
2445     SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
2446     while ( nIt->more() )
2447       nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
2448
2449     std::list< std::list< const SMDS_MeshNode*> > nodeGroups;
2450     SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
2451
2452     myCoincidentIDs.Clear();
2453     std::list< std::list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
2454     for ( ; groupIt != nodeGroups.end(); ++groupIt )
2455     {
2456       std::list< const SMDS_MeshNode*>& coincNodes = *groupIt;
2457       std::list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
2458       for ( ; n != coincNodes.end(); ++n )
2459         myCoincidentIDs.Add( (*n)->GetID() );
2460     }
2461   }
2462 }
2463
2464 //================================================================================
2465 /*
2466   Class       : CoincidentElements
2467   Description : Predicate of Coincident Elements
2468   Note        : This class is suitable only for visualization of Coincident Elements
2469 */
2470 //================================================================================
2471
2472 CoincidentElements::CoincidentElements()
2473 {
2474   myMesh = 0;
2475 }
2476
2477 void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
2478 {
2479   myMesh = theMesh;
2480 }
2481
2482 bool CoincidentElements::IsSatisfy( long theElementId )
2483 {
2484   if ( !myMesh ) return false;
2485
2486   if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
2487   {
2488     if ( e->GetType() != GetType() ) return false;
2489     std::set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
2490     const int nbNodes = e->NbNodes();
2491     SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
2492     while ( invIt->more() )
2493     {
2494       const SMDS_MeshElement* e2 = invIt->next();
2495       if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
2496
2497       bool sameNodes = true;
2498       for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
2499         sameNodes = ( elemNodes.count( e2->GetNode( i )));
2500       if ( sameNodes )
2501         return true;
2502     }
2503   }
2504   return false;
2505 }
2506
2507 SMDSAbs_ElementType CoincidentElements1D::GetType() const
2508 {
2509   return SMDSAbs_Edge;
2510 }
2511 SMDSAbs_ElementType CoincidentElements2D::GetType() const
2512 {
2513   return SMDSAbs_Face;
2514 }
2515 SMDSAbs_ElementType CoincidentElements3D::GetType() const
2516 {
2517   return SMDSAbs_Volume;
2518 }
2519
2520
2521 //================================================================================
2522 /*
2523   Class       : FreeBorders
2524   Description : Predicate for free borders
2525 */
2526 //================================================================================
2527
2528 FreeBorders::FreeBorders()
2529 {
2530   myMesh = 0;
2531 }
2532
2533 void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
2534 {
2535   myMesh = theMesh;
2536 }
2537
2538 bool FreeBorders::IsSatisfy( long theId )
2539 {
2540   return getNbMultiConnection( myMesh, theId ) == 1;
2541 }
2542
2543 SMDSAbs_ElementType FreeBorders::GetType() const
2544 {
2545   return SMDSAbs_Edge;
2546 }
2547
2548
2549 //================================================================================
2550 /*
2551   Class       : FreeEdges
2552   Description : Predicate for free Edges
2553 */
2554 //================================================================================
2555
2556 FreeEdges::FreeEdges()
2557 {
2558   myMesh = 0;
2559 }
2560
2561 void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
2562 {
2563   myMesh = theMesh;
2564 }
2565
2566 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId  )
2567 {
2568   TColStd_MapOfInteger aMap;
2569   for ( int i = 0; i < 2; i++ )
2570   {
2571     SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face);
2572     while( anElemIter->more() )
2573     {
2574       if ( const SMDS_MeshElement* anElem = anElemIter->next())
2575       {
2576         const int anId = anElem->GetID();
2577         if ( anId != theFaceId && !aMap.Add( anId ))
2578           return false;
2579       }
2580     }
2581   }
2582   return true;
2583 }
2584
2585 bool FreeEdges::IsSatisfy( long theId )
2586 {
2587   if ( myMesh == 0 )
2588     return false;
2589
2590   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2591   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
2592     return false;
2593
2594   SMDS_NodeIteratorPtr anIter = aFace->interlacedNodesIterator();
2595   if ( !anIter )
2596     return false;
2597
2598   int i = 0, nbNodes = aFace->NbNodes();
2599   std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
2600   while( anIter->more() )
2601     if ( ! ( aNodes[ i++ ] = anIter->next() ))
2602       return false;
2603   aNodes[ nbNodes ] = aNodes[ 0 ];
2604
2605   for ( i = 0; i < nbNodes; i++ )
2606     if ( IsFreeEdge( &aNodes[ i ], theId ) )
2607       return true;
2608
2609   return false;
2610 }
2611
2612 SMDSAbs_ElementType FreeEdges::GetType() const
2613 {
2614   return SMDSAbs_Face;
2615 }
2616
2617 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
2618   myElemId(theElemId)
2619 {
2620   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
2621   if(thePntId1 > thePntId2){
2622     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
2623   }
2624 }
2625
2626 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
2627   if(myPntId[0] < x.myPntId[0]) return true;
2628   if(myPntId[0] == x.myPntId[0])
2629     if(myPntId[1] < x.myPntId[1]) return true;
2630   return false;
2631 }
2632
2633 inline void UpdateBorders(const FreeEdges::Border& theBorder,
2634                           FreeEdges::TBorders& theRegistry,
2635                           FreeEdges::TBorders& theContainer)
2636 {
2637   if(theRegistry.find(theBorder) == theRegistry.end()){
2638     theRegistry.insert(theBorder);
2639     theContainer.insert(theBorder);
2640   }else{
2641     theContainer.erase(theBorder);
2642   }
2643 }
2644
2645 void FreeEdges::GetBoreders(TBorders& theBorders)
2646 {
2647   TBorders aRegistry;
2648   for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
2649   {
2650     const SMDS_MeshFace* anElem = anIter->next();
2651     long anElemId = anElem->GetID();
2652     SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
2653     if ( !aNodesIter->more() ) continue;
2654     long aNodeId[2] = {0,0};
2655     aNodeId[0] = anElem->GetNode( anElem->NbNodes()-1 )->GetID();
2656     for ( ; aNodesIter->more(); )
2657     {
2658       aNodeId[1] = aNodesIter->next()->GetID();
2659       Border aBorder( anElemId, aNodeId[0], aNodeId[1] );
2660       UpdateBorders( aBorder, aRegistry, theBorders );
2661       aNodeId[0] = aNodeId[1];
2662     }
2663   }
2664 }
2665
2666 //================================================================================
2667 /*
2668   Class       : FreeNodes
2669   Description : Predicate for free nodes
2670 */
2671 //================================================================================
2672
2673 FreeNodes::FreeNodes()
2674 {
2675   myMesh = 0;
2676 }
2677
2678 void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
2679 {
2680   myMesh = theMesh;
2681 }
2682
2683 bool FreeNodes::IsSatisfy( long theNodeId )
2684 {
2685   const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
2686   if (!aNode)
2687     return false;
2688
2689   return (aNode->NbInverseElements() < 1);
2690 }
2691
2692 SMDSAbs_ElementType FreeNodes::GetType() const
2693 {
2694   return SMDSAbs_Node;
2695 }
2696
2697
2698 //================================================================================
2699 /*
2700   Class       : FreeFaces
2701   Description : Predicate for free faces
2702 */
2703 //================================================================================
2704
2705 FreeFaces::FreeFaces()
2706 {
2707   myMesh = 0;
2708 }
2709
2710 void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
2711 {
2712   myMesh = theMesh;
2713 }
2714
2715 bool FreeFaces::IsSatisfy( long theId )
2716 {
2717   if (!myMesh) return false;
2718   // check that faces nodes refers to less than two common volumes
2719   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2720   if ( !aFace || aFace->GetType() != SMDSAbs_Face )
2721     return false;
2722
2723   int nbNode = aFace->NbNodes();
2724
2725   // collect volumes to check that number of volumes with count equal nbNode not less than 2
2726   typedef std::map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
2727   typedef std::map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
2728   TMapOfVolume mapOfVol;
2729
2730   SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
2731   while ( nodeItr->more() )
2732   {
2733     const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
2734     if ( !aNode ) continue;
2735     SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
2736     while ( volItr->more() )
2737     {
2738       SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
2739       TItrMapOfVolume    itr = mapOfVol.insert( std::make_pair( aVol, 0 )).first;
2740       (*itr).second++;
2741     }
2742   }
2743   int nbVol = 0;
2744   TItrMapOfVolume volItr = mapOfVol.begin();
2745   TItrMapOfVolume volEnd = mapOfVol.end();
2746   for ( ; volItr != volEnd; ++volItr )
2747     if ( (*volItr).second >= nbNode )
2748        nbVol++;
2749   // face is not free if number of volumes constructed on their nodes more than one
2750   return (nbVol < 2);
2751 }
2752
2753 SMDSAbs_ElementType FreeFaces::GetType() const
2754 {
2755   return SMDSAbs_Face;
2756 }
2757
2758 //================================================================================
2759 /*
2760   Class       : LinearOrQuadratic
2761   Description : Predicate to verify whether a mesh element is linear
2762 */
2763 //================================================================================
2764
2765 LinearOrQuadratic::LinearOrQuadratic()
2766 {
2767   myMesh = 0;
2768 }
2769
2770 void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
2771 {
2772   myMesh = theMesh;
2773 }
2774
2775 bool LinearOrQuadratic::IsSatisfy( long theId )
2776 {
2777   if (!myMesh) return false;
2778   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2779   if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
2780     return false;
2781   return (!anElem->IsQuadratic());
2782 }
2783
2784 void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
2785 {
2786   myType = theType;
2787 }
2788
2789 SMDSAbs_ElementType LinearOrQuadratic::GetType() const
2790 {
2791   return myType;
2792 }
2793
2794 //================================================================================
2795 /*
2796   Class       : GroupColor
2797   Description : Functor for check color of group to which mesh element belongs to
2798 */
2799 //================================================================================
2800
2801 GroupColor::GroupColor()
2802 {
2803 }
2804
2805 bool GroupColor::IsSatisfy( long theId )
2806 {
2807   return myIDs.count( theId );
2808 }
2809
2810 void GroupColor::SetType( SMDSAbs_ElementType theType )
2811 {
2812   myType = theType;
2813 }
2814
2815 SMDSAbs_ElementType GroupColor::GetType() const
2816 {
2817   return myType;
2818 }
2819
2820 static bool isEqual( const Quantity_Color& theColor1,
2821                      const Quantity_Color& theColor2 )
2822 {
2823   // tolerance to compare colors
2824   const double tol = 5*1e-3;
2825   return ( fabs( theColor1.Red()   - theColor2.Red() )   < tol &&
2826            fabs( theColor1.Green() - theColor2.Green() ) < tol &&
2827            fabs( theColor1.Blue()  - theColor2.Blue() )  < tol );
2828 }
2829
2830 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
2831 {
2832   myIDs.clear();
2833
2834   const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
2835   if ( !aMesh )
2836     return;
2837
2838   int nbGrp = aMesh->GetNbGroups();
2839   if ( !nbGrp )
2840     return;
2841
2842   // iterates on groups and find necessary elements ids
2843   const std::set<SMESHDS_GroupBase*>&       aGroups = aMesh->GetGroups();
2844   std::set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
2845   for (; GrIt != aGroups.end(); GrIt++)
2846   {
2847     SMESHDS_GroupBase* aGrp = (*GrIt);
2848     if ( !aGrp )
2849       continue;
2850     // check type and color of group
2851     if ( !isEqual( myColor, aGrp->GetColor() ))
2852       continue;
2853
2854     // IPAL52867 (prevent infinite recursion via GroupOnFilter)
2855     if ( SMESHDS_GroupOnFilter * gof = dynamic_cast< SMESHDS_GroupOnFilter* >( aGrp ))
2856       if ( gof->GetPredicate().get() == this )
2857         continue;
2858
2859     SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
2860     if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
2861       // add elements IDS into control
2862       int aSize = aGrp->Extent();
2863       for (int i = 0; i < aSize; i++)
2864         myIDs.insert( aGrp->GetID(i+1) );
2865     }
2866   }
2867 }
2868
2869 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
2870 {
2871   Kernel_Utils::Localizer loc;
2872   TCollection_AsciiString aStr = theStr;
2873   aStr.RemoveAll( ' ' );
2874   aStr.RemoveAll( '\t' );
2875   for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
2876     aStr.Remove( aPos, 2 );
2877   Standard_Real clr[3];
2878   clr[0] = clr[1] = clr[2] = 0.;
2879   for ( int i = 0; i < 3; i++ ) {
2880     TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
2881     if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
2882       clr[i] = tmpStr.RealValue();
2883   }
2884   myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
2885 }
2886
2887 //=======================================================================
2888 // name    : GetRangeStr
2889 // Purpose : Get range as a string.
2890 //           Example: "1,2,3,50-60,63,67,70-"
2891 //=======================================================================
2892
2893 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
2894 {
2895   theResStr.Clear();
2896   theResStr += TCollection_AsciiString( myColor.Red() );
2897   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
2898   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
2899 }
2900
2901 //================================================================================
2902 /*
2903   Class       : ElemGeomType
2904   Description : Predicate to check element geometry type
2905 */
2906 //================================================================================
2907
2908 ElemGeomType::ElemGeomType()
2909 {
2910   myMesh = 0;
2911   myType = SMDSAbs_All;
2912   myGeomType = SMDSGeom_TRIANGLE;
2913 }
2914
2915 void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
2916 {
2917   myMesh = theMesh;
2918 }
2919
2920 bool ElemGeomType::IsSatisfy( long theId )
2921 {
2922   if (!myMesh) return false;
2923   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2924   if ( !anElem )
2925     return false;
2926   const SMDSAbs_ElementType anElemType = anElem->GetType();
2927   if ( myType != SMDSAbs_All && anElemType != myType )
2928     return false;
2929   bool isOk = ( anElem->GetGeomType() == myGeomType );
2930   return isOk;
2931 }
2932
2933 void ElemGeomType::SetType( SMDSAbs_ElementType theType )
2934 {
2935   myType = theType;
2936 }
2937
2938 SMDSAbs_ElementType ElemGeomType::GetType() const
2939 {
2940   return myType;
2941 }
2942
2943 void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
2944 {
2945   myGeomType = theType;
2946 }
2947
2948 SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
2949 {
2950   return myGeomType;
2951 }
2952
2953 //================================================================================
2954 /*
2955   Class       : ElemEntityType
2956   Description : Predicate to check element entity type
2957 */
2958 //================================================================================
2959
2960 ElemEntityType::ElemEntityType():
2961   myMesh( 0 ),
2962   myType( SMDSAbs_All ),
2963   myEntityType( SMDSEntity_0D )
2964 {
2965 }
2966
2967 void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh )
2968 {
2969   myMesh = theMesh;
2970 }
2971
2972 bool ElemEntityType::IsSatisfy( long theId )
2973 {
2974   if ( !myMesh ) return false;
2975   if ( myType == SMDSAbs_Node )
2976     return myMesh->FindNode( theId );
2977   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2978   return ( anElem &&
2979            myEntityType == anElem->GetEntityType() );
2980 }
2981
2982 void ElemEntityType::SetType( SMDSAbs_ElementType theType )
2983 {
2984   myType = theType;
2985 }
2986
2987 SMDSAbs_ElementType ElemEntityType::GetType() const
2988 {
2989   return myType;
2990 }
2991
2992 void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType )
2993 {
2994   myEntityType = theEntityType;
2995 }
2996
2997 SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const
2998 {
2999   return myEntityType;
3000 }
3001
3002 //================================================================================
3003 /*!
3004  * \brief Class ConnectedElements
3005  */
3006 //================================================================================
3007
3008 ConnectedElements::ConnectedElements():
3009   myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {}
3010
3011 SMDSAbs_ElementType ConnectedElements::GetType() const
3012 { return myType; }
3013
3014 int ConnectedElements::GetNode() const
3015 { return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
3016
3017 std::vector<double> ConnectedElements::GetPoint() const
3018 { return myXYZ; }
3019
3020 void ConnectedElements::clearOkIDs()
3021 { myOkIDsReady = false; myOkIDs.clear(); }
3022
3023 void ConnectedElements::SetType( SMDSAbs_ElementType theType )
3024 {
3025   if ( myType != theType || myMeshModifTracer.IsMeshModified() )
3026     clearOkIDs();
3027   myType = theType;
3028 }
3029
3030 void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh )
3031 {
3032   myMeshModifTracer.SetMesh( theMesh );
3033   if ( myMeshModifTracer.IsMeshModified() )
3034   {
3035     clearOkIDs();
3036     if ( !myXYZ.empty() )
3037       SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh
3038   }
3039 }
3040
3041 void ConnectedElements::SetNode( int nodeID )
3042 {
3043   myNodeID = nodeID;
3044   myXYZ.clear();
3045
3046   bool isSameDomain = false;
3047   if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() )
3048     if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID ))
3049     {
3050       SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType );
3051       while ( !isSameDomain && eIt->more() )
3052         isSameDomain = IsSatisfy( eIt->next()->GetID() );
3053     }
3054   if ( !isSameDomain )
3055     clearOkIDs();
3056 }
3057
3058 void ConnectedElements::SetPoint( double x, double y, double z )
3059 {
3060   myXYZ.resize(3);
3061   myXYZ[0] = x;
3062   myXYZ[1] = y;
3063   myXYZ[2] = z;
3064   myNodeID = 0;
3065
3066   bool isSameDomain = false;
3067
3068   // find myNodeID by myXYZ if possible
3069   if ( myMeshModifTracer.GetMesh() )
3070   {
3071     SMESHUtils::Deleter<SMESH_ElementSearcher> searcher
3072       ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
3073
3074     std::vector< const SMDS_MeshElement* > foundElems;
3075     searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
3076
3077     if ( !foundElems.empty() )
3078     {
3079       myNodeID = foundElems[0]->GetNode(0)->GetID();
3080       if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() )
3081         isSameDomain = IsSatisfy( foundElems[0]->GetID() );
3082     }
3083   }
3084   if ( !isSameDomain )
3085     clearOkIDs();
3086 }
3087
3088 bool ConnectedElements::IsSatisfy( long theElementId )
3089 {
3090   // Here we do NOT check if the mesh has changed, we do it in Set...() only!!!
3091
3092   if ( !myOkIDsReady )
3093   {
3094     if ( !myMeshModifTracer.GetMesh() )
3095       return false;
3096     const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID );
3097     if ( !node0 )
3098       return false;
3099
3100     std::list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
3101     std::set< int > checkedNodeIDs;
3102     // algo:
3103     // foreach node in nodeQueue:
3104     //   foreach element sharing a node:
3105     //     add ID of an element of myType to myOkIDs;
3106     //     push all element nodes absent from checkedNodeIDs to nodeQueue;
3107     while ( !nodeQueue.empty() )
3108     {
3109       const SMDS_MeshNode* node = nodeQueue.front();
3110       nodeQueue.pop_front();
3111
3112       // loop on elements sharing the node
3113       SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
3114       while ( eIt->more() )
3115       {
3116         // keep elements of myType
3117         const SMDS_MeshElement* element = eIt->next();
3118         if ( myType == SMDSAbs_All || element->GetType() == myType )
3119           myOkIDs.insert( myOkIDs.end(), element->GetID() );
3120
3121         // enqueue nodes of the element
3122         SMDS_ElemIteratorPtr nIt = element->nodesIterator();
3123         while ( nIt->more() )
3124         {
3125           const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() );
3126           if ( checkedNodeIDs.insert( n->GetID() ).second )
3127             nodeQueue.push_back( n );
3128         }
3129       }
3130     }
3131     if ( myType == SMDSAbs_Node )
3132       std::swap( myOkIDs, checkedNodeIDs );
3133
3134     size_t totalNbElems = myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType );
3135     if ( myOkIDs.size() == totalNbElems )
3136       myOkIDs.clear();
3137
3138     myOkIDsReady = true;
3139   }
3140
3141   return myOkIDs.empty() ? true : myOkIDs.count( theElementId );
3142 }
3143
3144 //================================================================================
3145 /*!
3146  * \brief Class CoplanarFaces
3147  */
3148 //================================================================================
3149
3150 namespace
3151 {
3152   inline bool isLessAngle( const gp_Vec& v1, const gp_Vec& v2, const double cos )
3153   {
3154     double dot = v1 * v2; // cos * |v1| * |v2|
3155     double l1  = v1.SquareMagnitude();
3156     double l2  = v2.SquareMagnitude();
3157     return (( dot * cos >= 0 ) &&
3158             ( dot * dot ) / l1 / l2 >= ( cos * cos ));
3159   }
3160 }
3161 CoplanarFaces::CoplanarFaces()
3162   : myFaceID(0), myToler(0)
3163 {
3164 }
3165 void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
3166 {
3167   myMeshModifTracer.SetMesh( theMesh );
3168   if ( myMeshModifTracer.IsMeshModified() )
3169   {
3170     // Build a set of coplanar face ids
3171
3172     myCoplanarIDs.Clear();
3173
3174     if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler )
3175       return;
3176
3177     const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID );
3178     if ( !face || face->GetType() != SMDSAbs_Face )
3179       return;
3180
3181     bool normOK;
3182     gp_Vec myNorm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
3183     if (!normOK)
3184       return;
3185
3186     const double cosTol = Cos( myToler * M_PI / 180. );
3187     NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks;
3188
3189     std::list< std::pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
3190     faceQueue.push_back( std::make_pair( face, myNorm ));
3191     while ( !faceQueue.empty() )
3192     {
3193       face   = faceQueue.front().first;
3194       myNorm = faceQueue.front().second;
3195       faceQueue.pop_front();
3196
3197       for ( int i = 0, nbN = face->NbCornerNodes(); i < nbN; ++i )
3198       {
3199         const SMDS_MeshNode*  n1 = face->GetNode( i );
3200         const SMDS_MeshNode*  n2 = face->GetNode(( i+1 )%nbN);
3201         if ( !checkedLinks.Add( SMESH_TLink( n1, n2 )))
3202           continue;
3203         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face);
3204         while ( fIt->more() )
3205         {
3206           const SMDS_MeshElement* f = fIt->next();
3207           if ( f->GetNodeIndex( n2 ) > -1 )
3208           {
3209             gp_Vec norm = getNormale( static_cast<const SMDS_MeshFace*>(f), &normOK );
3210             if (!normOK || isLessAngle( myNorm, norm, cosTol))
3211             {
3212               myCoplanarIDs.Add( f->GetID() );
3213               faceQueue.push_back( std::make_pair( f, norm ));
3214             }
3215           }
3216         }
3217       }
3218     }
3219   }
3220 }
3221 bool CoplanarFaces::IsSatisfy( long theElementId )
3222 {
3223   return myCoplanarIDs.Contains( theElementId );
3224 }
3225
3226 /*
3227  *Class       : RangeOfIds
3228   *Description : Predicate for Range of Ids.
3229   *              Range may be specified with two ways.
3230   *              1. Using AddToRange method
3231   *              2. With SetRangeStr method. Parameter of this method is a string
3232   *                 like as "1,2,3,50-60,63,67,70-"
3233 */
3234
3235 //=======================================================================
3236 // name    : RangeOfIds
3237 // Purpose : Constructor
3238 //=======================================================================
3239 RangeOfIds::RangeOfIds()
3240 {
3241   myMesh = 0;
3242   myType = SMDSAbs_All;
3243 }
3244
3245 //=======================================================================
3246 // name    : SetMesh
3247 // Purpose : Set mesh
3248 //=======================================================================
3249 void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh )
3250 {
3251   myMesh = theMesh;
3252 }
3253
3254 //=======================================================================
3255 // name    : AddToRange
3256 // Purpose : Add ID to the range
3257 //=======================================================================
3258 bool RangeOfIds::AddToRange( long theEntityId )
3259 {
3260   myIds.Add( theEntityId );
3261   return true;
3262 }
3263
3264 //=======================================================================
3265 // name    : GetRangeStr
3266 // Purpose : Get range as a string.
3267 //           Example: "1,2,3,50-60,63,67,70-"
3268 //=======================================================================
3269 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
3270 {
3271   theResStr.Clear();
3272
3273   TColStd_SequenceOfInteger     anIntSeq;
3274   TColStd_SequenceOfAsciiString aStrSeq;
3275
3276   TColStd_MapIteratorOfMapOfInteger anIter( myIds );
3277   for ( ; anIter.More(); anIter.Next() )
3278   {
3279     int anId = anIter.Key();
3280     TCollection_AsciiString aStr( anId );
3281     anIntSeq.Append( anId );
3282     aStrSeq.Append( aStr );
3283   }
3284
3285   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
3286   {
3287     int aMinId = myMin( i );
3288     int aMaxId = myMax( i );
3289
3290     TCollection_AsciiString aStr;
3291     if ( aMinId != IntegerFirst() )
3292       aStr += aMinId;
3293
3294     aStr += "-";
3295
3296     if ( aMaxId != IntegerLast() )
3297       aStr += aMaxId;
3298
3299     // find position of the string in result sequence and insert string in it
3300     if ( anIntSeq.Length() == 0 )
3301     {
3302       anIntSeq.Append( aMinId );
3303       aStrSeq.Append( aStr );
3304     }
3305     else
3306     {
3307       if ( aMinId < anIntSeq.First() )
3308       {
3309         anIntSeq.Prepend( aMinId );
3310         aStrSeq.Prepend( aStr );
3311       }
3312       else if ( aMinId > anIntSeq.Last() )
3313       {
3314         anIntSeq.Append( aMinId );
3315         aStrSeq.Append( aStr );
3316       }
3317       else
3318         for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
3319           if ( aMinId < anIntSeq( j ) )
3320           {
3321             anIntSeq.InsertBefore( j, aMinId );
3322             aStrSeq.InsertBefore( j, aStr );
3323             break;
3324           }
3325     }
3326   }
3327
3328   if ( aStrSeq.Length() == 0 )
3329     return;
3330
3331   theResStr = aStrSeq( 1 );
3332   for ( int j = 2, k = aStrSeq.Length(); j <= k; j++  )
3333   {
3334     theResStr += ",";
3335     theResStr += aStrSeq( j );
3336   }
3337 }
3338
3339 //=======================================================================
3340 // name    : SetRangeStr
3341 // Purpose : Define range with string
3342 //           Example of entry string: "1,2,3,50-60,63,67,70-"
3343 //=======================================================================
3344 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
3345 {
3346   myMin.Clear();
3347   myMax.Clear();
3348   myIds.Clear();
3349
3350   TCollection_AsciiString aStr = theStr;
3351   for ( int i = 1; i <= aStr.Length(); ++i )
3352   {
3353     char c = aStr.Value( i );
3354     if ( !isdigit( c ) && c != ',' && c != '-' )
3355       aStr.SetValue( i, ',');
3356   }
3357   aStr.RemoveAll( ' ' );
3358
3359   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
3360   int i = 1;
3361   while ( tmpStr != "" )
3362   {
3363     tmpStr = aStr.Token( ",", i++ );
3364     int aPos = tmpStr.Search( '-' );
3365
3366     if ( aPos == -1 )
3367     {
3368       if ( tmpStr.IsIntegerValue() )
3369         myIds.Add( tmpStr.IntegerValue() );
3370       else
3371         return false;
3372     }
3373     else
3374     {
3375       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
3376       TCollection_AsciiString aMinStr = tmpStr;
3377
3378       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
3379       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
3380
3381       if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) ||
3382            (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
3383         return false;
3384
3385       myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
3386       myMax.Append( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
3387     }
3388   }
3389
3390   return true;
3391 }
3392
3393 //=======================================================================
3394 // name    : GetType
3395 // Purpose : Get type of supported entities
3396 //=======================================================================
3397 SMDSAbs_ElementType RangeOfIds::GetType() const
3398 {
3399   return myType;
3400 }
3401
3402 //=======================================================================
3403 // name    : SetType
3404 // Purpose : Set type of supported entities
3405 //=======================================================================
3406 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
3407 {
3408   myType = theType;
3409 }
3410
3411 //=======================================================================
3412 // name    : IsSatisfy
3413 // Purpose : Verify whether entity satisfies to this rpedicate
3414 //=======================================================================
3415 bool RangeOfIds::IsSatisfy( long theId )
3416 {
3417   if ( !myMesh )
3418     return false;
3419
3420   if ( myType == SMDSAbs_Node )
3421   {
3422     if ( myMesh->FindNode( theId ) == 0 )
3423       return false;
3424   }
3425   else
3426   {
3427     const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3428     if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All ))
3429       return false;
3430   }
3431
3432   if ( myIds.Contains( theId ) )
3433     return true;
3434
3435   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
3436     if ( theId >= myMin( i ) && theId <= myMax( i ) )
3437       return true;
3438
3439   return false;
3440 }
3441
3442 /*
3443   Class       : Comparator
3444   Description : Base class for comparators
3445 */
3446 Comparator::Comparator():
3447   myMargin(0)
3448 {}
3449
3450 Comparator::~Comparator()
3451 {}
3452
3453 void Comparator::SetMesh( const SMDS_Mesh* theMesh )
3454 {
3455   if ( myFunctor )
3456     myFunctor->SetMesh( theMesh );
3457 }
3458
3459 void Comparator::SetMargin( double theValue )
3460 {
3461   myMargin = theValue;
3462 }
3463
3464 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
3465 {
3466   myFunctor = theFunct;
3467 }
3468
3469 SMDSAbs_ElementType Comparator::GetType() const
3470 {
3471   return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
3472 }
3473
3474 double Comparator::GetMargin()
3475 {
3476   return myMargin;
3477 }
3478
3479
3480 /*
3481   Class       : LessThan
3482   Description : Comparator "<"
3483 */
3484 bool LessThan::IsSatisfy( long theId )
3485 {
3486   return myFunctor && myFunctor->GetValue( theId ) < myMargin;
3487 }
3488
3489
3490 /*
3491   Class       : MoreThan
3492   Description : Comparator ">"
3493 */
3494 bool MoreThan::IsSatisfy( long theId )
3495 {
3496   return myFunctor && myFunctor->GetValue( theId ) > myMargin;
3497 }
3498
3499
3500 /*
3501   Class       : EqualTo
3502   Description : Comparator "="
3503 */
3504 EqualTo::EqualTo():
3505   myToler(Precision::Confusion())
3506 {}
3507
3508 bool EqualTo::IsSatisfy( long theId )
3509 {
3510   return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
3511 }
3512
3513 void EqualTo::SetTolerance( double theToler )
3514 {
3515   myToler = theToler;
3516 }
3517
3518 double EqualTo::GetTolerance()
3519 {
3520   return myToler;
3521 }
3522
3523 /*
3524   Class       : LogicalNOT
3525   Description : Logical NOT predicate
3526 */
3527 LogicalNOT::LogicalNOT()
3528 {}
3529
3530 LogicalNOT::~LogicalNOT()
3531 {}
3532
3533 bool LogicalNOT::IsSatisfy( long theId )
3534 {
3535   return myPredicate && !myPredicate->IsSatisfy( theId );
3536 }
3537
3538 void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
3539 {
3540   if ( myPredicate )
3541     myPredicate->SetMesh( theMesh );
3542 }
3543
3544 void LogicalNOT::SetPredicate( PredicatePtr thePred )
3545 {
3546   myPredicate = thePred;
3547 }
3548
3549 SMDSAbs_ElementType LogicalNOT::GetType() const
3550 {
3551   return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
3552 }
3553
3554
3555 /*
3556   Class       : LogicalBinary
3557   Description : Base class for binary logical predicate
3558 */
3559 LogicalBinary::LogicalBinary()
3560 {}
3561
3562 LogicalBinary::~LogicalBinary()
3563 {}
3564
3565 void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
3566 {
3567   if ( myPredicate1 )
3568     myPredicate1->SetMesh( theMesh );
3569
3570   if ( myPredicate2 )
3571     myPredicate2->SetMesh( theMesh );
3572 }
3573
3574 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
3575 {
3576   myPredicate1 = thePredicate;
3577 }
3578
3579 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
3580 {
3581   myPredicate2 = thePredicate;
3582 }
3583
3584 SMDSAbs_ElementType LogicalBinary::GetType() const
3585 {
3586   if ( !myPredicate1 || !myPredicate2 )
3587     return SMDSAbs_All;
3588
3589   SMDSAbs_ElementType aType1 = myPredicate1->GetType();
3590   SMDSAbs_ElementType aType2 = myPredicate2->GetType();
3591
3592   return aType1 == aType2 ? aType1 : SMDSAbs_All;
3593 }
3594
3595
3596 /*
3597   Class       : LogicalAND
3598   Description : Logical AND
3599 */
3600 bool LogicalAND::IsSatisfy( long theId )
3601 {
3602   return
3603     myPredicate1 &&
3604     myPredicate2 &&
3605     myPredicate1->IsSatisfy( theId ) &&
3606     myPredicate2->IsSatisfy( theId );
3607 }
3608
3609
3610 /*
3611   Class       : LogicalOR
3612   Description : Logical OR
3613 */
3614 bool LogicalOR::IsSatisfy( long theId )
3615 {
3616   return
3617     myPredicate1 &&
3618     myPredicate2 &&
3619     (myPredicate1->IsSatisfy( theId ) ||
3620     myPredicate2->IsSatisfy( theId ));
3621 }
3622
3623
3624 /*
3625                               FILTER
3626 */
3627
3628 // #ifdef WITH_TBB
3629 // #include <tbb/parallel_for.h>
3630 // #include <tbb/enumerable_thread_specific.h>
3631
3632 // namespace Parallel
3633 // {
3634 //   typedef tbb::enumerable_thread_specific< TIdSequence > TIdSeq;
3635
3636 //   struct Predicate
3637 //   {
3638 //     const SMDS_Mesh* myMesh;
3639 //     PredicatePtr     myPredicate;
3640 //     TIdSeq &         myOKIds;
3641 //     Predicate( const SMDS_Mesh* m, PredicatePtr p, TIdSeq & ids ):
3642 //       myMesh(m), myPredicate(p->Duplicate()), myOKIds(ids) {}
3643 //     void operator() ( const tbb::blocked_range<size_t>& r ) const
3644 //     {
3645 //       for ( size_t i = r.begin(); i != r.end(); ++i )
3646 //         if ( myPredicate->IsSatisfy( i ))
3647 //           myOKIds.local().push_back();
3648 //     }
3649 //   }
3650 // }
3651 // #endif
3652
3653 Filter::Filter()
3654 {}
3655
3656 Filter::~Filter()
3657 {}
3658
3659 void Filter::SetPredicate( PredicatePtr thePredicate )
3660 {
3661   myPredicate = thePredicate;
3662 }
3663
3664 void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
3665                             PredicatePtr         thePredicate,
3666                             TIdSequence&         theSequence,
3667                             SMDS_ElemIteratorPtr theElements )
3668 {
3669   theSequence.clear();
3670
3671   if ( !theMesh || !thePredicate )
3672     return;
3673
3674   thePredicate->SetMesh( theMesh );
3675
3676   if ( !theElements )
3677     theElements = theMesh->elementsIterator( thePredicate->GetType() );
3678
3679   if ( theElements ) {
3680     while ( theElements->more() ) {
3681       const SMDS_MeshElement* anElem = theElements->next();
3682       if ( thePredicate->GetType() == SMDSAbs_All ||
3683            thePredicate->GetType() == anElem->GetType() )
3684       {
3685         long anId = anElem->GetID();
3686         if ( thePredicate->IsSatisfy( anId ) )
3687           theSequence.push_back( anId );
3688       }
3689     }
3690   }
3691 }
3692
3693 void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
3694                             Filter::TIdSequence& theSequence,
3695                             SMDS_ElemIteratorPtr theElements )
3696 {
3697   GetElementsId(theMesh,myPredicate,theSequence,theElements);
3698 }
3699
3700 /*
3701                               ManifoldPart
3702 */
3703
3704 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
3705
3706 /*
3707    Internal class Link
3708 */
3709
3710 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
3711                           SMDS_MeshNode* theNode2 )
3712 {
3713   myNode1 = theNode1;
3714   myNode2 = theNode2;
3715 }
3716
3717 ManifoldPart::Link::~Link()
3718 {
3719   myNode1 = 0;
3720   myNode2 = 0;
3721 }
3722
3723 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
3724 {
3725   if ( myNode1 == theLink.myNode1 &&
3726        myNode2 == theLink.myNode2 )
3727     return true;
3728   else if ( myNode1 == theLink.myNode2 &&
3729             myNode2 == theLink.myNode1 )
3730     return true;
3731   else
3732     return false;
3733 }
3734
3735 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
3736 {
3737   if(myNode1 < x.myNode1) return true;
3738   if(myNode1 == x.myNode1)
3739     if(myNode2 < x.myNode2) return true;
3740   return false;
3741 }
3742
3743 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
3744                             const ManifoldPart::Link& theLink2 )
3745 {
3746   return theLink1.IsEqual( theLink2 );
3747 }
3748
3749 ManifoldPart::ManifoldPart()
3750 {
3751   myMesh = 0;
3752   myAngToler = Precision::Angular();
3753   myIsOnlyManifold = true;
3754 }
3755
3756 ManifoldPart::~ManifoldPart()
3757 {
3758   myMesh = 0;
3759 }
3760
3761 void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
3762 {
3763   myMesh = theMesh;
3764   process();
3765 }
3766
3767 SMDSAbs_ElementType ManifoldPart::GetType() const
3768 { return SMDSAbs_Face; }
3769
3770 bool ManifoldPart::IsSatisfy( long theElementId )
3771 {
3772   return myMapIds.Contains( theElementId );
3773 }
3774
3775 void ManifoldPart::SetAngleTolerance( const double theAngToler )
3776 { myAngToler = theAngToler; }
3777
3778 double ManifoldPart::GetAngleTolerance() const
3779 { return myAngToler; }
3780
3781 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
3782 { myIsOnlyManifold = theIsOnly; }
3783
3784 void ManifoldPart::SetStartElem( const long  theStartId )
3785 { myStartElemId = theStartId; }
3786
3787 bool ManifoldPart::process()
3788 {
3789   myMapIds.Clear();
3790   myMapBadGeomIds.Clear();
3791
3792   myAllFacePtr.clear();
3793   myAllFacePtrIntDMap.clear();
3794   if ( !myMesh )
3795     return false;
3796
3797   // collect all faces into own map
3798   SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
3799   for (; anFaceItr->more(); )
3800   {
3801     SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
3802     myAllFacePtr.push_back( aFacePtr );
3803     myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
3804   }
3805
3806   SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
3807   if ( !aStartFace )
3808     return false;
3809
3810   // the map of non manifold links and bad geometry
3811   TMapOfLink aMapOfNonManifold;
3812   TColStd_MapOfInteger aMapOfTreated;
3813
3814   // begin cycle on faces from start index and run on vector till the end
3815   //  and from begin to start index to cover whole vector
3816   const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
3817   bool isStartTreat = false;
3818   for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
3819   {
3820     if ( fi == aStartIndx )
3821       isStartTreat = true;
3822     // as result next time when fi will be equal to aStartIndx
3823
3824     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
3825     if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
3826       continue;
3827
3828     aMapOfTreated.Add( aFacePtr->GetID() );
3829     TColStd_MapOfInteger aResFaces;
3830     if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
3831                          aMapOfNonManifold, aResFaces ) )
3832       continue;
3833     TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
3834     for ( ; anItr.More(); anItr.Next() )
3835     {
3836       int aFaceId = anItr.Key();
3837       aMapOfTreated.Add( aFaceId );
3838       myMapIds.Add( aFaceId );
3839     }
3840
3841     if ( fi == int( myAllFacePtr.size() - 1 ))
3842       fi = 0;
3843   } // end run on vector of faces
3844   return !myMapIds.IsEmpty();
3845 }
3846
3847 static void getLinks( const SMDS_MeshFace* theFace,
3848                       ManifoldPart::TVectorOfLink& theLinks )
3849 {
3850   int aNbNode = theFace->NbNodes();
3851   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
3852   int i = 1;
3853   SMDS_MeshNode* aNode = 0;
3854   for ( ; aNodeItr->more() && i <= aNbNode; )
3855   {
3856
3857     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
3858     if ( i == 1 )
3859       aNode = aN1;
3860     i++;
3861     SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
3862     i++;
3863     ManifoldPart::Link aLink( aN1, aN2 );
3864     theLinks.push_back( aLink );
3865   }
3866 }
3867
3868 bool ManifoldPart::findConnected
3869                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
3870                   SMDS_MeshFace*                           theStartFace,
3871                   ManifoldPart::TMapOfLink&                theNonManifold,
3872                   TColStd_MapOfInteger&                    theResFaces )
3873 {
3874   theResFaces.Clear();
3875   if ( !theAllFacePtrInt.size() )
3876     return false;
3877
3878   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
3879   {
3880     myMapBadGeomIds.Add( theStartFace->GetID() );
3881     return false;
3882   }
3883
3884   ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
3885   ManifoldPart::TVectorOfLink aSeqOfBoundary;
3886   theResFaces.Add( theStartFace->GetID() );
3887   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
3888
3889   expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3890                  aDMapLinkFace, theNonManifold, theStartFace );
3891
3892   bool isDone = false;
3893   while ( !isDone && aMapOfBoundary.size() != 0 )
3894   {
3895     bool isToReset = false;
3896     ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
3897     for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
3898     {
3899       ManifoldPart::Link aLink = *pLink;
3900       if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
3901         continue;
3902       // each link could be treated only once
3903       aMapToSkip.insert( aLink );
3904
3905       ManifoldPart::TVectorOfFacePtr aFaces;
3906       // find next
3907       if ( myIsOnlyManifold &&
3908            (theNonManifold.find( aLink ) != theNonManifold.end()) )
3909         continue;
3910       else
3911       {
3912         getFacesByLink( aLink, aFaces );
3913         // filter the element to keep only indicated elements
3914         ManifoldPart::TVectorOfFacePtr aFiltered;
3915         ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3916         for ( ; pFace != aFaces.end(); ++pFace )
3917         {
3918           SMDS_MeshFace* aFace = *pFace;
3919           if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
3920             aFiltered.push_back( aFace );
3921         }
3922         aFaces = aFiltered;
3923         if ( aFaces.size() < 2 )  // no neihgbour faces
3924           continue;
3925         else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
3926         {
3927           theNonManifold.insert( aLink );
3928           continue;
3929         }
3930       }
3931
3932       // compare normal with normals of neighbor element
3933       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
3934       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3935       for ( ; pFace != aFaces.end(); ++pFace )
3936       {
3937         SMDS_MeshFace* aNextFace = *pFace;
3938         if ( aPrevFace == aNextFace )
3939           continue;
3940         int anNextFaceID = aNextFace->GetID();
3941         if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
3942          // should not be with non manifold restriction. probably bad topology
3943           continue;
3944         // check if face was treated and skipped
3945         if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
3946              !isInPlane( aPrevFace, aNextFace ) )
3947           continue;
3948         // add new element to connected and extend the boundaries.
3949         theResFaces.Add( anNextFaceID );
3950         expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3951                         aDMapLinkFace, theNonManifold, aNextFace );
3952         isToReset = true;
3953       }
3954     }
3955     isDone = !isToReset;
3956   }
3957
3958   return !theResFaces.IsEmpty();
3959 }
3960
3961 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
3962                               const SMDS_MeshFace* theFace2 )
3963 {
3964   gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
3965   gp_XYZ aNorm2XYZ = getNormale( theFace2 );
3966   if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
3967   {
3968     myMapBadGeomIds.Add( theFace2->GetID() );
3969     return false;
3970   }
3971   if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
3972     return true;
3973
3974   return false;
3975 }
3976
3977 void ManifoldPart::expandBoundary
3978                    ( ManifoldPart::TMapOfLink&            theMapOfBoundary,
3979                      ManifoldPart::TVectorOfLink&         theSeqOfBoundary,
3980                      ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
3981                      ManifoldPart::TMapOfLink&            theNonManifold,
3982                      SMDS_MeshFace*                       theNextFace ) const
3983 {
3984   ManifoldPart::TVectorOfLink aLinks;
3985   getLinks( theNextFace, aLinks );
3986   int aNbLink = (int)aLinks.size();
3987   for ( int i = 0; i < aNbLink; i++ )
3988   {
3989     ManifoldPart::Link aLink = aLinks[ i ];
3990     if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
3991       continue;
3992     if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
3993     {
3994       if ( myIsOnlyManifold )
3995       {
3996         // remove from boundary
3997         theMapOfBoundary.erase( aLink );
3998         ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
3999         for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
4000         {
4001           ManifoldPart::Link aBoundLink = *pLink;
4002           if ( aBoundLink.IsEqual( aLink ) )
4003           {
4004             theSeqOfBoundary.erase( pLink );
4005             break;
4006           }
4007         }
4008       }
4009     }
4010     else
4011     {
4012       theMapOfBoundary.insert( aLink );
4013       theSeqOfBoundary.push_back( aLink );
4014       theDMapLinkFacePtr[ aLink ] = theNextFace;
4015     }
4016   }
4017 }
4018
4019 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
4020                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
4021 {
4022
4023   // take all faces that shared first node
4024   SMDS_ElemIteratorPtr anItr = theLink.myNode1->GetInverseElementIterator( SMDSAbs_Face );
4025   SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > faces( anItr ), facesEnd;
4026   std::set<const SMDS_MeshElement *> aSetOfFaces( faces, facesEnd );
4027
4028   // take all faces that shared second node
4029   anItr = theLink.myNode2->GetInverseElementIterator( SMDSAbs_Face );
4030   // find the common part of two sets
4031   for ( ; anItr->more(); )
4032   {
4033     const SMDS_MeshElement* aFace = anItr->next();
4034     if ( aSetOfFaces.count( aFace ))
4035       theFaces.push_back( (SMDS_MeshFace*) aFace );
4036   }
4037 }
4038
4039 /*
4040   Class       : BelongToMeshGroup
4041   Description : Verify whether a mesh element is included into a mesh group
4042 */
4043 BelongToMeshGroup::BelongToMeshGroup(): myGroup( 0 )
4044 {
4045 }
4046
4047 void BelongToMeshGroup::SetGroup( SMESHDS_GroupBase* g )
4048 {
4049   myGroup = g;
4050 }
4051
4052 void BelongToMeshGroup::SetStoreName( const std::string& sn )
4053 {
4054   myStoreName = sn;
4055 }
4056
4057 void BelongToMeshGroup::SetMesh( const SMDS_Mesh* theMesh )
4058 {
4059   if ( myGroup && myGroup->GetMesh() != theMesh )
4060   {
4061     myGroup = 0;
4062   }
4063   if ( !myGroup && !myStoreName.empty() )
4064   {
4065     if ( const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh))
4066     {
4067       const std::set<SMESHDS_GroupBase*>& grps = aMesh->GetGroups();
4068       std::set<SMESHDS_GroupBase*>::const_iterator g = grps.begin();
4069       for ( ; g != grps.end() && !myGroup; ++g )
4070         if ( *g && myStoreName == (*g)->GetStoreName() )
4071           myGroup = *g;
4072     }
4073   }
4074   if ( myGroup )
4075   {
4076     myGroup->IsEmpty(); // make GroupOnFilter update its predicate
4077   }
4078 }
4079
4080 bool BelongToMeshGroup::IsSatisfy( long theElementId )
4081 {
4082   return myGroup ? myGroup->Contains( theElementId ) : false;
4083 }
4084
4085 SMDSAbs_ElementType BelongToMeshGroup::GetType() const
4086 {
4087   return myGroup ? myGroup->GetType() : SMDSAbs_All;
4088 }
4089
4090 //================================================================================
4091 //  ElementsOnSurface
4092 //================================================================================
4093
4094 ElementsOnSurface::ElementsOnSurface()
4095 {
4096   myIds.Clear();
4097   myType = SMDSAbs_All;
4098   mySurf.Nullify();
4099   myToler = Precision::Confusion();
4100   myUseBoundaries = false;
4101 }
4102
4103 ElementsOnSurface::~ElementsOnSurface()
4104 {
4105 }
4106
4107 void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
4108 {
4109   myMeshModifTracer.SetMesh( theMesh );
4110   if ( myMeshModifTracer.IsMeshModified())
4111     process();
4112 }
4113
4114 bool ElementsOnSurface::IsSatisfy( long theElementId )
4115 {
4116   return myIds.Contains( theElementId );
4117 }
4118
4119 SMDSAbs_ElementType ElementsOnSurface::GetType() const
4120 { return myType; }
4121
4122 void ElementsOnSurface::SetTolerance( const double theToler )
4123 {
4124   if ( myToler != theToler )
4125   {
4126     myToler = theToler;
4127     process();
4128   }
4129 }
4130
4131 double ElementsOnSurface::GetTolerance() const
4132 { return myToler; }
4133
4134 void ElementsOnSurface::SetUseBoundaries( bool theUse )
4135 {
4136   if ( myUseBoundaries != theUse ) {
4137     myUseBoundaries = theUse;
4138     SetSurface( mySurf, myType );
4139   }
4140 }
4141
4142 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
4143                                     const SMDSAbs_ElementType theType )
4144 {
4145   myIds.Clear();
4146   myType = theType;
4147   mySurf.Nullify();
4148   if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
4149     return;
4150   mySurf = TopoDS::Face( theShape );
4151   BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
4152   Standard_Real
4153     u1 = SA.FirstUParameter(),
4154     u2 = SA.LastUParameter(),
4155     v1 = SA.FirstVParameter(),
4156     v2 = SA.LastVParameter();
4157   Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
4158   myProjector.Init( surf, u1,u2, v1,v2 );
4159   process();
4160 }
4161
4162 void ElementsOnSurface::process()
4163 {
4164   myIds.Clear();
4165   if ( mySurf.IsNull() )
4166     return;
4167
4168   if ( !myMeshModifTracer.GetMesh() )
4169     return;
4170
4171   myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
4172
4173   SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType );
4174   for(; anIter->more(); )
4175     process( anIter->next() );
4176 }
4177
4178 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
4179 {
4180   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
4181   bool isSatisfy = true;
4182   for ( ; aNodeItr->more(); )
4183   {
4184     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
4185     if ( !isOnSurface( aNode ) )
4186     {
4187       isSatisfy = false;
4188       break;
4189     }
4190   }
4191   if ( isSatisfy )
4192     myIds.Add( theElemPtr->GetID() );
4193 }
4194
4195 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
4196 {
4197   if ( mySurf.IsNull() )
4198     return false;
4199
4200   gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
4201   //  double aToler2 = myToler * myToler;
4202 //   if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
4203 //   {
4204 //     gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
4205 //     if ( aPln.SquareDistance( aPnt ) > aToler2 )
4206 //       return false;
4207 //   }
4208 //   else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
4209 //   {
4210 //     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
4211 //     double aRad = aCyl.Radius();
4212 //     gp_Ax3 anAxis = aCyl.Position();
4213 //     gp_XYZ aLoc = aCyl.Location().XYZ();
4214 //     double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4215 //     double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4216 //     if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
4217 //       return false;
4218 //   }
4219 //   else
4220 //     return false;
4221   myProjector.Perform( aPnt );
4222   bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
4223
4224   return isOn;
4225 }
4226
4227
4228 //================================================================================
4229 //  ElementsOnShape
4230 //================================================================================
4231
4232 namespace {
4233   const int theIsCheckedFlag = 0x0000100;
4234 }
4235
4236 struct ElementsOnShape::Classifier
4237 {
4238   Classifier() { mySolidClfr = 0; myFlags = 0; }
4239   ~Classifier();
4240   void Init(const TopoDS_Shape& s, double tol, const Bnd_B3d* box = 0 );
4241   bool IsOut(const gp_Pnt& p)        { return SetChecked( true ), (this->*myIsOutFun)( p ); }
4242   TopAbs_ShapeEnum ShapeType() const { return myShape.ShapeType(); }
4243   const TopoDS_Shape& Shape() const  { return myShape; }
4244   const Bnd_B3d* GetBndBox() const   { return & myBox; }
4245   bool IsChecked()                   { return myFlags & theIsCheckedFlag; }
4246   bool IsSetFlag( int flag ) const   { return myFlags & flag; }
4247   void SetChecked( bool is ) { is ? SetFlag( theIsCheckedFlag ) : UnsetFlag( theIsCheckedFlag ); }
4248   void SetFlag  ( int flag ) { myFlags |= flag; }
4249   void UnsetFlag( int flag ) { myFlags &= ~flag; }
4250
4251 private:
4252   bool isOutOfSolid (const gp_Pnt& p);
4253   bool isOutOfBox   (const gp_Pnt& p);
4254   bool isOutOfFace  (const gp_Pnt& p);
4255   bool isOutOfEdge  (const gp_Pnt& p);
4256   bool isOutOfVertex(const gp_Pnt& p);
4257   bool isBox        (const TopoDS_Shape& s);
4258
4259   bool (Classifier::*          myIsOutFun)(const gp_Pnt& p);
4260   BRepClass3d_SolidClassifier* mySolidClfr; // ptr because of a run-time forbidden copy-constructor
4261   Bnd_B3d                      myBox;
4262   GeomAPI_ProjectPointOnSurf   myProjFace;
4263   GeomAPI_ProjectPointOnCurve  myProjEdge;
4264   gp_Pnt                       myVertexXYZ;
4265   TopoDS_Shape                 myShape;
4266   double                       myTol;
4267   int                          myFlags;
4268 };
4269
4270 struct ElementsOnShape::OctreeClassifier : public SMESH_Octree
4271 {
4272   OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers );
4273   OctreeClassifier( const OctreeClassifier*                           otherTree,
4274                     const std::vector< ElementsOnShape::Classifier >& clsOther,
4275                     std::vector< ElementsOnShape::Classifier >&       cls );
4276   void GetClassifiersAtPoint( const gp_XYZ& p,
4277                               std::vector< ElementsOnShape::Classifier* >& classifiers );
4278 protected:
4279   OctreeClassifier() {}
4280   SMESH_Octree* newChild() const { return new OctreeClassifier; }
4281   void          buildChildrenData();
4282   Bnd_B3d*      buildRootBox();
4283
4284   std::vector< ElementsOnShape::Classifier* > myClassifiers;
4285 };
4286
4287
4288 ElementsOnShape::ElementsOnShape():
4289   myOctree(0),
4290   myType(SMDSAbs_All),
4291   myToler(Precision::Confusion()),
4292   myAllNodesFlag(false)
4293 {
4294 }
4295
4296 ElementsOnShape::~ElementsOnShape()
4297 {
4298   clearClassifiers();
4299 }
4300
4301 Predicate* ElementsOnShape::clone() const
4302 {
4303   ElementsOnShape* cln = new ElementsOnShape();
4304   cln->SetAllNodes ( myAllNodesFlag );
4305   cln->SetTolerance( myToler );
4306   cln->SetMesh     ( myMeshModifTracer.GetMesh() );
4307   cln->myShape = myShape; // avoid creation of myClassifiers
4308   cln->SetShape    ( myShape, myType );
4309   cln->myClassifiers.resize( myClassifiers.size() );
4310   for ( size_t i = 0; i < myClassifiers.size(); ++i )
4311     cln->myClassifiers[ i ].Init( BRepBuilderAPI_Copy( myClassifiers[ i ].Shape()),
4312                                   myToler, myClassifiers[ i ].GetBndBox() );
4313   if ( myOctree ) // copy myOctree
4314   {
4315     cln->myOctree = new OctreeClassifier( myOctree, myClassifiers, cln->myClassifiers );
4316   }
4317   return cln;
4318 }
4319
4320 SMDSAbs_ElementType ElementsOnShape::GetType() const
4321 {
4322   return myType;
4323 }
4324
4325 void ElementsOnShape::SetTolerance (const double theToler)
4326 {
4327   if (myToler != theToler) {
4328     myToler = theToler;
4329     SetShape(myShape, myType);
4330   }
4331 }
4332
4333 double ElementsOnShape::GetTolerance() const
4334 {
4335   return myToler;
4336 }
4337
4338 void ElementsOnShape::SetAllNodes (bool theAllNodes)
4339 {
4340   myAllNodesFlag = theAllNodes;
4341 }
4342
4343 void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
4344 {
4345   myMeshModifTracer.SetMesh( theMesh );
4346   if ( myMeshModifTracer.IsMeshModified())
4347   {
4348     size_t nbNodes = theMesh ? theMesh->NbNodes() : 0;
4349     if ( myNodeIsChecked.size() == nbNodes )
4350     {
4351       std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4352     }
4353     else
4354     {
4355       SMESHUtils::FreeVector( myNodeIsChecked );
4356       SMESHUtils::FreeVector( myNodeIsOut );
4357       myNodeIsChecked.resize( nbNodes, false );
4358       myNodeIsOut.resize( nbNodes );
4359     }
4360   }
4361 }
4362
4363 bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut )
4364 {
4365   if ( n->GetID() >= (int) myNodeIsChecked.size() ||
4366        !myNodeIsChecked[ n->GetID() ])
4367     return false;
4368
4369   isOut = myNodeIsOut[ n->GetID() ];
4370   return true;
4371 }
4372
4373 void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool  isOut )
4374 {
4375   if ( n->GetID() < (int) myNodeIsChecked.size() )
4376   {
4377     myNodeIsChecked[ n->GetID() ] = true;
4378     myNodeIsOut    [ n->GetID() ] = isOut;
4379   }
4380 }
4381
4382 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
4383                                 const SMDSAbs_ElementType theType)
4384 {
4385   bool shapeChanges = ( myShape != theShape );
4386   myType  = theType;
4387   myShape = theShape;
4388   if ( myShape.IsNull() ) return;
4389
4390   if ( shapeChanges )
4391   {
4392     // find most complex shapes
4393     TopTools_IndexedMapOfShape shapesMap;
4394     TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
4395     TopExp_Explorer sub;
4396     for ( int i = 0; i < 4; ++i )
4397     {
4398       if ( shapesMap.IsEmpty() )
4399         for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() )
4400           shapesMap.Add( sub.Current() );
4401       if ( i > 0 )
4402         for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() )
4403           shapesMap.Add( sub.Current() );
4404     }
4405
4406     clearClassifiers();
4407     myClassifiers.resize( shapesMap.Extent() );
4408     for ( int i = 0; i < shapesMap.Extent(); ++i )
4409       myClassifiers[ i ].Init( shapesMap( i+1 ), myToler );
4410   }
4411
4412   if ( theType == SMDSAbs_Node )
4413   {
4414     SMESHUtils::FreeVector( myNodeIsChecked );
4415     SMESHUtils::FreeVector( myNodeIsOut );
4416   }
4417   else
4418   {
4419     std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4420   }
4421 }
4422
4423 void ElementsOnShape::clearClassifiers()
4424 {
4425   // for ( size_t i = 0; i < myClassifiers.size(); ++i )
4426   //   delete myClassifiers[ i ];
4427   myClassifiers.clear();
4428
4429   delete myOctree;
4430   myOctree = 0;
4431 }
4432
4433 bool ElementsOnShape::IsSatisfy( long elemId )
4434 {
4435   if ( myClassifiers.empty() )
4436     return false;
4437
4438   const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh();
4439   if ( myType == SMDSAbs_Node )
4440     return IsSatisfy( mesh->FindNode( elemId ));
4441   return IsSatisfy( mesh->FindElement( elemId ));
4442 }
4443
4444 bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem)
4445 {
4446   if ( !elem )
4447     return false;
4448
4449   bool isSatisfy = myAllNodesFlag, isNodeOut;
4450
4451   gp_XYZ centerXYZ (0, 0, 0);
4452
4453   if ( !myOctree && myClassifiers.size() > 5 )
4454   {
4455     myWorkClassifiers.resize( myClassifiers.size() );
4456     for ( size_t i = 0; i < myClassifiers.size(); ++i )
4457       myWorkClassifiers[ i ] = & myClassifiers[ i ];
4458     myOctree = new OctreeClassifier( myWorkClassifiers );
4459   }
4460
4461   SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
4462   while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
4463   {
4464     SMESH_TNodeXYZ aPnt( aNodeItr->next() );
4465     centerXYZ += aPnt;
4466
4467     isNodeOut = true;
4468     if ( !getNodeIsOut( aPnt._node, isNodeOut ))
4469     {
4470       if ( myOctree )
4471       {
4472         myWorkClassifiers.clear();
4473         myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers );
4474
4475         for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
4476           myWorkClassifiers[i]->SetChecked( false );
4477
4478         for ( size_t i = 0; i < myWorkClassifiers.size() && isNodeOut; ++i )
4479           if ( !myWorkClassifiers[i]->IsChecked() )
4480             isNodeOut = myWorkClassifiers[i]->IsOut( aPnt );
4481       }
4482       else
4483       {
4484         for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i )
4485           isNodeOut = myClassifiers[i].IsOut( aPnt );
4486       }
4487       setNodeIsOut( aPnt._node, isNodeOut );
4488     }
4489     isSatisfy = !isNodeOut;
4490   }
4491
4492   // Check the center point for volumes MantisBug 0020168
4493   if ( isSatisfy &&
4494        myAllNodesFlag &&
4495        myClassifiers[0].ShapeType() == TopAbs_SOLID )
4496   {
4497     centerXYZ /= elem->NbNodes();
4498     isSatisfy = false;
4499     if ( myOctree )
4500       for ( size_t i = 0; i < myWorkClassifiers.size() && !isSatisfy; ++i )
4501         isSatisfy = ! myWorkClassifiers[i]->IsOut( centerXYZ );
4502     else
4503       for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i )
4504         isSatisfy = ! myClassifiers[i].IsOut( centerXYZ );
4505   }
4506
4507   return isSatisfy;
4508 }
4509
4510 //================================================================================
4511 /*!
4512  * \brief Check and optionally return a satisfying shape
4513  */
4514 //================================================================================
4515
4516 bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node,
4517                                  TopoDS_Shape*        okShape)
4518 {
4519   if ( !node )
4520     return false;
4521
4522   if ( !myOctree && myClassifiers.size() > 5 )
4523   {
4524     myWorkClassifiers.resize( myClassifiers.size() );
4525     for ( size_t i = 0; i < myClassifiers.size(); ++i )
4526       myWorkClassifiers[ i ] = & myClassifiers[ i ];
4527     myOctree = new OctreeClassifier( myWorkClassifiers );
4528   }
4529
4530   bool isNodeOut = true;
4531
4532   if ( okShape || !getNodeIsOut( node, isNodeOut ))
4533   {
4534     SMESH_NodeXYZ aPnt = node;
4535     if ( myOctree )
4536     {
4537       myWorkClassifiers.clear();
4538       myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers );
4539
4540       for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
4541         myWorkClassifiers[i]->SetChecked( false );
4542
4543       for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
4544         if ( !myWorkClassifiers[i]->IsChecked() &&
4545              !myWorkClassifiers[i]->IsOut( aPnt ))
4546         {
4547           isNodeOut = false;
4548           if ( okShape )
4549             *okShape = myWorkClassifiers[i]->Shape();
4550           break;
4551         }
4552     }
4553     else
4554     {
4555       for ( size_t i = 0; i < myClassifiers.size(); ++i )
4556         if ( !myClassifiers[i].IsOut( aPnt ))
4557         {
4558           isNodeOut = false;
4559           if ( okShape )
4560             *okShape = myWorkClassifiers[i]->Shape();
4561           break;
4562         }
4563     }
4564     setNodeIsOut( node, isNodeOut );
4565   }
4566
4567   return !isNodeOut;
4568 }
4569
4570 void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape,
4571                                         double              theTol,
4572                                         const Bnd_B3d*      theBox )
4573 {
4574   myShape = theShape;
4575   myTol   = theTol;
4576   myFlags = 0;
4577
4578   bool isShapeBox = false;
4579   switch ( myShape.ShapeType() )
4580   {
4581   case TopAbs_SOLID:
4582   {
4583     if (( isShapeBox = isBox( theShape )))
4584     {
4585       myIsOutFun = & ElementsOnShape::Classifier::isOutOfBox;
4586     }
4587     else
4588     {
4589       mySolidClfr = new BRepClass3d_SolidClassifier(theShape);
4590       myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid;
4591     }
4592     break;
4593   }
4594   case TopAbs_FACE:
4595   {
4596     Standard_Real u1,u2,v1,v2;
4597     Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
4598     surf->Bounds( u1,u2,v1,v2 );
4599     myProjFace.Init(surf, u1,u2, v1,v2, myTol );
4600     myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace;
4601     break;
4602   }
4603   case TopAbs_EDGE:
4604   {
4605     Standard_Real u1, u2;
4606     Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2);
4607     myProjEdge.Init(curve, u1, u2);
4608     myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
4609     break;
4610   }
4611   case TopAbs_VERTEX:
4612   {
4613     myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
4614     myIsOutFun = & ElementsOnShape::Classifier::isOutOfVertex;
4615     break;
4616   }
4617   default:
4618     throw SALOME_Exception("Programmer error in usage of ElementsOnShape::Classifier");
4619   }
4620
4621   if ( !isShapeBox )
4622   {
4623     if ( theBox )
4624     {
4625       myBox = *theBox;
4626     }
4627     else
4628     {
4629       Bnd_Box box;
4630       BRepBndLib::Add( myShape, box );
4631       myBox.Clear();
4632       myBox.Add( box.CornerMin() );
4633       myBox.Add( box.CornerMax() );
4634       gp_XYZ halfSize = 0.5 * ( box.CornerMax().XYZ() - box.CornerMin().XYZ() );
4635       for ( int iDim = 1; iDim <= 3; ++iDim )
4636       {
4637         double x = halfSize.Coord( iDim );
4638         halfSize.SetCoord( iDim, x + Max( myTol, 1e-2 * x ));
4639       }
4640       myBox.SetHSize( halfSize );
4641     }
4642   }
4643 }
4644
4645 ElementsOnShape::Classifier::~Classifier()
4646 {
4647   delete mySolidClfr; mySolidClfr = 0;
4648 }
4649
4650 bool ElementsOnShape::Classifier::isOutOfSolid (const gp_Pnt& p)
4651 {
4652   if ( isOutOfBox( p )) return true;
4653   mySolidClfr->Perform( p, myTol );
4654   return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON );
4655 }
4656
4657 bool ElementsOnShape::Classifier::isOutOfBox (const gp_Pnt& p)
4658 {
4659   return myBox.IsOut( p.XYZ() );
4660 }
4661
4662 bool ElementsOnShape::Classifier::isOutOfFace  (const gp_Pnt& p)
4663 {
4664   if ( isOutOfBox( p )) return true;
4665   myProjFace.Perform( p );
4666   if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
4667   {
4668     // check relatively to the face
4669     Standard_Real u, v;
4670     myProjFace.LowerDistanceParameters(u, v);
4671     gp_Pnt2d aProjPnt (u, v);
4672     BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
4673     if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
4674       return false;
4675   }
4676   return true;
4677 }
4678
4679 bool ElementsOnShape::Classifier::isOutOfEdge  (const gp_Pnt& p)
4680 {
4681   if ( isOutOfBox( p )) return true;
4682   myProjEdge.Perform( p );
4683   return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
4684 }
4685
4686 bool ElementsOnShape::Classifier::isOutOfVertex(const gp_Pnt& p)
4687 {
4688   return ( myVertexXYZ.Distance( p ) > myTol );
4689 }
4690
4691 bool ElementsOnShape::Classifier::isBox (const TopoDS_Shape& theShape)
4692 {
4693   TopTools_IndexedMapOfShape vMap;
4694   TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
4695   if ( vMap.Extent() != 8 )
4696     return false;
4697
4698   myBox.Clear();
4699   for ( int i = 1; i <= 8; ++i )
4700     myBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))).XYZ() );
4701
4702   gp_XYZ pMin = myBox.CornerMin(), pMax = myBox.CornerMax();
4703   for ( int i = 1; i <= 8; ++i )
4704   {
4705     gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i )));
4706     for ( int iC = 1; iC <= 3; ++ iC )
4707     {
4708       double d1 = Abs( pMin.Coord( iC ) - p.Coord( iC ));
4709       double d2 = Abs( pMax.Coord( iC ) - p.Coord( iC ));
4710       if ( Min( d1, d2 ) > myTol )
4711         return false;
4712     }
4713   }
4714   myBox.Enlarge( myTol );
4715   return true;
4716 }
4717
4718 ElementsOnShape::
4719 OctreeClassifier::OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers )
4720   :SMESH_Octree( new SMESH_TreeLimit )
4721 {
4722   myClassifiers = classifiers;
4723   compute();
4724 }
4725
4726 ElementsOnShape::
4727 OctreeClassifier::OctreeClassifier( const OctreeClassifier*                           otherTree,
4728                                     const std::vector< ElementsOnShape::Classifier >& clsOther,
4729                                     std::vector< ElementsOnShape::Classifier >&       cls )
4730   :SMESH_Octree( new SMESH_TreeLimit )
4731 {
4732   myBox = new Bnd_B3d( *otherTree->getBox() );
4733
4734   if (( myIsLeaf = otherTree->isLeaf() ))
4735   {
4736     myClassifiers.resize( otherTree->myClassifiers.size() );
4737     for ( size_t i = 0; i < otherTree->myClassifiers.size(); ++i )
4738     {
4739       int ind = otherTree->myClassifiers[i] - & clsOther[0];
4740       myClassifiers[ i ] = & cls[ ind ];
4741     }
4742   }
4743   else if ( otherTree->myChildren )
4744   {
4745     myChildren = new SMESH_Tree< Bnd_B3d, 8 > * [ 8 ];
4746     for ( int i = 0; i < nbChildren(); i++ )
4747       myChildren[i] =
4748         new OctreeClassifier( static_cast<const OctreeClassifier*>( otherTree->myChildren[i]),
4749                               clsOther, cls );
4750   }
4751 }
4752
4753 void ElementsOnShape::
4754 OctreeClassifier::GetClassifiersAtPoint( const gp_XYZ& point,
4755                                          std::vector< ElementsOnShape::Classifier* >& result )
4756 {
4757   if ( getBox()->IsOut( point ))
4758     return;
4759
4760   if ( isLeaf() )
4761   {
4762     for ( size_t i = 0; i < myClassifiers.size(); ++i )
4763       if ( !myClassifiers[i]->GetBndBox()->IsOut( point ))
4764         result.push_back( myClassifiers[i] );
4765   }
4766   else
4767   {
4768     for (int i = 0; i < nbChildren(); i++)
4769       ((OctreeClassifier*) myChildren[i])->GetClassifiersAtPoint( point, result );
4770   }
4771 }
4772
4773 void ElementsOnShape::OctreeClassifier::buildChildrenData()
4774 {
4775   // distribute myClassifiers among myChildren
4776
4777   const int childFlag[8] = { 0x0000001,
4778                              0x0000002,
4779                              0x0000004,
4780                              0x0000008,
4781                              0x0000010,
4782                              0x0000020,
4783                              0x0000040,
4784                              0x0000080 };
4785   int nbInChild[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
4786
4787   for ( size_t i = 0; i < myClassifiers.size(); ++i )
4788   {
4789     for ( int j = 0; j < nbChildren(); j++ )
4790     {
4791       if ( !myClassifiers[i]->GetBndBox()->IsOut( *myChildren[j]->getBox() ))
4792       {
4793         myClassifiers[i]->SetFlag( childFlag[ j ]);
4794         ++nbInChild[ j ];
4795       }
4796     }
4797   }
4798
4799   for ( int j = 0; j < nbChildren(); j++ )
4800   {
4801     OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ j ]);
4802     child->myClassifiers.resize( nbInChild[ j ]);
4803     for ( size_t i = 0; nbInChild[ j ] && i < myClassifiers.size(); ++i )
4804     {
4805       if ( myClassifiers[ i ]->IsSetFlag( childFlag[ j ]))
4806       {
4807         --nbInChild[ j ];
4808         child->myClassifiers[ nbInChild[ j ]] = myClassifiers[ i ];
4809         myClassifiers[ i ]->UnsetFlag( childFlag[ j ]);
4810       }
4811     }
4812   }
4813   SMESHUtils::FreeVector( myClassifiers );
4814
4815   // define if a child isLeaf()
4816   for ( int i = 0; i < nbChildren(); i++ )
4817   {
4818     OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ i ]);
4819     child->myIsLeaf = ( child->myClassifiers.size() <= 5  );
4820   }
4821 }
4822
4823 Bnd_B3d* ElementsOnShape::OctreeClassifier::buildRootBox()
4824 {
4825   Bnd_B3d* box = new Bnd_B3d;
4826   for ( size_t i = 0; i < myClassifiers.size(); ++i )
4827     box->Add( *myClassifiers[i]->GetBndBox() );
4828   return box;
4829 }
4830
4831 /*
4832   Class       : BelongToGeom
4833   Description : Predicate for verifying whether entity belongs to
4834                 specified geometrical support
4835 */
4836
4837 BelongToGeom::BelongToGeom()
4838   : myMeshDS(NULL),
4839     myType(SMDSAbs_NbElementTypes),
4840     myIsSubshape(false),
4841     myTolerance(Precision::Confusion())
4842 {}
4843
4844 Predicate* BelongToGeom::clone() const
4845 {
4846   BelongToGeom* cln = new BelongToGeom( *this );
4847   cln->myElementsOnShapePtr.reset( static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ));
4848   return cln;
4849 }
4850
4851 void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
4852 {
4853   if ( myMeshDS != theMesh )
4854   {
4855     myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
4856     init();
4857   }
4858 }
4859
4860 void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
4861 {
4862   if ( myShape != theShape )
4863   {
4864     myShape = theShape;
4865     init();
4866   }
4867 }
4868
4869 static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
4870                         const TopoDS_Shape&               theShape)
4871 {
4872   if (theMap.Contains(theShape)) return true;
4873
4874   if (theShape.ShapeType() == TopAbs_COMPOUND ||
4875       theShape.ShapeType() == TopAbs_COMPSOLID)
4876   {
4877     TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
4878     for (; anIt.More(); anIt.Next())
4879     {
4880       if (!IsSubShape(theMap, anIt.Value())) {
4881         return false;
4882       }
4883     }
4884     return true;
4885   }
4886
4887   return false;
4888 }
4889
4890 void BelongToGeom::init()
4891 {
4892   if ( !myMeshDS || myShape.IsNull() ) return;
4893
4894   // is sub-shape of main shape?
4895   TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
4896   if (aMainShape.IsNull()) {
4897     myIsSubshape = false;
4898   }
4899   else {
4900     TopTools_IndexedMapOfShape aMap;
4901     TopExp::MapShapes( aMainShape, aMap );
4902     myIsSubshape = IsSubShape( aMap, myShape );
4903     if ( myIsSubshape )
4904     {
4905       aMap.Clear();
4906       TopExp::MapShapes( myShape, aMap );
4907       mySubShapesIDs.Clear();
4908       for ( int i = 1; i <= aMap.Extent(); ++i )
4909       {
4910         int subID = myMeshDS->ShapeToIndex( aMap( i ));
4911         if ( subID > 0 )
4912           mySubShapesIDs.Add( subID );
4913       }
4914     }
4915   }
4916
4917   //if (!myIsSubshape) // to be always ready to check an element not bound to geometry
4918   {
4919     if ( !myElementsOnShapePtr )
4920       myElementsOnShapePtr.reset( new ElementsOnShape() );
4921     myElementsOnShapePtr->SetTolerance( myTolerance );
4922     myElementsOnShapePtr->SetAllNodes( true ); // "belong", while false means "lays on"
4923     myElementsOnShapePtr->SetMesh( myMeshDS );
4924     myElementsOnShapePtr->SetShape( myShape, myType );
4925   }
4926 }
4927
4928 bool BelongToGeom::IsSatisfy (long theId)
4929 {
4930   if (myMeshDS == 0 || myShape.IsNull())
4931     return false;
4932
4933   if (!myIsSubshape)
4934   {
4935     return myElementsOnShapePtr->IsSatisfy(theId);
4936   }
4937
4938   // Case of sub-mesh
4939
4940   if (myType == SMDSAbs_Node)
4941   {
4942     if ( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ))
4943     {
4944       if ( aNode->getshapeId() < 1 )
4945         return myElementsOnShapePtr->IsSatisfy(theId);
4946       else
4947         return mySubShapesIDs.Contains( aNode->getshapeId() );
4948     }
4949   }
4950   else
4951   {
4952     if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
4953     {
4954       if ( myType == SMDSAbs_All || anElem->GetType() == myType )
4955       {
4956         if ( anElem->getshapeId() < 1 )
4957           return myElementsOnShapePtr->IsSatisfy(theId);
4958         else
4959           return mySubShapesIDs.Contains( anElem->getshapeId() );
4960       }
4961     }
4962   }
4963
4964   return false;
4965 }
4966
4967 void BelongToGeom::SetType (SMDSAbs_ElementType theType)
4968 {
4969   if ( myType != theType )
4970   {
4971     myType = theType;
4972     init();
4973   }
4974 }
4975
4976 SMDSAbs_ElementType BelongToGeom::GetType() const
4977 {
4978   return myType;
4979 }
4980
4981 TopoDS_Shape BelongToGeom::GetShape()
4982 {
4983   return myShape;
4984 }
4985
4986 const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const
4987 {
4988   return myMeshDS;
4989 }
4990
4991 void BelongToGeom::SetTolerance (double theTolerance)
4992 {
4993   myTolerance = theTolerance;
4994   init();
4995 }
4996
4997 double BelongToGeom::GetTolerance()
4998 {
4999   return myTolerance;
5000 }
5001
5002 /*
5003   Class       : LyingOnGeom
5004   Description : Predicate for verifying whether entiy lying or partially lying on
5005   specified geometrical support
5006 */
5007
5008 LyingOnGeom::LyingOnGeom()
5009   : myMeshDS(NULL),
5010     myType(SMDSAbs_NbElementTypes),
5011     myIsSubshape(false),
5012     myTolerance(Precision::Confusion())
5013 {}
5014
5015 Predicate* LyingOnGeom::clone() const
5016 {
5017   LyingOnGeom* cln = new LyingOnGeom( *this );
5018   cln->myElementsOnShapePtr.reset( static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ));
5019   return cln;
5020 }
5021
5022 void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
5023 {
5024   if ( myMeshDS != theMesh )
5025   {
5026     myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
5027     init();
5028   }
5029 }
5030
5031 void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
5032 {
5033   if ( myShape != theShape )
5034   {
5035     myShape = theShape;
5036     init();
5037   }
5038 }
5039
5040 void LyingOnGeom::init()
5041 {
5042   if (!myMeshDS || myShape.IsNull()) return;
5043
5044   // is sub-shape of main shape?
5045   TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
5046   if (aMainShape.IsNull()) {
5047     myIsSubshape = false;
5048   }
5049   else {
5050     myIsSubshape = myMeshDS->IsGroupOfSubShapes( myShape );
5051   }
5052
5053   if (myIsSubshape)
5054   {
5055     TopTools_IndexedMapOfShape shapes;
5056     TopExp::MapShapes( myShape, shapes );
5057     mySubShapesIDs.Clear();
5058     for ( int i = 1; i <= shapes.Extent(); ++i )
5059     {
5060       int subID = myMeshDS->ShapeToIndex( shapes( i ));
5061       if ( subID > 0 )
5062         mySubShapesIDs.Add( subID );
5063     }
5064   }
5065   // else // to be always ready to check an element not bound to geometry
5066   {
5067     if ( !myElementsOnShapePtr )
5068       myElementsOnShapePtr.reset( new ElementsOnShape() );
5069     myElementsOnShapePtr->SetTolerance( myTolerance );
5070     myElementsOnShapePtr->SetAllNodes( false ); // lays on, while true means "belong"
5071     myElementsOnShapePtr->SetMesh( myMeshDS );
5072     myElementsOnShapePtr->SetShape( myShape, myType );
5073   }
5074 }
5075
5076 bool LyingOnGeom::IsSatisfy( long theId )
5077 {
5078   if ( myMeshDS == 0 || myShape.IsNull() )
5079     return false;
5080
5081   if (!myIsSubshape)
5082   {
5083     return myElementsOnShapePtr->IsSatisfy(theId);
5084   }
5085
5086   // Case of sub-mesh
5087
5088   const SMDS_MeshElement* elem =
5089     ( myType == SMDSAbs_Node ) ? myMeshDS->FindNode( theId ) : myMeshDS->FindElement( theId );
5090
5091   if ( mySubShapesIDs.Contains( elem->getshapeId() ))
5092     return true;
5093
5094   if (( elem->GetType() != SMDSAbs_Node ) &&
5095       ( myType == SMDSAbs_All || elem->GetType() == myType ))
5096   {
5097     SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator();
5098     while ( nodeItr->more() )
5099     {
5100       const SMDS_MeshElement* aNode = nodeItr->next();
5101       if ( mySubShapesIDs.Contains( aNode->getshapeId() ))
5102         return true;
5103     }
5104   }
5105
5106   return false;
5107 }
5108
5109 void LyingOnGeom::SetType( SMDSAbs_ElementType theType )
5110 {
5111   if ( myType != theType )
5112   {
5113     myType = theType;
5114     init();
5115   }
5116 }
5117
5118 SMDSAbs_ElementType LyingOnGeom::GetType() const
5119 {
5120   return myType;
5121 }
5122
5123 TopoDS_Shape LyingOnGeom::GetShape()
5124 {
5125   return myShape;
5126 }
5127
5128 const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const
5129 {
5130   return myMeshDS;
5131 }
5132
5133 void LyingOnGeom::SetTolerance (double theTolerance)
5134 {
5135   myTolerance = theTolerance;
5136   init();
5137 }
5138
5139 double LyingOnGeom::GetTolerance()
5140 {
5141   return myTolerance;
5142 }
5143
5144 TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0)
5145 {}
5146
5147 TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n), myElem(0)
5148 {}
5149
5150 TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t), myElem(0)
5151 {}
5152
5153 TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray), myElem(theSequenceOfXYZ.myElem)
5154 {}
5155
5156 template <class InputIterator>
5157 TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd), myElem(0)
5158 {}
5159
5160 TSequenceOfXYZ::~TSequenceOfXYZ()
5161 {}
5162
5163 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
5164 {
5165   myArray = theSequenceOfXYZ.myArray;
5166   myElem  = theSequenceOfXYZ.myElem;
5167   return *this;
5168 }
5169
5170 gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
5171 {
5172   return myArray[n-1];
5173 }
5174
5175 const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
5176 {
5177   return myArray[n-1];
5178 }
5179
5180 void TSequenceOfXYZ::clear()
5181 {
5182   myArray.clear();
5183 }
5184
5185 void TSequenceOfXYZ::reserve(size_type n)
5186 {
5187   myArray.reserve(n);
5188 }
5189
5190 void TSequenceOfXYZ::push_back(const gp_XYZ& v)
5191 {
5192   myArray.push_back(v);
5193 }
5194
5195 TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
5196 {
5197   return myArray.size();
5198 }
5199
5200 SMDSAbs_EntityType TSequenceOfXYZ::getElementEntity() const
5201 {
5202   return myElem ? myElem->GetEntityType() : SMDSEntity_Last;
5203 }
5204
5205 TMeshModifTracer::TMeshModifTracer():
5206   myMeshModifTime(0), myMesh(0)
5207 {
5208 }
5209 void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh )
5210 {
5211   if ( theMesh != myMesh )
5212     myMeshModifTime = 0;
5213   myMesh = theMesh;
5214 }
5215 bool TMeshModifTracer::IsMeshModified()
5216 {
5217   bool modified = false;
5218   if ( myMesh )
5219   {
5220     modified = ( myMeshModifTime != myMesh->GetMTime() );
5221     myMeshModifTime = myMesh->GetMTime();
5222   }
5223   return modified;
5224 }