src/Controls/SMESH_Controls.cxx

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