X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=bf8c097c0037c7ec85b048bd7f02d66fc3f6c5ff;hp=769096df4344b7e69a9d61ac59e8d7d6d903b5f9;hb=1746a461949c030eb46ccb860e586ade2e5de5e3;hpb=193c49c87753b6ccabb2b5e6dc935aa480d2d43e diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 769096df4..bf8c097c0 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-2016 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 @@ -39,11 +39,14 @@ #include #include +#include +#include #include #include #include #include #include +#include #include #include #include @@ -174,7 +177,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; @@ -311,12 +314,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 || @@ -335,7 +338,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 )); } @@ -521,148 +524,164 @@ 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: { // 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 ) { @@ -701,7 +720,7 @@ double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) aMin = getAngle(P( P.size() ), P( 1 ), P( 2 )); aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 ))); - for ( int i = 2; i < P.size(); i++ ) + for ( size_t i = 2; i < P.size(); i++ ) { double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) ); aMin = Min(aMin,A0); @@ -1476,10 +1495,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); } @@ -1815,10 +1834,10 @@ void Length2D::GetValues(TValues& theValues) dynamic_cast(anElem); // use special nodes iterator SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator(); - long aNodeId[4]; + long aNodeId[4] = { 0,0,0,0 }; gp_Pnt P[4]; - double aLength; + double aLength = 0; const SMDS_MeshElement* aNode; if(anIter->more()){ aNode = anIter->next(); @@ -1852,7 +1871,7 @@ void Length2D::GetValues(TValues& theValues) } else { SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator(); - long aNodeId[2]; + long aNodeId[2] = {0,0}; gp_Pnt P[3]; double aLength; @@ -1940,7 +1959,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++) { @@ -2021,7 +2040,7 @@ void MultiConnection2D::GetValues(MValues& theValues) (anElem)->interlacedNodesElemIterator(); else aNodesIter = anElem->nodesIterator(); - long aNodeId[3]; + long aNodeId[3] = {0,0,0}; //int aNbConnects=0; const SMDS_MeshNode* aNode0; @@ -2097,6 +2116,42 @@ SMDSAbs_ElementType BallDiameter::GetType() const return SMDSAbs_Ball; } +//================================================================================ +/* + Class : NodeConnectivityNumber + Description : Functor returning number of elements connected to a node +*/ +//================================================================================ + +double NodeConnectivityNumber::GetValue( long theId ) +{ + double nb = 0; + + if ( const SMDS_MeshNode* node = myMesh->FindNode( theId )) + { + SMDSAbs_ElementType type; + if ( myMesh->NbVolumes() > 0 ) + type = SMDSAbs_Volume; + else if ( myMesh->NbFaces() > 0 ) + type = SMDSAbs_Face; + else if ( myMesh->NbEdges() > 0 ) + type = SMDSAbs_Edge; + else + return 0; + nb = node->NbInverseElements( type ); + } + return nb; +} + +double NodeConnectivityNumber::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; +} + +SMDSAbs_ElementType NodeConnectivityNumber::GetType() const +{ + return SMDSAbs_Node; +} /* PREDICATES @@ -2146,7 +2201,7 @@ bool BareBorderVolume::IsSatisfy(long theElementId ) if ( myTool.IsFreeFace( iF )) { const SMDS_MeshNode** n = myTool.GetFaceNodes(iF); - vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF)); + std::vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF)); if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false)) return true; } @@ -2285,15 +2340,15 @@ void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh ) while ( nIt->more() ) nodesToCheck.insert( nodesToCheck.end(), nIt->next() ); - list< list< const SMDS_MeshNode*> > nodeGroups; + std::list< std::list< const SMDS_MeshNode*> > nodeGroups; SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler ); myCoincidentIDs.Clear(); - list< list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin(); + std::list< std::list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin(); for ( ; groupIt != nodeGroups.end(); ++groupIt ) { - list< const SMDS_MeshNode*>& coincNodes = *groupIt; - list< const SMDS_MeshNode*>::iterator n = coincNodes.begin(); + std::list< const SMDS_MeshNode*>& coincNodes = *groupIt; + std::list< const SMDS_MeshNode*>::iterator n = coincNodes.begin(); for ( ; n != coincNodes.end(); ++n ) myCoincidentIDs.Add( (*n)->GetID() ); } @@ -2325,7 +2380,7 @@ bool CoincidentElements::IsSatisfy( long theElementId ) if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId )) { if ( e->GetType() != GetType() ) return false; - set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() ); + std::set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() ); const int nbNodes = e->NbNodes(); SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() ); while ( invIt->more() ) @@ -2494,7 +2549,7 @@ void FreeEdges::GetBoreders(TBorders& theBorders) interlacedNodesElemIterator(); else aNodesIter = anElem->nodesIterator(); - long aNodeId[2]; + long aNodeId[2] = {0,0}; const SMDS_MeshElement* aNode; if(aNodesIter->more()){ aNode = aNodesIter->next(); @@ -2572,18 +2627,20 @@ bool FreeFaces::IsSatisfy( long theId ) int nbNode = aFace->NbNodes(); // collect volumes to check that number of volumes with count equal nbNode not less than 2 - typedef map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters - typedef map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator + typedef std::map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters + typedef std::map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator TMapOfVolume mapOfVol; SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator(); - while ( nodeItr->more() ) { + while ( nodeItr->more() ) + { const SMDS_MeshNode* aNode = static_cast(nodeItr->next()); if ( !aNode ) continue; SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume); - while ( volItr->more() ) { + while ( volItr->more() ) + { SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next(); - TItrMapOfVolume itr = mapOfVol.insert(make_pair(aVol, 0)).first; + TItrMapOfVolume itr = mapOfVol.insert( std::make_pair( aVol, 0 )).first; (*itr).second++; } } @@ -2687,8 +2744,8 @@ void GroupColor::SetMesh( const SMDS_Mesh* theMesh ) return; // iterates on groups and find necessary elements ids - const std::set& aGroups = aMesh->GetGroups(); - set::const_iterator GrIt = aGroups.begin(); + const std::set& aGroups = aMesh->GetGroups(); + std::set::const_iterator GrIt = aGroups.begin(); for (; GrIt != aGroups.end(); GrIt++) { SMESHDS_GroupBase* aGrp = (*GrIt); @@ -2915,10 +2972,10 @@ void ConnectedElements::SetPoint( double x, double y, double z ) // find myNodeID by myXYZ if possible if ( myMeshModifTracer.GetMesh() ) { - auto_ptr searcher + SMESHUtils::Deleter searcher ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() )); - vector< const SMDS_MeshElement* > foundElems; + std::vector< const SMDS_MeshElement* > foundElems; searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems ); if ( !foundElems.empty() ) @@ -2944,7 +3001,7 @@ bool ConnectedElements::IsSatisfy( long theElementId ) if ( !node0 ) return false; - list< const SMDS_MeshNode* > nodeQueue( 1, node0 ); + std::list< const SMDS_MeshNode* > nodeQueue( 1, node0 ); std::set< int > checkedNodeIDs; // algo: // foreach node in nodeQueue: @@ -2994,6 +3051,17 @@ bool ConnectedElements::IsSatisfy( long theElementId ) */ //================================================================================ +namespace +{ + inline bool isLessAngle( const gp_Vec& v1, const gp_Vec& v2, const double cos ) + { + double dot = v1 * v2; // cos * |v1| * |v2| + double l1 = v1.SquareMagnitude(); + double l2 = v2.SquareMagnitude(); + return (( dot * cos >= 0 ) && + ( dot * dot ) / l1 / l2 >= ( cos * cos )); + } +} CoplanarFaces::CoplanarFaces() : myFaceID(0), myToler(0) { @@ -3005,7 +3073,7 @@ void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh ) { // Build a set of coplanar face ids - myCoplanarIDs.clear(); + myCoplanarIDs.Clear(); if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler ) return; @@ -3019,11 +3087,11 @@ void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh ) if (!normOK) return; - const double radianTol = myToler * M_PI / 180.; - std::set< SMESH_TLink > checkedLinks; + const double cosTol = Cos( myToler * M_PI / 180. ); + NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks; - std::list< pair< const SMDS_MeshElement*, gp_Vec > > faceQueue; - faceQueue.push_back( make_pair( face, myNorm )); + std::list< std::pair< const SMDS_MeshElement*, gp_Vec > > faceQueue; + faceQueue.push_back( std::make_pair( face, myNorm )); while ( !faceQueue.empty() ) { face = faceQueue.front().first; @@ -3034,7 +3102,7 @@ void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh ) { const SMDS_MeshNode* n1 = face->GetNode( i ); const SMDS_MeshNode* n2 = face->GetNode(( i+1 )%nbN); - if ( !checkedLinks.insert( SMESH_TLink( n1, n2 )).second ) + if ( !checkedLinks.Add( SMESH_TLink( n1, n2 ))) continue; SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face); while ( fIt->more() ) @@ -3043,10 +3111,10 @@ void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh ) if ( f->GetNodeIndex( n2 ) > -1 ) { gp_Vec norm = getNormale( static_cast(f), &normOK ); - if (!normOK || myNorm.Angle( norm ) <= radianTol) + if (!normOK || isLessAngle( myNorm, norm, cosTol)) { - myCoplanarIDs.insert( f->GetID() ); - faceQueue.push_back( make_pair( f, norm )); + myCoplanarIDs.Add( f->GetID() ); + faceQueue.push_back( std::make_pair( f, norm )); } } } @@ -3056,7 +3124,7 @@ void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh ) } bool CoplanarFaces::IsSatisfy( long theElementId ) { - return myCoplanarIDs.count( theElementId ); + return myCoplanarIDs.Contains( theElementId ); } /* @@ -3184,14 +3252,13 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr ) myIds.Clear(); TCollection_AsciiString aStr = theStr; - //aStr.RemoveAll( ' ' ); - //aStr.RemoveAll( '\t' ); for ( int i = 1; i <= aStr.Length(); ++i ) - if ( isspace( aStr.Value( i ))) + { + char c = aStr.Value( i ); + if ( !isdigit( c ) && c != ',' && c != '-' ) aStr.SetValue( i, ','); - - for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) ) - aStr.Remove( aPos, 1 ); + } + aStr.RemoveAll( ' ' ); TCollection_AsciiString tmpStr = aStr.Token( ",", 1 ); int i = 1; @@ -3667,7 +3734,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(); @@ -3920,9 +3987,9 @@ SMDSAbs_ElementType BelongToMeshGroup::GetType() const return myGroup ? myGroup->GetType() : SMDSAbs_All; } -/* - ElementsOnSurface -*/ +//================================================================================ +// ElementsOnSurface +//================================================================================ ElementsOnSurface::ElementsOnSurface() { @@ -4056,15 +4123,71 @@ 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; } + 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 isBox (const TopoDS_Shape& s); + + 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 ); +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) { } @@ -4073,6 +4196,25 @@ ElementsOnShape::~ElementsOnShape() clearClassifiers(); } +Predicate* ElementsOnShape::clone() const +{ + 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; @@ -4138,27 +4280,31 @@ 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() ); - } + 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 ) { @@ -4173,12 +4319,15 @@ 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 = @@ -4190,6 +4339,14 @@ bool ElementsOnShape::IsSatisfy (long elemId) gp_XYZ centerXYZ (0, 0, 0); + 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 ); + } + SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator(); while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) { @@ -4199,93 +4356,138 @@ bool ElementsOnShape::IsSatisfy (long elemId) 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 ) + 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(); -} - -bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p) -{ - return (this->*myIsOutFun)( p ); -} - -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(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; + 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); + Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2); myProjEdge.Init(curve, u1, u2); - myIsOutFun = & ElementsOnShape::TClassifier::isOutOfEdge; + 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; + 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 ); + } + } +} + +ElementsOnShape::Classifier::~Classifier() +{ + delete mySolidClfr; mySolidClfr = 0; } -bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p) +bool ElementsOnShape::Classifier::isOutOfSolid (const gp_Pnt& p) { - mySolidClfr.Perform( p, myTol ); - return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON ); + mySolidClfr->Perform( p, myTol ); + return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON ); } -bool ElementsOnShape::TClassifier::isOutOfBox (const gp_Pnt& p) +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) { myProjFace.Perform( p ); if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol ) @@ -4301,18 +4503,18 @@ 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) { 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 ); @@ -4339,6 +4541,118 @@ 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 ); + } +} + +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 ); + } +} + +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 @@ -4348,25 +4662,38 @@ 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 = new BelongToGeom( *this ); + cln->myElementsOnShapePtr.reset( static_cast( myElementsOnShapePtr->clone() )); + return cln; +} + void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh ) { - myMeshDS = dynamic_cast(theMesh); - init(); + if ( myMeshDS != theMesh ) + { + myMeshDS = dynamic_cast(theMesh); + init(); + } } 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; @@ -4388,7 +4715,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(); @@ -4397,38 +4724,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) @@ -4441,49 +4761,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 ) - { - 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() ) + if ( anElem->GetType() == myType ) { - 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() ); } } } @@ -4493,8 +4792,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 @@ -4515,8 +4817,7 @@ const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const void BelongToGeom::SetTolerance (double theTolerance) { myTolerance = theTolerance; - if (!myIsSubshape) - init(); + init(); } double BelongToGeom::GetTolerance() @@ -4527,26 +4828,39 @@ 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 = new LyingOnGeom( *this ); + cln->myElementsOnShapePtr.reset( static_cast( myElementsOnShapePtr->clone() )); + return cln; +} + void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh ) { - myMeshDS = dynamic_cast(theMesh); - init(); + if ( myMeshDS != theMesh ) + { + myMeshDS = dynamic_cast(theMesh); + init(); + } } void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape ) { - myShape = theShape; - init(); + if ( myShape != theShape ) + { + myShape = theShape; + init(); + } } void LyingOnGeom::init() @@ -4574,13 +4888,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 ); } } @@ -4602,7 +4917,7 @@ bool LyingOnGeom::IsSatisfy( long theId ) if ( mySubShapesIDs.Contains( elem->getshapeId() )) return true; - if ( elem->GetType() != SMDSAbs_Node ) + if ( elem->GetType() != SMDSAbs_Node && elem->GetType() == myType ) { SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator(); while ( nodeItr->more() ) @@ -4618,8 +4933,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 @@ -4640,8 +4958,7 @@ const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const void LyingOnGeom::SetTolerance (double theTolerance) { myTolerance = theTolerance; - if (!myIsSubshape) - init(); + init(); } double LyingOnGeom::GetTolerance() @@ -4649,39 +4966,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) {}