src/Controls/SMESH_Controls.cxx

   1 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 //  This library is free software; you can redistribute it and/or
   5 //  modify it under the terms of the GNU Lesser General Public
   6 //  License as published by the Free Software Foundation; either
   7 //  version 2.1 of the License.
   8 //
   9 //  This library is distributed in the hope that it will be useful,
  10 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 //  Lesser General Public License for more details.
  13 //
  14 //  You should have received a copy of the GNU Lesser General Public
  15 //  License along with this library; if not, write to the Free Software
  16 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  19
  20 #include "SMESH_Controls.hxx"
  21
  22 #include <set>
  23
  24 #include <BRep_Tool.hxx>
  25 #include <gp_Ax3.hxx>
  26 #include <gp_Cylinder.hxx>
  27 #include <gp_Dir.hxx>
  28 #include <gp_Pnt.hxx>
  29 #include <gp_Pln.hxx>
  30 #include <gp_Vec.hxx>
  31 #include <gp_XYZ.hxx>
  32 #include <Geom_Plane.hxx>
  33 #include <Geom_CylindricalSurface.hxx>
  34 #include <Precision.hxx>
  35 #include <TColgp_Array1OfXYZ.hxx>
  36 #include <TColgp_SequenceOfXYZ.hxx>
  37 #include <TColStd_MapOfInteger.hxx>
  38 #include <TColStd_SequenceOfAsciiString.hxx>
  39 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
  40 #include <TopAbs.hxx>
  41 #include <TopoDS.hxx>
  42 #include <TopoDS_Face.hxx>
  43 #include <TopoDS_Shape.hxx>
  44
  45 #include "SMDS_Mesh.hxx"
  46 #include "SMDS_Iterator.hxx"
  47 #include "SMDS_MeshElement.hxx"
  48 #include "SMDS_MeshNode.hxx"
  49
  50
  51
  52 /*
  53                             AUXILIARY METHODS
  54 */
  55
  56 static inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
  57 {
  58   gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
  59
  60   return v1.Magnitude() < gp::Resolution() ||
  61          v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
  62 }
  63
  64 static inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
  65 {
  66   gp_Vec aVec1( P2 - P1 );
  67   gp_Vec aVec2( P3 - P1 );
  68   return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
  69 }
  70
  71 static inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
  72 {
  73   return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
  74 }
  75
  76 static inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
  77 {
  78   double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
  79   return aDist;
  80 }
  81
  82 static int getNbMultiConnection( SMDS_Mesh* theMesh, const int theId )
  83 {
  84   if ( theMesh == 0 )
  85     return 0;
  86
  87   const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
  88   if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 )
  89     return 0;
  90
  91   TColStd_MapOfInteger aMap;
  92
  93   int aResult = 0;
  94   SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
  95   if ( anIter != 0 )
  96   {
  97     while( anIter->more() )
  98     {
  99       const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
 100       if ( aNode == 0 )
 101         return 0;
 102       SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
 103       while( anElemIter->more() )
 104       {
 105         const SMDS_MeshElement* anElem = anElemIter->next();
 106         if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge )
 107         {
 108           int anId = anElem->GetID();
 109
 110           if ( anIter->more() )              // i.e. first node
 111             aMap.Add( anId );
 112           else if ( aMap.Contains( anId ) )
 113             aResult++;
 114         }
 115       }
 116     }
 117   }
 118
 119   return aResult;
 120 }
 121
 122
 123 using namespace SMESH::Controls;
 124
 125 /*
 126                                 FUNCTORS
 127 */
 128
 129 /*
 130   Class       : NumericalFunctor
 131   Description : Base class for numerical functors
 132 */
 133 NumericalFunctor::NumericalFunctor():
 134   myMesh(NULL)
 135 {
 136   myPrecision = -1;
 137 }
 138
 139 void NumericalFunctor::SetMesh( SMDS_Mesh* theMesh )
 140 {
 141   myMesh = theMesh;
 142 }
 143
 144 bool NumericalFunctor::GetPoints(const int             theId,
 145                                  TColgp_SequenceOfXYZ& theRes ) const
 146 {
 147   theRes.Clear();
 148
 149   if ( myMesh == 0 )
 150     return false;
 151
 152   return GetPoints( myMesh->FindElement( theId ), theRes );
 153 }
 154
 155 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
 156                                  TColgp_SequenceOfXYZ&   theRes )
 157 {
 158   theRes.Clear();
 159
 160   if ( anElem == 0)
 161     return false;
 162
 163   // Get nodes of the element
 164   SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
 165   if ( anIter != 0 )
 166   {
 167     while( anIter->more() )
 168     {
 169       const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
 170       if ( aNode != 0 )
 171         theRes.Append( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
 172     }
 173   }
 174
 175   return true;
 176 }
 177
 178 long  NumericalFunctor::GetPrecision() const
 179 {
 180   return myPrecision;
 181 }
 182
 183 void  NumericalFunctor::SetPrecision( const long thePrecision )
 184 {
 185   myPrecision = thePrecision;
 186 }
 187
 188 double NumericalFunctor::GetValue( long theId )
 189 {
 190   TColgp_SequenceOfXYZ P;
 191   if ( GetPoints( theId, P ))
 192   {
 193     double aVal = GetValue( P );
 194     if ( myPrecision >= 0 )
 195     {
 196       double prec = pow( 10., (double)( myPrecision ) );
 197       aVal = floor( aVal * prec + 0.5 ) / prec;
 198     }
 199
 200     return aVal;
 201   }
 202
 203   return 0.;
 204 }
 205
 206 /*
 207   Class       : MinimumAngle
 208   Description : Functor for calculation of minimum angle
 209 */
 210
 211 double MinimumAngle::GetValue( const TColgp_SequenceOfXYZ& P )
 212 {
 213   double aMin;
 214
 215   if ( P.Length() == 3 )
 216   {
 217     double A0 = getAngle( P( 3 ), P( 1 ), P( 2 ) );
 218     double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
 219     double A2 = getAngle( P( 2 ), P( 3 ), P( 1 ) );
 220
 221     aMin = Min( A0, Min( A1, A2 ) );
 222   }
 223   else if ( P.Length() == 4 )
 224   {
 225     double A0 = getAngle( P( 4 ), P( 1 ), P( 2 ) );
 226     double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
 227     double A2 = getAngle( P( 2 ), P( 3 ), P( 4 ) );
 228     double A3 = getAngle( P( 3 ), P( 4 ), P( 1 ) );
 229
 230     aMin = Min( Min( A0, A1 ), Min( A2, A3 ) );
 231   }
 232   else
 233     return 0.;
 234
 235   return aMin * 180 / PI;
 236 }
 237
 238 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
 239 {
 240   const double aBestAngle = PI / nbNodes;
 241   return ( fabs( aBestAngle - Value ));
 242 }
 243
 244 SMDSAbs_ElementType MinimumAngle::GetType() const
 245 {
 246   return SMDSAbs_Face;
 247 }
 248
 249
 250 /*
 251   Class       : AspectRatio
 252   Description : Functor for calculating aspect ratio
 253 */
 254 double AspectRatio::GetValue( const TColgp_SequenceOfXYZ& P )
 255 {
 256   int nbNodes = P.Length();
 257
 258   if ( nbNodes != 3 && nbNodes != 4 )
 259     return 0;
 260
 261   // Compute lengths of the sides
 262
 263   double aLen[ nbNodes ];
 264   for ( int i = 0; i < nbNodes - 1; i++ )
 265     aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
 266   aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
 267
 268   // Compute aspect ratio
 269
 270   if ( nbNodes == 3 )
 271   {
 272     double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
 273     if ( anArea <= Precision::Confusion() )
 274       return 0.;
 275     double aMaxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
 276     static double aCoef = sqrt( 3. ) / 4;
 277
 278     return aCoef * aMaxLen * aMaxLen / anArea;
 279   }
 280   else
 281   {
 282     double aMinLen = Min( Min( aLen[ 0 ], aLen[ 1 ] ), Min( aLen[ 2 ], aLen[ 3 ] ) );
 283     if ( aMinLen <= Precision::Confusion() )
 284       return 0.;
 285     double aMaxLen = Max( Max( aLen[ 0 ], aLen[ 1 ] ), Max( aLen[ 2 ], aLen[ 3 ] ) );
 286
 287     return aMaxLen / aMinLen;
 288   }
 289 }
 290
 291 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
 292 {
 293   // the aspect ratio is in the range [1.0,infinity]
 294   // 1.0 = good
 295   // infinity = bad
 296   return Value / 1000.;
 297 }
 298
 299 SMDSAbs_ElementType AspectRatio::GetType() const
 300 {
 301   return SMDSAbs_Face;
 302 }
 303
 304
 305 /*
 306   Class       : AspectRatio3D
 307   Description : Functor for calculating aspect ratio
 308 */
 309
 310 static inline double getHalfPerimeter(double theTria[3]){
 311   return (theTria[0] + theTria[1] + theTria[2])/2.0;
 312 }
 313
 314 static inline double getArea(double theHalfPerim, double theTria[3]){
 315   return sqrt(theHalfPerim*
 316               (theHalfPerim-theTria[0])*
 317               (theHalfPerim-theTria[1])*
 318               (theHalfPerim-theTria[2]));
 319 }
 320
 321 static inline double getVolume(double theLen[6]){
 322   double a2 = theLen[0]*theLen[0];
 323   double b2 = theLen[1]*theLen[1];
 324   double c2 = theLen[2]*theLen[2];
 325   double d2 = theLen[3]*theLen[3];
 326   double e2 = theLen[4]*theLen[4];
 327   double f2 = theLen[5]*theLen[5];
 328   double P = 4.0*a2*b2*d2;
 329   double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
 330   double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
 331   return sqrt(P-Q+R)/12.0;
 332 }
 333
 334 static inline double getHeight( const gp_Pnt& P1, const gp_Pnt& P2,
 335                                 const gp_Pnt& P3, const gp_Pnt& P4)
 336 {
 337   gp_Vec aVec1( P2.XYZ() - P1.XYZ() );
 338   gp_Vec aVec2( P3.XYZ() - P1.XYZ() );
 339   gp_Vec aNorm = aVec1 ^ aVec2;
 340   aNorm /= aNorm.Magnitude();
 341   gp_Vec aVec3( P4.XYZ() - P1.XYZ() );
 342   double aDist = aVec1 * aVec2;
 343   return fabs( aDist );
 344 }
 345
 346 static inline double getMaxHeight( const TColgp_SequenceOfXYZ& P )
 347 {
 348   double aHeight = getHeight(P(1),P(2),P(3),P(4));
 349   aHeight = max(aHeight,getHeight(P(1),P(2),P(4),P(3)));
 350   aHeight = max(aHeight,getHeight(P(1),P(3),P(4),P(2)));
 351   aHeight = max(aHeight,getHeight(P(2),P(3),P(4),P(1)));
 352   return aHeight;
 353 }
 354
 355 double AspectRatio3D::GetValue( const TColgp_SequenceOfXYZ& P )
 356 {
 357   double aQuality = 0.0;
 358   int nbNodes = P.Length();
 359   switch(nbNodes){
 360   case 4:{
 361     double aLen[6] = {
 362       getDistance(P(1),P(2)), // a
 363       getDistance(P(2),P(3)), // b
 364       getDistance(P(3),P(1)), // c
 365       getDistance(P(2),P(4)), // d
 366       getDistance(P(3),P(4)), // e
 367       getDistance(P(1),P(4))  // f
 368     };
 369     double aTria[4][3] = {
 370       {aLen[0],aLen[1],aLen[2]}, // abc
 371       {aLen[0],aLen[3],aLen[5]}, // adf
 372       {aLen[1],aLen[3],aLen[4]}, // bde
 373       {aLen[2],aLen[4],aLen[5]}  // cef
 374     };
 375     double aHalfPerim = getHalfPerimeter(aTria[0]);
 376     double anArea = getArea(aHalfPerim,aTria[0]);
 377     aHalfPerim = getHalfPerimeter(aTria[1]);
 378     anArea += getArea(aHalfPerim,aTria[1]);
 379     aHalfPerim = getHalfPerimeter(aTria[2]);
 380     anArea += getArea(aHalfPerim,aTria[2]);
 381     double aVolume = getVolume(aLen);
 382     double aHeight = getMaxHeight(P);
 383     aQuality = 1.0/3.0*aHeight*anArea/aVolume;
 384     break;
 385   }
 386   }
 387   return aQuality;
 388 }
 389
 390 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
 391 {
 392   // the aspect ratio is in the range [1.0,infinity]
 393   // 1.0 = good
 394   // infinity = bad
 395   return Value / 1000.;
 396 }
 397
 398 SMDSAbs_ElementType AspectRatio3D::GetType() const
 399 {
 400   return SMDSAbs_Volume;
 401 }
 402
 403
 404 /*
 405   Class       : Warping
 406   Description : Functor for calculating warping
 407 */
 408 double Warping::GetValue( const TColgp_SequenceOfXYZ& P )
 409 {
 410   if ( P.Length() != 4 )
 411     return 0;
 412
 413   gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4;
 414
 415   double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
 416   double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
 417   double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
 418   double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
 419
 420   return Max( Max( A1, A2 ), Max( A3, A4 ) );
 421 }
 422
 423 double Warping::ComputeA( const gp_XYZ& thePnt1,
 424                           const gp_XYZ& thePnt2,
 425                           const gp_XYZ& thePnt3,
 426                           const gp_XYZ& theG ) const
 427 {
 428   double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
 429   double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
 430   double L = Min( aLen1, aLen2 ) * 0.5;
 431   if ( L < Precision::Confusion())
 432     return 0.;
 433
 434   gp_XYZ GI = ( thePnt2 - thePnt1 ) / 2. - theG;
 435   gp_XYZ GJ = ( thePnt3 - thePnt2 ) / 2. - theG;
 436   gp_XYZ N  = GI.Crossed( GJ );
 437
 438   if ( N.Modulus() < gp::Resolution() )
 439     return PI / 2;
 440
 441   N.Normalize();
 442
 443   double H = ( thePnt2 - theG ).Dot( N );
 444   return asin( fabs( H / L ) ) * 180 / PI;
 445 }
 446
 447 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
 448 {
 449   // the warp is in the range [0.0,PI/2]
 450   // 0.0 = good (no warp)
 451   // PI/2 = bad  (face pliee)
 452   return Value;
 453 }
 454
 455 SMDSAbs_ElementType Warping::GetType() const
 456 {
 457   return SMDSAbs_Face;
 458 }
 459
 460
 461 /*
 462   Class       : Taper
 463   Description : Functor for calculating taper
 464 */
 465 double Taper::GetValue( const TColgp_SequenceOfXYZ& P )
 466 {
 467   if ( P.Length() != 4 )
 468     return 0;
 469
 470   // Compute taper
 471   double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2;
 472   double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2;
 473   double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2;
 474   double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2;
 475
 476   double JA = 0.25 * ( J1 + J2 + J3 + J4 );
 477   if ( JA <= Precision::Confusion() )
 478     return 0.;
 479
 480   double T1 = fabs( ( J1 - JA ) / JA );
 481   double T2 = fabs( ( J2 - JA ) / JA );
 482   double T3 = fabs( ( J3 - JA ) / JA );
 483   double T4 = fabs( ( J4 - JA ) / JA );
 484
 485   return Max( Max( T1, T2 ), Max( T3, T4 ) );
 486 }
 487
 488 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
 489 {
 490   // the taper is in the range [0.0,1.0]
 491   // 0.0  = good (no taper)
 492   // 1.0 = bad  (les cotes opposes sont allignes)
 493   return Value;
 494 }
 495
 496 SMDSAbs_ElementType Taper::GetType() const
 497 {
 498   return SMDSAbs_Face;
 499 }
 500
 501
 502 /*
 503   Class       : Skew
 504   Description : Functor for calculating skew in degrees
 505 */
 506 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
 507 {
 508   gp_XYZ p12 = ( p2 + p1 ) / 2;
 509   gp_XYZ p23 = ( p3 + p2 ) / 2;
 510   gp_XYZ p31 = ( p3 + p1 ) / 2;
 511
 512   gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
 513
 514   return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
 515 }
 516
 517 double Skew::GetValue( const TColgp_SequenceOfXYZ& P )
 518 {
 519   if ( P.Length() != 3 && P.Length() != 4 )
 520     return 0;
 521
 522   // Compute skew
 523   static double PI2 = PI / 2;
 524   if ( P.Length() == 3 )
 525   {
 526     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
 527     double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
 528     double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
 529
 530     return Max( A0, Max( A1, A2 ) ) * 180 / PI;
 531   }
 532   else
 533   {
 534     gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2;
 535     gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2;
 536     gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2;
 537     gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2;
 538
 539     gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
 540     double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
 541       ? 0 : fabs( PI2 - v1.Angle( v2 ) );
 542
 543     return A * 180 / PI;
 544   }
 545 }
 546
 547 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
 548 {
 549   // the skew is in the range [0.0,PI/2].
 550   // 0.0 = good
 551   // PI/2 = bad
 552   return Value;
 553 }
 554
 555 SMDSAbs_ElementType Skew::GetType() const
 556 {
 557   return SMDSAbs_Face;
 558 }
 559
 560
 561 /*
 562   Class       : Area
 563   Description : Functor for calculating area
 564 */
 565 double Area::GetValue( const TColgp_SequenceOfXYZ& P )
 566 {
 567   if ( P.Length() == 3 )
 568     return getArea( P( 1 ), P( 2 ), P( 3 ) );
 569   else if ( P.Length() == 4 )
 570     return getArea( P( 1 ), P( 2 ), P( 3 ) ) + getArea( P( 1 ), P( 3 ), P( 4 ) );
 571   else
 572     return 0;
 573 }
 574
 575 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
 576 {
 577   return Value;
 578 }
 579
 580 SMDSAbs_ElementType Area::GetType() const
 581 {
 582   return SMDSAbs_Face;
 583 }
 584
 585
 586 /*
 587   Class       : Length
 588   Description : Functor for calculating length off edge
 589 */
 590 double Length::GetValue( const TColgp_SequenceOfXYZ& P )
 591 {
 592   return ( P.Length() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0 );
 593 }
 594
 595 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
 596 {
 597   return Value;
 598 }
 599
 600 SMDSAbs_ElementType Length::GetType() const
 601 {
 602   return SMDSAbs_Edge;
 603 }
 604
 605
 606 /*
 607   Class       : MultiConnection
 608   Description : Functor for calculating number of faces conneted to the edge
 609 */
 610 double MultiConnection::GetValue( const TColgp_SequenceOfXYZ& P )
 611 {
 612   return 0;
 613 }
 614 double MultiConnection::GetValue( long theId )
 615 {
 616   return getNbMultiConnection( myMesh, theId );
 617 }
 618
 619 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
 620 {
 621   return Value;
 622 }
 623
 624 SMDSAbs_ElementType MultiConnection::GetType() const
 625 {
 626   return SMDSAbs_Edge;
 627 }
 628
 629
 630 /*
 631                             PREDICATES
 632 */
 633
 634 /*
 635   Class       : FreeBorders
 636   Description : Predicate for free borders
 637 */
 638
 639 FreeBorders::FreeBorders()
 640 {
 641   myMesh = 0;
 642 }
 643
 644 void FreeBorders::SetMesh( SMDS_Mesh* theMesh )
 645 {
 646   myMesh = theMesh;
 647 }
 648
 649 bool FreeBorders::IsSatisfy( long theId )
 650 {
 651   return getNbMultiConnection( myMesh, theId ) == 1;
 652 }
 653
 654 SMDSAbs_ElementType FreeBorders::GetType() const
 655 {
 656   return SMDSAbs_Edge;
 657 }
 658
 659
 660 /*
 661   Class       : FreeEdges
 662   Description : Predicate for free Edges
 663 */
 664 FreeEdges::FreeEdges()
 665 {
 666   myMesh = 0;
 667 }
 668
 669 void FreeEdges::SetMesh( SMDS_Mesh* theMesh )
 670 {
 671   myMesh = theMesh;
 672 }
 673
 674 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId  )
 675 {
 676   TColStd_MapOfInteger aMap;
 677   for ( int i = 0; i < 2; i++ )
 678   {
 679     SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator();
 680     while( anElemIter->more() )
 681     {
 682       const SMDS_MeshElement* anElem = anElemIter->next();
 683       if ( anElem != 0 && anElem->GetType() == SMDSAbs_Face )
 684       {
 685         int anId = anElem->GetID();
 686
 687         if ( i == 0 )
 688           aMap.Add( anId );
 689         else if ( aMap.Contains( anId ) && anId != theFaceId )
 690           return false;
 691       }
 692     }
 693   }
 694   return true;
 695 }
 696
 697 bool FreeEdges::IsSatisfy( long theId )
 698 {
 699   if ( myMesh == 0 )
 700     return false;
 701
 702   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
 703   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
 704     return false;
 705
 706   int nbNodes = aFace->NbNodes();
 707   const SMDS_MeshNode* aNodes[ nbNodes ];
 708   int i = 0;
 709   SMDS_ElemIteratorPtr anIter = aFace->nodesIterator();
 710   if ( anIter != 0 )
 711   {
 712     while( anIter->more() )
 713     {
 714       const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
 715       if ( aNode == 0 )
 716         return false;
 717       aNodes[ i++ ] = aNode;
 718     }
 719   }
 720
 721   for ( int i = 0; i < nbNodes - 1; i++ )
 722     if ( IsFreeEdge( &aNodes[ i ], theId ) )
 723       return true;
 724
 725   aNodes[ 1 ] = aNodes[ nbNodes - 1 ];
 726
 727   return IsFreeEdge( &aNodes[ 0 ], theId );
 728
 729 }
 730
 731 SMDSAbs_ElementType FreeEdges::GetType() const
 732 {
 733   return SMDSAbs_Face;
 734 }
 735
 736 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
 737   myElemId(theElemId)
 738 {
 739   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
 740   if(thePntId1 > thePntId2){
 741     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
 742   }
 743 }
 744
 745 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
 746   if(myPntId[0] < x.myPntId[0]) return true;
 747   if(myPntId[0] == x.myPntId[0])
 748     if(myPntId[1] < x.myPntId[1]) return true;
 749   return false;
 750 }
 751
 752 inline void UpdateBorders(const FreeEdges::Border& theBorder,
 753                           FreeEdges::TBorders& theRegistry,
 754                           FreeEdges::TBorders& theContainer)
 755 {
 756   if(theRegistry.find(theBorder) == theRegistry.end()){
 757     theRegistry.insert(theBorder);
 758     theContainer.insert(theBorder);
 759   }else{
 760     theContainer.erase(theBorder);
 761   }
 762 }
 763
 764 void FreeEdges::GetBoreders(TBorders& theBorders)
 765 {
 766   TBorders aRegistry;
 767   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
 768   for(; anIter->more(); ){
 769     const SMDS_MeshFace* anElem = anIter->next();
 770     long anElemId = anElem->GetID();
 771     SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
 772     long aNodeId[2];
 773     const SMDS_MeshElement* aNode;
 774     if(aNodesIter->more()){
 775       aNode = aNodesIter->next();
 776       aNodeId[0] = aNodeId[1] = aNode->GetID();
 777     }
 778     for(; aNodesIter->more(); ){
 779       aNode = aNodesIter->next();
 780       long anId = aNode->GetID();
 781       Border aBorder(anElemId,aNodeId[1],anId);
 782       aNodeId[1] = anId;
 783       //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
 784       UpdateBorders(aBorder,aRegistry,theBorders);
 785     }
 786     Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
 787     //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
 788     UpdateBorders(aBorder,aRegistry,theBorders);
 789   }
 790   //std::cout<<"theBorders.size() = "<<theBorders.size()<<endl;
 791 }
 792
 793 /*
 794   Class       : RangeOfIds
 795   Description : Predicate for Range of Ids.
 796                 Range may be specified with two ways.
 797                 1. Using AddToRange method
 798                 2. With SetRangeStr method. Parameter of this method is a string
 799                    like as "1,2,3,50-60,63,67,70-"
 800 */
 801
 802 //=======================================================================
 803 // name    : RangeOfIds
 804 // Purpose : Constructor
 805 //=======================================================================
 806 RangeOfIds::RangeOfIds()
 807 {
 808   myMesh = 0;
 809   myType = SMDSAbs_All;
 810 }
 811
 812 //=======================================================================
 813 // name    : SetMesh
 814 // Purpose : Set mesh
 815 //=======================================================================
 816 void RangeOfIds::SetMesh( SMDS_Mesh* theMesh )
 817 {
 818   myMesh = theMesh;
 819 }
 820
 821 //=======================================================================
 822 // name    : AddToRange
 823 // Purpose : Add ID to the range
 824 //=======================================================================
 825 bool RangeOfIds::AddToRange( long theEntityId )
 826 {
 827   myIds.Add( theEntityId );
 828   return true;
 829 }
 830
 831 //=======================================================================
 832 // name    : GetRangeStr
 833 // Purpose : Get range as a string.
 834 //           Example: "1,2,3,50-60,63,67,70-"
 835 //=======================================================================
 836 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
 837 {
 838   theResStr.Clear();
 839
 840   TColStd_SequenceOfInteger     anIntSeq;
 841   TColStd_SequenceOfAsciiString aStrSeq;
 842
 843   TColStd_MapIteratorOfMapOfInteger anIter( myIds );
 844   for ( ; anIter.More(); anIter.Next() )
 845   {
 846     int anId = anIter.Key();
 847     TCollection_AsciiString aStr( anId );
 848     anIntSeq.Append( anId );
 849     aStrSeq.Append( aStr );
 850   }
 851
 852   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
 853   {
 854     int aMinId = myMin( i );
 855     int aMaxId = myMax( i );
 856
 857     TCollection_AsciiString aStr;
 858     if ( aMinId != IntegerFirst() )
 859       aStr += aMinId;
 860
 861     aStr += "-";
 862
 863     if ( aMaxId != IntegerLast() )
 864       aStr += aMaxId;
 865
 866     // find position of the string in result sequence and insert string in it
 867     if ( anIntSeq.Length() == 0 )
 868     {
 869       anIntSeq.Append( aMinId );
 870       aStrSeq.Append( aStr );
 871     }
 872     else
 873     {
 874       if ( aMinId < anIntSeq.First() )
 875       {
 876         anIntSeq.Prepend( aMinId );
 877         aStrSeq.Prepend( aStr );
 878       }
 879       else if ( aMinId > anIntSeq.Last() )
 880       {
 881         anIntSeq.Append( aMinId );
 882         aStrSeq.Append( aStr );
 883       }
 884       else
 885         for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
 886           if ( aMinId < anIntSeq( j ) )
 887           {
 888             anIntSeq.InsertBefore( j, aMinId );
 889             aStrSeq.InsertBefore( j, aStr );
 890             break;
 891           }
 892     }
 893   }
 894
 895   if ( aStrSeq.Length() == 0 )
 896     return;
 897
 898   theResStr = aStrSeq( 1 );
 899   for ( int j = 2, k = aStrSeq.Length(); j <= k; j++  )
 900   {
 901     theResStr += ",";
 902     theResStr += aStrSeq( j );
 903   }
 904 }
 905
 906 //=======================================================================
 907 // name    : SetRangeStr
 908 // Purpose : Define range with string
 909 //           Example of entry string: "1,2,3,50-60,63,67,70-"
 910 //=======================================================================
 911 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
 912 {
 913   myMin.Clear();
 914   myMax.Clear();
 915   myIds.Clear();
 916
 917   TCollection_AsciiString aStr = theStr;
 918   aStr.RemoveAll( ' ' );
 919   aStr.RemoveAll( '\t' );
 920
 921   for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
 922     aStr.Remove( aPos, 2 );
 923
 924   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
 925   int i = 1;
 926   while ( tmpStr != "" )
 927   {
 928     tmpStr = aStr.Token( ",", i++ );
 929     int aPos = tmpStr.Search( '-' );
 930
 931     if ( aPos == -1 )
 932     {
 933       if ( tmpStr.IsIntegerValue() )
 934         myIds.Add( tmpStr.IntegerValue() );
 935       else
 936         return false;
 937     }
 938     else
 939     {
 940       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
 941       TCollection_AsciiString aMinStr = tmpStr;
 942
 943       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
 944       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
 945
 946       if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() ||
 947            !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() )
 948         return false;
 949
 950       myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
 951       myMax.Append( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
 952     }
 953   }
 954
 955   return true;
 956 }
 957
 958 //=======================================================================
 959 // name    : GetType
 960 // Purpose : Get type of supported entities
 961 //=======================================================================
 962 SMDSAbs_ElementType RangeOfIds::GetType() const
 963 {
 964   return myType;
 965 }
 966
 967 //=======================================================================
 968 // name    : SetType
 969 // Purpose : Set type of supported entities
 970 //=======================================================================
 971 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
 972 {
 973   myType = theType;
 974 }
 975
 976 //=======================================================================
 977 // name    : IsSatisfy
 978 // Purpose : Verify whether entity satisfies to this rpedicate
 979 //=======================================================================
 980 bool RangeOfIds::IsSatisfy( long theId )
 981 {
 982   if ( !myMesh )
 983     return false;
 984
 985   if ( myType == SMDSAbs_Node )
 986   {
 987     if ( myMesh->FindNode( theId ) == 0 )
 988       return false;
 989   }
 990   else
 991   {
 992     const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
 993     if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All )
 994       return false;
 995   }
 996
 997   if ( myIds.Contains( theId ) )
 998     return true;
 999
1000   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
1001     if ( theId >= myMin( i ) && theId <= myMax( i ) )
1002       return true;
1003
1004   return false;
1005 }
1006
1007 /*
1008   Class       : Comparator
1009   Description : Base class for comparators
1010 */
1011 Comparator::Comparator():
1012   myMargin(0)
1013 {}
1014
1015 Comparator::~Comparator()
1016 {}
1017
1018 void Comparator::SetMesh( SMDS_Mesh* theMesh )
1019 {
1020   if ( myFunctor )
1021     myFunctor->SetMesh( theMesh );
1022 }
1023
1024 void Comparator::SetMargin( double theValue )
1025 {
1026   myMargin = theValue;
1027 }
1028
1029 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
1030 {
1031   myFunctor = theFunct;
1032 }
1033
1034 SMDSAbs_ElementType Comparator::GetType() const
1035 {
1036   return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
1037 }
1038
1039 double Comparator::GetMargin()
1040 {
1041   return myMargin;
1042 }
1043
1044
1045 /*
1046   Class       : LessThan
1047   Description : Comparator "<"
1048 */
1049 bool LessThan::IsSatisfy( long theId )
1050 {
1051   return myFunctor && myFunctor->GetValue( theId ) < myMargin;
1052 }
1053
1054
1055 /*
1056   Class       : MoreThan
1057   Description : Comparator ">"
1058 */
1059 bool MoreThan::IsSatisfy( long theId )
1060 {
1061   return myFunctor && myFunctor->GetValue( theId ) > myMargin;
1062 }
1063
1064
1065 /*
1066   Class       : EqualTo
1067   Description : Comparator "="
1068 */
1069 EqualTo::EqualTo():
1070   myToler(Precision::Confusion())
1071 {}
1072
1073 bool EqualTo::IsSatisfy( long theId )
1074 {
1075   return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
1076 }
1077
1078 void EqualTo::SetTolerance( double theToler )
1079 {
1080   myToler = theToler;
1081 }
1082
1083 double EqualTo::GetTolerance()
1084 {
1085   return myToler;
1086 }
1087
1088 /*
1089   Class       : LogicalNOT
1090   Description : Logical NOT predicate
1091 */
1092 LogicalNOT::LogicalNOT()
1093 {}
1094
1095 LogicalNOT::~LogicalNOT()
1096 {}
1097
1098 bool LogicalNOT::IsSatisfy( long theId )
1099 {
1100   return myPredicate && !myPredicate->IsSatisfy( theId );
1101 }
1102
1103 void LogicalNOT::SetMesh( SMDS_Mesh* theMesh )
1104 {
1105   if ( myPredicate )
1106     myPredicate->SetMesh( theMesh );
1107 }
1108
1109 void LogicalNOT::SetPredicate( PredicatePtr thePred )
1110 {
1111   myPredicate = thePred;
1112 }
1113
1114 SMDSAbs_ElementType LogicalNOT::GetType() const
1115 {
1116   return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
1117 }
1118
1119
1120 /*
1121   Class       : LogicalBinary
1122   Description : Base class for binary logical predicate
1123 */
1124 LogicalBinary::LogicalBinary()
1125 {}
1126
1127 LogicalBinary::~LogicalBinary()
1128 {}
1129
1130 void LogicalBinary::SetMesh( SMDS_Mesh* theMesh )
1131 {
1132   if ( myPredicate1 )
1133     myPredicate1->SetMesh( theMesh );
1134
1135   if ( myPredicate2 )
1136     myPredicate2->SetMesh( theMesh );
1137 }
1138
1139 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
1140 {
1141   myPredicate1 = thePredicate;
1142 }
1143
1144 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
1145 {
1146   myPredicate2 = thePredicate;
1147 }
1148
1149 SMDSAbs_ElementType LogicalBinary::GetType() const
1150 {
1151   if ( !myPredicate1 || !myPredicate2 )
1152     return SMDSAbs_All;
1153
1154   SMDSAbs_ElementType aType1 = myPredicate1->GetType();
1155   SMDSAbs_ElementType aType2 = myPredicate2->GetType();
1156
1157   return aType1 == aType2 ? aType1 : SMDSAbs_All;
1158 }
1159
1160
1161 /*
1162   Class       : LogicalAND
1163   Description : Logical AND
1164 */
1165 bool LogicalAND::IsSatisfy( long theId )
1166 {
1167   return
1168     myPredicate1 &&
1169     myPredicate2 &&
1170     myPredicate1->IsSatisfy( theId ) &&
1171     myPredicate2->IsSatisfy( theId );
1172 }
1173
1174
1175 /*
1176   Class       : LogicalOR
1177   Description : Logical OR
1178 */
1179 bool LogicalOR::IsSatisfy( long theId )
1180 {
1181   return
1182     myPredicate1 &&
1183     myPredicate2 &&
1184     myPredicate1->IsSatisfy( theId ) ||
1185     myPredicate2->IsSatisfy( theId );
1186 }
1187
1188
1189 /*
1190                               FILTER
1191 */
1192
1193 Filter::Filter()
1194 {}
1195
1196 Filter::~Filter()
1197 {}
1198
1199 void Filter::SetPredicate( PredicatePtr thePredicate )
1200 {
1201   myPredicate = thePredicate;
1202 }
1203
1204
1205 template<class TElement, class TIterator, class TPredicate>
1206 void FillSequence(const TIterator& theIterator,
1207                   TPredicate& thePredicate,
1208                   Filter::TIdSequence& theSequence)
1209 {
1210   if ( theIterator ) {
1211     while( theIterator->more() ) {
1212       TElement anElem = theIterator->next();
1213       long anId = anElem->GetID();
1214       if ( thePredicate->IsSatisfy( anId ) )
1215         theSequence.push_back( anId );
1216     }
1217   }
1218 }
1219
1220 Filter::TIdSequence
1221 Filter::GetElementsId( SMDS_Mesh* theMesh )
1222 {
1223   TIdSequence aSequence;
1224   if ( !theMesh || !myPredicate ) return aSequence;
1225
1226   myPredicate->SetMesh( theMesh );
1227
1228   SMDSAbs_ElementType aType = myPredicate->GetType();
1229   switch(aType){
1230   case SMDSAbs_Node:{
1231     FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),myPredicate,aSequence);
1232     break;
1233   }
1234   case SMDSAbs_Edge:{
1235     FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1236     break;
1237   }
1238   case SMDSAbs_Face:{
1239     FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1240     break;
1241   }
1242   case SMDSAbs_Volume:{
1243     FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1244     break;
1245   }
1246   case SMDSAbs_All:{
1247     FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1248     FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1249     FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1250     break;
1251   }
1252   }
1253   return aSequence;
1254 }
1255
1256 /*
1257                               ManifoldPart
1258 */
1259
1260 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
1261
1262 /*
1263    Internal class Link
1264 */
1265
1266 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
1267                           SMDS_MeshNode* theNode2 )
1268 {
1269   myNode1 = theNode1;
1270   myNode2 = theNode2;
1271 }
1272
1273 ManifoldPart::Link::~Link()
1274 {
1275   myNode1 = 0;
1276   myNode2 = 0;
1277 }
1278
1279 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
1280 {
1281   if ( myNode1 == theLink.myNode1 &&
1282        myNode2 == theLink.myNode2 )
1283     return true;
1284   else if ( myNode1 == theLink.myNode2 &&
1285             myNode2 == theLink.myNode1 )
1286     return true;
1287   else
1288     return false;
1289 }
1290
1291 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
1292 {
1293   if(myNode1 < x.myNode1) return true;
1294   if(myNode1 == x.myNode1)
1295     if(myNode2 < x.myNode2) return true;
1296   return false;
1297 }
1298
1299 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
1300                             const ManifoldPart::Link& theLink2 )
1301 {
1302   return theLink1.IsEqual( theLink2 );
1303 }
1304
1305 ManifoldPart::ManifoldPart()
1306 {
1307   myMesh = 0;
1308   myAngToler = Precision::Angular();
1309   myIsOnlyManifold = true;
1310 }
1311
1312 ManifoldPart::~ManifoldPart()
1313 {
1314   myMesh = 0;
1315 }
1316
1317 void ManifoldPart::SetMesh( SMDS_Mesh* theMesh )
1318 {
1319   myMesh = theMesh;
1320   process();
1321 }
1322
1323 SMDSAbs_ElementType ManifoldPart::GetType() const
1324 { return SMDSAbs_Face; }
1325
1326 bool ManifoldPart::IsSatisfy( long theElementId )
1327 {
1328   return myMapIds.Contains( theElementId );
1329 }
1330
1331 void ManifoldPart::SetAngleTolerance( const double theAngToler )
1332 { myAngToler = theAngToler; }
1333
1334 double ManifoldPart::GetAngleTolerance() const
1335 { return myAngToler; }
1336
1337 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
1338 { myIsOnlyManifold = theIsOnly; }
1339
1340 void ManifoldPart::SetStartElem( const long  theStartId )
1341 { myStartElemId = theStartId; }
1342
1343 bool ManifoldPart::process()
1344 {
1345   myMapIds.Clear();
1346   myMapBadGeomIds.Clear();
1347
1348   myAllFacePtr.clear();
1349   myAllFacePtrIntDMap.clear();
1350   if ( !myMesh )
1351     return false;
1352
1353   // collect all faces into own map
1354   SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
1355   for (; anFaceItr->more(); )
1356   {
1357     SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
1358     myAllFacePtr.push_back( aFacePtr );
1359     myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
1360   }
1361
1362   SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
1363   if ( !aStartFace )
1364     return false;
1365
1366   // the map of non manifold links and bad geometry
1367   TMapOfLink aMapOfNonManifold;
1368   TColStd_MapOfInteger aMapOfTreated;
1369
1370   // begin cycle on faces from start index and run on vector till the end
1371   //  and from begin to start index to cover whole vector
1372   const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
1373   bool isStartTreat = false;
1374   for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
1375   {
1376     if ( fi == aStartIndx )
1377       isStartTreat = true;
1378     // as result next time when fi will be equal to aStartIndx
1379
1380     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
1381     if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
1382       continue;
1383
1384     aMapOfTreated.Add( aFacePtr->GetID() );
1385     TColStd_MapOfInteger aResFaces;
1386     if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
1387                          aMapOfNonManifold, aResFaces ) )
1388       continue;
1389     TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
1390     for ( ; anItr.More(); anItr.Next() )
1391     {
1392       int aFaceId = anItr.Key();
1393       aMapOfTreated.Add( aFaceId );
1394       myMapIds.Add( aFaceId );
1395     }
1396
1397     if ( fi == ( myAllFacePtr.size() - 1 ) )
1398       fi = 0;
1399   } // end run on vector of faces
1400   return !myMapIds.IsEmpty();
1401 }
1402
1403 static void getLinks( const SMDS_MeshFace* theFace,
1404                       ManifoldPart::TVectorOfLink& theLinks )
1405 {
1406   int aNbNode = theFace->NbNodes();
1407   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1408   int i = 1;
1409   SMDS_MeshNode* aNode = 0;
1410   for ( ; aNodeItr->more() && i <= aNbNode; )
1411   {
1412
1413     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
1414     if ( i == 1 )
1415       aNode = aN1;
1416     i++;
1417     SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
1418     i++;
1419     ManifoldPart::Link aLink( aN1, aN2 );
1420     theLinks.push_back( aLink );
1421   }
1422 }
1423
1424 static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
1425 {
1426   gp_XYZ n;
1427   int aNbNode = theFace->NbNodes();
1428   TColgp_Array1OfXYZ anArrOfXYZ(1,4);
1429   gp_XYZ p1, p2, p3, p4;
1430   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1431   int i = 1;
1432   for ( ; aNodeItr->more() && i <= 4; i++ )
1433   {
1434     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
1435     anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
1436   }
1437
1438   gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1);
1439   gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1);
1440   n  = q1 ^ q2;
1441   if ( aNbNode > 3 )
1442   {
1443     gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
1444     n += q2 ^ q3;
1445   }
1446   double len = n.Modulus();
1447   if ( len > 0 )
1448     n /= len;
1449
1450   return n;
1451 }
1452
1453 bool ManifoldPart::findConnected
1454                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
1455                   SMDS_MeshFace*                           theStartFace,
1456                   ManifoldPart::TMapOfLink&                theNonManifold,
1457                   TColStd_MapOfInteger&                    theResFaces )
1458 {
1459   theResFaces.Clear();
1460   if ( !theAllFacePtrInt.size() )
1461     return false;
1462
1463   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
1464   {
1465     myMapBadGeomIds.Add( theStartFace->GetID() );
1466     return false;
1467   }
1468
1469   ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
1470   ManifoldPart::TVectorOfLink aSeqOfBoundary;
1471   theResFaces.Add( theStartFace->GetID() );
1472   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
1473
1474   expandBoundary( aMapOfBoundary, aSeqOfBoundary,
1475                  aDMapLinkFace, theNonManifold, theStartFace );
1476
1477   bool isDone = false;
1478   while ( !isDone && aMapOfBoundary.size() != 0 )
1479   {
1480     bool isToReset = false;
1481     ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
1482     for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
1483     {
1484       ManifoldPart::Link aLink = *pLink;
1485       if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
1486         continue;
1487       // each link could be treated only once
1488       aMapToSkip.insert( aLink );
1489
1490       ManifoldPart::TVectorOfFacePtr aFaces;
1491       // find next
1492       if ( myIsOnlyManifold &&
1493            (theNonManifold.find( aLink ) != theNonManifold.end()) )
1494         continue;
1495       else
1496       {
1497         getFacesByLink( aLink, aFaces );
1498         // filter the element to keep only indicated elements
1499         ManifoldPart::TVectorOfFacePtr aFiltered;
1500         ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
1501         for ( ; pFace != aFaces.end(); ++pFace )
1502         {
1503           SMDS_MeshFace* aFace = *pFace;
1504           if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
1505             aFiltered.push_back( aFace );
1506         }
1507         aFaces = aFiltered;
1508         if ( aFaces.size() < 2 )  // no neihgbour faces
1509           continue;
1510         else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
1511         {
1512           theNonManifold.insert( aLink );
1513           continue;
1514         }
1515       }
1516
1517       // compare normal with normals of neighbor element
1518       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
1519       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
1520       for ( ; pFace != aFaces.end(); ++pFace )
1521       {
1522         SMDS_MeshFace* aNextFace = *pFace;
1523         if ( aPrevFace == aNextFace )
1524           continue;
1525         int anNextFaceID = aNextFace->GetID();
1526         if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
1527          // should not be with non manifold restriction. probably bad topology
1528           continue;
1529         // check if face was treated and skipped
1530         if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
1531              !isInPlane( aPrevFace, aNextFace ) )
1532           continue;
1533         // add new element to connected and extend the boundaries.
1534         theResFaces.Add( anNextFaceID );
1535         expandBoundary( aMapOfBoundary, aSeqOfBoundary,
1536                         aDMapLinkFace, theNonManifold, aNextFace );
1537         isToReset = true;
1538       }
1539     }
1540     isDone = !isToReset;
1541   }
1542
1543   return !theResFaces.IsEmpty();
1544 }
1545
1546 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
1547                               const SMDS_MeshFace* theFace2 )
1548 {
1549   gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
1550   gp_XYZ aNorm2XYZ = getNormale( theFace2 );
1551   if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
1552   {
1553     myMapBadGeomIds.Add( theFace2->GetID() );
1554     return false;
1555   }
1556   if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
1557     return true;
1558
1559   return false;
1560 }
1561
1562 void ManifoldPart::expandBoundary
1563                    ( ManifoldPart::TMapOfLink&            theMapOfBoundary,
1564                      ManifoldPart::TVectorOfLink&         theSeqOfBoundary,
1565                      ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
1566                      ManifoldPart::TMapOfLink&            theNonManifold,
1567                      SMDS_MeshFace*                       theNextFace ) const
1568 {
1569   ManifoldPart::TVectorOfLink aLinks;
1570   getLinks( theNextFace, aLinks );
1571   int aNbLink = aLinks.size();
1572   for ( int i = 0; i < aNbLink; i++ )
1573   {
1574     ManifoldPart::Link aLink = aLinks[ i ];
1575     if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
1576       continue;
1577     if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
1578     {
1579       if ( myIsOnlyManifold )
1580       {
1581         // remove from boundary
1582         theMapOfBoundary.erase( aLink );
1583         ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
1584         for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
1585         {
1586           ManifoldPart::Link aBoundLink = *pLink;
1587           if ( aBoundLink.IsEqual( aLink ) )
1588           {
1589             theSeqOfBoundary.erase( pLink );
1590             break;
1591           }
1592         }
1593       }
1594     }
1595     else
1596     {
1597       theMapOfBoundary.insert( aLink );
1598       theSeqOfBoundary.push_back( aLink );
1599       theDMapLinkFacePtr[ aLink ] = theNextFace;
1600     }
1601   }
1602 }
1603
1604 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
1605                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
1606 {
1607   SMDS_Mesh::SetOfFaces aSetOfFaces;
1608   // take all faces that shared first node
1609   SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
1610   for ( ; anItr->more(); )
1611   {
1612     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
1613     if ( !aFace )
1614       continue;
1615     aSetOfFaces.insert( aFace );
1616   }
1617   // take all faces that shared second node
1618   anItr = theLink.myNode2->facesIterator();
1619   // find the common part of two sets
1620   for ( ; anItr->more(); )
1621   {
1622     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
1623     if ( aSetOfFaces.find( aFace ) != aSetOfFaces.end() )
1624       theFaces.push_back( aFace );
1625   }
1626 }
1627
1628
1629 /*
1630    ElementsOnSurface
1631 */
1632
1633 ElementsOnSurface::ElementsOnSurface()
1634 {
1635   myMesh = 0;
1636   myIds.Clear();
1637   myType = SMDSAbs_All;
1638   mySurf.Nullify();
1639   myToler = Precision::Confusion();
1640 }
1641
1642 ElementsOnSurface::~ElementsOnSurface()
1643 {
1644   myMesh = 0;
1645 }
1646
1647 void ElementsOnSurface::SetMesh( SMDS_Mesh* theMesh )
1648 {
1649   if ( myMesh == theMesh )
1650     return;
1651   myMesh = theMesh;
1652   myIds.Clear();
1653   process();
1654 }
1655
1656 bool ElementsOnSurface::IsSatisfy( long theElementId )
1657 {
1658   return myIds.Contains( theElementId );
1659 }
1660
1661 SMDSAbs_ElementType ElementsOnSurface::GetType() const
1662 { return myType; }
1663
1664 void ElementsOnSurface::SetTolerance( const double theToler )
1665 { myToler = theToler; }
1666
1667 double ElementsOnSurface::GetTolerance() const
1668 {
1669   return myToler;
1670 }
1671
1672 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
1673                                     const SMDSAbs_ElementType theType )
1674 {
1675   myType = theType;
1676   mySurf.Nullify();
1677   if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
1678   {
1679     mySurf.Nullify();
1680     return;
1681   }
1682   TopoDS_Face aFace = TopoDS::Face( theShape );
1683   mySurf = BRep_Tool::Surface( aFace );
1684 }
1685
1686 void ElementsOnSurface::process()
1687 {
1688   myIds.Clear();
1689   if ( mySurf.IsNull() )
1690     return;
1691
1692   if ( myMesh == 0 )
1693     return;
1694
1695   if ( myType == SMDSAbs_Face || myType == SMDSAbs_All )
1696   {
1697     SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1698     for(; anIter->more(); )
1699       process( anIter->next() );
1700   }
1701
1702   if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All )
1703   {
1704     SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
1705     for(; anIter->more(); )
1706       process( anIter->next() );
1707   }
1708
1709   if ( myType == SMDSAbs_Node )
1710   {
1711     SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
1712     for(; anIter->more(); )
1713       process( anIter->next() );
1714   }
1715 }
1716
1717 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
1718 {
1719   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
1720   bool isSatisfy = true;
1721   for ( ; aNodeItr->more(); )
1722   {
1723     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
1724     if ( !isOnSurface( aNode ) )
1725     {
1726       isSatisfy = false;
1727       break;
1728     }
1729   }
1730   if ( isSatisfy )
1731     myIds.Add( theElemPtr->GetID() );
1732 }
1733
1734 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) const
1735 {
1736   if ( mySurf.IsNull() )
1737     return false;
1738
1739   gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
1740   double aToler2 = myToler * myToler;
1741   if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
1742   {
1743     gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
1744     if ( aPln.SquareDistance( aPnt ) > aToler2 )
1745       return false;
1746   }
1747   else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
1748   {
1749     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
1750     double aRad = aCyl.Radius();
1751     gp_Ax3 anAxis = aCyl.Position();
1752     gp_XYZ aLoc = aCyl.Location().XYZ();
1753     double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
1754     double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
1755     if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
1756       return false;
1757   }
1758   else
1759     return false;
1760
1761   return true;
1762 }