X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=0474d7a57c31299632904359f51441550daea167;hp=13a2a0008a8ecfa2845b68459e4467f1b35775dc;hb=HEAD;hpb=90e0893c89027a434e832eb72eeac85391b909e9 diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 13a2a0008..578557775 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2024 CEA, EDF, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -28,15 +28,15 @@ #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_MeshAlgos.hxx" #include "SMESH_OctreeNode.hxx" +#include "SMESH_Comment.hxx" +#include #include #include @@ -104,7 +104,7 @@ namespace { gp_Vec v1( P1 - P2 ), v2( P3 - P2 ); double dot = v1 * v2, len1 = v1.SquareMagnitude(), len2 = v2.SquareMagnitude(); - return ( len1 < gp::Resolution() || len2 < gp::Resolution() ? -1 : + return ( dot < 0 || len1 < gp::Resolution() || len2 < gp::Resolution() ? -1 : dot * dot / len1 / len2 ); } @@ -128,7 +128,7 @@ namespace { return aDist; } - int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId ) + int getNbMultiConnection( const SMDS_Mesh* theMesh, const smIdType theId ) { if ( theMesh == 0 ) return 0; @@ -144,13 +144,13 @@ namespace { // Case 1 Case 2 // | | | | | // | | | | | - // +-----+------+ +-----+------+ + // +-----+------+ +-----+------+ // | | | | // | | | | // 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 ); @@ -226,7 +226,7 @@ void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh ) myMesh = theMesh; } -bool NumericalFunctor::GetPoints(const int theId, +bool NumericalFunctor::GetPoints(const smIdType theId, TSequenceOfXYZ& theRes ) const { theRes.clear(); @@ -235,7 +235,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 ); @@ -253,26 +253,7 @@ 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 ) { SMESH_NodeXYZ p; while( anIter->more() ) { @@ -313,6 +294,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 @@ -324,12 +323,12 @@ double NumericalFunctor::Round( const double & aVal ) */ //================================================================================ -void NumericalFunctor::GetHistogram(int nbIntervals, - std::vector& nbEvents, - std::vector& funValues, - const std::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 || @@ -348,7 +347,7 @@ void NumericalFunctor::GetHistogram(int nbIntervals, } else { - std::vector::const_iterator id = elements.begin(); + std::vector::const_iterator id = elements.begin(); for ( ; id != elements.end(); ++id ) values.insert( GetValue( *id )); } @@ -494,7 +493,7 @@ double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P ) // } // } // { // polygons - + // } if( myPrecision >= 0 ) @@ -736,7 +735,7 @@ double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) double A0 = getCos2( P( i-1 ), P( i ), P( i+1 ) ); aMaxCos2 = Max( aMaxCos2, A0 ); } - if ( aMaxCos2 <= 0 ) + if ( aMaxCos2 < 0 ) return 0; // all nodes coincide double cos = sqrt( aMaxCos2 ); @@ -768,19 +767,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; } @@ -853,7 +842,7 @@ 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 @@ -864,10 +853,10 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) 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. ); + 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 ] ) ) ); @@ -897,24 +886,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; @@ -922,6 +911,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] @@ -1004,6 +998,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 ) @@ -1015,8 +1105,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 @@ -1028,6 +1118,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; @@ -1035,12 +1130,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; } @@ -1084,7 +1179,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 )}; @@ -1103,7 +1198,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 )}; @@ -1128,9 +1223,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 )}; @@ -1265,7 +1364,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 )}; @@ -1279,7 +1378,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 ); @@ -1288,7 +1387,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 )); } @@ -1317,23 +1416,14 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const */ //================================================================================ -double Warping::GetValue( const TSequenceOfXYZ& P ) +bool Warping::IsApplicable( const SMDS_MeshElement* element ) const { - if ( P.size() != 4 ) - return 0; - - gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.; - - double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G ); - double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G ); - double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G ); - double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G ); - - double val = Max( Max( A1, A2 ), Max( A3, A4 ) ); - - const double eps = 0.1; // val is in degrees + return NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4; +} - return val < eps ? 0. : val; +double Warping::GetValue( const TSequenceOfXYZ& P ) +{ + return ComputeValue(P); } double Warping::ComputeA( const gp_XYZ& thePnt1, @@ -1360,6 +1450,25 @@ double Warping::ComputeA( const gp_XYZ& thePnt1, return asin( fabs( H / L ) ) * 180. / M_PI; } +double Warping::ComputeValue(const TSequenceOfXYZ& thePoints) const +{ + if (thePoints.size() != 4) + return 0; + + gp_XYZ G = (thePoints(1) + thePoints(2) + thePoints(3) + thePoints(4)) / 4.; + + double A1 = ComputeA(thePoints(1), thePoints(2), thePoints(3), G); + double A2 = ComputeA(thePoints(2), thePoints(3), thePoints(4), G); + double A3 = ComputeA(thePoints(3), thePoints(4), thePoints(1), G); + double A4 = ComputeA(thePoints(4), thePoints(1), thePoints(2), G); + + double val = Max(Max(A1, A2), Max(A3, A4)); + + const double eps = 0.1; // val is in degrees + + return val < eps ? 0. : val; +} + double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const { // the warp is in the range [0.0,PI/2] @@ -1374,6 +1483,93 @@ SMDSAbs_ElementType Warping::GetType() const } +//================================================================================ +/* + Class : Warping3D + Description : Functor for calculating warping +*/ +//================================================================================ + +bool Warping3D::IsApplicable(const SMDS_MeshElement* element) const +{ + return NumericalFunctor::IsApplicable(element);//&& element->NbNodes() == 4; +} + +double Warping3D::GetValue(long theId) +{ + double aVal = 0; + myCurrElement = myMesh->FindElement(theId); + if (myCurrElement) + { + WValues aValues; + ProcessVolumeELement(aValues); + for (const auto& aValue: aValues) + { + aVal = Max(aVal, aValue.myWarp); + } + } + return aVal; +} + +double Warping3D::GetValue(const TSequenceOfXYZ& P) +{ + return ComputeValue(P); +} + +SMDSAbs_ElementType Warping3D::GetType() const +{ + return SMDSAbs_Volume; +} + +bool Warping3D::Value::operator<(const Warping3D::Value& x) const +{ + if (myPntIds.size() != x.myPntIds.size()) + return myPntIds.size() < x.myPntIds.size(); + + for (int anInd = 0; anInd < myPntIds.size(); ++anInd) + if (myPntIds[anInd] != x.myPntIds[anInd]) + return myPntIds[anInd] != x.myPntIds[anInd]; + + return false; +} + +// Compute value on each face of volume +void Warping3D::ProcessVolumeELement(WValues& theValues) +{ + SMDS_VolumeTool aVTool(myCurrElement); + double aCoord[3]; + for (int aFaceID = 0; aFaceID < aVTool.NbFaces(); ++aFaceID) + { + TSequenceOfXYZ aPoints; + std::set aNodes; + std::vector aNodeIds; + const SMDS_MeshNode** aNodesPtr = aVTool.GetFaceNodes(aFaceID); + + if (aNodesPtr) + { + for (int i = 0; i < aVTool.NbFaceNodes(aFaceID); ++i) + { + aNodesPtr[i]->GetXYZ(aCoord); + aPoints.push_back(gp_XYZ{ aCoord[0], aCoord[1], aCoord[2] }); + aNodeIds.push_back(aNodesPtr[i]->GetID()); + } + double aWarp = GetValue(aPoints); + Value aVal{ aWarp, aNodeIds }; + + theValues.push_back(aVal); + } + } +} + +void Warping3D::GetValues(WValues& theValues) +{ + for (SMDS_VolumeIteratorPtr anIter = myMesh->volumesIterator(); anIter->more(); ) + { + myCurrElement = anIter->next(); + ProcessVolumeELement(theValues); + } +} + //================================================================================ /* Class : Taper @@ -1381,6 +1577,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 ) @@ -1439,6 +1640,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 ) @@ -1553,13 +1759,36 @@ 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 */ //================================================================================ +Length2D::Length2D( SMDSAbs_ElementType type ): + myType ( type ) +{ +} + +bool Length2D::IsApplicable( const SMDS_MeshElement* element ) const +{ + return ( NumericalFunctor::IsApplicable( element ) && + element->GetEntityType() != SMDSEntity_Polyhedra ); +} + double Length2D::GetValue( const TSequenceOfXYZ& P ) { double aVal = 0; @@ -1804,7 +2033,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): @@ -1826,92 +2055,96 @@ 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(); + smIdType 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] = { 0,0,0,0 }; - gp_Pnt P[4]; - - double aLength = 0; - 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] = {0,0}; - 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); - } + else { + SMDS_NodeIteratorPtr aNodesIter = anElem->nodeIterator(); + smIdType aNodeId[2] = {0,0}; + gp_Pnt P[3]; + + 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(); + smIdType anId = aNode->GetID(); - aLength = P[0].Distance(P[1]); + P[2] = SMESH_NodeXYZ( aNode ); - Value aValue(aLength,aNodeId[0],aNodeId[1]); - theValues.insert(aValue); + 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 : Functor for calculating number of faces conneted to the edge + Description : computes distance between a face center and an underlying surface */ //================================================================================ @@ -1947,8 +2180,7 @@ double Deflection2D::GetValue( const TSequenceOfXYZ& P ) { gc += P(i+1); - if ( const SMDS_FacePosition* fPos = dynamic_cast - ( P.getElement()->GetNode( i )->GetPosition() )) + if ( SMDS_FacePositionPtr fPos = P.getElement()->GetNode( i )->GetPosition() ) { uv.ChangeCoord(1) += fPos->GetUParameter(); uv.ChangeCoord(2) += fPos->GetVParameter(); @@ -2006,7 +2238,7 @@ double Deflection2D::GetBadRate( double Value, int /*nbNodes*/ ) const */ //================================================================================ -double MultiConnection::GetValue( const TSequenceOfXYZ& P ) +double MultiConnection::GetValue( const TSequenceOfXYZ& /*P*/ ) { return 0; } @@ -2033,7 +2265,7 @@ SMDSAbs_ElementType MultiConnection::GetType() const */ //================================================================================ -double MultiConnection2D::GetValue( const TSequenceOfXYZ& P ) +double MultiConnection2D::GetValue( const TSequenceOfXYZ& /*P*/ ) { return 0; } @@ -2053,7 +2285,7 @@ double MultiConnection2D::GetValue( long theElementId ) if (!anIter) break; const SMDS_MeshNode *aNode, *aNode0 = 0; - TColStd_MapOfInteger aMap, aMapPrev; + NCollection_Map< smIdType, smIdHasher > aMap, aMapPrev; for (i = 0; i <= len; i++) { aMapPrev = aMap; @@ -2075,7 +2307,7 @@ double MultiConnection2D::GetValue( long theElementId ) while (anElemIter->more()) { const SMDS_MeshElement* anElem = anElemIter->next(); if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) { - int anId = anElem->GetID(); + smIdType anId = anElem->GetID(); aMap.Add(anId); if (aMapPrev.Contains(anId)) { @@ -2124,59 +2356,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] = {0,0,0}; + 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: "<::max() ) + aStr << aMaxId; // find position of the string in result sequence and insert string in it if ( anIntSeq.Length() == 0 ) { anIntSeq.Append( aMinId ); - aStrSeq.Append( aStr ); + aStrSeq.Append( (const char*)aStr ); } else { if ( aMinId < anIntSeq.First() ) { anIntSeq.Prepend( aMinId ); - aStrSeq.Prepend( aStr ); + aStrSeq.Prepend( (const char*)aStr ); } else if ( aMinId > anIntSeq.Last() ) { anIntSeq.Append( aMinId ); - aStrSeq.Append( aStr ); + aStrSeq.Append( (const char*)aStr ); } else for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ ) if ( aMinId < anIntSeq( j ) ) { anIntSeq.InsertBefore( j, aMinId ); - aStrSeq.InsertBefore( j, aStr ); + aStrSeq.InsertBefore( j, (const char*)aStr ); break; } } @@ -3324,13 +3588,14 @@ void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr ) if ( aStrSeq.Length() == 0 ) return; - - theResStr = aStrSeq( 1 ); + std::string aResStr; + aResStr = aStrSeq( 1 ); for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ ) { - theResStr += ","; - theResStr += aStrSeq( j ); + aResStr += ","; + aResStr += aStrSeq( j ); } + theResStr = aResStr.c_str(); } //======================================================================= @@ -3340,8 +3605,8 @@ void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr ) //======================================================================= bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr ) { - myMin.Clear(); - myMax.Clear(); + myMin.clear(); + myMax.clear(); myIds.Clear(); TCollection_AsciiString aStr = theStr; @@ -3379,8 +3644,8 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr ) (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) ) return false; - myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() ); - myMax.Append( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() ); + myMin.push_back( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() ); + myMax.push_back( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() ); } } @@ -3429,8 +3694,8 @@ bool RangeOfIds::IsSatisfy( long theId ) if ( myIds.Contains( theId ) ) return true; - for ( int i = 1, n = myMin.Length(); i <= n; i++ ) - if ( theId >= myMin( i ) && theId <= myMax( i ) ) + for ( size_t i = 0; i < myMin.size(); i++ ) + if ( theId >= myMin[i] && theId <= myMax[i] ) return true; return false; @@ -3658,9 +3923,10 @@ void Filter::SetPredicate( PredicatePtr thePredicate ) myPredicate = thePredicate; } -void Filter::GetElementsId( const SMDS_Mesh* theMesh, - PredicatePtr thePredicate, - TIdSequence& theSequence ) +void Filter::GetElementsId( const SMDS_Mesh* theMesh, + PredicatePtr thePredicate, + TIdSequence& theSequence, + SMDS_ElemIteratorPtr theElements ) { theSequence.clear(); @@ -3669,21 +3935,28 @@ void Filter::GetElementsId( const SMDS_Mesh* theMesh, thePredicate->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); } /* @@ -3798,7 +4071,7 @@ bool ManifoldPart::process() // the map of non manifold links and bad geometry TMapOfLink aMapOfNonManifold; - TColStd_MapOfInteger aMapOfTreated; + TIDsMap aMapOfTreated; // begin cycle on faces from start index and run on vector till the end // and from begin to start index to cover whole vector @@ -3811,18 +4084,18 @@ bool ManifoldPart::process() // as result next time when fi will be equal to aStartIndx SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ]; - if ( aMapOfTreated.Contains( aFacePtr->GetID() ) ) + if ( aMapOfTreated.Contains( aFacePtr->GetID()) ) continue; aMapOfTreated.Add( aFacePtr->GetID() ); - TColStd_MapOfInteger aResFaces; + TIDsMap aResFaces; if ( !findConnected( myAllFacePtrIntDMap, aFacePtr, aMapOfNonManifold, aResFaces ) ) continue; - TColStd_MapIteratorOfMapOfInteger anItr( aResFaces ); + TIDsMap::Iterator anItr( aResFaces ); for ( ; anItr.More(); anItr.Next() ) { - int aFaceId = anItr.Key(); + smIdType aFaceId = anItr.Key(); aMapOfTreated.Add( aFaceId ); myMapIds.Add( aFaceId ); } @@ -3858,7 +4131,7 @@ bool ManifoldPart::findConnected ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt, SMDS_MeshFace* theStartFace, ManifoldPart::TMapOfLink& theNonManifold, - TColStd_MapOfInteger& theResFaces ) + TIDsMap& theResFaces ) { theResFaces.Clear(); if ( !theAllFacePtrInt.size() ) @@ -3926,7 +4199,7 @@ bool ManifoldPart::findConnected SMDS_MeshFace* aNextFace = *pFace; if ( aPrevFace == aNextFace ) continue; - int anNextFaceID = aNextFace->GetID(); + smIdType anNextFaceID = aNextFace->GetID(); if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) ) // should not be with non manifold restriction. probably bad topology continue; @@ -4008,24 +4281,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 ); } } @@ -4115,8 +4384,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 @@ -4159,7 +4430,9 @@ void ElementsOnSurface::process() if ( !myMeshModifTracer.GetMesh() ) return; - myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType )); + int nbElems = FromSmIdType( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType )); + if ( nbElems > 0 ) + myIds.ReSize( nbElems ); SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType ); for(; anIter->more(); ) @@ -4226,18 +4499,20 @@ namespace { struct ElementsOnShape::Classifier { - Classifier() { mySolidClfr = 0; myFlags = 0; } + Classifier(): mySolidClfr(0), myProjFace(0), myProjEdge(0), myFlags(0) { myU = myV = 1e100; } ~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; } + void GetParams( double & u, double & v ) const { u = myU; v = myV; } private: bool isOutOfSolid (const gp_Pnt& p); @@ -4245,16 +4520,20 @@ private: 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 + BRepClass3d_SolidClassifier* mySolidClfr; Bnd_B3d myBox; - GeomAPI_ProjectPointOnSurf myProjFace; - GeomAPI_ProjectPointOnCurve myProjEdge; + GeomAPI_ProjectPointOnSurf* myProjFace; + GeomAPI_ProjectPointOnCurve* myProjEdge; gp_Pnt myVertexXYZ; TopoDS_Shape myShape; double myTol; + double myU, myV; // result of isOutOfFace() and isOutOfEdge() int myFlags; }; @@ -4266,6 +4545,8 @@ struct ElementsOnShape::OctreeClassifier : public SMESH_Octree 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; } @@ -4291,6 +4572,22 @@ ElementsOnShape::~ElementsOnShape() 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 + { + + if (SALOME::VerbosityActivated()) + std::cout << "Avoid ElementsOnShape::clone(), too large: " << size << " bytes " << std::endl; + + return 0; + } + ElementsOnShape* cln = new ElementsOnShape(); cln->SetAllNodes ( myAllNodesFlag ); cln->SetTolerance( myToler ); @@ -4315,9 +4612,12 @@ SMDSAbs_ElementType ElementsOnShape::GetType() const void ElementsOnShape::SetTolerance (const double theToler) { - if (myToler != theToler) { + if (myToler != theToler) + { myToler = theToler; - SetShape(myShape, myType); + TopoDS_Shape s = myShape; + myShape.Nullify(); + SetShape( s, myType ); } } @@ -4447,12 +4747,13 @@ bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem) for ( size_t i = 0; i < myClassifiers.size(); ++i ) myWorkClassifiers[ i ] = & myClassifiers[ i ]; myOctree = new OctreeClassifier( myWorkClassifiers ); + + SMESHUtils::FreeVector( myWorkClassifiers ); } - SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator(); - while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) + for ( int i = 0, nb = elem->NbNodes(); i < nb && (isSatisfy == myAllNodesFlag); ++i ) { - SMESH_TNodeXYZ aPnt( aNodeItr->next() ); + SMESH_TNodeXYZ aPnt( elem->GetNode( i )); centerXYZ += aPnt; isNodeOut = true; @@ -4488,16 +4789,28 @@ bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem) centerXYZ /= elem->NbNodes(); isSatisfy = false; 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; } +//================================================================================ +/*! + * \brief Check and optionally return a satisfying shape + */ +//================================================================================ + bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node, TopoDS_Shape* okShape) { @@ -4532,6 +4845,7 @@ bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node, isNodeOut = false; if ( okShape ) *okShape = myWorkClassifiers[i]->Shape(); + myWorkClassifiers[i]->GetParams( myU, myV ); break; } } @@ -4542,7 +4856,8 @@ bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node, { isNodeOut = false; if ( okShape ) - *okShape = myWorkClassifiers[i]->Shape(); + *okShape = myClassifiers[i].Shape(); + myClassifiers[i].GetParams( myU, myV ); break; } } @@ -4571,7 +4886,7 @@ void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape, } else { - mySolidClfr = new BRepClass3d_SolidClassifier(theShape); + mySolidClfr = new BRepClass3d_SolidClassifier( prepareSolid( theShape )); myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid; } break; @@ -4580,17 +4895,29 @@ void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape, { 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::Classifier::isOutOfFace; + if ( surf.IsNull() ) + myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone; + else + { + surf->Bounds( u1,u2,v1,v2 ); + myProjFace = new GeomAPI_ProjectPointOnSurf; + myProjFace->Init( surf, u1,u2, v1,v2, myTol ); + myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace; + } break; } case TopAbs_EDGE: { Standard_Real u1, u2; Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2); - myProjEdge.Init(curve, u1, u2); - myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge; + if ( curve.IsNull() ) + myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone; + else + { + myProjEdge = new GeomAPI_ProjectPointOnCurve; + myProjEdge->Init( curve, u1, u2 ); + myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge; + } break; } case TopAbs_VERTEX: @@ -4612,7 +4939,15 @@ void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape, else { Bnd_Box box; - BRepBndLib::Add( myShape, 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() ); @@ -4630,28 +4965,50 @@ void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape, ElementsOnShape::Classifier::~Classifier() { delete mySolidClfr; mySolidClfr = 0; + delete myProjFace; myProjFace = 0; + delete myProjEdge; myProjEdge = 0; +} + +TopoDS_Shape ElementsOnShape::Classifier::prepareSolid( const TopoDS_Shape& theSolid ) +{ + // 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::Classifier::isOutOfSolid (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) +bool ElementsOnShape::Classifier::isOutOfBox( const gp_Pnt& p ) { return myBox.IsOut( p.XYZ() ); } -bool ElementsOnShape::Classifier::isOutOfFace (const gp_Pnt& p) +bool ElementsOnShape::Classifier::isOutOfFace( const gp_Pnt& p ) { - myProjFace.Perform( p ); - if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol ) + if ( isOutOfBox( p )) return true; + myProjFace->Perform( p ); + if ( myProjFace->IsDone() && myProjFace->LowerDistance() <= myTol ) { // check relatively to the face - Standard_Real u, v; - myProjFace.LowerDistanceParameters(u, v); - gp_Pnt2d aProjPnt (u, v); + myProjFace->LowerDistanceParameters( myU, myV ); + gp_Pnt2d aProjPnt( myU, myV ); BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol ); if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON ) return false; @@ -4659,18 +5016,22 @@ bool ElementsOnShape::Classifier::isOutOfFace (const gp_Pnt& p) return true; } -bool ElementsOnShape::Classifier::isOutOfEdge (const gp_Pnt& p) +bool ElementsOnShape::Classifier::isOutOfEdge( const gp_Pnt& p ) { - myProjEdge.Perform( p ); - return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol ); + if ( isOutOfBox( p )) return true; + myProjEdge->Perform( p ); + bool isOn = ( myProjEdge->NbPoints() > 0 && myProjEdge->LowerDistance() <= myTol ); + if ( isOn ) + myU = myProjEdge->LowerDistanceParameter(); + return !isOn; } -bool ElementsOnShape::Classifier::isOutOfVertex(const gp_Pnt& p) +bool ElementsOnShape::Classifier::isOutOfVertex( const gp_Pnt& p ) { return ( myVertexXYZ.Distance( p ) > myTol ); } -bool ElementsOnShape::Classifier::isBox (const TopoDS_Shape& theShape) +bool ElementsOnShape::Classifier::isBox(const TopoDS_Shape& theShape ) { TopTools_IndexedMapOfShape vMap; TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap ); @@ -4752,6 +5113,19 @@ OctreeClassifier::GetClassifiersAtPoint( const gp_XYZ& point, } } +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 @@ -4798,7 +5172,8 @@ void ElementsOnShape::OctreeClassifier::buildChildrenData() for ( int i = 0; i < nbChildren(); i++ ) { OctreeClassifier* child = static_cast( myChildren[ i ]); - child->myIsLeaf = ( child->myClassifiers.size() <= 5 ); + child->myIsLeaf = ( child->myClassifiers.size() <= 5 || + child->maxSize() < child->myClassifiers[0]->Tolerance() ); } } @@ -4825,8 +5200,13 @@ BelongToGeom::BelongToGeom() Predicate* BelongToGeom::clone() const { - BelongToGeom* cln = new BelongToGeom( *this ); - cln->myElementsOnShapePtr.reset( static_cast( myElementsOnShapePtr->clone() )); + BelongToGeom* cln = 0; + if ( myElementsOnShapePtr ) + if ( ElementsOnShape* eos = static_cast( myElementsOnShapePtr->clone() )) + { + cln = new BelongToGeom( *this ); + cln->myElementsOnShapePtr.reset( eos ); + } return cln; } @@ -4837,6 +5217,8 @@ void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh ) myMeshDS = dynamic_cast(theMesh); init(); } + if ( myElementsOnShapePtr ) + myElementsOnShapePtr->SetMesh( myMeshDS ); } void BelongToGeom::SetGeom( const TopoDS_Shape& theShape ) @@ -4933,7 +5315,7 @@ bool BelongToGeom::IsSatisfy (long theId) { if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId )) { - if ( anElem->GetType() == myType ) + if ( myType == SMDSAbs_All || anElem->GetType() == myType ) { if ( anElem->getshapeId() < 1 ) return myElementsOnShapePtr->IsSatisfy(theId); @@ -4996,8 +5378,13 @@ LyingOnGeom::LyingOnGeom() Predicate* LyingOnGeom::clone() const { - LyingOnGeom* cln = new LyingOnGeom( *this ); - cln->myElementsOnShapePtr.reset( static_cast( myElementsOnShapePtr->clone() )); + LyingOnGeom* cln = 0; + if ( myElementsOnShapePtr ) + if ( ElementsOnShape* eos = static_cast( myElementsOnShapePtr->clone() )) + { + cln = new LyingOnGeom( *this ); + cln->myElementsOnShapePtr.reset( eos ); + } return cln; } @@ -5008,6 +5395,8 @@ void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh ) myMeshDS = dynamic_cast(theMesh); init(); } + if ( myElementsOnShapePtr ) + myElementsOnShapePtr->SetMesh( myMeshDS ); } void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape ) @@ -5073,7 +5462,8 @@ bool LyingOnGeom::IsSatisfy( long theId ) if ( mySubShapesIDs.Contains( elem->getshapeId() )) return true; - if ( elem->GetType() != SMDSAbs_Node && elem->GetType() == myType ) + if (( elem->GetType() != SMDSAbs_Node ) && + ( myType == SMDSAbs_All || elem->GetType() == myType )) { SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator(); while ( nodeItr->more() )