Salome HOME
22806: EDF SMESH: Regression: Prism_3D error
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
index 239aa41e4be9f406af25451bce4fb401b16db0d8..6311da713a53cdff7413d67c76f886b75066d35a 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2014  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
@@ -6,7 +6,7 @@
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -18,6 +18,7 @@
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
 
 #include "SMESH_ControlsDef.hxx"
 
@@ -32,6 +33,9 @@
 #include "SMESHDS_GroupBase.hxx"
 #include "SMESHDS_Mesh.hxx"
 #include "SMESH_OctreeNode.hxx"
+#include "SMESH_MeshAlgos.hxx"
+
+#include <Basics_Utils.hxx>
 
 #include <BRepAdaptor_Surface.hxx>
 #include <BRepClass_FaceClassifier.hxx>
@@ -45,6 +49,7 @@
 #include <TColStd_SequenceOfAsciiString.hxx>
 #include <TColgp_Array1OfXYZ.hxx>
 #include <TopAbs.hxx>
+#include <TopExp.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Edge.hxx>
 #include <TopoDS_Face.hxx>
@@ -64,7 +69,6 @@
 #include <set>
 #include <limits>
 
-
 /*
                             AUXILIARY METHODS
 */
@@ -210,10 +214,13 @@ using namespace SMESH::Controls;
  *                               FUNCTORS
  */
 
+//================================================================================
 /*
   Class       : NumericalFunctor
   Description : Base class for numerical functors
 */
+//================================================================================
+
 NumericalFunctor::NumericalFunctor():
   myMesh(NULL)
 {
@@ -300,7 +307,7 @@ double NumericalFunctor::GetValue( long theId )
   myCurrElement = myMesh->FindElement( theId );
 
   TSequenceOfXYZ P;
-  if ( GetPoints( theId, P ))
+  if ( GetPoints( theId, P )) // elem type is checked here
     aVal = Round( GetValue( P ));
 
   return aVal;
@@ -321,6 +328,7 @@ double NumericalFunctor::Round( const double & aVal )
  *  \param minmax - boundaries of diapason of values to divide into intervals
  */
 //================================================================================
+
 void NumericalFunctor::GetHistogram(int                  nbIntervals,
                                     std::vector<int>&    nbEvents,
                                     std::vector<double>& funValues,
@@ -403,9 +411,11 @@ void NumericalFunctor::GetHistogram(int                  nbIntervals,
 }
 
 //=======================================================================
-//function : GetValue
-//purpose  : 
-//=======================================================================
+/*
+  Class       : Volume
+  Description : Functor calculating volume of a 3D element
+*/
+//================================================================================
 
 double Volume::GetValue( long theElementId )
 {
@@ -417,21 +427,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;
@@ -442,6 +442,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)
@@ -508,6 +510,7 @@ SMDSAbs_ElementType MaxElementLength2D::GetType() const
   Class       : MaxElementLength3D
   Description : Functor calculating maximum length of 3D element
 */
+//================================================================================
 
 double MaxElementLength3D::GetValue( long theElementId )
 {
@@ -683,6 +686,7 @@ SMDSAbs_ElementType MaxElementLength3D::GetType() const
   Class       : MinimumAngle
   Description : Functor for calculation of minimum angle
 */
+//================================================================================
 
 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
 {
@@ -715,10 +719,13 @@ SMDSAbs_ElementType MinimumAngle::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : AspectRatio
   Description : Functor for calculating aspect ratio
 */
+//================================================================================
+
 double AspectRatio::GetValue( long theId )
 {
   double aVal = 0;
@@ -899,10 +906,13 @@ SMDSAbs_ElementType AspectRatio::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : AspectRatio3D
   Description : Functor for calculating aspect ratio
 */
+//================================================================================
+
 namespace{
 
   inline double getHalfPerimeter(double theTria[3]){
@@ -1269,10 +1279,13 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Warping
   Description : Functor for calculating warping
 */
+//================================================================================
+
 double Warping::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
@@ -1285,7 +1298,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,
@@ -1326,10 +1343,13 @@ SMDSAbs_ElementType Warping::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Taper
   Description : Functor for calculating taper
 */
+//================================================================================
+
 double Taper::GetValue( const TSequenceOfXYZ& P )
 {
   if ( P.size() != 4 )
@@ -1350,7 +1370,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
@@ -1366,11 +1390,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.;
@@ -1388,7 +1414,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 ) ) );
@@ -1408,11 +1434,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;
+
+    const double eps = 0.1; // val is in degrees
 
-    return A * 180. / M_PI;
+    return val < eps ? 0. : val;
   }
 }
 
@@ -1430,10 +1456,13 @@ SMDSAbs_ElementType Skew::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : Area
   Description : Functor for calculating area
 */
+//================================================================================
+
 double Area::GetValue( const TSequenceOfXYZ& P )
 {
   double val = 0.0;
@@ -1463,11 +1492,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() ) {
@@ -1488,10 +1519,12 @@ SMDSAbs_ElementType Length::GetType() const
   return SMDSAbs_Edge;
 }
 
+//================================================================================
 /*
   Class       : Length2D
   Description : Functor for calculating length of edge
 */
+//================================================================================
 
 double Length2D::GetValue( long theElementId)
 {
@@ -1799,10 +1832,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;
@@ -1823,10 +1859,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;
@@ -1915,6 +1954,7 @@ bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) cons
 }
 
 void MultiConnection2D::GetValues(MValues& theValues){
+  if ( !myMesh ) return;
   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
   for(; anIter->more(); ){
     const SMDS_MeshFace* anElem = anIter->next();
@@ -1970,10 +2010,13 @@ void MultiConnection2D::GetValues(MValues& theValues){
 
 }
 
+//================================================================================
 /*
   Class       : BallDiameter
   Description : Functor returning diameter of a ball element
 */
+//================================================================================
+
 double BallDiameter::GetValue( long theId )
 {
   double diameter = 0;
@@ -2002,10 +2045,12 @@ SMDSAbs_ElementType BallDiameter::GetType() const
                             PREDICATES
 */
 
+//================================================================================
 /*
   Class       : BadOrientedVolume
   Description : Predicate bad oriented volumes
 */
+//================================================================================
 
 BadOrientedVolume::BadOrientedVolume()
 {
@@ -2052,9 +2097,11 @@ bool BareBorderVolume::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : BareBorderFace
 */
+//================================================================================
 
 bool BareBorderFace::IsSatisfy(long theElementId )
 {
@@ -2092,9 +2139,11 @@ bool BareBorderFace::IsSatisfy(long theElementId )
   return ok;
 }
 
+//================================================================================
 /*
   Class       : OverConstrainedVolume
 */
+//================================================================================
 
 bool OverConstrainedVolume::IsSatisfy(long theElementId )
 {
@@ -2112,9 +2161,11 @@ bool OverConstrainedVolume::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : OverConstrainedFace
 */
+//================================================================================
 
 bool OverConstrainedFace::IsSatisfy(long theElementId )
 {
@@ -2145,10 +2196,12 @@ bool OverConstrainedFace::IsSatisfy(long theElementId )
   return false;
 }
 
+//================================================================================
 /*
   Class       : CoincidentNodes
   Description : Predicate of Coincident nodes
 */
+//================================================================================
 
 CoincidentNodes::CoincidentNodes()
 {
@@ -2190,11 +2243,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()
 {
@@ -2245,10 +2300,12 @@ SMDSAbs_ElementType CoincidentElements3D::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeBorders
   Description : Predicate for free borders
 */
+//================================================================================
 
 FreeBorders::FreeBorders()
 {
@@ -2271,10 +2328,13 @@ SMDSAbs_ElementType FreeBorders::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeEdges
   Description : Predicate for free Edges
 */
+//================================================================================
+
 FreeEdges::FreeEdges()
 {
   myMesh = 0;
@@ -2406,11 +2466,12 @@ void FreeEdges::GetBoreders(TBorders& theBorders)
   }
 }
 
-
+//================================================================================
 /*
   Class       : FreeNodes
   Description : Predicate for free nodes
 */
+//================================================================================
 
 FreeNodes::FreeNodes()
 {
@@ -2437,10 +2498,12 @@ SMDSAbs_ElementType FreeNodes::GetType() const
 }
 
 
+//================================================================================
 /*
   Class       : FreeFaces
   Description : Predicate for free faces
 */
+//================================================================================
 
 FreeFaces::FreeFaces()
 {
@@ -2493,10 +2556,12 @@ SMDSAbs_ElementType FreeFaces::GetType() const
   return SMDSAbs_Face;
 }
 
+//================================================================================
 /*
   Class       : LinearOrQuadratic
   Description : Predicate to verify whether a mesh element is linear
 */
+//================================================================================
 
 LinearOrQuadratic::LinearOrQuadratic()
 {
@@ -2527,10 +2592,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()
 {
@@ -2599,6 +2666,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' );
@@ -2619,6 +2687,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();
@@ -2627,10 +2696,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()
 {
@@ -2677,6 +2748,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
@@ -2873,11 +3135,14 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
   myIds.Clear();
 
   TCollection_AsciiString aStr = theStr;
-  aStr.RemoveAll( ' ' );
-  aStr.RemoveAll( '\t' );
+  //aStr.RemoveAll( ' ' );
+  //aStr.RemoveAll( '\t' );
+  for ( int i = 1; i <= aStr.Length(); ++i )
+    if ( isspace( aStr.Value( i )))
+      aStr.SetValue( i, ',');
 
   for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
-    aStr.Remove( aPos, 2 );
+    aStr.Remove( aPos, 1 );
 
   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
   int i = 1;
@@ -3148,6 +3413,31 @@ bool LogicalOR::IsSatisfy( long theId )
                               FILTER
 */
 
+// #ifdef WITH_TBB
+// #include <tbb/parallel_for.h>
+// #include <tbb/enumerable_thread_specific.h>
+
+// namespace Parallel
+// {
+//   typedef tbb::enumerable_thread_specific< TIdSequence > TIdSeq;
+
+//   struct Predicate
+//   {
+//     const SMDS_Mesh* myMesh;
+//     PredicatePtr     myPredicate;
+//     TIdSeq &         myOKIds;
+//     Predicate( const SMDS_Mesh* m, PredicatePtr p, TIdSeq & ids ):
+//       myMesh(m), myPredicate(p->Duplicate()), myOKIds(ids) {}
+//     void operator() ( const tbb::blocked_range<size_t>& r ) const
+//     {
+//       for ( size_t i = r.begin(); i != r.end(); ++i )
+//         if ( myPredicate->IsSatisfy( i ))
+//           myOKIds.local().push_back();
+//     }
+//   }
+// }
+// #endif
+
 Filter::Filter()
 {}
 
@@ -3798,8 +4088,15 @@ void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double th
   switch ( myShape.ShapeType() )
   {
   case TopAbs_SOLID: {
-    mySolidClfr.Load(theShape);
-    myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid;
+    if ( isBox( theShape ))
+    {
+      myIsOutFun = & ElementsOnShape::TClassifier::isOutOfBox;
+    }
+    else
+    {
+      mySolidClfr.Load(theShape);
+      myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid;
+    }
     break;
   }
   case TopAbs_FACE:  {
@@ -3833,6 +4130,11 @@ bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p)
   return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON );
 }
 
+bool ElementsOnShape::TClassifier::isOutOfBox (const gp_Pnt& p)
+{
+  return myBox.IsOut( p.XYZ() );
+}
+
 bool ElementsOnShape::TClassifier::isOutOfFace  (const gp_Pnt& p)
 {
   myProjFace.Perform( p );
@@ -3860,6 +4162,390 @@ bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p)
   return ( myVertexXYZ.Distance( p ) > myTol );
 }
 
+bool ElementsOnShape::TClassifier::isBox (const TopoDS_Shape& theShape)
+{
+  TopTools_IndexedMapOfShape vMap;
+  TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
+  if ( vMap.Extent() != 8 )
+    return false;
+
+  myBox.Clear();
+  for ( int i = 1; i <= 8; ++i )
+    myBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))).XYZ() );
+
+  gp_XYZ pMin = myBox.CornerMin(), pMax = myBox.CornerMax();
+  for ( int i = 1; i <= 8; ++i )
+  {
+    gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i )));
+    for ( int iC = 1; iC <= 3; ++ iC )
+    {
+      double d1 = Abs( pMin.Coord( iC ) - p.Coord( iC ));
+      double d2 = Abs( pMax.Coord( iC ) - p.Coord( iC ));
+      if ( Min( d1, d2 ) > myTol )
+        return false;
+    }
+  }
+  myBox.Enlarge( myTol );
+  return true;
+}
+
+
+/*
+  Class       : BelongToGeom
+  Description : Predicate for verifying whether entity belongs to
+                specified geometrical support
+*/
+
+BelongToGeom::BelongToGeom()
+  : myMeshDS(NULL),
+    myType(SMDSAbs_All),
+    myIsSubshape(false),
+    myTolerance(Precision::Confusion())
+{}
+
+void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
+{
+  myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
+  init();
+}
+
+void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
+{
+  myShape = theShape;
+  init();
+}
+
+static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
+                        const TopoDS_Shape& theShape)
+{
+  if (theMap.Contains(theShape)) return true;
+
+  if (theShape.ShapeType() == TopAbs_COMPOUND ||
+      theShape.ShapeType() == TopAbs_COMPSOLID)
+  {
+    TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
+    for (; anIt.More(); anIt.Next())
+    {
+      if (!IsSubShape(theMap, anIt.Value())) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  return false;
+}
+
+void BelongToGeom::init()
+{
+  if (!myMeshDS || myShape.IsNull()) return;
+
+  // is sub-shape of main shape?
+  TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
+  if (aMainShape.IsNull()) {
+    myIsSubshape = false;
+  }
+  else {
+    TopTools_IndexedMapOfShape aMap;
+    TopExp::MapShapes(aMainShape, aMap);
+    myIsSubshape = IsSubShape(aMap, myShape);
+  }
+
+  //if (!myIsSubshape) // to be always ready to check an element not bound to geometry
+  {
+    myElementsOnShapePtr.reset(new ElementsOnShape());
+    myElementsOnShapePtr->SetTolerance(myTolerance);
+    myElementsOnShapePtr->SetAllNodes(true); // "belong", while false means "lays on"
+    myElementsOnShapePtr->SetMesh(myMeshDS);
+    myElementsOnShapePtr->SetShape(myShape, myType);
+  }
+}
+
+static bool IsContains( const SMESHDS_Mesh*     theMeshDS,
+                        const TopoDS_Shape&     theShape,
+                        const SMDS_MeshElement* theElem,
+                        TopAbs_ShapeEnum        theFindShapeEnum,
+                        TopAbs_ShapeEnum        theAvoidShapeEnum = TopAbs_SHAPE )
+{
+  TopExp_Explorer anExp( theShape,theFindShapeEnum,theAvoidShapeEnum );
+
+  while( anExp.More() )
+  {
+    const TopoDS_Shape& aShape = anExp.Current();
+    if( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ) ){
+      if( aSubMesh->Contains( theElem ) )
+        return true;
+    }
+    anExp.Next();
+  }
+  return false;
+}
+
+bool BelongToGeom::IsSatisfy (long theId)
+{
+  if (myMeshDS == 0 || myShape.IsNull())
+    return false;
+
+  if (!myIsSubshape)
+  {
+    return myElementsOnShapePtr->IsSatisfy(theId);
+  }
+
+  // Case of submesh
+  if (myType == SMDSAbs_Node)
+  {
+    if( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ) )
+    {
+      if ( aNode->getshapeId() < 1 )
+        return myElementsOnShapePtr->IsSatisfy(theId);
+
+      const SMDS_PositionPtr& aPosition = aNode->GetPosition();
+      SMDS_TypeOfPosition aTypeOfPosition = aPosition->GetTypeOfPosition();
+      switch( aTypeOfPosition )
+      {
+      case SMDS_TOP_VERTEX : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_VERTEX ));
+      case SMDS_TOP_EDGE   : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_EDGE ));
+      case SMDS_TOP_FACE   : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_FACE ));
+      case SMDS_TOP_3DSPACE: return ( IsContains( myMeshDS,myShape,aNode,TopAbs_SOLID ) ||
+                                      IsContains( myMeshDS,myShape,aNode,TopAbs_SHELL ));
+      }
+    }
+  }
+  else
+  {
+    if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
+    {
+      if ( anElem->getshapeId() < 1 )
+        return myElementsOnShapePtr->IsSatisfy(theId);
+
+      if( myType == SMDSAbs_All )
+      {
+        return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ) ||
+                 IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ) ||
+                 IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
+                 IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
+      }
+      else if( myType == anElem->GetType() )
+      {
+        switch( myType )
+        {
+        case SMDSAbs_Edge  : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ));
+        case SMDSAbs_Face  : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ));
+        case SMDSAbs_Volume: return ( IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
+                                      IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
+        }
+      }
+    }
+  }
+
+  return false;
+}
+
+void BelongToGeom::SetType (SMDSAbs_ElementType theType)
+{
+  myType = theType;
+  init();
+}
+
+SMDSAbs_ElementType BelongToGeom::GetType() const
+{
+  return myType;
+}
+
+TopoDS_Shape BelongToGeom::GetShape()
+{
+  return myShape;
+}
+
+const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const
+{
+  return myMeshDS;
+}
+
+void BelongToGeom::SetTolerance (double theTolerance)
+{
+  myTolerance = theTolerance;
+  if (!myIsSubshape)
+    init();
+}
+
+double BelongToGeom::GetTolerance()
+{
+  return myTolerance;
+}
+
+/*
+  Class       : LyingOnGeom
+  Description : Predicate for verifying whether entiy lying or partially lying on
+                specified geometrical support
+*/
+
+LyingOnGeom::LyingOnGeom()
+  : myMeshDS(NULL),
+    myType(SMDSAbs_All),
+    myIsSubshape(false),
+    myTolerance(Precision::Confusion())
+{}
+
+void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
+{
+  myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
+  init();
+}
+
+void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
+{
+  myShape = theShape;
+  init();
+}
+
+void LyingOnGeom::init()
+{
+  if (!myMeshDS || myShape.IsNull()) return;
+
+  // is sub-shape of main shape?
+  TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
+  if (aMainShape.IsNull()) {
+    myIsSubshape = false;
+  }
+  else {
+    TopTools_IndexedMapOfShape aMap;
+    TopExp::MapShapes(aMainShape, aMap);
+    myIsSubshape = IsSubShape(aMap, myShape);
+  }
+
+  if (!myIsSubshape)
+  {
+    myElementsOnShapePtr.reset(new ElementsOnShape());
+    myElementsOnShapePtr->SetTolerance(myTolerance);
+    myElementsOnShapePtr->SetAllNodes(false); // lays on, while true means "belong"
+    myElementsOnShapePtr->SetMesh(myMeshDS);
+    myElementsOnShapePtr->SetShape(myShape, myType);
+  }
+}
+
+bool LyingOnGeom::IsSatisfy( long theId )
+{
+  if ( myMeshDS == 0 || myShape.IsNull() )
+    return false;
+
+  if (!myIsSubshape)
+  {
+    return myElementsOnShapePtr->IsSatisfy(theId);
+  }
+
+  // Case of submesh
+  if( myType == SMDSAbs_Node )
+  {
+    if( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ) )
+    {
+      const SMDS_PositionPtr& aPosition = aNode->GetPosition();
+      SMDS_TypeOfPosition aTypeOfPosition = aPosition->GetTypeOfPosition();
+      switch( aTypeOfPosition )
+      {
+      case SMDS_TOP_VERTEX : return IsContains( myMeshDS,myShape,aNode,TopAbs_VERTEX );
+      case SMDS_TOP_EDGE   : return IsContains( myMeshDS,myShape,aNode,TopAbs_EDGE );
+      case SMDS_TOP_FACE   : return IsContains( myMeshDS,myShape,aNode,TopAbs_FACE );
+      case SMDS_TOP_3DSPACE: return IsContains( myMeshDS,myShape,aNode,TopAbs_SHELL );
+      }
+    }
+  }
+  else
+  {
+    if( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ) )
+    {
+      if( myType == SMDSAbs_All )
+      {
+        return Contains( myMeshDS,myShape,anElem,TopAbs_EDGE ) ||
+               Contains( myMeshDS,myShape,anElem,TopAbs_FACE ) ||
+               Contains( myMeshDS,myShape,anElem,TopAbs_SHELL )||
+               Contains( myMeshDS,myShape,anElem,TopAbs_SOLID );
+      }
+      else if( myType == anElem->GetType() )
+      {
+        switch( myType )
+        {
+        case SMDSAbs_Edge  : return Contains( myMeshDS,myShape,anElem,TopAbs_EDGE );
+        case SMDSAbs_Face  : return Contains( myMeshDS,myShape,anElem,TopAbs_FACE );
+        case SMDSAbs_Volume: return Contains( myMeshDS,myShape,anElem,TopAbs_SHELL )||
+                                    Contains( myMeshDS,myShape,anElem,TopAbs_SOLID );
+        }
+      }
+    }
+  }
+
+  return false;
+}
+
+void LyingOnGeom::SetType( SMDSAbs_ElementType theType )
+{
+  myType = theType;
+  init();
+}
+
+SMDSAbs_ElementType LyingOnGeom::GetType() const
+{
+  return myType;
+}
+
+TopoDS_Shape LyingOnGeom::GetShape()
+{
+  return myShape;
+}
+
+const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const
+{
+  return myMeshDS;
+}
+
+void LyingOnGeom::SetTolerance (double theTolerance)
+{
+  myTolerance = theTolerance;
+  if (!myIsSubshape)
+    init();
+}
+
+double LyingOnGeom::GetTolerance()
+{
+  return myTolerance;
+}
+
+bool LyingOnGeom::Contains( const SMESHDS_Mesh*     theMeshDS,
+                            const TopoDS_Shape&     theShape,
+                            const SMDS_MeshElement* theElem,
+                            TopAbs_ShapeEnum        theFindShapeEnum,
+                            TopAbs_ShapeEnum        theAvoidShapeEnum )
+{
+  if (IsContains(theMeshDS, theShape, theElem, theFindShapeEnum, theAvoidShapeEnum))
+    return true;
+
+  TopTools_IndexedMapOfShape aSubShapes;
+  TopExp::MapShapes( theShape, aSubShapes );
+
+  for (int i = 1; i <= aSubShapes.Extent(); i++)
+  {
+    const TopoDS_Shape& aShape = aSubShapes.FindKey(i);
+
+    if( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ) ){
+      if( aSubMesh->Contains( theElem ) )
+        return true;
+
+      SMDS_NodeIteratorPtr aNodeIt = aSubMesh->GetNodes();
+      while ( aNodeIt->more() )
+      {
+        const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(aNodeIt->next());
+        SMDS_ElemIteratorPtr anElemIt = aNode->GetInverseElementIterator();
+        while ( anElemIt->more() )
+        {
+          const SMDS_MeshElement* anElement = static_cast<const SMDS_MeshElement*>(anElemIt->next());
+          if (anElement == theElem)
+            return true;
+        }
+      }
+    }
+  }
+  return false;
+}
 
 TSequenceOfXYZ::TSequenceOfXYZ()
 {}