src/Controls/SMESH_Controls.cxx

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