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