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