src/Controls/SMESH_Controls.cxx

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