src/Controls/SMESH_Controls.cxx

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