src/Controls/SMESH_Controls.cxx

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