src/Controls/SMESH_Controls.cxx

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