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