src/Controls/SMESH_Controls.cxx

   1 //  Copyright (C) 2007-2010  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.
  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 <set>
  26
  27 #include <BRepAdaptor_Surface.hxx>
  28 #include <BRepClass_FaceClassifier.hxx>
  29 #include <BRep_Tool.hxx>
  30
  31 #include <TopAbs.hxx>
  32 #include <TopoDS.hxx>
  33 #include <TopoDS_Edge.hxx>
  34 #include <TopoDS_Face.hxx>
  35 #include <TopoDS_Shape.hxx>
  36 #include <TopoDS_Vertex.hxx>
  37 #include <TopoDS_Iterator.hxx>
  38
  39 #include <Geom_CylindricalSurface.hxx>
  40 #include <Geom_Plane.hxx>
  41 #include <Geom_Surface.hxx>
  42
  43 #include <Precision.hxx>
  44 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
  45 #include <TColStd_MapOfInteger.hxx>
  46 #include <TColStd_SequenceOfAsciiString.hxx>
  47 #include <TColgp_Array1OfXYZ.hxx>
  48
  49 #include <gp_Ax3.hxx>
  50 #include <gp_Cylinder.hxx>
  51 #include <gp_Dir.hxx>
  52 #include <gp_Pln.hxx>
  53 #include <gp_Pnt.hxx>
  54 #include <gp_Vec.hxx>
  55 #include <gp_XYZ.hxx>
  56
  57 #include "SMDS_Mesh.hxx"
  58 #include "SMDS_Iterator.hxx"
  59 #include "SMDS_MeshElement.hxx"
  60 #include "SMDS_MeshNode.hxx"
  61 #include "SMDS_VolumeTool.hxx"
  62 #include "SMDS_QuadraticFaceOfNodes.hxx"
  63 #include "SMDS_QuadraticEdge.hxx"
  64
  65 #include "SMESHDS_Mesh.hxx"
  66 #include "SMESHDS_GroupBase.hxx"
  67
  68 /*
  69                             AUXILIARY METHODS
  70 */
  71
  72 namespace{
  73   inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
  74   {
  75     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
  76
  77     return v1.Magnitude() < gp::Resolution() ||
  78       v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
  79   }
  80
  81   inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
  82   {
  83     gp_Vec aVec1( P2 - P1 );
  84     gp_Vec aVec2( P3 - P1 );
  85     return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
  86   }
  87
  88   inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
  89   {
  90     return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
  91   }
  92
  93
  94
  95   inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
  96   {
  97     double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
  98     return aDist;
  99   }
 100
 101   int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId )
 102   {
 103     if ( theMesh == 0 )
 104       return 0;
 105
 106     const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
 107     if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */)
 108       return 0;
 109
 110     // for each pair of nodes in anEdge (there are 2 pairs in a quadratic edge)
 111     // count elements containing both nodes of the pair.
 112     // Note that there may be such cases for a quadratic edge (a horizontal line):
 113     //
 114     //  Case 1          Case 2
 115     //  |     |      |        |      |
 116     //  |     |      |        |      |
 117     //  +-----+------+  +-----+------+
 118     //  |            |  |            |
 119     //  |            |  |            |
 120     // result sould be 2 in both cases
 121     //
 122     int aResult0 = 0, aResult1 = 0;
 123      // last node, it is a medium one in a quadratic edge
 124     const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
 125     const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 );
 126     const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 );
 127     if ( aNode1 == aLastNode ) aNode1 = 0;
 128
 129     SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator();
 130     while( anElemIter->more() ) {
 131       const SMDS_MeshElement* anElem = anElemIter->next();
 132       if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
 133         SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
 134         while ( anIter->more() ) {
 135           if ( const SMDS_MeshElement* anElemNode = anIter->next() ) {
 136             if ( anElemNode == aNode0 ) {
 137               aResult0++;
 138               if ( !aNode1 ) break; // not a quadratic edge
 139             }
 140             else if ( anElemNode == aNode1 )
 141               aResult1++;
 142           }
 143         }
 144       }
 145     }
 146     int aResult = std::max ( aResult0, aResult1 );
 147
 148 //     TColStd_MapOfInteger aMap;
 149
 150 //     SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
 151 //     if ( anIter != 0 ) {
 152 //       while( anIter->more() ) {
 153 //      const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
 154 //      if ( aNode == 0 )
 155 //        return 0;
 156 //      SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
 157 //      while( anElemIter->more() ) {
 158 //        const SMDS_MeshElement* anElem = anElemIter->next();
 159 //        if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
 160 //          int anId = anElem->GetID();
 161
 162 //          if ( anIter->more() )              // i.e. first node
 163 //            aMap.Add( anId );
 164 //          else if ( aMap.Contains( anId ) )
 165 //            aResult++;
 166 //        }
 167 //      }
 168 //       }
 169 //     }
 170
 171     return aResult;
 172   }
 173
 174 }
 175
 176
 177
 178 using namespace SMESH::Controls;
 179
 180 /*
 181                                 FUNCTORS
 182 */
 183
 184 /*
 185   Class       : NumericalFunctor
 186   Description : Base class for numerical functors
 187 */
 188 NumericalFunctor::NumericalFunctor():
 189   myMesh(NULL)
 190 {
 191   myPrecision = -1;
 192 }
 193
 194 void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh )
 195 {
 196   myMesh = theMesh;
 197 }
 198
 199 bool NumericalFunctor::GetPoints(const int theId,
 200                                  TSequenceOfXYZ& theRes ) const
 201 {
 202   theRes.clear();
 203
 204   if ( myMesh == 0 )
 205     return false;
 206
 207   return GetPoints( myMesh->FindElement( theId ), theRes );
 208 }
 209
 210 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
 211                                  TSequenceOfXYZ&         theRes )
 212 {
 213   theRes.clear();
 214
 215   if ( anElem == 0)
 216     return false;
 217
 218   theRes.reserve( anElem->NbNodes() );
 219
 220   // Get nodes of the element
 221   SMDS_ElemIteratorPtr anIter;
 222
 223   if ( anElem->IsQuadratic() ) {
 224     switch ( anElem->GetType() ) {
 225     case SMDSAbs_Edge:
 226       anIter = static_cast<const SMDS_QuadraticEdge*>
 227         (anElem)->interlacedNodesElemIterator();
 228       break;
 229     case SMDSAbs_Face:
 230       anIter = static_cast<const SMDS_QuadraticFaceOfNodes*>
 231         (anElem)->interlacedNodesElemIterator();
 232       break;
 233     default:
 234       anIter = anElem->nodesIterator();
 235       //return false;
 236     }
 237   }
 238   else {
 239     anIter = anElem->nodesIterator();
 240   }
 241
 242   if ( anIter ) {
 243     while( anIter->more() ) {
 244       if ( const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( anIter->next() ))
 245         theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
 246     }
 247   }
 248
 249   return true;
 250 }
 251
 252 long  NumericalFunctor::GetPrecision() const
 253 {
 254   return myPrecision;
 255 }
 256
 257 void  NumericalFunctor::SetPrecision( const long thePrecision )
 258 {
 259   myPrecision = thePrecision;
 260 }
 261
 262 double NumericalFunctor::GetValue( long theId )
 263 {
 264   myCurrElement = myMesh->FindElement( theId );
 265   TSequenceOfXYZ P;
 266   if ( GetPoints( theId, P ))
 267   {
 268     double aVal = GetValue( P );
 269     if ( myPrecision >= 0 )
 270     {
 271       double prec = pow( 10., (double)( myPrecision ) );
 272       aVal = floor( aVal * prec + 0.5 ) / prec;
 273     }
 274     return aVal;
 275   }
 276
 277   return 0.;
 278 }
 279
 280 //================================================================================
 281 /*!
 282  * \brief Return histogram of functor values
 283  *  \param nbIntervals - number of intervals
 284  *  \param nbEvents - number of mesh elements having values within i-th interval
 285  *  \param funValues - boundaries of intervals
 286  *  \param elements - elements to check vulue of; empty list means "of all"
 287  */
 288 //================================================================================
 289
 290 void NumericalFunctor::GetHistogram(int                  nbIntervals,
 291                                     std::vector<int>&    nbEvents,
 292                                     std::vector<double>& funValues,
 293                                     const vector<int>&   elements)
 294 {
 295   if ( nbIntervals < 1 ||
 296        !myMesh ||
 297        !myMesh->GetMeshInfo().NbElements( GetType() ))
 298     return;
 299   nbEvents.resize( nbIntervals, 0 );
 300   funValues.resize( nbIntervals+1 );
 301
 302   // get all values sorted
 303   std::multiset< double > values;
 304   if ( elements.empty() )
 305   {
 306     SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType());
 307     while ( elemIt->more() )
 308       values.insert( GetValue( elemIt->next()->GetID() ));
 309   }
 310   else
 311   {
 312     vector<int>::const_iterator id = elements.begin();
 313     for ( ; id != elements.end(); ++id )
 314       values.insert( GetValue( *id ));
 315   }
 316
 317   // case nbIntervals == 1
 318   funValues[0] = *values.begin();
 319   funValues[nbIntervals] = *values.rbegin();
 320   if ( nbIntervals == 1 )
 321   {
 322     nbEvents[0] = values.size();
 323     return;
 324   }
 325   // case of 1 value
 326   if (funValues.front() == funValues.back())
 327   {
 328     nbEvents.resize( 1 );
 329     nbEvents[0] = values.size();
 330     funValues[1] = funValues.back();
 331     funValues.resize( 2 );
 332   }
 333   // generic case
 334   std::multiset< double >::iterator min = values.begin(), max;
 335   for ( int i = 0; i < nbIntervals; ++i )
 336   {
 337     double r = (i+1) / double( nbIntervals );
 338     funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
 339     if ( min != values.end() && *min <= funValues[i+1] )
 340     {
 341       max = values.upper_bound( funValues[i+1] ); // greater than funValues[i+1], or end()
 342       nbEvents[i] = std::distance( min, max );
 343       min = max;
 344     }
 345   }
 346 }
 347
 348 //=======================================================================
 349 //function : GetValue
 350 //purpose  :
 351 //=======================================================================
 352
 353 double Volume::GetValue( long theElementId )
 354 {
 355   if ( theElementId && myMesh ) {
 356     SMDS_VolumeTool aVolumeTool;
 357     if ( aVolumeTool.Set( myMesh->FindElement( theElementId )))
 358       return aVolumeTool.GetSize();
 359   }
 360   return 0;
 361 }
 362
 363 //=======================================================================
 364 //function : GetBadRate
 365 //purpose  : meaningless as it is not quality control functor
 366 //=======================================================================
 367
 368 double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
 369 {
 370   return Value;
 371 }
 372
 373 //=======================================================================
 374 //function : GetType
 375 //purpose  :
 376 //=======================================================================
 377
 378 SMDSAbs_ElementType Volume::GetType() const
 379 {
 380   return SMDSAbs_Volume;
 381 }
 382
 383
 384 /*
 385   Class       : MaxElementLength2D
 386   Description : Functor calculating maximum length of 2D element
 387 */
 388
 389 double MaxElementLength2D::GetValue( long theElementId )
 390 {
 391   TSequenceOfXYZ P;
 392   if( GetPoints( theElementId, P ) ) {
 393     double aVal = 0;
 394     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
 395     SMDSAbs_ElementType aType = aElem->GetType();
 396     int len = P.size();
 397     switch( aType ) {
 398     case SMDSAbs_Face:
 399       if( len == 3 ) { // triangles
 400         double L1 = getDistance(P( 1 ),P( 2 ));
 401         double L2 = getDistance(P( 2 ),P( 3 ));
 402         double L3 = getDistance(P( 3 ),P( 1 ));
 403         aVal = Max(L1,Max(L2,L3));
 404         break;
 405       }
 406       else if( len == 4 ) { // quadrangles
 407         double L1 = getDistance(P( 1 ),P( 2 ));
 408         double L2 = getDistance(P( 2 ),P( 3 ));
 409         double L3 = getDistance(P( 3 ),P( 4 ));
 410         double L4 = getDistance(P( 4 ),P( 1 ));
 411         double D1 = getDistance(P( 1 ),P( 3 ));
 412         double D2 = getDistance(P( 2 ),P( 4 ));
 413         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
 414         break;
 415       }
 416       else if( len == 6 ) { // quadratic triangles
 417         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
 418         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
 419         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
 420         aVal = Max(L1,Max(L2,L3));
 421         break;
 422       }
 423       else if( len == 8 ) { // quadratic quadrangles
 424         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
 425         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
 426         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
 427         double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
 428         double D1 = getDistance(P( 1 ),P( 5 ));
 429         double D2 = getDistance(P( 3 ),P( 7 ));
 430         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
 431         break;
 432       }
 433     }
 434
 435     if( myPrecision >= 0 )
 436     {
 437       double prec = pow( 10., (double)myPrecision );
 438       aVal = floor( aVal * prec + 0.5 ) / prec;
 439     }
 440     return aVal;
 441   }
 442   return 0.;
 443 }
 444
 445 double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const
 446 {
 447   return Value;
 448 }
 449
 450 SMDSAbs_ElementType MaxElementLength2D::GetType() const
 451 {
 452   return SMDSAbs_Face;
 453 }
 454
 455 /*
 456   Class       : MaxElementLength3D
 457   Description : Functor calculating maximum length of 3D element
 458 */
 459
 460 double MaxElementLength3D::GetValue( long theElementId )
 461 {
 462   TSequenceOfXYZ P;
 463   if( GetPoints( theElementId, P ) ) {
 464     double aVal = 0;
 465     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
 466     SMDSAbs_ElementType aType = aElem->GetType();
 467     int len = P.size();
 468     switch( aType ) {
 469     case SMDSAbs_Volume:
 470       if( len == 4 ) { // tetras
 471         double L1 = getDistance(P( 1 ),P( 2 ));
 472         double L2 = getDistance(P( 2 ),P( 3 ));
 473         double L3 = getDistance(P( 3 ),P( 1 ));
 474         double L4 = getDistance(P( 1 ),P( 4 ));
 475         double L5 = getDistance(P( 2 ),P( 4 ));
 476         double L6 = getDistance(P( 3 ),P( 4 ));
 477         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 478         break;
 479       }
 480       else if( len == 5 ) { // pyramids
 481         double L1 = getDistance(P( 1 ),P( 2 ));
 482         double L2 = getDistance(P( 2 ),P( 3 ));
 483         double L3 = getDistance(P( 3 ),P( 4 ));
 484         double L4 = getDistance(P( 4 ),P( 1 ));
 485         double L5 = getDistance(P( 1 ),P( 5 ));
 486         double L6 = getDistance(P( 2 ),P( 5 ));
 487         double L7 = getDistance(P( 3 ),P( 5 ));
 488         double L8 = getDistance(P( 4 ),P( 5 ));
 489         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 490         aVal = Max(aVal,Max(L7,L8));
 491         break;
 492       }
 493       else if( len == 6 ) { // pentas
 494         double L1 = getDistance(P( 1 ),P( 2 ));
 495         double L2 = getDistance(P( 2 ),P( 3 ));
 496         double L3 = getDistance(P( 3 ),P( 1 ));
 497         double L4 = getDistance(P( 4 ),P( 5 ));
 498         double L5 = getDistance(P( 5 ),P( 6 ));
 499         double L6 = getDistance(P( 6 ),P( 4 ));
 500         double L7 = getDistance(P( 1 ),P( 4 ));
 501         double L8 = getDistance(P( 2 ),P( 5 ));
 502         double L9 = getDistance(P( 3 ),P( 6 ));
 503         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 504         aVal = Max(aVal,Max(Max(L7,L8),L9));
 505         break;
 506       }
 507       else if( len == 8 ) { // hexas
 508         double L1 = getDistance(P( 1 ),P( 2 ));
 509         double L2 = getDistance(P( 2 ),P( 3 ));
 510         double L3 = getDistance(P( 3 ),P( 4 ));
 511         double L4 = getDistance(P( 4 ),P( 1 ));
 512         double L5 = getDistance(P( 5 ),P( 6 ));
 513         double L6 = getDistance(P( 6 ),P( 7 ));
 514         double L7 = getDistance(P( 7 ),P( 8 ));
 515         double L8 = getDistance(P( 8 ),P( 5 ));
 516         double L9 = getDistance(P( 1 ),P( 5 ));
 517         double L10= getDistance(P( 2 ),P( 6 ));
 518         double L11= getDistance(P( 3 ),P( 7 ));
 519         double L12= getDistance(P( 4 ),P( 8 ));
 520         double D1 = getDistance(P( 1 ),P( 7 ));
 521         double D2 = getDistance(P( 2 ),P( 8 ));
 522         double D3 = getDistance(P( 3 ),P( 5 ));
 523         double D4 = getDistance(P( 4 ),P( 6 ));
 524         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 525         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
 526         aVal = Max(aVal,Max(L11,L12));
 527         aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
 528         break;
 529       }
 530       else if( len == 10 ) { // quadratic tetras
 531         double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
 532         double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
 533         double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
 534         double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
 535         double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
 536         double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
 537         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 538         break;
 539       }
 540       else if( len == 13 ) { // quadratic pyramids
 541         double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
 542         double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
 543         double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
 544         double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
 545         double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
 546         double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
 547         double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
 548         double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
 549         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 550         aVal = Max(aVal,Max(L7,L8));
 551         break;
 552       }
 553       else if( len == 15 ) { // quadratic pentas
 554         double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
 555         double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
 556         double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
 557         double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
 558         double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
 559         double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
 560         double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
 561         double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
 562         double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
 563         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 564         aVal = Max(aVal,Max(Max(L7,L8),L9));
 565         break;
 566       }
 567       else if( len == 20 ) { // quadratic hexas
 568         double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
 569         double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
 570         double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
 571         double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
 572         double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
 573         double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
 574         double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
 575         double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
 576         double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
 577         double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
 578         double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
 579         double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
 580         double D1 = getDistance(P( 1 ),P( 7 ));
 581         double D2 = getDistance(P( 2 ),P( 8 ));
 582         double D3 = getDistance(P( 3 ),P( 5 ));
 583         double D4 = getDistance(P( 4 ),P( 6 ));
 584         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
 585         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
 586         aVal = Max(aVal,Max(L11,L12));
 587         aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
 588         break;
 589       }
 590       else if( len > 1 && aElem->IsPoly() ) { // polys
 591         // get the maximum distance between all pairs of nodes
 592         for( int i = 1; i <= len; i++ ) {
 593           for( int j = 1; j <= len; j++ ) {
 594             if( j > i ) { // optimization of the loop
 595               double D = getDistance( P(i), P(j) );
 596               aVal = Max( aVal, D );
 597             }
 598           }
 599         }
 600       }
 601     }
 602
 603     if( myPrecision >= 0 )
 604     {
 605       double prec = pow( 10., (double)myPrecision );
 606       aVal = floor( aVal * prec + 0.5 ) / prec;
 607     }
 608     return aVal;
 609   }
 610   return 0.;
 611 }
 612
 613 double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const
 614 {
 615   return Value;
 616 }
 617
 618 SMDSAbs_ElementType MaxElementLength3D::GetType() const
 619 {
 620   return SMDSAbs_Volume;
 621 }
 622
 623
 624 /*
 625   Class       : MinimumAngle
 626   Description : Functor for calculation of minimum angle
 627 */
 628
 629 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
 630 {
 631   double aMin;
 632
 633   if (P.size() <3)
 634     return 0.;
 635
 636   aMin = getAngle(P( P.size() ), P( 1 ), P( 2 ));
 637   aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 )));
 638
 639   for (int i=2; i<P.size();i++){
 640       double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) );
 641     aMin = Min(aMin,A0);
 642   }
 643
 644   return aMin * 180.0 / PI;
 645 }
 646
 647 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
 648 {
 649   //const double aBestAngle = PI / nbNodes;
 650   const double aBestAngle = 180.0 - ( 360.0 / double(nbNodes) );
 651   return ( fabs( aBestAngle - Value ));
 652 }
 653
 654 SMDSAbs_ElementType MinimumAngle::GetType() const
 655 {
 656   return SMDSAbs_Face;
 657 }
 658
 659
 660 /*
 661   Class       : AspectRatio
 662   Description : Functor for calculating aspect ratio
 663 */
 664 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
 665 {
 666   // According to "Mesh quality control" by Nadir Bouhamau referring to
 667   // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis.
 668   // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4
 669   // PAL10872
 670
 671   int nbNodes = P.size();
 672
 673   if ( nbNodes < 3 )
 674     return 0;
 675
 676   // Compute aspect ratio
 677
 678   if ( nbNodes == 3 ) {
 679     // Compute lengths of the sides
 680     std::vector< double > aLen (nbNodes);
 681     for ( int i = 0; i < nbNodes - 1; i++ )
 682       aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
 683     aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
 684     // Q = alfa * h * p / S, where
 685     //
 686     // alfa = sqrt( 3 ) / 6
 687     // h - length of the longest edge
 688     // p - half perimeter
 689     // S - triangle surface
 690     const double alfa = sqrt( 3. ) / 6.;
 691     double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
 692     double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
 693     double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
 694     if ( anArea <= Precision::Confusion() )
 695       return 0.;
 696     return alfa * maxLen * half_perimeter / anArea;
 697   }
 698   else if ( nbNodes == 6 ) { // quadratic triangles
 699     // Compute lengths of the sides
 700     std::vector< double > aLen (3);
 701     aLen[0] = getDistance( P(1), P(3) );
 702     aLen[1] = getDistance( P(3), P(5) );
 703     aLen[2] = getDistance( P(5), P(1) );
 704     // Q = alfa * h * p / S, where
 705     //
 706     // alfa = sqrt( 3 ) / 6
 707     // h - length of the longest edge
 708     // p - half perimeter
 709     // S - triangle surface
 710     const double alfa = sqrt( 3. ) / 6.;
 711     double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
 712     double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
 713     double anArea = getArea( P(1), P(3), P(5) );
 714     if ( anArea <= Precision::Confusion() )
 715       return 0.;
 716     return alfa * maxLen * half_perimeter / anArea;
 717   }
 718   else if( nbNodes == 4 ) { // quadrangle
 719     // Compute lengths of the sides
 720     std::vector< double > aLen (4);
 721     aLen[0] = getDistance( P(1), P(2) );
 722     aLen[1] = getDistance( P(2), P(3) );
 723     aLen[2] = getDistance( P(3), P(4) );
 724     aLen[3] = getDistance( P(4), P(1) );
 725     // Compute lengths of the diagonals
 726     std::vector< double > aDia (2);
 727     aDia[0] = getDistance( P(1), P(3) );
 728     aDia[1] = getDistance( P(2), P(4) );
 729     // Compute areas of all triangles which can be built
 730     // taking three nodes of the quadrangle
 731     std::vector< double > anArea (4);
 732     anArea[0] = getArea( P(1), P(2), P(3) );
 733     anArea[1] = getArea( P(1), P(2), P(4) );
 734     anArea[2] = getArea( P(1), P(3), P(4) );
 735     anArea[3] = getArea( P(2), P(3), P(4) );
 736     // Q = alpha * L * C1 / C2, where
 737     //
 738     // alpha = sqrt( 1/32 )
 739     // L = max( L1, L2, L3, L4, D1, D2 )
 740     // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
 741     // C2 = min( S1, S2, S3, S4 )
 742     // Li - lengths of the edges
 743     // Di - lengths of the diagonals
 744     // Si - areas of the triangles
 745     const double alpha = sqrt( 1 / 32. );
 746     double L = Max( aLen[ 0 ],
 747                  Max( aLen[ 1 ],
 748                    Max( aLen[ 2 ],
 749                      Max( aLen[ 3 ],
 750                        Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
 751     double C1 = sqrt( ( aLen[0] * aLen[0] +
 752                         aLen[1] * aLen[1] +
 753                         aLen[2] * aLen[2] +
 754                         aLen[3] * aLen[3] ) / 4. );
 755     double C2 = Min( anArea[ 0 ],
 756                   Min( anArea[ 1 ],
 757                     Min( anArea[ 2 ], anArea[ 3 ] ) ) );
 758     if ( C2 <= Precision::Confusion() )
 759       return 0.;
 760     return alpha * L * C1 / C2;
 761   }
 762   else if( nbNodes == 8 ){ // nbNodes==8 - quadratic quadrangle
 763     // Compute lengths of the sides
 764     std::vector< double > aLen (4);
 765     aLen[0] = getDistance( P(1), P(3) );
 766     aLen[1] = getDistance( P(3), P(5) );
 767     aLen[2] = getDistance( P(5), P(7) );
 768     aLen[3] = getDistance( P(7), P(1) );
 769     // Compute lengths of the diagonals
 770     std::vector< double > aDia (2);
 771     aDia[0] = getDistance( P(1), P(5) );
 772     aDia[1] = getDistance( P(3), P(7) );
 773     // Compute areas of all triangles which can be built
 774     // taking three nodes of the quadrangle
 775     std::vector< double > anArea (4);
 776     anArea[0] = getArea( P(1), P(3), P(5) );
 777     anArea[1] = getArea( P(1), P(3), P(7) );
 778     anArea[2] = getArea( P(1), P(5), P(7) );
 779     anArea[3] = getArea( P(3), P(5), P(7) );
 780     // Q = alpha * L * C1 / C2, where
 781     //
 782     // alpha = sqrt( 1/32 )
 783     // L = max( L1, L2, L3, L4, D1, D2 )
 784     // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
 785     // C2 = min( S1, S2, S3, S4 )
 786     // Li - lengths of the edges
 787     // Di - lengths of the diagonals
 788     // Si - areas of the triangles
 789     const double alpha = sqrt( 1 / 32. );
 790     double L = Max( aLen[ 0 ],
 791                  Max( aLen[ 1 ],
 792                    Max( aLen[ 2 ],
 793                      Max( aLen[ 3 ],
 794                        Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
 795     double C1 = sqrt( ( aLen[0] * aLen[0] +
 796                         aLen[1] * aLen[1] +
 797                         aLen[2] * aLen[2] +
 798                         aLen[3] * aLen[3] ) / 4. );
 799     double C2 = Min( anArea[ 0 ],
 800                   Min( anArea[ 1 ],
 801                     Min( anArea[ 2 ], anArea[ 3 ] ) ) );
 802     if ( C2 <= Precision::Confusion() )
 803       return 0.;
 804     return alpha * L * C1 / C2;
 805   }
 806   return 0;
 807 }
 808
 809 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
 810 {
 811   // the aspect ratio is in the range [1.0,infinity]
 812   // 1.0 = good
 813   // infinity = bad
 814   return Value / 1000.;
 815 }
 816
 817 SMDSAbs_ElementType AspectRatio::GetType() const
 818 {
 819   return SMDSAbs_Face;
 820 }
 821
 822
 823 /*
 824   Class       : AspectRatio3D
 825   Description : Functor for calculating aspect ratio
 826 */
 827 namespace{
 828
 829   inline double getHalfPerimeter(double theTria[3]){
 830     return (theTria[0] + theTria[1] + theTria[2])/2.0;
 831   }
 832
 833   inline double getArea(double theHalfPerim, double theTria[3]){
 834     return sqrt(theHalfPerim*
 835                 (theHalfPerim-theTria[0])*
 836                 (theHalfPerim-theTria[1])*
 837                 (theHalfPerim-theTria[2]));
 838   }
 839
 840   inline double getVolume(double theLen[6]){
 841     double a2 = theLen[0]*theLen[0];
 842     double b2 = theLen[1]*theLen[1];
 843     double c2 = theLen[2]*theLen[2];
 844     double d2 = theLen[3]*theLen[3];
 845     double e2 = theLen[4]*theLen[4];
 846     double f2 = theLen[5]*theLen[5];
 847     double P = 4.0*a2*b2*d2;
 848     double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
 849     double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
 850     return sqrt(P-Q+R)/12.0;
 851   }
 852
 853   inline double getVolume2(double theLen[6]){
 854     double a2 = theLen[0]*theLen[0];
 855     double b2 = theLen[1]*theLen[1];
 856     double c2 = theLen[2]*theLen[2];
 857     double d2 = theLen[3]*theLen[3];
 858     double e2 = theLen[4]*theLen[4];
 859     double f2 = theLen[5]*theLen[5];
 860
 861     double P = a2*e2*(b2+c2+d2+f2-a2-e2);
 862     double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
 863     double R = c2*d2*(a2+b2+e2+f2-c2-d2);
 864     double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
 865
 866     return sqrt(P+Q+R-S)/12.0;
 867   }
 868
 869   inline double getVolume(const TSequenceOfXYZ& P){
 870     gp_Vec aVec1( P( 2 ) - P( 1 ) );
 871     gp_Vec aVec2( P( 3 ) - P( 1 ) );
 872     gp_Vec aVec3( P( 4 ) - P( 1 ) );
 873     gp_Vec anAreaVec( aVec1 ^ aVec2 );
 874     return fabs(aVec3 * anAreaVec) / 6.0;
 875   }
 876
 877   inline double getMaxHeight(double theLen[6])
 878   {
 879     double aHeight = std::max(theLen[0],theLen[1]);
 880     aHeight = std::max(aHeight,theLen[2]);
 881     aHeight = std::max(aHeight,theLen[3]);
 882     aHeight = std::max(aHeight,theLen[4]);
 883     aHeight = std::max(aHeight,theLen[5]);
 884     return aHeight;
 885   }
 886
 887 }
 888
 889 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
 890 {
 891   double aQuality = 0.0;
 892   if(myCurrElement->IsPoly()) return aQuality;
 893
 894   int nbNodes = P.size();
 895
 896   if(myCurrElement->IsQuadratic()) {
 897     if(nbNodes==10) nbNodes=4; // quadratic tetrahedron
 898     else if(nbNodes==13) nbNodes=5; // quadratic pyramid
 899     else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
 900     else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
 901     else return aQuality;
 902   }
 903
 904   switch(nbNodes){
 905   case 4:{
 906     double aLen[6] = {
 907       getDistance(P( 1 ),P( 2 )), // a
 908       getDistance(P( 2 ),P( 3 )), // b
 909       getDistance(P( 3 ),P( 1 )), // c
 910       getDistance(P( 2 ),P( 4 )), // d
 911       getDistance(P( 3 ),P( 4 )), // e
 912       getDistance(P( 1 ),P( 4 ))  // f
 913     };
 914     double aTria[4][3] = {
 915       {aLen[0],aLen[1],aLen[2]}, // abc
 916       {aLen[0],aLen[3],aLen[5]}, // adf
 917       {aLen[1],aLen[3],aLen[4]}, // bde
 918       {aLen[2],aLen[4],aLen[5]}  // cef
 919     };
 920     double aSumArea = 0.0;
 921     double aHalfPerimeter = getHalfPerimeter(aTria[0]);
 922     double anArea = getArea(aHalfPerimeter,aTria[0]);
 923     aSumArea += anArea;
 924     aHalfPerimeter = getHalfPerimeter(aTria[1]);
 925     anArea = getArea(aHalfPerimeter,aTria[1]);
 926     aSumArea += anArea;
 927     aHalfPerimeter = getHalfPerimeter(aTria[2]);
 928     anArea = getArea(aHalfPerimeter,aTria[2]);
 929     aSumArea += anArea;
 930     aHalfPerimeter = getHalfPerimeter(aTria[3]);
 931     anArea = getArea(aHalfPerimeter,aTria[3]);
 932     aSumArea += anArea;
 933     double aVolume = getVolume(P);
 934     //double aVolume = getVolume(aLen);
 935     double aHeight = getMaxHeight(aLen);
 936     static double aCoeff = sqrt(2.0)/12.0;
 937     if ( aVolume > DBL_MIN )
 938       aQuality = aCoeff*aHeight*aSumArea/aVolume;
 939     break;
 940   }
 941   case 5:{
 942     {
 943       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
 944       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 945     }
 946     {
 947       gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
 948       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 949     }
 950     {
 951       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
 952       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 953     }
 954     {
 955       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
 956       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 957     }
 958     break;
 959   }
 960   case 6:{
 961     {
 962       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
 963       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 964     }
 965     {
 966       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
 967       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 968     }
 969     {
 970       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
 971       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 972     }
 973     {
 974       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
 975       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 976     }
 977     {
 978       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
 979       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 980     }
 981     {
 982       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
 983       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 984     }
 985     break;
 986   }
 987   case 8:{
 988     {
 989       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
 990       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 991     }
 992     {
 993       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
 994       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 995     }
 996     {
 997       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
 998       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
 999     }
1000     {
1001       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
1002       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1003     }
1004     {
1005       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
1006       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1007     }
1008     {
1009       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
1010       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1011     }
1012     {
1013       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
1014       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1015     }
1016     {
1017       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
1018       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1019     }
1020     {
1021       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
1022       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1023     }
1024     {
1025       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
1026       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1027     }
1028     {
1029       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
1030       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1031     }
1032     {
1033       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
1034       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1035     }
1036     {
1037       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
1038       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1039     }
1040     {
1041       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
1042       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1043     }
1044     {
1045       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
1046       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1047     }
1048     {
1049       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
1050       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1051     }
1052     {
1053       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
1054       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1055     }
1056     {
1057       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
1058       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1059     }
1060     {
1061       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
1062       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1063     }
1064     {
1065       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
1066       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1067     }
1068     {
1069       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
1070       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1071     }
1072     {
1073       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1074       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1075     }
1076     {
1077       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
1078       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1079     }
1080     {
1081       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
1082       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1083     }
1084     {
1085       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1086       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1087     }
1088     {
1089       gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
1090       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1091     }
1092     {
1093       gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
1094       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1095     }
1096     {
1097       gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
1098       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1099     }
1100     {
1101       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
1102       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1103     }
1104     {
1105       gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
1106       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1107     }
1108     {
1109       gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
1110       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1111     }
1112     {
1113       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
1114       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1115     }
1116     {
1117       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
1118       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1119     }
1120     break;
1121   }
1122   }
1123   if ( nbNodes > 4 ) {
1124     // avaluate aspect ratio of quadranle faces
1125     AspectRatio aspect2D;
1126     SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
1127     int nbFaces = SMDS_VolumeTool::NbFaces( type );
1128     TSequenceOfXYZ points(4);
1129     for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume
1130       if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
1131         continue;
1132       const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
1133       for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face
1134         points( p + 1 ) = P( pInd[ p ] + 1 );
1135       aQuality = std::max( aQuality, aspect2D.GetValue( points ));
1136     }
1137   }
1138   return aQuality;
1139 }
1140
1141 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
1142 {
1143   // the aspect ratio is in the range [1.0,infinity]
1144   // 1.0 = good
1145   // infinity = bad
1146   return Value / 1000.;
1147 }
1148
1149 SMDSAbs_ElementType AspectRatio3D::GetType() const
1150 {
1151   return SMDSAbs_Volume;
1152 }
1153
1154
1155 /*
1156   Class       : Warping
1157   Description : Functor for calculating warping
1158 */
1159 double Warping::GetValue( const TSequenceOfXYZ& P )
1160 {
1161   if ( P.size() != 4 )
1162     return 0;
1163
1164   gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.;
1165
1166   double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
1167   double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
1168   double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
1169   double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
1170
1171   return Max( Max( A1, A2 ), Max( A3, A4 ) );
1172 }
1173
1174 double Warping::ComputeA( const gp_XYZ& thePnt1,
1175                           const gp_XYZ& thePnt2,
1176                           const gp_XYZ& thePnt3,
1177                           const gp_XYZ& theG ) const
1178 {
1179   double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
1180   double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
1181   double L = Min( aLen1, aLen2 ) * 0.5;
1182   if ( L < Precision::Confusion())
1183     return 0.;
1184
1185   gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
1186   gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
1187   gp_XYZ N  = GI.Crossed( GJ );
1188
1189   if ( N.Modulus() < gp::Resolution() )
1190     return PI / 2;
1191
1192   N.Normalize();
1193
1194   double H = ( thePnt2 - theG ).Dot( N );
1195   return asin( fabs( H / L ) ) * 180. / PI;
1196 }
1197
1198 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
1199 {
1200   // the warp is in the range [0.0,PI/2]
1201   // 0.0 = good (no warp)
1202   // PI/2 = bad  (face pliee)
1203   return Value;
1204 }
1205
1206 SMDSAbs_ElementType Warping::GetType() const
1207 {
1208   return SMDSAbs_Face;
1209 }
1210
1211
1212 /*
1213   Class       : Taper
1214   Description : Functor for calculating taper
1215 */
1216 double Taper::GetValue( const TSequenceOfXYZ& P )
1217 {
1218   if ( P.size() != 4 )
1219     return 0.;
1220
1221   // Compute taper
1222   double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2.;
1223   double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2.;
1224   double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2.;
1225   double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2.;
1226
1227   double JA = 0.25 * ( J1 + J2 + J3 + J4 );
1228   if ( JA <= Precision::Confusion() )
1229     return 0.;
1230
1231   double T1 = fabs( ( J1 - JA ) / JA );
1232   double T2 = fabs( ( J2 - JA ) / JA );
1233   double T3 = fabs( ( J3 - JA ) / JA );
1234   double T4 = fabs( ( J4 - JA ) / JA );
1235
1236   return Max( Max( T1, T2 ), Max( T3, T4 ) );
1237 }
1238
1239 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
1240 {
1241   // the taper is in the range [0.0,1.0]
1242   // 0.0  = good (no taper)
1243   // 1.0 = bad  (les cotes opposes sont allignes)
1244   return Value;
1245 }
1246
1247 SMDSAbs_ElementType Taper::GetType() const
1248 {
1249   return SMDSAbs_Face;
1250 }
1251
1252
1253 /*
1254   Class       : Skew
1255   Description : Functor for calculating skew in degrees
1256 */
1257 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
1258 {
1259   gp_XYZ p12 = ( p2 + p1 ) / 2.;
1260   gp_XYZ p23 = ( p3 + p2 ) / 2.;
1261   gp_XYZ p31 = ( p3 + p1 ) / 2.;
1262
1263   gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
1264
1265   return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
1266 }
1267
1268 double Skew::GetValue( const TSequenceOfXYZ& P )
1269 {
1270   if ( P.size() != 3 && P.size() != 4 )
1271     return 0.;
1272
1273   // Compute skew
1274   static double PI2 = PI / 2.;
1275   if ( P.size() == 3 )
1276   {
1277     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
1278     double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
1279     double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
1280
1281     return Max( A0, Max( A1, A2 ) ) * 180. / PI;
1282   }
1283   else
1284   {
1285     gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.;
1286     gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.;
1287     gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.;
1288     gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.;
1289
1290     gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
1291     double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
1292       ? 0. : fabs( PI2 - v1.Angle( v2 ) );
1293
1294     //BUG SWP12743
1295     if ( A < Precision::Angular() )
1296       return 0.;
1297
1298     return A * 180. / PI;
1299   }
1300 }
1301
1302 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
1303 {
1304   // the skew is in the range [0.0,PI/2].
1305   // 0.0 = good
1306   // PI/2 = bad
1307   return Value;
1308 }
1309
1310 SMDSAbs_ElementType Skew::GetType() const
1311 {
1312   return SMDSAbs_Face;
1313 }
1314
1315
1316 /*
1317   Class       : Area
1318   Description : Functor for calculating area
1319 */
1320 double Area::GetValue( const TSequenceOfXYZ& P )
1321 {
1322   double val = 0.0;
1323   if ( P.size() > 2 ) {
1324     gp_Vec aVec1( P(2) - P(1) );
1325     gp_Vec aVec2( P(3) - P(1) );
1326     gp_Vec SumVec = aVec1 ^ aVec2;
1327     for (int i=4; i<=P.size(); i++) {
1328       gp_Vec aVec1( P(i-1) - P(1) );
1329       gp_Vec aVec2( P(i) - P(1) );
1330       gp_Vec tmp = aVec1 ^ aVec2;
1331       SumVec.Add(tmp);
1332     }
1333     val = SumVec.Magnitude() * 0.5;
1334   }
1335   return val;
1336 }
1337
1338 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
1339 {
1340   // meaningless as it is not a quality control functor
1341   return Value;
1342 }
1343
1344 SMDSAbs_ElementType Area::GetType() const
1345 {
1346   return SMDSAbs_Face;
1347 }
1348
1349
1350 /*
1351   Class       : Length
1352   Description : Functor for calculating length off edge
1353 */
1354 double Length::GetValue( const TSequenceOfXYZ& P )
1355 {
1356   switch ( P.size() ) {
1357   case 2:  return getDistance( P( 1 ), P( 2 ) );
1358   case 3:  return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) );
1359   default: return 0.;
1360   }
1361 }
1362
1363 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
1364 {
1365   // meaningless as it is not quality control functor
1366   return Value;
1367 }
1368
1369 SMDSAbs_ElementType Length::GetType() const
1370 {
1371   return SMDSAbs_Edge;
1372 }
1373
1374 /*
1375   Class       : Length2D
1376   Description : Functor for calculating length of edge
1377 */
1378
1379 double Length2D::GetValue( long theElementId)
1380 {
1381   TSequenceOfXYZ P;
1382
1383   //cout<<"Length2D::GetValue"<<endl;
1384   if (GetPoints(theElementId,P)){
1385     //for(int jj=1; jj<=P.size(); jj++)
1386     //  cout<<"jj="<<jj<<" P("<<P(jj).X()<<","<<P(jj).Y()<<","<<P(jj).Z()<<")"<<endl;
1387
1388     double  aVal;// = GetValue( P );
1389     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
1390     SMDSAbs_ElementType aType = aElem->GetType();
1391
1392     int len = P.size();
1393
1394     switch (aType){
1395     case SMDSAbs_All:
1396     case SMDSAbs_Node:
1397     case SMDSAbs_Edge:
1398       if (len == 2){
1399         aVal = getDistance( P( 1 ), P( 2 ) );
1400         break;
1401       }
1402       else if (len == 3){ // quadratic edge
1403         aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
1404         break;
1405       }
1406     case SMDSAbs_Face:
1407       if (len == 3){ // triangles
1408         double L1 = getDistance(P( 1 ),P( 2 ));
1409         double L2 = getDistance(P( 2 ),P( 3 ));
1410         double L3 = getDistance(P( 3 ),P( 1 ));
1411         aVal = Max(L1,Max(L2,L3));
1412         break;
1413       }
1414       else if (len == 4){ // quadrangles
1415         double L1 = getDistance(P( 1 ),P( 2 ));
1416         double L2 = getDistance(P( 2 ),P( 3 ));
1417         double L3 = getDistance(P( 3 ),P( 4 ));
1418         double L4 = getDistance(P( 4 ),P( 1 ));
1419         aVal = Max(Max(L1,L2),Max(L3,L4));
1420         break;
1421       }
1422       if (len == 6){ // quadratic triangles
1423         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1424         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1425         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
1426         aVal = Max(L1,Max(L2,L3));
1427         //cout<<"L1="<<L1<<" L2="<<L2<<"L3="<<L3<<" aVal="<<aVal<<endl;
1428         break;
1429       }
1430       else if (len == 8){ // quadratic quadrangles
1431         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1432         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1433         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
1434         double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
1435         aVal = Max(Max(L1,L2),Max(L3,L4));
1436         break;
1437       }
1438     case SMDSAbs_Volume:
1439       if (len == 4){ // tetraidrs
1440         double L1 = getDistance(P( 1 ),P( 2 ));
1441         double L2 = getDistance(P( 2 ),P( 3 ));
1442         double L3 = getDistance(P( 3 ),P( 1 ));
1443         double L4 = getDistance(P( 1 ),P( 4 ));
1444         double L5 = getDistance(P( 2 ),P( 4 ));
1445         double L6 = getDistance(P( 3 ),P( 4 ));
1446         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1447         break;
1448       }
1449       else if (len == 5){ // piramids
1450         double L1 = getDistance(P( 1 ),P( 2 ));
1451         double L2 = getDistance(P( 2 ),P( 3 ));
1452         double L3 = getDistance(P( 3 ),P( 4 ));
1453         double L4 = getDistance(P( 4 ),P( 1 ));
1454         double L5 = getDistance(P( 1 ),P( 5 ));
1455         double L6 = getDistance(P( 2 ),P( 5 ));
1456         double L7 = getDistance(P( 3 ),P( 5 ));
1457         double L8 = getDistance(P( 4 ),P( 5 ));
1458
1459         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1460         aVal = Max(aVal,Max(L7,L8));
1461         break;
1462       }
1463       else if (len == 6){ // pentaidres
1464         double L1 = getDistance(P( 1 ),P( 2 ));
1465         double L2 = getDistance(P( 2 ),P( 3 ));
1466         double L3 = getDistance(P( 3 ),P( 1 ));
1467         double L4 = getDistance(P( 4 ),P( 5 ));
1468         double L5 = getDistance(P( 5 ),P( 6 ));
1469         double L6 = getDistance(P( 6 ),P( 4 ));
1470         double L7 = getDistance(P( 1 ),P( 4 ));
1471         double L8 = getDistance(P( 2 ),P( 5 ));
1472         double L9 = getDistance(P( 3 ),P( 6 ));
1473
1474         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1475         aVal = Max(aVal,Max(Max(L7,L8),L9));
1476         break;
1477       }
1478       else if (len == 8){ // hexaider
1479         double L1 = getDistance(P( 1 ),P( 2 ));
1480         double L2 = getDistance(P( 2 ),P( 3 ));
1481         double L3 = getDistance(P( 3 ),P( 4 ));
1482         double L4 = getDistance(P( 4 ),P( 1 ));
1483         double L5 = getDistance(P( 5 ),P( 6 ));
1484         double L6 = getDistance(P( 6 ),P( 7 ));
1485         double L7 = getDistance(P( 7 ),P( 8 ));
1486         double L8 = getDistance(P( 8 ),P( 5 ));
1487         double L9 = getDistance(P( 1 ),P( 5 ));
1488         double L10= getDistance(P( 2 ),P( 6 ));
1489         double L11= getDistance(P( 3 ),P( 7 ));
1490         double L12= getDistance(P( 4 ),P( 8 ));
1491
1492         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1493         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
1494         aVal = Max(aVal,Max(L11,L12));
1495         break;
1496
1497       }
1498
1499       if (len == 10){ // quadratic tetraidrs
1500         double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
1501         double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
1502         double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
1503         double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1504         double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
1505         double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
1506         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1507         break;
1508       }
1509       else if (len == 13){ // quadratic piramids
1510         double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
1511         double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
1512         double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1513         double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1514         double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1515         double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
1516         double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
1517         double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
1518         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1519         aVal = Max(aVal,Max(L7,L8));
1520         break;
1521       }
1522       else if (len == 15){ // quadratic pentaidres
1523         double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
1524         double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
1525         double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1526         double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1527         double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
1528         double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
1529         double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
1530         double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
1531         double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
1532         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1533         aVal = Max(aVal,Max(Max(L7,L8),L9));
1534         break;
1535       }
1536       else if (len == 20){ // quadratic hexaider
1537         double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
1538         double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
1539         double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
1540         double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
1541         double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
1542         double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
1543         double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
1544         double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
1545         double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
1546         double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
1547         double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
1548         double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
1549         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
1550         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
1551         aVal = Max(aVal,Max(L11,L12));
1552         break;
1553
1554       }
1555
1556     default: aVal=-1;
1557     }
1558
1559     if (aVal <0){
1560       return 0.;
1561     }
1562
1563     if ( myPrecision >= 0 )
1564     {
1565       double prec = pow( 10., (double)( myPrecision ) );
1566       aVal = floor( aVal * prec + 0.5 ) / prec;
1567     }
1568
1569     return aVal;
1570
1571   }
1572   return 0.;
1573 }
1574
1575 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1576 {
1577   // meaningless as it is not quality control functor
1578   return Value;
1579 }
1580
1581 SMDSAbs_ElementType Length2D::GetType() const
1582 {
1583   return SMDSAbs_Face;
1584 }
1585
1586 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
1587   myLength(theLength)
1588 {
1589   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
1590   if(thePntId1 > thePntId2){
1591     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
1592   }
1593 }
1594
1595 bool Length2D::Value::operator<(const Length2D::Value& x) const{
1596   if(myPntId[0] < x.myPntId[0]) return true;
1597   if(myPntId[0] == x.myPntId[0])
1598     if(myPntId[1] < x.myPntId[1]) return true;
1599   return false;
1600 }
1601
1602 void Length2D::GetValues(TValues& theValues){
1603   TValues aValues;
1604   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1605   for(; anIter->more(); ){
1606     const SMDS_MeshFace* anElem = anIter->next();
1607
1608     if(anElem->IsQuadratic()) {
1609       const SMDS_QuadraticFaceOfNodes* F =
1610         static_cast<const SMDS_QuadraticFaceOfNodes*>(anElem);
1611       // use special nodes iterator
1612       SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
1613       long aNodeId[4];
1614       gp_Pnt P[4];
1615
1616       double aLength;
1617       const SMDS_MeshElement* aNode;
1618       if(anIter->more()){
1619         aNode = anIter->next();
1620         const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1621         P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1622         aNodeId[0] = aNodeId[1] = aNode->GetID();
1623         aLength = 0;
1624       }
1625       for(; anIter->more(); ){
1626         const SMDS_MeshNode* N1 = static_cast<const SMDS_MeshNode*> (anIter->next());
1627         P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z());
1628         aNodeId[2] = N1->GetID();
1629         aLength = P[1].Distance(P[2]);
1630         if(!anIter->more()) break;
1631         const SMDS_MeshNode* N2 = static_cast<const SMDS_MeshNode*> (anIter->next());
1632         P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z());
1633         aNodeId[3] = N2->GetID();
1634         aLength += P[2].Distance(P[3]);
1635         Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1636         Value aValue2(aLength,aNodeId[2],aNodeId[3]);
1637         P[1] = P[3];
1638         aNodeId[1] = aNodeId[3];
1639         theValues.insert(aValue1);
1640         theValues.insert(aValue2);
1641       }
1642       aLength += P[2].Distance(P[0]);
1643       Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1644       Value aValue2(aLength,aNodeId[2],aNodeId[0]);
1645       theValues.insert(aValue1);
1646       theValues.insert(aValue2);
1647     }
1648     else {
1649       SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1650       long aNodeId[2];
1651       gp_Pnt P[3];
1652
1653       double aLength;
1654       const SMDS_MeshElement* aNode;
1655       if(aNodesIter->more()){
1656         aNode = aNodesIter->next();
1657         const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1658         P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1659         aNodeId[0] = aNodeId[1] = aNode->GetID();
1660         aLength = 0;
1661       }
1662       for(; aNodesIter->more(); ){
1663         aNode = aNodesIter->next();
1664         const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1665         long anId = aNode->GetID();
1666
1667         P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1668
1669         aLength = P[1].Distance(P[2]);
1670
1671         Value aValue(aLength,aNodeId[1],anId);
1672         aNodeId[1] = anId;
1673         P[1] = P[2];
1674         theValues.insert(aValue);
1675       }
1676
1677       aLength = P[0].Distance(P[1]);
1678
1679       Value aValue(aLength,aNodeId[0],aNodeId[1]);
1680       theValues.insert(aValue);
1681     }
1682   }
1683 }
1684
1685 /*
1686   Class       : MultiConnection
1687   Description : Functor for calculating number of faces conneted to the edge
1688 */
1689 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
1690 {
1691   return 0;
1692 }
1693 double MultiConnection::GetValue( long theId )
1694 {
1695   return getNbMultiConnection( myMesh, theId );
1696 }
1697
1698 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
1699 {
1700   // meaningless as it is not quality control functor
1701   return Value;
1702 }
1703
1704 SMDSAbs_ElementType MultiConnection::GetType() const
1705 {
1706   return SMDSAbs_Edge;
1707 }
1708
1709 /*
1710   Class       : MultiConnection2D
1711   Description : Functor for calculating number of faces conneted to the edge
1712 */
1713 double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
1714 {
1715   return 0;
1716 }
1717
1718 double MultiConnection2D::GetValue( long theElementId )
1719 {
1720   int aResult = 0;
1721
1722   const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
1723   SMDSAbs_ElementType aType = aFaceElem->GetType();
1724
1725   switch (aType) {
1726   case SMDSAbs_Face:
1727     {
1728       int i = 0, len = aFaceElem->NbNodes();
1729       SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
1730       if (!anIter) break;
1731
1732       const SMDS_MeshNode *aNode, *aNode0;
1733       TColStd_MapOfInteger aMap, aMapPrev;
1734
1735       for (i = 0; i <= len; i++) {
1736         aMapPrev = aMap;
1737         aMap.Clear();
1738
1739         int aNb = 0;
1740         if (anIter->more()) {
1741           aNode = (SMDS_MeshNode*)anIter->next();
1742         } else {
1743           if (i == len)
1744             aNode = aNode0;
1745           else
1746             break;
1747         }
1748         if (!aNode) break;
1749         if (i == 0) aNode0 = aNode;
1750
1751         SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
1752         while (anElemIter->more()) {
1753           const SMDS_MeshElement* anElem = anElemIter->next();
1754           if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
1755             int anId = anElem->GetID();
1756
1757             aMap.Add(anId);
1758             if (aMapPrev.Contains(anId)) {
1759               aNb++;
1760             }
1761           }
1762         }
1763         aResult = Max(aResult, aNb);
1764       }
1765     }
1766     break;
1767   default:
1768     aResult = 0;
1769   }
1770
1771   return aResult;
1772 }
1773
1774 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1775 {
1776   // meaningless as it is not quality control functor
1777   return Value;
1778 }
1779
1780 SMDSAbs_ElementType MultiConnection2D::GetType() const
1781 {
1782   return SMDSAbs_Face;
1783 }
1784
1785 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
1786 {
1787   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
1788   if(thePntId1 > thePntId2){
1789     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
1790   }
1791 }
1792
1793 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{
1794   if(myPntId[0] < x.myPntId[0]) return true;
1795   if(myPntId[0] == x.myPntId[0])
1796     if(myPntId[1] < x.myPntId[1]) return true;
1797   return false;
1798 }
1799
1800 void MultiConnection2D::GetValues(MValues& theValues){
1801   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1802   for(; anIter->more(); ){
1803     const SMDS_MeshFace* anElem = anIter->next();
1804     SMDS_ElemIteratorPtr aNodesIter;
1805     if ( anElem->IsQuadratic() )
1806       aNodesIter = static_cast<const SMDS_QuadraticFaceOfNodes*>
1807         (anElem)->interlacedNodesElemIterator();
1808     else
1809       aNodesIter = anElem->nodesIterator();
1810     long aNodeId[3];
1811
1812     //int aNbConnects=0;
1813     const SMDS_MeshNode* aNode0;
1814     const SMDS_MeshNode* aNode1;
1815     const SMDS_MeshNode* aNode2;
1816     if(aNodesIter->more()){
1817       aNode0 = (SMDS_MeshNode*) aNodesIter->next();
1818       aNode1 = aNode0;
1819       const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
1820       aNodeId[0] = aNodeId[1] = aNodes->GetID();
1821     }
1822     for(; aNodesIter->more(); ) {
1823       aNode2 = (SMDS_MeshNode*) aNodesIter->next();
1824       long anId = aNode2->GetID();
1825       aNodeId[2] = anId;
1826
1827       Value aValue(aNodeId[1],aNodeId[2]);
1828       MValues::iterator aItr = theValues.find(aValue);
1829       if (aItr != theValues.end()){
1830         aItr->second += 1;
1831         //aNbConnects = nb;
1832       }
1833       else {
1834         theValues[aValue] = 1;
1835         //aNbConnects = 1;
1836       }
1837       //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1838       aNodeId[1] = aNodeId[2];
1839       aNode1 = aNode2;
1840     }
1841     Value aValue(aNodeId[0],aNodeId[2]);
1842     MValues::iterator aItr = theValues.find(aValue);
1843     if (aItr != theValues.end()) {
1844       aItr->second += 1;
1845       //aNbConnects = nb;
1846     }
1847     else {
1848       theValues[aValue] = 1;
1849       //aNbConnects = 1;
1850     }
1851     //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1852   }
1853
1854 }
1855
1856 /*
1857                             PREDICATES
1858 */
1859
1860 /*
1861   Class       : BadOrientedVolume
1862   Description : Predicate bad oriented volumes
1863 */
1864
1865 BadOrientedVolume::BadOrientedVolume()
1866 {
1867   myMesh = 0;
1868 }
1869
1870 void BadOrientedVolume::SetMesh( const SMDS_Mesh* theMesh )
1871 {
1872   myMesh = theMesh;
1873 }
1874
1875 bool BadOrientedVolume::IsSatisfy( long theId )
1876 {
1877   if ( myMesh == 0 )
1878     return false;
1879
1880   SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
1881   return !vTool.IsForward();
1882 }
1883
1884 SMDSAbs_ElementType BadOrientedVolume::GetType() const
1885 {
1886   return SMDSAbs_Volume;
1887 }
1888
1889 /*
1890   Class       : BareBorderVolume
1891 */
1892
1893 bool BareBorderVolume::IsSatisfy(long theElementId )
1894 {
1895   SMDS_VolumeTool  myTool;
1896   if ( myTool.Set( myMesh->FindElement(theElementId)))
1897   {
1898     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
1899       if ( myTool.IsFreeFace( iF ))
1900       {
1901         const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
1902         vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
1903         if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
1904           return true;
1905       }
1906   }
1907   return false;
1908 }
1909
1910 /*
1911   Class       : BareBorderFace
1912 */
1913
1914 bool BareBorderFace::IsSatisfy(long theElementId )
1915 {
1916   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
1917     if ( face->GetType() == SMDSAbs_Face )
1918     {
1919       int nbN = face->NbCornerNodes();
1920       for ( int i = 0; i < nbN; ++i )
1921       {
1922         // check if a link is shared by another face
1923         const SMDS_MeshNode* n1 = face->GetNode( i );
1924         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
1925         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
1926         bool isShared = false;
1927         while ( !isShared && fIt->more() )
1928         {
1929           const SMDS_MeshElement* f = fIt->next();
1930           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
1931         }
1932         if ( !isShared )
1933         {
1934           myLinkNodes.resize( 2 + face->IsQuadratic());
1935           myLinkNodes[0] = n1;
1936           myLinkNodes[1] = n2;
1937           if ( face->IsQuadratic() )
1938             myLinkNodes[2] = face->GetNode( i+nbN );
1939           return !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
1940         }
1941       }
1942     }
1943   return false;
1944 }
1945
1946 /*
1947   Class       : OverConstrainedVolume
1948 */
1949
1950 bool OverConstrainedVolume::IsSatisfy(long theElementId )
1951 {
1952   // An element is over-constrained if it has N-1 free borders where
1953   // N is the number of edges/faces for a 2D/3D element.
1954   SMDS_VolumeTool  myTool;
1955   if ( myTool.Set( myMesh->FindElement(theElementId)))
1956   {
1957     int nbSharedFaces = 0;
1958     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
1959       if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
1960         break;
1961     return ( nbSharedFaces == 1 );
1962   }
1963   return false;
1964 }
1965
1966 /*
1967   Class       : OverConstrainedFace
1968 */
1969
1970 bool OverConstrainedFace::IsSatisfy(long theElementId )
1971 {
1972   // An element is over-constrained if it has N-1 free borders where
1973   // N is the number of edges/faces for a 2D/3D element.
1974   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
1975     if ( face->GetType() == SMDSAbs_Face )
1976     {
1977       int nbSharedBorders = 0;
1978       int nbN = face->NbCornerNodes();
1979       for ( int i = 0; i < nbN; ++i )
1980       {
1981         // check if a link is shared by another face
1982         const SMDS_MeshNode* n1 = face->GetNode( i );
1983         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
1984         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
1985         bool isShared = false;
1986         while ( !isShared && fIt->more() )
1987         {
1988           const SMDS_MeshElement* f = fIt->next();
1989           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
1990         }
1991         if ( isShared && ++nbSharedBorders > 1 )
1992           break;
1993       }
1994       return ( nbSharedBorders == 1 );
1995     }
1996   return false;
1997 }
1998
1999 /*
2000   Class       : FreeBorders
2001   Description : Predicate for free borders
2002 */
2003
2004 FreeBorders::FreeBorders()
2005 {
2006   myMesh = 0;
2007 }
2008
2009 void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
2010 {
2011   myMesh = theMesh;
2012 }
2013
2014 bool FreeBorders::IsSatisfy( long theId )
2015 {
2016   return getNbMultiConnection( myMesh, theId ) == 1;
2017 }
2018
2019 SMDSAbs_ElementType FreeBorders::GetType() const
2020 {
2021   return SMDSAbs_Edge;
2022 }
2023
2024
2025 /*
2026   Class       : FreeEdges
2027   Description : Predicate for free Edges
2028 */
2029 FreeEdges::FreeEdges()
2030 {
2031   myMesh = 0;
2032 }
2033
2034 void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
2035 {
2036   myMesh = theMesh;
2037 }
2038
2039 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId  )
2040 {
2041   TColStd_MapOfInteger aMap;
2042   for ( int i = 0; i < 2; i++ )
2043   {
2044     SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator();
2045     while( anElemIter->more() )
2046     {
2047       const SMDS_MeshElement* anElem = anElemIter->next();
2048       if ( anElem != 0 && anElem->GetType() == SMDSAbs_Face )
2049       {
2050         int anId = anElem->GetID();
2051
2052         if ( i == 0 )
2053           aMap.Add( anId );
2054         else if ( aMap.Contains( anId ) && anId != theFaceId )
2055           return false;
2056       }
2057     }
2058   }
2059   return true;
2060 }
2061
2062 bool FreeEdges::IsSatisfy( long theId )
2063 {
2064   if ( myMesh == 0 )
2065     return false;
2066
2067   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2068   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
2069     return false;
2070
2071   SMDS_ElemIteratorPtr anIter;
2072   if ( aFace->IsQuadratic() ) {
2073     anIter = static_cast<const SMDS_QuadraticFaceOfNodes*>
2074       (aFace)->interlacedNodesElemIterator();
2075   }
2076   else {
2077     anIter = aFace->nodesIterator();
2078   }
2079   if ( anIter == 0 )
2080     return false;
2081
2082   int i = 0, nbNodes = aFace->NbNodes();
2083   std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
2084   while( anIter->more() )
2085   {
2086     const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
2087     if ( aNode == 0 )
2088       return false;
2089     aNodes[ i++ ] = aNode;
2090   }
2091   aNodes[ nbNodes ] = aNodes[ 0 ];
2092
2093   for ( i = 0; i < nbNodes; i++ )
2094     if ( IsFreeEdge( &aNodes[ i ], theId ) )
2095       return true;
2096
2097   return false;
2098 }
2099
2100 SMDSAbs_ElementType FreeEdges::GetType() const
2101 {
2102   return SMDSAbs_Face;
2103 }
2104
2105 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
2106   myElemId(theElemId)
2107 {
2108   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
2109   if(thePntId1 > thePntId2){
2110     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
2111   }
2112 }
2113
2114 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
2115   if(myPntId[0] < x.myPntId[0]) return true;
2116   if(myPntId[0] == x.myPntId[0])
2117     if(myPntId[1] < x.myPntId[1]) return true;
2118   return false;
2119 }
2120
2121 inline void UpdateBorders(const FreeEdges::Border& theBorder,
2122                           FreeEdges::TBorders& theRegistry,
2123                           FreeEdges::TBorders& theContainer)
2124 {
2125   if(theRegistry.find(theBorder) == theRegistry.end()){
2126     theRegistry.insert(theBorder);
2127     theContainer.insert(theBorder);
2128   }else{
2129     theContainer.erase(theBorder);
2130   }
2131 }
2132
2133 void FreeEdges::GetBoreders(TBorders& theBorders)
2134 {
2135   TBorders aRegistry;
2136   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
2137   for(; anIter->more(); ){
2138     const SMDS_MeshFace* anElem = anIter->next();
2139     long anElemId = anElem->GetID();
2140     SMDS_ElemIteratorPtr aNodesIter;
2141     if ( anElem->IsQuadratic() )
2142       aNodesIter = static_cast<const SMDS_QuadraticFaceOfNodes*>(anElem)->
2143         interlacedNodesElemIterator();
2144     else
2145       aNodesIter = anElem->nodesIterator();
2146     long aNodeId[2];
2147     const SMDS_MeshElement* aNode;
2148     if(aNodesIter->more()){
2149       aNode = aNodesIter->next();
2150       aNodeId[0] = aNodeId[1] = aNode->GetID();
2151     }
2152     for(; aNodesIter->more(); ){
2153       aNode = aNodesIter->next();
2154       long anId = aNode->GetID();
2155       Border aBorder(anElemId,aNodeId[1],anId);
2156       aNodeId[1] = anId;
2157       //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
2158       UpdateBorders(aBorder,aRegistry,theBorders);
2159     }
2160     Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
2161     //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
2162     UpdateBorders(aBorder,aRegistry,theBorders);
2163   }
2164   //std::cout<<"theBorders.size() = "<<theBorders.size()<<endl;
2165 }
2166
2167
2168 /*
2169   Class       : FreeNodes
2170   Description : Predicate for free nodes
2171 */
2172
2173 FreeNodes::FreeNodes()
2174 {
2175   myMesh = 0;
2176 }
2177
2178 void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
2179 {
2180   myMesh = theMesh;
2181 }
2182
2183 bool FreeNodes::IsSatisfy( long theNodeId )
2184 {
2185   const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
2186   if (!aNode)
2187     return false;
2188
2189   return (aNode->NbInverseElements() < 1);
2190 }
2191
2192 SMDSAbs_ElementType FreeNodes::GetType() const
2193 {
2194   return SMDSAbs_Node;
2195 }
2196
2197
2198 /*
2199   Class       : FreeFaces
2200   Description : Predicate for free faces
2201 */
2202
2203 FreeFaces::FreeFaces()
2204 {
2205   myMesh = 0;
2206 }
2207
2208 void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
2209 {
2210   myMesh = theMesh;
2211 }
2212
2213 bool FreeFaces::IsSatisfy( long theId )
2214 {
2215   if (!myMesh) return false;
2216   // check that faces nodes refers to less than two common volumes
2217   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2218   if ( !aFace || aFace->GetType() != SMDSAbs_Face )
2219     return false;
2220
2221   int nbNode = aFace->NbNodes();
2222
2223   // collect volumes check that number of volumss with count equal nbNode not less than 2
2224   typedef map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
2225   typedef map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
2226   TMapOfVolume mapOfVol;
2227
2228   SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
2229   while ( nodeItr->more() ) {
2230     const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
2231     if ( !aNode ) continue;
2232     SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
2233     while ( volItr->more() ) {
2234       SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
2235       TItrMapOfVolume itr = mapOfVol.insert(make_pair(aVol, 0)).first;
2236       (*itr).second++;
2237     }
2238   }
2239   int nbVol = 0;
2240   TItrMapOfVolume volItr = mapOfVol.begin();
2241   TItrMapOfVolume volEnd = mapOfVol.end();
2242   for ( ; volItr != volEnd; ++volItr )
2243     if ( (*volItr).second >= nbNode )
2244        nbVol++;
2245   // face is not free if number of volumes constructed on thier nodes more than one
2246   return (nbVol < 2);
2247 }
2248
2249 SMDSAbs_ElementType FreeFaces::GetType() const
2250 {
2251   return SMDSAbs_Face;
2252 }
2253
2254 /*
2255   Class       : LinearOrQuadratic
2256   Description : Predicate to verify whether a mesh element is linear
2257 */
2258
2259 LinearOrQuadratic::LinearOrQuadratic()
2260 {
2261   myMesh = 0;
2262 }
2263
2264 void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
2265 {
2266   myMesh = theMesh;
2267 }
2268
2269 bool LinearOrQuadratic::IsSatisfy( long theId )
2270 {
2271   if (!myMesh) return false;
2272   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2273   if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
2274     return false;
2275   return (!anElem->IsQuadratic());
2276 }
2277
2278 void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
2279 {
2280   myType = theType;
2281 }
2282
2283 SMDSAbs_ElementType LinearOrQuadratic::GetType() const
2284 {
2285   return myType;
2286 }
2287
2288 /*
2289   Class       : GroupColor
2290   Description : Functor for check color of group to whic mesh element belongs to
2291 */
2292
2293 GroupColor::GroupColor()
2294 {
2295 }
2296
2297 bool GroupColor::IsSatisfy( long theId )
2298 {
2299   return (myIDs.find( theId ) != myIDs.end());
2300 }
2301
2302 void GroupColor::SetType( SMDSAbs_ElementType theType )
2303 {
2304   myType = theType;
2305 }
2306
2307 SMDSAbs_ElementType GroupColor::GetType() const
2308 {
2309   return myType;
2310 }
2311
2312 static bool isEqual( const Quantity_Color& theColor1,
2313                      const Quantity_Color& theColor2 )
2314 {
2315   // tolerance to compare colors
2316   const double tol = 5*1e-3;
2317   return ( fabs( theColor1.Red() - theColor2.Red() ) < tol &&
2318            fabs( theColor1.Green() - theColor2.Green() ) < tol &&
2319            fabs( theColor1.Blue() - theColor2.Blue() ) < tol );
2320 }
2321
2322
2323 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
2324 {
2325   myIDs.clear();
2326
2327   const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
2328   if ( !aMesh )
2329     return;
2330
2331   int nbGrp = aMesh->GetNbGroups();
2332   if ( !nbGrp )
2333     return;
2334
2335   // iterates on groups and find necessary elements ids
2336   const std::set<SMESHDS_GroupBase*>& aGroups = aMesh->GetGroups();
2337   set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
2338   for (; GrIt != aGroups.end(); GrIt++) {
2339     SMESHDS_GroupBase* aGrp = (*GrIt);
2340     if ( !aGrp )
2341       continue;
2342     // check type and color of group
2343     if ( !isEqual( myColor, aGrp->GetColor() ) )
2344       continue;
2345     if ( myType != SMDSAbs_All && myType != (SMDSAbs_ElementType)aGrp->GetType() )
2346       continue;
2347
2348     SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
2349     if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
2350       // add elements IDS into control
2351       int aSize = aGrp->Extent();
2352       for (int i = 0; i < aSize; i++)
2353         myIDs.insert( aGrp->GetID(i+1) );
2354     }
2355   }
2356 }
2357
2358 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
2359 {
2360   TCollection_AsciiString aStr = theStr;
2361   aStr.RemoveAll( ' ' );
2362   aStr.RemoveAll( '\t' );
2363   for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
2364     aStr.Remove( aPos, 2 );
2365   Standard_Real clr[3];
2366   clr[0] = clr[1] = clr[2] = 0.;
2367   for ( int i = 0; i < 3; i++ ) {
2368     TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
2369     if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
2370       clr[i] = tmpStr.RealValue();
2371   }
2372   myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
2373 }
2374
2375 //=======================================================================
2376 // name    : GetRangeStr
2377 // Purpose : Get range as a string.
2378 //           Example: "1,2,3,50-60,63,67,70-"
2379 //=======================================================================
2380 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
2381 {
2382   theResStr.Clear();
2383   theResStr += TCollection_AsciiString( myColor.Red() );
2384   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
2385   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
2386 }
2387
2388 /*
2389   Class       : ElemGeomType
2390   Description : Predicate to check element geometry type
2391 */
2392
2393 ElemGeomType::ElemGeomType()
2394 {
2395   myMesh = 0;
2396   myType = SMDSAbs_All;
2397   myGeomType = SMDSGeom_TRIANGLE;
2398 }
2399
2400 void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
2401 {
2402   myMesh = theMesh;
2403 }
2404
2405 bool ElemGeomType::IsSatisfy( long theId )
2406 {
2407   if (!myMesh) return false;
2408   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2409   if ( !anElem )
2410     return false;
2411   const SMDSAbs_ElementType anElemType = anElem->GetType();
2412   if ( myType != SMDSAbs_All && anElemType != myType )
2413     return false;
2414   const int aNbNode = anElem->NbNodes();
2415   bool isOk = false;
2416   switch( anElemType )
2417   {
2418   case SMDSAbs_Node:
2419     isOk = (myGeomType == SMDSGeom_POINT);
2420     break;
2421
2422   case SMDSAbs_Edge:
2423     isOk = (myGeomType == SMDSGeom_EDGE);
2424     break;
2425
2426   case SMDSAbs_Face:
2427     if ( myGeomType == SMDSGeom_TRIANGLE )
2428       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 6 : aNbNode == 3));
2429     else if ( myGeomType == SMDSGeom_QUADRANGLE )
2430       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 8 : aNbNode == 4));
2431     else if ( myGeomType == SMDSGeom_POLYGON )
2432       isOk = anElem->IsPoly();
2433     break;
2434
2435   case SMDSAbs_Volume:
2436     if ( myGeomType == SMDSGeom_TETRA )
2437       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 10 : aNbNode == 4));
2438     else if ( myGeomType == SMDSGeom_PYRAMID )
2439       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 13 : aNbNode == 5));
2440     else if ( myGeomType == SMDSGeom_PENTA )
2441       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 15 : aNbNode == 6));
2442     else if ( myGeomType == SMDSGeom_HEXA )
2443       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 20 : aNbNode == 8));
2444      else if ( myGeomType == SMDSGeom_POLYHEDRA )
2445       isOk = anElem->IsPoly();
2446     break;
2447     default: break;
2448   }
2449   return isOk;
2450 }
2451
2452 void ElemGeomType::SetType( SMDSAbs_ElementType theType )
2453 {
2454   myType = theType;
2455 }
2456
2457 SMDSAbs_ElementType ElemGeomType::GetType() const
2458 {
2459   return myType;
2460 }
2461
2462 void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
2463 {
2464   myGeomType = theType;
2465 }
2466
2467 SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
2468 {
2469   return myGeomType;
2470 }
2471
2472 /*
2473   Class       : RangeOfIds
2474   Description : Predicate for Range of Ids.
2475                 Range may be specified with two ways.
2476                 1. Using AddToRange method
2477                 2. With SetRangeStr method. Parameter of this method is a string
2478                    like as "1,2,3,50-60,63,67,70-"
2479 */
2480
2481 //=======================================================================
2482 // name    : RangeOfIds
2483 // Purpose : Constructor
2484 //=======================================================================
2485 RangeOfIds::RangeOfIds()
2486 {
2487   myMesh = 0;
2488   myType = SMDSAbs_All;
2489 }
2490
2491 //=======================================================================
2492 // name    : SetMesh
2493 // Purpose : Set mesh
2494 //=======================================================================
2495 void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh )
2496 {
2497   myMesh = theMesh;
2498 }
2499
2500 //=======================================================================
2501 // name    : AddToRange
2502 // Purpose : Add ID to the range
2503 //=======================================================================
2504 bool RangeOfIds::AddToRange( long theEntityId )
2505 {
2506   myIds.Add( theEntityId );
2507   return true;
2508 }
2509
2510 //=======================================================================
2511 // name    : GetRangeStr
2512 // Purpose : Get range as a string.
2513 //           Example: "1,2,3,50-60,63,67,70-"
2514 //=======================================================================
2515 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
2516 {
2517   theResStr.Clear();
2518
2519   TColStd_SequenceOfInteger     anIntSeq;
2520   TColStd_SequenceOfAsciiString aStrSeq;
2521
2522   TColStd_MapIteratorOfMapOfInteger anIter( myIds );
2523   for ( ; anIter.More(); anIter.Next() )
2524   {
2525     int anId = anIter.Key();
2526     TCollection_AsciiString aStr( anId );
2527     anIntSeq.Append( anId );
2528     aStrSeq.Append( aStr );
2529   }
2530
2531   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
2532   {
2533     int aMinId = myMin( i );
2534     int aMaxId = myMax( i );
2535
2536     TCollection_AsciiString aStr;
2537     if ( aMinId != IntegerFirst() )
2538       aStr += aMinId;
2539
2540     aStr += "-";
2541
2542     if ( aMaxId != IntegerLast() )
2543       aStr += aMaxId;
2544
2545     // find position of the string in result sequence and insert string in it
2546     if ( anIntSeq.Length() == 0 )
2547     {
2548       anIntSeq.Append( aMinId );
2549       aStrSeq.Append( aStr );
2550     }
2551     else
2552     {
2553       if ( aMinId < anIntSeq.First() )
2554       {
2555         anIntSeq.Prepend( aMinId );
2556         aStrSeq.Prepend( aStr );
2557       }
2558       else if ( aMinId > anIntSeq.Last() )
2559       {
2560         anIntSeq.Append( aMinId );
2561         aStrSeq.Append( aStr );
2562       }
2563       else
2564         for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
2565           if ( aMinId < anIntSeq( j ) )
2566           {
2567             anIntSeq.InsertBefore( j, aMinId );
2568             aStrSeq.InsertBefore( j, aStr );
2569             break;
2570           }
2571     }
2572   }
2573
2574   if ( aStrSeq.Length() == 0 )
2575     return;
2576
2577   theResStr = aStrSeq( 1 );
2578   for ( int j = 2, k = aStrSeq.Length(); j <= k; j++  )
2579   {
2580     theResStr += ",";
2581     theResStr += aStrSeq( j );
2582   }
2583 }
2584
2585 //=======================================================================
2586 // name    : SetRangeStr
2587 // Purpose : Define range with string
2588 //           Example of entry string: "1,2,3,50-60,63,67,70-"
2589 //=======================================================================
2590 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
2591 {
2592   myMin.Clear();
2593   myMax.Clear();
2594   myIds.Clear();
2595
2596   TCollection_AsciiString aStr = theStr;
2597   aStr.RemoveAll( ' ' );
2598   aStr.RemoveAll( '\t' );
2599
2600   for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
2601     aStr.Remove( aPos, 2 );
2602
2603   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
2604   int i = 1;
2605   while ( tmpStr != "" )
2606   {
2607     tmpStr = aStr.Token( ",", i++ );
2608     int aPos = tmpStr.Search( '-' );
2609
2610     if ( aPos == -1 )
2611     {
2612       if ( tmpStr.IsIntegerValue() )
2613         myIds.Add( tmpStr.IntegerValue() );
2614       else
2615         return false;
2616     }
2617     else
2618     {
2619       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
2620       TCollection_AsciiString aMinStr = tmpStr;
2621
2622       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
2623       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
2624
2625       if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() ||
2626            !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() )
2627         return false;
2628
2629       myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
2630       myMax.Append( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
2631     }
2632   }
2633
2634   return true;
2635 }
2636
2637 //=======================================================================
2638 // name    : GetType
2639 // Purpose : Get type of supported entities
2640 //=======================================================================
2641 SMDSAbs_ElementType RangeOfIds::GetType() const
2642 {
2643   return myType;
2644 }
2645
2646 //=======================================================================
2647 // name    : SetType
2648 // Purpose : Set type of supported entities
2649 //=======================================================================
2650 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
2651 {
2652   myType = theType;
2653 }
2654
2655 //=======================================================================
2656 // name    : IsSatisfy
2657 // Purpose : Verify whether entity satisfies to this rpedicate
2658 //=======================================================================
2659 bool RangeOfIds::IsSatisfy( long theId )
2660 {
2661   if ( !myMesh )
2662     return false;
2663
2664   if ( myType == SMDSAbs_Node )
2665   {
2666     if ( myMesh->FindNode( theId ) == 0 )
2667       return false;
2668   }
2669   else
2670   {
2671     const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2672     if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All )
2673       return false;
2674   }
2675
2676   if ( myIds.Contains( theId ) )
2677     return true;
2678
2679   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
2680     if ( theId >= myMin( i ) && theId <= myMax( i ) )
2681       return true;
2682
2683   return false;
2684 }
2685
2686 /*
2687   Class       : Comparator
2688   Description : Base class for comparators
2689 */
2690 Comparator::Comparator():
2691   myMargin(0)
2692 {}
2693
2694 Comparator::~Comparator()
2695 {}
2696
2697 void Comparator::SetMesh( const SMDS_Mesh* theMesh )
2698 {
2699   if ( myFunctor )
2700     myFunctor->SetMesh( theMesh );
2701 }
2702
2703 void Comparator::SetMargin( double theValue )
2704 {
2705   myMargin = theValue;
2706 }
2707
2708 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
2709 {
2710   myFunctor = theFunct;
2711 }
2712
2713 SMDSAbs_ElementType Comparator::GetType() const
2714 {
2715   return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
2716 }
2717
2718 double Comparator::GetMargin()
2719 {
2720   return myMargin;
2721 }
2722
2723
2724 /*
2725   Class       : LessThan
2726   Description : Comparator "<"
2727 */
2728 bool LessThan::IsSatisfy( long theId )
2729 {
2730   return myFunctor && myFunctor->GetValue( theId ) < myMargin;
2731 }
2732
2733
2734 /*
2735   Class       : MoreThan
2736   Description : Comparator ">"
2737 */
2738 bool MoreThan::IsSatisfy( long theId )
2739 {
2740   return myFunctor && myFunctor->GetValue( theId ) > myMargin;
2741 }
2742
2743
2744 /*
2745   Class       : EqualTo
2746   Description : Comparator "="
2747 */
2748 EqualTo::EqualTo():
2749   myToler(Precision::Confusion())
2750 {}
2751
2752 bool EqualTo::IsSatisfy( long theId )
2753 {
2754   return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
2755 }
2756
2757 void EqualTo::SetTolerance( double theToler )
2758 {
2759   myToler = theToler;
2760 }
2761
2762 double EqualTo::GetTolerance()
2763 {
2764   return myToler;
2765 }
2766
2767 /*
2768   Class       : LogicalNOT
2769   Description : Logical NOT predicate
2770 */
2771 LogicalNOT::LogicalNOT()
2772 {}
2773
2774 LogicalNOT::~LogicalNOT()
2775 {}
2776
2777 bool LogicalNOT::IsSatisfy( long theId )
2778 {
2779   return myPredicate && !myPredicate->IsSatisfy( theId );
2780 }
2781
2782 void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
2783 {
2784   if ( myPredicate )
2785     myPredicate->SetMesh( theMesh );
2786 }
2787
2788 void LogicalNOT::SetPredicate( PredicatePtr thePred )
2789 {
2790   myPredicate = thePred;
2791 }
2792
2793 SMDSAbs_ElementType LogicalNOT::GetType() const
2794 {
2795   return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
2796 }
2797
2798
2799 /*
2800   Class       : LogicalBinary
2801   Description : Base class for binary logical predicate
2802 */
2803 LogicalBinary::LogicalBinary()
2804 {}
2805
2806 LogicalBinary::~LogicalBinary()
2807 {}
2808
2809 void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
2810 {
2811   if ( myPredicate1 )
2812     myPredicate1->SetMesh( theMesh );
2813
2814   if ( myPredicate2 )
2815     myPredicate2->SetMesh( theMesh );
2816 }
2817
2818 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
2819 {
2820   myPredicate1 = thePredicate;
2821 }
2822
2823 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
2824 {
2825   myPredicate2 = thePredicate;
2826 }
2827
2828 SMDSAbs_ElementType LogicalBinary::GetType() const
2829 {
2830   if ( !myPredicate1 || !myPredicate2 )
2831     return SMDSAbs_All;
2832
2833   SMDSAbs_ElementType aType1 = myPredicate1->GetType();
2834   SMDSAbs_ElementType aType2 = myPredicate2->GetType();
2835
2836   return aType1 == aType2 ? aType1 : SMDSAbs_All;
2837 }
2838
2839
2840 /*
2841   Class       : LogicalAND
2842   Description : Logical AND
2843 */
2844 bool LogicalAND::IsSatisfy( long theId )
2845 {
2846   return
2847     myPredicate1 &&
2848     myPredicate2 &&
2849     myPredicate1->IsSatisfy( theId ) &&
2850     myPredicate2->IsSatisfy( theId );
2851 }
2852
2853
2854 /*
2855   Class       : LogicalOR
2856   Description : Logical OR
2857 */
2858 bool LogicalOR::IsSatisfy( long theId )
2859 {
2860   return
2861     myPredicate1 &&
2862     myPredicate2 &&
2863     myPredicate1->IsSatisfy( theId ) ||
2864     myPredicate2->IsSatisfy( theId );
2865 }
2866
2867
2868 /*
2869                               FILTER
2870 */
2871
2872 Filter::Filter()
2873 {}
2874
2875 Filter::~Filter()
2876 {}
2877
2878 void Filter::SetPredicate( PredicatePtr thePredicate )
2879 {
2880   myPredicate = thePredicate;
2881 }
2882
2883 template<class TElement, class TIterator, class TPredicate>
2884 inline void FillSequence(const TIterator& theIterator,
2885                          TPredicate& thePredicate,
2886                          Filter::TIdSequence& theSequence)
2887 {
2888   if ( theIterator ) {
2889     while( theIterator->more() ) {
2890       TElement anElem = theIterator->next();
2891       long anId = anElem->GetID();
2892       if ( thePredicate->IsSatisfy( anId ) )
2893         theSequence.push_back( anId );
2894     }
2895   }
2896 }
2897
2898 void
2899 Filter::
2900 GetElementsId( const SMDS_Mesh* theMesh,
2901                PredicatePtr thePredicate,
2902                TIdSequence& theSequence )
2903 {
2904   theSequence.clear();
2905
2906   if ( !theMesh || !thePredicate )
2907     return;
2908
2909   thePredicate->SetMesh( theMesh );
2910
2911   SMDSAbs_ElementType aType = thePredicate->GetType();
2912   switch(aType){
2913   case SMDSAbs_Node:
2914     FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),thePredicate,theSequence);
2915     break;
2916   case SMDSAbs_Edge:
2917     FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
2918     break;
2919   case SMDSAbs_Face:
2920     FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
2921     break;
2922   case SMDSAbs_Volume:
2923     FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),thePredicate,theSequence);
2924     break;
2925   case SMDSAbs_All:
2926     FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
2927     FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
2928     FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),thePredicate,theSequence);
2929     break;
2930   }
2931 }
2932
2933 void
2934 Filter::GetElementsId( const SMDS_Mesh* theMesh,
2935                        Filter::TIdSequence& theSequence )
2936 {
2937   GetElementsId(theMesh,myPredicate,theSequence);
2938 }
2939
2940 /*
2941                               ManifoldPart
2942 */
2943
2944 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
2945
2946 /*
2947    Internal class Link
2948 */
2949
2950 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
2951                           SMDS_MeshNode* theNode2 )
2952 {
2953   myNode1 = theNode1;
2954   myNode2 = theNode2;
2955 }
2956
2957 ManifoldPart::Link::~Link()
2958 {
2959   myNode1 = 0;
2960   myNode2 = 0;
2961 }
2962
2963 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
2964 {
2965   if ( myNode1 == theLink.myNode1 &&
2966        myNode2 == theLink.myNode2 )
2967     return true;
2968   else if ( myNode1 == theLink.myNode2 &&
2969             myNode2 == theLink.myNode1 )
2970     return true;
2971   else
2972     return false;
2973 }
2974
2975 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
2976 {
2977   if(myNode1 < x.myNode1) return true;
2978   if(myNode1 == x.myNode1)
2979     if(myNode2 < x.myNode2) return true;
2980   return false;
2981 }
2982
2983 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
2984                             const ManifoldPart::Link& theLink2 )
2985 {
2986   return theLink1.IsEqual( theLink2 );
2987 }
2988
2989 ManifoldPart::ManifoldPart()
2990 {
2991   myMesh = 0;
2992   myAngToler = Precision::Angular();
2993   myIsOnlyManifold = true;
2994 }
2995
2996 ManifoldPart::~ManifoldPart()
2997 {
2998   myMesh = 0;
2999 }
3000
3001 void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
3002 {
3003   myMesh = theMesh;
3004   process();
3005 }
3006
3007 SMDSAbs_ElementType ManifoldPart::GetType() const
3008 { return SMDSAbs_Face; }
3009
3010 bool ManifoldPart::IsSatisfy( long theElementId )
3011 {
3012   return myMapIds.Contains( theElementId );
3013 }
3014
3015 void ManifoldPart::SetAngleTolerance( const double theAngToler )
3016 { myAngToler = theAngToler; }
3017
3018 double ManifoldPart::GetAngleTolerance() const
3019 { return myAngToler; }
3020
3021 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
3022 { myIsOnlyManifold = theIsOnly; }
3023
3024 void ManifoldPart::SetStartElem( const long  theStartId )
3025 { myStartElemId = theStartId; }
3026
3027 bool ManifoldPart::process()
3028 {
3029   myMapIds.Clear();
3030   myMapBadGeomIds.Clear();
3031
3032   myAllFacePtr.clear();
3033   myAllFacePtrIntDMap.clear();
3034   if ( !myMesh )
3035     return false;
3036
3037   // collect all faces into own map
3038   SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
3039   for (; anFaceItr->more(); )
3040   {
3041     SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
3042     myAllFacePtr.push_back( aFacePtr );
3043     myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
3044   }
3045
3046   SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
3047   if ( !aStartFace )
3048     return false;
3049
3050   // the map of non manifold links and bad geometry
3051   TMapOfLink aMapOfNonManifold;
3052   TColStd_MapOfInteger aMapOfTreated;
3053
3054   // begin cycle on faces from start index and run on vector till the end
3055   //  and from begin to start index to cover whole vector
3056   const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
3057   bool isStartTreat = false;
3058   for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
3059   {
3060     if ( fi == aStartIndx )
3061       isStartTreat = true;
3062     // as result next time when fi will be equal to aStartIndx
3063
3064     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
3065     if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
3066       continue;
3067
3068     aMapOfTreated.Add( aFacePtr->GetID() );
3069     TColStd_MapOfInteger aResFaces;
3070     if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
3071                          aMapOfNonManifold, aResFaces ) )
3072       continue;
3073     TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
3074     for ( ; anItr.More(); anItr.Next() )
3075     {
3076       int aFaceId = anItr.Key();
3077       aMapOfTreated.Add( aFaceId );
3078       myMapIds.Add( aFaceId );
3079     }
3080
3081     if ( fi == ( myAllFacePtr.size() - 1 ) )
3082       fi = 0;
3083   } // end run on vector of faces
3084   return !myMapIds.IsEmpty();
3085 }
3086
3087 static void getLinks( const SMDS_MeshFace* theFace,
3088                       ManifoldPart::TVectorOfLink& theLinks )
3089 {
3090   int aNbNode = theFace->NbNodes();
3091   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
3092   int i = 1;
3093   SMDS_MeshNode* aNode = 0;
3094   for ( ; aNodeItr->more() && i <= aNbNode; )
3095   {
3096
3097     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
3098     if ( i == 1 )
3099       aNode = aN1;
3100     i++;
3101     SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
3102     i++;
3103     ManifoldPart::Link aLink( aN1, aN2 );
3104     theLinks.push_back( aLink );
3105   }
3106 }
3107
3108 static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
3109 {
3110   gp_XYZ n;
3111   int aNbNode = theFace->NbNodes();
3112   TColgp_Array1OfXYZ anArrOfXYZ(1,4);
3113   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
3114   int i = 1;
3115   for ( ; aNodeItr->more() && i <= 4; i++ ) {
3116     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
3117     anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
3118   }
3119
3120   gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1);
3121   gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1);
3122   n  = q1 ^ q2;
3123   if ( aNbNode > 3 ) {
3124     gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
3125     n += q2 ^ q3;
3126   }
3127   double len = n.Modulus();
3128   if ( len > 0 )
3129     n /= len;
3130
3131   return n;
3132 }
3133
3134 bool ManifoldPart::findConnected
3135                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
3136                   SMDS_MeshFace*                           theStartFace,
3137                   ManifoldPart::TMapOfLink&                theNonManifold,
3138                   TColStd_MapOfInteger&                    theResFaces )
3139 {
3140   theResFaces.Clear();
3141   if ( !theAllFacePtrInt.size() )
3142     return false;
3143
3144   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
3145   {
3146     myMapBadGeomIds.Add( theStartFace->GetID() );
3147     return false;
3148   }
3149
3150   ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
3151   ManifoldPart::TVectorOfLink aSeqOfBoundary;
3152   theResFaces.Add( theStartFace->GetID() );
3153   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
3154
3155   expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3156                  aDMapLinkFace, theNonManifold, theStartFace );
3157
3158   bool isDone = false;
3159   while ( !isDone && aMapOfBoundary.size() != 0 )
3160   {
3161     bool isToReset = false;
3162     ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
3163     for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
3164     {
3165       ManifoldPart::Link aLink = *pLink;
3166       if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
3167         continue;
3168       // each link could be treated only once
3169       aMapToSkip.insert( aLink );
3170
3171       ManifoldPart::TVectorOfFacePtr aFaces;
3172       // find next
3173       if ( myIsOnlyManifold &&
3174            (theNonManifold.find( aLink ) != theNonManifold.end()) )
3175         continue;
3176       else
3177       {
3178         getFacesByLink( aLink, aFaces );
3179         // filter the element to keep only indicated elements
3180         ManifoldPart::TVectorOfFacePtr aFiltered;
3181         ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3182         for ( ; pFace != aFaces.end(); ++pFace )
3183         {
3184           SMDS_MeshFace* aFace = *pFace;
3185           if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
3186             aFiltered.push_back( aFace );
3187         }
3188         aFaces = aFiltered;
3189         if ( aFaces.size() < 2 )  // no neihgbour faces
3190           continue;
3191         else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
3192         {
3193           theNonManifold.insert( aLink );
3194           continue;
3195         }
3196       }
3197
3198       // compare normal with normals of neighbor element
3199       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
3200       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3201       for ( ; pFace != aFaces.end(); ++pFace )
3202       {
3203         SMDS_MeshFace* aNextFace = *pFace;
3204         if ( aPrevFace == aNextFace )
3205           continue;
3206         int anNextFaceID = aNextFace->GetID();
3207         if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
3208          // should not be with non manifold restriction. probably bad topology
3209           continue;
3210         // check if face was treated and skipped
3211         if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
3212              !isInPlane( aPrevFace, aNextFace ) )
3213           continue;
3214         // add new element to connected and extend the boundaries.
3215         theResFaces.Add( anNextFaceID );
3216         expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3217                         aDMapLinkFace, theNonManifold, aNextFace );
3218         isToReset = true;
3219       }
3220     }
3221     isDone = !isToReset;
3222   }
3223
3224   return !theResFaces.IsEmpty();
3225 }
3226
3227 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
3228                               const SMDS_MeshFace* theFace2 )
3229 {
3230   gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
3231   gp_XYZ aNorm2XYZ = getNormale( theFace2 );
3232   if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
3233   {
3234     myMapBadGeomIds.Add( theFace2->GetID() );
3235     return false;
3236   }
3237   if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
3238     return true;
3239
3240   return false;
3241 }
3242
3243 void ManifoldPart::expandBoundary
3244                    ( ManifoldPart::TMapOfLink&            theMapOfBoundary,
3245                      ManifoldPart::TVectorOfLink&         theSeqOfBoundary,
3246                      ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
3247                      ManifoldPart::TMapOfLink&            theNonManifold,
3248                      SMDS_MeshFace*                       theNextFace ) const
3249 {
3250   ManifoldPart::TVectorOfLink aLinks;
3251   getLinks( theNextFace, aLinks );
3252   int aNbLink = (int)aLinks.size();
3253   for ( int i = 0; i < aNbLink; i++ )
3254   {
3255     ManifoldPart::Link aLink = aLinks[ i ];
3256     if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
3257       continue;
3258     if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
3259     {
3260       if ( myIsOnlyManifold )
3261       {
3262         // remove from boundary
3263         theMapOfBoundary.erase( aLink );
3264         ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
3265         for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
3266         {
3267           ManifoldPart::Link aBoundLink = *pLink;
3268           if ( aBoundLink.IsEqual( aLink ) )
3269           {
3270             theSeqOfBoundary.erase( pLink );
3271             break;
3272           }
3273         }
3274       }
3275     }
3276     else
3277     {
3278       theMapOfBoundary.insert( aLink );
3279       theSeqOfBoundary.push_back( aLink );
3280       theDMapLinkFacePtr[ aLink ] = theNextFace;
3281     }
3282   }
3283 }
3284
3285 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
3286                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
3287 {
3288   SMDS_Mesh::SetOfFaces aSetOfFaces;
3289   // take all faces that shared first node
3290   SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
3291   for ( ; anItr->more(); )
3292   {
3293     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3294     if ( !aFace )
3295       continue;
3296     aSetOfFaces.Add( aFace );
3297   }
3298   // take all faces that shared second node
3299   anItr = theLink.myNode2->facesIterator();
3300   // find the common part of two sets
3301   for ( ; anItr->more(); )
3302   {
3303     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3304     if ( aSetOfFaces.Contains( aFace ) )
3305       theFaces.push_back( aFace );
3306   }
3307 }
3308
3309
3310 /*
3311    ElementsOnSurface
3312 */
3313
3314 ElementsOnSurface::ElementsOnSurface()
3315 {
3316   myMesh = 0;
3317   myIds.Clear();
3318   myType = SMDSAbs_All;
3319   mySurf.Nullify();
3320   myToler = Precision::Confusion();
3321   myUseBoundaries = false;
3322 }
3323
3324 ElementsOnSurface::~ElementsOnSurface()
3325 {
3326   myMesh = 0;
3327 }
3328
3329 void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
3330 {
3331   if ( myMesh == theMesh )
3332     return;
3333   myMesh = theMesh;
3334   process();
3335 }
3336
3337 bool ElementsOnSurface::IsSatisfy( long theElementId )
3338 {
3339   return myIds.Contains( theElementId );
3340 }
3341
3342 SMDSAbs_ElementType ElementsOnSurface::GetType() const
3343 { return myType; }
3344
3345 void ElementsOnSurface::SetTolerance( const double theToler )
3346 {
3347   if ( myToler != theToler )
3348     myIds.Clear();
3349   myToler = theToler;
3350 }
3351
3352 double ElementsOnSurface::GetTolerance() const
3353 { return myToler; }
3354
3355 void ElementsOnSurface::SetUseBoundaries( bool theUse )
3356 {
3357   if ( myUseBoundaries != theUse ) {
3358     myUseBoundaries = theUse;
3359     SetSurface( mySurf, myType );
3360   }
3361 }
3362
3363 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
3364                                     const SMDSAbs_ElementType theType )
3365 {
3366   myIds.Clear();
3367   myType = theType;
3368   mySurf.Nullify();
3369   if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
3370     return;
3371   mySurf = TopoDS::Face( theShape );
3372   BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
3373   Standard_Real
3374     u1 = SA.FirstUParameter(),
3375     u2 = SA.LastUParameter(),
3376     v1 = SA.FirstVParameter(),
3377     v2 = SA.LastVParameter();
3378   Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
3379   myProjector.Init( surf, u1,u2, v1,v2 );
3380   process();
3381 }
3382
3383 void ElementsOnSurface::process()
3384 {
3385   myIds.Clear();
3386   if ( mySurf.IsNull() )
3387     return;
3388
3389   if ( myMesh == 0 )
3390     return;
3391
3392   if ( myType == SMDSAbs_Face || myType == SMDSAbs_All )
3393   {
3394     myIds.ReSize( myMesh->NbFaces() );
3395     SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
3396     for(; anIter->more(); )
3397       process( anIter->next() );
3398   }
3399
3400   if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All )
3401   {
3402     myIds.ReSize( myIds.Extent() + myMesh->NbEdges() );
3403     SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
3404     for(; anIter->more(); )
3405       process( anIter->next() );
3406   }
3407
3408   if ( myType == SMDSAbs_Node )
3409   {
3410     myIds.ReSize( myMesh->NbNodes() );
3411     SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
3412     for(; anIter->more(); )
3413       process( anIter->next() );
3414   }
3415 }
3416
3417 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
3418 {
3419   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
3420   bool isSatisfy = true;
3421   for ( ; aNodeItr->more(); )
3422   {
3423     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
3424     if ( !isOnSurface( aNode ) )
3425     {
3426       isSatisfy = false;
3427       break;
3428     }
3429   }
3430   if ( isSatisfy )
3431     myIds.Add( theElemPtr->GetID() );
3432 }
3433
3434 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
3435 {
3436   if ( mySurf.IsNull() )
3437     return false;
3438
3439   gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
3440   //  double aToler2 = myToler * myToler;
3441 //   if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
3442 //   {
3443 //     gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
3444 //     if ( aPln.SquareDistance( aPnt ) > aToler2 )
3445 //       return false;
3446 //   }
3447 //   else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
3448 //   {
3449 //     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
3450 //     double aRad = aCyl.Radius();
3451 //     gp_Ax3 anAxis = aCyl.Position();
3452 //     gp_XYZ aLoc = aCyl.Location().XYZ();
3453 //     double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
3454 //     double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
3455 //     if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
3456 //       return false;
3457 //   }
3458 //   else
3459 //     return false;
3460   myProjector.Perform( aPnt );
3461   bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
3462
3463   return isOn;
3464 }
3465
3466
3467 /*
3468   ElementsOnShape
3469 */
3470
3471 ElementsOnShape::ElementsOnShape()
3472   : myMesh(0),
3473     myType(SMDSAbs_All),
3474     myToler(Precision::Confusion()),
3475     myAllNodesFlag(false)
3476 {
3477   myCurShapeType = TopAbs_SHAPE;
3478 }
3479
3480 ElementsOnShape::~ElementsOnShape()
3481 {
3482 }
3483
3484 void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
3485 {
3486   if (myMesh != theMesh) {
3487     myMesh = theMesh;
3488     SetShape(myShape, myType);
3489   }
3490 }
3491
3492 bool ElementsOnShape::IsSatisfy (long theElementId)
3493 {
3494   return myIds.Contains(theElementId);
3495 }
3496
3497 SMDSAbs_ElementType ElementsOnShape::GetType() const
3498 {
3499   return myType;
3500 }
3501
3502 void ElementsOnShape::SetTolerance (const double theToler)
3503 {
3504   if (myToler != theToler) {
3505     myToler = theToler;
3506     SetShape(myShape, myType);
3507   }
3508 }
3509
3510 double ElementsOnShape::GetTolerance() const
3511 {
3512   return myToler;
3513 }
3514
3515 void ElementsOnShape::SetAllNodes (bool theAllNodes)
3516 {
3517   if (myAllNodesFlag != theAllNodes) {
3518     myAllNodesFlag = theAllNodes;
3519     SetShape(myShape, myType);
3520   }
3521 }
3522
3523 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
3524                                 const SMDSAbs_ElementType theType)
3525 {
3526   myType = theType;
3527   myShape = theShape;
3528   myIds.Clear();
3529
3530   if (myMesh == 0) return;
3531
3532   switch (myType)
3533   {
3534   case SMDSAbs_All:
3535     myIds.ReSize(myMesh->NbEdges() + myMesh->NbFaces() + myMesh->NbVolumes());
3536     break;
3537   case SMDSAbs_Node:
3538     myIds.ReSize(myMesh->NbNodes());
3539     break;
3540   case SMDSAbs_Edge:
3541     myIds.ReSize(myMesh->NbEdges());
3542     break;
3543   case SMDSAbs_Face:
3544     myIds.ReSize(myMesh->NbFaces());
3545     break;
3546   case SMDSAbs_Volume:
3547     myIds.ReSize(myMesh->NbVolumes());
3548     break;
3549   default:
3550     break;
3551   }
3552
3553   myShapesMap.Clear();
3554   addShape(myShape);
3555 }
3556
3557 void ElementsOnShape::addShape (const TopoDS_Shape& theShape)
3558 {
3559   if (theShape.IsNull() || myMesh == 0)
3560     return;
3561
3562   if (!myShapesMap.Add(theShape)) return;
3563
3564   myCurShapeType = theShape.ShapeType();
3565   switch (myCurShapeType)
3566   {
3567   case TopAbs_COMPOUND:
3568   case TopAbs_COMPSOLID:
3569   case TopAbs_SHELL:
3570   case TopAbs_WIRE:
3571     {
3572       TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
3573       for (; anIt.More(); anIt.Next()) addShape(anIt.Value());
3574     }
3575     break;
3576   case TopAbs_SOLID:
3577     {
3578       myCurSC.Load(theShape);
3579       process();
3580     }
3581     break;
3582   case TopAbs_FACE:
3583     {
3584       TopoDS_Face aFace = TopoDS::Face(theShape);
3585       BRepAdaptor_Surface SA (aFace, true);
3586       Standard_Real
3587         u1 = SA.FirstUParameter(),
3588         u2 = SA.LastUParameter(),
3589         v1 = SA.FirstVParameter(),
3590         v2 = SA.LastVParameter();
3591       Handle(Geom_Surface) surf = BRep_Tool::Surface(aFace);
3592       myCurProjFace.Init(surf, u1,u2, v1,v2);
3593       myCurFace = aFace;
3594       process();
3595     }
3596     break;
3597   case TopAbs_EDGE:
3598     {
3599       TopoDS_Edge anEdge = TopoDS::Edge(theShape);
3600       Standard_Real u1, u2;
3601       Handle(Geom_Curve) curve = BRep_Tool::Curve(anEdge, u1, u2);
3602       myCurProjEdge.Init(curve, u1, u2);
3603       process();
3604     }
3605     break;
3606   case TopAbs_VERTEX:
3607     {
3608       TopoDS_Vertex aV = TopoDS::Vertex(theShape);
3609       myCurPnt = BRep_Tool::Pnt(aV);
3610       process();
3611     }
3612     break;
3613   default:
3614     break;
3615   }
3616 }
3617
3618 void ElementsOnShape::process()
3619 {
3620   if (myShape.IsNull() || myMesh == 0)
3621     return;
3622
3623   if (myType == SMDSAbs_Node)
3624   {
3625     SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
3626     while (anIter->more())
3627       process(anIter->next());
3628   }
3629   else
3630   {
3631     if (myType == SMDSAbs_Edge || myType == SMDSAbs_All)
3632     {
3633       SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
3634       while (anIter->more())
3635         process(anIter->next());
3636     }
3637
3638     if (myType == SMDSAbs_Face || myType == SMDSAbs_All)
3639     {
3640       SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
3641       while (anIter->more()) {
3642         process(anIter->next());
3643       }
3644     }
3645
3646     if (myType == SMDSAbs_Volume || myType == SMDSAbs_All)
3647     {
3648       SMDS_VolumeIteratorPtr anIter = myMesh->volumesIterator();
3649       while (anIter->more())
3650         process(anIter->next());
3651     }
3652   }
3653 }
3654
3655 void ElementsOnShape::process (const SMDS_MeshElement* theElemPtr)
3656 {
3657   if (myShape.IsNull())
3658     return;
3659
3660   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
3661   bool isSatisfy = myAllNodesFlag;
3662
3663   gp_XYZ centerXYZ (0, 0, 0);
3664
3665   while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
3666   {
3667     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
3668     gp_Pnt aPnt (aNode->X(), aNode->Y(), aNode->Z());
3669     centerXYZ += aPnt.XYZ();
3670
3671     switch (myCurShapeType)
3672     {
3673     case TopAbs_SOLID:
3674       {
3675         myCurSC.Perform(aPnt, myToler);
3676         isSatisfy = (myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON);
3677       }
3678       break;
3679     case TopAbs_FACE:
3680       {
3681         myCurProjFace.Perform(aPnt);
3682         isSatisfy = (myCurProjFace.IsDone() && myCurProjFace.LowerDistance() <= myToler);
3683         if (isSatisfy)
3684         {
3685           // check relatively the face
3686           Quantity_Parameter u, v;
3687           myCurProjFace.LowerDistanceParameters(u, v);
3688           gp_Pnt2d aProjPnt (u, v);
3689           BRepClass_FaceClassifier aClsf (myCurFace, aProjPnt, myToler);
3690           isSatisfy = (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON);
3691         }
3692       }
3693       break;
3694     case TopAbs_EDGE:
3695       {
3696         myCurProjEdge.Perform(aPnt);
3697         isSatisfy = (myCurProjEdge.NbPoints() > 0 && myCurProjEdge.LowerDistance() <= myToler);
3698       }
3699       break;
3700     case TopAbs_VERTEX:
3701       {
3702         isSatisfy = (aPnt.Distance(myCurPnt) <= myToler);
3703       }
3704       break;
3705     default:
3706       {
3707         isSatisfy = false;
3708       }
3709     }
3710   }
3711
3712   if (isSatisfy && myCurShapeType == TopAbs_SOLID) { // Check the center point for volumes MantisBug 0020168
3713     centerXYZ /= theElemPtr->NbNodes();
3714     gp_Pnt aCenterPnt (centerXYZ);
3715     myCurSC.Perform(aCenterPnt, myToler);
3716     if ( !(myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON))
3717       isSatisfy = false;
3718   }
3719
3720   if (isSatisfy)
3721     myIds.Add(theElemPtr->GetID());
3722 }
3723
3724 TSequenceOfXYZ::TSequenceOfXYZ()
3725 {}
3726
3727 TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n)
3728 {}
3729
3730 TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t)
3731 {}
3732
3733 TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray)
3734 {}
3735
3736 template <class InputIterator>
3737 TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd)
3738 {}
3739
3740 TSequenceOfXYZ::~TSequenceOfXYZ()
3741 {}
3742
3743 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
3744 {
3745   myArray = theSequenceOfXYZ.myArray;
3746   return *this;
3747 }
3748
3749 gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
3750 {
3751   return myArray[n-1];
3752 }
3753
3754 const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
3755 {
3756   return myArray[n-1];
3757 }
3758
3759 void TSequenceOfXYZ::clear()
3760 {
3761   myArray.clear();
3762 }
3763
3764 void TSequenceOfXYZ::reserve(size_type n)
3765 {
3766   myArray.reserve(n);
3767 }
3768
3769 void TSequenceOfXYZ::push_back(const gp_XYZ& v)
3770 {
3771   myArray.push_back(v);
3772 }
3773
3774 TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
3775 {
3776   return myArray.size();
3777 }