src/Controls/SMESH_Controls.cxx

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