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