X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=1ef3a5f2f4018efa8b7c14ca2454928338ecaa43;hp=c70528f279b3217a6b56f15a13edecac5c20f051;hb=24825596b3595ab0c0b1988368a886864d1080ef;hpb=23d90107acec5e54ded86d9f309fe5cb42720b78 diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index c70528f27..1ef3a5f2f 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2021 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 @@ -23,27 +23,34 @@ #include "SMESH_ControlsDef.hxx" #include "SMDS_BallElement.hxx" +#include "SMDS_FacePosition.hxx" #include "SMDS_Iterator.hxx" #include "SMDS_Mesh.hxx" #include "SMDS_MeshElement.hxx" #include "SMDS_MeshNode.hxx" -#include "SMDS_QuadraticEdge.hxx" -#include "SMDS_QuadraticFaceOfNodes.hxx" #include "SMDS_VolumeTool.hxx" #include "SMESHDS_GroupBase.hxx" +#include "SMESHDS_GroupOnFilter.hxx" #include "SMESHDS_Mesh.hxx" -#include "SMESH_OctreeNode.hxx" #include "SMESH_MeshAlgos.hxx" +#include "SMESH_OctreeNode.hxx" +#include #include #include +#include +#include +#include #include #include +#include #include #include #include +#include #include +#include #include #include #include @@ -91,6 +98,15 @@ namespace { v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 ); } + inline double getCos2( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) + { + gp_Vec v1( P1 - P2 ), v2( P3 - P2 ); + double dot = v1 * v2, len1 = v1.SquareMagnitude(), len2 = v2.SquareMagnitude(); + + return ( dot < 0 || len1 < gp::Resolution() || len2 < gp::Resolution() ? -1 : + dot * dot / len1 / len2 ); + } + inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) { gp_Vec aVec1( P2 - P1 ); @@ -127,16 +143,16 @@ namespace { // Case 1 Case 2 // | | | | | // | | | | | - // +-----+------+ +-----+------+ + // +-----+------+ +-----+------+ // | | | | // | | | | - // result sould be 2 in both cases + // result should be 2 in both cases // int aResult0 = 0, aResult1 = 0; - // last node, it is a medium one in a quadratic edge + // last node, it is a medium one in a quadratic edge const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 ); - const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 ); - const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 ); + const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 ); + const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 ); if ( aNode1 == aLastNode ) aNode1 = 0; SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator(); @@ -173,7 +189,7 @@ namespace { n += q2 ^ q3; } double len = n.Modulus(); - bool zeroLen = ( len <= numeric_limits::min()); + bool zeroLen = ( len <= std::numeric_limits::min()); if ( !zeroLen ) n /= len; @@ -218,7 +234,7 @@ bool NumericalFunctor::GetPoints(const int theId, return false; const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); - if ( !anElem || anElem->GetType() != this->GetType() ) + if ( !IsApplicable( anElem )) return false; return GetPoints( anElem, theRes ); @@ -236,34 +252,12 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, theRes.setElement( anElem ); // Get nodes of the element - SMDS_ElemIteratorPtr anIter; - - if ( anElem->IsQuadratic() ) { - switch ( anElem->GetType() ) { - case SMDSAbs_Edge: - anIter = dynamic_cast - (anElem)->interlacedNodesElemIterator(); - break; - case SMDSAbs_Face: - anIter = dynamic_cast - (anElem)->interlacedNodesElemIterator(); - break; - default: - anIter = anElem->nodesIterator(); - } - } - else { - anIter = anElem->nodesIterator(); - } - + SMDS_NodeIteratorPtr anIter= anElem->interlacedNodesIterator(); if ( anIter ) { - double xyz[3]; + SMESH_NodeXYZ p; while( anIter->more() ) { - if ( const SMDS_MeshNode* aNode = static_cast( anIter->next() )) - { - aNode->GetXYZ( xyz ); - theRes.push_back( gp_XYZ( xyz[0], xyz[1], xyz[2] )); - } + if ( p.Set( anIter->next() )) + theRes.push_back( p ); } } @@ -299,6 +293,24 @@ double NumericalFunctor::Round( const double & aVal ) return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal; } +//================================================================================ +/*! + * \brief Return true if a value can be computed for a given element. + * Some NumericalFunctor's are meaningful for elements of a certain + * geometry only. + */ +//================================================================================ + +bool NumericalFunctor::IsApplicable( const SMDS_MeshElement* element ) const +{ + return element && element->GetType() == this->GetType(); +} + +bool NumericalFunctor::IsApplicable( long theElementId ) const +{ + return IsApplicable( myMesh->FindElement( theElementId )); +} + //================================================================================ /*! * \brief Return histogram of functor values @@ -310,12 +322,12 @@ double NumericalFunctor::Round( const double & aVal ) */ //================================================================================ -void NumericalFunctor::GetHistogram(int nbIntervals, - std::vector& nbEvents, - std::vector& funValues, - const vector& elements, - const double* minmax, - const bool isLogarithmic) +void NumericalFunctor::GetHistogram(int nbIntervals, + std::vector& nbEvents, + std::vector& funValues, + const std::vector& elements, + const double* minmax, + const bool isLogarithmic) { if ( nbIntervals < 1 || !myMesh || @@ -334,7 +346,7 @@ void NumericalFunctor::GetHistogram(int nbIntervals, } else { - vector::const_iterator id = elements.begin(); + std::vector::const_iterator id = elements.begin(); for ( ; id != elements.end(); ++id ) values.insert( GetValue( *id )); } @@ -480,7 +492,7 @@ double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P ) // } // } // { // polygons - + // } if( myPrecision >= 0 ) @@ -520,148 +532,165 @@ double MaxElementLength3D::GetValue( long theElementId ) if( GetPoints( theElementId, P ) ) { double aVal = 0; const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId ); - SMDSAbs_ElementType aType = aElem->GetType(); + SMDSAbs_EntityType aType = aElem->GetEntityType(); int len = P.size(); - switch( aType ) { - case SMDSAbs_Volume: - if( len == 4 ) { // tetras - 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 ) { // pyramids - 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( 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 ) { // pentas - 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 ) { // hexas - 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 )); - double D1 = getDistance(P( 1 ),P( 7 )); - double D2 = getDistance(P( 2 ),P( 8 )); - double D3 = getDistance(P( 3 ),P( 5 )); - double D4 = getDistance(P( 4 ),P( 6 )); - 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)); - aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); - break; - } - else if( len == 12 ) { // hexagonal prism - for ( int i1 = 1; i1 < 12; ++i1 ) - for ( int i2 = i1+1; i1 <= 12; ++i1 ) - aVal = Max( aVal, getDistance(P( i1 ),P( i2 ))); - break; - } - else if( len == 10 ) { // quadratic tetras - double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 )); - double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); - double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 )); - double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 )); - aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); - break; - } - else if( len == 13 ) { // quadratic pyramids - double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); - double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); - double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); - double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 )); - double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 )); - double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 )); - aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); - aVal = Max(aVal,Max(L7,L8)); - break; - } - else if( len == 15 ) { // quadratic pentas - double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); - double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); - double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 )); - double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 )); - double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 )); - double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 )); - double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),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 == 20 || len == 27 ) { // quadratic hexas - double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 )); - double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 )); - double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 )); - double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 )); - double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 )); - double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 )); - double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 )); - double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 )); - double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 )); - double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 )); - double D1 = getDistance(P( 1 ),P( 7 )); - double D2 = getDistance(P( 2 ),P( 8 )); - double D3 = getDistance(P( 3 ),P( 5 )); - double D4 = getDistance(P( 4 ),P( 6 )); - 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)); - aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); - break; - } - else if( len > 1 && aElem->IsPoly() ) { // polys - // get the maximum distance between all pairs of nodes - for( int i = 1; i <= len; i++ ) { - for( int j = 1; j <= len; j++ ) { - if( j > i ) { // optimization of the loop - double D = getDistance( P(i), P(j) ); - aVal = Max( aVal, D ); - } + switch ( aType ) { + case SMDSEntity_Tetra: { // tetras + 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; + } + case SMDSEntity_Pyramid: { // pyramids + 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( 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; + } + case SMDSEntity_Penta: { // pentas + 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; + } + case SMDSEntity_Hexa: { // hexas + 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 )); + double D1 = getDistance(P( 1 ),P( 7 )); + double D2 = getDistance(P( 2 ),P( 8 )); + double D3 = getDistance(P( 3 ),P( 5 )); + double D4 = getDistance(P( 4 ),P( 6 )); + 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)); + aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); + break; + } + case SMDSEntity_Hexagonal_Prism: { // hexagonal prism + for ( int i1 = 1; i1 < 12; ++i1 ) + for ( int i2 = i1+1; i1 <= 12; ++i1 ) + aVal = Max( aVal, getDistance(P( i1 ),P( i2 ))); + break; + } + case SMDSEntity_Quad_Tetra: { // quadratic tetras + double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 )); + double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 )); + double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + break; + } + case SMDSEntity_Quad_Pyramid: { // quadratic pyramids + double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 )); + double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 )); + double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(L7,L8)); + break; + } + case SMDSEntity_Quad_Penta: + case SMDSEntity_BiQuad_Penta: { // quadratic pentas + double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 )); + double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 )); + double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 )); + double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),L9)); + break; + } + case SMDSEntity_Quad_Hexa: + case SMDSEntity_TriQuad_Hexa: { // quadratic hexas + double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 )); + double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 )); + double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 )); + double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 )); + double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 )); + double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 )); + double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 )); + double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 )); + double D1 = getDistance(P( 1 ),P( 7 )); + double D2 = getDistance(P( 2 ),P( 8 )); + double D3 = getDistance(P( 3 ),P( 5 )); + double D4 = getDistance(P( 4 ),P( 6 )); + 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)); + aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); + break; + } + case SMDSEntity_Quad_Polyhedra: + case SMDSEntity_Polyhedra: { // polys + // get the maximum distance between all pairs of nodes + for( int i = 1; i <= len; i++ ) { + for( int j = 1; j <= len; j++ ) { + if( j > i ) { // optimization of the loop + double D = getDistance( P(i), P(j) ); + aVal = Max( aVal, D ); } } } + break; } + case SMDSEntity_Node: + case SMDSEntity_0D: + case SMDSEntity_Edge: + case SMDSEntity_Quad_Edge: + case SMDSEntity_Triangle: + case SMDSEntity_Quad_Triangle: + case SMDSEntity_BiQuad_Triangle: + case SMDSEntity_Quadrangle: + case SMDSEntity_Quad_Quadrangle: + case SMDSEntity_BiQuad_Quadrangle: + case SMDSEntity_Polygon: + case SMDSEntity_Quad_Polygon: + case SMDSEntity_Ball: + case SMDSEntity_Last: return 0; + } // switch ( aType ) if( myPrecision >= 0 ) { @@ -692,21 +721,25 @@ SMDSAbs_ElementType MaxElementLength3D::GetType() const double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) { - double aMin; - - if (P.size() <3) + if ( P.size() < 3 ) return 0.; - aMin = getAngle(P( P.size() ), P( 1 ), P( 2 )); - aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 ))); + double aMaxCos2; - for ( int i = 2; i < P.size(); i++ ) + aMaxCos2 = getCos2( P( P.size() ), P( 1 ), P( 2 )); + aMaxCos2 = Max( aMaxCos2, getCos2( P( P.size()-1 ), P( P.size() ), P( 1 ))); + + for ( size_t i = 2; i < P.size(); i++ ) { - double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) ); - aMin = Min(aMin,A0); + double A0 = getCos2( P( i-1 ), P( i ), P( i+1 ) ); + aMaxCos2 = Max( aMaxCos2, A0 ); } + if ( aMaxCos2 < 0 ) + return 0; // all nodes coincide - return aMin * 180.0 / M_PI; + double cos = sqrt( aMaxCos2 ); + if ( cos >= 1 ) return 0; + return acos( cos ) * 180.0 / M_PI; } double MinimumAngle::GetBadRate( double Value, int nbNodes ) const @@ -733,19 +766,9 @@ double AspectRatio::GetValue( long theId ) { double aVal = 0; myCurrElement = myMesh->FindElement( theId ); - if ( myCurrElement && myCurrElement->GetVtkType() == VTK_QUAD ) - { - // issue 21723 - vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid(); - if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() )) - aVal = Round( vtkMeshQuality::QuadAspectRatio( avtkCell )); - } - else - { - TSequenceOfXYZ P; - if ( GetPoints( myCurrElement, P )) - aVal = Round( GetValue( P )); - } + TSequenceOfXYZ P; + if ( GetPoints( myCurrElement, P )) + aVal = Round( GetValue( P )); return aVal; } @@ -765,58 +788,51 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) if ( nbNodes == 3 ) { // Compute lengths of the sides - std::vector< double > aLen (nbNodes); - for ( int i = 0; i < nbNodes - 1; i++ ) - aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) ); - aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) ); + double aLen1 = getDistance( P( 1 ), P( 2 )); + double aLen2 = getDistance( P( 2 ), P( 3 )); + double aLen3 = getDistance( P( 3 ), P( 1 )); // Q = alfa * h * p / S, where // // alfa = sqrt( 3 ) / 6 // h - length of the longest edge // p - half perimeter // S - triangle surface - const double alfa = sqrt( 3. ) / 6.; - double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) ); - double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.; - double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) ); + const double alfa = sqrt( 3. ) / 6.; + double maxLen = Max( aLen1, Max( aLen2, aLen3 )); + double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.; + double anArea = getArea( P( 1 ), P( 2 ), P( 3 )); if ( anArea <= theEps ) return theInf; return alfa * maxLen * half_perimeter / anArea; } else if ( nbNodes == 6 ) { // quadratic triangles // Compute lengths of the sides - std::vector< double > aLen (3); - aLen[0] = getDistance( P(1), P(3) ); - aLen[1] = getDistance( P(3), P(5) ); - aLen[2] = getDistance( P(5), P(1) ); - // Q = alfa * h * p / S, where - // - // alfa = sqrt( 3 ) / 6 - // h - length of the longest edge - // p - half perimeter - // S - triangle surface - const double alfa = sqrt( 3. ) / 6.; - double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) ); - double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.; - double anArea = getArea( P(1), P(3), P(5) ); + double aLen1 = getDistance( P( 1 ), P( 3 )); + double aLen2 = getDistance( P( 3 ), P( 5 )); + double aLen3 = getDistance( P( 5 ), P( 1 )); + // algo same as for the linear triangle + const double alfa = sqrt( 3. ) / 6.; + double maxLen = Max( aLen1, Max( aLen2, aLen3 )); + double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.; + double anArea = getArea( P( 1 ), P( 3 ), P( 5 )); if ( anArea <= theEps ) return theInf; return alfa * maxLen * half_perimeter / anArea; } else if( nbNodes == 4 ) { // quadrangle // Compute lengths of the sides - std::vector< double > aLen (4); + double aLen[4]; aLen[0] = getDistance( P(1), P(2) ); aLen[1] = getDistance( P(2), P(3) ); aLen[2] = getDistance( P(3), P(4) ); aLen[3] = getDistance( P(4), P(1) ); // Compute lengths of the diagonals - std::vector< double > aDia (2); + double aDia[2]; aDia[0] = getDistance( P(1), P(3) ); aDia[1] = getDistance( P(2), P(4) ); // Compute areas of all triangles which can be built // taking three nodes of the quadrangle - std::vector< double > anArea (4); + double anArea[4]; anArea[0] = getArea( P(1), P(2), P(3) ); anArea[1] = getArea( P(1), P(2), P(4) ); anArea[2] = getArea( P(1), P(3), P(4) ); @@ -825,42 +841,42 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) // // alpha = sqrt( 1/32 ) // L = max( L1, L2, L3, L4, D1, D2 ) - // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 ) + // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 ) // C2 = min( S1, S2, S3, S4 ) // Li - lengths of the edges // Di - lengths of the diagonals // Si - areas of the triangles const double alpha = sqrt( 1 / 32. ); double L = Max( aLen[ 0 ], - Max( aLen[ 1 ], - Max( aLen[ 2 ], - Max( aLen[ 3 ], - Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) ); - double C1 = sqrt( ( aLen[0] * aLen[0] + - aLen[1] * aLen[1] + - aLen[2] * aLen[2] + - aLen[3] * aLen[3] ) / 4. ); + Max( aLen[ 1 ], + Max( aLen[ 2 ], + Max( aLen[ 3 ], + Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) ); + double C1 = sqrt( aLen[0] * aLen[0] + + aLen[1] * aLen[1] + + aLen[2] * aLen[2] + + aLen[3] * aLen[3] ); double C2 = Min( anArea[ 0 ], - Min( anArea[ 1 ], - Min( anArea[ 2 ], anArea[ 3 ] ) ) ); + Min( anArea[ 1 ], + Min( anArea[ 2 ], anArea[ 3 ] ) ) ); if ( C2 <= theEps ) return theInf; return alpha * L * C1 / C2; } else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle // Compute lengths of the sides - std::vector< double > aLen (4); + double aLen[4]; aLen[0] = getDistance( P(1), P(3) ); aLen[1] = getDistance( P(3), P(5) ); aLen[2] = getDistance( P(5), P(7) ); aLen[3] = getDistance( P(7), P(1) ); // Compute lengths of the diagonals - std::vector< double > aDia (2); + double aDia[2]; aDia[0] = getDistance( P(1), P(5) ); aDia[1] = getDistance( P(3), P(7) ); // Compute areas of all triangles which can be built // taking three nodes of the quadrangle - std::vector< double > anArea (4); + double anArea[4]; anArea[0] = getArea( P(1), P(3), P(5) ); anArea[1] = getArea( P(1), P(3), P(7) ); anArea[2] = getArea( P(1), P(5), P(7) ); @@ -869,24 +885,24 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) // // alpha = sqrt( 1/32 ) // L = max( L1, L2, L3, L4, D1, D2 ) - // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 ) + // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 ) // C2 = min( S1, S2, S3, S4 ) // Li - lengths of the edges // Di - lengths of the diagonals // Si - areas of the triangles const double alpha = sqrt( 1 / 32. ); double L = Max( aLen[ 0 ], - Max( aLen[ 1 ], - Max( aLen[ 2 ], - Max( aLen[ 3 ], - Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) ); - double C1 = sqrt( ( aLen[0] * aLen[0] + - aLen[1] * aLen[1] + - aLen[2] * aLen[2] + - aLen[3] * aLen[3] ) / 4. ); + Max( aLen[ 1 ], + Max( aLen[ 2 ], + Max( aLen[ 3 ], + Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) ); + double C1 = sqrt( aLen[0] * aLen[0] + + aLen[1] * aLen[1] + + aLen[2] * aLen[2] + + aLen[3] * aLen[3] ); double C2 = Min( anArea[ 0 ], - Min( anArea[ 1 ], - Min( anArea[ 2 ], anArea[ 3 ] ) ) ); + Min( anArea[ 1 ], + Min( anArea[ 2 ], anArea[ 3 ] ) ) ); if ( C2 <= theEps ) return theInf; return alpha * L * C1 / C2; @@ -894,6 +910,11 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) return 0; } +bool AspectRatio::IsApplicable( const SMDS_MeshElement* element ) const +{ + return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() ); +} + double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const { // the aspect ratio is in the range [1.0,infinity] @@ -976,6 +997,102 @@ namespace{ return aHeight; } + //================================================================================ + /*! + * \brief Standard quality of a tetrahedron but not normalized + */ + //================================================================================ + + double tetQualityByHomardMethod( const gp_XYZ & p1, + const gp_XYZ & p2, + const gp_XYZ & p3, + const gp_XYZ & p4 ) + { + gp_XYZ edgeVec[6]; + edgeVec[0] = ( p1 - p2 ); + edgeVec[1] = ( p2 - p3 ); + edgeVec[2] = ( p3 - p1 ); + edgeVec[3] = ( p4 - p1 ); + edgeVec[4] = ( p4 - p2 ); + edgeVec[5] = ( p4 - p3 ); + + double maxEdgeLen2 = edgeVec[0].SquareModulus(); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[1].SquareModulus() ); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[2].SquareModulus() ); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[3].SquareModulus() ); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[4].SquareModulus() ); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[5].SquareModulus() ); + double maxEdgeLen = Sqrt( maxEdgeLen2 ); + + gp_XYZ cross01 = edgeVec[0] ^ edgeVec[1]; + double sumArea = ( cross01 ).Modulus(); // actually double area + sumArea += ( edgeVec[0] ^ edgeVec[3] ).Modulus(); + sumArea += ( edgeVec[1] ^ edgeVec[4] ).Modulus(); + sumArea += ( edgeVec[2] ^ edgeVec[5] ).Modulus(); + + double sixVolume = Abs( cross01 * edgeVec[4] ); // 6 * volume + double quality = maxEdgeLen * sumArea / sixVolume; // not normalized!!! + return quality; + } + + //================================================================================ + /*! + * \brief HOMARD method of hexahedron quality + * 1. Decompose the hexa into 24 tetra: each face is splitted into 4 triangles by + * adding the diagonals and every triangle is connected to the center of the hexa. + * 2. Compute the quality of every tetra with the same formula as for the standard quality, + * except that the factor for the normalization is not the same because the final goal + * is to have a quality equal to 1 for a perfect cube. So the formula is: + * qual = max(lengthes of 6 edges) * (sum of surfaces of 4 faces) / (7.6569*6*volume) + * 3. The quality of the hexa is the highest value of the qualities of the 24 tetra + */ + //================================================================================ + + double hexQualityByHomardMethod( const TSequenceOfXYZ& P ) + { + gp_XYZ quadCenter[6]; + quadCenter[0] = ( P(1) + P(2) + P(3) + P(4) ) / 4.; + quadCenter[1] = ( P(5) + P(6) + P(7) + P(8) ) / 4.; + quadCenter[2] = ( P(1) + P(2) + P(6) + P(5) ) / 4.; + quadCenter[3] = ( P(2) + P(3) + P(7) + P(6) ) / 4.; + quadCenter[4] = ( P(3) + P(4) + P(8) + P(7) ) / 4.; + quadCenter[5] = ( P(1) + P(4) + P(8) + P(5) ) / 4.; + + gp_XYZ hexCenter = ( P(1) + P(2) + P(3) + P(4) + P(5) + P(6) + P(7) + P(8) ) / 8.; + + // quad 1 ( 1 2 3 4 ) + double quality = tetQualityByHomardMethod( P(1), P(2), quadCenter[0], hexCenter ); + quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[0], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[0], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(4), P(1), quadCenter[0], hexCenter )); + // quad 2 ( 5 6 7 8 ) + quality = Max( quality, tetQualityByHomardMethod( P(5), P(6), quadCenter[1], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(6), P(7), quadCenter[1], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(7), P(8), quadCenter[1], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[1], hexCenter )); + // quad 3 ( 1 2 6 5 ) + quality = Max( quality, tetQualityByHomardMethod( P(1), P(2), quadCenter[2], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(2), P(6), quadCenter[2], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(6), P(5), quadCenter[2], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[2], hexCenter )); + // quad 4 ( 2 3 7 6 ) + quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[3], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(3), P(7), quadCenter[3], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(7), P(6), quadCenter[3], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(6), P(2), quadCenter[3], hexCenter )); + // quad 5 ( 3 4 8 7 ) + quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[4], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[4], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(8), P(7), quadCenter[4], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(7), P(3), quadCenter[4], hexCenter )); + // quad 6 ( 1 4 8 5 ) + quality = Max( quality, tetQualityByHomardMethod( P(1), P(4), quadCenter[5], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[5], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[5], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[5], hexCenter )); + + return quality / 7.65685424949; + } } double AspectRatio3D::GetValue( long theId ) @@ -987,8 +1104,8 @@ double AspectRatio3D::GetValue( long theId ) // Action from CoTech | ACTION 31.3: // EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with // those of ParaView. The library used by ParaView for those calculations can be reused in SMESH. - vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid(); - if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() )) + vtkUnstructuredGrid* grid = const_cast( myMesh )->GetGrid(); + if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->GetVtkID() )) aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell )); } else @@ -1000,6 +1117,11 @@ double AspectRatio3D::GetValue( long theId ) return aVal; } +bool AspectRatio3D::IsApplicable( const SMDS_MeshElement* element ) const +{ + return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() ); +} + double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) { double aQuality = 0.0; @@ -1007,12 +1129,12 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) int nbNodes = P.size(); - if(myCurrElement->IsQuadratic()) { - if(nbNodes==10) nbNodes=4; // quadratic tetrahedron + if( myCurrElement->IsQuadratic() ) { + if (nbNodes==10) nbNodes=4; // quadratic tetrahedron else if(nbNodes==13) nbNodes=5; // quadratic pyramid else if(nbNodes==15) nbNodes=6; // quadratic pentahedron else if(nbNodes==20) nbNodes=8; // quadratic hexahedron - else if(nbNodes==27) nbNodes=8; // quadratic hexahedron + else if(nbNodes==27) nbNodes=8; // tri-quadratic hexahedron else return aQuality; } @@ -1056,7 +1178,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) case 5:{ { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )}; - aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )}; @@ -1075,7 +1197,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) case 6:{ { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )}; - aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )}; @@ -1100,9 +1222,13 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) break; } case 8:{ + + return hexQualityByHomardMethod( P ); // bos #23982 + + { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )}; - aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )}; @@ -1237,7 +1363,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) case 12: { gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )}; - aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality); + aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])); } { gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )}; @@ -1251,7 +1377,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) } // switch(nbNodes) if ( nbNodes > 4 ) { - // avaluate aspect ratio of quadranle faces + // evaluate aspect ratio of quadrangle faces AspectRatio aspect2D; SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes ); int nbFaces = SMDS_VolumeTool::NbFaces( type ); @@ -1260,7 +1386,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 ) continue; const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true ); - for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face + for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadrangle face points( p + 1 ) = P( pInd[ p ] + 1 ); aQuality = std::max( aQuality, aspect2D.GetValue( points )); } @@ -1289,6 +1415,11 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const */ //================================================================================ +bool Warping::IsApplicable( const SMDS_MeshElement* element ) const +{ + return NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4; +} + double Warping::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 4 ) @@ -1353,6 +1484,11 @@ SMDSAbs_ElementType Warping::GetType() const */ //================================================================================ +bool Taper::IsApplicable( const SMDS_MeshElement* element ) const +{ + return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4 ); +} + double Taper::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 4 ) @@ -1411,6 +1547,11 @@ static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 ); } +bool Skew::IsApplicable( const SMDS_MeshElement* element ) const +{ + return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() <= 4 ); +} + double Skew::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 3 && P.size() != 4 ) @@ -1475,10 +1616,10 @@ double Area::GetValue( const TSequenceOfXYZ& P ) gp_Vec aVec2( P(3) - P(1) ); gp_Vec SumVec = aVec1 ^ aVec2; - for (int i=4; i<=P.size(); i++) + for (size_t i=4; i<=P.size(); i++) { gp_Vec aVec1( P(i-1) - P(1) ); - gp_Vec aVec2( P(i) - P(1) ); + gp_Vec aVec2( P(i ) - P(1) ); gp_Vec tmp = aVec1 ^ aVec2; SumVec.Add(tmp); } @@ -1525,253 +1666,270 @@ SMDSAbs_ElementType Length::GetType() const return SMDSAbs_Edge; } +//================================================================================ +/* + Class : Length3D + Description : Functor for calculating minimal length of element edge +*/ +//================================================================================ + +Length3D::Length3D(): + Length2D ( SMDSAbs_Volume ) +{ +} + //================================================================================ /* Class : Length2D - Description : Functor for calculating minimal length of edge + Description : Functor for calculating minimal length of element edge */ //================================================================================ -double Length2D::GetValue( long theElementId ) +Length2D::Length2D( SMDSAbs_ElementType type ): + myType ( type ) { - TSequenceOfXYZ P; +} - if ( GetPoints( theElementId, P )) - { - double aVal = 0; - int len = P.size(); - SMDSAbs_EntityType aType = P.getElementEntity(); +bool Length2D::IsApplicable( const SMDS_MeshElement* element ) const +{ + return ( NumericalFunctor::IsApplicable( element ) && + element->GetEntityType() != SMDSEntity_Polyhedra ); +} - switch (aType) { - case SMDSEntity_Edge: - if (len == 2) - aVal = getDistance( P( 1 ), P( 2 ) ); - break; - case SMDSEntity_Quad_Edge: - if (len == 3) // quadratic edge - aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 )); - break; - case SMDSEntity_Triangle: - 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 = Min(L1,Min(L2,L3)); - } - break; - case SMDSEntity_Quadrangle: - 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 = Min(Min(L1,L2),Min(L3,L4)); - } - break; - case SMDSEntity_Quad_Triangle: - case SMDSEntity_BiQuad_Triangle: - 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 = Min(L1,Min(L2,L3)); - } - break; - case SMDSEntity_Quad_Quadrangle: - case SMDSEntity_BiQuad_Quadrangle: - if (len >= 8){ // quadratic quadrangles - 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( 7 )); - double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 )); - aVal = Min(Min(L1,L2),Min(L3,L4)); - } - break; - case SMDSEntity_Tetra: - if (len == 4){ // tetrahedra - 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); - } - break; - case SMDSEntity_Pyramid: - if (len == 5){ // piramids - 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( 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); - aVal = Min(aVal,Min(L7,L8)); - } - break; - case SMDSEntity_Penta: - 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); - aVal = Min(aVal,Min(Min(L7,L8),L9)); - } - break; - case SMDSEntity_Hexa: - if (len == 8){ // hexahedron - 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); - aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10))); - aVal = Min(aVal,Min(L11,L12)); - } - break; - case SMDSEntity_Quad_Tetra: - if (len == 10){ // quadratic tetraidrs - double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 )); - double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); - double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 )); - double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 )); - aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); - } - break; - case SMDSEntity_Quad_Pyramid: - if (len == 13){ // quadratic piramids - double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); - double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); - double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); - double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 )); - double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 )); - double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 )); - aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); - aVal = Min(aVal,Min(L7,L8)); - } - break; - case SMDSEntity_Quad_Penta: - if (len == 15){ // quadratic pentaidres - double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); - double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); - double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 )); - double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 )); - double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 )); - double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 )); - double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 )); - aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); - aVal = Min(aVal,Min(Min(L7,L8),L9)); - } - break; - case SMDSEntity_Quad_Hexa: - case SMDSEntity_TriQuad_Hexa: - if (len >= 20) { // quadratic hexaider - double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 )); - double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 )); - double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 )); - double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 )); - double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 )); - double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 )); - double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 )); - double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 )); - double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 )); - double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 )); - aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); - aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10))); - aVal = Min(aVal,Min(L11,L12)); - } - break; - case SMDSEntity_Polygon: - if ( len > 1 ) { - aVal = getDistance( P(1), P( P.size() )); - for ( size_t i = 1; i < P.size(); ++i ) - aVal = Min( aVal, getDistance( P( i ), P( i+1 ))); - } - break; - case SMDSEntity_Quad_Polygon: - if ( len > 2 ) { - aVal = getDistance( P(1), P( P.size() )) + getDistance( P(P.size()), P( P.size()-1 )); - for ( size_t i = 1; i < P.size()-1; i += 2 ) - aVal = Min( aVal, getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 ))); - } - break; - case SMDSEntity_Hexagonal_Prism: - if (len == 12) { // hexagonal prism - 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( 5 )); - double L5 = getDistance(P( 5 ),P( 6 )); - double L6 = getDistance(P( 6 ),P( 1 )); - - double L7 = getDistance(P( 7 ), P( 8 )); - double L8 = getDistance(P( 8 ), P( 9 )); - double L9 = getDistance(P( 9 ), P( 10 )); - double L10= getDistance(P( 10 ),P( 11 )); - double L11= getDistance(P( 11 ),P( 12 )); - double L12= getDistance(P( 12 ),P( 7 )); - - double L13 = getDistance(P( 1 ),P( 7 )); - double L14 = getDistance(P( 2 ),P( 8 )); - double L15 = getDistance(P( 3 ),P( 9 )); - double L16 = getDistance(P( 4 ),P( 10 )); - double L17 = getDistance(P( 5 ),P( 11 )); - double L18 = getDistance(P( 6 ),P( 12 )); - aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); - aVal = Min(aVal, Min(Min(Min(L7,L8),Min(L9,L10)),Min(L11,L12))); - aVal = Min(aVal, Min(Min(Min(L13,L14),Min(L15,L16)),Min(L17,L18))); - } - break; - case SMDSEntity_Polyhedra: - { +double Length2D::GetValue( const TSequenceOfXYZ& P ) +{ + double aVal = 0; + int len = P.size(); + SMDSAbs_EntityType aType = P.getElementEntity(); + + switch (aType) { + case SMDSEntity_Edge: + if (len == 2) + aVal = getDistance( P( 1 ), P( 2 ) ); + break; + case SMDSEntity_Quad_Edge: + if (len == 3) // quadratic edge + aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 )); + break; + case SMDSEntity_Triangle: + 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 = Min(L1,Min(L2,L3)); } break; - default: - return 0; + case SMDSEntity_Quadrangle: + 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 = Min(Min(L1,L2),Min(L3,L4)); } - - if (aVal < 0 ) { - return 0.; + break; + case SMDSEntity_Quad_Triangle: + case SMDSEntity_BiQuad_Triangle: + 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 = Min(L1,Min(L2,L3)); } - - if ( myPrecision >= 0 ) - { - double prec = pow( 10., (double)( myPrecision ) ); - aVal = floor( aVal * prec + 0.5 ) / prec; + break; + case SMDSEntity_Quad_Quadrangle: + case SMDSEntity_BiQuad_Quadrangle: + if (len >= 8){ // quadratic quadrangles + 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( 7 )); + double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 )); + aVal = Min(Min(L1,L2),Min(L3,L4)); + } + break; + case SMDSEntity_Tetra: + if (len == 4){ // tetrahedra + 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + } + break; + case SMDSEntity_Pyramid: + if (len == 5){ // pyramid + 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( 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(L7,L8)); + } + break; + case SMDSEntity_Penta: + if (len == 6) { // pentahedron + 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(Min(L7,L8),L9)); + } + break; + case SMDSEntity_Hexa: + if (len == 8){ // hexahedron + 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 = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10))); + aVal = Min(aVal,Min(L11,L12)); + } + break; + case SMDSEntity_Quad_Tetra: + if (len == 10){ // quadratic tetrahedron + double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 )); + double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 )); + double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + } + break; + case SMDSEntity_Quad_Pyramid: + if (len == 13){ // quadratic pyramid + double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 )); + double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 )); + double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(L7,L8)); + } + break; + case SMDSEntity_Quad_Penta: + case SMDSEntity_BiQuad_Penta: + if (len >= 15){ // quadratic pentahedron + double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 )); + double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 )); + double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 )); + double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(Min(L7,L8),L9)); + } + break; + case SMDSEntity_Quad_Hexa: + case SMDSEntity_TriQuad_Hexa: + if (len >= 20) { // quadratic hexahedron + double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 )); + double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 )); + double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 )); + double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 )); + double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 )); + double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 )); + double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 )); + double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10))); + aVal = Min(aVal,Min(L11,L12)); + } + break; + case SMDSEntity_Polygon: + if ( len > 1 ) { + aVal = getDistance( P(1), P( P.size() )); + for ( size_t i = 1; i < P.size(); ++i ) + aVal = Min( aVal, getDistance( P( i ), P( i+1 ))); + } + break; + case SMDSEntity_Quad_Polygon: + if ( len > 2 ) { + aVal = getDistance( P(1), P( P.size() )) + getDistance( P(P.size()), P( P.size()-1 )); + for ( size_t i = 1; i < P.size()-1; i += 2 ) + aVal = Min( aVal, getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 ))); + } + break; + case SMDSEntity_Hexagonal_Prism: + if (len == 12) { // hexagonal prism + 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( 5 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 1 )); + + double L7 = getDistance(P( 7 ), P( 8 )); + double L8 = getDistance(P( 8 ), P( 9 )); + double L9 = getDistance(P( 9 ), P( 10 )); + double L10= getDistance(P( 10 ),P( 11 )); + double L11= getDistance(P( 11 ),P( 12 )); + double L12= getDistance(P( 12 ),P( 7 )); + + double L13 = getDistance(P( 1 ),P( 7 )); + double L14 = getDistance(P( 2 ),P( 8 )); + double L15 = getDistance(P( 3 ),P( 9 )); + double L16 = getDistance(P( 4 ),P( 10 )); + double L17 = getDistance(P( 5 ),P( 11 )); + double L18 = getDistance(P( 6 ),P( 12 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal, Min(Min(Min(L7,L8),Min(L9,L10)),Min(L11,L12))); + aVal = Min(aVal, Min(Min(Min(L13,L14),Min(L15,L16)),Min(L17,L18))); } + break; + case SMDSEntity_Polyhedra: + { + } + break; + default: + return 0; + } - return aVal; + if (aVal < 0 ) { + return 0.; + } + if ( myPrecision >= 0 ) + { + double prec = pow( 10., (double)( myPrecision ) ); + aVal = floor( aVal * prec + 0.5 ) / prec; } - return 0.; + + return aVal; } double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -1782,7 +1940,7 @@ double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const SMDSAbs_ElementType Length2D::GetType() const { - return SMDSAbs_Face; + return myType; } Length2D::Value::Value(double theLength,long thePntId1, long thePntId2): @@ -1804,86 +1962,180 @@ bool Length2D::Value::operator<(const Length2D::Value& x) const void Length2D::GetValues(TValues& theValues) { - TValues aValues; - SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); - for(; anIter->more(); ){ - const SMDS_MeshFace* anElem = anIter->next(); + if ( myType == SMDSAbs_Face ) + { + for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); ) + { + const SMDS_MeshFace* anElem = anIter->next(); + if ( anElem->IsQuadratic() ) + { + // use special nodes iterator + SMDS_NodeIteratorPtr anIter = anElem->interlacedNodesIterator(); + long aNodeId[4] = { 0,0,0,0 }; + gp_Pnt P[4]; - if(anElem->IsQuadratic()) { - const SMDS_VtkFace* F = - dynamic_cast(anElem); - // use special nodes iterator - SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator(); - long aNodeId[4]; - gp_Pnt P[4]; - - double aLength; - const SMDS_MeshElement* aNode; - if(anIter->more()){ - aNode = anIter->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(; anIter->more(); ){ - const SMDS_MeshNode* N1 = static_cast (anIter->next()); - P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z()); - aNodeId[2] = N1->GetID(); - aLength = P[1].Distance(P[2]); - if(!anIter->more()) break; - const SMDS_MeshNode* N2 = static_cast (anIter->next()); - P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z()); - aNodeId[3] = N2->GetID(); - aLength += P[2].Distance(P[3]); + double aLength = 0; + if ( anIter->more() ) + { + const SMDS_MeshNode* aNode = anIter->next(); + P[0] = P[1] = SMESH_NodeXYZ( aNode ); + aNodeId[0] = aNodeId[1] = aNode->GetID(); + aLength = 0; + } + for ( ; anIter->more(); ) + { + const SMDS_MeshNode* N1 = anIter->next(); + P[2] = SMESH_NodeXYZ( N1 ); + aNodeId[2] = N1->GetID(); + aLength = P[1].Distance(P[2]); + if(!anIter->more()) break; + const SMDS_MeshNode* N2 = anIter->next(); + P[3] = SMESH_NodeXYZ( N2 ); + aNodeId[3] = N2->GetID(); + aLength += P[2].Distance(P[3]); + Value aValue1(aLength,aNodeId[1],aNodeId[2]); + Value aValue2(aLength,aNodeId[2],aNodeId[3]); + P[1] = P[3]; + aNodeId[1] = aNodeId[3]; + theValues.insert(aValue1); + theValues.insert(aValue2); + } + aLength += P[2].Distance(P[0]); Value aValue1(aLength,aNodeId[1],aNodeId[2]); - Value aValue2(aLength,aNodeId[2],aNodeId[3]); - P[1] = P[3]; - aNodeId[1] = aNodeId[3]; + Value aValue2(aLength,aNodeId[2],aNodeId[0]); theValues.insert(aValue1); theValues.insert(aValue2); } - aLength += P[2].Distance(P[0]); - Value aValue1(aLength,aNodeId[1],aNodeId[2]); - Value aValue2(aLength,aNodeId[2],aNodeId[0]); - theValues.insert(aValue1); - theValues.insert(aValue2); - } - else { - 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]); + else { + SMDS_NodeIteratorPtr aNodesIter = anElem->nodeIterator(); + long aNodeId[2] = {0,0}; + gp_Pnt P[3]; - Value aValue(aLength,aNodeId[0],aNodeId[1]); - theValues.insert(aValue); - } + double aLength; + const SMDS_MeshElement* aNode; + if ( aNodesIter->more()) + { + aNode = aNodesIter->next(); + P[0] = P[1] = SMESH_NodeXYZ( aNode ); + aNodeId[0] = aNodeId[1] = aNode->GetID(); + aLength = 0; + } + for( ; aNodesIter->more(); ) + { + aNode = aNodesIter->next(); + long anId = aNode->GetID(); + + P[2] = SMESH_NodeXYZ( aNode ); + + 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); + } + } } + else + { + // not implemented + } +} + +//================================================================================ +/* + Class : Deflection2D + Description : computes distance between a face center and an underlying surface +*/ +//================================================================================ + +double Deflection2D::GetValue( const TSequenceOfXYZ& P ) +{ + if ( myMesh && P.getElement() ) + { + // get underlying surface + if ( myShapeIndex != P.getElement()->getshapeId() ) + { + mySurface.Nullify(); + myShapeIndex = P.getElement()->getshapeId(); + const TopoDS_Shape& S = + static_cast< const SMESHDS_Mesh* >( myMesh )->IndexToShape( myShapeIndex ); + if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE ) + { + mySurface = new ShapeAnalysis_Surface( BRep_Tool::Surface( TopoDS::Face( S ))); + + GeomLib_IsPlanarSurface isPlaneCheck( mySurface->Surface() ); + if ( isPlaneCheck.IsPlanar() ) + myPlane.reset( new gp_Pln( isPlaneCheck.Plan() )); + else + myPlane.reset(); + } + } + // project gravity center to the surface + if ( !mySurface.IsNull() ) + { + gp_XYZ gc(0,0,0); + gp_XY uv(0,0); + int nbUV = 0; + for ( size_t i = 0; i < P.size(); ++i ) + { + gc += P(i+1); + + if ( SMDS_FacePositionPtr fPos = P.getElement()->GetNode( i )->GetPosition() ) + { + uv.ChangeCoord(1) += fPos->GetUParameter(); + uv.ChangeCoord(2) += fPos->GetVParameter(); + ++nbUV; + } + } + gc /= P.size(); + if ( nbUV ) uv /= nbUV; + + double maxLen = MaxElementLength2D().GetValue( P ); + double tol = 1e-3 * maxLen; + double dist; + if ( myPlane ) + { + dist = myPlane->Distance( gc ); + if ( dist < tol ) + dist = 0; + } + else + { + if ( uv.X() != 0 && uv.Y() != 0 ) // faster way + mySurface->NextValueOfUV( uv, gc, tol, 0.5 * maxLen ); + else + mySurface->ValueOfUV( gc, tol ); + dist = mySurface->Gap(); + } + return Round( dist ); + } + } + return 0; +} + +void Deflection2D::SetMesh( const SMDS_Mesh* theMesh ) +{ + NumericalFunctor::SetMesh( dynamic_cast( theMesh )); + myShapeIndex = -100; + myPlane.reset(); +} + +SMDSAbs_ElementType Deflection2D::GetType() const +{ + return SMDSAbs_Face; +} + +double Deflection2D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + // meaningless as it is not quality control functor + return Value; } //================================================================================ @@ -1893,7 +2145,7 @@ void Length2D::GetValues(TValues& theValues) */ //================================================================================ -double MultiConnection::GetValue( const TSequenceOfXYZ& P ) +double MultiConnection::GetValue( const TSequenceOfXYZ& /*P*/ ) { return 0; } @@ -1920,7 +2172,7 @@ SMDSAbs_ElementType MultiConnection::GetType() const */ //================================================================================ -double MultiConnection2D::GetValue( const TSequenceOfXYZ& P ) +double MultiConnection2D::GetValue( const TSequenceOfXYZ& /*P*/ ) { return 0; } @@ -1939,7 +2191,7 @@ double MultiConnection2D::GetValue( long theElementId ) SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator(); if (!anIter) break; - const SMDS_MeshNode *aNode, *aNode0; + const SMDS_MeshNode *aNode, *aNode0 = 0; TColStd_MapOfInteger aMap, aMapPrev; for (i = 0; i <= len; i++) { @@ -2011,59 +2263,24 @@ 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(); - SMDS_ElemIteratorPtr aNodesIter; - if ( anElem->IsQuadratic() ) - aNodesIter = dynamic_cast - (anElem)->interlacedNodesElemIterator(); - else - aNodesIter = anElem->nodesIterator(); - long aNodeId[3]; + for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); ) + { + const SMDS_MeshFace* anElem = anIter->next(); + SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator(); - //int aNbConnects=0; - const SMDS_MeshNode* aNode0; - const SMDS_MeshNode* aNode1; + const SMDS_MeshNode* aNode1 = anElem->GetNode( anElem->NbNodes() - 1 ); 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: "<SetMesh( theMesh ); - SMDS_ElemIteratorPtr elemIt = theMesh->elementsIterator( thePredicate->GetType() ); - if ( elemIt ) { - while ( elemIt->more() ) { - const SMDS_MeshElement* anElem = elemIt->next(); - long anId = anElem->GetID(); - if ( thePredicate->IsSatisfy( anId ) ) - theSequence.push_back( anId ); + if ( !theElements ) + theElements = theMesh->elementsIterator( thePredicate->GetType() ); + + if ( theElements ) { + while ( theElements->more() ) { + const SMDS_MeshElement* anElem = theElements->next(); + if ( thePredicate->GetType() == SMDSAbs_All || + thePredicate->GetType() == anElem->GetType() ) + { + long anId = anElem->GetID(); + if ( thePredicate->IsSatisfy( anId ) ) + theSequence.push_back( anId ); + } } } } void Filter::GetElementsId( const SMDS_Mesh* theMesh, - Filter::TIdSequence& theSequence ) + Filter::TIdSequence& theSequence, + SMDS_ElemIteratorPtr theElements ) { - GetElementsId(theMesh,myPredicate,theSequence); + GetElementsId(theMesh,myPredicate,theSequence,theElements); } /* @@ -3673,7 +3946,7 @@ bool ManifoldPart::process() myMapIds.Add( aFaceId ); } - if ( fi == ( myAllFacePtr.size() - 1 ) ) + if ( fi == int( myAllFacePtr.size() - 1 )) fi = 0; } // end run on vector of faces return !myMapIds.IsEmpty(); @@ -3854,24 +4127,20 @@ void ManifoldPart::expandBoundary void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink, ManifoldPart::TVectorOfFacePtr& theFaces ) const { - std::set 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.insert( aFace ); - } + SMDS_ElemIteratorPtr anItr = theLink.myNode1->GetInverseElementIterator( SMDSAbs_Face ); + SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > faces( anItr ), facesEnd; + std::set aSetOfFaces( faces, facesEnd ); + // take all faces that shared second node - anItr = theLink.myNode2->facesIterator(); + anItr = theLink.myNode2->GetInverseElementIterator( SMDSAbs_Face ); // find the common part of two sets for ( ; anItr->more(); ) { - SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next(); - if ( aSetOfFaces.count( aFace ) ) - theFaces.push_back( aFace ); + const SMDS_MeshElement* aFace = anItr->next(); + if ( aSetOfFaces.count( aFace )) + theFaces.push_back( (SMDS_MeshFace*) aFace ); } } @@ -3926,9 +4195,9 @@ SMDSAbs_ElementType BelongToMeshGroup::GetType() const return myGroup ? myGroup->GetType() : SMDSAbs_All; } -/* - ElementsOnSurface -*/ +//================================================================================ +// ElementsOnSurface +//================================================================================ ElementsOnSurface::ElementsOnSurface() { @@ -3961,8 +4230,10 @@ SMDSAbs_ElementType ElementsOnSurface::GetType() const void ElementsOnSurface::SetTolerance( const double theToler ) { if ( myToler != theToler ) - myIds.Clear(); - myToler = theToler; + { + myToler = theToler; + process(); + } } double ElementsOnSurface::GetTolerance() const @@ -4062,15 +4333,77 @@ bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) } -/* - ElementsOnShape -*/ +//================================================================================ +// ElementsOnShape +//================================================================================ -ElementsOnShape::ElementsOnShape() - : //myMesh(0), - myType(SMDSAbs_All), - myToler(Precision::Confusion()), - myAllNodesFlag(false) +namespace { + const int theIsCheckedFlag = 0x0000100; +} + +struct ElementsOnShape::Classifier +{ + Classifier() { mySolidClfr = 0; myFlags = 0; } + ~Classifier(); + void Init(const TopoDS_Shape& s, double tol, const Bnd_B3d* box = 0 ); + bool IsOut(const gp_Pnt& p) { return SetChecked( true ), (this->*myIsOutFun)( p ); } + TopAbs_ShapeEnum ShapeType() const { return myShape.ShapeType(); } + const TopoDS_Shape& Shape() const { return myShape; } + const Bnd_B3d* GetBndBox() const { return & myBox; } + double Tolerance() const { return myTol; } + bool IsChecked() { return myFlags & theIsCheckedFlag; } + bool IsSetFlag( int flag ) const { return myFlags & flag; } + void SetChecked( bool is ) { is ? SetFlag( theIsCheckedFlag ) : UnsetFlag( theIsCheckedFlag ); } + void SetFlag ( int flag ) { myFlags |= flag; } + void UnsetFlag( int flag ) { myFlags &= ~flag; } + +private: + bool isOutOfSolid (const gp_Pnt& p); + bool isOutOfBox (const gp_Pnt& p); + bool isOutOfFace (const gp_Pnt& p); + bool isOutOfEdge (const gp_Pnt& p); + bool isOutOfVertex(const gp_Pnt& p); + bool isOutOfNone (const gp_Pnt& /*p*/) { return true; } + bool isBox (const TopoDS_Shape& s); + + TopoDS_Shape prepareSolid( const TopoDS_Shape& theSolid ); + + bool (Classifier::* myIsOutFun)(const gp_Pnt& p); + BRepClass3d_SolidClassifier* mySolidClfr; // ptr because of a run-time forbidden copy-constructor + Bnd_B3d myBox; + GeomAPI_ProjectPointOnSurf myProjFace; + GeomAPI_ProjectPointOnCurve myProjEdge; + gp_Pnt myVertexXYZ; + TopoDS_Shape myShape; + double myTol; + int myFlags; +}; + +struct ElementsOnShape::OctreeClassifier : public SMESH_Octree +{ + OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers ); + OctreeClassifier( const OctreeClassifier* otherTree, + const std::vector< ElementsOnShape::Classifier >& clsOther, + std::vector< ElementsOnShape::Classifier >& cls ); + void GetClassifiersAtPoint( const gp_XYZ& p, + std::vector< ElementsOnShape::Classifier* >& classifiers ); + size_t GetSize(); + +protected: + OctreeClassifier() {} + SMESH_Octree* newChild() const { return new OctreeClassifier; } + void buildChildrenData(); + Bnd_B3d* buildRootBox(); + + std::vector< ElementsOnShape::Classifier* > myClassifiers; +}; + + +ElementsOnShape::ElementsOnShape(): + myOctree(0), + myType(SMDSAbs_All), + myToler(Precision::Confusion()), + myAllNodesFlag(false) { } @@ -4079,6 +4412,40 @@ ElementsOnShape::~ElementsOnShape() clearClassifiers(); } +Predicate* ElementsOnShape::clone() const +{ + size_t size = sizeof( *this ); + if ( myOctree ) + size += myOctree->GetSize(); + if ( !myClassifiers.empty() ) + size += sizeof( myClassifiers[0] ) * myClassifiers.size(); + if ( !myWorkClassifiers.empty() ) + size += sizeof( myWorkClassifiers[0] ) * myWorkClassifiers.size(); + if ( size > 1e+9 ) // 1G + { +#ifdef _DEBUG_ + std::cout << "Avoid ElementsOnShape::clone(), too large: " << size << " bytes " << std::endl; +#endif + return 0; + } + + ElementsOnShape* cln = new ElementsOnShape(); + cln->SetAllNodes ( myAllNodesFlag ); + cln->SetTolerance( myToler ); + cln->SetMesh ( myMeshModifTracer.GetMesh() ); + cln->myShape = myShape; // avoid creation of myClassifiers + cln->SetShape ( myShape, myType ); + cln->myClassifiers.resize( myClassifiers.size() ); + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + cln->myClassifiers[ i ].Init( BRepBuilderAPI_Copy( myClassifiers[ i ].Shape()), + myToler, myClassifiers[ i ].GetBndBox() ); + if ( myOctree ) // copy myOctree + { + cln->myOctree = new OctreeClassifier( myOctree, myClassifiers, cln->myClassifiers ); + } + return cln; +} + SMDSAbs_ElementType ElementsOnShape::GetType() const { return myType; @@ -4144,27 +4511,32 @@ void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool isOut ) void ElementsOnShape::SetShape (const TopoDS_Shape& theShape, const SMDSAbs_ElementType theType) { + bool shapeChanges = ( myShape != theShape ); myType = theType; myShape = theShape; if ( myShape.IsNull() ) 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 ( shapeChanges ) { - 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() ); - } + // find most complex shapes + 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() ); + } - clearClassifiers(); - myClassifiers.resize( shapesMap.Extent() ); - for ( int i = 0; i < shapesMap.Extent(); ++i ) - myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler ); + clearClassifiers(); + myClassifiers.resize( shapesMap.Extent() ); + for ( int i = 0; i < shapesMap.Extent(); ++i ) + myClassifiers[ i ].Init( shapesMap( i+1 ), myToler ); + } if ( theType == SMDSAbs_Node ) { @@ -4179,125 +4551,295 @@ void ElementsOnShape::SetShape (const TopoDS_Shape& theShape, void ElementsOnShape::clearClassifiers() { - for ( size_t i = 0; i < myClassifiers.size(); ++i ) - delete myClassifiers[ i ]; + // for ( size_t i = 0; i < myClassifiers.size(); ++i ) + // delete myClassifiers[ i ]; myClassifiers.clear(); + + delete myOctree; + myOctree = 0; } -bool ElementsOnShape::IsSatisfy (long elemId) +bool ElementsOnShape::IsSatisfy( long elemId ) { - const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh(); - const SMDS_MeshElement* elem = - ( myType == SMDSAbs_Node ? mesh->FindNode( elemId ) : mesh->FindElement( elemId )); - if ( !elem || myClassifiers.empty() ) + if ( myClassifiers.empty() ) + return false; + + const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh(); + if ( myType == SMDSAbs_Node ) + return IsSatisfy( mesh->FindNode( elemId )); + return IsSatisfy( mesh->FindElement( elemId )); +} + +bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem) +{ + if ( !elem ) return false; bool isSatisfy = myAllNodesFlag, isNodeOut; gp_XYZ centerXYZ (0, 0, 0); - SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator(); - while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) + if ( !myOctree && myClassifiers.size() > 5 ) { - SMESH_TNodeXYZ aPnt( aNodeItr->next() ); + myWorkClassifiers.resize( myClassifiers.size() ); + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + myWorkClassifiers[ i ] = & myClassifiers[ i ]; + myOctree = new OctreeClassifier( myWorkClassifiers ); + + SMESHUtils::FreeVector( myWorkClassifiers ); + } + + for ( int i = 0, nb = elem->NbNodes(); i < nb && (isSatisfy == myAllNodesFlag); ++i ) + { + SMESH_TNodeXYZ aPnt( elem->GetNode( i )); centerXYZ += aPnt; isNodeOut = true; if ( !getNodeIsOut( aPnt._node, isNodeOut )) { - for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i ) - isNodeOut = myClassifiers[i]->IsOut( aPnt ); + if ( myOctree ) + { + myWorkClassifiers.clear(); + myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers ); + + for ( size_t i = 0; i < myWorkClassifiers.size(); ++i ) + myWorkClassifiers[i]->SetChecked( false ); + for ( size_t i = 0; i < myWorkClassifiers.size() && isNodeOut; ++i ) + if ( !myWorkClassifiers[i]->IsChecked() ) + isNodeOut = myWorkClassifiers[i]->IsOut( aPnt ); + } + else + { + for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i ) + isNodeOut = myClassifiers[i].IsOut( aPnt ); + } setNodeIsOut( aPnt._node, isNodeOut ); } isSatisfy = !isNodeOut; } // Check the center point for volumes MantisBug 0020168 - if (isSatisfy && - myAllNodesFlag && - myClassifiers[0]->ShapeType() == TopAbs_SOLID) + 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 ); + if ( myOctree ) + { + myWorkClassifiers.clear(); + myOctree->GetClassifiersAtPoint( centerXYZ, myWorkClassifiers ); + for ( size_t i = 0; i < myWorkClassifiers.size() && !isSatisfy; ++i ) + isSatisfy = ! myWorkClassifiers[i]->IsOut( centerXYZ ); + } + else + { + 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(); -} +//================================================================================ +/*! + * \brief Check and optionally return a satisfying shape + */ +//================================================================================ -bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p) +bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node, + TopoDS_Shape* okShape) { - return (this->*myIsOutFun)( p ); + if ( !node ) + return false; + + if ( !myOctree && myClassifiers.size() > 5 ) + { + myWorkClassifiers.resize( myClassifiers.size() ); + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + myWorkClassifiers[ i ] = & myClassifiers[ i ]; + myOctree = new OctreeClassifier( myWorkClassifiers ); + } + + bool isNodeOut = true; + + if ( okShape || !getNodeIsOut( node, isNodeOut )) + { + SMESH_NodeXYZ aPnt = node; + if ( myOctree ) + { + myWorkClassifiers.clear(); + myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers ); + + for ( size_t i = 0; i < myWorkClassifiers.size(); ++i ) + myWorkClassifiers[i]->SetChecked( false ); + + for ( size_t i = 0; i < myWorkClassifiers.size(); ++i ) + if ( !myWorkClassifiers[i]->IsChecked() && + !myWorkClassifiers[i]->IsOut( aPnt )) + { + isNodeOut = false; + if ( okShape ) + *okShape = myWorkClassifiers[i]->Shape(); + break; + } + } + else + { + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + if ( !myClassifiers[i].IsOut( aPnt )) + { + isNodeOut = false; + if ( okShape ) + *okShape = myClassifiers[i].Shape(); + break; + } + } + setNodeIsOut( node, isNodeOut ); + } + + return !isNodeOut; } -void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double theTol) +void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape, + double theTol, + const Bnd_B3d* theBox ) { myShape = theShape; myTol = theTol; + myFlags = 0; + + bool isShapeBox = false; switch ( myShape.ShapeType() ) { - case TopAbs_SOLID: { - if ( isBox( theShape )) + case TopAbs_SOLID: + { + if (( isShapeBox = isBox( theShape ))) { - myIsOutFun = & ElementsOnShape::TClassifier::isOutOfBox; + myIsOutFun = & ElementsOnShape::Classifier::isOutOfBox; } else { - mySolidClfr.Load(theShape); - myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid; + mySolidClfr = new BRepClass3d_SolidClassifier( prepareSolid( theShape )); + myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid; } break; } - case TopAbs_FACE: { + 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; + if ( surf.IsNull() ) + myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone; + else + { + surf->Bounds( u1,u2,v1,v2 ); + myProjFace.Init(surf, u1,u2, v1,v2, myTol ); + myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace; + } break; } - case TopAbs_EDGE: { + 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; + Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2); + if ( curve.IsNull() ) + myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone; + else + { + myProjEdge.Init(curve, u1, u2); + myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge; + } break; } - case TopAbs_VERTEX:{ + case TopAbs_VERTEX: + { myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) ); - myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex; + myIsOutFun = & ElementsOnShape::Classifier::isOutOfVertex; break; } default: - throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier"); + throw SALOME_Exception("Programmer error in usage of ElementsOnShape::Classifier"); + } + + if ( !isShapeBox ) + { + if ( theBox ) + { + myBox = *theBox; + } + else + { + Bnd_Box box; + if ( myShape.ShapeType() == TopAbs_FACE ) + { + BRepAdaptor_Surface SA( TopoDS::Face( myShape ), /*useBoundaries=*/false ); + if ( SA.GetType() == GeomAbs_BSplineSurface ) + BRepBndLib::AddOptimal( myShape, box, + /*useTriangulation=*/true, /*useShapeTolerance=*/true ); + } + if ( box.IsVoid() ) + BRepBndLib::Add( myShape, box ); + myBox.Clear(); + myBox.Add( box.CornerMin() ); + myBox.Add( box.CornerMax() ); + gp_XYZ halfSize = 0.5 * ( box.CornerMax().XYZ() - box.CornerMin().XYZ() ); + for ( int iDim = 1; iDim <= 3; ++iDim ) + { + double x = halfSize.Coord( iDim ); + halfSize.SetCoord( iDim, x + Max( myTol, 1e-2 * x )); + } + myBox.SetHSize( halfSize ); + } } } -bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p) +ElementsOnShape::Classifier::~Classifier() +{ + delete mySolidClfr; mySolidClfr = 0; +} + +TopoDS_Shape ElementsOnShape::Classifier::prepareSolid( const TopoDS_Shape& theSolid ) { - mySolidClfr.Perform( p, myTol ); - return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON ); + // try to limit tolerance of theSolid down to myTol (issue #19026) + + // check if tolerance of theSolid is more than myTol + bool tolIsOk = true; // max tolerance is at VERTEXes + for ( TopExp_Explorer exp( theSolid, TopAbs_VERTEX ); exp.More() && tolIsOk; exp.Next() ) + tolIsOk = ( myTol >= BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current() ))); + if ( tolIsOk ) + return theSolid; + + // make a copy to prevent the original shape from changes + TopoDS_Shape resultShape = BRepBuilderAPI_Copy( theSolid ); + + if ( !GEOMUtils::FixShapeTolerance( resultShape, TopAbs_SHAPE, myTol )) + return theSolid; + return resultShape; } -bool ElementsOnShape::TClassifier::isOutOfBox (const gp_Pnt& p) +bool ElementsOnShape::Classifier::isOutOfSolid( const gp_Pnt& p ) +{ + if ( isOutOfBox( p )) return true; + mySolidClfr->Perform( p, myTol ); + return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON ); +} + +bool ElementsOnShape::Classifier::isOutOfBox( const gp_Pnt& p ) { return myBox.IsOut( p.XYZ() ); } -bool ElementsOnShape::TClassifier::isOutOfFace (const gp_Pnt& p) +bool ElementsOnShape::Classifier::isOutOfFace( const gp_Pnt& p ) { + if ( isOutOfBox( p )) return true; myProjFace.Perform( p ); if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol ) { // check relatively to the face - Quantity_Parameter u, v; + Standard_Real u, v; myProjFace.LowerDistanceParameters(u, v); gp_Pnt2d aProjPnt (u, v); BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol ); @@ -4307,18 +4849,19 @@ bool ElementsOnShape::TClassifier::isOutOfFace (const gp_Pnt& p) return true; } -bool ElementsOnShape::TClassifier::isOutOfEdge (const gp_Pnt& p) +bool ElementsOnShape::Classifier::isOutOfEdge( const gp_Pnt& p ) { + if ( isOutOfBox( p )) return true; myProjEdge.Perform( p ); return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol ); } -bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p) +bool ElementsOnShape::Classifier::isOutOfVertex( const gp_Pnt& p ) { return ( myVertexXYZ.Distance( p ) > myTol ); } -bool ElementsOnShape::TClassifier::isBox (const TopoDS_Shape& theShape) +bool ElementsOnShape::Classifier::isBox(const TopoDS_Shape& theShape ) { TopTools_IndexedMapOfShape vMap; TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap ); @@ -4345,6 +4888,132 @@ bool ElementsOnShape::TClassifier::isBox (const TopoDS_Shape& theShape) return true; } +ElementsOnShape:: +OctreeClassifier::OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers ) + :SMESH_Octree( new SMESH_TreeLimit ) +{ + myClassifiers = classifiers; + compute(); +} + +ElementsOnShape:: +OctreeClassifier::OctreeClassifier( const OctreeClassifier* otherTree, + const std::vector< ElementsOnShape::Classifier >& clsOther, + std::vector< ElementsOnShape::Classifier >& cls ) + :SMESH_Octree( new SMESH_TreeLimit ) +{ + myBox = new Bnd_B3d( *otherTree->getBox() ); + + if (( myIsLeaf = otherTree->isLeaf() )) + { + myClassifiers.resize( otherTree->myClassifiers.size() ); + for ( size_t i = 0; i < otherTree->myClassifiers.size(); ++i ) + { + int ind = otherTree->myClassifiers[i] - & clsOther[0]; + myClassifiers[ i ] = & cls[ ind ]; + } + } + else if ( otherTree->myChildren ) + { + myChildren = new SMESH_Tree< Bnd_B3d, 8 > * [ 8 ]; + for ( int i = 0; i < nbChildren(); i++ ) + myChildren[i] = + new OctreeClassifier( static_cast( otherTree->myChildren[i]), + clsOther, cls ); + } +} + +void ElementsOnShape:: +OctreeClassifier::GetClassifiersAtPoint( const gp_XYZ& point, + std::vector< ElementsOnShape::Classifier* >& result ) +{ + if ( getBox()->IsOut( point )) + return; + + if ( isLeaf() ) + { + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + if ( !myClassifiers[i]->GetBndBox()->IsOut( point )) + result.push_back( myClassifiers[i] ); + } + else + { + for (int i = 0; i < nbChildren(); i++) + ((OctreeClassifier*) myChildren[i])->GetClassifiersAtPoint( point, result ); + } +} + +size_t ElementsOnShape::OctreeClassifier::GetSize() +{ + size_t res = sizeof( *this ); + if ( !myClassifiers.empty() ) + res += sizeof( myClassifiers[0] ) * myClassifiers.size(); + + if ( !isLeaf() ) + for (int i = 0; i < nbChildren(); i++) + res += ((OctreeClassifier*) myChildren[i])->GetSize(); + + return res; +} + +void ElementsOnShape::OctreeClassifier::buildChildrenData() +{ + // distribute myClassifiers among myChildren + + const int childFlag[8] = { 0x0000001, + 0x0000002, + 0x0000004, + 0x0000008, + 0x0000010, + 0x0000020, + 0x0000040, + 0x0000080 }; + int nbInChild[8] = { 0, 0, 0, 0, 0, 0, 0, 0 }; + + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + { + for ( int j = 0; j < nbChildren(); j++ ) + { + if ( !myClassifiers[i]->GetBndBox()->IsOut( *myChildren[j]->getBox() )) + { + myClassifiers[i]->SetFlag( childFlag[ j ]); + ++nbInChild[ j ]; + } + } + } + + for ( int j = 0; j < nbChildren(); j++ ) + { + OctreeClassifier* child = static_cast( myChildren[ j ]); + child->myClassifiers.resize( nbInChild[ j ]); + for ( size_t i = 0; nbInChild[ j ] && i < myClassifiers.size(); ++i ) + { + if ( myClassifiers[ i ]->IsSetFlag( childFlag[ j ])) + { + --nbInChild[ j ]; + child->myClassifiers[ nbInChild[ j ]] = myClassifiers[ i ]; + myClassifiers[ i ]->UnsetFlag( childFlag[ j ]); + } + } + } + SMESHUtils::FreeVector( myClassifiers ); + + // define if a child isLeaf() + for ( int i = 0; i < nbChildren(); i++ ) + { + OctreeClassifier* child = static_cast( myChildren[ i ]); + child->myIsLeaf = ( child->myClassifiers.size() <= 5 || + child->maxSize() < child->myClassifiers[0]->Tolerance() ); + } +} + +Bnd_B3d* ElementsOnShape::OctreeClassifier::buildRootBox() +{ + Bnd_B3d* box = new Bnd_B3d; + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + box->Add( *myClassifiers[i]->GetBndBox() ); + return box; +} /* Class : BelongToGeom @@ -4354,25 +5023,45 @@ bool ElementsOnShape::TClassifier::isBox (const TopoDS_Shape& theShape) BelongToGeom::BelongToGeom() : myMeshDS(NULL), - myType(SMDSAbs_All), + myType(SMDSAbs_NbElementTypes), myIsSubshape(false), myTolerance(Precision::Confusion()) {} +Predicate* BelongToGeom::clone() const +{ + BelongToGeom* cln = 0; + if ( myElementsOnShapePtr ) + if ( ElementsOnShape* eos = static_cast( myElementsOnShapePtr->clone() )) + { + cln = new BelongToGeom( *this ); + cln->myElementsOnShapePtr.reset( eos ); + } + return cln; +} + void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh ) { - myMeshDS = dynamic_cast(theMesh); - init(); + if ( myMeshDS != theMesh ) + { + myMeshDS = dynamic_cast(theMesh); + init(); + } + if ( myElementsOnShapePtr ) + myElementsOnShapePtr->SetMesh( myMeshDS ); } void BelongToGeom::SetGeom( const TopoDS_Shape& theShape ) { - myShape = theShape; - init(); + if ( myShape != theShape ) + { + myShape = theShape; + init(); + } } static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap, - const TopoDS_Shape& theShape) + const TopoDS_Shape& theShape) { if (theMap.Contains(theShape)) return true; @@ -4394,7 +5083,7 @@ static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap, void BelongToGeom::init() { - if (!myMeshDS || myShape.IsNull()) return; + if ( !myMeshDS || myShape.IsNull() ) return; // is sub-shape of main shape? TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh(); @@ -4403,38 +5092,31 @@ void BelongToGeom::init() } else { TopTools_IndexedMapOfShape aMap; - TopExp::MapShapes(aMainShape, aMap); - myIsSubshape = IsSubShape(aMap, myShape); + TopExp::MapShapes( aMainShape, aMap ); + myIsSubshape = IsSubShape( aMap, myShape ); + if ( myIsSubshape ) + { + aMap.Clear(); + TopExp::MapShapes( myShape, aMap ); + mySubShapesIDs.Clear(); + for ( int i = 1; i <= aMap.Extent(); ++i ) + { + int subID = myMeshDS->ShapeToIndex( aMap( i )); + if ( subID > 0 ) + mySubShapesIDs.Add( subID ); + } + } } //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(); + if ( !myElementsOnShapePtr ) + myElementsOnShapePtr.reset( new ElementsOnShape() ); + myElementsOnShapePtr->SetTolerance( myTolerance ); + myElementsOnShapePtr->SetAllNodes( true ); // "belong", while false means "lays on" + myElementsOnShapePtr->SetMesh( myMeshDS ); + myElementsOnShapePtr->SetShape( myShape, myType ); } - return false; } bool BelongToGeom::IsSatisfy (long theId) @@ -4447,49 +5129,28 @@ bool BelongToGeom::IsSatisfy (long theId) return myElementsOnShapePtr->IsSatisfy(theId); } - // Case of submesh + // Case of sub-mesh + if (myType == SMDSAbs_Node) { - if( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ) ) + 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 + return mySubShapesIDs.Contains( aNode->getshapeId() ); } } else { if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId )) { - if ( anElem->getshapeId() < 1 ) - return myElementsOnShapePtr->IsSatisfy(theId); - - if( myType == SMDSAbs_All ) + if ( myType == SMDSAbs_All || anElem->GetType() == myType ) { - 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 )); - } + if ( anElem->getshapeId() < 1 ) + return myElementsOnShapePtr->IsSatisfy(theId); + else + return mySubShapesIDs.Contains( anElem->getshapeId() ); } } } @@ -4499,8 +5160,11 @@ bool BelongToGeom::IsSatisfy (long theId) void BelongToGeom::SetType (SMDSAbs_ElementType theType) { - myType = theType; - init(); + if ( myType != theType ) + { + myType = theType; + init(); + } } SMDSAbs_ElementType BelongToGeom::GetType() const @@ -4521,8 +5185,7 @@ const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const void BelongToGeom::SetTolerance (double theTolerance) { myTolerance = theTolerance; - if (!myIsSubshape) - init(); + init(); } double BelongToGeom::GetTolerance() @@ -4533,26 +5196,46 @@ double BelongToGeom::GetTolerance() /* Class : LyingOnGeom Description : Predicate for verifying whether entiy lying or partially lying on - specified geometrical support + specified geometrical support */ LyingOnGeom::LyingOnGeom() : myMeshDS(NULL), - myType(SMDSAbs_All), + myType(SMDSAbs_NbElementTypes), myIsSubshape(false), myTolerance(Precision::Confusion()) {} +Predicate* LyingOnGeom::clone() const +{ + LyingOnGeom* cln = 0; + if ( myElementsOnShapePtr ) + if ( ElementsOnShape* eos = static_cast( myElementsOnShapePtr->clone() )) + { + cln = new LyingOnGeom( *this ); + cln->myElementsOnShapePtr.reset( eos ); + } + return cln; +} + void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh ) { - myMeshDS = dynamic_cast(theMesh); - init(); + if ( myMeshDS != theMesh ) + { + myMeshDS = dynamic_cast(theMesh); + init(); + } + if ( myElementsOnShapePtr ) + myElementsOnShapePtr->SetMesh( myMeshDS ); } void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape ) { - myShape = theShape; - init(); + if ( myShape != theShape ) + { + myShape = theShape; + init(); + } } void LyingOnGeom::init() @@ -4580,13 +5263,14 @@ void LyingOnGeom::init() mySubShapesIDs.Add( subID ); } } - else + // else // to be always ready to check an element not bound to geometry { - myElementsOnShapePtr.reset(new ElementsOnShape()); - myElementsOnShapePtr->SetTolerance(myTolerance); - myElementsOnShapePtr->SetAllNodes(false); // lays on, while true means "belong" - myElementsOnShapePtr->SetMesh(myMeshDS); - myElementsOnShapePtr->SetShape(myShape, myType); + if ( !myElementsOnShapePtr ) + myElementsOnShapePtr.reset( new ElementsOnShape() ); + myElementsOnShapePtr->SetTolerance( myTolerance ); + myElementsOnShapePtr->SetAllNodes( false ); // lays on, while true means "belong" + myElementsOnShapePtr->SetMesh( myMeshDS ); + myElementsOnShapePtr->SetShape( myShape, myType ); } } @@ -4608,7 +5292,8 @@ bool LyingOnGeom::IsSatisfy( long theId ) if ( mySubShapesIDs.Contains( elem->getshapeId() )) return true; - if ( elem->GetType() != SMDSAbs_Node ) + if (( elem->GetType() != SMDSAbs_Node ) && + ( myType == SMDSAbs_All || elem->GetType() == myType )) { SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator(); while ( nodeItr->more() ) @@ -4624,8 +5309,11 @@ bool LyingOnGeom::IsSatisfy( long theId ) void LyingOnGeom::SetType( SMDSAbs_ElementType theType ) { - myType = theType; - init(); + if ( myType != theType ) + { + myType = theType; + init(); + } } SMDSAbs_ElementType LyingOnGeom::GetType() const @@ -4646,8 +5334,7 @@ const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const void LyingOnGeom::SetTolerance (double theTolerance) { myTolerance = theTolerance; - if (!myIsSubshape) - init(); + init(); } double LyingOnGeom::GetTolerance() @@ -4655,39 +5342,6 @@ 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_MapOfShape aSubShapes; - // TopExp_Explorer exp( theShape, theFindShapeEnum, theAvoidShapeEnum ); - // for ( ; exp.More(); exp.Next() ) - // { - // const TopoDS_Shape& aShape = exp.Current(); - // if ( !aSubShapes.Add( aShape )) continue; - - // if ( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape )) - // { - // if ( aSubMesh->Contains( theElem )) - // return true; - - // SMDS_ElemIteratorPtr nodeItr = theElem->nodesIterator(); - // while ( nodeItr->more() ) - // { - // const SMDS_MeshElement* aNode = nodeItr->next(); - // if ( aSubMesh->Contains( aNode )) - // return true; - // } - // } - // } - return false; -} - TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0) {}