src/Controls/SMESH_Controls.cxx

   1 // Copyright (C) 2007-2022  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 all its nodes are on the boundary.
2473   // A node is on the boundary if it is connected to one or more faces.
2474   SMDS_VolumeTool myTool;
2475   if (myTool.Set(myMesh->FindElement(theElementId)))
2476   {
2477     auto nodes = myTool.GetNodes();
2478
2479     for (int i = 0; i < myTool.NbNodes(); ++i)
2480     {
2481       auto node = nodes[i];
2482       if (node->NbInverseElements(SMDSAbs_Face) == 0)
2483       {
2484         return false;
2485       }
2486     }
2487     return true;
2488   }
2489   return false;
2490 }
2491
2492 //================================================================================
2493 /*
2494   Class       : OverConstrainedFace
2495 */
2496 //================================================================================
2497
2498 bool OverConstrainedFace::IsSatisfy(long theElementId )
2499 {
2500   // An element is over-constrained if all its nodes are on the boundary.
2501   // A node is on the boundary if it is connected to one or more faces.
2502   if (const SMDS_MeshElement *face = myMesh->FindElement(theElementId))
2503     if (face->GetType() == SMDSAbs_Face)
2504     {
2505       int nbN = face->NbCornerNodes();
2506       for (int i = 0; i < nbN; ++i)
2507       {
2508         const SMDS_MeshNode *n1 = face->GetNode(i);
2509         if (n1->NbInverseElements(SMDSAbs_Edge) == 0)
2510           return false;
2511       }
2512       return true;
2513     }
2514   return false;
2515 }
2516
2517 //================================================================================
2518 /*
2519   Class       : CoincidentNodes
2520   Description : Predicate of Coincident nodes
2521 */
2522 //================================================================================
2523
2524 CoincidentNodes::CoincidentNodes()
2525 {
2526   myToler = 1e-5;
2527 }
2528
2529 bool CoincidentNodes::IsSatisfy( long theElementId )
2530 {
2531   return myCoincidentIDs.Contains( theElementId );
2532 }
2533
2534 SMDSAbs_ElementType CoincidentNodes::GetType() const
2535 {
2536   return SMDSAbs_Node;
2537 }
2538
2539 void CoincidentNodes::SetTolerance( const double theToler )
2540 {
2541   if ( myToler != theToler )
2542   {
2543     SetMesh(0);
2544     myToler = theToler;
2545   }
2546 }
2547
2548 void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
2549 {
2550   myMeshModifTracer.SetMesh( theMesh );
2551   if ( myMeshModifTracer.IsMeshModified() )
2552   {
2553     TIDSortedNodeSet nodesToCheck;
2554     SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
2555     while ( nIt->more() )
2556       nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
2557
2558     std::list< std::list< const SMDS_MeshNode*> > nodeGroups;
2559     SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
2560
2561     myCoincidentIDs.Clear();
2562     std::list< std::list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
2563     for ( ; groupIt != nodeGroups.end(); ++groupIt )
2564     {
2565       std::list< const SMDS_MeshNode*>& coincNodes = *groupIt;
2566       std::list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
2567       for ( ; n != coincNodes.end(); ++n )
2568         myCoincidentIDs.Add( (*n)->GetID() );
2569     }
2570   }
2571 }
2572
2573 //================================================================================
2574 /*
2575   Class       : CoincidentElements
2576   Description : Predicate of Coincident Elements
2577   Note        : This class is suitable only for visualization of Coincident Elements
2578 */
2579 //================================================================================
2580
2581 CoincidentElements::CoincidentElements()
2582 {
2583   myMesh = 0;
2584 }
2585
2586 void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
2587 {
2588   myMesh = theMesh;
2589 }
2590
2591 bool CoincidentElements::IsSatisfy( long theElementId )
2592 {
2593   if ( !myMesh ) return false;
2594
2595   if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
2596   {
2597     if ( e->GetType() != GetType() ) return false;
2598     std::set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
2599     const int nbNodes = e->NbNodes();
2600     SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
2601     while ( invIt->more() )
2602     {
2603       const SMDS_MeshElement* e2 = invIt->next();
2604       if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
2605
2606       bool sameNodes = true;
2607       for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
2608         sameNodes = ( elemNodes.count( e2->GetNode( i )));
2609       if ( sameNodes )
2610         return true;
2611     }
2612   }
2613   return false;
2614 }
2615
2616 SMDSAbs_ElementType CoincidentElements1D::GetType() const
2617 {
2618   return SMDSAbs_Edge;
2619 }
2620 SMDSAbs_ElementType CoincidentElements2D::GetType() const
2621 {
2622   return SMDSAbs_Face;
2623 }
2624 SMDSAbs_ElementType CoincidentElements3D::GetType() const
2625 {
2626   return SMDSAbs_Volume;
2627 }
2628
2629
2630 //================================================================================
2631 /*
2632   Class       : FreeBorders
2633   Description : Predicate for free borders
2634 */
2635 //================================================================================
2636
2637 FreeBorders::FreeBorders()
2638 {
2639   myMesh = 0;
2640 }
2641
2642 void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
2643 {
2644   myMesh = theMesh;
2645 }
2646
2647 bool FreeBorders::IsSatisfy( long theId )
2648 {
2649   return getNbMultiConnection( myMesh, theId ) == 1;
2650 }
2651
2652 SMDSAbs_ElementType FreeBorders::GetType() const
2653 {
2654   return SMDSAbs_Edge;
2655 }
2656
2657
2658 //================================================================================
2659 /*
2660   Class       : FreeEdges
2661   Description : Predicate for free Edges
2662 */
2663 //================================================================================
2664
2665 FreeEdges::FreeEdges()
2666 {
2667   myMesh = 0;
2668 }
2669
2670 void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
2671 {
2672   myMesh = theMesh;
2673 }
2674
2675 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const smIdType theFaceId  )
2676 {
2677   SMDS_ElemIteratorPtr anElemIter = theNodes[ 0 ]->GetInverseElementIterator(SMDSAbs_Face);
2678   while( anElemIter->more() )
2679   {
2680     if ( const SMDS_MeshElement* anElem = anElemIter->next())
2681     {
2682       const smIdType anId = anElem->GetID();
2683       if ( anId != theFaceId && anElem->GetNodeIndex( theNodes[1] ) >= 0 )
2684         return false;
2685     }
2686   }
2687   return true;
2688 }
2689
2690 bool FreeEdges::IsSatisfy( long theId )
2691 {
2692   if ( myMesh == 0 )
2693     return false;
2694
2695   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2696   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
2697     return false;
2698
2699   SMDS_NodeIteratorPtr anIter = aFace->interlacedNodesIterator();
2700   if ( !anIter )
2701     return false;
2702
2703   int i = 0, nbNodes = aFace->NbNodes();
2704   std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
2705   while( anIter->more() )
2706     if ( ! ( aNodes[ i++ ] = anIter->next() ))
2707       return false;
2708   aNodes[ nbNodes ] = aNodes[ 0 ];
2709
2710   for ( i = 0; i < nbNodes; i++ )
2711     if ( IsFreeEdge( &aNodes[ i ], theId ) )
2712       return true;
2713
2714   return false;
2715 }
2716
2717 SMDSAbs_ElementType FreeEdges::GetType() const
2718 {
2719   return SMDSAbs_Face;
2720 }
2721
2722 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
2723   myElemId(theElemId)
2724 {
2725   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
2726   if(thePntId1 > thePntId2){
2727     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
2728   }
2729 }
2730
2731 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
2732   if(myPntId[0] < x.myPntId[0]) return true;
2733   if(myPntId[0] == x.myPntId[0])
2734     if(myPntId[1] < x.myPntId[1]) return true;
2735   return false;
2736 }
2737
2738 inline void UpdateBorders(const FreeEdges::Border& theBorder,
2739                           FreeEdges::TBorders& theRegistry,
2740                           FreeEdges::TBorders& theContainer)
2741 {
2742   if(theRegistry.find(theBorder) == theRegistry.end()){
2743     theRegistry.insert(theBorder);
2744     theContainer.insert(theBorder);
2745   }else{
2746     theContainer.erase(theBorder);
2747   }
2748 }
2749
2750 void FreeEdges::GetBoreders(TBorders& theBorders)
2751 {
2752   TBorders aRegistry;
2753   for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
2754   {
2755     const SMDS_MeshFace* anElem = anIter->next();
2756     long anElemId = anElem->GetID();
2757     SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
2758     if ( !aNodesIter->more() ) continue;
2759     long aNodeId[2] = {0,0};
2760     aNodeId[0] = anElem->GetNode( anElem->NbNodes()-1 )->GetID();
2761     for ( ; aNodesIter->more(); )
2762     {
2763       aNodeId[1] = aNodesIter->next()->GetID();
2764       Border aBorder( anElemId, aNodeId[0], aNodeId[1] );
2765       UpdateBorders( aBorder, aRegistry, theBorders );
2766       aNodeId[0] = aNodeId[1];
2767     }
2768   }
2769 }
2770
2771 //================================================================================
2772 /*
2773   Class       : FreeNodes
2774   Description : Predicate for free nodes
2775 */
2776 //================================================================================
2777
2778 FreeNodes::FreeNodes()
2779 {
2780   myMesh = 0;
2781 }
2782
2783 void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
2784 {
2785   myMesh = theMesh;
2786 }
2787
2788 bool FreeNodes::IsSatisfy( long theNodeId )
2789 {
2790   const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
2791   if (!aNode)
2792     return false;
2793
2794   return (aNode->NbInverseElements() < 1);
2795 }
2796
2797 SMDSAbs_ElementType FreeNodes::GetType() const
2798 {
2799   return SMDSAbs_Node;
2800 }
2801
2802
2803 //================================================================================
2804 /*
2805   Class       : FreeFaces
2806   Description : Predicate for free faces
2807 */
2808 //================================================================================
2809
2810 FreeFaces::FreeFaces()
2811 {
2812   myMesh = 0;
2813 }
2814
2815 void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
2816 {
2817   myMesh = theMesh;
2818 }
2819
2820 bool FreeFaces::IsSatisfy( long theId )
2821 {
2822   if (!myMesh) return false;
2823   // check that faces nodes refers to less than two common volumes
2824   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2825   if ( !aFace || aFace->GetType() != SMDSAbs_Face )
2826     return false;
2827
2828   int nbNode = aFace->NbNodes();
2829
2830   // collect volumes to check that number of volumes with count equal nbNode not less than 2
2831   typedef std::map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
2832   typedef std::map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
2833   TMapOfVolume mapOfVol;
2834
2835   SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
2836   while ( nodeItr->more() )
2837   {
2838     const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
2839     if ( !aNode ) continue;
2840     SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
2841     while ( volItr->more() )
2842     {
2843       SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
2844       TItrMapOfVolume    itr = mapOfVol.insert( std::make_pair( aVol, 0 )).first;
2845       (*itr).second++;
2846     }
2847   }
2848   int nbVol = 0;
2849   TItrMapOfVolume volItr = mapOfVol.begin();
2850   TItrMapOfVolume volEnd = mapOfVol.end();
2851   for ( ; volItr != volEnd; ++volItr )
2852     if ( (*volItr).second >= nbNode )
2853       nbVol++;
2854   // face is not free if number of volumes constructed on their nodes more than one
2855   return (nbVol < 2);
2856 }
2857
2858 SMDSAbs_ElementType FreeFaces::GetType() const
2859 {
2860   return SMDSAbs_Face;
2861 }
2862
2863 //================================================================================
2864 /*
2865   Class       : LinearOrQuadratic
2866   Description : Predicate to verify whether a mesh element is linear
2867 */
2868 //================================================================================
2869
2870 LinearOrQuadratic::LinearOrQuadratic()
2871 {
2872   myMesh = 0;
2873 }
2874
2875 void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
2876 {
2877   myMesh = theMesh;
2878 }
2879
2880 bool LinearOrQuadratic::IsSatisfy( long theId )
2881 {
2882   if (!myMesh) return false;
2883   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2884   if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
2885     return false;
2886   return (!anElem->IsQuadratic());
2887 }
2888
2889 void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
2890 {
2891   myType = theType;
2892 }
2893
2894 SMDSAbs_ElementType LinearOrQuadratic::GetType() const
2895 {
2896   return myType;
2897 }
2898
2899 //================================================================================
2900 /*
2901   Class       : GroupColor
2902   Description : Functor for check color of group to which mesh element belongs to
2903 */
2904 //================================================================================
2905
2906 GroupColor::GroupColor()
2907 {
2908 }
2909
2910 bool GroupColor::IsSatisfy( long theId )
2911 {
2912   return myIDs.count( theId );
2913 }
2914
2915 void GroupColor::SetType( SMDSAbs_ElementType theType )
2916 {
2917   myType = theType;
2918 }
2919
2920 SMDSAbs_ElementType GroupColor::GetType() const
2921 {
2922   return myType;
2923 }
2924
2925 static bool isEqual( const Quantity_Color& theColor1,
2926                      const Quantity_Color& theColor2 )
2927 {
2928   // tolerance to compare colors
2929   const double tol = 5*1e-3;
2930   return ( fabs( theColor1.Red()   - theColor2.Red() )   < tol &&
2931            fabs( theColor1.Green() - theColor2.Green() ) < tol &&
2932            fabs( theColor1.Blue()  - theColor2.Blue() )  < tol );
2933 }
2934
2935 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
2936 {
2937   myIDs.clear();
2938
2939   const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
2940   if ( !aMesh )
2941     return;
2942
2943   int nbGrp = aMesh->GetNbGroups();
2944   if ( !nbGrp )
2945     return;
2946
2947   // iterates on groups and find necessary elements ids
2948   const std::set<SMESHDS_GroupBase*>&       aGroups = aMesh->GetGroups();
2949   std::set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
2950   for (; GrIt != aGroups.end(); GrIt++)
2951   {
2952     SMESHDS_GroupBase* aGrp = (*GrIt);
2953     if ( !aGrp )
2954       continue;
2955     // check type and color of group
2956     if ( !isEqual( myColor, aGrp->GetColor() ))
2957       continue;
2958
2959     // IPAL52867 (prevent infinite recursion via GroupOnFilter)
2960     if ( SMESHDS_GroupOnFilter * gof = dynamic_cast< SMESHDS_GroupOnFilter* >( aGrp ))
2961       if ( gof->GetPredicate().get() == this )
2962         continue;
2963
2964     SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
2965     if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
2966       // add elements IDS into control
2967       smIdType aSize = aGrp->Extent();
2968       for (smIdType i = 0; i < aSize; i++)
2969         myIDs.insert( aGrp->GetID(i+1) );
2970     }
2971   }
2972 }
2973
2974 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
2975 {
2976   Kernel_Utils::Localizer loc;
2977   TCollection_AsciiString aStr = theStr;
2978   aStr.RemoveAll( ' ' );
2979   aStr.RemoveAll( '\t' );
2980   for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
2981     aStr.Remove( aPos, 2 );
2982   Standard_Real clr[3];
2983   clr[0] = clr[1] = clr[2] = 0.;
2984   for ( int i = 0; i < 3; i++ ) {
2985     TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
2986     if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
2987       clr[i] = tmpStr.RealValue();
2988   }
2989   myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
2990 }
2991
2992 //=======================================================================
2993 // name    : GetRangeStr
2994 // Purpose : Get range as a string.
2995 //           Example: "1,2,3,50-60,63,67,70-"
2996 //=======================================================================
2997
2998 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
2999 {
3000   theResStr.Clear();
3001   theResStr += TCollection_AsciiString( myColor.Red() );
3002   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
3003   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
3004 }
3005
3006 //================================================================================
3007 /*
3008   Class       : ElemGeomType
3009   Description : Predicate to check element geometry type
3010 */
3011 //================================================================================
3012
3013 ElemGeomType::ElemGeomType()
3014 {
3015   myMesh = 0;
3016   myType = SMDSAbs_All;
3017   myGeomType = SMDSGeom_TRIANGLE;
3018 }
3019
3020 void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
3021 {
3022   myMesh = theMesh;
3023 }
3024
3025 bool ElemGeomType::IsSatisfy( long theId )
3026 {
3027   if (!myMesh) return false;
3028   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3029   if ( !anElem )
3030     return false;
3031   const SMDSAbs_ElementType anElemType = anElem->GetType();
3032   if ( myType != SMDSAbs_All && anElemType != myType )
3033     return false;
3034   bool isOk = ( anElem->GetGeomType() == myGeomType );
3035   return isOk;
3036 }
3037
3038 void ElemGeomType::SetType( SMDSAbs_ElementType theType )
3039 {
3040   myType = theType;
3041 }
3042
3043 SMDSAbs_ElementType ElemGeomType::GetType() const
3044 {
3045   return myType;
3046 }
3047
3048 void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
3049 {
3050   myGeomType = theType;
3051 }
3052
3053 SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
3054 {
3055   return myGeomType;
3056 }
3057
3058 //================================================================================
3059 /*
3060   Class       : ElemEntityType
3061   Description : Predicate to check element entity type
3062 */
3063 //================================================================================
3064
3065 ElemEntityType::ElemEntityType():
3066   myMesh( 0 ),
3067   myType( SMDSAbs_All ),
3068   myEntityType( SMDSEntity_0D )
3069 {
3070 }
3071
3072 void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh )
3073 {
3074   myMesh = theMesh;
3075 }
3076
3077 bool ElemEntityType::IsSatisfy( long theId )
3078 {
3079   if ( !myMesh ) return false;
3080   if ( myType == SMDSAbs_Node )
3081     return myMesh->FindNode( theId );
3082   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3083   return ( anElem &&
3084            myEntityType == anElem->GetEntityType() );
3085 }
3086
3087 void ElemEntityType::SetType( SMDSAbs_ElementType theType )
3088 {
3089   myType = theType;
3090 }
3091
3092 SMDSAbs_ElementType ElemEntityType::GetType() const
3093 {
3094   return myType;
3095 }
3096
3097 void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType )
3098 {
3099   myEntityType = theEntityType;
3100 }
3101
3102 SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const
3103 {
3104   return myEntityType;
3105 }
3106
3107 //================================================================================
3108 /*!
3109  * \brief Class ConnectedElements
3110  */
3111 //================================================================================
3112
3113 ConnectedElements::ConnectedElements():
3114   myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {}
3115
3116 SMDSAbs_ElementType ConnectedElements::GetType() const
3117 { return myType; }
3118
3119 smIdType ConnectedElements::GetNode() const
3120 { return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
3121
3122 std::vector<double> ConnectedElements::GetPoint() const
3123 { return myXYZ; }
3124
3125 void ConnectedElements::clearOkIDs()
3126 { myOkIDsReady = false; myOkIDs.clear(); }
3127
3128 void ConnectedElements::SetType( SMDSAbs_ElementType theType )
3129 {
3130   if ( myType != theType || myMeshModifTracer.IsMeshModified() )
3131     clearOkIDs();
3132   myType = theType;
3133 }
3134
3135 void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh )
3136 {
3137   myMeshModifTracer.SetMesh( theMesh );
3138   if ( myMeshModifTracer.IsMeshModified() )
3139   {
3140     clearOkIDs();
3141     if ( !myXYZ.empty() )
3142       SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh
3143   }
3144 }
3145
3146 void ConnectedElements::SetNode( smIdType nodeID )
3147 {
3148   myNodeID = nodeID;
3149   myXYZ.clear();
3150
3151   bool isSameDomain = false;
3152   if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() )
3153     if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID ))
3154     {
3155       SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType );
3156       while ( !isSameDomain && eIt->more() )
3157         isSameDomain = IsSatisfy( eIt->next()->GetID() );
3158     }
3159   if ( !isSameDomain )
3160     clearOkIDs();
3161 }
3162
3163 void ConnectedElements::SetPoint( double x, double y, double z )
3164 {
3165   myXYZ.resize(3);
3166   myXYZ[0] = x;
3167   myXYZ[1] = y;
3168   myXYZ[2] = z;
3169   myNodeID = 0;
3170
3171   bool isSameDomain = false;
3172
3173   // find myNodeID by myXYZ if possible
3174   if ( myMeshModifTracer.GetMesh() )
3175   {
3176     SMESHUtils::Deleter<SMESH_ElementSearcher> searcher
3177       ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
3178
3179     std::vector< const SMDS_MeshElement* > foundElems;
3180     searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
3181
3182     if ( !foundElems.empty() )
3183     {
3184       myNodeID = foundElems[0]->GetNode(0)->GetID();
3185       if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() )
3186         isSameDomain = IsSatisfy( foundElems[0]->GetID() );
3187     }
3188   }
3189   if ( !isSameDomain )
3190     clearOkIDs();
3191 }
3192
3193 bool ConnectedElements::IsSatisfy( long theElementId )
3194 {
3195   // Here we do NOT check if the mesh has changed, we do it in Set...() only!!!
3196
3197   if ( !myOkIDsReady )
3198   {
3199     if ( !myMeshModifTracer.GetMesh() )
3200       return false;
3201     const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID );
3202     if ( !node0 )
3203       return false;
3204
3205     std::list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
3206     std::set< smIdType > checkedNodeIDs;
3207     // algo:
3208     // foreach node in nodeQueue:
3209     //   foreach element sharing a node:
3210     //     add ID of an element of myType to myOkIDs;
3211     //     push all element nodes absent from checkedNodeIDs to nodeQueue;
3212     while ( !nodeQueue.empty() )
3213     {
3214       const SMDS_MeshNode* node = nodeQueue.front();
3215       nodeQueue.pop_front();
3216
3217       // loop on elements sharing the node
3218       SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
3219       while ( eIt->more() )
3220       {
3221         // keep elements of myType
3222         const SMDS_MeshElement* element = eIt->next();
3223         if ( myType == SMDSAbs_All || element->GetType() == myType )
3224           myOkIDs.insert( myOkIDs.end(), element->GetID() );
3225
3226         // enqueue nodes of the element
3227         SMDS_ElemIteratorPtr nIt = element->nodesIterator();
3228         while ( nIt->more() )
3229         {
3230           const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() );
3231           if ( checkedNodeIDs.insert( n->GetID()).second )
3232             nodeQueue.push_back( n );
3233         }
3234       }
3235     }
3236     if ( myType == SMDSAbs_Node )
3237       std::swap( myOkIDs, checkedNodeIDs );
3238
3239     size_t totalNbElems = myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType );
3240     if ( myOkIDs.size() == totalNbElems )
3241       myOkIDs.clear();
3242
3243     myOkIDsReady = true;
3244   }
3245
3246   return myOkIDs.empty() ? true : myOkIDs.count( theElementId );
3247 }
3248
3249 //================================================================================
3250 /*!
3251  * \brief Class CoplanarFaces
3252  */
3253 //================================================================================
3254
3255 namespace
3256 {
3257   inline bool isLessAngle( const gp_Vec& v1, const gp_Vec& v2, const double cos )
3258   {
3259     double dot = v1 * v2; // cos * |v1| * |v2|
3260     double l1  = v1.SquareMagnitude();
3261     double l2  = v2.SquareMagnitude();
3262     return (( dot * cos >= 0 ) &&
3263             ( dot * dot ) / l1 / l2 >= ( cos * cos ));
3264   }
3265 }
3266 CoplanarFaces::CoplanarFaces()
3267   : myFaceID(0), myToler(0)
3268 {
3269 }
3270 void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
3271 {
3272   myMeshModifTracer.SetMesh( theMesh );
3273   if ( myMeshModifTracer.IsMeshModified() )
3274   {
3275     // Build a set of coplanar face ids
3276
3277     myCoplanarIDs.Clear();
3278
3279     if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler )
3280       return;
3281
3282     const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID );
3283     if ( !face || face->GetType() != SMDSAbs_Face )
3284       return;
3285
3286     bool normOK;
3287     gp_Vec myNorm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
3288     if (!normOK)
3289       return;
3290
3291     const double cosTol = Cos( myToler * M_PI / 180. );
3292     NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks;
3293
3294     std::list< std::pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
3295     faceQueue.push_back( std::make_pair( face, myNorm ));
3296     while ( !faceQueue.empty() )
3297     {
3298       face   = faceQueue.front().first;
3299       myNorm = faceQueue.front().second;
3300       faceQueue.pop_front();
3301
3302       for ( int i = 0, nbN = face->NbCornerNodes(); i < nbN; ++i )
3303       {
3304         const SMDS_MeshNode*  n1 = face->GetNode( i );
3305         const SMDS_MeshNode*  n2 = face->GetNode(( i+1 )%nbN);
3306         if ( !checkedLinks.Add( SMESH_TLink( n1, n2 )))
3307           continue;
3308         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face);
3309         while ( fIt->more() )
3310         {
3311           const SMDS_MeshElement* f = fIt->next();
3312           if ( f->GetNodeIndex( n2 ) > -1 )
3313           {
3314             gp_Vec norm = getNormale( static_cast<const SMDS_MeshFace*>(f), &normOK );
3315             if (!normOK || isLessAngle( myNorm, norm, cosTol))
3316             {
3317               myCoplanarIDs.Add( f->GetID() );
3318               faceQueue.push_back( std::make_pair( f, norm ));
3319             }
3320           }
3321         }
3322       }
3323     }
3324   }
3325 }
3326 bool CoplanarFaces::IsSatisfy( long theElementId )
3327 {
3328   return myCoplanarIDs.Contains( theElementId );
3329 }
3330
3331 /*
3332  *Class       : RangeOfIds
3333   *Description : Predicate for Range of Ids.
3334   *              Range may be specified with two ways.
3335   *              1. Using AddToRange method
3336   *              2. With SetRangeStr method. Parameter of this method is a string
3337   *                 like as "1,2,3,50-60,63,67,70-"
3338 */
3339
3340 //=======================================================================
3341 // name    : RangeOfIds
3342 // Purpose : Constructor
3343 //=======================================================================
3344 RangeOfIds::RangeOfIds()
3345 {
3346   myMesh = 0;
3347   myType = SMDSAbs_All;
3348 }
3349
3350 //=======================================================================
3351 // name    : SetMesh
3352 // Purpose : Set mesh
3353 //=======================================================================
3354 void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh )
3355 {
3356   myMesh = theMesh;
3357 }
3358
3359 //=======================================================================
3360 // name    : AddToRange
3361 // Purpose : Add ID to the range
3362 //=======================================================================
3363 bool RangeOfIds::AddToRange( long theEntityId )
3364 {
3365   myIds.Add( theEntityId );
3366   return true;
3367 }
3368
3369 //=======================================================================
3370 // name    : GetRangeStr
3371 // Purpose : Get range as a string.
3372 //           Example: "1,2,3,50-60,63,67,70-"
3373 //=======================================================================
3374 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
3375 {
3376   theResStr.Clear();
3377
3378   TIDsSeq                             anIntSeq;
3379   NCollection_Sequence< std::string > aStrSeq;
3380
3381   TIDsMap::Iterator anIter( myIds );
3382   for ( ; anIter.More(); anIter.Next() )
3383   {
3384     smIdType anId = anIter.Key();
3385     SMESH_Comment aStr( anId );
3386     anIntSeq.Append( anId );
3387     aStrSeq.Append( aStr );
3388   }
3389
3390   for ( smIdType i = 1, n = myMin.size(); i <= n; i++ )
3391   {
3392     smIdType aMinId = myMin[i];
3393     smIdType aMaxId = myMax[i];
3394
3395     SMESH_Comment aStr;
3396     if ( aMinId != IntegerFirst() )
3397       aStr << aMinId;
3398
3399     aStr << "-";
3400
3401     if ( aMaxId != std::numeric_limits<smIdType>::max() )
3402       aStr << aMaxId;
3403
3404     // find position of the string in result sequence and insert string in it
3405     if ( anIntSeq.Length() == 0 )
3406     {
3407       anIntSeq.Append( aMinId );
3408       aStrSeq.Append( (const char*)aStr );
3409     }
3410     else
3411     {
3412       if ( aMinId < anIntSeq.First() )
3413       {
3414         anIntSeq.Prepend( aMinId );
3415         aStrSeq.Prepend( (const char*)aStr );
3416       }
3417       else if ( aMinId > anIntSeq.Last() )
3418       {
3419         anIntSeq.Append( aMinId );
3420         aStrSeq.Append( (const char*)aStr );
3421       }
3422       else
3423         for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
3424           if ( aMinId < anIntSeq( j ) )
3425           {
3426             anIntSeq.InsertBefore( j, aMinId );
3427             aStrSeq.InsertBefore( j, (const char*)aStr );
3428             break;
3429           }
3430     }
3431   }
3432
3433   if ( aStrSeq.Length() == 0 )
3434     return;
3435   std::string aResStr;
3436   aResStr = aStrSeq( 1 );
3437   for ( int j = 2, k = aStrSeq.Length(); j <= k; j++  )
3438   {
3439     aResStr += ",";
3440     aResStr += aStrSeq( j );
3441   }
3442   theResStr = aResStr.c_str();
3443 }
3444
3445 //=======================================================================
3446 // name    : SetRangeStr
3447 // Purpose : Define range with string
3448 //           Example of entry string: "1,2,3,50-60,63,67,70-"
3449 //=======================================================================
3450 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
3451 {
3452   myMin.clear();
3453   myMax.clear();
3454   myIds.Clear();
3455
3456   TCollection_AsciiString aStr = theStr;
3457   for ( int i = 1; i <= aStr.Length(); ++i )
3458   {
3459     char c = aStr.Value( i );
3460     if ( !isdigit( c ) && c != ',' && c != '-' )
3461       aStr.SetValue( i, ',');
3462   }
3463   aStr.RemoveAll( ' ' );
3464
3465   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
3466   int i = 1;
3467   while ( tmpStr != "" )
3468   {
3469     tmpStr = aStr.Token( ",", i++ );
3470     int aPos = tmpStr.Search( '-' );
3471
3472     if ( aPos == -1 )
3473     {
3474       if ( tmpStr.IsIntegerValue() )
3475         myIds.Add( tmpStr.IntegerValue() );
3476       else
3477         return false;
3478     }
3479     else
3480     {
3481       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
3482       TCollection_AsciiString aMinStr = tmpStr;
3483
3484       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
3485       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
3486
3487       if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) ||
3488            (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
3489         return false;
3490
3491       myMin.push_back( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
3492       myMax.push_back( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
3493     }
3494   }
3495
3496   return true;
3497 }
3498
3499 //=======================================================================
3500 // name    : GetType
3501 // Purpose : Get type of supported entities
3502 //=======================================================================
3503 SMDSAbs_ElementType RangeOfIds::GetType() const
3504 {
3505   return myType;
3506 }
3507
3508 //=======================================================================
3509 // name    : SetType
3510 // Purpose : Set type of supported entities
3511 //=======================================================================
3512 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
3513 {
3514   myType = theType;
3515 }
3516
3517 //=======================================================================
3518 // name    : IsSatisfy
3519 // Purpose : Verify whether entity satisfies to this rpedicate
3520 //=======================================================================
3521 bool RangeOfIds::IsSatisfy( long theId )
3522 {
3523   if ( !myMesh )
3524     return false;
3525
3526   if ( myType == SMDSAbs_Node )
3527   {
3528     if ( myMesh->FindNode( theId ) == 0 )
3529       return false;
3530   }
3531   else
3532   {
3533     const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3534     if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All ))
3535       return false;
3536   }
3537
3538   if ( myIds.Contains( theId ) )
3539     return true;
3540
3541   for ( size_t i = 0; i < myMin.size(); i++ )
3542     if ( theId >= myMin[i] && theId <= myMax[i] )
3543       return true;
3544
3545   return false;
3546 }
3547
3548 /*
3549   Class       : Comparator
3550   Description : Base class for comparators
3551 */
3552 Comparator::Comparator():
3553   myMargin(0)
3554 {}
3555
3556 Comparator::~Comparator()
3557 {}
3558
3559 void Comparator::SetMesh( const SMDS_Mesh* theMesh )
3560 {
3561   if ( myFunctor )
3562     myFunctor->SetMesh( theMesh );
3563 }
3564
3565 void Comparator::SetMargin( double theValue )
3566 {
3567   myMargin = theValue;
3568 }
3569
3570 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
3571 {
3572   myFunctor = theFunct;
3573 }
3574
3575 SMDSAbs_ElementType Comparator::GetType() const
3576 {
3577   return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
3578 }
3579
3580 double Comparator::GetMargin()
3581 {
3582   return myMargin;
3583 }
3584
3585
3586 /*
3587   Class       : LessThan
3588   Description : Comparator "<"
3589 */
3590 bool LessThan::IsSatisfy( long theId )
3591 {
3592   return myFunctor && myFunctor->GetValue( theId ) < myMargin;
3593 }
3594
3595
3596 /*
3597   Class       : MoreThan
3598   Description : Comparator ">"
3599 */
3600 bool MoreThan::IsSatisfy( long theId )
3601 {
3602   return myFunctor && myFunctor->GetValue( theId ) > myMargin;
3603 }
3604
3605
3606 /*
3607   Class       : EqualTo
3608   Description : Comparator "="
3609 */
3610 EqualTo::EqualTo():
3611   myToler(Precision::Confusion())
3612 {}
3613
3614 bool EqualTo::IsSatisfy( long theId )
3615 {
3616   return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
3617 }
3618
3619 void EqualTo::SetTolerance( double theToler )
3620 {
3621   myToler = theToler;
3622 }
3623
3624 double EqualTo::GetTolerance()
3625 {
3626   return myToler;
3627 }
3628
3629 /*
3630   Class       : LogicalNOT
3631   Description : Logical NOT predicate
3632 */
3633 LogicalNOT::LogicalNOT()
3634 {}
3635
3636 LogicalNOT::~LogicalNOT()
3637 {}
3638
3639 bool LogicalNOT::IsSatisfy( long theId )
3640 {
3641   return myPredicate && !myPredicate->IsSatisfy( theId );
3642 }
3643
3644 void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
3645 {
3646   if ( myPredicate )
3647     myPredicate->SetMesh( theMesh );
3648 }
3649
3650 void LogicalNOT::SetPredicate( PredicatePtr thePred )
3651 {
3652   myPredicate = thePred;
3653 }
3654
3655 SMDSAbs_ElementType LogicalNOT::GetType() const
3656 {
3657   return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
3658 }
3659
3660
3661 /*
3662   Class       : LogicalBinary
3663   Description : Base class for binary logical predicate
3664 */
3665 LogicalBinary::LogicalBinary()
3666 {}
3667
3668 LogicalBinary::~LogicalBinary()
3669 {}
3670
3671 void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
3672 {
3673   if ( myPredicate1 )
3674     myPredicate1->SetMesh( theMesh );
3675
3676   if ( myPredicate2 )
3677     myPredicate2->SetMesh( theMesh );
3678 }
3679
3680 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
3681 {
3682   myPredicate1 = thePredicate;
3683 }
3684
3685 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
3686 {
3687   myPredicate2 = thePredicate;
3688 }
3689
3690 SMDSAbs_ElementType LogicalBinary::GetType() const
3691 {
3692   if ( !myPredicate1 || !myPredicate2 )
3693     return SMDSAbs_All;
3694
3695   SMDSAbs_ElementType aType1 = myPredicate1->GetType();
3696   SMDSAbs_ElementType aType2 = myPredicate2->GetType();
3697
3698   return aType1 == aType2 ? aType1 : SMDSAbs_All;
3699 }
3700
3701
3702 /*
3703   Class       : LogicalAND
3704   Description : Logical AND
3705 */
3706 bool LogicalAND::IsSatisfy( long theId )
3707 {
3708   return
3709     myPredicate1 &&
3710     myPredicate2 &&
3711     myPredicate1->IsSatisfy( theId ) &&
3712     myPredicate2->IsSatisfy( theId );
3713 }
3714
3715
3716 /*
3717   Class       : LogicalOR
3718   Description : Logical OR
3719 */
3720 bool LogicalOR::IsSatisfy( long theId )
3721 {
3722   return
3723     myPredicate1 &&
3724     myPredicate2 &&
3725     (myPredicate1->IsSatisfy( theId ) ||
3726     myPredicate2->IsSatisfy( theId ));
3727 }
3728
3729
3730 /*
3731                               FILTER
3732 */
3733
3734 // #ifdef WITH_TBB
3735 // #include <tbb/parallel_for.h>
3736 // #include <tbb/enumerable_thread_specific.h>
3737
3738 // namespace Parallel
3739 // {
3740 //   typedef tbb::enumerable_thread_specific< TIdSequence > TIdSeq;
3741
3742 //   struct Predicate
3743 //   {
3744 //     const SMDS_Mesh* myMesh;
3745 //     PredicatePtr     myPredicate;
3746 //     TIdSeq &         myOKIds;
3747 //     Predicate( const SMDS_Mesh* m, PredicatePtr p, TIdSeq & ids ):
3748 //       myMesh(m), myPredicate(p->Duplicate()), myOKIds(ids) {}
3749 //     void operator() ( const tbb::blocked_range<size_t>& r ) const
3750 //     {
3751 //       for ( size_t i = r.begin(); i != r.end(); ++i )
3752 //         if ( myPredicate->IsSatisfy( i ))
3753 //           myOKIds.local().push_back();
3754 //     }
3755 //   }
3756 // }
3757 // #endif
3758
3759 Filter::Filter()
3760 {}
3761
3762 Filter::~Filter()
3763 {}
3764
3765 void Filter::SetPredicate( PredicatePtr thePredicate )
3766 {
3767   myPredicate = thePredicate;
3768 }
3769
3770 void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
3771                             PredicatePtr         thePredicate,
3772                             TIdSequence&         theSequence,
3773                             SMDS_ElemIteratorPtr theElements )
3774 {
3775   theSequence.clear();
3776
3777   if ( !theMesh || !thePredicate )
3778     return;
3779
3780   thePredicate->SetMesh( theMesh );
3781
3782   if ( !theElements )
3783     theElements = theMesh->elementsIterator( thePredicate->GetType() );
3784
3785   if ( theElements ) {
3786     while ( theElements->more() ) {
3787       const SMDS_MeshElement* anElem = theElements->next();
3788       if ( thePredicate->GetType() == SMDSAbs_All ||
3789            thePredicate->GetType() == anElem->GetType() )
3790       {
3791         long anId = anElem->GetID();
3792         if ( thePredicate->IsSatisfy( anId ) )
3793           theSequence.push_back( anId );
3794       }
3795     }
3796   }
3797 }
3798
3799 void Filter::GetElementsId( const SMDS_Mesh*     theMesh,
3800                             Filter::TIdSequence& theSequence,
3801                             SMDS_ElemIteratorPtr theElements )
3802 {
3803   GetElementsId(theMesh,myPredicate,theSequence,theElements);
3804 }
3805
3806 /*
3807                               ManifoldPart
3808 */
3809
3810 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
3811
3812 /*
3813    Internal class Link
3814 */
3815
3816 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
3817                           SMDS_MeshNode* theNode2 )
3818 {
3819   myNode1 = theNode1;
3820   myNode2 = theNode2;
3821 }
3822
3823 ManifoldPart::Link::~Link()
3824 {
3825   myNode1 = 0;
3826   myNode2 = 0;
3827 }
3828
3829 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
3830 {
3831   if ( myNode1 == theLink.myNode1 &&
3832        myNode2 == theLink.myNode2 )
3833     return true;
3834   else if ( myNode1 == theLink.myNode2 &&
3835             myNode2 == theLink.myNode1 )
3836     return true;
3837   else
3838     return false;
3839 }
3840
3841 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
3842 {
3843   if(myNode1 < x.myNode1) return true;
3844   if(myNode1 == x.myNode1)
3845     if(myNode2 < x.myNode2) return true;
3846   return false;
3847 }
3848
3849 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
3850                             const ManifoldPart::Link& theLink2 )
3851 {
3852   return theLink1.IsEqual( theLink2 );
3853 }
3854
3855 ManifoldPart::ManifoldPart()
3856 {
3857   myMesh = 0;
3858   myAngToler = Precision::Angular();
3859   myIsOnlyManifold = true;
3860 }
3861
3862 ManifoldPart::~ManifoldPart()
3863 {
3864   myMesh = 0;
3865 }
3866
3867 void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
3868 {
3869   myMesh = theMesh;
3870   process();
3871 }
3872
3873 SMDSAbs_ElementType ManifoldPart::GetType() const
3874 { return SMDSAbs_Face; }
3875
3876 bool ManifoldPart::IsSatisfy( long theElementId )
3877 {
3878   return myMapIds.Contains( theElementId );
3879 }
3880
3881 void ManifoldPart::SetAngleTolerance( const double theAngToler )
3882 { myAngToler = theAngToler; }
3883
3884 double ManifoldPart::GetAngleTolerance() const
3885 { return myAngToler; }
3886
3887 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
3888 { myIsOnlyManifold = theIsOnly; }
3889
3890 void ManifoldPart::SetStartElem( const long  theStartId )
3891 { myStartElemId = theStartId; }
3892
3893 bool ManifoldPart::process()
3894 {
3895   myMapIds.Clear();
3896   myMapBadGeomIds.Clear();
3897
3898   myAllFacePtr.clear();
3899   myAllFacePtrIntDMap.clear();
3900   if ( !myMesh )
3901     return false;
3902
3903   // collect all faces into own map
3904   SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
3905   for (; anFaceItr->more(); )
3906   {
3907     SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
3908     myAllFacePtr.push_back( aFacePtr );
3909     myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
3910   }
3911
3912   SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
3913   if ( !aStartFace )
3914     return false;
3915
3916   // the map of non manifold links and bad geometry
3917   TMapOfLink aMapOfNonManifold;
3918   TIDsMap    aMapOfTreated;
3919
3920   // begin cycle on faces from start index and run on vector till the end
3921   //  and from begin to start index to cover whole vector
3922   const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
3923   bool isStartTreat = false;
3924   for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
3925   {
3926     if ( fi == aStartIndx )
3927       isStartTreat = true;
3928     // as result next time when fi will be equal to aStartIndx
3929
3930     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
3931     if ( aMapOfTreated.Contains( aFacePtr->GetID()) )
3932       continue;
3933
3934     aMapOfTreated.Add( aFacePtr->GetID() );
3935     TIDsMap aResFaces;
3936     if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
3937                          aMapOfNonManifold, aResFaces ) )
3938       continue;
3939     TIDsMap::Iterator anItr( aResFaces );
3940     for ( ; anItr.More(); anItr.Next() )
3941     {
3942       smIdType aFaceId = anItr.Key();
3943       aMapOfTreated.Add( aFaceId );
3944       myMapIds.Add( aFaceId );
3945     }
3946
3947     if ( fi == int( myAllFacePtr.size() - 1 ))
3948       fi = 0;
3949   } // end run on vector of faces
3950   return !myMapIds.IsEmpty();
3951 }
3952
3953 static void getLinks( const SMDS_MeshFace* theFace,
3954                       ManifoldPart::TVectorOfLink& theLinks )
3955 {
3956   int aNbNode = theFace->NbNodes();
3957   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
3958   int i = 1;
3959   SMDS_MeshNode* aNode = 0;
3960   for ( ; aNodeItr->more() && i <= aNbNode; )
3961   {
3962
3963     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
3964     if ( i == 1 )
3965       aNode = aN1;
3966     i++;
3967     SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
3968     i++;
3969     ManifoldPart::Link aLink( aN1, aN2 );
3970     theLinks.push_back( aLink );
3971   }
3972 }
3973
3974 bool ManifoldPart::findConnected
3975                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
3976                   SMDS_MeshFace*                           theStartFace,
3977                   ManifoldPart::TMapOfLink&                theNonManifold,
3978                   TIDsMap&                                 theResFaces )
3979 {
3980   theResFaces.Clear();
3981   if ( !theAllFacePtrInt.size() )
3982     return false;
3983
3984   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
3985   {
3986     myMapBadGeomIds.Add( theStartFace->GetID() );
3987     return false;
3988   }
3989
3990   ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
3991   ManifoldPart::TVectorOfLink aSeqOfBoundary;
3992   theResFaces.Add( theStartFace->GetID() );
3993   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
3994
3995   expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3996                  aDMapLinkFace, theNonManifold, theStartFace );
3997
3998   bool isDone = false;
3999   while ( !isDone && aMapOfBoundary.size() != 0 )
4000   {
4001     bool isToReset = false;
4002     ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
4003     for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
4004     {
4005       ManifoldPart::Link aLink = *pLink;
4006       if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
4007         continue;
4008       // each link could be treated only once
4009       aMapToSkip.insert( aLink );
4010
4011       ManifoldPart::TVectorOfFacePtr aFaces;
4012       // find next
4013       if ( myIsOnlyManifold &&
4014            (theNonManifold.find( aLink ) != theNonManifold.end()) )
4015         continue;
4016       else
4017       {
4018         getFacesByLink( aLink, aFaces );
4019         // filter the element to keep only indicated elements
4020         ManifoldPart::TVectorOfFacePtr aFiltered;
4021         ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
4022         for ( ; pFace != aFaces.end(); ++pFace )
4023         {
4024           SMDS_MeshFace* aFace = *pFace;
4025           if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
4026             aFiltered.push_back( aFace );
4027         }
4028         aFaces = aFiltered;
4029         if ( aFaces.size() < 2 )  // no neihgbour faces
4030           continue;
4031         else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
4032         {
4033           theNonManifold.insert( aLink );
4034           continue;
4035         }
4036       }
4037
4038       // compare normal with normals of neighbor element
4039       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
4040       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
4041       for ( ; pFace != aFaces.end(); ++pFace )
4042       {
4043         SMDS_MeshFace* aNextFace = *pFace;
4044         if ( aPrevFace == aNextFace )
4045           continue;
4046         smIdType anNextFaceID = aNextFace->GetID();
4047         if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
4048          // should not be with non manifold restriction. probably bad topology
4049           continue;
4050         // check if face was treated and skipped
4051         if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
4052              !isInPlane( aPrevFace, aNextFace ) )
4053           continue;
4054         // add new element to connected and extend the boundaries.
4055         theResFaces.Add( anNextFaceID );
4056         expandBoundary( aMapOfBoundary, aSeqOfBoundary,
4057                         aDMapLinkFace, theNonManifold, aNextFace );
4058         isToReset = true;
4059       }
4060     }
4061     isDone = !isToReset;
4062   }
4063
4064   return !theResFaces.IsEmpty();
4065 }
4066
4067 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
4068                               const SMDS_MeshFace* theFace2 )
4069 {
4070   gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
4071   gp_XYZ aNorm2XYZ = getNormale( theFace2 );
4072   if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
4073   {
4074     myMapBadGeomIds.Add( theFace2->GetID() );
4075     return false;
4076   }
4077   if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
4078     return true;
4079
4080   return false;
4081 }
4082
4083 void ManifoldPart::expandBoundary
4084                    ( ManifoldPart::TMapOfLink&            theMapOfBoundary,
4085                      ManifoldPart::TVectorOfLink&         theSeqOfBoundary,
4086                      ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
4087                      ManifoldPart::TMapOfLink&            theNonManifold,
4088                      SMDS_MeshFace*                       theNextFace ) const
4089 {
4090   ManifoldPart::TVectorOfLink aLinks;
4091   getLinks( theNextFace, aLinks );
4092   int aNbLink = (int)aLinks.size();
4093   for ( int i = 0; i < aNbLink; i++ )
4094   {
4095     ManifoldPart::Link aLink = aLinks[ i ];
4096     if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
4097       continue;
4098     if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
4099     {
4100       if ( myIsOnlyManifold )
4101       {
4102         // remove from boundary
4103         theMapOfBoundary.erase( aLink );
4104         ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
4105         for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
4106         {
4107           ManifoldPart::Link aBoundLink = *pLink;
4108           if ( aBoundLink.IsEqual( aLink ) )
4109           {
4110             theSeqOfBoundary.erase( pLink );
4111             break;
4112           }
4113         }
4114       }
4115     }
4116     else
4117     {
4118       theMapOfBoundary.insert( aLink );
4119       theSeqOfBoundary.push_back( aLink );
4120       theDMapLinkFacePtr[ aLink ] = theNextFace;
4121     }
4122   }
4123 }
4124
4125 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
4126                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
4127 {
4128
4129   // take all faces that shared first node
4130   SMDS_ElemIteratorPtr anItr = theLink.myNode1->GetInverseElementIterator( SMDSAbs_Face );
4131   SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > faces( anItr ), facesEnd;
4132   std::set<const SMDS_MeshElement *> aSetOfFaces( faces, facesEnd );
4133
4134   // take all faces that shared second node
4135   anItr = theLink.myNode2->GetInverseElementIterator( SMDSAbs_Face );
4136   // find the common part of two sets
4137   for ( ; anItr->more(); )
4138   {
4139     const SMDS_MeshElement* aFace = anItr->next();
4140     if ( aSetOfFaces.count( aFace ))
4141       theFaces.push_back( (SMDS_MeshFace*) aFace );
4142   }
4143 }
4144
4145 /*
4146   Class       : BelongToMeshGroup
4147   Description : Verify whether a mesh element is included into a mesh group
4148 */
4149 BelongToMeshGroup::BelongToMeshGroup(): myGroup( 0 )
4150 {
4151 }
4152
4153 void BelongToMeshGroup::SetGroup( SMESHDS_GroupBase* g )
4154 {
4155   myGroup = g;
4156 }
4157
4158 void BelongToMeshGroup::SetStoreName( const std::string& sn )
4159 {
4160   myStoreName = sn;
4161 }
4162
4163 void BelongToMeshGroup::SetMesh( const SMDS_Mesh* theMesh )
4164 {
4165   if ( myGroup && myGroup->GetMesh() != theMesh )
4166   {
4167     myGroup = 0;
4168   }
4169   if ( !myGroup && !myStoreName.empty() )
4170   {
4171     if ( const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh))
4172     {
4173       const std::set<SMESHDS_GroupBase*>& grps = aMesh->GetGroups();
4174       std::set<SMESHDS_GroupBase*>::const_iterator g = grps.begin();
4175       for ( ; g != grps.end() && !myGroup; ++g )
4176         if ( *g && myStoreName == (*g)->GetStoreName() )
4177           myGroup = *g;
4178     }
4179   }
4180   if ( myGroup )
4181   {
4182     myGroup->IsEmpty(); // make GroupOnFilter update its predicate
4183   }
4184 }
4185
4186 bool BelongToMeshGroup::IsSatisfy( long theElementId )
4187 {
4188   return myGroup ? myGroup->Contains( theElementId ) : false;
4189 }
4190
4191 SMDSAbs_ElementType BelongToMeshGroup::GetType() const
4192 {
4193   return myGroup ? myGroup->GetType() : SMDSAbs_All;
4194 }
4195
4196 //================================================================================
4197 //  ElementsOnSurface
4198 //================================================================================
4199
4200 ElementsOnSurface::ElementsOnSurface()
4201 {
4202   myIds.Clear();
4203   myType = SMDSAbs_All;
4204   mySurf.Nullify();
4205   myToler = Precision::Confusion();
4206   myUseBoundaries = false;
4207 }
4208
4209 ElementsOnSurface::~ElementsOnSurface()
4210 {
4211 }
4212
4213 void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
4214 {
4215   myMeshModifTracer.SetMesh( theMesh );
4216   if ( myMeshModifTracer.IsMeshModified())
4217     process();
4218 }
4219
4220 bool ElementsOnSurface::IsSatisfy( long theElementId )
4221 {
4222   return myIds.Contains( theElementId );
4223 }
4224
4225 SMDSAbs_ElementType ElementsOnSurface::GetType() const
4226 { return myType; }
4227
4228 void ElementsOnSurface::SetTolerance( const double theToler )
4229 {
4230   if ( myToler != theToler )
4231   {
4232     myToler = theToler;
4233     process();
4234   }
4235 }
4236
4237 double ElementsOnSurface::GetTolerance() const
4238 { return myToler; }
4239
4240 void ElementsOnSurface::SetUseBoundaries( bool theUse )
4241 {
4242   if ( myUseBoundaries != theUse ) {
4243     myUseBoundaries = theUse;
4244     SetSurface( mySurf, myType );
4245   }
4246 }
4247
4248 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
4249                                     const SMDSAbs_ElementType theType )
4250 {
4251   myIds.Clear();
4252   myType = theType;
4253   mySurf.Nullify();
4254   if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
4255     return;
4256   mySurf = TopoDS::Face( theShape );
4257   BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
4258   Standard_Real
4259     u1 = SA.FirstUParameter(),
4260     u2 = SA.LastUParameter(),
4261     v1 = SA.FirstVParameter(),
4262     v2 = SA.LastVParameter();
4263   Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
4264   myProjector.Init( surf, u1,u2, v1,v2 );
4265   process();
4266 }
4267
4268 void ElementsOnSurface::process()
4269 {
4270   myIds.Clear();
4271   if ( mySurf.IsNull() )
4272     return;
4273
4274   if ( !myMeshModifTracer.GetMesh() )
4275     return;
4276
4277   int nbElems = FromSmIdType<int>( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
4278   if ( nbElems > 0 )
4279     myIds.ReSize( nbElems );
4280
4281   SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType );
4282   for(; anIter->more(); )
4283     process( anIter->next() );
4284 }
4285
4286 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
4287 {
4288   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
4289   bool isSatisfy = true;
4290   for ( ; aNodeItr->more(); )
4291   {
4292     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
4293     if ( !isOnSurface( aNode ) )
4294     {
4295       isSatisfy = false;
4296       break;
4297     }
4298   }
4299   if ( isSatisfy )
4300     myIds.Add( theElemPtr->GetID() );
4301 }
4302
4303 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
4304 {
4305   if ( mySurf.IsNull() )
4306     return false;
4307
4308   gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
4309   //  double aToler2 = myToler * myToler;
4310 //   if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
4311 //   {
4312 //     gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
4313 //     if ( aPln.SquareDistance( aPnt ) > aToler2 )
4314 //       return false;
4315 //   }
4316 //   else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
4317 //   {
4318 //     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
4319 //     double aRad = aCyl.Radius();
4320 //     gp_Ax3 anAxis = aCyl.Position();
4321 //     gp_XYZ aLoc = aCyl.Location().XYZ();
4322 //     double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4323 //     double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4324 //     if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
4325 //       return false;
4326 //   }
4327 //   else
4328 //     return false;
4329   myProjector.Perform( aPnt );
4330   bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
4331
4332   return isOn;
4333 }
4334
4335
4336 //================================================================================
4337 //  ElementsOnShape
4338 //================================================================================
4339
4340 namespace {
4341   const int theIsCheckedFlag = 0x0000100;
4342 }
4343
4344 struct ElementsOnShape::Classifier
4345 {
4346   Classifier(): mySolidClfr(0), myProjFace(0), myProjEdge(0), myFlags(0) { myU = myV = 1e100; }
4347   ~Classifier();
4348   void Init(const TopoDS_Shape& s, double tol, const Bnd_B3d* box = 0 );
4349   bool IsOut(const gp_Pnt& p)        { return SetChecked( true ), (this->*myIsOutFun)( p ); }
4350   TopAbs_ShapeEnum ShapeType() const { return myShape.ShapeType(); }
4351   const TopoDS_Shape& Shape() const  { return myShape; }
4352   const Bnd_B3d* GetBndBox() const   { return & myBox; }
4353   double Tolerance() const           { return myTol; }
4354   bool IsChecked()                   { return myFlags & theIsCheckedFlag; }
4355   bool IsSetFlag( int flag ) const   { return myFlags & flag; }
4356   void SetChecked( bool is ) { is ? SetFlag( theIsCheckedFlag ) : UnsetFlag( theIsCheckedFlag ); }
4357   void SetFlag  ( int flag ) { myFlags |= flag; }
4358   void UnsetFlag( int flag ) { myFlags &= ~flag; }
4359   void GetParams( double & u, double & v ) const { u = myU; v = myV; }
4360
4361 private:
4362   bool isOutOfSolid (const gp_Pnt& p);
4363   bool isOutOfBox   (const gp_Pnt& p);
4364   bool isOutOfFace  (const gp_Pnt& p);
4365   bool isOutOfEdge  (const gp_Pnt& p);
4366   bool isOutOfVertex(const gp_Pnt& p);
4367   bool isOutOfNone  (const gp_Pnt& /*p*/) { return true; }
4368   bool isBox        (const TopoDS_Shape& s);
4369
4370   TopoDS_Shape prepareSolid( const TopoDS_Shape& theSolid );
4371
4372   bool (Classifier::*          myIsOutFun)(const gp_Pnt& p);
4373   BRepClass3d_SolidClassifier* mySolidClfr;
4374   Bnd_B3d                      myBox;
4375   GeomAPI_ProjectPointOnSurf*  myProjFace;
4376   GeomAPI_ProjectPointOnCurve* myProjEdge;
4377   gp_Pnt                       myVertexXYZ;
4378   TopoDS_Shape                 myShape;
4379   double                       myTol;
4380   double                       myU, myV; // result of isOutOfFace() and isOutOfEdge()
4381   int                          myFlags;
4382 };
4383
4384 struct ElementsOnShape::OctreeClassifier : public SMESH_Octree
4385 {
4386   OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers );
4387   OctreeClassifier( const OctreeClassifier*                           otherTree,
4388                     const std::vector< ElementsOnShape::Classifier >& clsOther,
4389                     std::vector< ElementsOnShape::Classifier >&       cls );
4390   void GetClassifiersAtPoint( const gp_XYZ& p,
4391                               std::vector< ElementsOnShape::Classifier* >& classifiers );
4392   size_t GetSize();
4393
4394 protected:
4395   OctreeClassifier() {}
4396   SMESH_Octree* newChild() const { return new OctreeClassifier; }
4397   void          buildChildrenData();
4398   Bnd_B3d*      buildRootBox();
4399
4400   std::vector< ElementsOnShape::Classifier* > myClassifiers;
4401 };
4402
4403
4404 ElementsOnShape::ElementsOnShape():
4405   myOctree(0),
4406   myType(SMDSAbs_All),
4407   myToler(Precision::Confusion()),
4408   myAllNodesFlag(false)
4409 {
4410 }
4411
4412 ElementsOnShape::~ElementsOnShape()
4413 {
4414   clearClassifiers();
4415 }
4416
4417 Predicate* ElementsOnShape::clone() const
4418 {
4419   size_t size = sizeof( *this );
4420   if ( myOctree )
4421     size += myOctree->GetSize();
4422   if ( !myClassifiers.empty() )
4423     size += sizeof( myClassifiers[0] ) * myClassifiers.size();
4424   if ( !myWorkClassifiers.empty() )
4425     size += sizeof( myWorkClassifiers[0] ) * myWorkClassifiers.size();
4426   if ( size > 1e+9 ) // 1G
4427   {
4428
4429   if (SALOME::VerbosityActivated())
4430     std::cout << "Avoid ElementsOnShape::clone(), too large: " << size << " bytes " << std::endl;
4431
4432     return 0;
4433   }
4434
4435   ElementsOnShape* cln = new ElementsOnShape();
4436   cln->SetAllNodes ( myAllNodesFlag );
4437   cln->SetTolerance( myToler );
4438   cln->SetMesh     ( myMeshModifTracer.GetMesh() );
4439   cln->myShape = myShape; // avoid creation of myClassifiers
4440   cln->SetShape    ( myShape, myType );
4441   cln->myClassifiers.resize( myClassifiers.size() );
4442   for ( size_t i = 0; i < myClassifiers.size(); ++i )
4443     cln->myClassifiers[ i ].Init( BRepBuilderAPI_Copy( myClassifiers[ i ].Shape()),
4444                                   myToler, myClassifiers[ i ].GetBndBox() );
4445   if ( myOctree ) // copy myOctree
4446   {
4447     cln->myOctree = new OctreeClassifier( myOctree, myClassifiers, cln->myClassifiers );
4448   }
4449   return cln;
4450 }
4451
4452 SMDSAbs_ElementType ElementsOnShape::GetType() const
4453 {
4454   return myType;
4455 }
4456
4457 void ElementsOnShape::SetTolerance (const double theToler)
4458 {
4459   if (myToler != theToler)
4460   {
4461     myToler = theToler;
4462     TopoDS_Shape s = myShape;
4463     myShape.Nullify();
4464     SetShape( s, myType );
4465   }
4466 }
4467
4468 double ElementsOnShape::GetTolerance() const
4469 {
4470   return myToler;
4471 }
4472
4473 void ElementsOnShape::SetAllNodes (bool theAllNodes)
4474 {
4475   myAllNodesFlag = theAllNodes;
4476 }
4477
4478 void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
4479 {
4480   myMeshModifTracer.SetMesh( theMesh );
4481   if ( myMeshModifTracer.IsMeshModified())
4482   {
4483     size_t nbNodes = theMesh ? theMesh->NbNodes() : 0;
4484     if ( myNodeIsChecked.size() == nbNodes )
4485     {
4486       std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4487     }
4488     else
4489     {
4490       SMESHUtils::FreeVector( myNodeIsChecked );
4491       SMESHUtils::FreeVector( myNodeIsOut );
4492       myNodeIsChecked.resize( nbNodes, false );
4493       myNodeIsOut.resize( nbNodes );
4494     }
4495   }
4496 }
4497
4498 bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut )
4499 {
4500   if ( n->GetID() >= (int) myNodeIsChecked.size() ||
4501        !myNodeIsChecked[ n->GetID() ])
4502     return false;
4503
4504   isOut = myNodeIsOut[ n->GetID() ];
4505   return true;
4506 }
4507
4508 void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool  isOut )
4509 {
4510   if ( n->GetID() < (int) myNodeIsChecked.size() )
4511   {
4512     myNodeIsChecked[ n->GetID() ] = true;
4513     myNodeIsOut    [ n->GetID() ] = isOut;
4514   }
4515 }
4516
4517 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
4518                                 const SMDSAbs_ElementType theType)
4519 {
4520   bool shapeChanges = ( myShape != theShape );
4521   myType  = theType;
4522   myShape = theShape;
4523   if ( myShape.IsNull() ) return;
4524
4525   if ( shapeChanges )
4526   {
4527     // find most complex shapes
4528     TopTools_IndexedMapOfShape shapesMap;
4529     TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
4530     TopExp_Explorer sub;
4531     for ( int i = 0; i < 4; ++i )
4532     {
4533       if ( shapesMap.IsEmpty() )
4534         for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() )
4535           shapesMap.Add( sub.Current() );
4536       if ( i > 0 )
4537         for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() )
4538           shapesMap.Add( sub.Current() );
4539     }
4540
4541     clearClassifiers();
4542     myClassifiers.resize( shapesMap.Extent() );
4543     for ( int i = 0; i < shapesMap.Extent(); ++i )
4544       myClassifiers[ i ].Init( shapesMap( i+1 ), myToler );
4545   }
4546
4547   if ( theType == SMDSAbs_Node )
4548   {
4549     SMESHUtils::FreeVector( myNodeIsChecked );
4550     SMESHUtils::FreeVector( myNodeIsOut );
4551   }
4552   else
4553   {
4554     std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4555   }
4556 }
4557
4558 void ElementsOnShape::clearClassifiers()
4559 {
4560   // for ( size_t i = 0; i < myClassifiers.size(); ++i )
4561   //   delete myClassifiers[ i ];
4562   myClassifiers.clear();
4563
4564   delete myOctree;
4565   myOctree = 0;
4566 }
4567
4568 bool ElementsOnShape::IsSatisfy( long elemId )
4569 {
4570   if ( myClassifiers.empty() )
4571     return false;
4572
4573   const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh();
4574   if ( myType == SMDSAbs_Node )
4575     return IsSatisfy( mesh->FindNode( elemId ));
4576   return IsSatisfy( mesh->FindElement( elemId ));
4577 }
4578
4579 bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem)
4580 {
4581   if ( !elem )
4582     return false;
4583
4584   bool isSatisfy = myAllNodesFlag, isNodeOut;
4585
4586   gp_XYZ centerXYZ (0, 0, 0);
4587
4588   if ( !myOctree && myClassifiers.size() > 5 )
4589   {
4590     myWorkClassifiers.resize( myClassifiers.size() );
4591     for ( size_t i = 0; i < myClassifiers.size(); ++i )
4592       myWorkClassifiers[ i ] = & myClassifiers[ i ];
4593     myOctree = new OctreeClassifier( myWorkClassifiers );
4594
4595     SMESHUtils::FreeVector( myWorkClassifiers );
4596   }
4597
4598   for ( int i = 0, nb = elem->NbNodes(); i < nb  && (isSatisfy == myAllNodesFlag); ++i )
4599   {
4600     SMESH_TNodeXYZ aPnt( elem->GetNode( i ));
4601     centerXYZ += aPnt;
4602
4603     isNodeOut = true;
4604     if ( !getNodeIsOut( aPnt._node, isNodeOut ))
4605     {
4606       if ( myOctree )
4607       {
4608         myWorkClassifiers.clear();
4609         myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers );
4610
4611         for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
4612           myWorkClassifiers[i]->SetChecked( false );
4613
4614         for ( size_t i = 0; i < myWorkClassifiers.size() && isNodeOut; ++i )
4615           if ( !myWorkClassifiers[i]->IsChecked() )
4616             isNodeOut = myWorkClassifiers[i]->IsOut( aPnt );
4617       }
4618       else
4619       {
4620         for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i )
4621           isNodeOut = myClassifiers[i].IsOut( aPnt );
4622       }
4623       setNodeIsOut( aPnt._node, isNodeOut );
4624     }
4625     isSatisfy = !isNodeOut;
4626   }
4627
4628   // Check the center point for volumes MantisBug 0020168
4629   if ( isSatisfy &&
4630        myAllNodesFlag &&
4631        myClassifiers[0].ShapeType() == TopAbs_SOLID )
4632   {
4633     centerXYZ /= elem->NbNodes();
4634     isSatisfy = false;
4635     if ( myOctree )
4636     {
4637       myWorkClassifiers.clear();
4638       myOctree->GetClassifiersAtPoint( centerXYZ, myWorkClassifiers );
4639       for ( size_t i = 0; i < myWorkClassifiers.size() && !isSatisfy; ++i )
4640         isSatisfy = ! myWorkClassifiers[i]->IsOut( centerXYZ );
4641     }
4642     else
4643     {
4644       for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i )
4645         isSatisfy = ! myClassifiers[i].IsOut( centerXYZ );
4646     }
4647   }
4648
4649   return isSatisfy;
4650 }
4651
4652 //================================================================================
4653 /*!
4654  * \brief Check and optionally return a satisfying shape
4655  */
4656 //================================================================================
4657
4658 bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node,
4659                                  TopoDS_Shape*        okShape)
4660 {
4661   if ( !node )
4662     return false;
4663
4664   if ( !myOctree && myClassifiers.size() > 5 )
4665   {
4666     myWorkClassifiers.resize( myClassifiers.size() );
4667     for ( size_t i = 0; i < myClassifiers.size(); ++i )
4668       myWorkClassifiers[ i ] = & myClassifiers[ i ];
4669     myOctree = new OctreeClassifier( myWorkClassifiers );
4670   }
4671
4672   bool isNodeOut = true;
4673
4674   if ( okShape || !getNodeIsOut( node, isNodeOut ))
4675   {
4676     SMESH_NodeXYZ aPnt = node;
4677     if ( myOctree )
4678     {
4679       myWorkClassifiers.clear();
4680       myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers );
4681
4682       for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
4683         myWorkClassifiers[i]->SetChecked( false );
4684
4685       for ( size_t i = 0; i < myWorkClassifiers.size(); ++i )
4686         if ( !myWorkClassifiers[i]->IsChecked() &&
4687              !myWorkClassifiers[i]->IsOut( aPnt ))
4688         {
4689           isNodeOut = false;
4690           if ( okShape )
4691             *okShape = myWorkClassifiers[i]->Shape();
4692           myWorkClassifiers[i]->GetParams( myU, myV );
4693           break;
4694         }
4695     }
4696     else
4697     {
4698       for ( size_t i = 0; i < myClassifiers.size(); ++i )
4699         if ( !myClassifiers[i].IsOut( aPnt ))
4700         {
4701           isNodeOut = false;
4702           if ( okShape )
4703             *okShape = myClassifiers[i].Shape();
4704           myClassifiers[i].GetParams( myU, myV );
4705           break;
4706         }
4707     }
4708     setNodeIsOut( node, isNodeOut );
4709   }
4710
4711   return !isNodeOut;
4712 }
4713
4714 void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape,
4715                                         double              theTol,
4716                                         const Bnd_B3d*      theBox )
4717 {
4718   myShape = theShape;
4719   myTol   = theTol;
4720   myFlags = 0;
4721
4722   bool isShapeBox = false;
4723   switch ( myShape.ShapeType() )
4724   {
4725   case TopAbs_SOLID:
4726   {
4727     if (( isShapeBox = isBox( theShape )))
4728     {
4729       myIsOutFun = & ElementsOnShape::Classifier::isOutOfBox;
4730     }
4731     else
4732     {
4733       mySolidClfr = new BRepClass3d_SolidClassifier( prepareSolid( theShape ));
4734       myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid;
4735     }
4736     break;
4737   }
4738   case TopAbs_FACE:
4739   {
4740     Standard_Real u1,u2,v1,v2;
4741     Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
4742     if ( surf.IsNull() )
4743       myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
4744     else
4745     {
4746       surf->Bounds( u1,u2,v1,v2 );
4747       myProjFace = new GeomAPI_ProjectPointOnSurf;
4748       myProjFace->Init( surf, u1,u2, v1,v2, myTol );
4749       myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace;
4750     }
4751     break;
4752   }
4753   case TopAbs_EDGE:
4754   {
4755     Standard_Real u1, u2;
4756     Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2);
4757     if ( curve.IsNull() )
4758       myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
4759     else
4760     {
4761       myProjEdge = new GeomAPI_ProjectPointOnCurve;
4762       myProjEdge->Init( curve, u1, u2 );
4763       myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
4764     }
4765     break;
4766   }
4767   case TopAbs_VERTEX:
4768   {
4769     myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
4770     myIsOutFun = & ElementsOnShape::Classifier::isOutOfVertex;
4771     break;
4772   }
4773   default:
4774     throw SALOME_Exception("Programmer error in usage of ElementsOnShape::Classifier");
4775   }
4776
4777   if ( !isShapeBox )
4778   {
4779     if ( theBox )
4780     {
4781       myBox = *theBox;
4782     }
4783     else
4784     {
4785       Bnd_Box box;
4786       if ( myShape.ShapeType() == TopAbs_FACE )
4787       {
4788         BRepAdaptor_Surface SA( TopoDS::Face( myShape ), /*useBoundaries=*/false );
4789         if ( SA.GetType() == GeomAbs_BSplineSurface )
4790           BRepBndLib::AddOptimal( myShape, box,
4791                                   /*useTriangulation=*/true, /*useShapeTolerance=*/true );
4792       }
4793       if ( box.IsVoid() )
4794         BRepBndLib::Add( myShape, box );
4795       myBox.Clear();
4796       myBox.Add( box.CornerMin() );
4797       myBox.Add( box.CornerMax() );
4798       gp_XYZ halfSize = 0.5 * ( box.CornerMax().XYZ() - box.CornerMin().XYZ() );
4799       for ( int iDim = 1; iDim <= 3; ++iDim )
4800       {
4801         double x = halfSize.Coord( iDim );
4802         halfSize.SetCoord( iDim, x + Max( myTol, 1e-2 * x ));
4803       }
4804       myBox.SetHSize( halfSize );
4805     }
4806   }
4807 }
4808
4809 ElementsOnShape::Classifier::~Classifier()
4810 {
4811   delete mySolidClfr; mySolidClfr = 0;
4812   delete myProjFace;  myProjFace = 0;
4813   delete myProjEdge;  myProjEdge = 0;
4814 }
4815
4816 TopoDS_Shape ElementsOnShape::Classifier::prepareSolid( const TopoDS_Shape& theSolid )
4817 {
4818   // try to limit tolerance of theSolid down to myTol (issue #19026)
4819
4820   // check if tolerance of theSolid is more than myTol
4821   bool tolIsOk = true; // max tolerance is at VERTEXes
4822   for ( TopExp_Explorer exp( theSolid, TopAbs_VERTEX ); exp.More() &&  tolIsOk; exp.Next() )
4823     tolIsOk = ( myTol >= BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current() )));
4824   if ( tolIsOk )
4825     return theSolid;
4826
4827   // make a copy to prevent the original shape from changes
4828   TopoDS_Shape resultShape = BRepBuilderAPI_Copy( theSolid );
4829
4830   if ( !GEOMUtils::FixShapeTolerance( resultShape, TopAbs_SHAPE, myTol ))
4831     return theSolid;
4832   return resultShape;
4833 }
4834
4835 bool ElementsOnShape::Classifier::isOutOfSolid( const gp_Pnt& p )
4836 {
4837   if ( isOutOfBox( p )) return true;
4838   mySolidClfr->Perform( p, myTol );
4839   return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON );
4840 }
4841
4842 bool ElementsOnShape::Classifier::isOutOfBox( const gp_Pnt& p )
4843 {
4844   return myBox.IsOut( p.XYZ() );
4845 }
4846
4847 bool ElementsOnShape::Classifier::isOutOfFace( const gp_Pnt& p )
4848 {
4849   if ( isOutOfBox( p )) return true;
4850   myProjFace->Perform( p );
4851   if ( myProjFace->IsDone() && myProjFace->LowerDistance() <= myTol )
4852   {
4853     // check relatively to the face
4854     myProjFace->LowerDistanceParameters( myU, myV );
4855     gp_Pnt2d aProjPnt( myU, myV );
4856     BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
4857     if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
4858       return false;
4859   }
4860   return true;
4861 }
4862
4863 bool ElementsOnShape::Classifier::isOutOfEdge( const gp_Pnt& p )
4864 {
4865   if ( isOutOfBox( p )) return true;
4866   myProjEdge->Perform( p );
4867   bool isOn = ( myProjEdge->NbPoints() > 0 && myProjEdge->LowerDistance() <= myTol );
4868   if ( isOn )
4869     myU = myProjEdge->LowerDistanceParameter();
4870   return !isOn;
4871 }
4872
4873 bool ElementsOnShape::Classifier::isOutOfVertex( const gp_Pnt& p )
4874 {
4875   return ( myVertexXYZ.Distance( p ) > myTol );
4876 }
4877
4878 bool ElementsOnShape::Classifier::isBox(const TopoDS_Shape& theShape )
4879 {
4880   TopTools_IndexedMapOfShape vMap;
4881   TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
4882   if ( vMap.Extent() != 8 )
4883     return false;
4884
4885   myBox.Clear();
4886   for ( int i = 1; i <= 8; ++i )
4887     myBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))).XYZ() );
4888
4889   gp_XYZ pMin = myBox.CornerMin(), pMax = myBox.CornerMax();
4890   for ( int i = 1; i <= 8; ++i )
4891   {
4892     gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i )));
4893     for ( int iC = 1; iC <= 3; ++ iC )
4894     {
4895       double d1 = Abs( pMin.Coord( iC ) - p.Coord( iC ));
4896       double d2 = Abs( pMax.Coord( iC ) - p.Coord( iC ));
4897       if ( Min( d1, d2 ) > myTol )
4898         return false;
4899     }
4900   }
4901   myBox.Enlarge( myTol );
4902   return true;
4903 }
4904
4905 ElementsOnShape::
4906 OctreeClassifier::OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers )
4907   :SMESH_Octree( new SMESH_TreeLimit )
4908 {
4909   myClassifiers = classifiers;
4910   compute();
4911 }
4912
4913 ElementsOnShape::
4914 OctreeClassifier::OctreeClassifier( const OctreeClassifier*                           otherTree,
4915                                     const std::vector< ElementsOnShape::Classifier >& clsOther,
4916                                     std::vector< ElementsOnShape::Classifier >&       cls )
4917   :SMESH_Octree( new SMESH_TreeLimit )
4918 {
4919   myBox = new Bnd_B3d( *otherTree->getBox() );
4920
4921   if (( myIsLeaf = otherTree->isLeaf() ))
4922   {
4923     myClassifiers.resize( otherTree->myClassifiers.size() );
4924     for ( size_t i = 0; i < otherTree->myClassifiers.size(); ++i )
4925     {
4926       int ind = otherTree->myClassifiers[i] - & clsOther[0];
4927       myClassifiers[ i ] = & cls[ ind ];
4928     }
4929   }
4930   else if ( otherTree->myChildren )
4931   {
4932     myChildren = new SMESH_Tree< Bnd_B3d, 8 > * [ 8 ];
4933     for ( int i = 0; i < nbChildren(); i++ )
4934       myChildren[i] =
4935         new OctreeClassifier( static_cast<const OctreeClassifier*>( otherTree->myChildren[i]),
4936                               clsOther, cls );
4937   }
4938 }
4939
4940 void ElementsOnShape::
4941 OctreeClassifier::GetClassifiersAtPoint( const gp_XYZ& point,
4942                                          std::vector< ElementsOnShape::Classifier* >& result )
4943 {
4944   if ( getBox()->IsOut( point ))
4945     return;
4946
4947   if ( isLeaf() )
4948   {
4949     for ( size_t i = 0; i < myClassifiers.size(); ++i )
4950       if ( !myClassifiers[i]->GetBndBox()->IsOut( point ))
4951         result.push_back( myClassifiers[i] );
4952   }
4953   else
4954   {
4955     for (int i = 0; i < nbChildren(); i++)
4956       ((OctreeClassifier*) myChildren[i])->GetClassifiersAtPoint( point, result );
4957   }
4958 }
4959
4960 size_t ElementsOnShape::OctreeClassifier::GetSize()
4961 {
4962   size_t res = sizeof( *this );
4963   if ( !myClassifiers.empty() )
4964     res += sizeof( myClassifiers[0] ) * myClassifiers.size();
4965
4966   if ( !isLeaf() )
4967     for (int i = 0; i < nbChildren(); i++)
4968       res += ((OctreeClassifier*) myChildren[i])->GetSize();
4969
4970   return res;
4971 }
4972
4973 void ElementsOnShape::OctreeClassifier::buildChildrenData()
4974 {
4975   // distribute myClassifiers among myChildren
4976
4977   const int childFlag[8] = { 0x0000001,
4978                              0x0000002,
4979                              0x0000004,
4980                              0x0000008,
4981                              0x0000010,
4982                              0x0000020,
4983                              0x0000040,
4984                              0x0000080 };
4985   int nbInChild[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
4986
4987   for ( size_t i = 0; i < myClassifiers.size(); ++i )
4988   {
4989     for ( int j = 0; j < nbChildren(); j++ )
4990     {
4991       if ( !myClassifiers[i]->GetBndBox()->IsOut( *myChildren[j]->getBox() ))
4992       {
4993         myClassifiers[i]->SetFlag( childFlag[ j ]);
4994         ++nbInChild[ j ];
4995       }
4996     }
4997   }
4998
4999   for ( int j = 0; j < nbChildren(); j++ )
5000   {
5001     OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ j ]);
5002     child->myClassifiers.resize( nbInChild[ j ]);
5003     for ( size_t i = 0; nbInChild[ j ] && i < myClassifiers.size(); ++i )
5004     {
5005       if ( myClassifiers[ i ]->IsSetFlag( childFlag[ j ]))
5006       {
5007         --nbInChild[ j ];
5008         child->myClassifiers[ nbInChild[ j ]] = myClassifiers[ i ];
5009         myClassifiers[ i ]->UnsetFlag( childFlag[ j ]);
5010       }
5011     }
5012   }
5013   SMESHUtils::FreeVector( myClassifiers );
5014
5015   // define if a child isLeaf()
5016   for ( int i = 0; i < nbChildren(); i++ )
5017   {
5018     OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ i ]);
5019     child->myIsLeaf = ( child->myClassifiers.size() <= 5  ||
5020                         child->maxSize() < child->myClassifiers[0]->Tolerance() );
5021   }
5022 }
5023
5024 Bnd_B3d* ElementsOnShape::OctreeClassifier::buildRootBox()
5025 {
5026   Bnd_B3d* box = new Bnd_B3d;
5027   for ( size_t i = 0; i < myClassifiers.size(); ++i )
5028     box->Add( *myClassifiers[i]->GetBndBox() );
5029   return box;
5030 }
5031
5032 /*
5033   Class       : BelongToGeom
5034   Description : Predicate for verifying whether entity belongs to
5035                 specified geometrical support
5036 */
5037
5038 BelongToGeom::BelongToGeom()
5039   : myMeshDS(NULL),
5040     myType(SMDSAbs_NbElementTypes),
5041     myIsSubshape(false),
5042     myTolerance(Precision::Confusion())
5043 {}
5044
5045 Predicate* BelongToGeom::clone() const
5046 {
5047   BelongToGeom* cln = 0;
5048   if ( myElementsOnShapePtr )
5049     if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
5050     {
5051       cln = new BelongToGeom( *this );
5052       cln->myElementsOnShapePtr.reset( eos );
5053     }
5054   return cln;
5055 }
5056
5057 void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
5058 {
5059   if ( myMeshDS != theMesh )
5060   {
5061     myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
5062     init();
5063   }
5064   if ( myElementsOnShapePtr )
5065     myElementsOnShapePtr->SetMesh( myMeshDS );
5066 }
5067
5068 void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
5069 {
5070   if ( myShape != theShape )
5071   {
5072     myShape = theShape;
5073     init();
5074   }
5075 }
5076
5077 static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
5078                         const TopoDS_Shape&               theShape)
5079 {
5080   if (theMap.Contains(theShape)) return true;
5081
5082   if (theShape.ShapeType() == TopAbs_COMPOUND ||
5083       theShape.ShapeType() == TopAbs_COMPSOLID)
5084   {
5085     TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
5086     for (; anIt.More(); anIt.Next())
5087     {
5088       if (!IsSubShape(theMap, anIt.Value())) {
5089         return false;
5090       }
5091     }
5092     return true;
5093   }
5094
5095   return false;
5096 }
5097
5098 void BelongToGeom::init()
5099 {
5100   if ( !myMeshDS || myShape.IsNull() ) return;
5101
5102   // is sub-shape of main shape?
5103   TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
5104   if (aMainShape.IsNull()) {
5105     myIsSubshape = false;
5106   }
5107   else {
5108     TopTools_IndexedMapOfShape aMap;
5109     TopExp::MapShapes( aMainShape, aMap );
5110     myIsSubshape = IsSubShape( aMap, myShape );
5111     if ( myIsSubshape )
5112     {
5113       aMap.Clear();
5114       TopExp::MapShapes( myShape, aMap );
5115       mySubShapesIDs.Clear();
5116       for ( int i = 1; i <= aMap.Extent(); ++i )
5117       {
5118         int subID = myMeshDS->ShapeToIndex( aMap( i ));
5119         if ( subID > 0 )
5120           mySubShapesIDs.Add( subID );
5121       }
5122     }
5123   }
5124
5125   //if (!myIsSubshape) // to be always ready to check an element not bound to geometry
5126   {
5127     if ( !myElementsOnShapePtr )
5128       myElementsOnShapePtr.reset( new ElementsOnShape() );
5129     myElementsOnShapePtr->SetTolerance( myTolerance );
5130     myElementsOnShapePtr->SetAllNodes( true ); // "belong", while false means "lays on"
5131     myElementsOnShapePtr->SetMesh( myMeshDS );
5132     myElementsOnShapePtr->SetShape( myShape, myType );
5133   }
5134 }
5135
5136 bool BelongToGeom::IsSatisfy (long theId)
5137 {
5138   if (myMeshDS == 0 || myShape.IsNull())
5139     return false;
5140
5141   if (!myIsSubshape)
5142   {
5143     return myElementsOnShapePtr->IsSatisfy(theId);
5144   }
5145
5146   // Case of sub-mesh
5147
5148   if (myType == SMDSAbs_Node)
5149   {
5150     if ( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ))
5151     {
5152       if ( aNode->getshapeId() < 1 )
5153         return myElementsOnShapePtr->IsSatisfy(theId);
5154       else
5155         return mySubShapesIDs.Contains( aNode->getshapeId() );
5156     }
5157   }
5158   else
5159   {
5160     if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
5161     {
5162       if ( myType == SMDSAbs_All || anElem->GetType() == myType )
5163       {
5164         if ( anElem->getshapeId() < 1 )
5165           return myElementsOnShapePtr->IsSatisfy(theId);
5166         else
5167           return mySubShapesIDs.Contains( anElem->getshapeId() );
5168       }
5169     }
5170   }
5171
5172   return false;
5173 }
5174
5175 void BelongToGeom::SetType (SMDSAbs_ElementType theType)
5176 {
5177   if ( myType != theType )
5178   {
5179     myType = theType;
5180     init();
5181   }
5182 }
5183
5184 SMDSAbs_ElementType BelongToGeom::GetType() const
5185 {
5186   return myType;
5187 }
5188
5189 TopoDS_Shape BelongToGeom::GetShape()
5190 {
5191   return myShape;
5192 }
5193
5194 const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const
5195 {
5196   return myMeshDS;
5197 }
5198
5199 void BelongToGeom::SetTolerance (double theTolerance)
5200 {
5201   myTolerance = theTolerance;
5202   init();
5203 }
5204
5205 double BelongToGeom::GetTolerance()
5206 {
5207   return myTolerance;
5208 }
5209
5210 /*
5211   Class       : LyingOnGeom
5212   Description : Predicate for verifying whether entiy lying or partially lying on
5213   specified geometrical support
5214 */
5215
5216 LyingOnGeom::LyingOnGeom()
5217   : myMeshDS(NULL),
5218     myType(SMDSAbs_NbElementTypes),
5219     myIsSubshape(false),
5220     myTolerance(Precision::Confusion())
5221 {}
5222
5223 Predicate* LyingOnGeom::clone() const
5224 {
5225   LyingOnGeom* cln = 0;
5226   if ( myElementsOnShapePtr )
5227     if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
5228     {
5229       cln = new LyingOnGeom( *this );
5230       cln->myElementsOnShapePtr.reset( eos );
5231     }
5232   return cln;
5233 }
5234
5235 void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
5236 {
5237   if ( myMeshDS != theMesh )
5238   {
5239     myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
5240     init();
5241   }
5242   if ( myElementsOnShapePtr )
5243     myElementsOnShapePtr->SetMesh( myMeshDS );
5244 }
5245
5246 void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
5247 {
5248   if ( myShape != theShape )
5249   {
5250     myShape = theShape;
5251     init();
5252   }
5253 }
5254
5255 void LyingOnGeom::init()
5256 {
5257   if (!myMeshDS || myShape.IsNull()) return;
5258
5259   // is sub-shape of main shape?
5260   TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
5261   if (aMainShape.IsNull()) {
5262     myIsSubshape = false;
5263   }
5264   else {
5265     myIsSubshape = myMeshDS->IsGroupOfSubShapes( myShape );
5266   }
5267
5268   if (myIsSubshape)
5269   {
5270     TopTools_IndexedMapOfShape shapes;
5271     TopExp::MapShapes( myShape, shapes );
5272     mySubShapesIDs.Clear();
5273     for ( int i = 1; i <= shapes.Extent(); ++i )
5274     {
5275       int subID = myMeshDS->ShapeToIndex( shapes( i ));
5276       if ( subID > 0 )
5277         mySubShapesIDs.Add( subID );
5278     }
5279   }
5280   // else // to be always ready to check an element not bound to geometry
5281   {
5282     if ( !myElementsOnShapePtr )
5283       myElementsOnShapePtr.reset( new ElementsOnShape() );
5284     myElementsOnShapePtr->SetTolerance( myTolerance );
5285     myElementsOnShapePtr->SetAllNodes( false ); // lays on, while true means "belong"
5286     myElementsOnShapePtr->SetMesh( myMeshDS );
5287     myElementsOnShapePtr->SetShape( myShape, myType );
5288   }
5289 }
5290
5291 bool LyingOnGeom::IsSatisfy( long theId )
5292 {
5293   if ( myMeshDS == 0 || myShape.IsNull() )
5294     return false;
5295
5296   if (!myIsSubshape)
5297   {
5298     return myElementsOnShapePtr->IsSatisfy(theId);
5299   }
5300
5301   // Case of sub-mesh
5302
5303   const SMDS_MeshElement* elem =
5304     ( myType == SMDSAbs_Node ) ? myMeshDS->FindNode( theId ) : myMeshDS->FindElement( theId );
5305
5306   if ( mySubShapesIDs.Contains( elem->getshapeId() ))
5307     return true;
5308
5309   if (( elem->GetType() != SMDSAbs_Node ) &&
5310       ( myType == SMDSAbs_All || elem->GetType() == myType ))
5311   {
5312     SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator();
5313     while ( nodeItr->more() )
5314     {
5315       const SMDS_MeshElement* aNode = nodeItr->next();
5316       if ( mySubShapesIDs.Contains( aNode->getshapeId() ))
5317         return true;
5318     }
5319   }
5320
5321   return false;
5322 }
5323
5324 void LyingOnGeom::SetType( SMDSAbs_ElementType theType )
5325 {
5326   if ( myType != theType )
5327   {
5328     myType = theType;
5329     init();
5330   }
5331 }
5332
5333 SMDSAbs_ElementType LyingOnGeom::GetType() const
5334 {
5335   return myType;
5336 }
5337
5338 TopoDS_Shape LyingOnGeom::GetShape()
5339 {
5340   return myShape;
5341 }
5342
5343 const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const
5344 {
5345   return myMeshDS;
5346 }
5347
5348 void LyingOnGeom::SetTolerance (double theTolerance)
5349 {
5350   myTolerance = theTolerance;
5351   init();
5352 }
5353
5354 double LyingOnGeom::GetTolerance()
5355 {
5356   return myTolerance;
5357 }
5358
5359 TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0)
5360 {}
5361
5362 TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n), myElem(0)
5363 {}
5364
5365 TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t), myElem(0)
5366 {}
5367
5368 TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray), myElem(theSequenceOfXYZ.myElem)
5369 {}
5370
5371 template <class InputIterator>
5372 TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd), myElem(0)
5373 {}
5374
5375 TSequenceOfXYZ::~TSequenceOfXYZ()
5376 {}
5377
5378 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
5379 {
5380   myArray = theSequenceOfXYZ.myArray;
5381   myElem  = theSequenceOfXYZ.myElem;
5382   return *this;
5383 }
5384
5385 gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
5386 {
5387   return myArray[n-1];
5388 }
5389
5390 const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
5391 {
5392   return myArray[n-1];
5393 }
5394
5395 void TSequenceOfXYZ::clear()
5396 {
5397   myArray.clear();
5398 }
5399
5400 void TSequenceOfXYZ::reserve(size_type n)
5401 {
5402   myArray.reserve(n);
5403 }
5404
5405 void TSequenceOfXYZ::push_back(const gp_XYZ& v)
5406 {
5407   myArray.push_back(v);
5408 }
5409
5410 TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
5411 {
5412   return myArray.size();
5413 }
5414
5415 SMDSAbs_EntityType TSequenceOfXYZ::getElementEntity() const
5416 {
5417   return myElem ? myElem->GetEntityType() : SMDSEntity_Last;
5418 }
5419
5420 TMeshModifTracer::TMeshModifTracer():
5421   myMeshModifTime(0), myMesh(0)
5422 {
5423 }
5424 void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh )
5425 {
5426   if ( theMesh != myMesh )
5427     myMeshModifTime = 0;
5428   myMesh = theMesh;
5429 }
5430 bool TMeshModifTracer::IsMeshModified()
5431 {
5432   bool modified = false;
5433   if ( myMesh )
5434   {
5435     modified = ( myMeshModifTime != myMesh->GetMTime() );
5436     myMeshModifTime = myMesh->GetMTime();
5437   }
5438   return modified;
5439 }