src/Controls/SMESH_Controls.cxx

   1 // Copyright (C) 2007-2021  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, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "SMESH_ControlsDef.hxx"
  24
  25 #include "SMDS_BallElement.hxx"
  26 #include "SMDS_FacePosition.hxx"
  27 #include "SMDS_Iterator.hxx"
  28 #include "SMDS_Mesh.hxx"
  29 #include "SMDS_MeshElement.hxx"
  30 #include "SMDS_MeshNode.hxx"
  31 #include "SMDS_VolumeTool.hxx"
  32 #include "SMESHDS_GroupBase.hxx"
  33 #include "SMESHDS_GroupOnFilter.hxx"
  34 #include "SMESHDS_Mesh.hxx"
  35 #include "SMESH_MeshAlgos.hxx"
  36 #include "SMESH_OctreeNode.hxx"
  37 #include "SMESH_Comment.hxx"
  38
  39 #include <GEOMUtils.hxx>
  40 #include <Basics_Utils.hxx>
  41
  42 #include <BRepAdaptor_Surface.hxx>
  43 #include <BRepBndLib.hxx>
  44 #include <BRepBuilderAPI_Copy.hxx>
  45 #include <BRepClass3d_SolidClassifier.hxx>
  46 #include <BRepClass_FaceClassifier.hxx>
  47 #include <BRep_Tool.hxx>
  48 #include <GeomLib_IsPlanarSurface.hxx>
  49 #include <Geom_CylindricalSurface.hxx>
  50 #include <Geom_Plane.hxx>
  51 #include <Geom_Surface.hxx>
  52 #include <NCollection_Map.hxx>
  53 #include <Precision.hxx>
  54 #include <ShapeAnalysis_Surface.hxx>
  55 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
  56 #include <TColStd_MapOfInteger.hxx>
  57 #include <TColStd_SequenceOfAsciiString.hxx>
  58 #include <TColgp_Array1OfXYZ.hxx>
  59 #include <TopAbs.hxx>
  60 #include <TopExp.hxx>
  61 #include <TopoDS.hxx>
  62 #include <TopoDS_Edge.hxx>
  63 #include <TopoDS_Face.hxx>
  64 #include <TopoDS_Iterator.hxx>
  65 #include <TopoDS_Shape.hxx>
  66 #include <TopoDS_Vertex.hxx>
  67 #include <gp_Ax3.hxx>
  68 #include <gp_Cylinder.hxx>
  69 #include <gp_Dir.hxx>
  70 #include <gp_Pln.hxx>
  71 #include <gp_Pnt.hxx>
  72 #include <gp_Vec.hxx>
  73 #include <gp_XYZ.hxx>
  74
  75 #include <vtkMeshQuality.h>
  76
  77 #include <set>
  78 #include <limits>
  79
  80 /*
  81                             AUXILIARY METHODS
  82 */
  83
  84 namespace {
  85
  86   const double theEps = 1e-100;
  87   const double theInf = 1e+100;
  88
  89   inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
  90   {
  91     return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() );
  92   }
  93
  94   inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
  95   {
  96     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
  97
  98     return v1.Magnitude() < gp::Resolution() ||
  99       v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
 100   }
 101
 102   inline double getCos2( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
 103   {
 104     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
 105     double dot = v1 * v2, len1 = v1.SquareMagnitude(), len2 = v2.SquareMagnitude();
 106
 107     return ( dot < 0 || len1 < gp::Resolution() || len2 < gp::Resolution() ? -1 :
 108              dot * dot / len1 / len2 );
 109   }
 110
 111   inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
 112   {
 113     gp_Vec aVec1( P2 - P1 );
 114     gp_Vec aVec2( P3 - P1 );
 115     return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
 116   }
 117
 118   inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
 119   {
 120     return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
 121   }
 122
 123
 124
 125   inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
 126   {
 127     double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
 128     return aDist;
 129   }
 130
 131   int getNbMultiConnection( const SMDS_Mesh* theMesh, const smIdType theId )
 132   {
 133     if ( theMesh == 0 )
 134       return 0;
 135
 136     const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
 137     if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */)
 138       return 0;
 139
 140     // for each pair of nodes in anEdge (there are 2 pairs in a quadratic edge)
 141     // count elements containing both nodes of the pair.
 142     // Note that there may be such cases for a quadratic edge (a horizontal line):
 143     //
 144     //  Case 1          Case 2
 145     //  |     |      |        |      |
 146     //  |     |      |        |      |
 147     //  +-----+------+  +-----+------+
 148     //  |            |  |            |
 149     //  |            |  |            |
 150     // result should be 2 in both cases
 151     //
 152     int aResult0 = 0, aResult1 = 0;
 153     // last node, it is a medium one in a quadratic edge
 154     const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
 155     const SMDS_MeshNode*    aNode0 = anEdge->GetNode( 0 );
 156     const SMDS_MeshNode*    aNode1 = anEdge->GetNode( 1 );
 157     if ( aNode1 == aLastNode ) aNode1 = 0;
 158
 159     SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator();
 160     while( anElemIter->more() ) {
 161       const SMDS_MeshElement* anElem = anElemIter->next();
 162       if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
 163         SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
 164         while ( anIter->more() ) {
 165           if ( const SMDS_MeshElement* anElemNode = anIter->next() ) {
 166             if ( anElemNode == aNode0 ) {
 167               aResult0++;
 168               if ( !aNode1 ) break; // not a quadratic edge
 169             }
 170             else if ( anElemNode == aNode1 )
 171               aResult1++;
 172           }
 173         }
 174       }
 175     }
 176     int aResult = std::max ( aResult0, aResult1 );
 177
 178     return aResult;
 179   }
 180
 181   gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 )
 182   {
 183     int aNbNode = theFace->NbNodes();
 184
 185     gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0));
 186     gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0));
 187     gp_XYZ n  = q1 ^ q2;
 188     if ( aNbNode > 3 ) {
 189       gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0));
 190       n += q2 ^ q3;
 191     }
 192     double len = n.Modulus();
 193     bool zeroLen = ( len <= std::numeric_limits<double>::min());
 194     if ( !zeroLen )
 195       n /= len;
 196
 197     if (ok) *ok = !zeroLen;
 198
 199     return n;
 200   }
 201 }
 202
 203
 204
 205 using namespace SMESH::Controls;
 206
 207 /*
 208  *                               FUNCTORS
 209  */
 210
 211 //================================================================================
 212 /*
 213   Class       : NumericalFunctor
 214   Description : Base class for numerical functors
 215 */
 216 //================================================================================
 217
 218 NumericalFunctor::NumericalFunctor():
 219   myMesh(NULL)
 220 {
 221   myPrecision = -1;
 222 }
 223
 224 void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh )
 225 {
 226   myMesh = theMesh;
 227 }
 228
 229 bool NumericalFunctor::GetPoints(const smIdType       theId,
 230                                  TSequenceOfXYZ& theRes ) const
 231 {
 232   theRes.clear();
 233
 234   if ( myMesh == 0 )
 235     return false;
 236
 237   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
 238   if ( !IsApplicable( anElem ))
 239     return false;
 240
 241   return GetPoints( anElem, theRes );
 242 }
 243
 244 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
 245                                  TSequenceOfXYZ&         theRes )
 246 {
 247   theRes.clear();
 248
 249   if ( anElem == 0 )
 250     return false;
 251
 252   theRes.reserve( anElem->NbNodes() );
 253   theRes.setElement( anElem );
 254
 255   // Get nodes of the element
 256   SMDS_NodeIteratorPtr anIter= anElem->interlacedNodesIterator();
 257   if ( anIter ) {
 258     SMESH_NodeXYZ p;
 259     while( anIter->more() ) {
 260       if ( p.Set( anIter->next() ))
 261         theRes.push_back( p );
 262     }
 263   }
 264
 265   return true;
 266 }
 267
 268 long  NumericalFunctor::GetPrecision() const
 269 {
 270   return myPrecision;
 271 }
 272
 273 void  NumericalFunctor::SetPrecision( const long thePrecision )
 274 {
 275   myPrecision = thePrecision;
 276   myPrecisionValue = pow( 10., (double)( myPrecision ) );
 277 }
 278
 279 double NumericalFunctor::GetValue( long theId )
 280 {
 281   double aVal = 0;
 282
 283   myCurrElement = myMesh->FindElement( theId );
 284
 285   TSequenceOfXYZ P;
 286   if ( GetPoints( theId, P )) // elem type is checked here
 287     aVal = Round( GetValue( P ));
 288
 289   return aVal;
 290 }
 291
 292 double NumericalFunctor::Round( const double & aVal )
 293 {
 294   return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal;
 295 }
 296
 297 //================================================================================
 298 /*!
 299  * \brief Return true if a value can be computed for a given element.
 300  *        Some NumericalFunctor's are meaningful for elements of a certain
 301  *        geometry only.
 302  */
 303 //================================================================================
 304
 305 bool NumericalFunctor::IsApplicable( const SMDS_MeshElement* element ) const
 306 {
 307   return element && element->GetType() == this->GetType();
 308 }
 309
 310 bool NumericalFunctor::IsApplicable( long theElementId ) const
 311 {
 312   return IsApplicable( myMesh->FindElement( theElementId ));
 313 }
 314
 315 //================================================================================
 316 /*!
 317  * \brief Return histogram of functor values
 318  *  \param nbIntervals - number of intervals
 319  *  \param nbEvents - number of mesh elements having values within i-th interval
 320  *  \param funValues - boundaries of intervals
 321  *  \param elements - elements to check vulue of; empty list means "of all"
 322  *  \param minmax - boundaries of diapason of values to divide into intervals
 323  */
 324 //================================================================================
 325
 326 void NumericalFunctor::GetHistogram(int                          nbIntervals,
 327                                     std::vector<int>&            nbEvents,
 328                                     std::vector<double>&         funValues,
 329                                     const std::vector<smIdType>& elements,
 330                                     const double*                minmax,
 331                                     const bool                   isLogarithmic)
 332 {
 333   if ( nbIntervals < 1 ||
 334        !myMesh ||
 335        !myMesh->GetMeshInfo().NbElements( GetType() ))
 336     return;
 337   nbEvents.resize( nbIntervals, 0 );
 338   funValues.resize( nbIntervals+1 );
 339
 340   // get all values sorted
 341   std::multiset< double > values;
 342   if ( elements.empty() )
 343   {
 344     SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator( GetType() );
 345     while ( elemIt->more() )
 346       values.insert( GetValue( elemIt->next()->GetID() ));
 347   }
 348   else
 349   {
 350     std::vector<smIdType>::const_iterator id = elements.begin();
 351     for ( ; id != elements.end(); ++id )
 352       values.insert( GetValue( *id ));
 353   }
 354
 355   if ( minmax )
 356   {
 357     funValues[0] = minmax[0];
 358     funValues[nbIntervals] = minmax[1];
 359   }
 360   else
 361   {
 362     funValues[0] = *values.begin();
 363     funValues[nbIntervals] = *values.rbegin();
 364   }
 365   // case nbIntervals == 1
 366   if ( nbIntervals == 1 )
 367   {
 368     nbEvents[0] = values.size();
 369     return;
 370   }
 371   // case of 1 value
 372   if (funValues.front() == funValues.back())
 373   {
 374     nbEvents.resize( 1 );
 375     nbEvents[0] = values.size();
 376     funValues[1] = funValues.back();
 377     funValues.resize( 2 );
 378   }
 379   // generic case
 380   std::multiset< double >::iterator min = values.begin(), max;
 381   for ( int i = 0; i < nbIntervals; ++i )
 382   {
 383     // find end value of i-th interval
 384     double r = (i+1) / double(nbIntervals);
 385     if (isLogarithmic && funValues.front() > 1e-07 && funValues.back() > 1e-07) {
 386       double logmin = log10(funValues.front());
 387       double lval = logmin + r * (log10(funValues.back()) - logmin);
 388       funValues[i+1] = pow(10.0, lval);
 389     }
 390     else {
 391       funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
 392     }
 393
 394     // count values in the i-th interval if there are any
 395     if ( min != values.end() && *min <= funValues[i+1] )
 396     {
 397       // find the first value out of the interval
 398       max = values.upper_bound( funValues[i+1] ); // max is greater than funValues[i+1], or end()
 399       nbEvents[i] = std::distance( min, max );
 400       min = max;
 401     }
 402   }
 403   // add values larger than minmax[1]
 404   nbEvents.back() += std::distance( min, values.end() );
 405 }
 406
 407 //=======================================================================
 408 /*
 409   Class       : Volume
 410   Description : Functor calculating volume of a 3D element
 411 */
 412 //================================================================================
 413
 414 double Volume::GetValue( long theElementId )
 415 {
 416   if ( theElementId && myMesh ) {
 417     SMDS_VolumeTool aVolumeTool;
 418     if ( aVolumeTool.Set( myMesh->FindElement( theElementId )))
 419       return aVolumeTool.GetSize();
 420   }
 421   return 0;
 422 }
 423
 424 double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
 425 {
 426   return Value;
 427 }
 428
 429 SMDSAbs_ElementType Volume::GetType() const
 430 {
 431   return SMDSAbs_Volume;
 432 }
 433
 434 //=======================================================================
 435 /*
 436   Class       : MaxElementLength2D
 437   Description : Functor calculating maximum length of 2D element
 438 */
 439 //================================================================================
 440
 441 double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P )
 442 {
 443   if(P.size() == 0)
 444     return 0.;
 445   double aVal = 0;
 446   int len = P.size();
 447   if( len == 3 ) { // triangles
 448     double L1 = getDistance(P( 1 ),P( 2 ));
 449     double L2 = getDistance(P( 2 ),P( 3 ));
 450     double L3 = getDistance(P( 3 ),P( 1 ));
 451     aVal = Max(L1,Max(L2,L3));
 452   }
 453   else if( len == 4 ) { // quadrangles
 454     double L1 = getDistance(P( 1 ),P( 2 ));
 455     double L2 = getDistance(P( 2 ),P( 3 ));
 456     double L3 = getDistance(P( 3 ),P( 4 ));
 457     double L4 = getDistance(P( 4 ),P( 1 ));
 458     double D1 = getDistance(P( 1 ),P( 3 ));
 459     double D2 = getDistance(P( 2 ),P( 4 ));
 460     aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
 461   }
 462   else if( len == 6 ) { // quadratic triangles
 463     double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
 464     double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
 465     double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
 466     aVal = Max(L1,Max(L2,L3));
 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   }
 477   // Diagonals are undefined for concave polygons
 478   // else if ( P.getElementEntity() == SMDSEntity_Quad_Polygon && P.size() > 2 ) // quad polygon
 479   // {
 480   //   // sides
 481   //   aVal = getDistance( P( 1 ), P( P.size() )) + getDistance( P( P.size() ), P( P.size()-1 ));
 482   //   for ( size_t i = 1; i < P.size()-1; i += 2 )
 483   //   {
 484   //     double L = getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 ));
 485   //     aVal = Max( aVal, L );
 486   //   }
 487   //   // diagonals
 488   //   for ( int i = P.size()-5; i > 0; i -= 2 )
 489   //     for ( int j = i + 4; j < P.size() + i - 2; i += 2 )
 490   //     {
 491   //       double D = getDistance( P( i ), P( j ));
 492   //       aVal = Max( aVal, D );
 493   //     }
 494   // }
 495   // { // polygons
 496
 497   // }
 498
 499   if( myPrecision >= 0 )
 500   {
 501     double prec = pow( 10., (double)myPrecision );
 502     aVal = floor( aVal * prec + 0.5 ) / prec;
 503   }
 504   return aVal;
 505 }
 506
 507 double MaxElementLength2D::GetValue( long theElementId )
 508 {
 509   TSequenceOfXYZ P;
 510   return GetPoints( theElementId, P ) ? GetValue(P) : 0.0;
 511 }
 512
 513 double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const
 514 {
 515   return Value;
 516 }
 517
 518 SMDSAbs_ElementType MaxElementLength2D::GetType() const
 519 {
 520   return SMDSAbs_Face;
 521 }
 522
 523 //=======================================================================
 524 /*
 525   Class       : MaxElementLength3D
 526   Description : Functor calculating maximum length of 3D element
 527 */
 528 //================================================================================
 529
 530 double MaxElementLength3D::GetValue( long theElementId )
 531 {
 532   TSequenceOfXYZ P;
 533   if( GetPoints( theElementId, P ) ) {
 534     double aVal = 0;
 535     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
 536     SMDSAbs_EntityType      aType = aElem->GetEntityType();
 537     int len = P.size();
 538     switch ( aType ) {
 539     case SMDSEntity_Tetra: { // tetras
 540       double L1 = getDistance(P( 1 ),P( 2 ));
 541       double L2 = getDistance(P( 2 ),P( 3 ));
 542       double L3 = getDistance(P( 3 ),P( 1 ));
 543       double L4 = getDistance(P( 1 ),P( 4 ));
 544       double L5 = getDistance(P( 2 ),P( 4 ));
 545       double L6 = getDistance(P( 3 ),P( 4 ));
 546       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 547       break;
 548     }
 549     case SMDSEntity_Pyramid: { // pyramids
 550       double L1 = getDistance(P( 1 ),P( 2 ));
 551       double L2 = getDistance(P( 2 ),P( 3 ));
 552       double L3 = getDistance(P( 3 ),P( 4 ));
 553       double L4 = getDistance(P( 4 ),P( 1 ));
 554       double L5 = getDistance(P( 1 ),P( 5 ));
 555       double L6 = getDistance(P( 2 ),P( 5 ));
 556       double L7 = getDistance(P( 3 ),P( 5 ));
 557       double L8 = getDistance(P( 4 ),P( 5 ));
 558       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 559       aVal = Max(aVal,Max(L7,L8));
 560       break;
 561     }
 562     case SMDSEntity_Penta: { // pentas
 563       double L1 = getDistance(P( 1 ),P( 2 ));
 564       double L2 = getDistance(P( 2 ),P( 3 ));
 565       double L3 = getDistance(P( 3 ),P( 1 ));
 566       double L4 = getDistance(P( 4 ),P( 5 ));
 567       double L5 = getDistance(P( 5 ),P( 6 ));
 568       double L6 = getDistance(P( 6 ),P( 4 ));
 569       double L7 = getDistance(P( 1 ),P( 4 ));
 570       double L8 = getDistance(P( 2 ),P( 5 ));
 571       double L9 = getDistance(P( 3 ),P( 6 ));
 572       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 573       aVal = Max(aVal,Max(Max(L7,L8),L9));
 574       break;
 575     }
 576     case SMDSEntity_Hexa: { // hexas
 577       double L1 = getDistance(P( 1 ),P( 2 ));
 578       double L2 = getDistance(P( 2 ),P( 3 ));
 579       double L3 = getDistance(P( 3 ),P( 4 ));
 580       double L4 = getDistance(P( 4 ),P( 1 ));
 581       double L5 = getDistance(P( 5 ),P( 6 ));
 582       double L6 = getDistance(P( 6 ),P( 7 ));
 583       double L7 = getDistance(P( 7 ),P( 8 ));
 584       double L8 = getDistance(P( 8 ),P( 5 ));
 585       double L9 = getDistance(P( 1 ),P( 5 ));
 586       double L10= getDistance(P( 2 ),P( 6 ));
 587       double L11= getDistance(P( 3 ),P( 7 ));
 588       double L12= getDistance(P( 4 ),P( 8 ));
 589       double D1 = getDistance(P( 1 ),P( 7 ));
 590       double D2 = getDistance(P( 2 ),P( 8 ));
 591       double D3 = getDistance(P( 3 ),P( 5 ));
 592       double D4 = getDistance(P( 4 ),P( 6 ));
 593       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 594       aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
 595       aVal = Max(aVal,Max(L11,L12));
 596       aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
 597       break;
 598     }
 599     case SMDSEntity_Hexagonal_Prism: { // hexagonal prism
 600       for ( int i1 = 1; i1 < 12; ++i1 )
 601         for ( int i2 = i1+1; i1 <= 12; ++i1 )
 602           aVal = Max( aVal, getDistance(P( i1 ),P( i2 )));
 603       break;
 604     }
 605     case SMDSEntity_Quad_Tetra: { // quadratic tetras
 606       double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
 607       double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
 608       double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
 609       double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
 610       double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
 611       double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
 612       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 613       break;
 614     }
 615     case SMDSEntity_Quad_Pyramid: { // quadratic pyramids
 616       double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
 617       double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
 618       double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
 619       double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
 620       double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
 621       double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
 622       double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
 623       double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
 624       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 625       aVal = Max(aVal,Max(L7,L8));
 626       break;
 627     }
 628     case SMDSEntity_Quad_Penta:
 629     case SMDSEntity_BiQuad_Penta: { // quadratic pentas
 630       double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
 631       double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
 632       double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
 633       double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
 634       double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
 635       double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
 636       double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
 637       double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
 638       double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
 639       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 640       aVal = Max(aVal,Max(Max(L7,L8),L9));
 641       break;
 642     }
 643     case SMDSEntity_Quad_Hexa:
 644     case SMDSEntity_TriQuad_Hexa: { // quadratic hexas
 645       double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
 646       double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
 647       double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
 648       double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
 649       double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
 650       double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
 651       double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
 652       double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
 653       double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
 654       double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
 655       double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
 656       double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
 657       double D1 = getDistance(P( 1 ),P( 7 ));
 658       double D2 = getDistance(P( 2 ),P( 8 ));
 659       double D3 = getDistance(P( 3 ),P( 5 ));
 660       double D4 = getDistance(P( 4 ),P( 6 ));
 661       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 662       aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
 663       aVal = Max(aVal,Max(L11,L12));
 664       aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
 665       break;
 666     }
 667     case SMDSEntity_Quad_Polyhedra:
 668     case SMDSEntity_Polyhedra: { // polys
 669       // get the maximum distance between all pairs of nodes
 670       for( int i = 1; i <= len; i++ ) {
 671         for( int j = 1; j <= len; j++ ) {
 672           if( j > i ) { // optimization of the loop
 673             double D = getDistance( P(i), P(j) );
 674             aVal = Max( aVal, D );
 675           }
 676         }
 677       }
 678       break;
 679     }
 680     case SMDSEntity_Node:
 681     case SMDSEntity_0D:
 682     case SMDSEntity_Edge:
 683     case SMDSEntity_Quad_Edge:
 684     case SMDSEntity_Triangle:
 685     case SMDSEntity_Quad_Triangle:
 686     case SMDSEntity_BiQuad_Triangle:
 687     case SMDSEntity_Quadrangle:
 688     case SMDSEntity_Quad_Quadrangle:
 689     case SMDSEntity_BiQuad_Quadrangle:
 690     case SMDSEntity_Polygon:
 691     case SMDSEntity_Quad_Polygon:
 692     case SMDSEntity_Ball:
 693     case SMDSEntity_Last: return 0;
 694     } // switch ( aType )
 695
 696     if( myPrecision >= 0 )
 697     {
 698       double prec = pow( 10., (double)myPrecision );
 699       aVal = floor( aVal * prec + 0.5 ) / prec;
 700     }
 701     return aVal;
 702   }
 703   return 0.;
 704 }
 705
 706 double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const
 707 {
 708   return Value;
 709 }
 710
 711 SMDSAbs_ElementType MaxElementLength3D::GetType() const
 712 {
 713   return SMDSAbs_Volume;
 714 }
 715
 716 //=======================================================================
 717 /*
 718   Class       : MinimumAngle
 719   Description : Functor for calculation of minimum angle
 720 */
 721 //================================================================================
 722
 723 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
 724 {
 725   if ( P.size() < 3 )
 726     return 0.;
 727
 728   double aMaxCos2;
 729
 730   aMaxCos2 = getCos2( P( P.size() ), P( 1 ), P( 2 ));
 731   aMaxCos2 = Max( aMaxCos2, getCos2( P( P.size()-1 ), P( P.size() ), P( 1 )));
 732
 733   for ( size_t i = 2; i < P.size(); i++ )
 734   {
 735     double A0 = getCos2( P( i-1 ), P( i ), P( i+1 ) );
 736     aMaxCos2 = Max( aMaxCos2, A0 );
 737   }
 738   if ( aMaxCos2 < 0 )
 739     return 0; // all nodes coincide
 740
 741   double cos = sqrt( aMaxCos2 );
 742   if ( cos >=  1 ) return 0;
 743   return acos( cos ) * 180.0 / M_PI;
 744 }
 745
 746 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
 747 {
 748   //const double aBestAngle = PI / nbNodes;
 749   const double aBestAngle = 180.0 - ( 360.0 / double(nbNodes) );
 750   return ( fabs( aBestAngle - Value ));
 751 }
 752
 753 SMDSAbs_ElementType MinimumAngle::GetType() const
 754 {
 755   return SMDSAbs_Face;
 756 }
 757
 758
 759 //================================================================================
 760 /*
 761   Class       : AspectRatio
 762   Description : Functor for calculating aspect ratio
 763 */
 764 //================================================================================
 765
 766 double AspectRatio::GetValue( long theId )
 767 {
 768   double aVal = 0;
 769   myCurrElement = myMesh->FindElement( theId );
 770   TSequenceOfXYZ P;
 771   if ( GetPoints( myCurrElement, P ))
 772     aVal = Round( GetValue( P ));
 773   return aVal;
 774 }
 775
 776 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
 777 {
 778   // According to "Mesh quality control" by Nadir Bouhamau referring to
 779   // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis.
 780   // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4
 781   // PAL10872
 782
 783   int nbNodes = P.size();
 784
 785   if ( nbNodes < 3 )
 786     return 0;
 787
 788   // Compute aspect ratio
 789
 790   if ( nbNodes == 3 ) {
 791     // Compute lengths of the sides
 792     double aLen1 = getDistance( P( 1 ), P( 2 ));
 793     double aLen2 = getDistance( P( 2 ), P( 3 ));
 794     double aLen3 = getDistance( P( 3 ), P( 1 ));
 795     // Q = alfa * h * p / S, where
 796     //
 797     // alfa = sqrt( 3 ) / 6
 798     // h - length of the longest edge
 799     // p - half perimeter
 800     // S - triangle surface
 801     const double     alfa = sqrt( 3. ) / 6.;
 802     double         maxLen = Max( aLen1, Max( aLen2, aLen3 ));
 803     double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.;
 804     double         anArea = getArea( P( 1 ), P( 2 ), P( 3 ));
 805     if ( anArea <= theEps  )
 806       return theInf;
 807     return alfa * maxLen * half_perimeter / anArea;
 808   }
 809   else if ( nbNodes == 6 ) { // quadratic triangles
 810     // Compute lengths of the sides
 811     double aLen1 = getDistance( P( 1 ), P( 3 ));
 812     double aLen2 = getDistance( P( 3 ), P( 5 ));
 813     double aLen3 = getDistance( P( 5 ), P( 1 ));
 814     // algo same as for the linear triangle
 815     const double     alfa = sqrt( 3. ) / 6.;
 816     double         maxLen = Max( aLen1, Max( aLen2, aLen3 ));
 817     double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.;
 818     double         anArea = getArea( P( 1 ), P( 3 ), P( 5 ));
 819     if ( anArea <= theEps )
 820       return theInf;
 821     return alfa * maxLen * half_perimeter / anArea;
 822   }
 823   else if( nbNodes == 4 ) { // quadrangle
 824     // Compute lengths of the sides
 825     double aLen[4];
 826     aLen[0] = getDistance( P(1), P(2) );
 827     aLen[1] = getDistance( P(2), P(3) );
 828     aLen[2] = getDistance( P(3), P(4) );
 829     aLen[3] = getDistance( P(4), P(1) );
 830     // Compute lengths of the diagonals
 831     double aDia[2];
 832     aDia[0] = getDistance( P(1), P(3) );
 833     aDia[1] = getDistance( P(2), P(4) );
 834     // Compute areas of all triangles which can be built
 835     // taking three nodes of the quadrangle
 836     double anArea[4];
 837     anArea[0] = getArea( P(1), P(2), P(3) );
 838     anArea[1] = getArea( P(1), P(2), P(4) );
 839     anArea[2] = getArea( P(1), P(3), P(4) );
 840     anArea[3] = getArea( P(2), P(3), P(4) );
 841     // Q = alpha * L * C1 / C2, where
 842     //
 843     // alpha = sqrt( 1/32 )
 844     // L = max( L1, L2, L3, L4, D1, D2 )
 845     // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
 846     // C2 = min( S1, S2, S3, S4 )
 847     // Li - lengths of the edges
 848     // Di - lengths of the diagonals
 849     // Si - areas of the triangles
 850     const double alpha = sqrt( 1 / 32. );
 851     double L = Max( aLen[ 0 ],
 852                     Max( aLen[ 1 ],
 853                          Max( aLen[ 2 ],
 854                               Max( aLen[ 3 ],
 855                                    Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
 856     double C1 = sqrt( aLen[0] * aLen[0] +
 857                       aLen[1] * aLen[1] +
 858                       aLen[2] * aLen[2] +
 859                       aLen[3] * aLen[3] );
 860     double C2 = Min( anArea[ 0 ],
 861                      Min( anArea[ 1 ],
 862                           Min( anArea[ 2 ], anArea[ 3 ] ) ) );
 863     if ( C2 <= theEps )
 864       return theInf;
 865     return alpha * L * C1 / C2;
 866   }
 867   else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle
 868     // Compute lengths of the sides
 869     double aLen[4];
 870     aLen[0] = getDistance( P(1), P(3) );
 871     aLen[1] = getDistance( P(3), P(5) );
 872     aLen[2] = getDistance( P(5), P(7) );
 873     aLen[3] = getDistance( P(7), P(1) );
 874     // Compute lengths of the diagonals
 875     double aDia[2];
 876     aDia[0] = getDistance( P(1), P(5) );
 877     aDia[1] = getDistance( P(3), P(7) );
 878     // Compute areas of all triangles which can be built
 879     // taking three nodes of the quadrangle
 880     double anArea[4];
 881     anArea[0] = getArea( P(1), P(3), P(5) );
 882     anArea[1] = getArea( P(1), P(3), P(7) );
 883     anArea[2] = getArea( P(1), P(5), P(7) );
 884     anArea[3] = getArea( P(3), P(5), P(7) );
 885     // Q = alpha * L * C1 / C2, where
 886     //
 887     // alpha = sqrt( 1/32 )
 888     // L = max( L1, L2, L3, L4, D1, D2 )
 889     // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
 890     // C2 = min( S1, S2, S3, S4 )
 891     // Li - lengths of the edges
 892     // Di - lengths of the diagonals
 893     // Si - areas of the triangles
 894     const double alpha = sqrt( 1 / 32. );
 895     double L = Max( aLen[ 0 ],
 896                     Max( aLen[ 1 ],
 897                          Max( aLen[ 2 ],
 898                               Max( aLen[ 3 ],
 899                                    Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
 900     double C1 = sqrt( aLen[0] * aLen[0] +
 901                       aLen[1] * aLen[1] +
 902                       aLen[2] * aLen[2] +
 903                       aLen[3] * aLen[3] );
 904     double C2 = Min( anArea[ 0 ],
 905                      Min( anArea[ 1 ],
 906                           Min( anArea[ 2 ], anArea[ 3 ] ) ) );
 907     if ( C2 <= theEps )
 908       return theInf;
 909     return alpha * L * C1 / C2;
 910   }
 911   return 0;
 912 }
 913
 914 bool AspectRatio::IsApplicable( const SMDS_MeshElement* element ) const
 915 {
 916   return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() );
 917 }
 918
 919 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
 920 {
 921   // the aspect ratio is in the range [1.0,infinity]
 922   // < 1.0 = very bad, zero area
 923   // 1.0 = good
 924   // infinity = bad
 925   return ( Value < 0.9 ) ? 1000 : Value / 1000.;
 926 }
 927
 928 SMDSAbs_ElementType AspectRatio::GetType() const
 929 {
 930   return SMDSAbs_Face;
 931 }
 932
 933
 934 //================================================================================
 935 /*
 936   Class       : AspectRatio3D
 937   Description : Functor for calculating aspect ratio
 938 */
 939 //================================================================================
 940
 941 namespace{
 942
 943   inline double getHalfPerimeter(double theTria[3]){
 944     return (theTria[0] + theTria[1] + theTria[2])/2.0;
 945   }
 946
 947   inline double getArea(double theHalfPerim, double theTria[3]){
 948     return sqrt(theHalfPerim*
 949                 (theHalfPerim-theTria[0])*
 950                 (theHalfPerim-theTria[1])*
 951                 (theHalfPerim-theTria[2]));
 952   }
 953
 954   inline double getVolume(double theLen[6]){
 955     double a2 = theLen[0]*theLen[0];
 956     double b2 = theLen[1]*theLen[1];
 957     double c2 = theLen[2]*theLen[2];
 958     double d2 = theLen[3]*theLen[3];
 959     double e2 = theLen[4]*theLen[4];
 960     double f2 = theLen[5]*theLen[5];
 961     double P = 4.0*a2*b2*d2;
 962     double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
 963     double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
 964     return sqrt(P-Q+R)/12.0;
 965   }
 966
 967   inline double getVolume2(double theLen[6]){
 968     double a2 = theLen[0]*theLen[0];
 969     double b2 = theLen[1]*theLen[1];
 970     double c2 = theLen[2]*theLen[2];
 971     double d2 = theLen[3]*theLen[3];
 972     double e2 = theLen[4]*theLen[4];
 973     double f2 = theLen[5]*theLen[5];
 974
 975     double P = a2*e2*(b2+c2+d2+f2-a2-e2);
 976     double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
 977     double R = c2*d2*(a2+b2+e2+f2-c2-d2);
 978     double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
 979
 980     return sqrt(P+Q+R-S)/12.0;
 981   }
 982
 983   inline double getVolume(const TSequenceOfXYZ& P){
 984     gp_Vec aVec1( P( 2 ) - P( 1 ) );
 985     gp_Vec aVec2( P( 3 ) - P( 1 ) );
 986     gp_Vec aVec3( P( 4 ) - P( 1 ) );
 987     gp_Vec anAreaVec( aVec1 ^ aVec2 );
 988     return fabs(aVec3 * anAreaVec) / 6.0;
 989   }
 990
 991   inline double getMaxHeight(double theLen[6])
 992   {
 993     double aHeight = std::max(theLen[0],theLen[1]);
 994     aHeight = std::max(aHeight,theLen[2]);
 995     aHeight = std::max(aHeight,theLen[3]);
 996     aHeight = std::max(aHeight,theLen[4]);
 997     aHeight = std::max(aHeight,theLen[5]);
 998     return aHeight;
 999   }
1000
1001   //================================================================================
1002   /*!
1003    * \brief Standard quality of a tetrahedron but not normalized
1004    */
1005   //================================================================================
1006
1007   double tetQualityByHomardMethod( const gp_XYZ & p1,
1008                                    const gp_XYZ & p2,
1009                                    const gp_XYZ & p3,
1010                                    const gp_XYZ & p4 )
1011   {
1012     gp_XYZ edgeVec[6];
1013     edgeVec[0] = ( p1 - p2 );
1014     edgeVec[1] = ( p2 - p3 );
1015     edgeVec[2] = ( p3 - p1 );
1016     edgeVec[3] = ( p4 - p1 );
1017     edgeVec[4] = ( p4 - p2 );
1018     edgeVec[5] = ( p4 - p3 );
1019
1020     double maxEdgeLen2            = edgeVec[0].SquareModulus();
1021     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[1].SquareModulus() );
1022     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[2].SquareModulus() );
1023     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[3].SquareModulus() );
1024     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[4].SquareModulus() );
1025     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[5].SquareModulus() );
1026     double maxEdgeLen = Sqrt( maxEdgeLen2 );
1027
1028     gp_XYZ cross01 = edgeVec[0] ^ edgeVec[1];
1029     double sumArea = ( cross01                 ).Modulus(); // actually double area
1030     sumArea       += ( edgeVec[0] ^ edgeVec[3] ).Modulus();
1031     sumArea       += ( edgeVec[1] ^ edgeVec[4] ).Modulus();
1032     sumArea       += ( edgeVec[2] ^ edgeVec[5] ).Modulus();
1033
1034     double sixVolume = Abs( cross01 * edgeVec[4] ); // 6 * volume
1035     double quality   = maxEdgeLen * sumArea / sixVolume; // not normalized!!!
1036     return quality;
1037   }
1038
1039   //================================================================================
1040   /*!
1041    * \brief HOMARD method of hexahedron quality
1042    * 1. Decompose the hexa into 24 tetra: each face is splitted into 4 triangles by
1043    *    adding the diagonals and every triangle is connected to the center of the hexa.
1044    * 2. Compute the quality of every tetra with the same formula as for the standard quality,
1045    *    except that the factor for the normalization is not the same because the final goal
1046    *    is to have a quality equal to 1 for a perfect cube. So the formula is:
1047    *    qual = max(lengthes of 6 edges) * (sum of surfaces of 4 faces) / (7.6569*6*volume)
1048    * 3. The quality of the hexa is the highest value of the qualities of the 24 tetra
1049    */
1050   //================================================================================
1051
1052   double hexQualityByHomardMethod( const TSequenceOfXYZ& P )
1053   {
1054     gp_XYZ quadCenter[6];
1055     quadCenter[0] = ( P(1) + P(2) + P(3) + P(4) ) / 4.;
1056     quadCenter[1] = ( P(5) + P(6) + P(7) + P(8) ) / 4.;
1057     quadCenter[2] = ( P(1) + P(2) + P(6) + P(5) ) / 4.;
1058     quadCenter[3] = ( P(2) + P(3) + P(7) + P(6) ) / 4.;
1059     quadCenter[4] = ( P(3) + P(4) + P(8) + P(7) ) / 4.;
1060     quadCenter[5] = ( P(1) + P(4) + P(8) + P(5) ) / 4.;
1061
1062     gp_XYZ hexCenter = ( P(1) + P(2) + P(3) + P(4) + P(5) + P(6) + P(7) + P(8) ) / 8.;
1063
1064     // quad 1 ( 1 2 3 4 )
1065     double quality =        tetQualityByHomardMethod( P(1), P(2), quadCenter[0], hexCenter );
1066     quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[0], hexCenter ));
1067     quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[0], hexCenter ));
1068     quality = Max( quality, tetQualityByHomardMethod( P(4), P(1), quadCenter[0], hexCenter ));
1069     // quad 2 ( 5 6 7 8 )
1070     quality = Max( quality, tetQualityByHomardMethod( P(5), P(6), quadCenter[1], hexCenter ));
1071     quality = Max( quality, tetQualityByHomardMethod( P(6), P(7), quadCenter[1], hexCenter ));
1072     quality = Max( quality, tetQualityByHomardMethod( P(7), P(8), quadCenter[1], hexCenter ));
1073     quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[1], hexCenter ));
1074     // quad 3 ( 1 2 6 5 )
1075     quality = Max( quality, tetQualityByHomardMethod( P(1), P(2), quadCenter[2], hexCenter ));
1076     quality = Max( quality, tetQualityByHomardMethod( P(2), P(6), quadCenter[2], hexCenter ));
1077     quality = Max( quality, tetQualityByHomardMethod( P(6), P(5), quadCenter[2], hexCenter ));
1078     quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[2], hexCenter ));
1079     // quad 4 ( 2 3 7 6 )
1080     quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[3], hexCenter ));
1081     quality = Max( quality, tetQualityByHomardMethod( P(3), P(7), quadCenter[3], hexCenter ));
1082     quality = Max( quality, tetQualityByHomardMethod( P(7), P(6), quadCenter[3], hexCenter ));
1083     quality = Max( quality, tetQualityByHomardMethod( P(6), P(2), quadCenter[3], hexCenter ));
1084     // quad 5 ( 3 4 8 7 )
1085     quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[4], hexCenter ));
1086     quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[4], hexCenter ));
1087     quality = Max( quality, tetQualityByHomardMethod( P(8), P(7), quadCenter[4], hexCenter ));
1088     quality = Max( quality, tetQualityByHomardMethod( P(7), P(3), quadCenter[4], hexCenter ));
1089     // quad 6 ( 1 4 8 5 )
1090     quality = Max( quality, tetQualityByHomardMethod( P(1), P(4), quadCenter[5], hexCenter ));
1091     quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[5], hexCenter ));
1092     quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[5], hexCenter ));
1093     quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[5], hexCenter ));
1094
1095     return quality / 7.65685424949;
1096   }
1097 }
1098
1099 double AspectRatio3D::GetValue( long theId )
1100 {
1101   double aVal = 0;
1102   myCurrElement = myMesh->FindElement( theId );
1103   if ( myCurrElement && myCurrElement->GetVtkType() == VTK_TETRA )
1104   {
1105     // Action from CoTech | ACTION 31.3:
1106     // EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with
1107     // those of ParaView. The library used by ParaView for those calculations can be reused in SMESH.
1108     vtkUnstructuredGrid* grid = const_cast<SMDS_Mesh*>( myMesh )->GetGrid();
1109     if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->GetVtkID() ))
1110       aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
1111   }
1112   else
1113   {
1114     TSequenceOfXYZ P;
1115     if ( GetPoints( myCurrElement, P ))
1116       aVal = Round( GetValue( P ));
1117   }
1118   return aVal;
1119 }
1120
1121 bool AspectRatio3D::IsApplicable( const SMDS_MeshElement* element ) const
1122 {
1123   return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() );
1124 }
1125
1126 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
1127 {
1128   double aQuality = 0.0;
1129   if(myCurrElement->IsPoly()) return aQuality;
1130
1131   int nbNodes = P.size();
1132
1133   if( myCurrElement->IsQuadratic() ) {
1134     if     (nbNodes==10) nbNodes=4; // quadratic tetrahedron
1135     else if(nbNodes==13) nbNodes=5; // quadratic pyramid
1136     else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
1137     else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
1138     else if(nbNodes==27) nbNodes=8; // tri-quadratic hexahedron
1139     else return aQuality;
1140   }
1141
1142   switch(nbNodes) {
1143   case 4:{
1144     double aLen[6] = {
1145       getDistance(P( 1 ),P( 2 )), // a
1146       getDistance(P( 2 ),P( 3 )), // b
1147       getDistance(P( 3 ),P( 1 )), // c
1148       getDistance(P( 2 ),P( 4 )), // d
1149       getDistance(P( 3 ),P( 4 )), // e
1150       getDistance(P( 1 ),P( 4 ))  // f
1151     };
1152     double aTria[4][3] = {
1153       {aLen[0],aLen[1],aLen[2]}, // abc
1154       {aLen[0],aLen[3],aLen[5]}, // adf
1155       {aLen[1],aLen[3],aLen[4]}, // bde
1156       {aLen[2],aLen[4],aLen[5]}  // cef
1157     };
1158     double aSumArea = 0.0;
1159     double aHalfPerimeter = getHalfPerimeter(aTria[0]);
1160     double anArea = getArea(aHalfPerimeter,aTria[0]);
1161     aSumArea += anArea;
1162     aHalfPerimeter = getHalfPerimeter(aTria[1]);
1163     anArea = getArea(aHalfPerimeter,aTria[1]);
1164     aSumArea += anArea;
1165     aHalfPerimeter = getHalfPerimeter(aTria[2]);
1166     anArea = getArea(aHalfPerimeter,aTria[2]);
1167     aSumArea += anArea;
1168     aHalfPerimeter = getHalfPerimeter(aTria[3]);
1169     anArea = getArea(aHalfPerimeter,aTria[3]);
1170     aSumArea += anArea;
1171     double aVolume = getVolume(P);
1172     //double aVolume = getVolume(aLen);
1173     double aHeight = getMaxHeight(aLen);
1174     static double aCoeff = sqrt(2.0)/12.0;
1175     if ( aVolume > DBL_MIN )
1176       aQuality = aCoeff*aHeight*aSumArea/aVolume;
1177     break;
1178   }
1179   case 5:{
1180     {
1181       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
1182       aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
1183     }
1184     {
1185       gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
1186       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1187     }
1188     {
1189       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
1190       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1191     }
1192     {
1193       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
1194       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1195     }
1196     break;
1197   }
1198   case 6:{
1199     {
1200       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
1201       aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
1202     }
1203     {
1204       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
1205       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1206     }
1207     {
1208       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
1209       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1210     }
1211     {
1212       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1213       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1214     }
1215     {
1216       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
1217       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1218     }
1219     {
1220       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
1221       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1222     }
1223     break;
1224   }
1225   case 8:{
1226
1227     return hexQualityByHomardMethod( P ); // bos #23982
1228
1229
1230     {
1231       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1232       aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
1233     }
1234     {
1235       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
1236       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1237     }
1238     {
1239       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
1240       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1241     }
1242     {
1243       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
1244       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1245     }
1246     {
1247       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
1248       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1249     }
1250     {
1251       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
1252       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1253     }
1254     {
1255       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
1256       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1257     }
1258     {
1259       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
1260       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1261     }
1262     {
1263       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
1264       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1265     }
1266     {
1267       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
1268       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1269     }
1270     {
1271       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
1272       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1273     }
1274     {
1275       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
1276       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1277     }
1278     {
1279       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
1280       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1281     }
1282     {
1283       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
1284       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1285     }
1286     {
1287       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
1288       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1289     }
1290     {
1291       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
1292       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1293     }
1294     {
1295       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
1296       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1297     }
1298     {
1299       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
1300       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1301     }
1302     {
1303       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
1304       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1305     }
1306     {
1307       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
1308       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1309     }
1310     {
1311       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
1312       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1313     }
1314     {
1315       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1316       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1317     }
1318     {
1319       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
1320       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1321     }
1322     {
1323       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
1324       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1325     }
1326     {
1327       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1328       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1329     }
1330     {
1331       gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
1332       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1333     }
1334     {
1335       gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
1336       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1337     }
1338     {
1339       gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
1340       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1341     }
1342     {
1343       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
1344       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1345     }
1346     {
1347       gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
1348       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1349     }
1350     {
1351       gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
1352       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1353     }
1354     {
1355       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
1356       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1357     }
1358     {
1359       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
1360       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1361     }
1362     break;
1363   }
1364   case 12:
1365     {
1366       gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
1367       aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8]));
1368     }
1369     {
1370       gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
1371       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1372     }
1373     {
1374       gp_XYZ aXYZ[8] = {P( 3 ),P( 4 ),P( 6 ),P( 1 ),P( 9 ),P( 10 ),P( 12 ),P( 7 )};
1375       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1376     }
1377     break;
1378   } // switch(nbNodes)
1379
1380   if ( nbNodes > 4 ) {
1381     // evaluate aspect ratio of quadrangle faces
1382     AspectRatio aspect2D;
1383     SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
1384     int nbFaces = SMDS_VolumeTool::NbFaces( type );
1385     TSequenceOfXYZ points(4);
1386     for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume
1387       if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
1388         continue;
1389       const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
1390       for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadrangle face
1391         points( p + 1 ) = P( pInd[ p ] + 1 );
1392       aQuality = std::max( aQuality, aspect2D.GetValue( points ));
1393     }
1394   }
1395   return aQuality;
1396 }
1397
1398 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
1399 {
1400   // the aspect ratio is in the range [1.0,infinity]
1401   // 1.0 = good
1402   // infinity = bad
1403   return Value / 1000.;
1404 }
1405
1406 SMDSAbs_ElementType AspectRatio3D::GetType() const
1407 {
1408   return SMDSAbs_Volume;
1409 }
1410
1411
1412 //================================================================================
1413 /*
1414   Class       : Warping
1415   Description : Functor for calculating warping
1416 */
1417 //================================================================================
1418
1419 bool Warping::IsApplicable( const SMDS_MeshElement* element ) const
1420 {
1421   return NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4;
1422 }
1423
1424 double Warping::GetValue( const TSequenceOfXYZ& P )
1425 {
1426   if ( P.size() != 4 )
1427     return 0;
1428
1429   gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.;
1430
1431   double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
1432   double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
1433   double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
1434   double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
1435
1436   double val = Max( Max( A1, A2 ), Max( A3, A4 ) );
1437
1438   const double eps = 0.1; // val is in degrees
1439
1440   return val < eps ? 0. : val;
1441 }
1442
1443 double Warping::ComputeA( const gp_XYZ& thePnt1,
1444                           const gp_XYZ& thePnt2,
1445                           const gp_XYZ& thePnt3,
1446                           const gp_XYZ& theG ) const
1447 {
1448   double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
1449   double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
1450   double L = Min( aLen1, aLen2 ) * 0.5;
1451   if ( L < theEps )
1452     return theInf;
1453
1454   gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
1455   gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
1456   gp_XYZ N  = GI.Crossed( GJ );
1457
1458   if ( N.Modulus() < gp::Resolution() )
1459     return M_PI / 2;
1460
1461   N.Normalize();
1462
1463   double H = ( thePnt2 - theG ).Dot( N );
1464   return asin( fabs( H / L ) ) * 180. / M_PI;
1465 }
1466
1467 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
1468 {
1469   // the warp is in the range [0.0,PI/2]
1470   // 0.0 = good (no warp)
1471   // PI/2 = bad  (face pliee)
1472   return Value;
1473 }
1474
1475 SMDSAbs_ElementType Warping::GetType() const
1476 {
1477   return SMDSAbs_Face;
1478 }
1479
1480
1481 //================================================================================
1482 /*
1483   Class       : Taper
1484   Description : Functor for calculating taper
1485 */
1486 //================================================================================
1487
1488 bool Taper::IsApplicable( const SMDS_MeshElement* element ) const
1489 {
1490   return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4 );
1491 }
1492
1493 double Taper::GetValue( const TSequenceOfXYZ& P )
1494 {
1495   if ( P.size() != 4 )
1496     return 0.;
1497
1498   // Compute taper
1499   double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) );
1500   double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) );
1501   double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) );
1502   double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) );
1503
1504   double JA = 0.25 * ( J1 + J2 + J3 + J4 );
1505   if ( JA <= theEps )
1506     return theInf;
1507
1508   double T1 = fabs( ( J1 - JA ) / JA );
1509   double T2 = fabs( ( J2 - JA ) / JA );
1510   double T3 = fabs( ( J3 - JA ) / JA );
1511   double T4 = fabs( ( J4 - JA ) / JA );
1512
1513   double val = Max( Max( T1, T2 ), Max( T3, T4 ) );
1514
1515   const double eps = 0.01;
1516
1517   return val < eps ? 0. : val;
1518 }
1519
1520 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
1521 {
1522   // the taper is in the range [0.0,1.0]
1523   // 0.0 = good (no taper)
1524   // 1.0 = bad  (les cotes opposes sont allignes)
1525   return Value;
1526 }
1527
1528 SMDSAbs_ElementType Taper::GetType() const
1529 {
1530   return SMDSAbs_Face;
1531 }
1532
1533 //================================================================================
1534 /*
1535   Class       : Skew
1536   Description : Functor for calculating skew in degrees
1537 */
1538 //================================================================================
1539
1540 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
1541 {
1542   gp_XYZ p12 = ( p2 + p1 ) / 2.;
1543   gp_XYZ p23 = ( p3 + p2 ) / 2.;
1544   gp_XYZ p31 = ( p3 + p1 ) / 2.;
1545
1546   gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
1547
1548   return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
1549 }
1550
1551 bool Skew::IsApplicable( const SMDS_MeshElement* element ) const
1552 {
1553   return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() <= 4 );
1554 }
1555
1556 double Skew::GetValue( const TSequenceOfXYZ& P )
1557 {
1558   if ( P.size() != 3 && P.size() != 4 )
1559     return 0.;
1560
1561   // Compute skew
1562   const double PI2 = M_PI / 2.;
1563   if ( P.size() == 3 )
1564   {
1565     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
1566     double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
1567     double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
1568
1569     return Max( A0, Max( A1, A2 ) ) * 180. / M_PI;
1570   }
1571   else
1572   {
1573     gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.;
1574     gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.;
1575     gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.;
1576     gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.;
1577
1578     gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
1579     double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
1580       ? 0. : fabs( PI2 - v1.Angle( v2 ) );
1581
1582     double val = A * 180. / M_PI;
1583
1584     const double eps = 0.1; // val is in degrees
1585
1586     return val < eps ? 0. : val;
1587   }
1588 }
1589
1590 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
1591 {
1592   // the skew is in the range [0.0,PI/2].
1593   // 0.0 = good
1594   // PI/2 = bad
1595   return Value;
1596 }
1597
1598 SMDSAbs_ElementType Skew::GetType() const
1599 {
1600   return SMDSAbs_Face;
1601 }
1602
1603
1604 //================================================================================
1605 /*
1606   Class       : Area
1607   Description : Functor for calculating area
1608 */
1609 //================================================================================
1610
1611 double Area::GetValue( const TSequenceOfXYZ& P )
1612 {
1613   double val = 0.0;
1614   if ( P.size() > 2 )
1615   {
1616     gp_Vec aVec1( P(2) - P(1) );
1617     gp_Vec aVec2( P(3) - P(1) );
1618     gp_Vec SumVec = aVec1 ^ aVec2;
1619
1620     for (size_t i=4; i<=P.size(); i++)
1621     {
1622       gp_Vec aVec1( P(i-1) - P(1) );
1623       gp_Vec aVec2( P(i  ) - P(1) );
1624       gp_Vec tmp = aVec1 ^ aVec2;
1625       SumVec.Add(tmp);
1626     }
1627     val = SumVec.Magnitude() * 0.5;
1628   }
1629   return val;
1630 }
1631
1632 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
1633 {
1634   // meaningless as it is not a quality control functor
1635   return Value;
1636 }
1637
1638 SMDSAbs_ElementType Area::GetType() const
1639 {
1640   return SMDSAbs_Face;
1641 }
1642
1643 //================================================================================
1644 /*
1645   Class       : Length
1646   Description : Functor for calculating length of edge
1647 */
1648 //================================================================================
1649
1650 double Length::GetValue( const TSequenceOfXYZ& P )
1651 {
1652   switch ( P.size() ) {
1653   case 2:  return getDistance( P( 1 ), P( 2 ) );
1654   case 3:  return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) );
1655   default: return 0.;
1656   }
1657 }
1658
1659 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
1660 {
1661   // meaningless as it is not quality control functor
1662   return Value;
1663 }
1664
1665 SMDSAbs_ElementType Length::GetType() const
1666 {
1667   return SMDSAbs_Edge;
1668 }
1669
1670 //================================================================================
1671 /*
1672   Class       : Length3D
1673   Description : Functor for calculating minimal length of element edge
1674 */
1675 //================================================================================
1676
1677 Length3D::Length3D():
1678   Length2D ( SMDSAbs_Volume )
1679 {
1680 }
1681
1682 //================================================================================
1683 /*
1684   Class       : Length2D
1685   Description : Functor for calculating minimal length of element edge
1686 */
1687 //================================================================================
1688
1689 Length2D::Length2D( SMDSAbs_ElementType type ):
1690   myType ( type )
1691 {
1692 }
1693
1694 bool Length2D::IsApplicable( const SMDS_MeshElement* element ) const
1695 {
1696   return ( NumericalFunctor::IsApplicable( element ) &&
1697            element->GetEntityType() != SMDSEntity_Polyhedra );
1698 }
1699
1700 double Length2D::GetValue( const TSequenceOfXYZ& P )
1701 {
1702   double aVal = 0;
1703   int len = P.size();
1704   SMDSAbs_EntityType aType = P.getElementEntity();
1705
1706   switch (aType) {
1707   case SMDSEntity_Edge:
1708     if (len == 2)
1709       aVal = getDistance( P( 1 ), P( 2 ) );
1710     break;
1711   case SMDSEntity_Quad_Edge:
1712     if (len == 3) // quadratic edge
1713       aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
1714     break;
1715   case SMDSEntity_Triangle:
1716     if (len == 3){ // triangles
1717       double L1 = getDistance(P( 1 ),P( 2 ));
1718       double L2 = getDistance(P( 2 ),P( 3 ));
1719       double L3 = getDistance(P( 3 ),P( 1 ));
1720       aVal = Min(L1,Min(L2,L3));
1721     }
1722     break;
1723   case SMDSEntity_Quadrangle:
1724     if (len == 4){ // quadrangles
1725       double L1 = getDistance(P( 1 ),P( 2 ));
1726       double L2 = getDistance(P( 2 ),P( 3 ));
1727       double L3 = getDistance(P( 3 ),P( 4 ));
1728       double L4 = getDistance(P( 4 ),P( 1 ));
1729       aVal = Min(Min(L1,L2),Min(L3,L4));
1730     }
1731     break;
1732   case SMDSEntity_Quad_Triangle:
1733   case SMDSEntity_BiQuad_Triangle:
1734     if (len >= 6){ // quadratic triangles
1735       double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1736       double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1737       double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
1738       aVal = Min(L1,Min(L2,L3));
1739     }
1740     break;
1741   case SMDSEntity_Quad_Quadrangle:
1742   case SMDSEntity_BiQuad_Quadrangle:
1743     if (len >= 8){ // quadratic quadrangles
1744       double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1745       double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1746       double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
1747       double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
1748       aVal = Min(Min(L1,L2),Min(L3,L4));
1749     }
1750     break;
1751   case SMDSEntity_Tetra:
1752     if (len == 4){ // tetrahedra
1753       double L1 = getDistance(P( 1 ),P( 2 ));
1754       double L2 = getDistance(P( 2 ),P( 3 ));
1755       double L3 = getDistance(P( 3 ),P( 1 ));
1756       double L4 = getDistance(P( 1 ),P( 4 ));
1757       double L5 = getDistance(P( 2 ),P( 4 ));
1758       double L6 = getDistance(P( 3 ),P( 4 ));
1759       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1760     }
1761     break;
1762   case SMDSEntity_Pyramid:
1763     if (len == 5){ // pyramid
1764       double L1 = getDistance(P( 1 ),P( 2 ));
1765       double L2 = getDistance(P( 2 ),P( 3 ));
1766       double L3 = getDistance(P( 3 ),P( 4 ));
1767       double L4 = getDistance(P( 4 ),P( 1 ));
1768       double L5 = getDistance(P( 1 ),P( 5 ));
1769       double L6 = getDistance(P( 2 ),P( 5 ));
1770       double L7 = getDistance(P( 3 ),P( 5 ));
1771       double L8 = getDistance(P( 4 ),P( 5 ));
1772
1773       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1774       aVal = Min(aVal,Min(L7,L8));
1775     }
1776     break;
1777   case SMDSEntity_Penta:
1778     if (len == 6) { // pentahedron
1779       double L1 = getDistance(P( 1 ),P( 2 ));
1780       double L2 = getDistance(P( 2 ),P( 3 ));
1781       double L3 = getDistance(P( 3 ),P( 1 ));
1782       double L4 = getDistance(P( 4 ),P( 5 ));
1783       double L5 = getDistance(P( 5 ),P( 6 ));
1784       double L6 = getDistance(P( 6 ),P( 4 ));
1785       double L7 = getDistance(P( 1 ),P( 4 ));
1786       double L8 = getDistance(P( 2 ),P( 5 ));
1787       double L9 = getDistance(P( 3 ),P( 6 ));
1788
1789       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1790       aVal = Min(aVal,Min(Min(L7,L8),L9));
1791     }
1792     break;
1793   case SMDSEntity_Hexa:
1794     if (len == 8){ // hexahedron
1795       double L1 = getDistance(P( 1 ),P( 2 ));
1796       double L2 = getDistance(P( 2 ),P( 3 ));
1797       double L3 = getDistance(P( 3 ),P( 4 ));
1798       double L4 = getDistance(P( 4 ),P( 1 ));
1799       double L5 = getDistance(P( 5 ),P( 6 ));
1800       double L6 = getDistance(P( 6 ),P( 7 ));
1801       double L7 = getDistance(P( 7 ),P( 8 ));
1802       double L8 = getDistance(P( 8 ),P( 5 ));
1803       double L9 = getDistance(P( 1 ),P( 5 ));
1804       double L10= getDistance(P( 2 ),P( 6 ));
1805       double L11= getDistance(P( 3 ),P( 7 ));
1806       double L12= getDistance(P( 4 ),P( 8 ));
1807
1808       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1809       aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
1810       aVal = Min(aVal,Min(L11,L12));
1811     }
1812     break;
1813   case SMDSEntity_Quad_Tetra:
1814     if (len == 10){ // quadratic tetrahedron
1815       double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
1816       double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
1817       double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
1818       double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1819       double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
1820       double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
1821       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1822     }
1823     break;
1824   case SMDSEntity_Quad_Pyramid:
1825     if (len == 13){ // quadratic pyramid
1826       double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
1827       double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
1828       double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1829       double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1830       double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1831       double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
1832       double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
1833       double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
1834       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1835       aVal = Min(aVal,Min(L7,L8));
1836     }
1837     break;
1838   case SMDSEntity_Quad_Penta:
1839   case SMDSEntity_BiQuad_Penta:
1840     if (len >= 15){ // quadratic pentahedron
1841       double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
1842       double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
1843       double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1844       double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1845       double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
1846       double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
1847       double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
1848       double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
1849       double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
1850       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1851       aVal = Min(aVal,Min(Min(L7,L8),L9));
1852     }
1853     break;
1854   case SMDSEntity_Quad_Hexa:
1855   case SMDSEntity_TriQuad_Hexa:
1856     if (len >= 20) { // quadratic hexahedron
1857       double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
1858       double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
1859       double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
1860       double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
1861       double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
1862       double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
1863       double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
1864       double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
1865       double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
1866       double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
1867       double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
1868       double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
1869       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1870       aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
1871       aVal = Min(aVal,Min(L11,L12));
1872     }
1873     break;
1874   case SMDSEntity_Polygon:
1875     if ( len > 1 ) {
1876       aVal = getDistance( P(1), P( P.size() ));
1877       for ( size_t i = 1; i < P.size(); ++i )
1878         aVal = Min( aVal, getDistance( P( i ), P( i+1 )));
1879     }
1880     break;
1881   case SMDSEntity_Quad_Polygon:
1882     if ( len > 2 ) {
1883       aVal = getDistance( P(1), P( P.size() )) + getDistance( P(P.size()), P( P.size()-1 ));
1884       for ( size_t i = 1; i < P.size()-1; i += 2 )
1885         aVal = Min( aVal, getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 )));
1886     }
1887     break;
1888   case SMDSEntity_Hexagonal_Prism:
1889     if (len == 12) { // hexagonal prism
1890       double L1 = getDistance(P( 1 ),P( 2 ));
1891       double L2 = getDistance(P( 2 ),P( 3 ));
1892       double L3 = getDistance(P( 3 ),P( 4 ));
1893       double L4 = getDistance(P( 4 ),P( 5 ));
1894       double L5 = getDistance(P( 5 ),P( 6 ));
1895       double L6 = getDistance(P( 6 ),P( 1 ));
1896
1897       double L7 = getDistance(P( 7 ), P( 8 ));
1898       double L8 = getDistance(P( 8 ), P( 9 ));
1899       double L9 = getDistance(P( 9 ), P( 10 ));
1900       double L10= getDistance(P( 10 ),P( 11 ));
1901       double L11= getDistance(P( 11 ),P( 12 ));
1902       double L12= getDistance(P( 12 ),P( 7 ));
1903
1904       double L13 = getDistance(P( 1 ),P( 7 ));
1905       double L14 = getDistance(P( 2 ),P( 8 ));
1906       double L15 = getDistance(P( 3 ),P( 9 ));
1907       double L16 = getDistance(P( 4 ),P( 10 ));
1908       double L17 = getDistance(P( 5 ),P( 11 ));
1909       double L18 = getDistance(P( 6 ),P( 12 ));
1910       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1911       aVal = Min(aVal, Min(Min(Min(L7,L8),Min(L9,L10)),Min(L11,L12)));
1912       aVal = Min(aVal, Min(Min(Min(L13,L14),Min(L15,L16)),Min(L17,L18)));
1913     }
1914     break;
1915   case SMDSEntity_Polyhedra:
1916   {
1917   }
1918   break;
1919   default:
1920     return 0;
1921   }
1922
1923   if (aVal < 0 ) {
1924     return 0.;
1925   }
1926
1927   if ( myPrecision >= 0 )
1928   {
1929     double prec = pow( 10., (double)( myPrecision ) );
1930     aVal = floor( aVal * prec + 0.5 ) / prec;
1931   }
1932
1933   return aVal;
1934 }
1935
1936 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1937 {
1938   // meaningless as it is not a quality control functor
1939   return Value;
1940 }
1941
1942 SMDSAbs_ElementType Length2D::GetType() const
1943 {
1944   return myType;
1945 }
1946
1947 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
1948   myLength(theLength)
1949 {
1950   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
1951   if(thePntId1 > thePntId2){
1952     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
1953   }
1954 }
1955
1956 bool Length2D::Value::operator<(const Length2D::Value& x) const
1957 {
1958   if(myPntId[0] < x.myPntId[0]) return true;
1959   if(myPntId[0] == x.myPntId[0])
1960     if(myPntId[1] < x.myPntId[1]) return true;
1961   return false;
1962 }
1963
1964 void Length2D::GetValues(TValues& theValues)
1965 {
1966   if ( myType == SMDSAbs_Face )
1967   {
1968     for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
1969     {
1970       const SMDS_MeshFace* anElem = anIter->next();
1971       if ( anElem->IsQuadratic() )
1972       {
1973         // use special nodes iterator
1974         SMDS_NodeIteratorPtr anIter = anElem->interlacedNodesIterator();
1975         smIdType aNodeId[4] = { 0,0,0,0 };
1976         gp_Pnt P[4];
1977
1978         double aLength = 0;
1979         if ( anIter->more() )
1980         {
1981           const SMDS_MeshNode* aNode = anIter->next();
1982           P[0] = P[1] = SMESH_NodeXYZ( aNode );
1983           aNodeId[0] = aNodeId[1] = aNode->GetID();
1984           aLength = 0;
1985         }
1986         for ( ; anIter->more(); )
1987         {
1988           const SMDS_MeshNode* N1 = anIter->next();
1989           P[2] = SMESH_NodeXYZ( N1 );
1990           aNodeId[2] = N1->GetID();
1991           aLength = P[1].Distance(P[2]);
1992           if(!anIter->more()) break;
1993           const SMDS_MeshNode* N2 = anIter->next();
1994           P[3] = SMESH_NodeXYZ( N2 );
1995           aNodeId[3] = N2->GetID();
1996           aLength += P[2].Distance(P[3]);
1997           Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1998           Value aValue2(aLength,aNodeId[2],aNodeId[3]);
1999           P[1] = P[3];
2000           aNodeId[1] = aNodeId[3];
2001           theValues.insert(aValue1);
2002           theValues.insert(aValue2);
2003         }
2004         aLength += P[2].Distance(P[0]);
2005         Value aValue1(aLength,aNodeId[1],aNodeId[2]);
2006         Value aValue2(aLength,aNodeId[2],aNodeId[0]);
2007         theValues.insert(aValue1);
2008         theValues.insert(aValue2);
2009       }
2010       else {
2011         SMDS_NodeIteratorPtr aNodesIter = anElem->nodeIterator();
2012         smIdType aNodeId[2] = {0,0};
2013         gp_Pnt P[3];
2014
2015         double aLength;
2016         const SMDS_MeshElement* aNode;
2017         if ( aNodesIter->more())
2018         {
2019           aNode = aNodesIter->next();
2020           P[0] = P[1] = SMESH_NodeXYZ( aNode );
2021           aNodeId[0] = aNodeId[1] = aNode->GetID();
2022           aLength = 0;
2023         }
2024         for( ; aNodesIter->more(); )
2025         {
2026           aNode = aNodesIter->next();
2027           smIdType anId = aNode->GetID();
2028
2029           P[2] = SMESH_NodeXYZ( aNode );
2030
2031           aLength = P[1].Distance(P[2]);
2032
2033           Value aValue(aLength,aNodeId[1],anId);
2034           aNodeId[1] = anId;
2035           P[1] = P[2];
2036           theValues.insert(aValue);
2037         }
2038
2039         aLength = P[0].Distance(P[1]);
2040
2041         Value aValue(aLength,aNodeId[0],aNodeId[1]);
2042         theValues.insert(aValue);
2043       }
2044     }
2045   }
2046   else
2047   {
2048     // not implemented
2049   }
2050 }
2051
2052 //================================================================================
2053 /*
2054   Class       : Deflection2D
2055   Description : computes distance between a face center and an underlying surface
2056 */
2057 //================================================================================
2058
2059 double Deflection2D::GetValue( const TSequenceOfXYZ& P )
2060 {
2061   if ( myMesh && P.getElement() )
2062   {
2063     // get underlying surface
2064     if ( myShapeIndex != P.getElement()->getshapeId() )
2065     {
2066       mySurface.Nullify();
2067       myShapeIndex = P.getElement()->getshapeId();
2068       const TopoDS_Shape& S =
2069         static_cast< const SMESHDS_Mesh* >( myMesh )->IndexToShape( myShapeIndex );
2070       if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE )
2071       {
2072         mySurface = new ShapeAnalysis_Surface( BRep_Tool::Surface( TopoDS::Face( S )));
2073
2074         GeomLib_IsPlanarSurface isPlaneCheck( mySurface->Surface() );
2075         if ( isPlaneCheck.IsPlanar() )
2076           myPlane.reset( new gp_Pln( isPlaneCheck.Plan() ));
2077         else
2078           myPlane.reset();
2079       }
2080     }
2081     // project gravity center to the surface
2082     if ( !mySurface.IsNull() )
2083     {
2084       gp_XYZ gc(0,0,0);
2085       gp_XY  uv(0,0);
2086       int nbUV = 0;
2087       for ( size_t i = 0; i < P.size(); ++i )
2088       {
2089         gc += P(i+1);
2090
2091         if ( SMDS_FacePositionPtr fPos = P.getElement()->GetNode( i )->GetPosition() )
2092         {
2093           uv.ChangeCoord(1) += fPos->GetUParameter();
2094           uv.ChangeCoord(2) += fPos->GetVParameter();
2095           ++nbUV;
2096         }
2097       }
2098       gc /= P.size();
2099       if ( nbUV ) uv /= nbUV;
2100
2101       double maxLen = MaxElementLength2D().GetValue( P );
2102       double    tol = 1e-3 * maxLen;
2103       double dist;
2104       if ( myPlane )
2105       {
2106         dist = myPlane->Distance( gc );
2107         if ( dist < tol )
2108           dist = 0;
2109       }
2110       else
2111       {
2112         if ( uv.X() != 0 && uv.Y() != 0 ) // faster way
2113           mySurface->NextValueOfUV( uv, gc, tol, 0.5 * maxLen );
2114         else
2115           mySurface->ValueOfUV( gc, tol );
2116         dist = mySurface->Gap();
2117       }
2118       return Round( dist );
2119     }
2120   }
2121   return 0;
2122 }
2123
2124 void Deflection2D::SetMesh( const SMDS_Mesh* theMesh )
2125 {
2126   NumericalFunctor::SetMesh( dynamic_cast<const SMESHDS_Mesh* >( theMesh ));
2127   myShapeIndex = -100;
2128   myPlane.reset();
2129 }
2130
2131 SMDSAbs_ElementType Deflection2D::GetType() const
2132 {
2133   return SMDSAbs_Face;
2134 }
2135
2136 double Deflection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
2137 {
2138   // meaningless as it is not quality control functor
2139   return Value;
2140 }
2141
2142 //================================================================================
2143 /*
2144   Class       : MultiConnection
2145   Description : Functor for calculating number of faces conneted to the edge
2146 */
2147 //================================================================================
2148
2149 double MultiConnection::GetValue( const TSequenceOfXYZ& /*P*/ )
2150 {
2151   return 0;
2152 }
2153 double MultiConnection::GetValue( long theId )
2154 {
2155   return getNbMultiConnection( myMesh, theId );
2156 }
2157
2158 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
2159 {
2160   // meaningless as it is not quality control functor
2161   return Value;
2162 }
2163
2164 SMDSAbs_ElementType MultiConnection::GetType() const
2165 {
2166   return SMDSAbs_Edge;
2167 }
2168
2169 //================================================================================
2170 /*
2171   Class       : MultiConnection2D
2172   Description : Functor for calculating number of faces conneted to the edge
2173 */
2174 //================================================================================
2175
2176 double MultiConnection2D::GetValue( const TSequenceOfXYZ& /*P*/ )
2177 {
2178   return 0;
2179 }
2180
2181 double MultiConnection2D::GetValue( long theElementId )
2182 {
2183   int aResult = 0;
2184
2185   const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
2186   SMDSAbs_ElementType aType = aFaceElem->GetType();
2187
2188   switch (aType) {
2189   case SMDSAbs_Face:
2190     {
2191       int i = 0, len = aFaceElem->NbNodes();
2192       SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
2193       if (!anIter) break;
2194
2195       const SMDS_MeshNode *aNode, *aNode0 = 0;
2196       NCollection_Map< smIdType, smIdHasher > aMap, aMapPrev;
2197
2198       for (i = 0; i <= len; i++) {
2199         aMapPrev = aMap;
2200         aMap.Clear();
2201
2202         int aNb = 0;
2203         if (anIter->more()) {
2204           aNode = (SMDS_MeshNode*)anIter->next();
2205         } else {
2206           if (i == len)
2207             aNode = aNode0;
2208           else
2209             break;
2210         }
2211         if (!aNode) break;
2212         if (i == 0) aNode0 = aNode;
2213
2214         SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
2215         while (anElemIter->more()) {
2216           const SMDS_MeshElement* anElem = anElemIter->next();
2217           if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
2218             smIdType anId = anElem->GetID();
2219
2220             aMap.Add(anId);
2221             if (aMapPrev.Contains(anId)) {
2222               aNb++;
2223             }
2224           }
2225         }
2226         aResult = Max(aResult, aNb);
2227       }
2228     }
2229     break;
2230   default:
2231     aResult = 0;
2232   }
2233
2234   return aResult;
2235 }
2236
2237 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
2238 {
2239   // meaningless as it is not quality control functor
2240   return Value;
2241 }
2242
2243 SMDSAbs_ElementType MultiConnection2D::GetType() const
2244 {
2245   return SMDSAbs_Face;
2246 }
2247
2248 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
2249 {
2250   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
2251   if(thePntId1 > thePntId2){
2252     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
2253   }
2254 }
2255
2256 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const
2257 {
2258   if(myPntId[0] < x.myPntId[0]) return true;
2259   if(myPntId[0] == x.myPntId[0])
2260     if(myPntId[1] < x.myPntId[1]) return true;
2261   return false;
2262 }
2263
2264 void MultiConnection2D::GetValues(MValues& theValues)
2265 {
2266   if ( !myMesh ) return;
2267   for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
2268   {
2269     const SMDS_MeshFace*     anElem = anIter->next();
2270     SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
2271
2272     const SMDS_MeshNode* aNode1 = anElem->GetNode( anElem->NbNodes() - 1 );
2273     const SMDS_MeshNode* aNode2;
2274     for ( ; aNodesIter->more(); )
2275     {
2276       aNode2 = aNodesIter->next();
2277
2278       Value aValue ( aNode1->GetID(), aNode2->GetID() );
2279       MValues::iterator aItr = theValues.insert( std::make_pair( aValue, 0 )).first;
2280       aItr->second++;
2281       aNode1 = aNode2;
2282     }
2283   }
2284   return;
2285 }
2286
2287 //================================================================================
2288 /*
2289   Class       : BallDiameter
2290   Description : Functor returning diameter of a ball element
2291 */
2292 //================================================================================
2293
2294 double BallDiameter::GetValue( long theId )
2295 {
2296   double diameter = 0;
2297
2298   if ( const SMDS_BallElement* ball =
2299        myMesh->DownCast< SMDS_BallElement >( myMesh->FindElement( theId )))
2300   {
2301     diameter = ball->GetDiameter();
2302   }
2303   return diameter;
2304 }
2305
2306 double BallDiameter::GetBadRate( double Value, int /*nbNodes*/ ) const
2307 {
2308   // meaningless as it is not a quality control functor
2309   return Value;
2310 }
2311
2312 SMDSAbs_ElementType BallDiameter::GetType() const
2313 {
2314   return SMDSAbs_Ball;
2315 }
2316
2317 //================================================================================
2318 /*
2319   Class       : NodeConnectivityNumber
2320   Description : Functor returning number of elements connected to a node
2321 */
2322 //================================================================================
2323
2324 double NodeConnectivityNumber::GetValue( long theId )
2325 {
2326   double nb = 0;
2327
2328   if ( const SMDS_MeshNode* node = myMesh->FindNode( theId ))
2329   {
2330     SMDSAbs_ElementType type;
2331     if ( myMesh->NbVolumes() > 0 )
2332       type = SMDSAbs_Volume;
2333     else if ( myMesh->NbFaces() > 0 )
2334       type = SMDSAbs_Face;
2335     else if ( myMesh->NbEdges() > 0 )
2336       type = SMDSAbs_Edge;
2337     else
2338       return 0;
2339     nb = node->NbInverseElements( type );
2340   }
2341   return nb;
2342 }
2343
2344 double NodeConnectivityNumber::GetBadRate( double Value, int /*nbNodes*/ ) const
2345 {
2346   return Value;
2347 }
2348
2349 SMDSAbs_ElementType NodeConnectivityNumber::GetType() const
2350 {
2351   return SMDSAbs_Node;
2352 }
2353
2354 /*
2355                             PREDICATES
2356 */
2357
2358 //================================================================================
2359 /*
2360   Class       : BadOrientedVolume
2361   Description : Predicate bad oriented volumes
2362 */
2363 //================================================================================
2364
2365 BadOrientedVolume::BadOrientedVolume()
2366 {
2367   myMesh = 0;
2368 }
2369
2370 void BadOrientedVolume::SetMesh( const SMDS_Mesh* theMesh )
2371 {
2372   myMesh = theMesh;
2373 }
2374
2375 bool BadOrientedVolume::IsSatisfy( long theId )
2376 {
2377   if ( myMesh == 0 )
2378     return false;
2379
2380   SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
2381
2382   bool isOk = true;
2383   if ( vTool.IsPoly() )
2384   {
2385     isOk = true;
2386     for ( int i = 0; i < vTool.NbFaces() && isOk; ++i )
2387       isOk = vTool.IsFaceExternal( i );
2388   }
2389   else
2390   {
2391     isOk = vTool.IsForward();
2392   }
2393   return !isOk;
2394 }
2395
2396 SMDSAbs_ElementType BadOrientedVolume::GetType() const
2397 {
2398   return SMDSAbs_Volume;
2399 }
2400
2401 /*
2402   Class       : BareBorderVolume
2403 */
2404
2405 bool BareBorderVolume::IsSatisfy(long theElementId )
2406 {
2407   SMDS_VolumeTool  myTool;
2408   if ( myTool.Set( myMesh->FindElement(theElementId)))
2409   {
2410     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2411       if ( myTool.IsFreeFace( iF ))
2412       {
2413         const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
2414         std::vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
2415         if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
2416           return true;
2417       }
2418   }
2419   return false;
2420 }
2421
2422 //================================================================================
2423 /*
2424   Class       : BareBorderFace
2425 */
2426 //================================================================================
2427
2428 bool BareBorderFace::IsSatisfy(long theElementId )
2429 {
2430   bool ok = false;
2431   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2432   {
2433     if ( face->GetType() == SMDSAbs_Face )
2434     {
2435       int nbN = face->NbCornerNodes();
2436       for ( int i = 0; i < nbN && !ok; ++i )
2437       {
2438         // check if a link is shared by another face
2439         const SMDS_MeshNode* n1 = face->GetNode( i );
2440         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2441         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2442         bool isShared = false;
2443         while ( !isShared && fIt->more() )
2444         {
2445           const SMDS_MeshElement* f = fIt->next();
2446           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2447         }
2448         if ( !isShared )
2449         {
2450           const int iQuad = face->IsQuadratic();
2451           myLinkNodes.resize( 2 + iQuad);
2452           myLinkNodes[0] = n1;
2453           myLinkNodes[1] = n2;
2454           if ( iQuad )
2455             myLinkNodes[2] = face->GetNode( i+nbN );
2456           ok = !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
2457         }
2458       }
2459     }
2460   }
2461   return ok;
2462 }
2463
2464 //================================================================================
2465 /*
2466   Class       : OverConstrainedVolume
2467 */
2468 //================================================================================
2469
2470 bool OverConstrainedVolume::IsSatisfy(long theElementId )
2471 {
2472   // An element is over-constrained if it has N-1 free borders where
2473   // N is the number of edges/faces for a 2D/3D element.
2474   SMDS_VolumeTool  myTool;
2475   if ( myTool.Set( myMesh->FindElement(theElementId)))
2476   {
2477     int nbSharedFaces = 0;
2478     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2479       if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
2480         break;
2481     return ( nbSharedFaces == 1 );
2482   }
2483   return false;
2484 }
2485
2486 //================================================================================
2487 /*
2488   Class       : OverConstrainedFace
2489 */
2490 //================================================================================
2491
2492 bool OverConstrainedFace::IsSatisfy(long theElementId )
2493 {
2494   // An element is over-constrained if it has N-1 free borders where
2495   // N is the number of edges/faces for a 2D/3D element.
2496   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2497     if ( face->GetType() == SMDSAbs_Face )
2498     {
2499       int nbSharedBorders = 0;
2500       int nbN = face->NbCornerNodes();
2501       for ( int i = 0; i < nbN; ++i )
2502       {
2503         // check if a link is shared by another face
2504         const SMDS_MeshNode* n1 = face->GetNode( i );
2505         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2506         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2507         bool isShared = false;
2508         while ( !isShared && fIt->more() )
2509         {
2510           const SMDS_MeshElement* f = fIt->next();
2511           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2512         }
2513         if ( isShared && ++nbSharedBorders > 1 )
2514           break;
2515       }
2516       return ( nbSharedBorders == 1 );
2517     }
2518   return false;
2519 }
2520
2521 //================================================================================
2522 /*
2523   Class       : CoincidentNodes
2524   Description : Predicate of Coincident nodes
2525 */
2526 //================================================================================
2527
2528 CoincidentNodes::CoincidentNodes()
2529 {
2530   myToler = 1e-5;
2531 }
2532
2533 bool CoincidentNodes::IsSatisfy( long theElementId )
2534 {
2535   return myCoincidentIDs.Contains( theElementId );
2536 }
2537
2538 SMDSAbs_ElementType CoincidentNodes::GetType() const
2539 {
2540   return SMDSAbs_Node;
2541 }
2542
2543 void CoincidentNodes::SetTolerance( const double theToler )
2544 {
2545   if ( myToler != theToler )
2546   {
2547     SetMesh(0);
2548     myToler = theToler;
2549   }
2550 }
2551
2552 void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
2553 {
2554   myMeshModifTracer.SetMesh( theMesh );
2555   if ( myMeshModifTracer.IsMeshModified() )
2556   {
2557     TIDSortedNodeSet nodesToCheck;
2558     SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
2559     while ( nIt->more() )
2560       nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
2561
2562     std::list< std::list< const SMDS_MeshNode*> > nodeGroups;
2563     SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
2564
2565     myCoincidentIDs.Clear();
2566     std::list< std::list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
2567     for ( ; groupIt != nodeGroups.end(); ++groupIt )
2568     {
2569       std::list< const SMDS_MeshNode*>& coincNodes = *groupIt;
2570       std::list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
2571       for ( ; n != coincNodes.end(); ++n )
2572         myCoincidentIDs.Add( (*n)->GetID() );
2573     }
2574   }
2575 }
2576
2577 //================================================================================
2578 /*
2579   Class       : CoincidentElements
2580   Description : Predicate of Coincident Elements
2581   Note        : This class is suitable only for visualization of Coincident Elements
2582 */
2583 //================================================================================
2584
2585 CoincidentElements::CoincidentElements()
2586 {
2587   myMesh = 0;
2588 }
2589
2590 void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
2591 {
2592   myMesh = theMesh;
2593 }
2594
2595 bool CoincidentElements::IsSatisfy( long theElementId )
2596 {
2597   if ( !myMesh ) return false;
2598
2599   if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
2600   {
2601     if ( e->GetType() != GetType() ) return false;
2602     std::set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
2603     const int nbNodes = e->NbNodes();
2604     SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
2605     while ( invIt->more() )
2606     {
2607       const SMDS_MeshElement* e2 = invIt->next();
2608       if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
2609
2610       bool sameNodes = true;
2611       for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
2612         sameNodes = ( elemNodes.count( e2->GetNode( i )));
2613       if ( sameNodes )
2614         return true;
2615     }
2616   }
2617   return false;
2618 }
2619
2620 SMDSAbs_ElementType CoincidentElements1D::GetType() const
2621 {
2622   return SMDSAbs_Edge;
2623 }
2624 SMDSAbs_ElementType CoincidentElements2D::GetType() const
2625 {
2626   return SMDSAbs_Face;
2627 }
2628 SMDSAbs_ElementType CoincidentElements3D::GetType() const
2629 {
2630   return SMDSAbs_Volume;
2631 }
2632
2633
2634 //================================================================================
2635 /*
2636   Class       : FreeBorders
2637   Description : Predicate for free borders
2638 */
2639 //================================================================================
2640
2641 FreeBorders::FreeBorders()
2642 {
2643   myMesh = 0;
2644 }
2645
2646 void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
2647 {
2648   myMesh = theMesh;
2649 }
2650
2651 bool FreeBorders::IsSatisfy( long theId )
2652 {
2653   return getNbMultiConnection( myMesh, theId ) == 1;
2654 }
2655
2656 SMDSAbs_ElementType FreeBorders::GetType() const
2657 {
2658   return SMDSAbs_Edge;
2659 }
2660
2661
2662 //================================================================================
2663 /*
2664   Class       : FreeEdges
2665   Description : Predicate for free Edges
2666 */
2667 //================================================================================
2668
2669 FreeEdges::FreeEdges()
2670 {
2671   myMesh = 0;
2672 }
2673
2674 void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
2675 {
2676   myMesh = theMesh;
2677 }
2678
2679 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const smIdType theFaceId  )
2680 {
2681   SMDS_ElemIteratorPtr anElemIter = theNodes[ 0 ]->GetInverseElementIterator(SMDSAbs_Face);
2682   while( anElemIter->more() )
2683   {
2684     if ( const SMDS_MeshElement* anElem = anElemIter->next())
2685     {
2686       const smIdType anId = anElem->GetID();
2687       if ( anId != theFaceId && anElem->GetNodeIndex( theNodes[1] ) >= 0 )
2688         return false;
2689     }
2690   }
2691   return true;
2692 }
2693
2694 bool FreeEdges::IsSatisfy( long theId )
2695 {
2696   if ( myMesh == 0 )
2697     return false;
2698
2699   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2700   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
2701     return false;
2702
2703   SMDS_NodeIteratorPtr anIter = aFace->interlacedNodesIterator();
2704   if ( !anIter )
2705     return false;
2706
2707   int i = 0, nbNodes = aFace->NbNodes();
2708   std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
2709   while( anIter->more() )
2710     if ( ! ( aNodes[ i++ ] = anIter->next() ))
2711       return false;
2712   aNodes[ nbNodes ] = aNodes[ 0 ];
2713
2714   for ( i = 0; i < nbNodes; i++ )
2715     if ( IsFreeEdge( &aNodes[ i ], theId ) )
2716       return true;
2717
2718   return false;
2719 }
2720
2721 SMDSAbs_ElementType FreeEdges::GetType() const
2722 {
2723   return SMDSAbs_Face;
2724 }
2725
2726 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
2727   myElemId(theElemId)
2728 {
2729   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
2730   if(thePntId1 > thePntId2){
2731     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
2732   }
2733 }
2734
2735 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
2736   if(myPntId[0] < x.myPntId[0]) return true;
2737   if(myPntId[0] == x.myPntId[0])
2738     if(myPntId[1] < x.myPntId[1]) return true;
2739   return false;
2740 }
2741
2742 inline void UpdateBorders(const FreeEdges::Border& theBorder,
2743                           FreeEdges::TBorders& theRegistry,
2744                           FreeEdges::TBorders& theContainer)
2745 {
2746   if(theRegistry.find(theBorder) == theRegistry.end()){
2747     theRegistry.insert(theBorder);
2748     theContainer.insert(theBorder);
2749   }else{
2750     theContainer.erase(theBorder);
2751   }
2752 }
2753
2754 void FreeEdges::GetBoreders(TBorders& theBorders)
2755 {
2756   TBorders aRegistry;
2757   for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
2758   {
2759     const SMDS_MeshFace* anElem = anIter->next();
2760     long anElemId = anElem->GetID();
2761     SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
2762     if ( !aNodesIter->more() ) continue;
2763     long aNodeId[2] = {0,0};
2764     aNodeId[0] = anElem->GetNode( anElem->NbNodes()-1 )->GetID();
2765     for ( ; aNodesIter->more(); )
2766     {
2767       aNodeId[1] = aNodesIter->next()->GetID();
2768       Border aBorder( anElemId, aNodeId[0], aNodeId[1] );
2769       UpdateBorders( aBorder, aRegistry, theBorders );
2770       aNodeId[0] = aNodeId[1];
2771     }
2772   }
2773 }
2774
2775 //================================================================================
2776 /*
2777   Class       : FreeNodes
2778   Description : Predicate for free nodes
2779 */
2780 //================================================================================
2781
2782 FreeNodes::FreeNodes()
2783 {
2784   myMesh = 0;
2785 }
2786
2787 void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
2788 {
2789   myMesh = theMesh;
2790 }
2791
2792 bool FreeNodes::IsSatisfy( long theNodeId )
2793 {
2794   const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
2795   if (!aNode)
2796     return false;
2797
2798   return (aNode->NbInverseElements() < 1);
2799 }
2800
2801 SMDSAbs_ElementType FreeNodes::GetType() const
2802 {
2803   return SMDSAbs_Node;
2804 }
2805
2806
2807 //================================================================================
2808 /*
2809   Class       : FreeFaces
2810   Description : Predicate for free faces
2811 */
2812 //================================================================================
2813
2814 FreeFaces::FreeFaces()
2815 {
2816   myMesh = 0;
2817 }
2818
2819 void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
2820 {
2821   myMesh = theMesh;
2822 }
2823
2824 bool FreeFaces::IsSatisfy( long theId )
2825 {
2826   if (!myMesh) return false;
2827   // check that faces nodes refers to less than two common volumes
2828   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2829   if ( !aFace || aFace->GetType() != SMDSAbs_Face )
2830     return false;
2831
2832   int nbNode = aFace->NbNodes();
2833
2834   // collect volumes to check that number of volumes with count equal nbNode not less than 2
2835   typedef std::map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
2836   typedef std::map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
2837   TMapOfVolume mapOfVol;
2838
2839   SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
2840   while ( nodeItr->more() )
2841   {
2842     const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
2843     if ( !aNode ) continue;
2844     SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
2845     while ( volItr->more() )
2846     {
2847       SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
2848       TItrMapOfVolume    itr = mapOfVol.insert( std::make_pair( aVol, 0 )).first;
2849       (*itr).second++;
2850     }
2851   }
2852   int nbVol = 0;
2853   TItrMapOfVolume volItr = mapOfVol.begin();
2854   TItrMapOfVolume volEnd = mapOfVol.end();
2855   for ( ; volItr != volEnd; ++volItr )
2856     if ( (*volItr).second >= nbNode )
2857       nbVol++;
2858   // face is not free if number of volumes constructed on their nodes more than one
2859   return (nbVol < 2);
2860 }
2861
2862 SMDSAbs_ElementType FreeFaces::GetType() const
2863 {
2864   return SMDSAbs_Face;
2865 }
2866
2867 //================================================================================
2868 /*
2869   Class       : LinearOrQuadratic
2870   Description : Predicate to verify whether a mesh element is linear
2871 */
2872 //================================================================================
2873
2874 LinearOrQuadratic::LinearOrQuadratic()
2875 {
2876   myMesh = 0;
2877 }
2878
2879 void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
2880 {
2881   myMesh = theMesh;
2882 }
2883
2884 bool LinearOrQuadratic::IsSatisfy( long theId )
2885 {
2886   if (!myMesh) return false;
2887   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2888   if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
2889     return false;
2890   return (!anElem->IsQuadratic());
2891 }
2892
2893 void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
2894 {
2895   myType = theType;
2896 }
2897
2898 SMDSAbs_ElementType LinearOrQuadratic::GetType() const
2899 {
2900   return myType;
2901 }
2902
2903 //================================================================================
2904 /*
2905   Class       : GroupColor
2906   Description : Functor for check color of group to which mesh element belongs to
2907 */
2908 //================================================================================
2909
2910 GroupColor::GroupColor()
2911 {
2912 }
2913
2914 bool GroupColor::IsSatisfy( long theId )
2915 {
2916   return myIDs.count( theId );
2917 }
2918
2919 void GroupColor::SetType( SMDSAbs_ElementType theType )
2920 {
2921   myType = theType;
2922 }
2923
2924 SMDSAbs_ElementType GroupColor::GetType() const
2925 {
2926   return myType;
2927 }
2928
2929 static bool isEqual( const Quantity_Color& theColor1,
2930                      const Quantity_Color& theColor2 )
2931 {
2932   // tolerance to compare colors
2933   const double tol = 5*1e-3;
2934   return ( fabs( theColor1.Red()   - theColor2.Red() )   < tol &&
2935            fabs( theColor1.Green() - theColor2.Green() ) < tol &&
2936            fabs( theColor1.Blue()  - theColor2.Blue() )  < tol );
2937 }
2938
2939 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
2940 {
2941   myIDs.clear();
2942
2943   const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
2944   if ( !aMesh )
2945     return;
2946
2947   int nbGrp = aMesh->GetNbGroups();
2948   if ( !nbGrp )
2949     return;
2950
2951   // iterates on groups and find necessary elements ids
2952   const std::set<SMESHDS_GroupBase*>&       aGroups = aMesh->GetGroups();
2953   std::set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
2954   for (; GrIt != aGroups.end(); GrIt++)
2955   {
2956     SMESHDS_GroupBase* aGrp = (*GrIt);
2957     if ( !aGrp )
2958       continue;
2959     // check type and color of group
2960     if ( !isEqual( myColor, aGrp->GetColor() ))
2961       continue;
2962
2963     // IPAL52867 (prevent infinite recursion via GroupOnFilter)
2964     if ( SMESHDS_GroupOnFilter * gof = dynamic_cast< SMESHDS_GroupOnFilter* >( aGrp ))
2965       if ( gof->GetPredicate().get() == this )
2966         continue;
2967
2968     SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
2969     if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
2970       // add elements IDS into control
2971       smIdType aSize = aGrp->Extent();
2972       for (smIdType i = 0; i < aSize; i++)
2973         myIDs.insert( aGrp->GetID(i+1) );
2974     }
2975   }
2976 }
2977
2978 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
2979 {
2980   Kernel_Utils::Localizer loc;
2981   TCollection_AsciiString aStr = theStr;
2982   aStr.RemoveAll( ' ' );
2983   aStr.RemoveAll( '\t' );
2984   for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
2985     aStr.Remove( aPos, 2 );
2986   Standard_Real clr[3];
2987   clr[0] = clr[1] = clr[2] = 0.;
2988   for ( int i = 0; i < 3; i++ ) {
2989     TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
2990     if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
2991       clr[i] = tmpStr.RealValue();
2992   }
2993   myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
2994 }
2995
2996 //=======================================================================
2997 // name    : GetRangeStr
2998 // Purpose : Get range as a string.
2999 //           Example: "1,2,3,50-60,63,67,70-"
3000 //=======================================================================
3001
3002 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
3003 {
3004   theResStr.Clear();
3005   theResStr += TCollection_AsciiString( myColor.Red() );
3006   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
3007   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
3008 }
3009
3010 //================================================================================
3011 /*
3012   Class       : ElemGeomType
3013   Description : Predicate to check element geometry type
3014 */
3015 //================================================================================
3016
3017 ElemGeomType::ElemGeomType()
3018 {
3019   myMesh = 0;
3020   myType = SMDSAbs_All;
3021   myGeomType = SMDSGeom_TRIANGLE;
3022 }
3023
3024 void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
3025 {
3026   myMesh = theMesh;
3027 }
3028
3029 bool ElemGeomType::IsSatisfy( long theId )
3030 {
3031   if (!myMesh) return false;
3032   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3033   if ( !anElem )
3034     return false;
3035   const SMDSAbs_ElementType anElemType = anElem->GetType();
3036   if ( myType != SMDSAbs_All && anElemType != myType )
3037     return false;
3038   bool isOk = ( anElem->GetGeomType() == myGeomType );
3039   return isOk;
3040 }
3041
3042 void ElemGeomType::SetType( SMDSAbs_ElementType theType )
3043 {
3044   myType = theType;
3045 }
3046
3047 SMDSAbs_ElementType ElemGeomType::GetType() const
3048 {
3049   return myType;
3050 }
3051
3052 void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
3053 {
3054   myGeomType = theType;
3055 }
3056
3057 SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
3058 {
3059   return myGeomType;
3060 }
3061
3062 //================================================================================
3063 /*
3064   Class       : ElemEntityType
3065   Description : Predicate to check element entity type
3066 */
3067 //================================================================================
3068
3069 ElemEntityType::ElemEntityType():
3070   myMesh( 0 ),
3071   myType( SMDSAbs_All ),
3072   myEntityType( SMDSEntity_0D )
3073 {
3074 }
3075
3076 void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh )
3077 {
3078   myMesh = theMesh;
3079 }
3080
3081 bool ElemEntityType::IsSatisfy( long theId )
3082 {
3083   if ( !myMesh ) return false;
3084   if ( myType == SMDSAbs_Node )
3085     return myMesh->FindNode( theId );
3086   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3087   return ( anElem &&
3088            myEntityType == anElem->GetEntityType() );
3089 }
3090
3091 void ElemEntityType::SetType( SMDSAbs_ElementType theType )
3092 {
3093   myType = theType;
3094 }
3095
3096 SMDSAbs_ElementType ElemEntityType::GetType() const
3097 {
3098   return myType;
3099 }
3100
3101 void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType )
3102 {
3103   myEntityType = theEntityType;
3104 }
3105
3106 SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const
3107 {
3108   return myEntityType;
3109 }
3110
3111 //================================================================================
3112 /*!
3113  * \brief Class ConnectedElements
3114  */
3115 //================================================================================
3116
3117 ConnectedElements::ConnectedElements():
3118   myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {}
3119
3120 SMDSAbs_ElementType ConnectedElements::GetType() const
3121 { return myType; }
3122
3123 smIdType ConnectedElements::GetNode() const
3124 { return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
3125
3126 std::vector<double> ConnectedElements::GetPoint() const
3127 { return myXYZ; }
3128
3129 void ConnectedElements::clearOkIDs()
3130 { myOkIDsReady = false; myOkIDs.clear(); }
3131
3132 void ConnectedElements::SetType( SMDSAbs_ElementType theType )
3133 {
3134   if ( myType != theType || myMeshModifTracer.IsMeshModified() )
3135     clearOkIDs();
3136   myType = theType;
3137 }
3138
3139 void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh )
3140 {
3141   myMeshModifTracer.SetMesh( theMesh );
3142   if ( myMeshModifTracer.IsMeshModified() )
3143   {
3144     clearOkIDs();
3145     if ( !myXYZ.empty() )
3146       SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh
3147   }
3148 }
3149
3150 void ConnectedElements::SetNode( smIdType nodeID )
3151 {
3152   myNodeID = nodeID;
3153   myXYZ.clear();
3154
3155   bool isSameDomain = false;
3156   if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() )
3157     if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID ))
3158     {
3159       SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType );
3160       while ( !isSameDomain && eIt->more() )
3161         isSameDomain = IsSatisfy( eIt->next()->GetID() );
3162     }
3163   if ( !isSameDomain )
3164     clearOkIDs();
3165 }
3166
3167 void ConnectedElements::SetPoint( double x, double y, double z )
3168 {
3169   myXYZ.resize(3);
3170   myXYZ[0] = x;
3171   myXYZ[1] = y;
3172   myXYZ[2] = z;
3173   myNodeID = 0;
3174
3175   bool isSameDomain = false;
3176
3177   // find myNodeID by myXYZ if possible
3178   if ( myMeshModifTracer.GetMesh() )
3179   {
3180     SMESHUtils::Deleter<SMESH_ElementSearcher> searcher
3181       ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
3182
3183     std::vector< const SMDS_MeshElement* > foundElems;
3184     searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
3185
3186     if ( !foundElems.empty() )
3187     {
3188       myNodeID = foundElems[0]->GetNode(0)->GetID();
3189       if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() )
3190         isSameDomain = IsSatisfy( foundElems[0]->GetID() );
3191     }
3192   }
3193   if ( !isSameDomain )
3194     clearOkIDs();
3195 }
3196
3197 bool ConnectedElements::IsSatisfy( long theElementId )
3198 {
3199   // Here we do NOT check if the mesh has changed, we do it in Set...() only!!!
3200
3201   if ( !myOkIDsReady )
3202   {
3203     if ( !myMeshModifTracer.GetMesh() )
3204       return false;
3205     const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID );
3206     if ( !node0 )
3207       return false;
3208
3209     std::list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
3210     std::set< smIdType > checkedNodeIDs;
3211     // algo:
3212     // foreach node in nodeQueue:
3213     //   foreach element sharing a node:
3214     //     add ID of an element of myType to myOkIDs;
3215     //     push all element nodes absent from checkedNodeIDs to nodeQueue;
3216     while ( !nodeQueue.empty() )
3217     {
3218       const SMDS_MeshNode* node = nodeQueue.front();
3219       nodeQueue.pop_front();
3220
3221       // loop on elements sharing the node
3222       SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
3223       while ( eIt->more() )
3224       {
3225         // keep elements of myType
3226         const SMDS_MeshElement* element = eIt->next();
3227         if ( myType == SMDSAbs_All || element->GetType() == myType )
3228           myOkIDs.insert( myOkIDs.end(), element->GetID() );
3229
3230         // enqueue nodes of the element
3231         SMDS_ElemIteratorPtr nIt = element->nodesIterator();
3232         while ( nIt->more() )
3233         {
3234           const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() );
3235           if ( checkedNodeIDs.insert( n->GetID()).second )
3236             nodeQueue.push_back( n );
3237         }
3238       }
3239     }
3240     if ( myType == SMDSAbs_Node )
3241       std::swap( myOkIDs, checkedNodeIDs );
3242
3243     size_t totalNbElems = myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType );
3244     if ( myOkIDs.size() == totalNbElems )
3245       myOkIDs.clear();
3246
3247     myOkIDsReady = true;
3248   }
3249
3250   return myOkIDs.empty() ? true : myOkIDs.count( theElementId );
3251 }
3252
3253 //================================================================================
3254 /*!
3255  * \brief Class CoplanarFaces
3256  */
3257 //================================================================================
3258
3259 namespace
3260 {
3261   inline bool isLessAngle( const gp_Vec& v1, const gp_Vec& v2, const double cos )
3262   {
3263     double dot = v1 * v2; // cos * |v1| * |v2|
3264     double l1  = v1.SquareMagnitude();
3265     double l2  = v2.SquareMagnitude();
3266     return (( dot * cos >= 0 ) &&
3267             ( dot * dot ) / l1 / l2 >= ( cos * cos ));
3268   }
3269 }
3270 CoplanarFaces::CoplanarFaces()
3271   : myFaceID(0), myToler(0)
3272 {
3273 }
3274 void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
3275 {
3276   myMeshModifTracer.SetMesh( theMesh );
3277   if ( myMeshModifTracer.IsMeshModified() )
3278   {
3279     // Build a set of coplanar face ids
3280
3281     myCoplanarIDs.Clear();
3282
3283     if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler )
3284       return;
3285
3286     const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID );
3287     if ( !face || face->GetType() != SMDSAbs_Face )
3288       return;
3289
3290     bool normOK;
3291     gp_Vec myNorm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
3292     if (!normOK)
3293       return;
3294
3295     const double cosTol = Cos( myToler * M_PI / 180. );
3296     NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks;
3297
3298     std::list< std::pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
3299     faceQueue.push_back( std::make_pair( face, myNorm ));
3300     while ( !faceQueue.empty() )
3301     {
3302       face   = faceQueue.front().first;
3303       myNorm = faceQueue.front().second;
3304       faceQueue.pop_front();
3305
3306       for ( int i = 0, nbN = face->NbCornerNodes(); i < nbN; ++i )
3307       {
3308         const SMDS_MeshNode*  n1 = face->GetNode( i );
3309         const SMDS_MeshNode*  n2 = face->GetNode(( i+1 )%nbN);
3310         if ( !checkedLinks.Add( SMESH_TLink( n1, n2 )))
3311           continue;
3312         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face);
3313         while ( fIt->more() )
3314         {
3315           const SMDS_MeshElement* f = fIt->next();
3316           if ( f->GetNodeIndex( n2 ) > -1 )
3317           {
3318             gp_Vec norm = getNormale( static_cast<const SMDS_MeshFace*>(f), &normOK );
3319             if (!normOK || isLessAngle( myNorm, norm, cosTol))
3320             {
3321               myCoplanarIDs.Add( f->GetID() );
3322               faceQueue.push_back( std::make_pair( f, norm ));
3323             }
3324           }
3325         }
3326       }
3327     }
3328   }
3329 }
3330 bool CoplanarFaces::IsSatisfy( long theElementId )
3331 {
3332   return myCoplanarIDs.Contains( theElementId );
3333 }
3334
3335 /*
3336  *Class       : RangeOfIds
3337   *Description : Predicate for Range of Ids.
3338   *              Range may be specified with two ways.
3339   *              1. Using AddToRange method
3340   *              2. With SetRangeStr method. Parameter of this method is a string
3341   *                 like as "1,2,3,50-60,63,67,70-"
3342 */
3343
3344 //=======================================================================
3345 // name    : RangeOfIds
3346 // Purpose : Constructor
3347 //=======================================================================
3348 RangeOfIds::RangeOfIds()
3349 {
3350   myMesh = 0;
3351   myType = SMDSAbs_All;
3352 }
3353
3354 //=======================================================================
3355 // name    : SetMesh
3356 // Purpose : Set mesh
3357 //=======================================================================
3358 void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh )
3359 {
3360   myMesh = theMesh;
3361 }
3362
3363 //=======================================================================
3364 // name    : AddToRange
3365 // Purpose : Add ID to the range
3366 //=======================================================================
3367 bool RangeOfIds::AddToRange( long theEntityId )
3368 {
3369   myIds.Add( theEntityId );
3370   return true;
3371 }
3372
3373 //=======================================================================
3374 // name    : GetRangeStr
3375 // Purpose : Get range as a string.
3376 //           Example: "1,2,3,50-60,63,67,70-"
3377 //=======================================================================
3378 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
3379 {
3380   theResStr.Clear();
3381
3382   TIDsSeq                             anIntSeq;
3383   NCollection_Sequence< std::string > aStrSeq;
3384
3385   TIDsMap::Iterator anIter( myIds );
3386   for ( ; anIter.More(); anIter.Next() )
3387   {
3388     smIdType anId = anIter.Key();
3389     SMESH_Comment aStr( anId );
3390     anIntSeq.Append( anId );
3391     aStrSeq.Append( aStr );
3392   }
3393
3394   for ( smIdType i = 1, n = myMin.size(); i <= n; i++ )
3395   {
3396     smIdType aMinId = myMin[i];
3397     smIdType aMaxId = myMax[i];
3398
3399     SMESH_Comment aStr;
3400     if ( aMinId != IntegerFirst() )
3401       aStr << aMinId;
3402
3403     aStr << "-";
3404
3405     if ( aMaxId != std::numeric_limits<smIdType>::max() )
3406       aStr << aMaxId;
3407
3408     // find position of the string in result sequence and insert string in it
3409     if ( anIntSeq.Length() == 0 )
3410     {
3411       anIntSeq.Append( aMinId );
3412       aStrSeq.Append( (const char*)aStr );
3413     }
3414     else
3415     {
3416       if ( aMinId < anIntSeq.First() )
3417       {
3418         anIntSeq.Prepend( aMinId );
3419         aStrSeq.Prepend( (const char*)aStr );
3420       }
3421       else if ( aMinId > anIntSeq.Last() )
3422       {
3423         anIntSeq.Append( aMinId );
3424         aStrSeq.Append( (const char*)aStr );
3425       }
3426       else
3427         for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
3428           if ( aMinId < anIntSeq( j ) )
3429           {
3430             anIntSeq.InsertBefore( j, aMinId );
3431             aStrSeq.InsertBefore( j, (const char*)aStr );
3432             break;
3433           }
3434     }
3435   }
3436
3437   if ( aStrSeq.Length() == 0 )
3438     return;
3439   std::string aResStr;
3440   aResStr = aStrSeq( 1 );
3441   for ( int j = 2, k = aStrSeq.Length(); j <= k; j++  )
3442   {
3443     aResStr += ",";
3444     aResStr += aStrSeq( j );
3445   }
3446   theResStr = aResStr.c_str();
3447 }
3448
3449 //=======================================================================
3450 // name    : SetRangeStr
3451 // Purpose : Define range with string
3452 //           Example of entry string: "1,2,3,50-60,63,67,70-"
3453 //=======================================================================
3454 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
3455 {
3456   myMin.clear();
3457   myMax.clear();
3458   myIds.Clear();
3459
3460   TCollection_AsciiString aStr = theStr;
3461   for ( int i = 1; i <= aStr.Length(); ++i )
3462   {
3463     char c = aStr.Value( i );
3464     if ( !isdigit( c ) && c != ',' && c != '-' )
3465       aStr.SetValue( i, ',');
3466   }
3467   aStr.RemoveAll( ' ' );
3468
3469   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
3470   int i = 1;
3471   while ( tmpStr != "" )
3472   {
3473     tmpStr = aStr.Token( ",", i++ );
3474     int aPos = tmpStr.Search( '-' );
3475
3476     if ( aPos == -1 )
3477     {
3478       if ( tmpStr.IsIntegerValue() )
3479         myIds.Add( tmpStr.IntegerValue() );
3480       else
3481         return false;
3482     }
3483     else
3484     {
3485       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
3486       TCollection_AsciiString aMinStr = tmpStr;
3487
3488       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
3489       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
3490
3491       if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) ||
3492            (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
3493         return false;
3494
3495       myMin.push_back( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
3496       myMax.push_back( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
3497     }
3498   }
3499
3500   return true;
3501 }
3502
3503 //=======================================================================
3504 // name    : GetType
3505 // Purpose : Get type of supported entities
3506 //=======================================================================
3507 SMDSAbs_ElementType RangeOfIds::GetType() const
3508 {
3509   return myType;
3510 }
3511
3512 //=======================================================================
3513 // name    : SetType
3514 // Purpose : Set type of supported entities
3515 //=======================================================================
3516 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
3517 {
3518   myType = theType;
3519 }
3520
3521 //=======================================================================
3522 // name    : IsSatisfy
3523 // Purpose : Verify whether entity satisfies to this rpedicate
3524 //=======================================================================
3525 bool RangeOfIds::IsSatisfy( long theId )
3526 {
3527   if ( !myMesh )
3528     return false;
3529
3530   if ( myType == SMDSAbs_Node )
3531   {
3532     if ( myMesh->FindNode( theId ) == 0 )
3533       return false;
3534   }
3535   else
3536   {
3537     const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3538     if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All ))
3539       return false;
3540   }
3541
3542   if ( myIds.Contains( theId ) )
3543     return true;
3544
3545   for ( size_t i = 0; i < myMin.size(); i++ )
3546     if ( theId >= myMin[i] && theId <= myMax[i] )
3547       return true;
3548
3549   return false;
3550 }
3551
3552 /*
3553   Class       : Comparator
3554   Description : Base class for comparators
3555 */
3556 Comparator::Comparator():
3557   myMargin(0)
3558 {}
3559
3560 Comparator::~Comparator()
3561 {}
3562
3563 void Comparator::SetMesh( const SMDS_Mesh* theMesh )
3564 {
3565   if ( myFunctor )
3566     myFunctor->SetMesh( theMesh );
3567 }
3568
3569 void Comparator::SetMargin( double theValue )
3570 {
3571   myMargin = theValue;
3572 }
3573
3574 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
3575 {
3576   myFunctor = theFunct;
3577 }
3578
3579 SMDSAbs_ElementType Comparator::GetType() const
3580 {
3581   return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
3582 }
3583
3584 double Comparator::GetMargin()
3585 {
3586   return myMargin;
3587 }
3588
3589
3590 /*
3591   Class       : LessThan
3592   Description : Comparator "<"
3593 */
3594 bool LessThan::IsSatisfy( long theId )
3595 {
3596   return myFunctor && myFunctor->GetValue( theId ) < myMargin;
3597 }
3598
3599
3600 /*
3601   Class       : MoreThan
3602   Description : Comparator ">"
3603 */
3604 bool MoreThan::IsSatisfy( long theId )
3605 {
3606   return myFunctor && myFunctor->GetValue( theId ) > myMargin;
3607 }
3608
3609
3610 /*
3611   Class       : EqualTo
3612   Description : Comparator "="
3613 */
3614 EqualTo::EqualTo():
3615   myToler(Precision::Confusion())
3616 {}
3617
3618 bool EqualTo::IsSatisfy( long theId )
3619 {
3620   return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
3621 }
3622
3623 void EqualTo::SetTolerance( double theToler )
3624 {
3625   myToler = theToler;
3626 }
3627
3628 double EqualTo::GetTolerance()
3629 {
3630   return myToler;
3631 }
3632
3633 /*
3634   Class       : LogicalNOT
3635   Description : Logical NOT predicate
3636 */
3637 LogicalNOT::LogicalNOT()
3638 {}
3639
3640 LogicalNOT::~LogicalNOT()
3641 {}
3642
3643 bool LogicalNOT::IsSatisfy( long theId )
3644 {
3645   return myPredicate && !myPredicate->IsSatisfy( theId );
3646 }
3647
3648 void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
3649 {
3650   if ( myPredicate )
3651     myPredicate->SetMesh( theMesh );
3652 }
3653
3654 void LogicalNOT::SetPredicate( PredicatePtr thePred )
3655 {
3656   myPredicate = thePred;
3657 }
3658
3659 SMDSAbs_ElementType LogicalNOT::GetType() const
3660 {
3661   return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
3662 }
3663
3664
3665 /*
3666   Class       : LogicalBinary
3667   Description : Base class for binary logical predicate
3668 */
3669 LogicalBinary::LogicalBinary()
3670 {}
3671
3672 LogicalBinary::~LogicalBinary()
3673 {}
3674
3675 void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
3676 {
3677   if ( myPredicate1 )
3678     myPredicate1->SetMesh( theMesh );
3679
3680   if ( myPredicate2 )
3681     myPredicate2->SetMesh( theMesh );
3682 }
3683
3684 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
3685 {
3686   myPredicate1 = thePredicate;
3687 }
3688
3689 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
3690 {
3691   myPredicate2 = thePredicate;
3692 }
3693
3694 SMDSAbs_ElementType LogicalBinary::GetType() const
3695 {
3696   if ( !myPredicate1 || !myPredicate2 )
3697     return SMDSAbs_All;
3698
3699   SMDSAbs_ElementType aType1 = myPredicate1->GetType();
3700   SMDSAbs_ElementType aType2 = myPredicate2->GetType();
3701
3702   return aType1 == aType2 ? aType1 : SMDSAbs_All;
3703 }
3704
3705
3706 /*
3707   Class       : LogicalAND
3708   Description : Logical AND
3709 */
3710 bool LogicalAND::IsSatisfy( long theId )
3711 {
3712   return
3713     myPredicate1 &&
3714     myPredicate2 &&
3715     myPredicate1->IsSatisfy( theId ) &&
3716     myPredicate2->IsSatisfy( theId );
3717 }
3718
3719
3720 /*
3721   Class       : LogicalOR
3722   Description : Logical OR
3723 */
3724 bool LogicalOR::IsSatisfy( long theId )
3725 {
3726   return
3727     myPredicate1 &&
3728     myPredicate2 &&
3729     (myPredicate1->IsSatisfy( theId ) ||
3730     myPredicate2->IsSatisfy( theId ));
3731 }
3732
3733
3734 /*
3735                               FILTER
3736 */
3737
3738 // #ifdef WITH_TBB
3739 // #include <tbb/parallel_for.h>
3740 // #include <tbb/enumerable_thread_specific.h>
3741
3742 // namespace Parallel
3743 // {
3744 //   typedef tbb::enumerable_thread_specific< TIdSequence > TIdSeq;
3745
3746 //   struct Predicate
3747 //   {
3748 //     const SMDS_Mesh* myMesh;
3749 //     PredicatePtr     myPredicate;
3750 //     TIdSeq &         myOKIds;
3751 //     Predicate( const SMDS_Mesh* m, PredicatePtr p, TIdSeq & ids ):
3752 //       myMesh(m), myPredicate(p->Duplicate()), myOKIds(ids) {}
3753 //     void operator() ( const tbb::blocked_range<size_t>& r ) const
3754 //     {
3755 //       for ( size_t i = r.begin(); i != r.end(); ++i )
3756 //         if ( myPredicate->IsSatisfy( i ))
3757 //           myOKIds.local().push_back();
3758 //     }
3759 //   }
3760 // }
3761 // #endif
3762
3763 Filter::Filter()
3764 {}
3765
3766 Filter::~Filter()
3767 {}
3768
3769 void Filter::SetPredicate( PredicatePtr thePredicate )
3770 {
3771   myPredicate = thePredicate;
3772 }
3773
3774 void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
3775                             PredicatePtr         thePredicate,
3776                             TIdSequence&         theSequence,
3777                             SMDS_ElemIteratorPtr theElements )
3778 {
3779   theSequence.clear();
3780
3781   if ( !theMesh || !thePredicate )
3782     return;
3783
3784   thePredicate->SetMesh( theMesh );
3785
3786   if ( !theElements )
3787     theElements = theMesh->elementsIterator( thePredicate->GetType() );
3788
3789   if ( theElements ) {
3790     while ( theElements->more() ) {
3791       const SMDS_MeshElement* anElem = theElements->next();
3792       if ( thePredicate->GetType() == SMDSAbs_All ||
3793            thePredicate->GetType() == anElem->GetType() )
3794       {
3795         long anId = anElem->GetID();
3796         if ( thePredicate->IsSatisfy( anId ) )
3797           theSequence.push_back( anId );
3798       }
3799     }
3800   }
3801 }
3802
3803 void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
3804                             Filter::TIdSequence& theSequence,
3805                             SMDS_ElemIteratorPtr theElements )
3806 {
3807   GetElementsId(theMesh,myPredicate,theSequence,theElements);
3808 }
3809
3810 /*
3811                               ManifoldPart
3812 */
3813
3814 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
3815
3816 /*
3817    Internal class Link
3818 */
3819
3820 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
3821                           SMDS_MeshNode* theNode2 )
3822 {
3823   myNode1 = theNode1;
3824   myNode2 = theNode2;
3825 }
3826
3827 ManifoldPart::Link::~Link()
3828 {
3829   myNode1 = 0;
3830   myNode2 = 0;
3831 }
3832
3833 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
3834 {
3835   if ( myNode1 == theLink.myNode1 &&
3836        myNode2 == theLink.myNode2 )
3837     return true;
3838   else if ( myNode1 == theLink.myNode2 &&
3839             myNode2 == theLink.myNode1 )
3840     return true;
3841   else
3842     return false;
3843 }
3844
3845 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
3846 {
3847   if(myNode1 < x.myNode1) return true;
3848   if(myNode1 == x.myNode1)
3849     if(myNode2 < x.myNode2) return true;
3850   return false;
3851 }
3852
3853 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
3854                             const ManifoldPart::Link& theLink2 )
3855 {
3856   return theLink1.IsEqual( theLink2 );
3857 }
3858
3859 ManifoldPart::ManifoldPart()
3860 {
3861   myMesh = 0;
3862   myAngToler = Precision::Angular();
3863   myIsOnlyManifold = true;
3864 }
3865
3866 ManifoldPart::~ManifoldPart()
3867 {
3868   myMesh = 0;
3869 }
3870
3871 void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
3872 {
3873   myMesh = theMesh;
3874   process();
3875 }
3876
3877 SMDSAbs_ElementType ManifoldPart::GetType() const
3878 { return SMDSAbs_Face; }
3879
3880 bool ManifoldPart::IsSatisfy( long theElementId )
3881 {
3882   return myMapIds.Contains( theElementId );
3883 }
3884
3885 void ManifoldPart::SetAngleTolerance( const double theAngToler )
3886 { myAngToler = theAngToler; }
3887
3888 double ManifoldPart::GetAngleTolerance() const
3889 { return myAngToler; }
3890
3891 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
3892 { myIsOnlyManifold = theIsOnly; }
3893
3894 void ManifoldPart::SetStartElem( const long  theStartId )
3895 { myStartElemId = theStartId; }
3896
3897 bool ManifoldPart::process()
3898 {
3899   myMapIds.Clear();
3900   myMapBadGeomIds.Clear();
3901
3902   myAllFacePtr.clear();
3903   myAllFacePtrIntDMap.clear();
3904   if ( !myMesh )
3905     return false;
3906
3907   // collect all faces into own map
3908   SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
3909   for (; anFaceItr->more(); )
3910   {
3911     SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
3912     myAllFacePtr.push_back( aFacePtr );
3913     myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
3914   }
3915
3916   SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
3917   if ( !aStartFace )
3918     return false;
3919
3920   // the map of non manifold links and bad geometry
3921   TMapOfLink aMapOfNonManifold;
3922   TIDsMap    aMapOfTreated;
3923
3924   // begin cycle on faces from start index and run on vector till the end
3925   //  and from begin to start index to cover whole vector
3926   const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
3927   bool isStartTreat = false;
3928   for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
3929   {
3930     if ( fi == aStartIndx )
3931       isStartTreat = true;
3932     // as result next time when fi will be equal to aStartIndx
3933
3934     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
3935     if ( aMapOfTreated.Contains( aFacePtr->GetID()) )
3936       continue;
3937
3938     aMapOfTreated.Add( aFacePtr->GetID() );
3939     TIDsMap aResFaces;
3940     if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
3941                          aMapOfNonManifold, aResFaces ) )
3942       continue;
3943     TIDsMap::Iterator anItr( aResFaces );
3944     for ( ; anItr.More(); anItr.Next() )
3945     {
3946       smIdType aFaceId = anItr.Key();
3947       aMapOfTreated.Add( aFaceId );
3948       myMapIds.Add( aFaceId );
3949     }
3950
3951     if ( fi == int( myAllFacePtr.size() - 1 ))
3952       fi = 0;
3953   } // end run on vector of faces
3954   return !myMapIds.IsEmpty();
3955 }
3956
3957 static void getLinks( const SMDS_MeshFace* theFace,
3958                       ManifoldPart::TVectorOfLink& theLinks )
3959 {
3960   int aNbNode = theFace->NbNodes();
3961   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
3962   int i = 1;
3963   SMDS_MeshNode* aNode = 0;
3964   for ( ; aNodeItr->more() && i <= aNbNode; )
3965   {
3966
3967     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
3968     if ( i == 1 )
3969       aNode = aN1;
3970     i++;
3971     SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
3972     i++;
3973     ManifoldPart::Link aLink( aN1, aN2 );
3974     theLinks.push_back( aLink );
3975   }
3976 }
3977
3978 bool ManifoldPart::findConnected
3979                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
3980                   SMDS_MeshFace*                           theStartFace,
3981                   ManifoldPart::TMapOfLink&                theNonManifold,
3982                   TIDsMap&                                 theResFaces )
3983 {
3984   theResFaces.Clear();
3985   if ( !theAllFacePtrInt.size() )
3986     return false;
3987
3988   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
3989   {
3990     myMapBadGeomIds.Add( theStartFace->GetID() );
3991     return false;
3992   }
3993
3994   ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
3995   ManifoldPart::TVectorOfLink aSeqOfBoundary;
3996   theResFaces.Add( theStartFace->GetID() );
3997   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
3998
3999   expandBoundary( aMapOfBoundary, aSeqOfBoundary,
4000                  aDMapLinkFace, theNonManifold, theStartFace );
4001
4002   bool isDone = false;
4003   while ( !isDone && aMapOfBoundary.size() != 0 )
4004   {
4005     bool isToReset = false;
4006     ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
4007     for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
4008     {
4009       ManifoldPart::Link aLink = *pLink;
4010       if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
4011         continue;
4012       // each link could be treated only once
4013       aMapToSkip.insert( aLink );
4014
4015       ManifoldPart::TVectorOfFacePtr aFaces;
4016       // find next
4017       if ( myIsOnlyManifold &&
4018            (theNonManifold.find( aLink ) != theNonManifold.end()) )
4019         continue;
4020       else
4021       {
4022         getFacesByLink( aLink, aFaces );
4023         // filter the element to keep only indicated elements
4024         ManifoldPart::TVectorOfFacePtr aFiltered;
4025         ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
4026         for ( ; pFace != aFaces.end(); ++pFace )
4027         {
4028           SMDS_MeshFace* aFace = *pFace;
4029           if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
4030             aFiltered.push_back( aFace );
4031         }
4032         aFaces = aFiltered;
4033         if ( aFaces.size() < 2 )  // no neihgbour faces
4034           continue;
4035         else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
4036         {
4037           theNonManifold.insert( aLink );
4038           continue;
4039         }
4040       }
4041
4042       // compare normal with normals of neighbor element
4043       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
4044       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
4045       for ( ; pFace != aFaces.end(); ++pFace )
4046       {
4047         SMDS_MeshFace* aNextFace = *pFace;
4048         if ( aPrevFace == aNextFace )
4049           continue;
4050         smIdType anNextFaceID = aNextFace->GetID();
4051         if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
4052          // should not be with non manifold restriction. probably bad topology
4053           continue;
4054         // check if face was treated and skipped
4055         if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
4056              !isInPlane( aPrevFace, aNextFace ) )
4057           continue;
4058         // add new element to connected and extend the boundaries.
4059         theResFaces.Add( anNextFaceID );
4060         expandBoundary( aMapOfBoundary, aSeqOfBoundary,
4061                         aDMapLinkFace, theNonManifold, aNextFace );
4062         isToReset = true;
4063       }
4064     }
4065     isDone = !isToReset;
4066   }
4067
4068   return !theResFaces.IsEmpty();
4069 }
4070
4071 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
4072                               const SMDS_MeshFace* theFace2 )
4073 {
4074   gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
4075   gp_XYZ aNorm2XYZ = getNormale( theFace2 );
4076   if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
4077   {
4078     myMapBadGeomIds.Add( theFace2->GetID() );
4079     return false;
4080   }
4081   if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
4082     return true;
4083
4084   return false;
4085 }
4086
4087 void ManifoldPart::expandBoundary
4088                    ( ManifoldPart::TMapOfLink&            theMapOfBoundary,
4089                      ManifoldPart::TVectorOfLink&         theSeqOfBoundary,
4090                      ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
4091                      ManifoldPart::TMapOfLink&            theNonManifold,
4092                      SMDS_MeshFace*                       theNextFace ) const
4093 {
4094   ManifoldPart::TVectorOfLink aLinks;
4095   getLinks( theNextFace, aLinks );
4096   int aNbLink = (int)aLinks.size();
4097   for ( int i = 0; i < aNbLink; i++ )
4098   {
4099     ManifoldPart::Link aLink = aLinks[ i ];
4100     if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
4101       continue;
4102     if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
4103     {
4104       if ( myIsOnlyManifold )
4105       {
4106         // remove from boundary
4107         theMapOfBoundary.erase( aLink );
4108         ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
4109         for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
4110         {
4111           ManifoldPart::Link aBoundLink = *pLink;
4112           if ( aBoundLink.IsEqual( aLink ) )
4113           {
4114             theSeqOfBoundary.erase( pLink );
4115             break;
4116           }
4117         }
4118       }
4119     }
4120     else
4121     {
4122       theMapOfBoundary.insert( aLink );
4123       theSeqOfBoundary.push_back( aLink );
4124       theDMapLinkFacePtr[ aLink ] = theNextFace;
4125     }
4126   }
4127 }
4128
4129 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
4130                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
4131 {
4132
4133   // take all faces that shared first node
4134   SMDS_ElemIteratorPtr anItr = theLink.myNode1->GetInverseElementIterator( SMDSAbs_Face );
4135   SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > faces( anItr ), facesEnd;
4136   std::set<const SMDS_MeshElement *> aSetOfFaces( faces, facesEnd );
4137
4138   // take all faces that shared second node
4139   anItr = theLink.myNode2->GetInverseElementIterator( SMDSAbs_Face );
4140   // find the common part of two sets
4141   for ( ; anItr->more(); )
4142   {
4143     const SMDS_MeshElement* aFace = anItr->next();
4144     if ( aSetOfFaces.count( aFace ))
4145       theFaces.push_back( (SMDS_MeshFace*) aFace );
4146   }
4147 }
4148
4149 /*
4150   Class       : BelongToMeshGroup
4151   Description : Verify whether a mesh element is included into a mesh group
4152 */
4153 BelongToMeshGroup::BelongToMeshGroup(): myGroup( 0 )
4154 {
4155 }
4156
4157 void BelongToMeshGroup::SetGroup( SMESHDS_GroupBase* g )
4158 {
4159   myGroup = g;
4160 }
4161
4162 void BelongToMeshGroup::SetStoreName( const std::string& sn )
4163 {
4164   myStoreName = sn;
4165 }
4166
4167 void BelongToMeshGroup::SetMesh( const SMDS_Mesh* theMesh )
4168 {
4169   if ( myGroup && myGroup->GetMesh() != theMesh )
4170   {
4171     myGroup = 0;
4172   }
4173   if ( !myGroup && !myStoreName.empty() )
4174   {
4175     if ( const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh))
4176     {
4177       const std::set<SMESHDS_GroupBase*>& grps = aMesh->GetGroups();
4178       std::set<SMESHDS_GroupBase*>::const_iterator g = grps.begin();
4179       for ( ; g != grps.end() && !myGroup; ++g )
4180         if ( *g && myStoreName == (*g)->GetStoreName() )
4181           myGroup = *g;
4182     }
4183   }
4184   if ( myGroup )
4185   {
4186     myGroup->IsEmpty(); // make GroupOnFilter update its predicate
4187   }
4188 }
4189
4190 bool BelongToMeshGroup::IsSatisfy( long theElementId )
4191 {
4192   return myGroup ? myGroup->Contains( theElementId ) : false;
4193 }
4194
4195 SMDSAbs_ElementType BelongToMeshGroup::GetType() const
4196 {
4197   return myGroup ? myGroup->GetType() : SMDSAbs_All;
4198 }
4199
4200 //================================================================================
4201 //  ElementsOnSurface
4202 //================================================================================
4203
4204 ElementsOnSurface::ElementsOnSurface()
4205 {
4206   myIds.Clear();
4207   myType = SMDSAbs_All;
4208   mySurf.Nullify();
4209   myToler = Precision::Confusion();
4210   myUseBoundaries = false;
4211 }
4212
4213 ElementsOnSurface::~ElementsOnSurface()
4214 {
4215 }
4216
4217 void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
4218 {
4219   myMeshModifTracer.SetMesh( theMesh );
4220   if ( myMeshModifTracer.IsMeshModified())
4221     process();
4222 }
4223
4224 bool ElementsOnSurface::IsSatisfy( long theElementId )
4225 {
4226   return myIds.Contains( theElementId );
4227 }
4228
4229 SMDSAbs_ElementType ElementsOnSurface::GetType() const
4230 { return myType; }
4231
4232 void ElementsOnSurface::SetTolerance( const double theToler )
4233 {
4234   if ( myToler != theToler )
4235   {
4236     myToler = theToler;
4237     process();
4238   }
4239 }
4240
4241 double ElementsOnSurface::GetTolerance() const
4242 { return myToler; }
4243
4244 void ElementsOnSurface::SetUseBoundaries( bool theUse )
4245 {
4246   if ( myUseBoundaries != theUse ) {
4247     myUseBoundaries = theUse;
4248     SetSurface( mySurf, myType );
4249   }
4250 }
4251
4252 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
4253                                     const SMDSAbs_ElementType theType )
4254 {
4255   myIds.Clear();
4256   myType = theType;
4257   mySurf.Nullify();
4258   if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
4259     return;
4260   mySurf = TopoDS::Face( theShape );
4261   BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
4262   Standard_Real
4263     u1 = SA.FirstUParameter(),
4264     u2 = SA.LastUParameter(),
4265     v1 = SA.FirstVParameter(),
4266     v2 = SA.LastVParameter();
4267   Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
4268   myProjector.Init( surf, u1,u2, v1,v2 );
4269   process();
4270 }
4271
4272 void ElementsOnSurface::process()
4273 {
4274   myIds.Clear();
4275   if ( mySurf.IsNull() )
4276     return;
4277
4278   if ( !myMeshModifTracer.GetMesh() )
4279     return;
4280
4281   int nbElems = FromSmIdType<int>( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
4282   if ( nbElems > 0 )
4283     myIds.ReSize( nbElems );
4284
4285   SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType );
4286   for(; anIter->more(); )
4287     process( anIter->next() );
4288 }
4289
4290 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
4291 {
4292   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
4293   bool isSatisfy = true;
4294   for ( ; aNodeItr->more(); )
4295   {
4296     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
4297     if ( !isOnSurface( aNode ) )
4298     {
4299       isSatisfy = false;
4300       break;
4301     }
4302   }
4303   if ( isSatisfy )
4304     myIds.Add( theElemPtr->GetID() );
4305 }
4306
4307 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
4308 {
4309   if ( mySurf.IsNull() )
4310     return false;
4311
4312   gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
4313   //  double aToler2 = myToler * myToler;
4314 //   if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
4315 //   {
4316 //     gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
4317 //     if ( aPln.SquareDistance( aPnt ) > aToler2 )
4318 //       return false;
4319 //   }
4320 //   else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
4321 //   {
4322 //     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
4323 //     double aRad = aCyl.Radius();
4324 //     gp_Ax3 anAxis = aCyl.Position();
4325 //     gp_XYZ aLoc = aCyl.Location().XYZ();
4326 //     double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4327 //     double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4328 //     if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
4329 //       return false;
4330 //   }
4331 //   else
4332 //     return false;
4333   myProjector.Perform( aPnt );
4334   bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
4335
4336   return isOn;
4337 }
4338
4339
4340 //================================================================================
4341 //  ElementsOnShape
4342 //================================================================================
4343
4344 namespace {
4345   const int theIsCheckedFlag = 0x0000100;
4346 }
4347
4348 struct ElementsOnShape::Classifier
4349 {
4350   Classifier() { mySolidClfr = 0; myFlags = 0; }
4351   ~Classifier();
4352   void Init(const TopoDS_Shape& s, double tol, const Bnd_B3d* box = 0 );
4353   bool IsOut(const gp_Pnt& p)        { return SetChecked( true ), (this->*myIsOutFun)( p ); }
4354   TopAbs_ShapeEnum ShapeType() const { return myShape.ShapeType(); }
4355   const TopoDS_Shape& Shape() const  { return myShape; }
4356   const Bnd_B3d* GetBndBox() const   { return & myBox; }
4357   double Tolerance() const           { return myTol; }
4358   bool IsChecked()                   { return myFlags & theIsCheckedFlag; }
4359   bool IsSetFlag( int flag ) const   { return myFlags & flag; }
4360   void SetChecked( bool is ) { is ? SetFlag( theIsCheckedFlag ) : UnsetFlag( theIsCheckedFlag ); }
4361   void SetFlag  ( int flag ) { myFlags |= flag; }
4362   void UnsetFlag( int flag ) { myFlags &= ~flag; }
4363
4364 private:
4365   bool isOutOfSolid (const gp_Pnt& p);
4366   bool isOutOfBox   (const gp_Pnt& p);
4367   bool isOutOfFace  (const gp_Pnt& p);
4368   bool isOutOfEdge  (const gp_Pnt& p);
4369   bool isOutOfVertex(const gp_Pnt& p);
4370   bool isOutOfNone  (const gp_Pnt& /*p*/) { return true; }
4371   bool isBox        (const TopoDS_Shape& s);
4372
4373   TopoDS_Shape prepareSolid( const TopoDS_Shape& theSolid );
4374
4375   bool (Classifier::*          myIsOutFun)(const gp_Pnt& p);
4376   BRepClass3d_SolidClassifier* mySolidClfr; // ptr because of a run-time forbidden copy-constructor
4377   Bnd_B3d                      myBox;
4378   GeomAPI_ProjectPointOnSurf   myProjFace;
4379   GeomAPI_ProjectPointOnCurve  myProjEdge;
4380   gp_Pnt                       myVertexXYZ;
4381   TopoDS_Shape                 myShape;
4382   double                       myTol;
4383   int                          myFlags;
4384 };
4385
4386 struct ElementsOnShape::OctreeClassifier : public SMESH_Octree
4387 {
4388   OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers );
4389   OctreeClassifier( const OctreeClassifier*                           otherTree,
4390                     const std::vector< ElementsOnShape::Classifier >& clsOther,
4391                     std::vector< ElementsOnShape::Classifier >&       cls );
4392   void GetClassifiersAtPoint( const gp_XYZ& p,
4393                               std::vector< ElementsOnShape::Classifier* >& classifiers );
4394   size_t GetSize();
4395
4396 protected:
4397   OctreeClassifier() {}
4398   SMESH_Octree* newChild() const { return new OctreeClassifier; }
4399   void          buildChildrenData();
4400   Bnd_B3d*      buildRootBox();
4401
4402   std::vector< ElementsOnShape::Classifier* > myClassifiers;
4403 };
4404
4405
4406 ElementsOnShape::ElementsOnShape():
4407   myOctree(0),
4408   myType(SMDSAbs_All),
4409   myToler(Precision::Confusion()),
4410   myAllNodesFlag(false)
4411 {
4412 }
4413
4414 ElementsOnShape::~ElementsOnShape()
4415 {
4416   clearClassifiers();
4417 }
4418
4419 Predicate* ElementsOnShape::clone() const
4420 {
4421   size_t size = sizeof( *this );
4422   if ( myOctree )
4423     size += myOctree->GetSize();
4424   if ( !myClassifiers.empty() )
4425     size += sizeof( myClassifiers[0] ) * myClassifiers.size();
4426   if ( !myWorkClassifiers.empty() )
4427     size += sizeof( myWorkClassifiers[0] ) * myWorkClassifiers.size();
4428   if ( size > 1e+9 ) // 1G
4429   {
4430 #ifdef _DEBUG_
4431     std::cout << "Avoid ElementsOnShape::clone(), too large: " << size << " bytes " << std::endl;
4432 #endif
4433     return 0;
4434   }
4435
4436   ElementsOnShape* cln = new ElementsOnShape();
4437   cln->SetAllNodes ( myAllNodesFlag );
4438   cln->SetTolerance( myToler );
4439   cln->SetMesh     ( myMeshModifTracer.GetMesh() );
4440   cln->myShape = myShape; // avoid creation of myClassifiers
4441   cln->SetShape    ( myShape, myType );
4442   cln->myClassifiers.resize( myClassifiers.size() );
4443   for ( size_t i = 0; i < myClassifiers.size(); ++i )
4444     cln->myClassifiers[ i ].Init( BRepBuilderAPI_Copy( myClassifiers[ i ].Shape()),
4445                                   myToler, myClassifiers[ i ].GetBndBox() );
4446   if ( myOctree ) // copy myOctree
4447   {
4448     cln->myOctree = new OctreeClassifier( myOctree, myClassifiers, cln->myClassifiers );
4449   }
4450   return cln;
4451 }
4452
4453 SMDSAbs_ElementType ElementsOnShape::GetType() const
4454 {
4455   return myType;
4456 }
4457
4458 void ElementsOnShape::SetTolerance (const double theToler)
4459 {
4460   if (myToler != theToler) {
4461     myToler = theToler;
4462     SetShape(myShape, myType);
4463   }
4464 }
4465
4466 double ElementsOnShape::GetTolerance() const
4467 {
4468   return myToler;
4469 }
4470
4471 void ElementsOnShape::SetAllNodes (bool theAllNodes)
4472 {
4473   myAllNodesFlag = theAllNodes;
4474 }
4475
4476 void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
4477 {
4478   myMeshModifTracer.SetMesh( theMesh );
4479   if ( myMeshModifTracer.IsMeshModified())
4480   {
4481     size_t nbNodes = theMesh ? theMesh->NbNodes() : 0;
4482     if ( myNodeIsChecked.size() == nbNodes )
4483     {
4484       std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4485     }
4486     else
4487     {
4488       SMESHUtils::FreeVector( myNodeIsChecked );
4489       SMESHUtils::FreeVector( myNodeIsOut );
4490       myNodeIsChecked.resize( nbNodes, false );
4491       myNodeIsOut.resize( nbNodes );
4492     }
4493   }
4494 }
4495
4496 bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut )
4497 {
4498   if ( n->GetID() >= (int) myNodeIsChecked.size() ||
4499        !myNodeIsChecked[ n->GetID() ])
4500     return false;
4501
4502   isOut = myNodeIsOut[ n->GetID() ];
4503   return true;
4504 }
4505
4506 void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool  isOut )
4507 {
4508   if ( n->GetID() < (int) myNodeIsChecked.size() )
4509   {
4510     myNodeIsChecked[ n->GetID() ] = true;
4511     myNodeIsOut    [ n->GetID() ] = isOut;
4512   }
4513 }
4514
4515 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
4516                                 const SMDSAbs_ElementType theType)
4517 {
4518   bool shapeChanges = ( myShape != theShape );
4519   myType  = theType;
4520   myShape = theShape;
4521   if ( myShape.IsNull() ) return;
4522
4523   if ( shapeChanges )
4524   {
4525     // find most complex shapes
4526     TopTools_IndexedMapOfShape shapesMap;
4527     TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
4528     TopExp_Explorer sub;
4529     for ( int i = 0; i < 4; ++i )
4530     {
4531       if ( shapesMap.IsEmpty() )
4532         for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() )
4533           shapesMap.Add( sub.Current() );
4534       if ( i > 0 )
4535         for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() )
4536           shapesMap.Add( sub.Current() );
4537     }
4538
4539     clearClassifiers();
4540     myClassifiers.resize( shapesMap.Extent() );
4541     for ( int i = 0; i < shapesMap.Extent(); ++i )
4542       myClassifiers[ i ].Init( shapesMap( i+1 ), myToler );
4543   }
4544
4545   if ( theType == SMDSAbs_Node )
4546   {
4547     SMESHUtils::FreeVector( myNodeIsChecked );
4548     SMESHUtils::FreeVector( myNodeIsOut );
4549   }
4550   else
4551   {
4552     std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4553   }
4554 }
4555
4556 void ElementsOnShape::clearClassifiers()
4557 {
4558   // for ( size_t i = 0; i < myClassifiers.size(); ++i )
4559   //   delete myClassifiers[ i ];
4560   myClassifiers.clear();
4561
4562   delete myOctree;
4563   myOctree = 0;
4564 }
4565
4566 bool ElementsOnShape::IsSatisfy( long elemId )
4567 {
4568   if ( myClassifiers.empty() )
4569     return false;
4570
4571   const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh();
4572   if ( myType == SMDSAbs_Node )
4573     return IsSatisfy( mesh->FindNode( elemId ));
4574   return IsSatisfy( mesh->FindElement( elemId ));
4575 }
4576
4577 bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem)
4578 {
4579   if ( !elem )
4580     return false;
4581
4582   bool isSatisfy = myAllNodesFlag, isNodeOut;
4583
4584   gp_XYZ centerXYZ (0, 0, 0);
4585
4586   if ( !myOctree && myClassifiers.size() > 5 )
4587   {
4588     myWorkClassifiers.resize( myClassifiers.size() );
4589     for ( size_t i = 0; i < myClassifiers.size(); ++i )
4590       myWorkClassifiers[ i ] = & myClassifiers[ i ];
4591     myOctree = new OctreeClassifier( myWorkClassifiers );
4592
4593     SMESHUtils::FreeVector( myWorkClassifiers );
4594   }
4595
4596   for ( int i = 0, nb = elem->NbNodes(); i < nb  && (isSatisfy == myAllNodesFlag); ++i )
4597   {
4598     SMESH_TNodeXYZ aPnt( elem->GetNode( i ));
4599     centerXYZ += aPnt;
4600
4601     isNodeOut = true;
4602     if ( !getNodeIsOut( aPnt._node, isNodeOut ))
4603     {
4604       if ( myOctree )
4605       {
4606         myWorkClassifiers.clear();
4607         myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers );
4608
4609         for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
4610           myWorkClassifiers[i]->SetChecked( false );
4611
4612         for ( size_t i = 0; i < myWorkClassifiers.size() && isNodeOut; ++i )
4613           if ( !myWorkClassifiers[i]->IsChecked() )
4614             isNodeOut = myWorkClassifiers[i]->IsOut( aPnt );
4615       }
4616       else
4617       {
4618         for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i )
4619           isNodeOut = myClassifiers[i].IsOut( aPnt );
4620       }
4621       setNodeIsOut( aPnt._node, isNodeOut );
4622     }
4623     isSatisfy = !isNodeOut;
4624   }
4625
4626   // Check the center point for volumes MantisBug 0020168
4627   if ( isSatisfy &&
4628        myAllNodesFlag &&
4629        myClassifiers[0].ShapeType() == TopAbs_SOLID )
4630   {
4631     centerXYZ /= elem->NbNodes();
4632     isSatisfy = false;
4633     if ( myOctree )
4634     {
4635       myWorkClassifiers.clear();
4636       myOctree->GetClassifiersAtPoint( centerXYZ, myWorkClassifiers );
4637       for ( size_t i = 0; i < myWorkClassifiers.size() && !isSatisfy; ++i )
4638         isSatisfy = ! myWorkClassifiers[i]->IsOut( centerXYZ );
4639     }
4640     else
4641     {
4642       for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i )
4643         isSatisfy = ! myClassifiers[i].IsOut( centerXYZ );
4644     }
4645   }
4646
4647   return isSatisfy;
4648 }
4649
4650 //================================================================================
4651 /*!
4652  * \brief Check and optionally return a satisfying shape
4653  */
4654 //================================================================================
4655
4656 bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node,
4657                                  TopoDS_Shape*        okShape)
4658 {
4659   if ( !node )
4660     return false;
4661
4662   if ( !myOctree && myClassifiers.size() > 5 )
4663   {
4664     myWorkClassifiers.resize( myClassifiers.size() );
4665     for ( size_t i = 0; i < myClassifiers.size(); ++i )
4666       myWorkClassifiers[ i ] = & myClassifiers[ i ];
4667     myOctree = new OctreeClassifier( myWorkClassifiers );
4668   }
4669
4670   bool isNodeOut = true;
4671
4672   if ( okShape || !getNodeIsOut( node, isNodeOut ))
4673   {
4674     SMESH_NodeXYZ aPnt = node;
4675     if ( myOctree )
4676     {
4677       myWorkClassifiers.clear();
4678       myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers );
4679
4680       for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
4681         myWorkClassifiers[i]->SetChecked( false );
4682
4683       for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
4684         if ( !myWorkClassifiers[i]->IsChecked() &&
4685              !myWorkClassifiers[i]->IsOut( aPnt ))
4686         {
4687           isNodeOut = false;
4688           if ( okShape )
4689             *okShape = myWorkClassifiers[i]->Shape();
4690           break;
4691         }
4692     }
4693     else
4694     {
4695       for ( size_t i = 0; i < myClassifiers.size(); ++i )
4696         if ( !myClassifiers[i].IsOut( aPnt ))
4697         {
4698           isNodeOut = false;
4699           if ( okShape )
4700             *okShape = myClassifiers[i].Shape();
4701           break;
4702         }
4703     }
4704     setNodeIsOut( node, isNodeOut );
4705   }
4706
4707   return !isNodeOut;
4708 }
4709
4710 void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape,
4711                                         double              theTol,
4712                                         const Bnd_B3d*      theBox )
4713 {
4714   myShape = theShape;
4715   myTol   = theTol;
4716   myFlags = 0;
4717
4718   bool isShapeBox = false;
4719   switch ( myShape.ShapeType() )
4720   {
4721   case TopAbs_SOLID:
4722   {
4723     if (( isShapeBox = isBox( theShape )))
4724     {
4725       myIsOutFun = & ElementsOnShape::Classifier::isOutOfBox;
4726     }
4727     else
4728     {
4729       mySolidClfr = new BRepClass3d_SolidClassifier( prepareSolid( theShape ));
4730       myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid;
4731     }
4732     break;
4733   }
4734   case TopAbs_FACE:
4735   {
4736     Standard_Real u1,u2,v1,v2;
4737     Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
4738     if ( surf.IsNull() )
4739       myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
4740     else
4741     {
4742       surf->Bounds( u1,u2,v1,v2 );
4743       myProjFace.Init(surf, u1,u2, v1,v2, myTol );
4744       myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace;
4745     }
4746     break;
4747   }
4748   case TopAbs_EDGE:
4749   {
4750     Standard_Real u1, u2;
4751     Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2);
4752     if ( curve.IsNull() )
4753       myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
4754     else
4755     {
4756       myProjEdge.Init(curve, u1, u2);
4757       myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
4758     }
4759     break;
4760   }
4761   case TopAbs_VERTEX:
4762   {
4763     myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
4764     myIsOutFun = & ElementsOnShape::Classifier::isOutOfVertex;
4765     break;
4766   }
4767   default:
4768     throw SALOME_Exception("Programmer error in usage of ElementsOnShape::Classifier");
4769   }
4770
4771   if ( !isShapeBox )
4772   {
4773     if ( theBox )
4774     {
4775       myBox = *theBox;
4776     }
4777     else
4778     {
4779       Bnd_Box box;
4780       if ( myShape.ShapeType() == TopAbs_FACE )
4781       {
4782         BRepAdaptor_Surface SA( TopoDS::Face( myShape ), /*useBoundaries=*/false );
4783         if ( SA.GetType() == GeomAbs_BSplineSurface )
4784           BRepBndLib::AddOptimal( myShape, box,
4785                                   /*useTriangulation=*/true, /*useShapeTolerance=*/true );
4786       }
4787       if ( box.IsVoid() )
4788         BRepBndLib::Add( myShape, box );
4789       myBox.Clear();
4790       myBox.Add( box.CornerMin() );
4791       myBox.Add( box.CornerMax() );
4792       gp_XYZ halfSize = 0.5 * ( box.CornerMax().XYZ() - box.CornerMin().XYZ() );
4793       for ( int iDim = 1; iDim <= 3; ++iDim )
4794       {
4795         double x = halfSize.Coord( iDim );
4796         halfSize.SetCoord( iDim, x + Max( myTol, 1e-2 * x ));
4797       }
4798       myBox.SetHSize( halfSize );
4799     }
4800   }
4801 }
4802
4803 ElementsOnShape::Classifier::~Classifier()
4804 {
4805   delete mySolidClfr; mySolidClfr = 0;
4806 }
4807
4808 TopoDS_Shape ElementsOnShape::Classifier::prepareSolid( const TopoDS_Shape& theSolid )
4809 {
4810   // try to limit tolerance of theSolid down to myTol (issue #19026)
4811
4812   // check if tolerance of theSolid is more than myTol
4813   bool tolIsOk = true; // max tolerance is at VERTEXes
4814   for ( TopExp_Explorer exp( theSolid, TopAbs_VERTEX ); exp.More() &&  tolIsOk; exp.Next() )
4815     tolIsOk = ( myTol >= BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current() )));
4816   if ( tolIsOk )
4817     return theSolid;
4818
4819   // make a copy to prevent the original shape from changes
4820   TopoDS_Shape resultShape = BRepBuilderAPI_Copy( theSolid );
4821
4822   if ( !GEOMUtils::FixShapeTolerance( resultShape, TopAbs_SHAPE, myTol ))
4823     return theSolid;
4824   return resultShape;
4825 }
4826
4827 bool ElementsOnShape::Classifier::isOutOfSolid( const gp_Pnt& p )
4828 {
4829   if ( isOutOfBox( p )) return true;
4830   mySolidClfr->Perform( p, myTol );
4831   return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON );
4832 }
4833
4834 bool ElementsOnShape::Classifier::isOutOfBox( const gp_Pnt& p )
4835 {
4836   return myBox.IsOut( p.XYZ() );
4837 }
4838
4839 bool ElementsOnShape::Classifier::isOutOfFace( const gp_Pnt& p )
4840 {
4841   if ( isOutOfBox( p )) return true;
4842   myProjFace.Perform( p );
4843   if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
4844   {
4845     // check relatively to the face
4846     Standard_Real u, v;
4847     myProjFace.LowerDistanceParameters(u, v);
4848     gp_Pnt2d aProjPnt (u, v);
4849     BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
4850     if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
4851       return false;
4852   }
4853   return true;
4854 }
4855
4856 bool ElementsOnShape::Classifier::isOutOfEdge( const gp_Pnt& p )
4857 {
4858   if ( isOutOfBox( p )) return true;
4859   myProjEdge.Perform( p );
4860   return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
4861 }
4862
4863 bool ElementsOnShape::Classifier::isOutOfVertex( const gp_Pnt& p )
4864 {
4865   return ( myVertexXYZ.Distance( p ) > myTol );
4866 }
4867
4868 bool ElementsOnShape::Classifier::isBox(const TopoDS_Shape& theShape )
4869 {
4870   TopTools_IndexedMapOfShape vMap;
4871   TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
4872   if ( vMap.Extent() != 8 )
4873     return false;
4874
4875   myBox.Clear();
4876   for ( int i = 1; i <= 8; ++i )
4877     myBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))).XYZ() );
4878
4879   gp_XYZ pMin = myBox.CornerMin(), pMax = myBox.CornerMax();
4880   for ( int i = 1; i <= 8; ++i )
4881   {
4882     gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i )));
4883     for ( int iC = 1; iC <= 3; ++ iC )
4884     {
4885       double d1 = Abs( pMin.Coord( iC ) - p.Coord( iC ));
4886       double d2 = Abs( pMax.Coord( iC ) - p.Coord( iC ));
4887       if ( Min( d1, d2 ) > myTol )
4888         return false;
4889     }
4890   }
4891   myBox.Enlarge( myTol );
4892   return true;
4893 }
4894
4895 ElementsOnShape::
4896 OctreeClassifier::OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers )
4897   :SMESH_Octree( new SMESH_TreeLimit )
4898 {
4899   myClassifiers = classifiers;
4900   compute();
4901 }
4902
4903 ElementsOnShape::
4904 OctreeClassifier::OctreeClassifier( const OctreeClassifier*                           otherTree,
4905                                     const std::vector< ElementsOnShape::Classifier >& clsOther,
4906                                     std::vector< ElementsOnShape::Classifier >&       cls )
4907   :SMESH_Octree( new SMESH_TreeLimit )
4908 {
4909   myBox = new Bnd_B3d( *otherTree->getBox() );
4910
4911   if (( myIsLeaf = otherTree->isLeaf() ))
4912   {
4913     myClassifiers.resize( otherTree->myClassifiers.size() );
4914     for ( size_t i = 0; i < otherTree->myClassifiers.size(); ++i )
4915     {
4916       int ind = otherTree->myClassifiers[i] - & clsOther[0];
4917       myClassifiers[ i ] = & cls[ ind ];
4918     }
4919   }
4920   else if ( otherTree->myChildren )
4921   {
4922     myChildren = new SMESH_Tree< Bnd_B3d, 8 > * [ 8 ];
4923     for ( int i = 0; i < nbChildren(); i++ )
4924       myChildren[i] =
4925         new OctreeClassifier( static_cast<const OctreeClassifier*>( otherTree->myChildren[i]),
4926                               clsOther, cls );
4927   }
4928 }
4929
4930 void ElementsOnShape::
4931 OctreeClassifier::GetClassifiersAtPoint( const gp_XYZ& point,
4932                                          std::vector< ElementsOnShape::Classifier* >& result )
4933 {
4934   if ( getBox()->IsOut( point ))
4935     return;
4936
4937   if ( isLeaf() )
4938   {
4939     for ( size_t i = 0; i < myClassifiers.size(); ++i )
4940       if ( !myClassifiers[i]->GetBndBox()->IsOut( point ))
4941         result.push_back( myClassifiers[i] );
4942   }
4943   else
4944   {
4945     for (int i = 0; i < nbChildren(); i++)
4946       ((OctreeClassifier*) myChildren[i])->GetClassifiersAtPoint( point, result );
4947   }
4948 }
4949
4950 size_t ElementsOnShape::OctreeClassifier::GetSize()
4951 {
4952   size_t res = sizeof( *this );
4953   if ( !myClassifiers.empty() )
4954     res += sizeof( myClassifiers[0] ) * myClassifiers.size();
4955
4956   if ( !isLeaf() )
4957     for (int i = 0; i < nbChildren(); i++)
4958       res += ((OctreeClassifier*) myChildren[i])->GetSize();
4959
4960   return res;
4961 }
4962
4963 void ElementsOnShape::OctreeClassifier::buildChildrenData()
4964 {
4965   // distribute myClassifiers among myChildren
4966
4967   const int childFlag[8] = { 0x0000001,
4968                              0x0000002,
4969                              0x0000004,
4970                              0x0000008,
4971                              0x0000010,
4972                              0x0000020,
4973                              0x0000040,
4974                              0x0000080 };
4975   int nbInChild[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
4976
4977   for ( size_t i = 0; i < myClassifiers.size(); ++i )
4978   {
4979     for ( int j = 0; j < nbChildren(); j++ )
4980     {
4981       if ( !myClassifiers[i]->GetBndBox()->IsOut( *myChildren[j]->getBox() ))
4982       {
4983         myClassifiers[i]->SetFlag( childFlag[ j ]);
4984         ++nbInChild[ j ];
4985       }
4986     }
4987   }
4988
4989   for ( int j = 0; j < nbChildren(); j++ )
4990   {
4991     OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ j ]);
4992     child->myClassifiers.resize( nbInChild[ j ]);
4993     for ( size_t i = 0; nbInChild[ j ] && i < myClassifiers.size(); ++i )
4994     {
4995       if ( myClassifiers[ i ]->IsSetFlag( childFlag[ j ]))
4996       {
4997         --nbInChild[ j ];
4998         child->myClassifiers[ nbInChild[ j ]] = myClassifiers[ i ];
4999         myClassifiers[ i ]->UnsetFlag( childFlag[ j ]);
5000       }
5001     }
5002   }
5003   SMESHUtils::FreeVector( myClassifiers );
5004
5005   // define if a child isLeaf()
5006   for ( int i = 0; i < nbChildren(); i++ )
5007   {
5008     OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ i ]);
5009     child->myIsLeaf = ( child->myClassifiers.size() <= 5  ||
5010                         child->maxSize() < child->myClassifiers[0]->Tolerance() );
5011   }
5012 }
5013
5014 Bnd_B3d* ElementsOnShape::OctreeClassifier::buildRootBox()
5015 {
5016   Bnd_B3d* box = new Bnd_B3d;
5017   for ( size_t i = 0; i < myClassifiers.size(); ++i )
5018     box->Add( *myClassifiers[i]->GetBndBox() );
5019   return box;
5020 }
5021
5022 /*
5023   Class       : BelongToGeom
5024   Description : Predicate for verifying whether entity belongs to
5025                 specified geometrical support
5026 */
5027
5028 BelongToGeom::BelongToGeom()
5029   : myMeshDS(NULL),
5030     myType(SMDSAbs_NbElementTypes),
5031     myIsSubshape(false),
5032     myTolerance(Precision::Confusion())
5033 {}
5034
5035 Predicate* BelongToGeom::clone() const
5036 {
5037   BelongToGeom* cln = 0;
5038   if ( myElementsOnShapePtr )
5039     if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
5040     {
5041       cln = new BelongToGeom( *this );
5042       cln->myElementsOnShapePtr.reset( eos );
5043     }
5044   return cln;
5045 }
5046
5047 void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
5048 {
5049   if ( myMeshDS != theMesh )
5050   {
5051     myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
5052     init();
5053   }
5054   if ( myElementsOnShapePtr )
5055     myElementsOnShapePtr->SetMesh( myMeshDS );
5056 }
5057
5058 void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
5059 {
5060   if ( myShape != theShape )
5061   {
5062     myShape = theShape;
5063     init();
5064   }
5065 }
5066
5067 static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
5068                         const TopoDS_Shape&               theShape)
5069 {
5070   if (theMap.Contains(theShape)) return true;
5071
5072   if (theShape.ShapeType() == TopAbs_COMPOUND ||
5073       theShape.ShapeType() == TopAbs_COMPSOLID)
5074   {
5075     TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
5076     for (; anIt.More(); anIt.Next())
5077     {
5078       if (!IsSubShape(theMap, anIt.Value())) {
5079         return false;
5080       }
5081     }
5082     return true;
5083   }
5084
5085   return false;
5086 }
5087
5088 void BelongToGeom::init()
5089 {
5090   if ( !myMeshDS || myShape.IsNull() ) return;
5091
5092   // is sub-shape of main shape?
5093   TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
5094   if (aMainShape.IsNull()) {
5095     myIsSubshape = false;
5096   }
5097   else {
5098     TopTools_IndexedMapOfShape aMap;
5099     TopExp::MapShapes( aMainShape, aMap );
5100     myIsSubshape = IsSubShape( aMap, myShape );
5101     if ( myIsSubshape )
5102     {
5103       aMap.Clear();
5104       TopExp::MapShapes( myShape, aMap );
5105       mySubShapesIDs.Clear();
5106       for ( int i = 1; i <= aMap.Extent(); ++i )
5107       {
5108         int subID = myMeshDS->ShapeToIndex( aMap( i ));
5109         if ( subID > 0 )
5110           mySubShapesIDs.Add( subID );
5111       }
5112     }
5113   }
5114
5115   //if (!myIsSubshape) // to be always ready to check an element not bound to geometry
5116   {
5117     if ( !myElementsOnShapePtr )
5118       myElementsOnShapePtr.reset( new ElementsOnShape() );
5119     myElementsOnShapePtr->SetTolerance( myTolerance );
5120     myElementsOnShapePtr->SetAllNodes( true ); // "belong", while false means "lays on"
5121     myElementsOnShapePtr->SetMesh( myMeshDS );
5122     myElementsOnShapePtr->SetShape( myShape, myType );
5123   }
5124 }
5125
5126 bool BelongToGeom::IsSatisfy (long theId)
5127 {
5128   if (myMeshDS == 0 || myShape.IsNull())
5129     return false;
5130
5131   if (!myIsSubshape)
5132   {
5133     return myElementsOnShapePtr->IsSatisfy(theId);
5134   }
5135
5136   // Case of sub-mesh
5137
5138   if (myType == SMDSAbs_Node)
5139   {
5140     if ( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ))
5141     {
5142       if ( aNode->getshapeId() < 1 )
5143         return myElementsOnShapePtr->IsSatisfy(theId);
5144       else
5145         return mySubShapesIDs.Contains( aNode->getshapeId() );
5146     }
5147   }
5148   else
5149   {
5150     if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
5151     {
5152       if ( myType == SMDSAbs_All || anElem->GetType() == myType )
5153       {
5154         if ( anElem->getshapeId() < 1 )
5155           return myElementsOnShapePtr->IsSatisfy(theId);
5156         else
5157           return mySubShapesIDs.Contains( anElem->getshapeId() );
5158       }
5159     }
5160   }
5161
5162   return false;
5163 }
5164
5165 void BelongToGeom::SetType (SMDSAbs_ElementType theType)
5166 {
5167   if ( myType != theType )
5168   {
5169     myType = theType;
5170     init();
5171   }
5172 }
5173
5174 SMDSAbs_ElementType BelongToGeom::GetType() const
5175 {
5176   return myType;
5177 }
5178
5179 TopoDS_Shape BelongToGeom::GetShape()
5180 {
5181   return myShape;
5182 }
5183
5184 const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const
5185 {
5186   return myMeshDS;
5187 }
5188
5189 void BelongToGeom::SetTolerance (double theTolerance)
5190 {
5191   myTolerance = theTolerance;
5192   init();
5193 }
5194
5195 double BelongToGeom::GetTolerance()
5196 {
5197   return myTolerance;
5198 }
5199
5200 /*
5201   Class       : LyingOnGeom
5202   Description : Predicate for verifying whether entiy lying or partially lying on
5203   specified geometrical support
5204 */
5205
5206 LyingOnGeom::LyingOnGeom()
5207   : myMeshDS(NULL),
5208     myType(SMDSAbs_NbElementTypes),
5209     myIsSubshape(false),
5210     myTolerance(Precision::Confusion())
5211 {}
5212
5213 Predicate* LyingOnGeom::clone() const
5214 {
5215   LyingOnGeom* cln = 0;
5216   if ( myElementsOnShapePtr )
5217     if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
5218     {
5219       cln = new LyingOnGeom( *this );
5220       cln->myElementsOnShapePtr.reset( eos );
5221     }
5222   return cln;
5223 }
5224
5225 void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
5226 {
5227   if ( myMeshDS != theMesh )
5228   {
5229     myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
5230     init();
5231   }
5232   if ( myElementsOnShapePtr )
5233     myElementsOnShapePtr->SetMesh( myMeshDS );
5234 }
5235
5236 void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
5237 {
5238   if ( myShape != theShape )
5239   {
5240     myShape = theShape;
5241     init();
5242   }
5243 }
5244
5245 void LyingOnGeom::init()
5246 {
5247   if (!myMeshDS || myShape.IsNull()) return;
5248
5249   // is sub-shape of main shape?
5250   TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
5251   if (aMainShape.IsNull()) {
5252     myIsSubshape = false;
5253   }
5254   else {
5255     myIsSubshape = myMeshDS->IsGroupOfSubShapes( myShape );
5256   }
5257
5258   if (myIsSubshape)
5259   {
5260     TopTools_IndexedMapOfShape shapes;
5261     TopExp::MapShapes( myShape, shapes );
5262     mySubShapesIDs.Clear();
5263     for ( int i = 1; i <= shapes.Extent(); ++i )
5264     {
5265       int subID = myMeshDS->ShapeToIndex( shapes( i ));
5266       if ( subID > 0 )
5267         mySubShapesIDs.Add( subID );
5268     }
5269   }
5270   // else // to be always ready to check an element not bound to geometry
5271   {
5272     if ( !myElementsOnShapePtr )
5273       myElementsOnShapePtr.reset( new ElementsOnShape() );
5274     myElementsOnShapePtr->SetTolerance( myTolerance );
5275     myElementsOnShapePtr->SetAllNodes( false ); // lays on, while true means "belong"
5276     myElementsOnShapePtr->SetMesh( myMeshDS );
5277     myElementsOnShapePtr->SetShape( myShape, myType );
5278   }
5279 }
5280
5281 bool LyingOnGeom::IsSatisfy( long theId )
5282 {
5283   if ( myMeshDS == 0 || myShape.IsNull() )
5284     return false;
5285
5286   if (!myIsSubshape)
5287   {
5288     return myElementsOnShapePtr->IsSatisfy(theId);
5289   }
5290
5291   // Case of sub-mesh
5292
5293   const SMDS_MeshElement* elem =
5294     ( myType == SMDSAbs_Node ) ? myMeshDS->FindNode( theId ) : myMeshDS->FindElement( theId );
5295
5296   if ( mySubShapesIDs.Contains( elem->getshapeId() ))
5297     return true;
5298
5299   if (( elem->GetType() != SMDSAbs_Node ) &&
5300       ( myType == SMDSAbs_All || elem->GetType() == myType ))
5301   {
5302     SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator();
5303     while ( nodeItr->more() )
5304     {
5305       const SMDS_MeshElement* aNode = nodeItr->next();
5306       if ( mySubShapesIDs.Contains( aNode->getshapeId() ))
5307         return true;
5308     }
5309   }
5310
5311   return false;
5312 }
5313
5314 void LyingOnGeom::SetType( SMDSAbs_ElementType theType )
5315 {
5316   if ( myType != theType )
5317   {
5318     myType = theType;
5319     init();
5320   }
5321 }
5322
5323 SMDSAbs_ElementType LyingOnGeom::GetType() const
5324 {
5325   return myType;
5326 }
5327
5328 TopoDS_Shape LyingOnGeom::GetShape()
5329 {
5330   return myShape;
5331 }
5332
5333 const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const
5334 {
5335   return myMeshDS;
5336 }
5337
5338 void LyingOnGeom::SetTolerance (double theTolerance)
5339 {
5340   myTolerance = theTolerance;
5341   init();
5342 }
5343
5344 double LyingOnGeom::GetTolerance()
5345 {
5346   return myTolerance;
5347 }
5348
5349 TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0)
5350 {}
5351
5352 TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n), myElem(0)
5353 {}
5354
5355 TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t), myElem(0)
5356 {}
5357
5358 TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray), myElem(theSequenceOfXYZ.myElem)
5359 {}
5360
5361 template <class InputIterator>
5362 TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd), myElem(0)
5363 {}
5364
5365 TSequenceOfXYZ::~TSequenceOfXYZ()
5366 {}
5367
5368 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
5369 {
5370   myArray = theSequenceOfXYZ.myArray;
5371   myElem  = theSequenceOfXYZ.myElem;
5372   return *this;
5373 }
5374
5375 gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
5376 {
5377   return myArray[n-1];
5378 }
5379
5380 const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
5381 {
5382   return myArray[n-1];
5383 }
5384
5385 void TSequenceOfXYZ::clear()
5386 {
5387   myArray.clear();
5388 }
5389
5390 void TSequenceOfXYZ::reserve(size_type n)
5391 {
5392   myArray.reserve(n);
5393 }
5394
5395 void TSequenceOfXYZ::push_back(const gp_XYZ& v)
5396 {
5397   myArray.push_back(v);
5398 }
5399
5400 TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
5401 {
5402   return myArray.size();
5403 }
5404
5405 SMDSAbs_EntityType TSequenceOfXYZ::getElementEntity() const
5406 {
5407   return myElem ? myElem->GetEntityType() : SMDSEntity_Last;
5408 }
5409
5410 TMeshModifTracer::TMeshModifTracer():
5411   myMeshModifTime(0), myMesh(0)
5412 {
5413 }
5414 void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh )
5415 {
5416   if ( theMesh != myMesh )
5417     myMeshModifTime = 0;
5418   myMesh = theMesh;
5419 }
5420 bool TMeshModifTracer::IsMeshModified()
5421 {
5422   bool modified = false;
5423   if ( myMesh )
5424   {
5425     modified = ( myMeshModifTime != myMesh->GetMTime() );
5426     myMeshModifTime = myMesh->GetMTime();
5427   }
5428   return modified;
5429 }