X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=3dd59cf5dadbdf63cc4afdab916c7009b4387666;hb=cca53b9163eee6588f78ba14dc33d3d36b60ad5a;hp=06bfd5304a9fb4422a87b325dee2880818a9cc41;hpb=c3bf92bd87b770fd81631a3853f7f5bb1ac6a4e8;p=modules%2Fsmesh.git diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 06bfd5304..3dd59cf5d 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -17,90 +17,110 @@ // // See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org +#include "SMESH_ControlsDef.hxx" + #include +#include +#include +#include +#include #include +#include #include #include +#include +#include #include -#include +#include #include +#include +#include +#include +#include +#include +#include #include "SMDS_Mesh.hxx" #include "SMDS_Iterator.hxx" #include "SMDS_MeshElement.hxx" #include "SMDS_MeshNode.hxx" +#include "SMDS_VolumeTool.hxx" -#include "SMESH_Controls.hxx" /* AUXILIARY METHODS */ -static inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) -{ - return gp_Vec( P1 - P2 ).Angle( gp_Vec( P3 - P2 ) ); -} - -static inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) -{ - gp_Vec aVec1( P2 - P1 ); - gp_Vec aVec2( P3 - P1 ); - return ( aVec1 ^ aVec2 ).Magnitude() * 0.5; -} - -static inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 ) -{ - return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() ); -} +namespace{ + inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) + { + gp_Vec v1( P1 - P2 ), v2( P3 - P2 ); + + return v1.Magnitude() < gp::Resolution() || + v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 ); + } -static inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 ) -{ - double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) ); - return aDist; -} + inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) + { + gp_Vec aVec1( P2 - P1 ); + gp_Vec aVec2( P3 - P1 ); + return ( aVec1 ^ aVec2 ).Magnitude() * 0.5; + } -static int getNbMultiConnection( SMDS_Mesh* theMesh, const int theId ) -{ - if ( theMesh == 0 ) - return 0; + inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 ) + { + return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() ); + } - const SMDS_MeshElement* anEdge = theMesh->FindElement( theId ); - if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 ) - return 0; - TColStd_MapOfInteger aMap; - int aResult = 0; - SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator(); - if ( anIter != 0 ) + inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 ) { - while( anIter->more() ) - { - const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); - if ( aNode == 0 ) - return 0; - SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); - while( anElemIter->more() ) - { - const SMDS_MeshElement* anElem = anElemIter->next(); - if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) - { - int anId = anElem->GetID(); + double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) ); + return aDist; + } - if ( anIter->more() ) // i.e. first node - aMap.Add( anId ); - else if ( aMap.Contains( anId ) ) - aResult++; - } + int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId ) + { + if ( theMesh == 0 ) + return 0; + + const SMDS_MeshElement* anEdge = theMesh->FindElement( theId ); + if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 ) + return 0; + + TColStd_MapOfInteger aMap; + + int aResult = 0; + SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator(); + if ( anIter != 0 ) { + while( anIter->more() ) { + const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); + if ( aNode == 0 ) + return 0; + SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); + while( anElemIter->more() ) { + const SMDS_MeshElement* anElem = anElemIter->next(); + if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) { + int anId = anElem->GetID(); + + if ( anIter->more() ) // i.e. first node + aMap.Add( anId ); + else if ( aMap.Contains( anId ) ) + aResult++; + } + } } } + + return aResult; } - return aResult; } + using namespace SMESH::Controls; /* @@ -113,76 +133,106 @@ using namespace SMESH::Controls; */ NumericalFunctor::NumericalFunctor(): myMesh(NULL) -{} +{ + myPrecision = -1; +} -void NumericalFunctor::SetMesh( SMDS_Mesh* theMesh ) +void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh ) { myMesh = theMesh; } -bool NumericalFunctor::getPoints( const int theId, - TColgp_SequenceOfXYZ& theRes ) const +bool NumericalFunctor::GetPoints(const int theId, + TSequenceOfXYZ& theRes ) const { - theRes.Clear(); + theRes.clear(); if ( myMesh == 0 ) return false; - // Get nodes of the face - const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); - if ( anElem == 0 || anElem->GetType() != GetType() ) - return false; + return GetPoints( myMesh->FindElement( theId ), theRes ); +} + +bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, + TSequenceOfXYZ& theRes ) +{ + theRes.clear(); - int nbNodes = anElem->NbNodes(); + if ( anElem == 0) + return false; + // Get nodes of the element SMDS_ElemIteratorPtr anIter = anElem->nodesIterator(); if ( anIter != 0 ) { while( anIter->more() ) { const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); - if ( aNode != 0 ) - theRes.Append( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); + if ( aNode != 0 ){ + theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); + } } } return true; } +long NumericalFunctor::GetPrecision() const +{ + return myPrecision; +} + +void NumericalFunctor::SetPrecision( const long thePrecision ) +{ + myPrecision = thePrecision; +} + +double NumericalFunctor::GetValue( long theId ) +{ + TSequenceOfXYZ P; + if ( GetPoints( theId, P )) + { + double aVal = GetValue( P ); + if ( myPrecision >= 0 ) + { + double prec = pow( 10., (double)( myPrecision ) ); + aVal = floor( aVal * prec + 0.5 ) / prec; + } + return aVal; + } + + return 0.; +} /* Class : MinimumAngle Description : Functor for calculation of minimum angle */ -double MinimumAngle::GetValue( long theId ) -{ - TColgp_SequenceOfXYZ P; - if ( !getPoints( theId, P ) || P.Length() != 3 && P.Length() != 4 ) - return 0; +double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) +{ double aMin; - if ( P.Length() == 3 ) - { - double A0 = getAngle( P( 3 ), P( 1 ), P( 2 ) ); - double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) ); - double A2 = getAngle( P( 2 ), P( 3 ), P( 1 ) ); + if (P.size() <3) + return 0.; - aMin = Min( A0, Min( A1, A2 ) ); - } - else - { - double A0 = getAngle( P( 4 ), P( 1 ), P( 2 ) ); - double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) ); - double A2 = getAngle( P( 2 ), P( 3 ), P( 4 ) ); - double A3 = getAngle( P( 3 ), P( 4 ), P( 1 ) ); - - aMin = Min( Min( A0, A1 ), Min( A2, A3 ) ); - } + aMin = getAngle(P( P.size() ), P( 1 ), P( 2 )); + aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 ))); + for (int i=2; i 3) + aArea = getArea( P( 1 ), P( 2 ), P( 3 ) ); else - return getArea( P( 1 ), P( 2 ), P( 3 ) ) + getArea( P( 1 ), P( 3 ), P( 4 ) ); + return 0; + + for (int i=4; i<=P.size(); i++) + aArea += getArea(P(1),P(i-1),P(i)); + return aArea; +} + +double Area::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; } SMDSAbs_ElementType Area::GetType() const @@ -383,51 +790,421 @@ SMDSAbs_ElementType Area::GetType() const Class : Length Description : Functor for calculating length off edge */ -double Length::GetValue( long theId ) -{ - TColgp_SequenceOfXYZ P; - return getPoints( theId, P ) && P.Length() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0; -} - -SMDSAbs_ElementType Length::GetType() const +double Length::GetValue( const TSequenceOfXYZ& P ) { - return SMDSAbs_Edge; + return ( P.size() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0 ); } - -/* - Class : MultiConnection - Description : Functor for calculating number of faces conneted to the edge -*/ -double MultiConnection::GetValue( long theId ) +double Length::GetBadRate( double Value, int /*nbNodes*/ ) const { - return getNbMultiConnection( myMesh, theId ); + return Value; } -SMDSAbs_ElementType MultiConnection::GetType() const +SMDSAbs_ElementType Length::GetType() const { return SMDSAbs_Edge; } - -/* - PREDICATES -*/ - /* - Class : FreeBorders - Description : Predicate for free borders + Class : Length2D + Description : Functor for calculating length of edge */ -FreeBorders::FreeBorders() +double Length2D::GetValue( long theElementId) { - myMesh = 0; -} + TSequenceOfXYZ P; -void FreeBorders::SetMesh( SMDS_Mesh* theMesh ) -{ - myMesh = theMesh; -} + if (GetPoints(theElementId,P)){ + + double aVal;// = GetValue( P ); + const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId ); + SMDSAbs_ElementType aType = aElem->GetType(); + + int len = P.size(); + + switch (aType){ + case SMDSAbs_All: + case SMDSAbs_Node: + case SMDSAbs_Edge: + if (len == 2){ + aVal = getDistance( P( 1 ), P( 2 ) ); + break; + } + case SMDSAbs_Face: + if (len == 3){ // triangles + 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)); + break; + } + else if (len == 4){ // quadrangles + double L1 = getDistance(P( 1 ),P( 2 )); + 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)); + break; + } + case SMDSAbs_Volume: + if (len == 4){ // tetraidrs + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + double L4 = getDistance(P( 1 ),P( 4 )); + double L5 = getDistance(P( 2 ),P( 4 )); + double L6 = getDistance(P( 3 ),P( 4 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + break; + } + else if (len == 5){ // piramids + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 1 )); + double L5 = getDistance(P( 1 ),P( 5 )); + double L6 = getDistance(P( 2 ),P( 5 )); + double L7 = getDistance(P( 3 ),P( 5 )); + double L8 = getDistance(P( 4 ),P( 5 )); + + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(L7,L8)); + break; + } + else if (len == 6){ // pentaidres + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 4 )); + double L7 = getDistance(P( 1 ),P( 4 )); + double L8 = getDistance(P( 2 ),P( 5 )); + double L9 = getDistance(P( 3 ),P( 6 )); + + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),L9)); + break; + } + else if (len == 8){ // hexaider + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 7 )); + double L7 = getDistance(P( 7 ),P( 8 )); + double L8 = getDistance(P( 8 ),P( 5 )); + double L9 = getDistance(P( 1 ),P( 5 )); + double L10= getDistance(P( 2 ),P( 6 )); + double L11= getDistance(P( 3 ),P( 7 )); + double L12= getDistance(P( 4 ),P( 8 )); + + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10))); + aVal = Max(aVal,Max(L11,L12)); + break; + + } + + default: aVal=-1; + } + + if (aVal <0){ + return 0.; + } + + if ( myPrecision >= 0 ) + { + double prec = pow( 10., (double)( myPrecision ) ); + aVal = floor( aVal * prec + 0.5 ) / prec; + } + + return aVal; + + } + return 0.; +} + +double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; +} + +SMDSAbs_ElementType Length2D::GetType() const +{ + return SMDSAbs_Face; +} + +Length2D::Value::Value(double theLength,long thePntId1, long thePntId2): + myLength(theLength) +{ + myPntId[0] = thePntId1; myPntId[1] = thePntId2; + if(thePntId1 > thePntId2){ + myPntId[1] = thePntId1; myPntId[0] = thePntId2; + } +} + +bool Length2D::Value::operator<(const Length2D::Value& x) const{ + if(myPntId[0] < x.myPntId[0]) return true; + if(myPntId[0] == x.myPntId[0]) + if(myPntId[1] < x.myPntId[1]) return true; + return false; +} + +void Length2D::GetValues(TValues& theValues){ + TValues aValues; + SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); + for(; anIter->more(); ){ + const SMDS_MeshFace* anElem = anIter->next(); + SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator(); + long aNodeId[2]; + gp_Pnt P[3]; + + double aLength; + const SMDS_MeshElement* aNode; + if(aNodesIter->more()){ + aNode = aNodesIter->next(); + const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode; + P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z()); + aNodeId[0] = aNodeId[1] = aNode->GetID(); + aLength = 0; + } + for(; aNodesIter->more(); ){ + aNode = aNodesIter->next(); + const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode; + long anId = aNode->GetID(); + + P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z()); + + aLength = P[1].Distance(P[2]); + + Value aValue(aLength,aNodeId[1],anId); + aNodeId[1] = anId; + P[1] = P[2]; + theValues.insert(aValue); + } + + aLength = P[0].Distance(P[1]); + + Value aValue(aLength,aNodeId[0],aNodeId[1]); + theValues.insert(aValue); + } +} + +/* + Class : MultiConnection + Description : Functor for calculating number of faces conneted to the edge +*/ +double MultiConnection::GetValue( const TSequenceOfXYZ& P ) +{ + return 0; +} +double MultiConnection::GetValue( long theId ) +{ + return getNbMultiConnection( myMesh, theId ); +} + +double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; +} + +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; +} + +double MultiConnection2D::GetValue( long theElementId ) +{ + TSequenceOfXYZ P; + int aResult = 0; + + if (GetPoints(theElementId,P)){ + const SMDS_MeshElement* anFaceElem = myMesh->FindElement( theElementId ); + SMDSAbs_ElementType aType = anFaceElem->GetType(); + + int len = P.size(); + + TColStd_MapOfInteger aMap; + int aResult = 0; + + switch (aType){ + case SMDSAbs_All: + case SMDSAbs_Node: + case SMDSAbs_Edge: + case SMDSAbs_Face: + if (len == 3){ // triangles + int Nb[3] = {0,0,0}; + + int i=0; + SMDS_ElemIteratorPtr anIter = anFaceElem->nodesIterator(); + if ( anIter != 0 ) { + while( anIter->more() ) { + const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); + if ( aNode == 0 ){ + break; + } + SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); + while( anElemIter->more() ) { + const SMDS_MeshElement* anElem = anElemIter->next(); + if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) { + int anId = anElem->GetID(); + + if ( anIter->more() ) // i.e. first node + aMap.Add( anId ); + else if ( aMap.Contains( anId ) ){ + Nb[i]++; + } + } + else if ( anElem != 0 && anElem->GetType() == SMDSAbs_Edge ) i++; + } + } + } + + aResult = Max(Max(Nb[0],Nb[1]),Nb[2]); + } + break; + case SMDSAbs_Volume: + default: aResult=0; + } + + } + return aResult;//getNbMultiConnection( myMesh, theId ); +} + +double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; +} + +SMDSAbs_ElementType MultiConnection2D::GetType() const +{ + return SMDSAbs_Face; +} + +MultiConnection2D::Value::Value(long thePntId1, long thePntId2) +{ + myPntId[0] = thePntId1; myPntId[1] = thePntId2; + if(thePntId1 > thePntId2){ + myPntId[1] = thePntId1; myPntId[0] = thePntId2; + } +} + +bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{ + if(myPntId[0] < x.myPntId[0]) return true; + if(myPntId[0] == x.myPntId[0]) + if(myPntId[1] < x.myPntId[1]) return true; + return false; +} + +void MultiConnection2D::GetValues(MValues& theValues){ + SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); + for(; anIter->more(); ){ + const SMDS_MeshFace* anElem = anIter->next(); + SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator(); + long aNodeId[3]; + + //int aNbConnects=0; + const SMDS_MeshNode* aNode0; + const SMDS_MeshNode* aNode1; + const SMDS_MeshNode* aNode2; + if(aNodesIter->more()){ + aNode0 = (SMDS_MeshNode*) aNodesIter->next(); + aNode1 = aNode0; + const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1; + aNodeId[0] = aNodeId[1] = aNodes->GetID(); + } + for(; aNodesIter->more(); ){ + aNode2 = (SMDS_MeshNode*) aNodesIter->next(); + long anId = aNode2->GetID(); + aNodeId[2] = anId; + + Value aValue(aNodeId[1],aNodeId[2]); + MValues::iterator aItr = theValues.find(aValue); + if (aItr != theValues.end()){ + aItr->second += 1; + //aNbConnects = nb; + } else { + theValues[aValue] = 1; + //aNbConnects = 1; + } + //cout << "NodeIds: "<GetID(); + Border aBorder(anElemId,aNodeId[1],anId); + aNodeId[1] = anId; + //std::cout< anIntSeq.Last() ) + { + anIntSeq.Append( aMinId ); + aStrSeq.Append( aStr ); + } + else + for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ ) + if ( aMinId < anIntSeq( j ) ) + { + anIntSeq.InsertBefore( j, aMinId ); + aStrSeq.InsertBefore( j, aStr ); + break; + } + } + } + + if ( aStrSeq.Length() == 0 ) + return; + + theResStr = aStrSeq( 1 ); + for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ ) + { + theResStr += ","; + theResStr += aStrSeq( j ); } - //std::cout<<"aContainer.size() = "<FindNode( theId ) == 0 ) + return false; + } + else + { + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All ) + return false; + } + + if ( myIds.Contains( theId ) ) + return true; + + for ( int i = 1, n = myMin.Length(); i <= n; i++ ) + if ( theId >= myMin( i ) && theId <= myMax( i ) ) + return true; + + return false; +} /* Class : Comparator @@ -560,7 +1575,7 @@ Comparator::Comparator(): Comparator::~Comparator() {} -void Comparator::SetMesh( SMDS_Mesh* theMesh ) +void Comparator::SetMesh( const SMDS_Mesh* theMesh ) { if ( myFunctor ) myFunctor->SetMesh( theMesh ); @@ -581,6 +1596,11 @@ SMDSAbs_ElementType Comparator::GetType() const return myFunctor ? myFunctor->GetType() : SMDSAbs_All; } +double Comparator::GetMargin() +{ + return myMargin; +} + /* Class : LessThan @@ -620,6 +1640,10 @@ void EqualTo::SetTolerance( double theToler ) myToler = theToler; } +double EqualTo::GetTolerance() +{ + return myToler; +} /* Class : LogicalNOT @@ -636,7 +1660,7 @@ bool LogicalNOT::IsSatisfy( long theId ) return myPredicate && !myPredicate->IsSatisfy( theId ); } -void LogicalNOT::SetMesh( SMDS_Mesh* theMesh ) +void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh ) { if ( myPredicate ) myPredicate->SetMesh( theMesh ); @@ -663,7 +1687,7 @@ LogicalBinary::LogicalBinary() LogicalBinary::~LogicalBinary() {} -void LogicalBinary::SetMesh( SMDS_Mesh* theMesh ) +void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh ) { if ( myPredicate1 ) myPredicate1->SetMesh( theMesh ); @@ -737,38 +1761,566 @@ void Filter::SetPredicate( PredicatePtr thePredicate ) myPredicate = thePredicate; } -Filter::TIdSequence -Filter::GetElementsId( SMDS_Mesh* theMesh ) +template +inline void FillSequence(const TIterator& theIterator, + TPredicate& thePredicate, + Filter::TIdSequence& theSequence) { - TIdSequence aSequence; - if ( !theMesh || !myPredicate ) return aSequence; + if ( theIterator ) { + while( theIterator->more() ) { + TElement anElem = theIterator->next(); + long anId = anElem->GetID(); + if ( thePredicate->IsSatisfy( anId ) ) + theSequence.push_back( anId ); + } + } +} - myPredicate->SetMesh( theMesh ); +void +Filter:: +GetElementsId( const SMDS_Mesh* theMesh, + PredicatePtr thePredicate, + TIdSequence& theSequence ) +{ + theSequence.clear(); - SMDSAbs_ElementType aType = myPredicate->GetType(); + if ( !theMesh || !thePredicate ) + return; + + thePredicate->SetMesh( theMesh ); + + SMDSAbs_ElementType aType = thePredicate->GetType(); switch(aType){ - case SMDSAbs_Edge:{ - SMDS_EdgeIteratorPtr anIter = theMesh->edgesIterator(); - if ( anIter != 0 ) { - while( anIter->more() ) { - const SMDS_MeshElement* anElem = anIter->next(); - long anId = anElem->GetID(); - if ( myPredicate->IsSatisfy( anId ) ) - aSequence.push_back( anId ); + case SMDSAbs_Node: + FillSequence(theMesh->nodesIterator(),thePredicate,theSequence); + break; + case SMDSAbs_Edge: + FillSequence(theMesh->edgesIterator(),thePredicate,theSequence); + break; + case SMDSAbs_Face: + FillSequence(theMesh->facesIterator(),thePredicate,theSequence); + break; + case SMDSAbs_Volume: + FillSequence(theMesh->volumesIterator(),thePredicate,theSequence); + break; + case SMDSAbs_All: + FillSequence(theMesh->edgesIterator(),thePredicate,theSequence); + FillSequence(theMesh->facesIterator(),thePredicate,theSequence); + FillSequence(theMesh->volumesIterator(),thePredicate,theSequence); + break; + } +} + +void +Filter::GetElementsId( const SMDS_Mesh* theMesh, + Filter::TIdSequence& theSequence ) +{ + GetElementsId(theMesh,myPredicate,theSequence); +} + +/* + ManifoldPart +*/ + +typedef std::set TMapOfFacePtr; + +/* + Internal class Link +*/ + +ManifoldPart::Link::Link( SMDS_MeshNode* theNode1, + SMDS_MeshNode* theNode2 ) +{ + myNode1 = theNode1; + myNode2 = theNode2; +} + +ManifoldPart::Link::~Link() +{ + myNode1 = 0; + myNode2 = 0; +} + +bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const +{ + if ( myNode1 == theLink.myNode1 && + myNode2 == theLink.myNode2 ) + return true; + else if ( myNode1 == theLink.myNode2 && + myNode2 == theLink.myNode1 ) + return true; + else + return false; +} + +bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const +{ + if(myNode1 < x.myNode1) return true; + if(myNode1 == x.myNode1) + if(myNode2 < x.myNode2) return true; + return false; +} + +bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1, + const ManifoldPart::Link& theLink2 ) +{ + return theLink1.IsEqual( theLink2 ); +} + +ManifoldPart::ManifoldPart() +{ + myMesh = 0; + myAngToler = Precision::Angular(); + myIsOnlyManifold = true; +} + +ManifoldPart::~ManifoldPart() +{ + myMesh = 0; +} + +void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; + process(); +} + +SMDSAbs_ElementType ManifoldPart::GetType() const +{ return SMDSAbs_Face; } + +bool ManifoldPart::IsSatisfy( long theElementId ) +{ + return myMapIds.Contains( theElementId ); +} + +void ManifoldPart::SetAngleTolerance( const double theAngToler ) +{ myAngToler = theAngToler; } + +double ManifoldPart::GetAngleTolerance() const +{ return myAngToler; } + +void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly ) +{ myIsOnlyManifold = theIsOnly; } + +void ManifoldPart::SetStartElem( const long theStartId ) +{ myStartElemId = theStartId; } + +bool ManifoldPart::process() +{ + myMapIds.Clear(); + myMapBadGeomIds.Clear(); + + myAllFacePtr.clear(); + myAllFacePtrIntDMap.clear(); + if ( !myMesh ) + return false; + + // collect all faces into own map + SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator(); + for (; anFaceItr->more(); ) + { + SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next(); + myAllFacePtr.push_back( aFacePtr ); + myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1; + } + + SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId ); + if ( !aStartFace ) + return false; + + // the map of non manifold links and bad geometry + TMapOfLink aMapOfNonManifold; + TColStd_MapOfInteger aMapOfTreated; + + // begin cycle on faces from start index and run on vector till the end + // and from begin to start index to cover whole vector + const int aStartIndx = myAllFacePtrIntDMap[aStartFace]; + bool isStartTreat = false; + for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ ) + { + if ( fi == aStartIndx ) + isStartTreat = true; + // as result next time when fi will be equal to aStartIndx + + SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ]; + if ( aMapOfTreated.Contains( aFacePtr->GetID() ) ) + continue; + + aMapOfTreated.Add( aFacePtr->GetID() ); + TColStd_MapOfInteger aResFaces; + if ( !findConnected( myAllFacePtrIntDMap, aFacePtr, + aMapOfNonManifold, aResFaces ) ) + continue; + TColStd_MapIteratorOfMapOfInteger anItr( aResFaces ); + for ( ; anItr.More(); anItr.Next() ) + { + int aFaceId = anItr.Key(); + aMapOfTreated.Add( aFaceId ); + myMapIds.Add( aFaceId ); + } + + if ( fi == ( myAllFacePtr.size() - 1 ) ) + fi = 0; + } // end run on vector of faces + return !myMapIds.IsEmpty(); +} + +static void getLinks( const SMDS_MeshFace* theFace, + ManifoldPart::TVectorOfLink& theLinks ) +{ + int aNbNode = theFace->NbNodes(); + SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator(); + int i = 1; + SMDS_MeshNode* aNode = 0; + for ( ; aNodeItr->more() && i <= aNbNode; ) + { + + SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next(); + if ( i == 1 ) + aNode = aN1; + i++; + SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next(); + i++; + ManifoldPart::Link aLink( aN1, aN2 ); + theLinks.push_back( aLink ); + } +} + +static gp_XYZ getNormale( const SMDS_MeshFace* theFace ) +{ + gp_XYZ n; + int aNbNode = theFace->NbNodes(); + TColgp_Array1OfXYZ anArrOfXYZ(1,4); + SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator(); + int i = 1; + for ( ; aNodeItr->more() && i <= 4; i++ ) + { + SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); + anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); + } + + gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1); + gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1); + n = q1 ^ q2; + if ( aNbNode > 3 ) + { + gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1); + n += q2 ^ q3; + } + double len = n.Modulus(); + if ( len > 0 ) + n /= len; + + return n; +} + +bool ManifoldPart::findConnected + ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt, + SMDS_MeshFace* theStartFace, + ManifoldPart::TMapOfLink& theNonManifold, + TColStd_MapOfInteger& theResFaces ) +{ + theResFaces.Clear(); + if ( !theAllFacePtrInt.size() ) + return false; + + if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() ) + { + myMapBadGeomIds.Add( theStartFace->GetID() ); + return false; + } + + ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip; + ManifoldPart::TVectorOfLink aSeqOfBoundary; + theResFaces.Add( theStartFace->GetID() ); + ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace; + + expandBoundary( aMapOfBoundary, aSeqOfBoundary, + aDMapLinkFace, theNonManifold, theStartFace ); + + bool isDone = false; + while ( !isDone && aMapOfBoundary.size() != 0 ) + { + bool isToReset = false; + ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin(); + for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink ) + { + ManifoldPart::Link aLink = *pLink; + if ( aMapToSkip.find( aLink ) != aMapToSkip.end() ) + continue; + // each link could be treated only once + aMapToSkip.insert( aLink ); + + ManifoldPart::TVectorOfFacePtr aFaces; + // find next + if ( myIsOnlyManifold && + (theNonManifold.find( aLink ) != theNonManifold.end()) ) + continue; + else + { + getFacesByLink( aLink, aFaces ); + // filter the element to keep only indicated elements + ManifoldPart::TVectorOfFacePtr aFiltered; + ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin(); + for ( ; pFace != aFaces.end(); ++pFace ) + { + SMDS_MeshFace* aFace = *pFace; + if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() ) + aFiltered.push_back( aFace ); + } + aFaces = aFiltered; + if ( aFaces.size() < 2 ) // no neihgbour faces + continue; + else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case + { + theNonManifold.insert( aLink ); + continue; + } + } + + // compare normal with normals of neighbor element + SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ]; + ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin(); + for ( ; pFace != aFaces.end(); ++pFace ) + { + SMDS_MeshFace* aNextFace = *pFace; + if ( aPrevFace == aNextFace ) + continue; + int anNextFaceID = aNextFace->GetID(); + if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) ) + // should not be with non manifold restriction. probably bad topology + continue; + // check if face was treated and skipped + if ( myMapBadGeomIds.Contains( anNextFaceID ) || + !isInPlane( aPrevFace, aNextFace ) ) + continue; + // add new element to connected and extend the boundaries. + theResFaces.Add( anNextFaceID ); + expandBoundary( aMapOfBoundary, aSeqOfBoundary, + aDMapLinkFace, theNonManifold, aNextFace ); + isToReset = true; } } + isDone = !isToReset; } - case SMDSAbs_Face:{ - SMDS_FaceIteratorPtr anIter = theMesh->facesIterator(); - if ( anIter != 0 ) { - while( anIter->more() ) { - const SMDS_MeshElement* anElem = anIter->next(); - long anId = anElem->GetID(); - if ( myPredicate->IsSatisfy( anId ) ) - aSequence.push_back( anId ); + + return !theResFaces.IsEmpty(); +} + +bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1, + const SMDS_MeshFace* theFace2 ) +{ + gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) ); + gp_XYZ aNorm2XYZ = getNormale( theFace2 ); + if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() ) + { + myMapBadGeomIds.Add( theFace2->GetID() ); + return false; + } + if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) ) + return true; + + return false; +} + +void ManifoldPart::expandBoundary + ( ManifoldPart::TMapOfLink& theMapOfBoundary, + ManifoldPart::TVectorOfLink& theSeqOfBoundary, + ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr, + ManifoldPart::TMapOfLink& theNonManifold, + SMDS_MeshFace* theNextFace ) const +{ + ManifoldPart::TVectorOfLink aLinks; + getLinks( theNextFace, aLinks ); + int aNbLink = aLinks.size(); + for ( int i = 0; i < aNbLink; i++ ) + { + ManifoldPart::Link aLink = aLinks[ i ]; + if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) ) + continue; + if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() ) + { + if ( myIsOnlyManifold ) + { + // remove from boundary + theMapOfBoundary.erase( aLink ); + ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin(); + for ( ; pLink != theSeqOfBoundary.end(); ++pLink ) + { + ManifoldPart::Link aBoundLink = *pLink; + if ( aBoundLink.IsEqual( aLink ) ) + { + theSeqOfBoundary.erase( pLink ); + break; + } + } } } + else + { + theMapOfBoundary.insert( aLink ); + theSeqOfBoundary.push_back( aLink ); + theDMapLinkFacePtr[ aLink ] = theNextFace; + } + } +} + +void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink, + ManifoldPart::TVectorOfFacePtr& theFaces ) const +{ + SMDS_Mesh::SetOfFaces aSetOfFaces; + // take all faces that shared first node + SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator(); + for ( ; anItr->more(); ) + { + SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next(); + if ( !aFace ) + continue; + aSetOfFaces.Add( aFace ); + } + // take all faces that shared second node + anItr = theLink.myNode2->facesIterator(); + // find the common part of two sets + for ( ; anItr->more(); ) + { + SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next(); + if ( aSetOfFaces.Contains( aFace ) ) + theFaces.push_back( aFace ); + } +} + + +/* + ElementsOnSurface +*/ + +ElementsOnSurface::ElementsOnSurface() +{ + myMesh = 0; + myIds.Clear(); + myType = SMDSAbs_All; + mySurf.Nullify(); + myToler = Precision::Confusion(); +} + +ElementsOnSurface::~ElementsOnSurface() +{ + myMesh = 0; +} + +void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh ) +{ + if ( myMesh == theMesh ) + return; + myMesh = theMesh; + myIds.Clear(); + process(); +} + +bool ElementsOnSurface::IsSatisfy( long theElementId ) +{ + return myIds.Contains( theElementId ); +} + +SMDSAbs_ElementType ElementsOnSurface::GetType() const +{ return myType; } + +void ElementsOnSurface::SetTolerance( const double theToler ) +{ myToler = theToler; } + +double ElementsOnSurface::GetTolerance() const +{ + return myToler; +} + +void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape, + const SMDSAbs_ElementType theType ) +{ + myType = theType; + mySurf.Nullify(); + if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE ) + { + mySurf.Nullify(); + return; + } + TopoDS_Face aFace = TopoDS::Face( theShape ); + mySurf = BRep_Tool::Surface( aFace ); +} + +void ElementsOnSurface::process() +{ + myIds.Clear(); + if ( mySurf.IsNull() ) + return; + + if ( myMesh == 0 ) + return; + + if ( myType == SMDSAbs_Face || myType == SMDSAbs_All ) + { + SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); + for(; anIter->more(); ) + process( anIter->next() ); + } + + if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All ) + { + SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator(); + for(; anIter->more(); ) + process( anIter->next() ); + } + + if ( myType == SMDSAbs_Node ) + { + SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator(); + for(; anIter->more(); ) + process( anIter->next() ); + } +} + +void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr ) +{ + SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator(); + bool isSatisfy = true; + for ( ; aNodeItr->more(); ) + { + SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); + if ( !isOnSurface( aNode ) ) + { + isSatisfy = false; + break; + } } + if ( isSatisfy ) + myIds.Add( theElemPtr->GetID() ); +} + +bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) const +{ + if ( mySurf.IsNull() ) + return false; + + gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() ); + double aToler2 = myToler * myToler; + if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane))) + { + gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln(); + if ( aPln.SquareDistance( aPnt ) > aToler2 ) + return false; } - return aSequence; + else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface))) + { + gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder(); + double aRad = aCyl.Radius(); + gp_Ax3 anAxis = aCyl.Position(); + gp_XYZ aLoc = aCyl.Location().XYZ(); + double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc ); + double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc ); + if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 ) + return false; + } + else + return false; + + return true; }