src/Controls/SMESH_Controls.cxx

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