src/Controls/SMESH_Controls.cxx

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