src/Controls/SMESH_Controls.cxx

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