src/Controls/SMESH_Controls.cxx

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