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