X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=42c668044def6d2b04df05bb2536ab948c218c0e;hp=3e44afe282f09f0a2c53c05b5466a46986c87f9d;hb=51b0f287860107c56c08e61d40a3aea9ba425483;hpb=1067ffa6e7e5c394e3a1b17219d8b355a57607cd diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 3e44afe28..42c668044 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -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 @@ -33,6 +33,9 @@ #include "SMESHDS_GroupBase.hxx" #include "SMESHDS_Mesh.hxx" #include "SMESH_OctreeNode.hxx" +#include "SMESH_MeshAlgos.hxx" + +#include #include #include @@ -46,6 +49,7 @@ #include #include #include +#include #include #include #include @@ -65,7 +69,6 @@ #include #include - /* AUXILIARY METHODS */ @@ -211,10 +214,13 @@ using namespace SMESH::Controls; * FUNCTORS */ +//================================================================================ /* Class : NumericalFunctor Description : Base class for numerical functors */ +//================================================================================ + NumericalFunctor::NumericalFunctor(): myMesh(NULL) { @@ -301,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; @@ -327,7 +333,8 @@ void NumericalFunctor::GetHistogram(int nbIntervals, std::vector& nbEvents, std::vector& funValues, const vector& elements, - const double* minmax) + const double* minmax, + const bool isLogarithmic) { if ( nbIntervals < 1 || !myMesh || @@ -380,8 +387,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 +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 ) { @@ -411,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; @@ -436,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) @@ -502,6 +510,7 @@ SMDSAbs_ElementType MaxElementLength2D::GetType() const Class : MaxElementLength3D Description : Functor calculating maximum length of 3D element */ +//================================================================================ double MaxElementLength3D::GetValue( long theElementId ) { @@ -677,6 +686,7 @@ SMDSAbs_ElementType MaxElementLength3D::GetType() const Class : MinimumAngle Description : Functor for calculation of minimum angle */ +//================================================================================ double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) { @@ -709,10 +719,13 @@ SMDSAbs_ElementType MinimumAngle::GetType() const } +//================================================================================ /* Class : AspectRatio Description : Functor for calculating aspect ratio */ +//================================================================================ + double AspectRatio::GetValue( long theId ) { double aVal = 0; @@ -893,10 +906,13 @@ SMDSAbs_ElementType AspectRatio::GetType() const } +//================================================================================ /* Class : AspectRatio3D Description : Functor for calculating aspect ratio */ +//================================================================================ + namespace{ inline double getHalfPerimeter(double theTria[3]){ @@ -1263,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 ) @@ -1279,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, @@ -1320,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 ) @@ -1344,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 @@ -1360,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.; @@ -1382,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 ) ) ); @@ -1402,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; - return A * 180. / M_PI; + const double eps = 0.1; // val is in degrees + + return val < eps ? 0. : val; } } @@ -1424,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; @@ -1457,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() ) { @@ -1482,12 +1519,14 @@ SMDSAbs_ElementType Length::GetType() const return SMDSAbs_Edge; } +//================================================================================ /* Class : Length2D Description : Functor for calculating length of edge */ +//================================================================================ -double Length2D::GetValue( long theElementId) +double Length2D::GetValue( long theElementId ) { TSequenceOfXYZ P; @@ -1519,7 +1558,7 @@ double Length2D::GetValue( long theElementId) double L1 = getDistance(P( 1 ),P( 2 )); double L2 = getDistance(P( 2 ),P( 3 )); double L3 = getDistance(P( 3 ),P( 1 )); - aVal = Max(L1,Max(L2,L3)); + aVal = Min(L1,Min(L2,L3)); break; } else if (len == 4){ // quadrangles @@ -1527,14 +1566,14 @@ double Length2D::GetValue( long theElementId) double L2 = getDistance(P( 2 ),P( 3 )); double L3 = getDistance(P( 3 ),P( 4 )); double L4 = getDistance(P( 4 ),P( 1 )); - aVal = Max(Max(L1,L2),Max(L3,L4)); + aVal = Min(Min(L1,L2),Min(L3,L4)); break; } if (len == 6){ // quadratic triangles double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 )); double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 )); double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 )); - aVal = Max(L1,Max(L2,L3)); + aVal = Min(L1,Min(L2,L3)); //cout<<"L1="<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 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 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 @@ -2867,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; @@ -3142,6 +3413,31 @@ bool LogicalOR::IsSatisfy( long theId ) FILTER */ +// #ifdef WITH_TBB +// #include +// #include + +// 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& r ) const +// { +// for ( size_t i = r.begin(); i != r.end(); ++i ) +// if ( myPredicate->IsSatisfy( i )) +// myOKIds.local().push_back(); +// } +// } +// } +// #endif + Filter::Filter() {} @@ -3671,23 +3967,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,167 +3994,608 @@ 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) +{ + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified()) + { + size_t nbNodes = theMesh ? theMesh->NbNodes() : 0; + if ( myNodeIsChecked.size() == nbNodes ) + { + std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false ); + } + else + { + SMESHUtils::FreeVector( myNodeIsChecked ); + SMESHUtils::FreeVector( myNodeIsOut ); + myNodeIsChecked.resize( nbNodes, false ); + myNodeIsOut.resize( nbNodes ); + } + } +} + +bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut ) +{ + if ( n->GetID() >= (int) myNodeIsChecked.size() || + !myNodeIsChecked[ n->GetID() ]) + return false; + + isOut = myNodeIsOut[ n->GetID() ]; + return true; +} + +void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool isOut ) +{ + if ( n->GetID() < (int) myNodeIsChecked.size() ) + { + myNodeIsChecked[ n->GetID() ] = true; + myNodeIsOut [ n->GetID() ] = isOut; } } void ElementsOnShape::SetShape (const TopoDS_Shape& theShape, const SMDSAbs_ElementType theType) { - myType = theType; + myType = theType; myShape = theShape; - myIds.Clear(); + if ( myShape.IsNull() ) return; - const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh(); - - if ( !myMesh ) return; + 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() ); + } - myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType )); + clearClassifiers(); + myClassifiers.resize( shapesMap.Extent() ); + for ( int i = 0; i < shapesMap.Extent(); ++i ) + myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler ); + + if ( theType == SMDSAbs_Node ) + { + SMESHUtils::FreeVector( myNodeIsChecked ); + SMESHUtils::FreeVector( myNodeIsOut ); + } + else + { + std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false ); + } +} - myShapesMap.Clear(); - addShape(myShape); +void ElementsOnShape::clearClassifiers() +{ + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + delete myClassifiers[ i ]; + myClassifiers.clear(); } -void ElementsOnShape::addShape (const TopoDS_Shape& theShape) +bool ElementsOnShape::IsSatisfy (long elemId) { - if (theShape.IsNull() || myMeshModifTracer.GetMesh() == 0) - return; + const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh(); + const SMDS_MeshElement* elem = + ( myType == SMDSAbs_Node ? mesh->FindNode( elemId ) : mesh->FindElement( elemId )); + if ( !elem || myClassifiers.empty() ) + return false; - if (!myShapesMap.Add(theShape)) return; + bool isSatisfy = myAllNodesFlag, isNodeOut; + + gp_XYZ centerXYZ (0, 0, 0); - myCurShapeType = theShape.ShapeType(); - switch (myCurShapeType) + SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator(); + while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) { - case TopAbs_COMPOUND: - case TopAbs_COMPSOLID: - case TopAbs_SHELL: - case TopAbs_WIRE: + SMESH_TNodeXYZ aPnt( aNodeItr->next() ); + centerXYZ += aPnt; + + isNodeOut = true; + if ( !getNodeIsOut( aPnt._node, isNodeOut )) { - TopoDS_Iterator anIt (theShape, Standard_True, Standard_True); - for (; anIt.More(); anIt.Next()) addShape(anIt.Value()); + for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i ) + isNodeOut = myClassifiers[i]->IsOut( aPnt ); + + setNodeIsOut( aPnt._node, isNodeOut ); } - break; - case TopAbs_SOLID: + isSatisfy = !isNodeOut; + } + + // Check the center point for volumes MantisBug 0020168 + if (isSatisfy && + myAllNodesFlag && + myClassifiers[0]->ShapeType() == TopAbs_SOLID) + { + centerXYZ /= elem->NbNodes(); + isSatisfy = false; + for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i ) + isSatisfy = ! myClassifiers[i]->IsOut( centerXYZ ); + } + + return isSatisfy; +} + +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: { + if ( isBox( theShape )) { - myCurSC.Load(theShape); - process(); + myIsOutFun = & ElementsOnShape::TClassifier::isOutOfBox; } - 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(); + else + { + mySolidClfr.Load(theShape); + myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid; } break; - case TopAbs_EDGE: + } + 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: { + 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:{ + myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) ); + myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex; + break; + } + default: + throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier"); + } +} + +bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p) +{ + mySolidClfr.Perform( p, myTol ); + 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 ); + if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol ) + { + // 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; +} + +bool ElementsOnShape::TClassifier::isOutOfEdge (const gp_Pnt& p) +{ + myProjEdge.Perform( p ); + return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol ); +} + +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 ) { - 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(); + 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; } - break; - case TopAbs_VERTEX: + } + 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(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()) { - TopoDS_Vertex aV = TopoDS::Vertex(theShape); - myCurPnt = BRep_Tool::Pnt(aV); - process(); + if (!IsSubShape(theMap, anIt.Value())) { + return false; + } } - break; - default: - break; + return true; } + + return false; } -void ElementsOnShape::process() +void BelongToGeom::init() { - const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh(); - if (myShape.IsNull() || myMesh == 0) - return; + 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); + } - SMDS_ElemIteratorPtr anIter = myMesh->elementsIterator(myType); - while (anIter->more()) - process(anIter->next()); + //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); + } } -void ElementsOnShape::process (const SMDS_MeshElement* theElemPtr) +static bool IsContains( const SMESHDS_Mesh* theMeshDS, + const TopoDS_Shape& theShape, + const SMDS_MeshElement* theElem, + TopAbs_ShapeEnum theFindShapeEnum, + TopAbs_ShapeEnum theAvoidShapeEnum = TopAbs_SHAPE ) { - if (myShape.IsNull()) - return; + TopExp_Explorer anExp( theShape,theFindShapeEnum,theAvoidShapeEnum ); - SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator(); - bool isSatisfy = myAllNodesFlag; + 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; +} - gp_XYZ centerXYZ (0, 0, 0); +bool BelongToGeom::IsSatisfy (long theId) +{ + if (myMeshDS == 0 || myShape.IsNull()) + return false; - while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) + if (!myIsSubshape) { - SMESH_TNodeXYZ aPnt ( aNodeItr->next() ); - centerXYZ += aPnt; + return myElementsOnShapePtr->IsSatisfy(theId); + } - switch (myCurShapeType) + // Case of submesh + if (myType == SMDSAbs_Node) + { + if( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ) ) { - case TopAbs_SOLID: + if ( aNode->getshapeId() < 1 ) + return myElementsOnShapePtr->IsSatisfy(theId); + + const SMDS_PositionPtr& aPosition = aNode->GetPosition(); + SMDS_TypeOfPosition aTypeOfPosition = aPosition->GetTypeOfPosition(); + switch( aTypeOfPosition ) { - myCurSC.Perform(aPnt, myToler); - isSatisfy = (myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON); + 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 )); } - break; - case TopAbs_FACE: + } + } + 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() ) { - myCurProjFace.Perform(aPnt); - isSatisfy = (myCurProjFace.IsDone() && myCurProjFace.LowerDistance() <= myToler); - if (isSatisfy) + switch( myType ) { - // 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); + 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 )); } } - break; - case TopAbs_EDGE: + } + } + + 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(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 ) { - myCurProjEdge.Perform(aPnt); - isSatisfy = (myCurProjEdge.NbPoints() > 0 && myCurProjEdge.LowerDistance() <= myToler); + 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 ); } - break; - case TopAbs_VERTEX: + } + } + else + { + if( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ) ) + { + if( myType == SMDSAbs_All ) { - isSatisfy = (myCurPnt.Distance(aPnt) <= myToler); + 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 ); } - break; - default: + else if( myType == anElem->GetType() ) { - isSatisfy = false; + 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 ); + } } } } - 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; - } + return false; +} - if (isSatisfy) - myIds.Add(theElemPtr->GetID()); +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(aNodeIt->next()); + SMDS_ElemIteratorPtr anElemIt = aNode->GetInverseElementIterator(); + while ( anElemIt->more() ) + { + const SMDS_MeshElement* anElement = static_cast(anElemIt->next()); + if (anElement == theElem) + return true; + } + } + } + } + return false; } TSequenceOfXYZ::TSequenceOfXYZ()