src/Controls/SMESH_Controls.cxx

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