src/Controls/SMESH_Controls.cxx

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