src/Controls/SMESH_Controls.cxx

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