Salome HOME
Make ElemEntityType pass nodes through as well
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
index 3e44afe282f09f0a2c53c05b5466a46986c87f9d..7686e655db6ba1215f01d19318b80dcf844a738a 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -33,6 +33,7 @@
 #include "SMESHDS_GroupBase.hxx"
 #include "SMESHDS_Mesh.hxx"
 #include "SMESH_OctreeNode.hxx"
+#include "SMESH_MeshAlgos.hxx"
 
 #include <BRepAdaptor_Surface.hxx>
 #include <BRepClass_FaceClassifier.hxx>
@@ -65,6 +66,7 @@
 #include <set>
 #include <limits>
 
+#include <Basics_Utils.hxx>
 
 /*
                             AUXILIARY METHODS
@@ -211,10 +213,13 @@ using namespace SMESH::Controls;
  *                               FUNCTORS
  */
 
+//================================================================================
 /*
   Class       : NumericalFunctor
   Description : Base class for numerical functors
 */
+//================================================================================
+
 NumericalFunctor::NumericalFunctor():
   myMesh(NULL)
 {
@@ -327,7 +332,8 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
                                     std::vector<int>&    nbEvents,
                                     std::vector<double>& funValues,
                                     const vector<int>&   elements,
-                                    const double*        minmax)
+                                    const double*        minmax,
+                                    const bool           isLogarithmic)
 {
   if ( nbIntervals < 1 ||
        !myMesh ||
@@ -380,8 +386,15 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
   for ( int i = 0; i < nbIntervals; ++i )
   {
     // find end value of i-th interval
-    double r = (i+1) / double( nbIntervals );
-    funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
+    double r = (i+1) / double(nbIntervals);
+    if (isLogarithmic && funValues.front() > 1e-07 && funValues.back() > 1e-07) {
+      double logmin = log10(funValues.front());
+      double lval = logmin + r * (log10(funValues.back()) - logmin);
+      funValues[i+1] = pow(10.0, lval);
+    }
+    else {
+      funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
+    }
 
     // count values in the i-th interval if there are any
     if ( min != values.end() && *min <= funValues[i+1] )
@@ -397,9 +410,11 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
 }
 
 //=======================================================================
-//function : GetValue
-//purpose  : 
-//=======================================================================
+/*
+  Class       : Volume
+  Description : Functor calculating volume of a 3D element
+*/
+//================================================================================
 
 double Volume::GetValue( long theElementId )
 {
@@ -411,21 +426,11 @@ double Volume::GetValue( long theElementId )
   return 0;
 }
 
-//=======================================================================
-//function : GetBadRate
-//purpose  : meaningless as it is not quality control functor
-//=======================================================================
-
 double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
 {
   return Value;
 }
 
-//=======================================================================
-//function : GetType
-//purpose  : 
-//=======================================================================
-
 SMDSAbs_ElementType Volume::GetType() const
 {
   return SMDSAbs_Volume;
@@ -436,6 +441,8 @@ SMDSAbs_ElementType Volume::GetType() const
   Class       : MaxElementLength2D
   Description : Functor calculating maximum length of 2D element
 */
+//================================================================================
+
 double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P )
 {
   if(P.size() == 0)
@@ -502,6 +509,7 @@ SMDSAbs_ElementType MaxElementLength2D::GetType() const
   Class       : MaxElementLength3D
   Description : Functor calculating maximum length of 3D element
 */
+//================================================================================
 
 double MaxElementLength3D::GetValue( long theElementId )
 {
@@ -677,6 +685,7 @@ SMDSAbs_ElementType MaxElementLength3D::GetType() const
   Class       : MinimumAngle
   Description : Functor for calculation of minimum angle
 */
+//================================================================================
 
 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
 {
@@ -709,10 +718,13 @@ SMDSAbs_ElementType MinimumAngle::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : AspectRatio
   Description : Functor for calculating aspect ratio
 */
+//================================================================================
+
 double AspectRatio::GetValue( long theId )
 {
   double aVal = 0;
@@ -893,10 +905,13 @@ SMDSAbs_ElementType AspectRatio::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : AspectRatio3D
   Description : Functor for calculating aspect ratio
 */
+//================================================================================
+
 namespace{
 
   inline double getHalfPerimeter(double theTria[3]){
@@ -1263,10 +1278,13 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Warping
   Description : Functor for calculating warping
 */
+//================================================================================
+
 double Warping::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
@@ -1279,7 +1297,11 @@ double Warping::GetValue( const TSequenceOfXYZ& P )
   double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
   double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
 
-  return Max( Max( A1, A2 ), Max( A3, A4 ) );
+  double val = Max( Max( A1, A2 ), Max( A3, A4 ) );
+
+  const double eps = 0.1; // val is in degrees
+
+  return val < eps ? 0. : val;
 }
 
 double Warping::ComputeA( const gp_XYZ& thePnt1,
@@ -1320,10 +1342,13 @@ SMDSAbs_ElementType Warping::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Taper
   Description : Functor for calculating taper
 */
+//================================================================================
+
 double Taper::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
@@ -1344,7 +1369,11 @@ double Taper::GetValue( const TSequenceOfXYZ& P )
   double T3 = fabs( ( J3 - JA ) / JA );
   double T4 = fabs( ( J4 - JA ) / JA );
 
-  return Max( Max( T1, T2 ), Max( T3, T4 ) );
+  double val = Max( Max( T1, T2 ), Max( T3, T4 ) );
+
+  const double eps = 0.01;
+
+  return val < eps ? 0. : val;
 }
 
 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
@@ -1360,11 +1389,13 @@ SMDSAbs_ElementType Taper::GetType() const
   return SMDSAbs_Face;
 }
 
-
+//================================================================================
 /*
   Class       : Skew
   Description : Functor for calculating skew in degrees
 */
+//================================================================================
+
 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
 {
   gp_XYZ p12 = ( p2 + p1 ) / 2.;
@@ -1382,7 +1413,7 @@ double Skew::GetValue( const TSequenceOfXYZ& P )
     return 0.;
 
   // Compute skew
-  static double PI2 = M_PI / 2.;
+  const double PI2 = M_PI / 2.;
   if ( P.size() == 3 )
   {
     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
@@ -1402,11 +1433,11 @@ double Skew::GetValue( const TSequenceOfXYZ& P )
     double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
       ? 0. : fabs( PI2 - v1.Angle( v2 ) );
 
-    //BUG SWP12743
-    if ( A < theEps )
-      return theInf;
+    double val = A * 180. / M_PI;
 
-    return A * 180. / M_PI;
+    const double eps = 0.1; // val is in degrees
+
+    return val < eps ? 0. : val;
   }
 }
 
@@ -1424,10 +1455,13 @@ SMDSAbs_ElementType Skew::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Area
   Description : Functor for calculating area
 */
+//================================================================================
+
 double Area::GetValue( const TSequenceOfXYZ& P )
 {
   double val = 0.0;
@@ -1457,11 +1491,13 @@ SMDSAbs_ElementType Area::GetType() const
   return SMDSAbs_Face;
 }
 
-
+//================================================================================
 /*
   Class       : Length
   Description : Functor for calculating length of edge
 */
+//================================================================================
+
 double Length::GetValue( const TSequenceOfXYZ& P )
 {
   switch ( P.size() ) {
@@ -1482,10 +1518,12 @@ SMDSAbs_ElementType Length::GetType() const
   return SMDSAbs_Edge;
 }
 
+//================================================================================
 /*
   Class       : Length2D
   Description : Functor for calculating length of edge
 */
+//================================================================================
 
 double Length2D::GetValue( long theElementId)
 {
@@ -1793,10 +1831,13 @@ void Length2D::GetValues(TValues& theValues){
   }
 }
 
+//================================================================================
 /*
   Class       : MultiConnection
   Description : Functor for calculating number of faces conneted to the edge
 */
+//================================================================================
+
 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
 {
   return 0;
@@ -1817,10 +1858,13 @@ SMDSAbs_ElementType MultiConnection::GetType() const
   return SMDSAbs_Edge;
 }
 
+//================================================================================
 /*
   Class       : MultiConnection2D
   Description : Functor for calculating number of faces conneted to the edge
 */
+//================================================================================
+
 double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
 {
   return 0;
@@ -1964,10 +2008,13 @@ void MultiConnection2D::GetValues(MValues& theValues){
 
 }
 
+//================================================================================
 /*
   Class       : BallDiameter
   Description : Functor returning diameter of a ball element
 */
+//================================================================================
+
 double BallDiameter::GetValue( long theId )
 {
   double diameter = 0;
@@ -1996,10 +2043,12 @@ SMDSAbs_ElementType BallDiameter::GetType() const
                             PREDICATES
 */
 
+//================================================================================
 /*
   Class       : BadOrientedVolume
   Description : Predicate bad oriented volumes
 */
+//================================================================================
 
 BadOrientedVolume::BadOrientedVolume()
 {
@@ -2046,9 +2095,11 @@ bool BareBorderVolume::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : BareBorderFace
 */
+//================================================================================
 
 bool BareBorderFace::IsSatisfy(long theElementId )
 {
@@ -2086,9 +2137,11 @@ bool BareBorderFace::IsSatisfy(long theElementId )
   return ok;
 }
 
+//================================================================================
 /*
   Class       : OverConstrainedVolume
 */
+//================================================================================
 
 bool OverConstrainedVolume::IsSatisfy(long theElementId )
 {
@@ -2106,9 +2159,11 @@ bool OverConstrainedVolume::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : OverConstrainedFace
 */
+//================================================================================
 
 bool OverConstrainedFace::IsSatisfy(long theElementId )
 {
@@ -2139,10 +2194,12 @@ bool OverConstrainedFace::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : CoincidentNodes
   Description : Predicate of Coincident nodes
 */
+//================================================================================
 
 CoincidentNodes::CoincidentNodes()
 {
@@ -2184,11 +2241,13 @@ void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
   }
 }
 
+//================================================================================
 /*
   Class       : CoincidentElements
   Description : Predicate of Coincident Elements
   Note        : This class is suitable only for visualization of Coincident Elements
 */
+//================================================================================
 
 CoincidentElements::CoincidentElements()
 {
@@ -2239,10 +2298,12 @@ SMDSAbs_ElementType CoincidentElements3D::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeBorders
   Description : Predicate for free borders
 */
+//================================================================================
 
 FreeBorders::FreeBorders()
 {
@@ -2265,10 +2326,13 @@ SMDSAbs_ElementType FreeBorders::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeEdges
   Description : Predicate for free Edges
 */
+//================================================================================
+
 FreeEdges::FreeEdges()
 {
   myMesh = 0;
@@ -2400,11 +2464,12 @@ void FreeEdges::GetBoreders(TBorders& theBorders)
   }
 }
 
-
+//================================================================================
 /*
   Class       : FreeNodes
   Description : Predicate for free nodes
 */
+//================================================================================
 
 FreeNodes::FreeNodes()
 {
@@ -2431,10 +2496,12 @@ SMDSAbs_ElementType FreeNodes::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeFaces
   Description : Predicate for free faces
 */
+//================================================================================
 
 FreeFaces::FreeFaces()
 {
@@ -2487,10 +2554,12 @@ SMDSAbs_ElementType FreeFaces::GetType() const
   return SMDSAbs_Face;
 }
 
+//================================================================================
 /*
   Class       : LinearOrQuadratic
   Description : Predicate to verify whether a mesh element is linear
 */
+//================================================================================
 
 LinearOrQuadratic::LinearOrQuadratic()
 {
@@ -2521,10 +2590,12 @@ SMDSAbs_ElementType LinearOrQuadratic::GetType() const
   return myType;
 }
 
+//================================================================================
 /*
   Class       : GroupColor
   Description : Functor for check color of group to whic mesh element belongs to
 */
+//================================================================================
 
 GroupColor::GroupColor()
 {
@@ -2593,6 +2664,7 @@ void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
 
 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
 {
+  Kernel_Utils::Localizer loc;
   TCollection_AsciiString aStr = theStr;
   aStr.RemoveAll( ' ' );
   aStr.RemoveAll( '\t' );
@@ -2613,6 +2685,7 @@ void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
 // Purpose : Get range as a string.
 //           Example: "1,2,3,50-60,63,67,70-"
 //=======================================================================
+
 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
 {
   theResStr.Clear();
@@ -2621,10 +2694,12 @@ void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
 }
 
+//================================================================================
 /*
   Class       : ElemGeomType
   Description : Predicate to check element geometry type
 */
+//================================================================================
 
 ElemGeomType::ElemGeomType()
 {
@@ -2671,6 +2746,197 @@ SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
   return myGeomType;
 }
 
+//================================================================================
+/*
+  Class       : ElemEntityType
+  Description : Predicate to check element entity type
+*/
+//================================================================================
+
+ElemEntityType::ElemEntityType():
+  myMesh( 0 ),
+  myType( SMDSAbs_All ),
+  myEntityType( SMDSEntity_0D )
+{
+}
+
+void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh )
+{
+  myMesh = theMesh;
+}
+
+bool ElemEntityType::IsSatisfy( long theId )
+{
+  if ( !myMesh ) return false;
+  if ( myType == SMDSAbs_Node )
+    return myMesh->FindNode( theId );
+  const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
+  return ( anElem &&
+           myEntityType == anElem->GetEntityType() );
+}
+
+void ElemEntityType::SetType( SMDSAbs_ElementType theType )
+{
+  myType = theType;
+}
+
+SMDSAbs_ElementType ElemEntityType::GetType() const
+{
+  return myType;
+}
+
+void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType )
+{
+  myEntityType = theEntityType;
+}
+
+SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const
+{
+  return myEntityType;
+}
+
+//================================================================================
+/*!
+ * \brief Class ConnectedElements
+ */
+//================================================================================
+
+ConnectedElements::ConnectedElements():
+  myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {}
+
+SMDSAbs_ElementType ConnectedElements::GetType() const
+{ return myType; }
+
+int ConnectedElements::GetNode() const
+{ return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
+
+std::vector<double> ConnectedElements::GetPoint() const
+{ return myXYZ; }
+
+void ConnectedElements::clearOkIDs()
+{ myOkIDsReady = false; myOkIDs.clear(); }
+
+void ConnectedElements::SetType( SMDSAbs_ElementType theType )
+{
+  if ( myType != theType || myMeshModifTracer.IsMeshModified() )
+    clearOkIDs();
+  myType = theType;
+}
+
+void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh )
+{
+  myMeshModifTracer.SetMesh( theMesh );
+  if ( myMeshModifTracer.IsMeshModified() )
+  {
+    clearOkIDs();
+    if ( !myXYZ.empty() )
+      SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh
+  }
+}
+
+void ConnectedElements::SetNode( int nodeID )
+{
+  myNodeID = nodeID;
+  myXYZ.clear();
+
+  bool isSameDomain = false;
+  if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() )
+    if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID ))
+    {
+      SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType );
+      while ( !isSameDomain && eIt->more() )
+        isSameDomain = IsSatisfy( eIt->next()->GetID() );
+    }
+  if ( !isSameDomain )
+    clearOkIDs();
+}
+
+void ConnectedElements::SetPoint( double x, double y, double z )
+{
+  myXYZ.resize(3);
+  myXYZ[0] = x;
+  myXYZ[1] = y;
+  myXYZ[2] = z;
+  myNodeID = 0;
+
+  bool isSameDomain = false;
+
+  // find myNodeID by myXYZ if possible
+  if ( myMeshModifTracer.GetMesh() )
+  {
+    auto_ptr<SMESH_ElementSearcher> searcher
+      ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
+
+    vector< const SMDS_MeshElement* > foundElems;
+    searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
+
+    if ( !foundElems.empty() )
+    {
+      myNodeID = foundElems[0]->GetNode(0)->GetID();
+      if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() )
+        isSameDomain = IsSatisfy( foundElems[0]->GetID() );
+    }
+  }
+  if ( !isSameDomain )
+    clearOkIDs();
+}
+
+bool ConnectedElements::IsSatisfy( long theElementId )
+{
+  // Here we do NOT check if the mesh has changed, we do it in Set...() only!!!
+
+  if ( !myOkIDsReady )
+  {
+    if ( !myMeshModifTracer.GetMesh() )
+      return false;
+    const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID );
+    if ( !node0 )
+      return false;
+
+    list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
+    std::set< int > checkedNodeIDs;
+    // algo:
+    // foreach node in nodeQueue:
+    //   foreach element sharing a node:
+    //     add ID of an element of myType to myOkIDs;
+    //     push all element nodes absent from checkedNodeIDs to nodeQueue;
+    while ( !nodeQueue.empty() )
+    {
+      const SMDS_MeshNode* node = nodeQueue.front();
+      nodeQueue.pop_front();
+
+      // loop on elements sharing the node
+      SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
+      while ( eIt->more() )
+      {
+        // keep elements of myType
+        const SMDS_MeshElement* element = eIt->next();
+        if ( element->GetType() == myType )
+          myOkIDs.insert( myOkIDs.end(), element->GetID() );
+
+        // enqueue nodes of the element
+        SMDS_ElemIteratorPtr nIt = element->nodesIterator();
+        while ( nIt->more() )
+        {
+          const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() );
+          if ( checkedNodeIDs.insert( n->GetID() ).second )
+            nodeQueue.push_back( n );
+        }
+      }
+    }
+    if ( myType == SMDSAbs_Node )
+      std::swap( myOkIDs, checkedNodeIDs );
+
+    size_t totalNbElems = myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType );
+    if ( myOkIDs.size() == totalNbElems )
+      myOkIDs.clear();
+
+    myOkIDsReady = true;
+  }
+
+  return myOkIDs.empty() ? true : myOkIDs.count( theElementId );
+}
+
 //================================================================================
 /*!
  * \brief Class CoplanarFaces
@@ -3671,23 +3937,11 @@ ElementsOnShape::ElementsOnShape()
     myToler(Precision::Confusion()),
     myAllNodesFlag(false)
 {
-  myCurShapeType = TopAbs_SHAPE;
 }
 
 ElementsOnShape::~ElementsOnShape()
 {
-}
-
-void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
-{
-  myMeshModifTracer.SetMesh( theMesh );
-  if ( myMeshModifTracer.IsMeshModified())
-    SetShape(myShape, myType);
-}
-
-bool ElementsOnShape::IsSatisfy (long theElementId)
-{
-  return myIds.Contains(theElementId);
+  clearClassifiers();
 }
 
 SMDSAbs_ElementType ElementsOnShape::GetType() const
@@ -3710,169 +3964,163 @@ double ElementsOnShape::GetTolerance() const
 
 void ElementsOnShape::SetAllNodes (bool theAllNodes)
 {
-  if (myAllNodesFlag != theAllNodes) {
-    myAllNodesFlag = theAllNodes;
-    SetShape(myShape, myType);
-  }
+  myAllNodesFlag = theAllNodes;
+}
+
+void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
+{
+  myMesh = theMesh;
 }
 
 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
                                 const SMDSAbs_ElementType theType)
 {
-  myType = theType;
+  myType  = theType;
   myShape = theShape;
-  myIds.Clear();
-
-  const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh();
+  if ( myShape.IsNull() ) return;
   
-  if ( !myMesh ) return;
-
-  myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
+  TopTools_IndexedMapOfShape shapesMap;
+  TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
+  TopExp_Explorer sub;
+  for ( int i = 0; i < 4; ++i )
+  {
+    if ( shapesMap.IsEmpty() )
+      for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() )
+        shapesMap.Add( sub.Current() );
+    if ( i > 0 )
+      for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() )
+        shapesMap.Add( sub.Current() );
+  }
 
-  myShapesMap.Clear();
-  addShape(myShape);
+  clearClassifiers();
+  myClassifiers.resize( shapesMap.Extent() );
+  for ( int i = 0; i < shapesMap.Extent(); ++i )
+    myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler );
 }
 
-void ElementsOnShape::addShape (const TopoDS_Shape& theShape)
+void ElementsOnShape::clearClassifiers()
 {
-  if (theShape.IsNull() || myMeshModifTracer.GetMesh() == 0)
-    return;
+  for ( size_t i = 0; i < myClassifiers.size(); ++i )
+    delete myClassifiers[ i ];
+  myClassifiers.clear();
+}
 
-  if (!myShapesMap.Add(theShape)) return;
+bool ElementsOnShape::IsSatisfy (long elemId)
+{
+  const SMDS_MeshElement* elem =
+    ( myType == SMDSAbs_Node ? myMesh->FindNode( elemId ) : myMesh->FindElement( elemId ));
+  if ( !elem || myClassifiers.empty() )
+    return false;
 
-  myCurShapeType = theShape.ShapeType();
-  switch (myCurShapeType)
+  for ( size_t i = 0; i < myClassifiers.size(); ++i )
   {
-  case TopAbs_COMPOUND:
-  case TopAbs_COMPSOLID:
-  case TopAbs_SHELL:
-  case TopAbs_WIRE:
+    SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
+    bool isSatisfy = myAllNodesFlag;
+    
+    gp_XYZ centerXYZ (0, 0, 0);
+
+    while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
     {
-      TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
-      for (; anIt.More(); anIt.Next()) addShape(anIt.Value());
+      SMESH_TNodeXYZ aPnt ( aNodeItr->next() );
+      centerXYZ += aPnt;
+      isSatisfy = ! myClassifiers[i]->IsOut( aPnt );
     }
-    break;
-  case TopAbs_SOLID:
+
+    // Check the center point for volumes MantisBug 0020168
+    if (isSatisfy &&
+        myAllNodesFlag &&
+        myClassifiers[i]->ShapeType() == TopAbs_SOLID)
     {
-      myCurSC.Load(theShape);
-      process();
+      centerXYZ /= elem->NbNodes();
+      isSatisfy = ! myClassifiers[i]->IsOut( centerXYZ );
     }
+    if ( isSatisfy )
+      return true;
+  }
+
+  return false;
+}
+
+TopAbs_ShapeEnum ElementsOnShape::TClassifier::ShapeType() const
+{
+  return myShape.ShapeType();
+}
+
+bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p)
+{
+  return (this->*myIsOutFun)( p );
+}
+
+void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double theTol)
+{
+  myShape = theShape;
+  myTol   = theTol;
+  switch ( myShape.ShapeType() )
+  {
+  case TopAbs_SOLID: {
+    mySolidClfr.Load(theShape);
+    myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid;
     break;
-  case TopAbs_FACE:
-    {
-      TopoDS_Face aFace = TopoDS::Face(theShape);
-      BRepAdaptor_Surface SA (aFace, true);
-      Standard_Real
-        u1 = SA.FirstUParameter(),
-        u2 = SA.LastUParameter(),
-        v1 = SA.FirstVParameter(),
-        v2 = SA.LastVParameter();
-      Handle(Geom_Surface) surf = BRep_Tool::Surface(aFace);
-      myCurProjFace.Init(surf, u1,u2, v1,v2);
-      myCurFace = aFace;
-      process();
-    }
+  }
+  case TopAbs_FACE:  {
+    Standard_Real u1,u2,v1,v2;
+    Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
+    surf->Bounds( u1,u2,v1,v2 );
+    myProjFace.Init(surf, u1,u2, v1,v2, myTol );
+    myIsOutFun = & ElementsOnShape::TClassifier::isOutOfFace;
     break;
-  case TopAbs_EDGE:
-    {
-      TopoDS_Edge anEdge = TopoDS::Edge(theShape);
-      Standard_Real u1, u2;
-      Handle(Geom_Curve) curve = BRep_Tool::Curve(anEdge, u1, u2);
-      myCurProjEdge.Init(curve, u1, u2);
-      process();
-    }
+  }
+  case TopAbs_EDGE:  {
+    Standard_Real u1, u2;
+    Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge(theShape), u1, u2);
+    myProjEdge.Init(curve, u1, u2);
+    myIsOutFun = & ElementsOnShape::TClassifier::isOutOfEdge;
     break;
-  case TopAbs_VERTEX:
-    {
-      TopoDS_Vertex aV = TopoDS::Vertex(theShape);
-      myCurPnt = BRep_Tool::Pnt(aV);
-      process();
-    }
+  }
+  case TopAbs_VERTEX:{
+    myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
+    myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex;
     break;
+  }
   default:
-    break;
+    throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier");
   }
 }
 
-void ElementsOnShape::process()
+bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p)
 {
-  const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh();
-  if (myShape.IsNull() || myMesh == 0)
-    return;
-
-  SMDS_ElemIteratorPtr anIter = myMesh->elementsIterator(myType);
-  while (anIter->more())
-    process(anIter->next());
+  mySolidClfr.Perform( p, myTol );
+  return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON );
 }
 
-void ElementsOnShape::process (const SMDS_MeshElement* theElemPtr)
+bool ElementsOnShape::TClassifier::isOutOfFace  (const gp_Pnt& p)
 {
-  if (myShape.IsNull())
-    return;
-
-  SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
-  bool isSatisfy = myAllNodesFlag;
-
-  gp_XYZ centerXYZ (0, 0, 0);
-
-  while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
+  myProjFace.Perform( p );
+  if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
   {
-    SMESH_TNodeXYZ aPnt ( aNodeItr->next() );
-    centerXYZ += aPnt;
-
-    switch (myCurShapeType)
-    {
-    case TopAbs_SOLID:
-      {
-        myCurSC.Perform(aPnt, myToler);
-        isSatisfy = (myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON);
-      }
-      break;
-    case TopAbs_FACE:
-      {
-        myCurProjFace.Perform(aPnt);
-        isSatisfy = (myCurProjFace.IsDone() && myCurProjFace.LowerDistance() <= myToler);
-        if (isSatisfy)
-        {
-          // check relatively the face
-          Quantity_Parameter u, v;
-          myCurProjFace.LowerDistanceParameters(u, v);
-          gp_Pnt2d aProjPnt (u, v);
-          BRepClass_FaceClassifier aClsf (myCurFace, aProjPnt, myToler);
-          isSatisfy = (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON);
-        }
-      }
-      break;
-    case TopAbs_EDGE:
-      {
-        myCurProjEdge.Perform(aPnt);
-        isSatisfy = (myCurProjEdge.NbPoints() > 0 && myCurProjEdge.LowerDistance() <= myToler);
-      }
-      break;
-    case TopAbs_VERTEX:
-      {
-        isSatisfy = (myCurPnt.Distance(aPnt) <= myToler);
-      }
-      break;
-    default:
-      {
-        isSatisfy = false;
-      }
-    }
+    // check relatively to the face
+    Quantity_Parameter u, v;
+    myProjFace.LowerDistanceParameters(u, v);
+    gp_Pnt2d aProjPnt (u, v);
+    BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
+    if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
+      return false;
   }
+  return true;
+}
 
-  if (isSatisfy && myCurShapeType == TopAbs_SOLID) { // Check the center point for volumes MantisBug 0020168
-    centerXYZ /= theElemPtr->NbNodes();
-    gp_Pnt aCenterPnt (centerXYZ);
-    myCurSC.Perform(aCenterPnt, myToler);
-    if ( !(myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON))
-      isSatisfy = false;
-  }
+bool ElementsOnShape::TClassifier::isOutOfEdge  (const gp_Pnt& p)
+{
+  myProjEdge.Perform( p );
+  return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
+}
 
-  if (isSatisfy)
-    myIds.Add(theElemPtr->GetID());
+bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p)
+{
+  return ( myVertexXYZ.Distance( p ) > myTol );
 }
 
+
 TSequenceOfXYZ::TSequenceOfXYZ()
 {}