src/Controls/SMESH_Controls.cxx

   1 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 //  This library is free software; you can redistribute it and/or
   5 //  modify it under the terms of the GNU Lesser General Public
   6 //  License as published by the Free Software Foundation; either
   7 //  version 2.1 of the License.
   8 //
   9 //  This library is distributed in the hope that it will be useful,
  10 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 //  Lesser General Public License for more details.
  13 //
  14 //  You should have received a copy of the GNU Lesser General Public
  15 //  License along with this library; if not, write to the Free Software
  16 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  19
  20 #include "SMESH_ControlsDef.hxx"
  21
  22 #include <set>
  23
  24 #include <BRep_Tool.hxx>
  25 #include <gp_Ax3.hxx>
  26 #include <gp_Cylinder.hxx>
  27 #include <gp_Dir.hxx>
  28 #include <gp_Pnt.hxx>
  29 #include <gp_Pln.hxx>
  30 #include <gp_Vec.hxx>
  31 #include <gp_XYZ.hxx>
  32 #include <Geom_Plane.hxx>
  33 #include <Geom_CylindricalSurface.hxx>
  34 #include <Precision.hxx>
  35 #include <TColgp_Array1OfXYZ.hxx>
  36 #include <TColStd_MapOfInteger.hxx>
  37 #include <TColStd_SequenceOfAsciiString.hxx>
  38 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
  39 #include <TopAbs.hxx>
  40 #include <TopoDS.hxx>
  41 #include <TopoDS_Face.hxx>
  42 #include <TopoDS_Shape.hxx>
  43
  44 #include "SMDS_Mesh.hxx"
  45 #include "SMDS_Iterator.hxx"
  46 #include "SMDS_MeshElement.hxx"
  47 #include "SMDS_MeshNode.hxx"
  48
  49
  50 /*
  51                             AUXILIARY METHODS
  52 */
  53
  54 namespace{
  55   inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
  56   {
  57     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
  58
  59     return v1.Magnitude() < gp::Resolution() ||
  60       v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
  61   }
  62
  63   inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
  64   {
  65     gp_Vec aVec1( P2 - P1 );
  66     gp_Vec aVec2( P3 - P1 );
  67     return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
  68   }
  69
  70   inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
  71   {
  72     return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
  73   }
  74
  75
  76
  77   inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
  78   {
  79     double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
  80     return aDist;
  81   }
  82
  83   int getNbMultiConnection( SMDS_Mesh* theMesh, const int theId )
  84   {
  85     if ( theMesh == 0 )
  86       return 0;
  87
  88     const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
  89     if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 )
  90       return 0;
  91
  92     TColStd_MapOfInteger aMap;
  93
  94     int aResult = 0;
  95     SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
  96     if ( anIter != 0 ) {
  97       while( anIter->more() ) {
  98         const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
  99         if ( aNode == 0 )
 100           return 0;
 101         SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
 102         while( anElemIter->more() ) {
 103           const SMDS_MeshElement* anElem = anElemIter->next();
 104           if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
 105             int anId = anElem->GetID();
 106
 107             if ( anIter->more() )              // i.e. first node
 108               aMap.Add( anId );
 109             else if ( aMap.Contains( anId ) )
 110               aResult++;
 111           }
 112         }
 113       }
 114     }
 115
 116     return aResult;
 117   }
 118
 119 }
 120
 121
 122
 123 using namespace SMESH::Controls;
 124
 125 /*
 126                                 FUNCTORS
 127 */
 128
 129 /*
 130   Class       : NumericalFunctor
 131   Description : Base class for numerical functors
 132 */
 133 NumericalFunctor::NumericalFunctor():
 134   myMesh(NULL)
 135 {
 136   myPrecision = -1;
 137 }
 138
 139 void NumericalFunctor::SetMesh( SMDS_Mesh* theMesh )
 140 {
 141   myMesh = theMesh;
 142 }
 143
 144 bool NumericalFunctor::GetPoints(const int theId,
 145                                  TSequenceOfXYZ& theRes ) const
 146 {
 147   theRes.clear();
 148
 149   if ( myMesh == 0 )
 150     return false;
 151
 152   return GetPoints( myMesh->FindElement( theId ), theRes );
 153 }
 154
 155 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
 156                                  TSequenceOfXYZ& theRes )
 157 {
 158   theRes.clear();
 159
 160   if ( anElem == 0)
 161     return false;
 162
 163   // Get nodes of the element
 164   SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
 165   if ( anIter != 0 )
 166   {
 167     while( anIter->more() )
 168     {
 169       const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
 170       if ( aNode != 0 ){
 171         theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
 172       }
 173     }
 174   }
 175
 176   return true;
 177 }
 178
 179 long  NumericalFunctor::GetPrecision() const
 180 {
 181   return myPrecision;
 182 }
 183
 184 void  NumericalFunctor::SetPrecision( const long thePrecision )
 185 {
 186   myPrecision = thePrecision;
 187 }
 188
 189 double NumericalFunctor::GetValue( long theId )
 190 {
 191   TSequenceOfXYZ P;
 192   if ( GetPoints( theId, P ))
 193   {
 194     double aVal = GetValue( P );
 195     if ( myPrecision >= 0 )
 196     {
 197       double prec = pow( 10., (double)( myPrecision ) );
 198       aVal = floor( aVal * prec + 0.5 ) / prec;
 199     }
 200     return aVal;
 201   }
 202
 203   return 0.;
 204 }
 205
 206 /*
 207   Class       : MinimumAngle
 208   Description : Functor for calculation of minimum angle
 209 */
 210
 211 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
 212 {
 213   double aMin;
 214
 215   if ( P.size() == 3 )
 216   {
 217     double A0 = getAngle( P( 3 ), P( 1 ), P( 2 ) );
 218     double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
 219     double A2 = getAngle( P( 2 ), P( 3 ), P( 1 ) );
 220
 221     aMin = Min( A0, Min( A1, A2 ) );
 222   }
 223   else if ( P.size() == 4 )
 224   {
 225     double A0 = getAngle( P( 4 ), P( 1 ), P( 2 ) );
 226     double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
 227     double A2 = getAngle( P( 2 ), P( 3 ), P( 4 ) );
 228     double A3 = getAngle( P( 3 ), P( 4 ), P( 1 ) );
 229
 230     aMin = Min( Min( A0, A1 ), Min( A2, A3 ) );
 231   }
 232   else
 233     return 0.;
 234
 235   return aMin * 180 / PI;
 236 }
 237
 238 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
 239 {
 240   const double aBestAngle = PI / nbNodes;
 241   return ( fabs( aBestAngle - Value ));
 242 }
 243
 244 SMDSAbs_ElementType MinimumAngle::GetType() const
 245 {
 246   return SMDSAbs_Face;
 247 }
 248
 249
 250 /*
 251   Class       : AspectRatio
 252   Description : Functor for calculating aspect ratio
 253 */
 254 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
 255 {
 256   int nbNodes = P.size();
 257
 258   if ( nbNodes != 3 && nbNodes != 4 )
 259     return 0;
 260
 261   // Compute lengths of the sides
 262
 263   double aLen[ nbNodes ];
 264   for ( int i = 0; i < nbNodes - 1; i++ )
 265     aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
 266   aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
 267
 268   // Compute aspect ratio
 269
 270   if ( nbNodes == 3 )
 271   {
 272     double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
 273     if ( anArea <= Precision::Confusion() )
 274       return 0.;
 275     double aMaxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
 276     static double aCoef = sqrt( 3. ) / 4;
 277
 278     return aCoef * aMaxLen * aMaxLen / anArea;
 279   }
 280   else
 281   {
 282     double aMinLen = Min( Min( aLen[ 0 ], aLen[ 1 ] ), Min( aLen[ 2 ], aLen[ 3 ] ) );
 283     if ( aMinLen <= Precision::Confusion() )
 284       return 0.;
 285     double aMaxLen = Max( Max( aLen[ 0 ], aLen[ 1 ] ), Max( aLen[ 2 ], aLen[ 3 ] ) );
 286
 287     return aMaxLen / aMinLen;
 288   }
 289 }
 290
 291 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
 292 {
 293   // the aspect ratio is in the range [1.0,infinity]
 294   // 1.0 = good
 295   // infinity = bad
 296   return Value / 1000.;
 297 }
 298
 299 SMDSAbs_ElementType AspectRatio::GetType() const
 300 {
 301   return SMDSAbs_Face;
 302 }
 303
 304
 305 /*
 306   Class       : AspectRatio3D
 307   Description : Functor for calculating aspect ratio
 308 */
 309 namespace{
 310
 311   inline double getHalfPerimeter(double theTria[3]){
 312     return (theTria[0] + theTria[1] + theTria[2])/2.0;
 313   }
 314
 315   inline double getArea(double theHalfPerim, double theTria[3]){
 316     return sqrt(theHalfPerim*
 317                 (theHalfPerim-theTria[0])*
 318                 (theHalfPerim-theTria[1])*
 319                 (theHalfPerim-theTria[2]));
 320   }
 321
 322   inline double getVolume(double theLen[6]){
 323     double a2 = theLen[0]*theLen[0];
 324     double b2 = theLen[1]*theLen[1];
 325     double c2 = theLen[2]*theLen[2];
 326     double d2 = theLen[3]*theLen[3];
 327     double e2 = theLen[4]*theLen[4];
 328     double f2 = theLen[5]*theLen[5];
 329     double P = 4.0*a2*b2*d2;
 330     double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
 331     double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
 332     return sqrt(P-Q+R)/12.0;
 333   }
 334
 335   inline double getVolume2(double theLen[6]){
 336     double a2 = theLen[0]*theLen[0];
 337     double b2 = theLen[1]*theLen[1];
 338     double c2 = theLen[2]*theLen[2];
 339     double d2 = theLen[3]*theLen[3];
 340     double e2 = theLen[4]*theLen[4];
 341     double f2 = theLen[5]*theLen[5];
 342
 343     double P = a2*e2*(b2+c2+d2+f2-a2-e2);
 344     double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
 345     double R = c2*d2*(a2+b2+e2+f2-c2-d2);
 346     double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
 347
 348     return sqrt(P+Q+R-S)/12.0;
 349   }
 350
 351   inline double getVolume(const TSequenceOfXYZ& P){
 352     gp_Vec aVec1( P( 2 ) - P( 1 ) );
 353     gp_Vec aVec2( P( 3 ) - P( 1 ) );
 354     gp_Vec aVec3( P( 4 ) - P( 1 ) );
 355     gp_Vec anAreaVec( aVec1 ^ aVec2 );
 356     return abs(aVec3 * anAreaVec) / 6.0;
 357   }
 358
 359   inline double getMaxHeight(double theLen[6])
 360   {
 361     double aHeight = max(theLen[0],theLen[1]);
 362     aHeight = max(aHeight,theLen[2]);
 363     aHeight = max(aHeight,theLen[3]);
 364     aHeight = max(aHeight,theLen[4]);
 365     aHeight = max(aHeight,theLen[5]);
 366     return aHeight;
 367   }
 368
 369 }
 370
 371 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
 372 {
 373   double aQuality = 0.0;
 374   int nbNodes = P.size();
 375   switch(nbNodes){
 376   case 4:{
 377     double aLen[6] = {
 378       getDistance(P( 1 ),P( 2 )), // a
 379       getDistance(P( 2 ),P( 3 )), // b
 380       getDistance(P( 3 ),P( 1 )), // c
 381       getDistance(P( 2 ),P( 4 )), // d
 382       getDistance(P( 3 ),P( 4 )), // e
 383       getDistance(P( 1 ),P( 4 ))  // f
 384     };
 385     double aTria[4][3] = {
 386       {aLen[0],aLen[1],aLen[2]}, // abc
 387       {aLen[0],aLen[3],aLen[5]}, // adf
 388       {aLen[1],aLen[3],aLen[4]}, // bde
 389       {aLen[2],aLen[4],aLen[5]}  // cef
 390     };
 391     double aSumArea = 0.0;
 392     double aHalfPerimeter = getHalfPerimeter(aTria[0]);
 393     double anArea = getArea(aHalfPerimeter,aTria[0]);
 394     aSumArea += anArea;
 395     aHalfPerimeter = getHalfPerimeter(aTria[1]);
 396     anArea = getArea(aHalfPerimeter,aTria[1]);
 397     aSumArea += anArea;
 398     aHalfPerimeter = getHalfPerimeter(aTria[2]);
 399     anArea = getArea(aHalfPerimeter,aTria[2]);
 400     aSumArea += anArea;
 401     aHalfPerimeter = getHalfPerimeter(aTria[3]);
 402     anArea = getArea(aHalfPerimeter,aTria[3]);
 403     aSumArea += anArea;
 404     double aVolume = getVolume(P);
 405     //double aVolume = getVolume(aLen);
 406     double aHeight = getMaxHeight(aLen);
 407     static double aCoeff = sqrt(6.0)/36.0;
 408     aQuality = aCoeff*aHeight*aSumArea/aVolume;
 409     break;
 410   }
 411   case 5:{
 412     {
 413       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
 414       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 415     }
 416     {
 417       gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
 418       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 419     }
 420     {
 421       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
 422       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 423     }
 424     {
 425       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
 426       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 427     }
 428     break;
 429   }
 430   case 6:{
 431     {
 432       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
 433       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 434     }
 435     {
 436       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
 437       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 438     }
 439     {
 440       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
 441       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 442     }
 443     {
 444       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
 445       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 446     }
 447     {
 448       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
 449       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 450     }
 451     {
 452       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
 453       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 454     }
 455     break;
 456   }
 457   case 8:{
 458     {
 459       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
 460       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 461     }
 462     {
 463       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
 464       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 465     }
 466     {
 467       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
 468       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 469     }
 470     {
 471       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
 472       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 473     }
 474     {
 475       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
 476       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 477     }
 478     {
 479       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
 480       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 481     }
 482     {
 483       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
 484       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 485     }
 486     {
 487       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
 488       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 489     }
 490     {
 491       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
 492       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 493     }
 494     {
 495       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
 496       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 497     }
 498     {
 499       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
 500       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 501     }
 502     {
 503       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
 504       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 505     }
 506     {
 507       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
 508       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 509     }
 510     {
 511       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
 512       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 513     }
 514     {
 515       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
 516       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 517     }
 518     {
 519       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
 520       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 521     }
 522     {
 523       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
 524       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 525     }
 526     {
 527       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
 528       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 529     }
 530     {
 531       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
 532       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 533     }
 534     {
 535       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
 536       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 537     }
 538     {
 539       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
 540       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 541     }
 542     {
 543       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
 544       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 545     }
 546     {
 547       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
 548       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 549     }
 550     {
 551       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
 552       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 553     }
 554     {
 555       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
 556       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 557     }
 558     {
 559       gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
 560       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 561     }
 562     {
 563       gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
 564       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 565     }
 566     {
 567       gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
 568       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 569     }
 570     {
 571       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
 572       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 573     }
 574     {
 575       gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
 576       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 577     }
 578     {
 579       gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
 580       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 581     }
 582     {
 583       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
 584       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 585     }
 586     {
 587       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
 588       aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 589     }
 590     break;
 591   }
 592   }
 593   return aQuality;
 594 }
 595
 596 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
 597 {
 598   // the aspect ratio is in the range [1.0,infinity]
 599   // 1.0 = good
 600   // infinity = bad
 601   return Value / 1000.;
 602 }
 603
 604 SMDSAbs_ElementType AspectRatio3D::GetType() const
 605 {
 606   return SMDSAbs_Volume;
 607 }
 608
 609
 610 /*
 611   Class       : Warping
 612   Description : Functor for calculating warping
 613 */
 614 double Warping::GetValue( const TSequenceOfXYZ& P )
 615 {
 616   if ( P.size() != 4 )
 617     return 0;
 618
 619   gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4;
 620
 621   double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
 622   double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
 623   double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
 624   double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
 625
 626   return Max( Max( A1, A2 ), Max( A3, A4 ) );
 627 }
 628
 629 double Warping::ComputeA( const gp_XYZ& thePnt1,
 630                           const gp_XYZ& thePnt2,
 631                           const gp_XYZ& thePnt3,
 632                           const gp_XYZ& theG ) const
 633 {
 634   double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
 635   double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
 636   double L = Min( aLen1, aLen2 ) * 0.5;
 637   if ( L < Precision::Confusion())
 638     return 0.;
 639
 640   gp_XYZ GI = ( thePnt2 - thePnt1 ) / 2. - theG;
 641   gp_XYZ GJ = ( thePnt3 - thePnt2 ) / 2. - theG;
 642   gp_XYZ N  = GI.Crossed( GJ );
 643
 644   if ( N.Modulus() < gp::Resolution() )
 645     return PI / 2;
 646
 647   N.Normalize();
 648
 649   double H = ( thePnt2 - theG ).Dot( N );
 650   return asin( fabs( H / L ) ) * 180 / PI;
 651 }
 652
 653 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
 654 {
 655   // the warp is in the range [0.0,PI/2]
 656   // 0.0 = good (no warp)
 657   // PI/2 = bad  (face pliee)
 658   return Value;
 659 }
 660
 661 SMDSAbs_ElementType Warping::GetType() const
 662 {
 663   return SMDSAbs_Face;
 664 }
 665
 666
 667 /*
 668   Class       : Taper
 669   Description : Functor for calculating taper
 670 */
 671 double Taper::GetValue( const TSequenceOfXYZ& P )
 672 {
 673   if ( P.size() != 4 )
 674     return 0;
 675
 676   // Compute taper
 677   double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2;
 678   double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2;
 679   double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2;
 680   double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2;
 681
 682   double JA = 0.25 * ( J1 + J2 + J3 + J4 );
 683   if ( JA <= Precision::Confusion() )
 684     return 0.;
 685
 686   double T1 = fabs( ( J1 - JA ) / JA );
 687   double T2 = fabs( ( J2 - JA ) / JA );
 688   double T3 = fabs( ( J3 - JA ) / JA );
 689   double T4 = fabs( ( J4 - JA ) / JA );
 690
 691   return Max( Max( T1, T2 ), Max( T3, T4 ) );
 692 }
 693
 694 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
 695 {
 696   // the taper is in the range [0.0,1.0]
 697   // 0.0  = good (no taper)
 698   // 1.0 = bad  (les cotes opposes sont allignes)
 699   return Value;
 700 }
 701
 702 SMDSAbs_ElementType Taper::GetType() const
 703 {
 704   return SMDSAbs_Face;
 705 }
 706
 707
 708 /*
 709   Class       : Skew
 710   Description : Functor for calculating skew in degrees
 711 */
 712 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
 713 {
 714   gp_XYZ p12 = ( p2 + p1 ) / 2;
 715   gp_XYZ p23 = ( p3 + p2 ) / 2;
 716   gp_XYZ p31 = ( p3 + p1 ) / 2;
 717
 718   gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
 719
 720   return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
 721 }
 722
 723 double Skew::GetValue( const TSequenceOfXYZ& P )
 724 {
 725   if ( P.size() != 3 && P.size() != 4 )
 726     return 0;
 727
 728   // Compute skew
 729   static double PI2 = PI / 2;
 730   if ( P.size() == 3 )
 731   {
 732     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
 733     double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
 734     double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
 735
 736     return Max( A0, Max( A1, A2 ) ) * 180 / PI;
 737   }
 738   else
 739   {
 740     gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2;
 741     gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2;
 742     gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2;
 743     gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2;
 744
 745     gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
 746     double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
 747       ? 0 : fabs( PI2 - v1.Angle( v2 ) );
 748
 749     return A * 180 / PI;
 750   }
 751 }
 752
 753 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
 754 {
 755   // the skew is in the range [0.0,PI/2].
 756   // 0.0 = good
 757   // PI/2 = bad
 758   return Value;
 759 }
 760
 761 SMDSAbs_ElementType Skew::GetType() const
 762 {
 763   return SMDSAbs_Face;
 764 }
 765
 766
 767 /*
 768   Class       : Area
 769   Description : Functor for calculating area
 770 */
 771 double Area::GetValue( const TSequenceOfXYZ& P )
 772 {
 773   if ( P.size() == 3 )
 774     return getArea( P( 1 ), P( 2 ), P( 3 ) );
 775   else if ( P.size() == 4 )
 776     return getArea( P( 1 ), P( 2 ), P( 3 ) ) + getArea( P( 1 ), P( 3 ), P( 4 ) );
 777   else
 778     return 0;
 779 }
 780
 781 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
 782 {
 783   return Value;
 784 }
 785
 786 SMDSAbs_ElementType Area::GetType() const
 787 {
 788   return SMDSAbs_Face;
 789 }
 790
 791
 792 /*
 793   Class       : Length
 794   Description : Functor for calculating length off edge
 795 */
 796 double Length::GetValue( const TSequenceOfXYZ& P )
 797 {
 798   return ( P.size() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0 );
 799 }
 800
 801 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
 802 {
 803   return Value;
 804 }
 805
 806 SMDSAbs_ElementType Length::GetType() const
 807 {
 808   return SMDSAbs_Edge;
 809 }
 810
 811 /*
 812   Class       : Length2D
 813   Description : Functor for calculating length of edge
 814 */
 815
 816 double Length2D::GetValue( long theElementId)
 817 {
 818   TSequenceOfXYZ P;
 819
 820   if (GetPoints(theElementId,P)){
 821
 822     double  aVal;// = GetValue( P );
 823     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
 824     SMDSAbs_ElementType aType = aElem->GetType();
 825
 826     int len = P.size();
 827
 828     switch (aType){
 829     case SMDSAbs_All:
 830     case SMDSAbs_Node:
 831     case SMDSAbs_Edge:
 832       if (len == 2){
 833         aVal = getDistance( P( 1 ), P( 2 ) );
 834         break;
 835       }
 836     case SMDSAbs_Face:
 837       if (len == 3){ // triangles
 838         double L1 = getDistance(P( 1 ),P( 2 ));
 839         double L2 = getDistance(P( 2 ),P( 3 ));
 840         double L3 = getDistance(P( 3 ),P( 1 ));
 841         aVal = Max(L1,Max(L2,L3));
 842         break;
 843       }
 844       else if (len == 4){ // quadrangles
 845         double L1 = getDistance(P( 1 ),P( 2 ));
 846         double L2 = getDistance(P( 2 ),P( 3 ));
 847         double L3 = getDistance(P( 3 ),P( 4 ));
 848         double L4 = getDistance(P( 4 ),P( 1 ));
 849         aVal = Max(Max(L1,L2),Max(L3,L4));
 850         break;
 851       }
 852     case SMDSAbs_Volume:
 853       if (len == 4){ // tetraidrs
 854         double L1 = getDistance(P( 1 ),P( 2 ));
 855         double L2 = getDistance(P( 2 ),P( 3 ));
 856         double L3 = getDistance(P( 3 ),P( 1 ));
 857         double L4 = getDistance(P( 1 ),P( 4 ));
 858         double L5 = getDistance(P( 2 ),P( 4 ));
 859         double L6 = getDistance(P( 3 ),P( 4 ));
 860         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 861         break;
 862       }
 863       else if (len == 5){ // piramids
 864         double L1 = getDistance(P( 1 ),P( 2 ));
 865         double L2 = getDistance(P( 2 ),P( 3 ));
 866         double L3 = getDistance(P( 3 ),P( 1 ));
 867         double L4 = getDistance(P( 4 ),P( 1 ));
 868         double L5 = getDistance(P( 1 ),P( 5 ));
 869         double L6 = getDistance(P( 2 ),P( 5 ));
 870         double L7 = getDistance(P( 3 ),P( 5 ));
 871         double L8 = getDistance(P( 4 ),P( 5 ));
 872
 873         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 874         aVal = Max(aVal,Max(L7,L8));
 875         break;
 876       }
 877       else if (len == 6){ // pentaidres
 878         double L1 = getDistance(P( 1 ),P( 2 ));
 879         double L2 = getDistance(P( 2 ),P( 3 ));
 880         double L3 = getDistance(P( 3 ),P( 1 ));
 881         double L4 = getDistance(P( 4 ),P( 5 ));
 882         double L5 = getDistance(P( 5 ),P( 6 ));
 883         double L6 = getDistance(P( 6 ),P( 4 ));
 884         double L7 = getDistance(P( 1 ),P( 4 ));
 885         double L8 = getDistance(P( 2 ),P( 5 ));
 886         double L9 = getDistance(P( 3 ),P( 6 ));
 887
 888         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 889         aVal = Max(aVal,Max(Max(L7,L8),L9));
 890         break;
 891       }
 892       else if (len == 8){ // hexaider
 893         double L1 = getDistance(P( 1 ),P( 2 ));
 894         double L2 = getDistance(P( 2 ),P( 3 ));
 895         double L3 = getDistance(P( 3 ),P( 4 ));
 896         double L4 = getDistance(P( 4 ),P( 1 ));
 897         double L5 = getDistance(P( 5 ),P( 6 ));
 898         double L6 = getDistance(P( 6 ),P( 7 ));
 899         double L7 = getDistance(P( 7 ),P( 8 ));
 900         double L8 = getDistance(P( 8 ),P( 5 ));
 901         double L9 = getDistance(P( 1 ),P( 5 ));
 902         double L10= getDistance(P( 2 ),P( 6 ));
 903         double L11= getDistance(P( 3 ),P( 7 ));
 904         double L12= getDistance(P( 4 ),P( 8 ));
 905
 906         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 907         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
 908         aVal = Max(aVal,Max(L11,L12));
 909         break;
 910
 911       }
 912
 913     default: aVal=-1;
 914     }
 915
 916     if (aVal <0){
 917       return 0.;
 918     }
 919
 920     if ( myPrecision >= 0 )
 921     {
 922       double prec = pow( 10., (double)( myPrecision ) );
 923       aVal = floor( aVal * prec + 0.5 ) / prec;
 924     }
 925
 926     return aVal;
 927
 928   }
 929   return 0.;
 930 }
 931
 932 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
 933 {
 934   return Value;
 935 }
 936
 937 SMDSAbs_ElementType Length2D::GetType() const
 938 {
 939   return SMDSAbs_Face;
 940 }
 941
 942 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
 943   myLength(theLength)
 944 {
 945   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
 946   if(thePntId1 > thePntId2){
 947     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
 948   }
 949 }
 950
 951 bool Length2D::Value::operator<(const Length2D::Value& x) const{
 952   if(myPntId[0] < x.myPntId[0]) return true;
 953   if(myPntId[0] == x.myPntId[0])
 954     if(myPntId[1] < x.myPntId[1]) return true;
 955   return false;
 956 }
 957
 958 void Length2D::GetValues(TValues& theValues){
 959   TValues aValues;
 960   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
 961   int i = 0;
 962   for(; anIter->more(); ){
 963     const SMDS_MeshFace* anElem = anIter->next();
 964     long anElemId = anElem->GetID();
 965     SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
 966     long aNodeId[2];
 967     gp_Pnt P[3];
 968
 969     double aLength;
 970     const SMDS_MeshElement* aNode;
 971     if(aNodesIter->more()){
 972       aNode = aNodesIter->next();
 973       const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
 974       P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
 975       aNodeId[0] = aNodeId[1] = aNode->GetID();
 976       aLength = 0;
 977     }
 978     for(; aNodesIter->more(); ){
 979       aNode = aNodesIter->next();
 980       const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
 981       long anId = aNode->GetID();
 982
 983       P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
 984
 985       aLength = P[1].Distance(P[2]);
 986
 987       Value aValue(aLength,aNodeId[1],anId);
 988       aNodeId[1] = anId;
 989       P[1] = P[2];
 990       theValues.insert(aValue);
 991     }
 992
 993     aLength = P[0].Distance(P[1]);
 994
 995     Value aValue(aLength,aNodeId[0],aNodeId[1]);
 996     theValues.insert(aValue);
 997   }
 998 }
 999
1000 /*
1001   Class       : MultiConnection
1002   Description : Functor for calculating number of faces conneted to the edge
1003 */
1004 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
1005 {
1006   return 0;
1007 }
1008 double MultiConnection::GetValue( long theId )
1009 {
1010   return getNbMultiConnection( myMesh, theId );
1011 }
1012
1013 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
1014 {
1015   return Value;
1016 }
1017
1018 SMDSAbs_ElementType MultiConnection::GetType() const
1019 {
1020   return SMDSAbs_Edge;
1021 }
1022
1023 /*
1024   Class       : MultiConnection2D
1025   Description : Functor for calculating number of faces conneted to the edge
1026 */
1027 double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
1028 {
1029   return 0;
1030 }
1031
1032 double MultiConnection2D::GetValue( long theElementId )
1033 {
1034   TSequenceOfXYZ P;
1035   int aResult = 0;
1036
1037   if (GetPoints(theElementId,P)){
1038     double  aVal;
1039     const SMDS_MeshElement* anFaceElem = myMesh->FindElement( theElementId );
1040     SMDSAbs_ElementType aType = anFaceElem->GetType();
1041
1042     int len = P.size();
1043
1044     TColStd_MapOfInteger aMap;
1045     int aResult = 0;
1046
1047     switch (aType){
1048     case SMDSAbs_All:
1049     case SMDSAbs_Node:
1050     case SMDSAbs_Edge:
1051     case SMDSAbs_Face:
1052       if (len == 3){ // triangles
1053         int Nb[3] = {0,0,0};
1054
1055         int i=0;
1056         SMDS_ElemIteratorPtr anIter = anFaceElem->nodesIterator();
1057         if ( anIter != 0 ) {
1058           while( anIter->more() ) {
1059             const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
1060             if ( aNode == 0 ){
1061               break;
1062             }
1063             SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
1064             while( anElemIter->more() ) {
1065               const SMDS_MeshElement* anElem = anElemIter->next();
1066               if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
1067                 int anId = anElem->GetID();
1068
1069                 if ( anIter->more() )              // i.e. first node
1070                   aMap.Add( anId );
1071                 else if ( aMap.Contains( anId ) ){
1072                   Nb[i]++;
1073                 }
1074               }
1075               else if ( anElem != 0 && anElem->GetType() == SMDSAbs_Edge ) i++;
1076             }
1077           }
1078         }
1079
1080         aResult = Max(Max(Nb[0],Nb[1]),Nb[2]);
1081       }
1082       break;
1083     case SMDSAbs_Volume:
1084     default: aResult=0;
1085     }
1086
1087   }
1088   return aResult;//getNbMultiConnection( myMesh, theId );
1089 }
1090
1091 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1092 {
1093   return Value;
1094 }
1095
1096 SMDSAbs_ElementType MultiConnection2D::GetType() const
1097 {
1098   return SMDSAbs_Face;
1099 }
1100
1101 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
1102 {
1103   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
1104   if(thePntId1 > thePntId2){
1105     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
1106   }
1107 }
1108
1109 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{
1110   if(myPntId[0] < x.myPntId[0]) return true;
1111   if(myPntId[0] == x.myPntId[0])
1112     if(myPntId[1] < x.myPntId[1]) return true;
1113   return false;
1114 }
1115
1116 void MultiConnection2D::GetValues(MValues& theValues){
1117   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1118   for(; anIter->more(); ){
1119     const SMDS_MeshFace* anElem = anIter->next();
1120     long anElemId = anElem->GetID();
1121     SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1122     long aNodeId[3];
1123
1124     //int aNbConnects=0;
1125     const SMDS_MeshNode* aNode0;
1126     const SMDS_MeshNode* aNode1;
1127     const SMDS_MeshNode* aNode2;
1128     if(aNodesIter->more()){
1129       aNode0 = (SMDS_MeshNode*) aNodesIter->next();
1130       aNode1 = aNode0;
1131       const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
1132       aNodeId[0] = aNodeId[1] = aNodes->GetID();
1133     }
1134     for(; aNodesIter->more(); ){
1135       aNode2 = (SMDS_MeshNode*) aNodesIter->next();
1136       long anId = aNode2->GetID();
1137       aNodeId[2] = anId;
1138
1139       Value aValue(aNodeId[1],aNodeId[2]);
1140       MValues::iterator aItr = theValues.find(aValue);
1141       if (aItr != theValues.end()){
1142         aItr->second += 1;
1143         //aNbConnects = nb;
1144       } else {
1145         theValues[aValue] = 1;
1146         //aNbConnects = 1;
1147       }
1148       //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1149       aNodeId[1] = aNodeId[2];
1150       aNode1 = aNode2;
1151     }
1152     Value aValue(aNodeId[0],aNodeId[2]);
1153     MValues::iterator aItr = theValues.find(aValue);
1154     if (aItr != theValues.end()){
1155       aItr->second += 1;
1156       //aNbConnects = nb;
1157     } else {
1158       theValues[aValue] = 1;
1159       //aNbConnects = 1;
1160     }
1161     //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1162   }
1163
1164 }
1165
1166 /*
1167                             PREDICATES
1168 */
1169
1170 /*
1171   Class       : FreeBorders
1172   Description : Predicate for free borders
1173 */
1174
1175 FreeBorders::FreeBorders()
1176 {
1177   myMesh = 0;
1178 }
1179
1180 void FreeBorders::SetMesh( SMDS_Mesh* theMesh )
1181 {
1182   myMesh = theMesh;
1183 }
1184
1185 bool FreeBorders::IsSatisfy( long theId )
1186 {
1187   return getNbMultiConnection( myMesh, theId ) == 1;
1188 }
1189
1190 SMDSAbs_ElementType FreeBorders::GetType() const
1191 {
1192   return SMDSAbs_Edge;
1193 }
1194
1195
1196 /*
1197   Class       : FreeEdges
1198   Description : Predicate for free Edges
1199 */
1200 FreeEdges::FreeEdges()
1201 {
1202   myMesh = 0;
1203 }
1204
1205 void FreeEdges::SetMesh( SMDS_Mesh* theMesh )
1206 {
1207   myMesh = theMesh;
1208 }
1209
1210 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId  )
1211 {
1212   TColStd_MapOfInteger aMap;
1213   for ( int i = 0; i < 2; i++ )
1214   {
1215     SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator();
1216     while( anElemIter->more() )
1217     {
1218       const SMDS_MeshElement* anElem = anElemIter->next();
1219       if ( anElem != 0 && anElem->GetType() == SMDSAbs_Face )
1220       {
1221         int anId = anElem->GetID();
1222
1223         if ( i == 0 )
1224           aMap.Add( anId );
1225         else if ( aMap.Contains( anId ) && anId != theFaceId )
1226           return false;
1227       }
1228     }
1229   }
1230   return true;
1231 }
1232
1233 bool FreeEdges::IsSatisfy( long theId )
1234 {
1235   if ( myMesh == 0 )
1236     return false;
1237
1238   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
1239   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
1240     return false;
1241
1242   int nbNodes = aFace->NbNodes();
1243   const SMDS_MeshNode* aNodes[ nbNodes ];
1244   int i = 0;
1245   SMDS_ElemIteratorPtr anIter = aFace->nodesIterator();
1246   if ( anIter != 0 )
1247   {
1248     while( anIter->more() )
1249     {
1250       const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
1251       if ( aNode == 0 )
1252         return false;
1253       aNodes[ i++ ] = aNode;
1254     }
1255   }
1256
1257   for ( int i = 0; i < nbNodes - 1; i++ )
1258     if ( IsFreeEdge( &aNodes[ i ], theId ) )
1259       return true;
1260
1261   aNodes[ 1 ] = aNodes[ nbNodes - 1 ];
1262
1263   return IsFreeEdge( &aNodes[ 0 ], theId );
1264
1265 }
1266
1267 SMDSAbs_ElementType FreeEdges::GetType() const
1268 {
1269   return SMDSAbs_Face;
1270 }
1271
1272 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
1273   myElemId(theElemId)
1274 {
1275   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
1276   if(thePntId1 > thePntId2){
1277     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
1278   }
1279 }
1280
1281 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
1282   if(myPntId[0] < x.myPntId[0]) return true;
1283   if(myPntId[0] == x.myPntId[0])
1284     if(myPntId[1] < x.myPntId[1]) return true;
1285   return false;
1286 }
1287
1288 inline void UpdateBorders(const FreeEdges::Border& theBorder,
1289                           FreeEdges::TBorders& theRegistry,
1290                           FreeEdges::TBorders& theContainer)
1291 {
1292   if(theRegistry.find(theBorder) == theRegistry.end()){
1293     theRegistry.insert(theBorder);
1294     theContainer.insert(theBorder);
1295   }else{
1296     theContainer.erase(theBorder);
1297   }
1298 }
1299
1300 void FreeEdges::GetBoreders(TBorders& theBorders)
1301 {
1302   TBorders aRegistry;
1303   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1304   for(; anIter->more(); ){
1305     const SMDS_MeshFace* anElem = anIter->next();
1306     long anElemId = anElem->GetID();
1307     SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1308     long aNodeId[2];
1309     const SMDS_MeshElement* aNode;
1310     if(aNodesIter->more()){
1311       aNode = aNodesIter->next();
1312       aNodeId[0] = aNodeId[1] = aNode->GetID();
1313     }
1314     for(; aNodesIter->more(); ){
1315       aNode = aNodesIter->next();
1316       long anId = aNode->GetID();
1317       Border aBorder(anElemId,aNodeId[1],anId);
1318       aNodeId[1] = anId;
1319       //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
1320       UpdateBorders(aBorder,aRegistry,theBorders);
1321     }
1322     Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
1323     //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
1324     UpdateBorders(aBorder,aRegistry,theBorders);
1325   }
1326   //std::cout<<"theBorders.size() = "<<theBorders.size()<<endl;
1327 }
1328
1329 /*
1330   Class       : RangeOfIds
1331   Description : Predicate for Range of Ids.
1332                 Range may be specified with two ways.
1333                 1. Using AddToRange method
1334                 2. With SetRangeStr method. Parameter of this method is a string
1335                    like as "1,2,3,50-60,63,67,70-"
1336 */
1337
1338 //=======================================================================
1339 // name    : RangeOfIds
1340 // Purpose : Constructor
1341 //=======================================================================
1342 RangeOfIds::RangeOfIds()
1343 {
1344   myMesh = 0;
1345   myType = SMDSAbs_All;
1346 }
1347
1348 //=======================================================================
1349 // name    : SetMesh
1350 // Purpose : Set mesh
1351 //=======================================================================
1352 void RangeOfIds::SetMesh( SMDS_Mesh* theMesh )
1353 {
1354   myMesh = theMesh;
1355 }
1356
1357 //=======================================================================
1358 // name    : AddToRange
1359 // Purpose : Add ID to the range
1360 //=======================================================================
1361 bool RangeOfIds::AddToRange( long theEntityId )
1362 {
1363   myIds.Add( theEntityId );
1364   return true;
1365 }
1366
1367 //=======================================================================
1368 // name    : GetRangeStr
1369 // Purpose : Get range as a string.
1370 //           Example: "1,2,3,50-60,63,67,70-"
1371 //=======================================================================
1372 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
1373 {
1374   theResStr.Clear();
1375
1376   TColStd_SequenceOfInteger     anIntSeq;
1377   TColStd_SequenceOfAsciiString aStrSeq;
1378
1379   TColStd_MapIteratorOfMapOfInteger anIter( myIds );
1380   for ( ; anIter.More(); anIter.Next() )
1381   {
1382     int anId = anIter.Key();
1383     TCollection_AsciiString aStr( anId );
1384     anIntSeq.Append( anId );
1385     aStrSeq.Append( aStr );
1386   }
1387
1388   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
1389   {
1390     int aMinId = myMin( i );
1391     int aMaxId = myMax( i );
1392
1393     TCollection_AsciiString aStr;
1394     if ( aMinId != IntegerFirst() )
1395       aStr += aMinId;
1396
1397     aStr += "-";
1398
1399     if ( aMaxId != IntegerLast() )
1400       aStr += aMaxId;
1401
1402     // find position of the string in result sequence and insert string in it
1403     if ( anIntSeq.Length() == 0 )
1404     {
1405       anIntSeq.Append( aMinId );
1406       aStrSeq.Append( aStr );
1407     }
1408     else
1409     {
1410       if ( aMinId < anIntSeq.First() )
1411       {
1412         anIntSeq.Prepend( aMinId );
1413         aStrSeq.Prepend( aStr );
1414       }
1415       else if ( aMinId > anIntSeq.Last() )
1416       {
1417         anIntSeq.Append( aMinId );
1418         aStrSeq.Append( aStr );
1419       }
1420       else
1421         for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
1422           if ( aMinId < anIntSeq( j ) )
1423           {
1424             anIntSeq.InsertBefore( j, aMinId );
1425             aStrSeq.InsertBefore( j, aStr );
1426             break;
1427           }
1428     }
1429   }
1430
1431   if ( aStrSeq.Length() == 0 )
1432     return;
1433
1434   theResStr = aStrSeq( 1 );
1435   for ( int j = 2, k = aStrSeq.Length(); j <= k; j++  )
1436   {
1437     theResStr += ",";
1438     theResStr += aStrSeq( j );
1439   }
1440 }
1441
1442 //=======================================================================
1443 // name    : SetRangeStr
1444 // Purpose : Define range with string
1445 //           Example of entry string: "1,2,3,50-60,63,67,70-"
1446 //=======================================================================
1447 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
1448 {
1449   myMin.Clear();
1450   myMax.Clear();
1451   myIds.Clear();
1452
1453   TCollection_AsciiString aStr = theStr;
1454   aStr.RemoveAll( ' ' );
1455   aStr.RemoveAll( '\t' );
1456
1457   for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
1458     aStr.Remove( aPos, 2 );
1459
1460   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
1461   int i = 1;
1462   while ( tmpStr != "" )
1463   {
1464     tmpStr = aStr.Token( ",", i++ );
1465     int aPos = tmpStr.Search( '-' );
1466
1467     if ( aPos == -1 )
1468     {
1469       if ( tmpStr.IsIntegerValue() )
1470         myIds.Add( tmpStr.IntegerValue() );
1471       else
1472         return false;
1473     }
1474     else
1475     {
1476       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
1477       TCollection_AsciiString aMinStr = tmpStr;
1478
1479       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
1480       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
1481
1482       if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() ||
1483            !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() )
1484         return false;
1485
1486       myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
1487       myMax.Append( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
1488     }
1489   }
1490
1491   return true;
1492 }
1493
1494 //=======================================================================
1495 // name    : GetType
1496 // Purpose : Get type of supported entities
1497 //=======================================================================
1498 SMDSAbs_ElementType RangeOfIds::GetType() const
1499 {
1500   return myType;
1501 }
1502
1503 //=======================================================================
1504 // name    : SetType
1505 // Purpose : Set type of supported entities
1506 //=======================================================================
1507 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
1508 {
1509   myType = theType;
1510 }
1511
1512 //=======================================================================
1513 // name    : IsSatisfy
1514 // Purpose : Verify whether entity satisfies to this rpedicate
1515 //=======================================================================
1516 bool RangeOfIds::IsSatisfy( long theId )
1517 {
1518   if ( !myMesh )
1519     return false;
1520
1521   if ( myType == SMDSAbs_Node )
1522   {
1523     if ( myMesh->FindNode( theId ) == 0 )
1524       return false;
1525   }
1526   else
1527   {
1528     const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
1529     if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All )
1530       return false;
1531   }
1532
1533   if ( myIds.Contains( theId ) )
1534     return true;
1535
1536   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
1537     if ( theId >= myMin( i ) && theId <= myMax( i ) )
1538       return true;
1539
1540   return false;
1541 }
1542
1543 /*
1544   Class       : Comparator
1545   Description : Base class for comparators
1546 */
1547 Comparator::Comparator():
1548   myMargin(0)
1549 {}
1550
1551 Comparator::~Comparator()
1552 {}
1553
1554 void Comparator::SetMesh( SMDS_Mesh* theMesh )
1555 {
1556   if ( myFunctor )
1557     myFunctor->SetMesh( theMesh );
1558 }
1559
1560 void Comparator::SetMargin( double theValue )
1561 {
1562   myMargin = theValue;
1563 }
1564
1565 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
1566 {
1567   myFunctor = theFunct;
1568 }
1569
1570 SMDSAbs_ElementType Comparator::GetType() const
1571 {
1572   return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
1573 }
1574
1575 double Comparator::GetMargin()
1576 {
1577   return myMargin;
1578 }
1579
1580
1581 /*
1582   Class       : LessThan
1583   Description : Comparator "<"
1584 */
1585 bool LessThan::IsSatisfy( long theId )
1586 {
1587   return myFunctor && myFunctor->GetValue( theId ) < myMargin;
1588 }
1589
1590
1591 /*
1592   Class       : MoreThan
1593   Description : Comparator ">"
1594 */
1595 bool MoreThan::IsSatisfy( long theId )
1596 {
1597   return myFunctor && myFunctor->GetValue( theId ) > myMargin;
1598 }
1599
1600
1601 /*
1602   Class       : EqualTo
1603   Description : Comparator "="
1604 */
1605 EqualTo::EqualTo():
1606   myToler(Precision::Confusion())
1607 {}
1608
1609 bool EqualTo::IsSatisfy( long theId )
1610 {
1611   return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
1612 }
1613
1614 void EqualTo::SetTolerance( double theToler )
1615 {
1616   myToler = theToler;
1617 }
1618
1619 double EqualTo::GetTolerance()
1620 {
1621   return myToler;
1622 }
1623
1624 /*
1625   Class       : LogicalNOT
1626   Description : Logical NOT predicate
1627 */
1628 LogicalNOT::LogicalNOT()
1629 {}
1630
1631 LogicalNOT::~LogicalNOT()
1632 {}
1633
1634 bool LogicalNOT::IsSatisfy( long theId )
1635 {
1636   return myPredicate && !myPredicate->IsSatisfy( theId );
1637 }
1638
1639 void LogicalNOT::SetMesh( SMDS_Mesh* theMesh )
1640 {
1641   if ( myPredicate )
1642     myPredicate->SetMesh( theMesh );
1643 }
1644
1645 void LogicalNOT::SetPredicate( PredicatePtr thePred )
1646 {
1647   myPredicate = thePred;
1648 }
1649
1650 SMDSAbs_ElementType LogicalNOT::GetType() const
1651 {
1652   return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
1653 }
1654
1655
1656 /*
1657   Class       : LogicalBinary
1658   Description : Base class for binary logical predicate
1659 */
1660 LogicalBinary::LogicalBinary()
1661 {}
1662
1663 LogicalBinary::~LogicalBinary()
1664 {}
1665
1666 void LogicalBinary::SetMesh( SMDS_Mesh* theMesh )
1667 {
1668   if ( myPredicate1 )
1669     myPredicate1->SetMesh( theMesh );
1670
1671   if ( myPredicate2 )
1672     myPredicate2->SetMesh( theMesh );
1673 }
1674
1675 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
1676 {
1677   myPredicate1 = thePredicate;
1678 }
1679
1680 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
1681 {
1682   myPredicate2 = thePredicate;
1683 }
1684
1685 SMDSAbs_ElementType LogicalBinary::GetType() const
1686 {
1687   if ( !myPredicate1 || !myPredicate2 )
1688     return SMDSAbs_All;
1689
1690   SMDSAbs_ElementType aType1 = myPredicate1->GetType();
1691   SMDSAbs_ElementType aType2 = myPredicate2->GetType();
1692
1693   return aType1 == aType2 ? aType1 : SMDSAbs_All;
1694 }
1695
1696
1697 /*
1698   Class       : LogicalAND
1699   Description : Logical AND
1700 */
1701 bool LogicalAND::IsSatisfy( long theId )
1702 {
1703   return
1704     myPredicate1 &&
1705     myPredicate2 &&
1706     myPredicate1->IsSatisfy( theId ) &&
1707     myPredicate2->IsSatisfy( theId );
1708 }
1709
1710
1711 /*
1712   Class       : LogicalOR
1713   Description : Logical OR
1714 */
1715 bool LogicalOR::IsSatisfy( long theId )
1716 {
1717   return
1718     myPredicate1 &&
1719     myPredicate2 &&
1720     myPredicate1->IsSatisfy( theId ) ||
1721     myPredicate2->IsSatisfy( theId );
1722 }
1723
1724
1725 /*
1726                               FILTER
1727 */
1728
1729 Filter::Filter()
1730 {}
1731
1732 Filter::~Filter()
1733 {}
1734
1735 void Filter::SetPredicate( PredicatePtr thePredicate )
1736 {
1737   myPredicate = thePredicate;
1738 }
1739
1740
1741 template<class TElement, class TIterator, class TPredicate>
1742 void FillSequence(const TIterator& theIterator,
1743                   TPredicate& thePredicate,
1744                   Filter::TIdSequence& theSequence)
1745 {
1746   if ( theIterator ) {
1747     while( theIterator->more() ) {
1748       TElement anElem = theIterator->next();
1749       long anId = anElem->GetID();
1750       if ( thePredicate->IsSatisfy( anId ) )
1751         theSequence.push_back( anId );
1752     }
1753   }
1754 }
1755
1756 Filter::TIdSequence
1757 Filter::GetElementsId( SMDS_Mesh* theMesh )
1758 {
1759   TIdSequence aSequence;
1760   if ( !theMesh || !myPredicate ) return aSequence;
1761
1762   myPredicate->SetMesh( theMesh );
1763
1764   SMDSAbs_ElementType aType = myPredicate->GetType();
1765   switch(aType){
1766   case SMDSAbs_Node:{
1767     FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),myPredicate,aSequence);
1768     break;
1769   }
1770   case SMDSAbs_Edge:{
1771     FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1772     break;
1773   }
1774   case SMDSAbs_Face:{
1775     FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1776     break;
1777   }
1778   case SMDSAbs_Volume:{
1779     FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1780     break;
1781   }
1782   case SMDSAbs_All:{
1783     FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1784     FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1785     FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1786     break;
1787   }
1788   }
1789   return aSequence;
1790 }
1791
1792 /*
1793                               ManifoldPart
1794 */
1795
1796 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
1797
1798 /*
1799    Internal class Link
1800 */
1801
1802 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
1803                           SMDS_MeshNode* theNode2 )
1804 {
1805   myNode1 = theNode1;
1806   myNode2 = theNode2;
1807 }
1808
1809 ManifoldPart::Link::~Link()
1810 {
1811   myNode1 = 0;
1812   myNode2 = 0;
1813 }
1814
1815 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
1816 {
1817   if ( myNode1 == theLink.myNode1 &&
1818        myNode2 == theLink.myNode2 )
1819     return true;
1820   else if ( myNode1 == theLink.myNode2 &&
1821             myNode2 == theLink.myNode1 )
1822     return true;
1823   else
1824     return false;
1825 }
1826
1827 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
1828 {
1829   if(myNode1 < x.myNode1) return true;
1830   if(myNode1 == x.myNode1)
1831     if(myNode2 < x.myNode2) return true;
1832   return false;
1833 }
1834
1835 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
1836                             const ManifoldPart::Link& theLink2 )
1837 {
1838   return theLink1.IsEqual( theLink2 );
1839 }
1840
1841 ManifoldPart::ManifoldPart()
1842 {
1843   myMesh = 0;
1844   myAngToler = Precision::Angular();
1845   myIsOnlyManifold = true;
1846 }
1847
1848 ManifoldPart::~ManifoldPart()
1849 {
1850   myMesh = 0;
1851 }
1852
1853 void ManifoldPart::SetMesh( SMDS_Mesh* theMesh )
1854 {
1855   myMesh = theMesh;
1856   process();
1857 }
1858
1859 SMDSAbs_ElementType ManifoldPart::GetType() const
1860 { return SMDSAbs_Face; }
1861
1862 bool ManifoldPart::IsSatisfy( long theElementId )
1863 {
1864   return myMapIds.Contains( theElementId );
1865 }
1866
1867 void ManifoldPart::SetAngleTolerance( const double theAngToler )
1868 { myAngToler = theAngToler; }
1869
1870 double ManifoldPart::GetAngleTolerance() const
1871 { return myAngToler; }
1872
1873 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
1874 { myIsOnlyManifold = theIsOnly; }
1875
1876 void ManifoldPart::SetStartElem( const long  theStartId )
1877 { myStartElemId = theStartId; }
1878
1879 bool ManifoldPart::process()
1880 {
1881   myMapIds.Clear();
1882   myMapBadGeomIds.Clear();
1883
1884   myAllFacePtr.clear();
1885   myAllFacePtrIntDMap.clear();
1886   if ( !myMesh )
1887     return false;
1888
1889   // collect all faces into own map
1890   SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
1891   for (; anFaceItr->more(); )
1892   {
1893     SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
1894     myAllFacePtr.push_back( aFacePtr );
1895     myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
1896   }
1897
1898   SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
1899   if ( !aStartFace )
1900     return false;
1901
1902   // the map of non manifold links and bad geometry
1903   TMapOfLink aMapOfNonManifold;
1904   TColStd_MapOfInteger aMapOfTreated;
1905
1906   // begin cycle on faces from start index and run on vector till the end
1907   //  and from begin to start index to cover whole vector
1908   const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
1909   bool isStartTreat = false;
1910   for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
1911   {
1912     if ( fi == aStartIndx )
1913       isStartTreat = true;
1914     // as result next time when fi will be equal to aStartIndx
1915
1916     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
1917     if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
1918       continue;
1919
1920     aMapOfTreated.Add( aFacePtr->GetID() );
1921     TColStd_MapOfInteger aResFaces;
1922     if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
1923                          aMapOfNonManifold, aResFaces ) )
1924       continue;
1925     TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
1926     for ( ; anItr.More(); anItr.Next() )
1927     {
1928       int aFaceId = anItr.Key();
1929       aMapOfTreated.Add( aFaceId );
1930       myMapIds.Add( aFaceId );
1931     }
1932
1933     if ( fi == ( myAllFacePtr.size() - 1 ) )
1934       fi = 0;
1935   } // end run on vector of faces
1936   return !myMapIds.IsEmpty();
1937 }
1938
1939 static void getLinks( const SMDS_MeshFace* theFace,
1940                       ManifoldPart::TVectorOfLink& theLinks )
1941 {
1942   int aNbNode = theFace->NbNodes();
1943   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1944   int i = 1;
1945   SMDS_MeshNode* aNode = 0;
1946   for ( ; aNodeItr->more() && i <= aNbNode; )
1947   {
1948
1949     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
1950     if ( i == 1 )
1951       aNode = aN1;
1952     i++;
1953     SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
1954     i++;
1955     ManifoldPart::Link aLink( aN1, aN2 );
1956     theLinks.push_back( aLink );
1957   }
1958 }
1959
1960 static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
1961 {
1962   gp_XYZ n;
1963   int aNbNode = theFace->NbNodes();
1964   TColgp_Array1OfXYZ anArrOfXYZ(1,4);
1965   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1966   int i = 1;
1967   for ( ; aNodeItr->more() && i <= 4; i++ )
1968   {
1969     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
1970     anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
1971   }
1972
1973   gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1);
1974   gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1);
1975   n  = q1 ^ q2;
1976   if ( aNbNode > 3 )
1977   {
1978     gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
1979     n += q2 ^ q3;
1980   }
1981   double len = n.Modulus();
1982   if ( len > 0 )
1983     n /= len;
1984
1985   return n;
1986 }
1987
1988 bool ManifoldPart::findConnected
1989                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
1990                   SMDS_MeshFace*                           theStartFace,
1991                   ManifoldPart::TMapOfLink&                theNonManifold,
1992                   TColStd_MapOfInteger&                    theResFaces )
1993 {
1994   theResFaces.Clear();
1995   if ( !theAllFacePtrInt.size() )
1996     return false;
1997
1998   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
1999   {
2000     myMapBadGeomIds.Add( theStartFace->GetID() );
2001     return false;
2002   }
2003
2004   ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
2005   ManifoldPart::TVectorOfLink aSeqOfBoundary;
2006   theResFaces.Add( theStartFace->GetID() );
2007   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
2008
2009   expandBoundary( aMapOfBoundary, aSeqOfBoundary,
2010                  aDMapLinkFace, theNonManifold, theStartFace );
2011
2012   bool isDone = false;
2013   while ( !isDone && aMapOfBoundary.size() != 0 )
2014   {
2015     bool isToReset = false;
2016     ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
2017     for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
2018     {
2019       ManifoldPart::Link aLink = *pLink;
2020       if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
2021         continue;
2022       // each link could be treated only once
2023       aMapToSkip.insert( aLink );
2024
2025       ManifoldPart::TVectorOfFacePtr aFaces;
2026       // find next
2027       if ( myIsOnlyManifold &&
2028            (theNonManifold.find( aLink ) != theNonManifold.end()) )
2029         continue;
2030       else
2031       {
2032         getFacesByLink( aLink, aFaces );
2033         // filter the element to keep only indicated elements
2034         ManifoldPart::TVectorOfFacePtr aFiltered;
2035         ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
2036         for ( ; pFace != aFaces.end(); ++pFace )
2037         {
2038           SMDS_MeshFace* aFace = *pFace;
2039           if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
2040             aFiltered.push_back( aFace );
2041         }
2042         aFaces = aFiltered;
2043         if ( aFaces.size() < 2 )  // no neihgbour faces
2044           continue;
2045         else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
2046         {
2047           theNonManifold.insert( aLink );
2048           continue;
2049         }
2050       }
2051
2052       // compare normal with normals of neighbor element
2053       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
2054       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
2055       for ( ; pFace != aFaces.end(); ++pFace )
2056       {
2057         SMDS_MeshFace* aNextFace = *pFace;
2058         if ( aPrevFace == aNextFace )
2059           continue;
2060         int anNextFaceID = aNextFace->GetID();
2061         if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
2062          // should not be with non manifold restriction. probably bad topology
2063           continue;
2064         // check if face was treated and skipped
2065         if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
2066              !isInPlane( aPrevFace, aNextFace ) )
2067           continue;
2068         // add new element to connected and extend the boundaries.
2069         theResFaces.Add( anNextFaceID );
2070         expandBoundary( aMapOfBoundary, aSeqOfBoundary,
2071                         aDMapLinkFace, theNonManifold, aNextFace );
2072         isToReset = true;
2073       }
2074     }
2075     isDone = !isToReset;
2076   }
2077
2078   return !theResFaces.IsEmpty();
2079 }
2080
2081 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
2082                               const SMDS_MeshFace* theFace2 )
2083 {
2084   gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
2085   gp_XYZ aNorm2XYZ = getNormale( theFace2 );
2086   if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
2087   {
2088     myMapBadGeomIds.Add( theFace2->GetID() );
2089     return false;
2090   }
2091   if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
2092     return true;
2093
2094   return false;
2095 }
2096
2097 void ManifoldPart::expandBoundary
2098                    ( ManifoldPart::TMapOfLink&            theMapOfBoundary,
2099                      ManifoldPart::TVectorOfLink&         theSeqOfBoundary,
2100                      ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
2101                      ManifoldPart::TMapOfLink&            theNonManifold,
2102                      SMDS_MeshFace*                       theNextFace ) const
2103 {
2104   ManifoldPart::TVectorOfLink aLinks;
2105   getLinks( theNextFace, aLinks );
2106   int aNbLink = aLinks.size();
2107   for ( int i = 0; i < aNbLink; i++ )
2108   {
2109     ManifoldPart::Link aLink = aLinks[ i ];
2110     if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
2111       continue;
2112     if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
2113     {
2114       if ( myIsOnlyManifold )
2115       {
2116         // remove from boundary
2117         theMapOfBoundary.erase( aLink );
2118         ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
2119         for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
2120         {
2121           ManifoldPart::Link aBoundLink = *pLink;
2122           if ( aBoundLink.IsEqual( aLink ) )
2123           {
2124             theSeqOfBoundary.erase( pLink );
2125             break;
2126           }
2127         }
2128       }
2129     }
2130     else
2131     {
2132       theMapOfBoundary.insert( aLink );
2133       theSeqOfBoundary.push_back( aLink );
2134       theDMapLinkFacePtr[ aLink ] = theNextFace;
2135     }
2136   }
2137 }
2138
2139 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
2140                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
2141 {
2142   SMDS_Mesh::SetOfFaces aSetOfFaces;
2143   // take all faces that shared first node
2144   SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
2145   for ( ; anItr->more(); )
2146   {
2147     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
2148     if ( !aFace )
2149       continue;
2150     aSetOfFaces.insert( aFace );
2151   }
2152   // take all faces that shared second node
2153   anItr = theLink.myNode2->facesIterator();
2154   // find the common part of two sets
2155   for ( ; anItr->more(); )
2156   {
2157     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
2158     if ( aSetOfFaces.find( aFace ) != aSetOfFaces.end() )
2159       theFaces.push_back( aFace );
2160   }
2161 }
2162
2163
2164 /*
2165    ElementsOnSurface
2166 */
2167
2168 ElementsOnSurface::ElementsOnSurface()
2169 {
2170   myMesh = 0;
2171   myIds.Clear();
2172   myType = SMDSAbs_All;
2173   mySurf.Nullify();
2174   myToler = Precision::Confusion();
2175 }
2176
2177 ElementsOnSurface::~ElementsOnSurface()
2178 {
2179   myMesh = 0;
2180 }
2181
2182 void ElementsOnSurface::SetMesh( SMDS_Mesh* theMesh )
2183 {
2184   if ( myMesh == theMesh )
2185     return;
2186   myMesh = theMesh;
2187   myIds.Clear();
2188   process();
2189 }
2190
2191 bool ElementsOnSurface::IsSatisfy( long theElementId )
2192 {
2193   return myIds.Contains( theElementId );
2194 }
2195
2196 SMDSAbs_ElementType ElementsOnSurface::GetType() const
2197 { return myType; }
2198
2199 void ElementsOnSurface::SetTolerance( const double theToler )
2200 { myToler = theToler; }
2201
2202 double ElementsOnSurface::GetTolerance() const
2203 {
2204   return myToler;
2205 }
2206
2207 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
2208                                     const SMDSAbs_ElementType theType )
2209 {
2210   myType = theType;
2211   mySurf.Nullify();
2212   if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
2213   {
2214     mySurf.Nullify();
2215     return;
2216   }
2217   TopoDS_Face aFace = TopoDS::Face( theShape );
2218   mySurf = BRep_Tool::Surface( aFace );
2219 }
2220
2221 void ElementsOnSurface::process()
2222 {
2223   myIds.Clear();
2224   if ( mySurf.IsNull() )
2225     return;
2226
2227   if ( myMesh == 0 )
2228     return;
2229
2230   if ( myType == SMDSAbs_Face || myType == SMDSAbs_All )
2231   {
2232     SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
2233     for(; anIter->more(); )
2234       process( anIter->next() );
2235   }
2236
2237   if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All )
2238   {
2239     SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
2240     for(; anIter->more(); )
2241       process( anIter->next() );
2242   }
2243
2244   if ( myType == SMDSAbs_Node )
2245   {
2246     SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
2247     for(; anIter->more(); )
2248       process( anIter->next() );
2249   }
2250 }
2251
2252 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
2253 {
2254   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
2255   bool isSatisfy = true;
2256   for ( ; aNodeItr->more(); )
2257   {
2258     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
2259     if ( !isOnSurface( aNode ) )
2260     {
2261       isSatisfy = false;
2262       break;
2263     }
2264   }
2265   if ( isSatisfy )
2266     myIds.Add( theElemPtr->GetID() );
2267 }
2268
2269 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) const
2270 {
2271   if ( mySurf.IsNull() )
2272     return false;
2273
2274   gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
2275   double aToler2 = myToler * myToler;
2276   if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
2277   {
2278     gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
2279     if ( aPln.SquareDistance( aPnt ) > aToler2 )
2280       return false;
2281   }
2282   else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
2283   {
2284     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
2285     double aRad = aCyl.Radius();
2286     gp_Ax3 anAxis = aCyl.Position();
2287     gp_XYZ aLoc = aCyl.Location().XYZ();
2288     double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
2289     double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
2290     if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
2291       return false;
2292   }
2293   else
2294     return false;
2295
2296   return true;
2297 }