src/Controls/SMESH_Controls.cxx

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