src/Controls/SMESH_Controls.cxx

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