src/Controls/SMESH_Controls.cxx

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