src/Controls/SMESH_Controls.cxx

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