src/Controls/SMESH_Controls.cxx

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