X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=f86b5638719dfb201d6e15ab0a6e4487a1b0f174;hp=8ba13a8f1f4e0727d548282be267cb69b4fbe571;hb=b54f78cc796b6ae41d77cc43b55c01b16a4c0445;hpb=0635c9fc80f67d1e5dc0e94ec85f487286a92070 diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 8ba13a8f1..f86b56387 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -1,50 +1,69 @@ -// Copyright (C) 2007-2008 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 +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS // -// This library is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 2.1 of the License. +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. // -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // + #include "SMESH_ControlsDef.hxx" -#include +#include "SMDS_BallElement.hxx" +#include "SMDS_FacePosition.hxx" +#include "SMDS_Iterator.hxx" +#include "SMDS_Mesh.hxx" +#include "SMDS_MeshElement.hxx" +#include "SMDS_MeshNode.hxx" +#include "SMDS_QuadraticEdge.hxx" +#include "SMDS_QuadraticFaceOfNodes.hxx" +#include "SMDS_VolumeTool.hxx" +#include "SMESHDS_GroupBase.hxx" +#include "SMESHDS_GroupOnFilter.hxx" +#include "SMESHDS_Mesh.hxx" +#include "SMESH_MeshAlgos.hxx" +#include "SMESH_OctreeNode.hxx" + +#include #include +#include +#include +#include #include #include - -#include -#include -#include -#include -#include -#include -#include - +#include #include #include #include - +#include #include +#include #include #include #include #include - +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include @@ -53,22 +72,25 @@ #include #include -#include "SMDS_Mesh.hxx" -#include "SMDS_Iterator.hxx" -#include "SMDS_MeshElement.hxx" -#include "SMDS_MeshNode.hxx" -#include "SMDS_VolumeTool.hxx" -#include "SMDS_QuadraticFaceOfNodes.hxx" -#include "SMDS_QuadraticEdge.hxx" +#include -#include "SMESHDS_Mesh.hxx" -#include "SMESHDS_GroupBase.hxx" +#include +#include /* AUXILIARY METHODS */ -namespace{ +namespace { + + const double theEps = 1e-100; + const double theInf = 1e+100; + + inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode ) + { + return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() ); + } + inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) { gp_Vec v1( P1 - P2 ), v2( P3 - P2 ); @@ -77,6 +99,15 @@ namespace{ v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 ); } + inline double getCos2( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) + { + gp_Vec v1( P1 - P2 ), v2( P3 - P2 ); + double dot = v1 * v2, len1 = v1.SquareMagnitude(), len2 = v2.SquareMagnitude(); + + return ( dot < 0 || len1 < gp::Resolution() || len2 < gp::Resolution() ? -1 : + dot * dot / len1 / len2 ); + } + inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) { gp_Vec aVec1( P2 - P1 ); @@ -116,13 +147,13 @@ namespace{ // +-----+------+ +-----+------+ // | | | | // | | | | - // result sould be 2 in both cases + // result should be 2 in both cases // int aResult0 = 0, aResult1 = 0; // last node, it is a medium one in a quadratic edge const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 ); - const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 ); - const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 ); + const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 ); + const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 ); if ( aNode1 == aLastNode ) aNode1 = 0; SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator(); @@ -144,32 +175,29 @@ namespace{ } int aResult = std::max ( aResult0, aResult1 ); -// TColStd_MapOfInteger aMap; - -// SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator(); -// if ( anIter != 0 ) { -// while( anIter->more() ) { -// const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); -// if ( aNode == 0 ) -// return 0; -// SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); -// while( anElemIter->more() ) { -// const SMDS_MeshElement* anElem = anElemIter->next(); -// if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) { -// int anId = anElem->GetID(); - -// if ( anIter->more() ) // i.e. first node -// aMap.Add( anId ); -// else if ( aMap.Contains( anId ) ) -// aResult++; -// } -// } -// } -// } - return aResult; } + gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 ) + { + int aNbNode = theFace->NbNodes(); + + gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0)); + gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0)); + gp_XYZ n = q1 ^ q2; + if ( aNbNode > 3 ) { + gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0)); + n += q2 ^ q3; + } + double len = n.Modulus(); + bool zeroLen = ( len <= std::numeric_limits::min()); + if ( !zeroLen ) + n /= len; + + if (ok) *ok = !zeroLen; + + return n; + } } @@ -177,13 +205,16 @@ namespace{ using namespace SMESH::Controls; /* - FUNCTORS -*/ + * FUNCTORS + */ +//================================================================================ /* Class : NumericalFunctor Description : Base class for numerical functors */ +//================================================================================ + NumericalFunctor::NumericalFunctor(): myMesh(NULL) { @@ -195,7 +226,7 @@ void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh ) myMesh = theMesh; } -bool NumericalFunctor::GetPoints(const int theId, +bool NumericalFunctor::GetPoints(const int theId, TSequenceOfXYZ& theRes ) const { theRes.clear(); @@ -203,7 +234,11 @@ bool NumericalFunctor::GetPoints(const int theId, if ( myMesh == 0 ) return false; - return GetPoints( myMesh->FindElement( theId ), theRes ); + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + if ( !anElem || anElem->GetType() != this->GetType() ) + return false; + + return GetPoints( anElem, theRes ); } bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, @@ -211,10 +246,11 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, { theRes.clear(); - if ( anElem == 0) + if ( anElem == 0 ) return false; theRes.reserve( anElem->NbNodes() ); + theRes.setElement( anElem ); // Get nodes of the element SMDS_ElemIteratorPtr anIter; @@ -222,16 +258,15 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, if ( anElem->IsQuadratic() ) { switch ( anElem->GetType() ) { case SMDSAbs_Edge: - anIter = static_cast + anIter = dynamic_cast (anElem)->interlacedNodesElemIterator(); break; case SMDSAbs_Face: - anIter = static_cast + anIter = dynamic_cast (anElem)->interlacedNodesElemIterator(); break; default: anIter = anElem->nodesIterator(); - //return false; } } else { @@ -239,9 +274,10 @@ bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, } if ( anIter ) { + SMESH_NodeXYZ p; while( anIter->more() ) { - if ( const SMDS_MeshNode* aNode = static_cast( anIter->next() )) - theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); + if ( p.Set( anIter->next() )) + theRes.push_back( p ); } } @@ -256,30 +292,125 @@ long NumericalFunctor::GetPrecision() const void NumericalFunctor::SetPrecision( const long thePrecision ) { myPrecision = thePrecision; + myPrecisionValue = pow( 10., (double)( myPrecision ) ); } double NumericalFunctor::GetValue( long theId ) { + double aVal = 0; + myCurrElement = myMesh->FindElement( theId ); + TSequenceOfXYZ P; - if ( GetPoints( theId, P )) + if ( GetPoints( theId, P )) // elem type is checked here + aVal = Round( GetValue( P )); + + return aVal; +} + +double NumericalFunctor::Round( const double & aVal ) +{ + return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal; +} + +//================================================================================ +/*! + * \brief Return histogram of functor values + * \param nbIntervals - number of intervals + * \param nbEvents - number of mesh elements having values within i-th interval + * \param funValues - boundaries of intervals + * \param elements - elements to check vulue of; empty list means "of all" + * \param minmax - boundaries of diapason of values to divide into intervals + */ +//================================================================================ + +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 || + !myMesh->GetMeshInfo().NbElements( GetType() )) + return; + nbEvents.resize( nbIntervals, 0 ); + funValues.resize( nbIntervals+1 ); + + // get all values sorted + std::multiset< double > values; + if ( elements.empty() ) + { + SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator( GetType() ); + while ( elemIt->more() ) + values.insert( GetValue( elemIt->next()->GetID() )); + } + else + { + std::vector::const_iterator id = elements.begin(); + for ( ; id != elements.end(); ++id ) + values.insert( GetValue( *id )); + } + + if ( minmax ) + { + funValues[0] = minmax[0]; + funValues[nbIntervals] = minmax[1]; + } + else + { + funValues[0] = *values.begin(); + funValues[nbIntervals] = *values.rbegin(); + } + // case nbIntervals == 1 + if ( nbIntervals == 1 ) + { + nbEvents[0] = values.size(); + return; + } + // case of 1 value + if (funValues.front() == funValues.back()) + { + nbEvents.resize( 1 ); + nbEvents[0] = values.size(); + funValues[1] = funValues.back(); + funValues.resize( 2 ); + } + // generic case + std::multiset< double >::iterator min = values.begin(), max; + for ( int i = 0; i < nbIntervals; ++i ) { - double aVal = GetValue( P ); - if ( myPrecision >= 0 ) + // find end value of i-th interval + double r = (i+1) / double(nbIntervals); + if (isLogarithmic && funValues.front() > 1e-07 && funValues.back() > 1e-07) { + double logmin = log10(funValues.front()); + double lval = logmin + r * (log10(funValues.back()) - logmin); + funValues[i+1] = pow(10.0, lval); + } + else { + funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r; + } + + // count values in the i-th interval if there are any + if ( min != values.end() && *min <= funValues[i+1] ) { - double prec = pow( 10., (double)( myPrecision ) ); - aVal = floor( aVal * prec + 0.5 ) / prec; + // find the first value out of the interval + max = values.upper_bound( funValues[i+1] ); // max is greater than funValues[i+1], or end() + nbEvents[i] = std::distance( min, max ); + min = max; } - return aVal; } - - return 0.; + // add values larger than minmax[1] + nbEvents.back() += std::distance( min, values.end() ); } //======================================================================= -//function : GetValue -//purpose : -//======================================================================= +/* + Class : Volume + Description : Functor calculating volume of a 3D element +*/ +//================================================================================ double Volume::GetValue( long theElementId ) { @@ -291,48 +422,326 @@ double Volume::GetValue( long theElementId ) return 0; } -//======================================================================= -//function : GetBadRate -//purpose : meaningless as it is not quality control functor -//======================================================================= - double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const { return Value; } +SMDSAbs_ElementType Volume::GetType() const +{ + return SMDSAbs_Volume; +} + //======================================================================= -//function : GetType -//purpose : +/* + Class : MaxElementLength2D + Description : Functor calculating maximum length of 2D element +*/ +//================================================================================ + +double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P ) +{ + if(P.size() == 0) + return 0.; + double aVal = 0; + int len = P.size(); + if( len == 3 ) { // triangles + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + aVal = Max(L1,Max(L2,L3)); + } + else if( len == 4 ) { // quadrangles + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + double D1 = getDistance(P( 1 ),P( 3 )); + double D2 = getDistance(P( 2 ),P( 4 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2)); + } + else if( len == 6 ) { // quadratic triangles + double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 )); + double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 )); + double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 )); + aVal = Max(L1,Max(L2,L3)); + } + else if( len == 8 || len == 9 ) { // quadratic quadrangles + double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 )); + double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 )); + double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 )); + double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 )); + double D1 = getDistance(P( 1 ),P( 5 )); + double D2 = getDistance(P( 3 ),P( 7 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2)); + } + // Diagonals are undefined for concave polygons + // else if ( P.getElementEntity() == SMDSEntity_Quad_Polygon && P.size() > 2 ) // quad polygon + // { + // // sides + // aVal = getDistance( P( 1 ), P( P.size() )) + getDistance( P( P.size() ), P( P.size()-1 )); + // for ( size_t i = 1; i < P.size()-1; i += 2 ) + // { + // double L = getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 )); + // aVal = Max( aVal, L ); + // } + // // diagonals + // for ( int i = P.size()-5; i > 0; i -= 2 ) + // for ( int j = i + 4; j < P.size() + i - 2; i += 2 ) + // { + // double D = getDistance( P( i ), P( j )); + // aVal = Max( aVal, D ); + // } + // } + // { // polygons + + // } + + if( myPrecision >= 0 ) + { + double prec = pow( 10., (double)myPrecision ); + aVal = floor( aVal * prec + 0.5 ) / prec; + } + return aVal; +} + +double MaxElementLength2D::GetValue( long theElementId ) +{ + TSequenceOfXYZ P; + return GetPoints( theElementId, P ) ? GetValue(P) : 0.0; +} + +double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; +} + +SMDSAbs_ElementType MaxElementLength2D::GetType() const +{ + return SMDSAbs_Face; +} + //======================================================================= +/* + Class : MaxElementLength3D + Description : Functor calculating maximum length of 3D element +*/ +//================================================================================ -SMDSAbs_ElementType Volume::GetType() const +double MaxElementLength3D::GetValue( long theElementId ) { - return SMDSAbs_Volume; + TSequenceOfXYZ P; + if( GetPoints( theElementId, P ) ) { + double aVal = 0; + const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId ); + SMDSAbs_EntityType aType = aElem->GetEntityType(); + int len = P.size(); + switch ( aType ) { + case SMDSEntity_Tetra: { // tetras + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + double L4 = getDistance(P( 1 ),P( 4 )); + double L5 = getDistance(P( 2 ),P( 4 )); + double L6 = getDistance(P( 3 ),P( 4 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + break; + } + case SMDSEntity_Pyramid: { // pyramids + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + double L5 = getDistance(P( 1 ),P( 5 )); + double L6 = getDistance(P( 2 ),P( 5 )); + double L7 = getDistance(P( 3 ),P( 5 )); + double L8 = getDistance(P( 4 ),P( 5 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(L7,L8)); + break; + } + case SMDSEntity_Penta: { // pentas + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 4 )); + double L7 = getDistance(P( 1 ),P( 4 )); + double L8 = getDistance(P( 2 ),P( 5 )); + double L9 = getDistance(P( 3 ),P( 6 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),L9)); + break; + } + case SMDSEntity_Hexa: { // hexas + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 7 )); + double L7 = getDistance(P( 7 ),P( 8 )); + double L8 = getDistance(P( 8 ),P( 5 )); + double L9 = getDistance(P( 1 ),P( 5 )); + double L10= getDistance(P( 2 ),P( 6 )); + double L11= getDistance(P( 3 ),P( 7 )); + double L12= getDistance(P( 4 ),P( 8 )); + double D1 = getDistance(P( 1 ),P( 7 )); + double D2 = getDistance(P( 2 ),P( 8 )); + double D3 = getDistance(P( 3 ),P( 5 )); + double D4 = getDistance(P( 4 ),P( 6 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10))); + aVal = Max(aVal,Max(L11,L12)); + aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); + break; + } + case SMDSEntity_Hexagonal_Prism: { // hexagonal prism + for ( int i1 = 1; i1 < 12; ++i1 ) + for ( int i2 = i1+1; i1 <= 12; ++i1 ) + aVal = Max( aVal, getDistance(P( i1 ),P( i2 ))); + break; + } + case SMDSEntity_Quad_Tetra: { // quadratic tetras + double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 )); + double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 )); + double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + break; + } + case SMDSEntity_Quad_Pyramid: { // quadratic pyramids + double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 )); + double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 )); + double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(L7,L8)); + break; + } + case SMDSEntity_Quad_Penta: + case SMDSEntity_BiQuad_Penta: { // quadratic pentas + double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 )); + double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 )); + double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 )); + double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),L9)); + break; + } + case SMDSEntity_Quad_Hexa: + case SMDSEntity_TriQuad_Hexa: { // quadratic hexas + double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 )); + double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 )); + double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 )); + double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 )); + double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 )); + double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 )); + double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 )); + double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 )); + double D1 = getDistance(P( 1 ),P( 7 )); + double D2 = getDistance(P( 2 ),P( 8 )); + double D3 = getDistance(P( 3 ),P( 5 )); + double D4 = getDistance(P( 4 ),P( 6 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10))); + aVal = Max(aVal,Max(L11,L12)); + aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); + break; + } + case SMDSEntity_Quad_Polyhedra: + case SMDSEntity_Polyhedra: { // polys + // get the maximum distance between all pairs of nodes + for( int i = 1; i <= len; i++ ) { + for( int j = 1; j <= len; j++ ) { + if( j > i ) { // optimization of the loop + double D = getDistance( P(i), P(j) ); + aVal = Max( aVal, D ); + } + } + } + break; + } + case SMDSEntity_Node: + case SMDSEntity_0D: + case SMDSEntity_Edge: + case SMDSEntity_Quad_Edge: + case SMDSEntity_Triangle: + case SMDSEntity_Quad_Triangle: + case SMDSEntity_BiQuad_Triangle: + case SMDSEntity_Quadrangle: + case SMDSEntity_Quad_Quadrangle: + case SMDSEntity_BiQuad_Quadrangle: + case SMDSEntity_Polygon: + case SMDSEntity_Quad_Polygon: + case SMDSEntity_Ball: + case SMDSEntity_Last: return 0; + } // switch ( aType ) + + if( myPrecision >= 0 ) + { + double prec = pow( 10., (double)myPrecision ); + aVal = floor( aVal * prec + 0.5 ) / prec; + } + return aVal; + } + return 0.; +} + +double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; } +SMDSAbs_ElementType MaxElementLength3D::GetType() const +{ + return SMDSAbs_Volume; +} +//======================================================================= /* Class : MinimumAngle Description : Functor for calculation of minimum angle */ +//================================================================================ double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) { - double aMin; - - if (P.size() <3) + if ( P.size() < 3 ) return 0.; - aMin = getAngle(P( P.size() ), P( 1 ), P( 2 )); - aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 ))); + double aMaxCos2; - for (int i=2; i= 1 ) return 0; + return acos( cos ) * 180.0 / M_PI; } double MinimumAngle::GetBadRate( double Value, int nbNodes ) const @@ -348,10 +757,33 @@ SMDSAbs_ElementType MinimumAngle::GetType() const } +//================================================================================ /* Class : AspectRatio Description : Functor for calculating aspect ratio */ +//================================================================================ + +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 )); + } + return aVal; +} + double AspectRatio::GetValue( const TSequenceOfXYZ& P ) { // According to "Mesh quality control" by Nadir Bouhamau referring to @@ -368,94 +800,135 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P ) if ( nbNodes == 3 ) { // Compute lengths of the sides - std::vector< double > aLen (nbNodes); - for ( int i = 0; i < nbNodes - 1; i++ ) - aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) ); - aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) ); + double aLen1 = getDistance( P( 1 ), P( 2 )); + double aLen2 = getDistance( P( 2 ), P( 3 )); + double aLen3 = getDistance( P( 3 ), P( 1 )); // Q = alfa * h * p / S, where // // alfa = sqrt( 3 ) / 6 // h - length of the longest edge // p - half perimeter // S - triangle surface - const double alfa = sqrt( 3. ) / 6.; - double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) ); - double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.; - double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) ); - if ( anArea <= Precision::Confusion() ) - return 0.; + const double alfa = sqrt( 3. ) / 6.; + double maxLen = Max( aLen1, Max( aLen2, aLen3 )); + double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.; + double anArea = getArea( P( 1 ), P( 2 ), P( 3 )); + if ( anArea <= theEps ) + return theInf; return alfa * maxLen * half_perimeter / anArea; } else if ( nbNodes == 6 ) { // quadratic triangles // Compute lengths of the sides - std::vector< double > aLen (3); - aLen[0] = getDistance( P(1), P(3) ); - aLen[1] = getDistance( P(3), P(5) ); - aLen[2] = getDistance( P(5), P(1) ); - // Q = alfa * h * p / S, where - // - // alfa = sqrt( 3 ) / 6 - // h - length of the longest edge - // p - half perimeter - // S - triangle surface - const double alfa = sqrt( 3. ) / 6.; - double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) ); - double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.; - double anArea = getArea( P(1), P(3), P(5) ); - if ( anArea <= Precision::Confusion() ) - return 0.; + double aLen1 = getDistance( P( 1 ), P( 3 )); + double aLen2 = getDistance( P( 3 ), P( 5 )); + double aLen3 = getDistance( P( 5 ), P( 1 )); + // algo same as for the linear triangle + const double alfa = sqrt( 3. ) / 6.; + double maxLen = Max( aLen1, Max( aLen2, aLen3 )); + double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.; + double anArea = getArea( P( 1 ), P( 3 ), P( 5 )); + if ( anArea <= theEps ) + return theInf; return alfa * maxLen * half_perimeter / anArea; } else if( nbNodes == 4 ) { // quadrangle - // return aspect ratio of the worst triange which can be built + // Compute lengths of the sides + double aLen[4]; + aLen[0] = getDistance( P(1), P(2) ); + aLen[1] = getDistance( P(2), P(3) ); + aLen[2] = getDistance( P(3), P(4) ); + aLen[3] = getDistance( P(4), P(1) ); + // Compute lengths of the diagonals + double aDia[2]; + aDia[0] = getDistance( P(1), P(3) ); + aDia[1] = getDistance( P(2), P(4) ); + // Compute areas of all triangles which can be built // taking three nodes of the quadrangle - TSequenceOfXYZ triaPnts(3); - // triangle on nodes 1 3 2 - triaPnts(1) = P(1); - triaPnts(2) = P(3); - triaPnts(3) = P(2); - double ar = GetValue( triaPnts ); - // triangle on nodes 1 3 4 - triaPnts(3) = P(4); - ar = Max ( ar, GetValue( triaPnts )); - // triangle on nodes 1 2 4 - triaPnts(2) = P(2); - ar = Max ( ar, GetValue( triaPnts )); - // triangle on nodes 3 2 4 - triaPnts(1) = P(3); - ar = Max ( ar, GetValue( triaPnts )); - - return ar; - } - else { // nbNodes==8 - quadratic quadrangle - // return aspect ratio of the worst triange which can be built + double anArea[4]; + anArea[0] = getArea( P(1), P(2), P(3) ); + anArea[1] = getArea( P(1), P(2), P(4) ); + anArea[2] = getArea( P(1), P(3), P(4) ); + anArea[3] = getArea( P(2), P(3), P(4) ); + // Q = alpha * L * C1 / C2, where + // + // alpha = sqrt( 1/32 ) + // L = max( L1, L2, L3, L4, D1, D2 ) + // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 ) + // 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. ); + double C2 = Min( anArea[ 0 ], + Min( anArea[ 1 ], + Min( anArea[ 2 ], anArea[ 3 ] ) ) ); + if ( C2 <= theEps ) + return theInf; + return alpha * L * C1 / C2; + } + else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle + // Compute lengths of the sides + double aLen[4]; + aLen[0] = getDistance( P(1), P(3) ); + aLen[1] = getDistance( P(3), P(5) ); + aLen[2] = getDistance( P(5), P(7) ); + aLen[3] = getDistance( P(7), P(1) ); + // Compute lengths of the diagonals + double aDia[2]; + aDia[0] = getDistance( P(1), P(5) ); + aDia[1] = getDistance( P(3), P(7) ); + // Compute areas of all triangles which can be built // taking three nodes of the quadrangle - TSequenceOfXYZ triaPnts(3); - // triangle on nodes 1 3 2 - triaPnts(1) = P(1); - triaPnts(2) = P(5); - triaPnts(3) = P(3); - double ar = GetValue( triaPnts ); - // triangle on nodes 1 3 4 - triaPnts(3) = P(7); - ar = Max ( ar, GetValue( triaPnts )); - // triangle on nodes 1 2 4 - triaPnts(2) = P(3); - ar = Max ( ar, GetValue( triaPnts )); - // triangle on nodes 3 2 4 - triaPnts(1) = P(5); - ar = Max ( ar, GetValue( triaPnts )); - - return ar; + double anArea[4]; + anArea[0] = getArea( P(1), P(3), P(5) ); + anArea[1] = getArea( P(1), P(3), P(7) ); + anArea[2] = getArea( P(1), P(5), P(7) ); + anArea[3] = getArea( P(3), P(5), P(7) ); + // Q = alpha * L * C1 / C2, where + // + // alpha = sqrt( 1/32 ) + // L = max( L1, L2, L3, L4, D1, D2 ) + // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 ) + // 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. ); + double C2 = Min( anArea[ 0 ], + Min( anArea[ 1 ], + Min( anArea[ 2 ], anArea[ 3 ] ) ) ); + if ( C2 <= theEps ) + return theInf; + return alpha * L * C1 / C2; } + return 0; } double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const { // the aspect ratio is in the range [1.0,infinity] + // < 1.0 = very bad, zero area // 1.0 = good // infinity = bad - return Value / 1000.; + return ( Value < 0.9 ) ? 1000 : Value / 1000.; } SMDSAbs_ElementType AspectRatio::GetType() const @@ -464,10 +937,13 @@ SMDSAbs_ElementType AspectRatio::GetType() const } +//================================================================================ /* Class : AspectRatio3D Description : Functor for calculating aspect ratio */ +//================================================================================ + namespace{ inline double getHalfPerimeter(double theTria[3]){ @@ -476,9 +952,9 @@ namespace{ inline double getArea(double theHalfPerim, double theTria[3]){ return sqrt(theHalfPerim* - (theHalfPerim-theTria[0])* - (theHalfPerim-theTria[1])* - (theHalfPerim-theTria[2])); + (theHalfPerim-theTria[0])* + (theHalfPerim-theTria[1])* + (theHalfPerim-theTria[2])); } inline double getVolume(double theLen[6]){ @@ -530,6 +1006,28 @@ namespace{ } +double AspectRatio3D::GetValue( long theId ) +{ + double aVal = 0; + myCurrElement = myMesh->FindElement( theId ); + if ( myCurrElement && myCurrElement->GetVtkType() == VTK_TETRA ) + { + // 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() )) + aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell )); + } + else + { + TSequenceOfXYZ P; + if ( GetPoints( myCurrElement, P )) + aVal = Round( GetValue( P )); + } + return aVal; +} + double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) { double aQuality = 0.0; @@ -542,10 +1040,11 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) 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 return aQuality; } - switch(nbNodes){ + switch(nbNodes) { case 4:{ double aLen[6] = { getDistance(P( 1 ),P( 2 )), // a @@ -763,7 +1262,22 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) } break; } - } + 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); + } + { + gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality); + } + { + gp_XYZ aXYZ[8] = {P( 3 ),P( 4 ),P( 6 ),P( 1 ),P( 9 ),P( 10 ),P( 12 ),P( 7 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality); + } + break; + } // switch(nbNodes) + if ( nbNodes > 4 ) { // avaluate aspect ratio of quadranle faces AspectRatio aspect2D; @@ -796,10 +1310,13 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const } +//================================================================================ /* Class : Warping Description : Functor for calculating warping */ +//================================================================================ + double Warping::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 4 ) @@ -812,7 +1329,11 @@ double Warping::GetValue( const TSequenceOfXYZ& P ) double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G ); double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G ); - return Max( Max( A1, A2 ), Max( A3, A4 ) ); + double val = Max( Max( A1, A2 ), Max( A3, A4 ) ); + + const double eps = 0.1; // val is in degrees + + return val < eps ? 0. : val; } double Warping::ComputeA( const gp_XYZ& thePnt1, @@ -823,20 +1344,20 @@ double Warping::ComputeA( const gp_XYZ& thePnt1, double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) ); double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) ); double L = Min( aLen1, aLen2 ) * 0.5; - if ( L < Precision::Confusion()) - return 0.; + if ( L < theEps ) + return theInf; gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG; gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG; gp_XYZ N = GI.Crossed( GJ ); if ( N.Modulus() < gp::Resolution() ) - return PI / 2; + return M_PI / 2; N.Normalize(); double H = ( thePnt2 - theG ).Dot( N ); - return asin( fabs( H / L ) ) * 180. / PI; + return asin( fabs( H / L ) ) * 180. / M_PI; } double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -853,37 +1374,44 @@ SMDSAbs_ElementType Warping::GetType() const } +//================================================================================ /* Class : Taper Description : Functor for calculating taper */ +//================================================================================ + double Taper::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 4 ) return 0.; // Compute taper - double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2.; - double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2.; - double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2.; - double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2.; + double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ); + double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ); + double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ); + double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ); double JA = 0.25 * ( J1 + J2 + J3 + J4 ); - if ( JA <= Precision::Confusion() ) - return 0.; + if ( JA <= theEps ) + return theInf; double T1 = fabs( ( J1 - JA ) / JA ); double T2 = fabs( ( J2 - JA ) / JA ); double T3 = fabs( ( J3 - JA ) / JA ); double T4 = fabs( ( J4 - JA ) / JA ); - return Max( Max( T1, T2 ), Max( T3, T4 ) ); + double val = Max( Max( T1, T2 ), Max( T3, T4 ) ); + + const double eps = 0.01; + + return val < eps ? 0. : val; } double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const { // the taper is in the range [0.0,1.0] - // 0.0 = good (no taper) + // 0.0 = good (no taper) // 1.0 = bad (les cotes opposes sont allignes) return Value; } @@ -893,11 +1421,13 @@ SMDSAbs_ElementType Taper::GetType() const return SMDSAbs_Face; } - +//================================================================================ /* Class : Skew Description : Functor for calculating skew in degrees */ +//================================================================================ + static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 ) { gp_XYZ p12 = ( p2 + p1 ) / 2.; @@ -915,14 +1445,14 @@ double Skew::GetValue( const TSequenceOfXYZ& P ) return 0.; // Compute skew - static double PI2 = PI / 2.; + const double PI2 = M_PI / 2.; if ( P.size() == 3 ) { double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) ); double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) ); double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) ); - return Max( A0, Max( A1, A2 ) ) * 180. / PI; + return Max( A0, Max( A1, A2 ) ) * 180. / M_PI; } else { @@ -935,11 +1465,11 @@ double Skew::GetValue( const TSequenceOfXYZ& P ) double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution() ? 0. : fabs( PI2 - v1.Angle( v2 ) ); - //BUG SWP12743 - if ( A < Precision::Angular() ) - return 0.; + double val = A * 180. / M_PI; - return A * 180. / PI; + const double eps = 0.1; // val is in degrees + + return val < eps ? 0. : val; } } @@ -957,22 +1487,32 @@ SMDSAbs_ElementType Skew::GetType() const } +//================================================================================ /* Class : Area Description : Functor for calculating area */ +//================================================================================ + double Area::GetValue( const TSequenceOfXYZ& P ) { - gp_Vec aVec1( P(2) - P(1) ); - gp_Vec aVec2( P(3) - P(1) ); - gp_Vec SumVec = aVec1 ^ aVec2; - for (int i=4; i<=P.size(); i++) { - gp_Vec aVec1( P(i-1) - P(1) ); - gp_Vec aVec2( P(i) - P(1) ); - gp_Vec tmp = aVec1 ^ aVec2; - SumVec.Add(tmp); + double val = 0.0; + if ( P.size() > 2 ) + { + gp_Vec aVec1( P(2) - P(1) ); + gp_Vec aVec2( P(3) - P(1) ); + gp_Vec SumVec = aVec1 ^ aVec2; + + 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 tmp = aVec1 ^ aVec2; + SumVec.Add(tmp); + } + val = SumVec.Magnitude() * 0.5; } - return SumVec.Magnitude() * 0.5; + return val; } double Area::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -986,11 +1526,13 @@ SMDSAbs_ElementType Area::GetType() const return SMDSAbs_Face; } - +//================================================================================ /* Class : Length - Description : Functor for calculating length off edge + Description : Functor for calculating length of edge */ +//================================================================================ + double Length::GetValue( const TSequenceOfXYZ& P ) { switch ( P.size() ) { @@ -1011,210 +1553,252 @@ SMDSAbs_ElementType Length::GetType() const return SMDSAbs_Edge; } +//================================================================================ /* Class : Length2D - Description : Functor for calculating length of edge + Description : Functor for calculating minimal length of edge */ +//================================================================================ -double Length2D::GetValue( long theElementId) +double Length2D::GetValue( const TSequenceOfXYZ& P ) { - TSequenceOfXYZ P; - - //cout<<"Length2D::GetValue"<FindElement( theElementId ); - SMDSAbs_ElementType aType = aElem->GetType(); + double aVal = 0; + int len = P.size(); + SMDSAbs_EntityType aType = P.getElementEntity(); - int len = P.size(); - - switch (aType){ - case SMDSAbs_All: - case SMDSAbs_Node: - case SMDSAbs_Edge: - if (len == 2){ - aVal = getDistance( P( 1 ), P( 2 ) ); - break; - } - else if (len == 3){ // quadratic edge - aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 )); - break; - } - case SMDSAbs_Face: - if (len == 3){ // triangles - double L1 = getDistance(P( 1 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 1 )); - aVal = Max(L1,Max(L2,L3)); - break; - } - else if (len == 4){ // quadrangles - double L1 = getDistance(P( 1 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 4 )); - double L4 = getDistance(P( 4 ),P( 1 )); - aVal = Max(Max(L1,L2),Max(L3,L4)); - break; - } - if (len == 6){ // quadratic triangles - double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 )); - double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 )); - double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 )); - aVal = Max(L1,Max(L2,L3)); - //cout<<"L1="<FindElement( theId ); + if ( !anElem ) + return false; const SMDSAbs_ElementType anElemType = anElem->GetType(); - if ( !anElem || (myType != SMDSAbs_All && anElemType != myType) ) + if ( myType != SMDSAbs_All && anElemType != myType ) return false; - const int aNbNode = anElem->NbNodes(); - bool isOk = false; - switch( anElemType ) - { - case SMDSAbs_Node: - isOk = (myGeomType == SMDSGeom_POINT); - break; - - case SMDSAbs_Edge: - isOk = (myGeomType == SMDSGeom_EDGE); - break; - - case SMDSAbs_Face: - if ( myGeomType == SMDSGeom_TRIANGLE ) - isOk = (!anElem->IsPoly() && aNbNode == 3); - else if ( myGeomType == SMDSGeom_QUADRANGLE ) - isOk = (!anElem->IsPoly() && aNbNode == 4); - else if ( myGeomType == SMDSGeom_POLYGON ) - isOk = anElem->IsPoly(); - break; - - case SMDSAbs_Volume: - if ( myGeomType == SMDSGeom_TETRA ) - isOk = (!anElem->IsPoly() && aNbNode == 4); - else if ( myGeomType == SMDSGeom_PYRAMID ) - isOk = (!anElem->IsPoly() && aNbNode == 5); - else if ( myGeomType == SMDSGeom_PENTA ) - isOk = (!anElem->IsPoly() && aNbNode == 6); - else if ( myGeomType == SMDSGeom_HEXA ) - isOk = (!anElem->IsPoly() && aNbNode == 8); - else if ( myGeomType == SMDSGeom_POLYHEDRA ) - isOk = anElem->IsPoly(); - break; - default: break; - } + bool isOk = ( anElem->GetGeomType() == myGeomType ); return isOk; } @@ -1999,13 +2947,286 @@ SMDSAbs_GeometryType ElemGeomType::GetGeomType() const return myGeomType; } +//================================================================================ +/* + Class : ElemEntityType + Description : Predicate to check element entity type +*/ +//================================================================================ + +ElemEntityType::ElemEntityType(): + myMesh( 0 ), + myType( SMDSAbs_All ), + myEntityType( SMDSEntity_0D ) +{ +} + +void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; +} + +bool ElemEntityType::IsSatisfy( long theId ) +{ + if ( !myMesh ) return false; + if ( myType == SMDSAbs_Node ) + return myMesh->FindNode( theId ); + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + return ( anElem && + myEntityType == anElem->GetEntityType() ); +} + +void ElemEntityType::SetType( SMDSAbs_ElementType theType ) +{ + myType = theType; +} + +SMDSAbs_ElementType ElemEntityType::GetType() const +{ + return myType; +} + +void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType ) +{ + myEntityType = theEntityType; +} + +SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const +{ + return myEntityType; +} + +//================================================================================ +/*! + * \brief Class ConnectedElements + */ +//================================================================================ + +ConnectedElements::ConnectedElements(): + myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {} + +SMDSAbs_ElementType ConnectedElements::GetType() const +{ return myType; } + +int ConnectedElements::GetNode() const +{ return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ + +std::vector ConnectedElements::GetPoint() const +{ return myXYZ; } + +void ConnectedElements::clearOkIDs() +{ myOkIDsReady = false; myOkIDs.clear(); } + +void ConnectedElements::SetType( SMDSAbs_ElementType theType ) +{ + if ( myType != theType || myMeshModifTracer.IsMeshModified() ) + clearOkIDs(); + myType = theType; +} + +void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified() ) + { + clearOkIDs(); + if ( !myXYZ.empty() ) + SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh + } +} + +void ConnectedElements::SetNode( int nodeID ) +{ + myNodeID = nodeID; + myXYZ.clear(); + + bool isSameDomain = false; + if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() ) + if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID )) + { + SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType ); + while ( !isSameDomain && eIt->more() ) + isSameDomain = IsSatisfy( eIt->next()->GetID() ); + } + if ( !isSameDomain ) + clearOkIDs(); +} + +void ConnectedElements::SetPoint( double x, double y, double z ) +{ + myXYZ.resize(3); + myXYZ[0] = x; + myXYZ[1] = y; + myXYZ[2] = z; + myNodeID = 0; + + bool isSameDomain = false; + + // find myNodeID by myXYZ if possible + if ( myMeshModifTracer.GetMesh() ) + { + SMESHUtils::Deleter searcher + ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() )); + + std::vector< const SMDS_MeshElement* > foundElems; + searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems ); + + if ( !foundElems.empty() ) + { + myNodeID = foundElems[0]->GetNode(0)->GetID(); + if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() ) + isSameDomain = IsSatisfy( foundElems[0]->GetID() ); + } + } + if ( !isSameDomain ) + clearOkIDs(); +} + +bool ConnectedElements::IsSatisfy( long theElementId ) +{ + // Here we do NOT check if the mesh has changed, we do it in Set...() only!!! + + if ( !myOkIDsReady ) + { + if ( !myMeshModifTracer.GetMesh() ) + return false; + const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID ); + if ( !node0 ) + return false; + + std::list< const SMDS_MeshNode* > nodeQueue( 1, node0 ); + std::set< int > checkedNodeIDs; + // algo: + // foreach node in nodeQueue: + // foreach element sharing a node: + // add ID of an element of myType to myOkIDs; + // push all element nodes absent from checkedNodeIDs to nodeQueue; + while ( !nodeQueue.empty() ) + { + const SMDS_MeshNode* node = nodeQueue.front(); + nodeQueue.pop_front(); + + // loop on elements sharing the node + SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(); + while ( eIt->more() ) + { + // keep elements of myType + const SMDS_MeshElement* element = eIt->next(); + if ( element->GetType() == myType ) + myOkIDs.insert( myOkIDs.end(), element->GetID() ); + + // enqueue nodes of the element + SMDS_ElemIteratorPtr nIt = element->nodesIterator(); + while ( nIt->more() ) + { + const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() ); + if ( checkedNodeIDs.insert( n->GetID() ).second ) + nodeQueue.push_back( n ); + } + } + } + if ( myType == SMDSAbs_Node ) + std::swap( myOkIDs, checkedNodeIDs ); + + size_t totalNbElems = myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ); + if ( myOkIDs.size() == totalNbElems ) + myOkIDs.clear(); + + myOkIDsReady = true; + } + + return myOkIDs.empty() ? true : myOkIDs.count( theElementId ); +} + +//================================================================================ +/*! + * \brief Class CoplanarFaces + */ +//================================================================================ + +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) +{ +} +void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified() ) + { + // Build a set of coplanar face ids + + myCoplanarIDs.Clear(); + + if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler ) + return; + + const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID ); + if ( !face || face->GetType() != SMDSAbs_Face ) + return; + + bool normOK; + gp_Vec myNorm = getNormale( static_cast(face), &normOK ); + if (!normOK) + return; + + const double cosTol = Cos( myToler * M_PI / 180. ); + NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks; + + 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; + myNorm = faceQueue.front().second; + faceQueue.pop_front(); + + for ( int i = 0, nbN = face->NbCornerNodes(); i < nbN; ++i ) + { + const SMDS_MeshNode* n1 = face->GetNode( i ); + const SMDS_MeshNode* n2 = face->GetNode(( i+1 )%nbN); + if ( !checkedLinks.Add( SMESH_TLink( n1, n2 ))) + continue; + SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() ) + { + const SMDS_MeshElement* f = fIt->next(); + if ( f->GetNodeIndex( n2 ) > -1 ) + { + gp_Vec norm = getNormale( static_cast(f), &normOK ); + if (!normOK || isLessAngle( myNorm, norm, cosTol)) + { + myCoplanarIDs.Add( f->GetID() ); + faceQueue.push_back( std::make_pair( f, norm )); + } + } + } + } + } + } +} +bool CoplanarFaces::IsSatisfy( long theElementId ) +{ + return myCoplanarIDs.Contains( theElementId ); +} + /* - Class : RangeOfIds - Description : Predicate for Range of Ids. - Range may be specified with two ways. - 1. Using AddToRange method - 2. With SetRangeStr method. Parameter of this method is a string - like as "1,2,3,50-60,63,67,70-" + *Class : RangeOfIds + *Description : Predicate for Range of Ids. + * Range may be specified with two ways. + * 1. Using AddToRange method + * 2. With SetRangeStr method. Parameter of this method is a string + * like as "1,2,3,50-60,63,67,70-" */ //======================================================================= @@ -2124,11 +3345,13 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr ) myIds.Clear(); TCollection_AsciiString aStr = theStr; + for ( int i = 1; i <= aStr.Length(); ++i ) + { + char c = aStr.Value( i ); + if ( !isdigit( c ) && c != ',' && c != '-' ) + aStr.SetValue( i, ','); + } aStr.RemoveAll( ' ' ); - aStr.RemoveAll( '\t' ); - - for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) ) - aStr.Remove( aPos, 2 ); TCollection_AsciiString tmpStr = aStr.Token( ",", 1 ); int i = 1; @@ -2152,8 +3375,8 @@ bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr ) while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' ); while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' ); - if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() || - !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() ) + if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) || + (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) ) return false; myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() ); @@ -2199,7 +3422,7 @@ bool RangeOfIds::IsSatisfy( long theId ) else { const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); - if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All ) + if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All )) return false; } @@ -2390,8 +3613,8 @@ bool LogicalOR::IsSatisfy( long theId ) return myPredicate1 && myPredicate2 && - myPredicate1->IsSatisfy( theId ) || - myPredicate2->IsSatisfy( theId ); + (myPredicate1->IsSatisfy( theId ) || + myPredicate2->IsSatisfy( theId )); } @@ -2399,6 +3622,31 @@ bool LogicalOR::IsSatisfy( long theId ) FILTER */ +// #ifdef WITH_TBB +// #include +// #include + +// namespace Parallel +// { +// typedef tbb::enumerable_thread_specific< TIdSequence > TIdSeq; + +// struct Predicate +// { +// const SMDS_Mesh* myMesh; +// PredicatePtr myPredicate; +// TIdSeq & myOKIds; +// Predicate( const SMDS_Mesh* m, PredicatePtr p, TIdSeq & ids ): +// myMesh(m), myPredicate(p->Duplicate()), myOKIds(ids) {} +// void operator() ( const tbb::blocked_range& r ) const +// { +// for ( size_t i = r.begin(); i != r.end(); ++i ) +// if ( myPredicate->IsSatisfy( i )) +// myOKIds.local().push_back(); +// } +// } +// } +// #endif + Filter::Filter() {} @@ -2410,26 +3658,9 @@ void Filter::SetPredicate( PredicatePtr thePredicate ) myPredicate = thePredicate; } -template -inline void FillSequence(const TIterator& theIterator, - TPredicate& thePredicate, - Filter::TIdSequence& theSequence) -{ - if ( theIterator ) { - while( theIterator->more() ) { - TElement anElem = theIterator->next(); - long anId = anElem->GetID(); - if ( thePredicate->IsSatisfy( anId ) ) - theSequence.push_back( anId ); - } - } -} - -void -Filter:: -GetElementsId( const SMDS_Mesh* theMesh, - PredicatePtr thePredicate, - TIdSequence& theSequence ) +void Filter::GetElementsId( const SMDS_Mesh* theMesh, + PredicatePtr thePredicate, + TIdSequence& theSequence ) { theSequence.clear(); @@ -2438,31 +3669,19 @@ GetElementsId( const SMDS_Mesh* theMesh, thePredicate->SetMesh( theMesh ); - SMDSAbs_ElementType aType = thePredicate->GetType(); - switch(aType){ - case SMDSAbs_Node: - FillSequence(theMesh->nodesIterator(),thePredicate,theSequence); - break; - case SMDSAbs_Edge: - FillSequence(theMesh->edgesIterator(),thePredicate,theSequence); - break; - case SMDSAbs_Face: - FillSequence(theMesh->facesIterator(),thePredicate,theSequence); - break; - case SMDSAbs_Volume: - FillSequence(theMesh->volumesIterator(),thePredicate,theSequence); - break; - case SMDSAbs_All: - FillSequence(theMesh->edgesIterator(),thePredicate,theSequence); - FillSequence(theMesh->facesIterator(),thePredicate,theSequence); - FillSequence(theMesh->volumesIterator(),thePredicate,theSequence); - break; + 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 ); + } } } -void -Filter::GetElementsId( const SMDS_Mesh* theMesh, - Filter::TIdSequence& theSequence ) +void Filter::GetElementsId( const SMDS_Mesh* theMesh, + Filter::TIdSequence& theSequence ) { GetElementsId(theMesh,myPredicate,theSequence); } @@ -2608,7 +3827,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(); @@ -2635,32 +3854,6 @@ static void getLinks( const SMDS_MeshFace* theFace, } } -static gp_XYZ getNormale( const SMDS_MeshFace* theFace ) -{ - gp_XYZ n; - int aNbNode = theFace->NbNodes(); - TColgp_Array1OfXYZ anArrOfXYZ(1,4); - SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator(); - int i = 1; - for ( ; aNodeItr->more() && i <= 4; i++ ) { - SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); - anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); - } - - gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1); - gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1); - n = q1 ^ q2; - if ( aNbNode > 3 ) { - gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1); - n += q2 ^ q3; - } - double len = n.Modulus(); - if ( len > 0 ) - n /= len; - - return n; -} - bool ManifoldPart::findConnected ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt, SMDS_MeshFace* theStartFace, @@ -2815,7 +4008,7 @@ void ManifoldPart::expandBoundary void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink, ManifoldPart::TVectorOfFacePtr& theFaces ) const { - SMDS_Mesh::SetOfFaces aSetOfFaces; + std::set aSetOfFaces; // take all faces that shared first node SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator(); for ( ; anItr->more(); ) @@ -2823,7 +4016,7 @@ void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink, SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next(); if ( !aFace ) continue; - aSetOfFaces.Add( aFace ); + aSetOfFaces.insert( aFace ); } // take all faces that shared second node anItr = theLink.myNode2->facesIterator(); @@ -2831,19 +4024,68 @@ void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink, for ( ; anItr->more(); ) { SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next(); - if ( aSetOfFaces.Contains( aFace ) ) + if ( aSetOfFaces.count( aFace ) ) theFaces.push_back( aFace ); } } - /* - ElementsOnSurface + Class : BelongToMeshGroup + Description : Verify whether a mesh element is included into a mesh group */ +BelongToMeshGroup::BelongToMeshGroup(): myGroup( 0 ) +{ +} + +void BelongToMeshGroup::SetGroup( SMESHDS_GroupBase* g ) +{ + myGroup = g; +} + +void BelongToMeshGroup::SetStoreName( const std::string& sn ) +{ + myStoreName = sn; +} + +void BelongToMeshGroup::SetMesh( const SMDS_Mesh* theMesh ) +{ + if ( myGroup && myGroup->GetMesh() != theMesh ) + { + myGroup = 0; + } + if ( !myGroup && !myStoreName.empty() ) + { + if ( const SMESHDS_Mesh* aMesh = dynamic_cast(theMesh)) + { + const std::set& grps = aMesh->GetGroups(); + std::set::const_iterator g = grps.begin(); + for ( ; g != grps.end() && !myGroup; ++g ) + if ( *g && myStoreName == (*g)->GetStoreName() ) + myGroup = *g; + } + } + if ( myGroup ) + { + myGroup->IsEmpty(); // make GroupOnFilter update its predicate + } +} + +bool BelongToMeshGroup::IsSatisfy( long theElementId ) +{ + return myGroup ? myGroup->Contains( theElementId ) : false; +} + +SMDSAbs_ElementType BelongToMeshGroup::GetType() const +{ + return myGroup ? myGroup->GetType() : SMDSAbs_All; +} + +//================================================================================ +// ElementsOnSurface +//================================================================================ ElementsOnSurface::ElementsOnSurface() { - myMesh = 0; myIds.Clear(); myType = SMDSAbs_All; mySurf.Nullify(); @@ -2853,15 +4095,13 @@ ElementsOnSurface::ElementsOnSurface() ElementsOnSurface::~ElementsOnSurface() { - myMesh = 0; } void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh ) { - if ( myMesh == theMesh ) - return; - myMesh = theMesh; - process(); + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified()) + process(); } bool ElementsOnSurface::IsSatisfy( long theElementId ) @@ -2916,32 +4156,14 @@ void ElementsOnSurface::process() if ( mySurf.IsNull() ) return; - if ( myMesh == 0 ) + if ( !myMeshModifTracer.GetMesh() ) return; - if ( myType == SMDSAbs_Face || myType == SMDSAbs_All ) - { - myIds.ReSize( myMesh->NbFaces() ); - SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); - for(; anIter->more(); ) - process( anIter->next() ); - } - - if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All ) - { - myIds.ReSize( myIds.Extent() + myMesh->NbEdges() ); - SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator(); - for(; anIter->more(); ) - process( anIter->next() ); - } + myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType )); - if ( myType == SMDSAbs_Node ) - { - myIds.ReSize( myMesh->NbNodes() ); - SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator(); - for(; anIter->more(); ) - process( anIter->next() ); - } + SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType ); + for(; anIter->more(); ) + process( anIter->next() ); } void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr ) @@ -2994,34 +4216,96 @@ bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) } -/* - ElementsOnShape -*/ +//================================================================================ +// ElementsOnShape +//================================================================================ -ElementsOnShape::ElementsOnShape() - : myMesh(0), - myType(SMDSAbs_All), - myToler(Precision::Confusion()), - myAllNodesFlag(false) -{ - myCurShapeType = TopAbs_SHAPE; +namespace { + const int theIsCheckedFlag = 0x0000100; } -ElementsOnShape::~ElementsOnShape() +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) { } -void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh) +ElementsOnShape::~ElementsOnShape() { - if (myMesh != theMesh) { - myMesh = theMesh; - SetShape(myShape, myType); - } + clearClassifiers(); } -bool ElementsOnShape::IsSatisfy (long theElementId) +Predicate* ElementsOnShape::clone() const { - return myIds.Contains(theElementId); + 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 @@ -3044,198 +4328,878 @@ double ElementsOnShape::GetTolerance() const void ElementsOnShape::SetAllNodes (bool theAllNodes) { - if (myAllNodesFlag != theAllNodes) { - myAllNodesFlag = theAllNodes; - SetShape(myShape, myType); + myAllNodesFlag = theAllNodes; +} + +void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh) +{ + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified()) + { + size_t nbNodes = theMesh ? theMesh->NbNodes() : 0; + if ( myNodeIsChecked.size() == nbNodes ) + { + std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false ); + } + else + { + SMESHUtils::FreeVector( myNodeIsChecked ); + SMESHUtils::FreeVector( myNodeIsOut ); + myNodeIsChecked.resize( nbNodes, false ); + myNodeIsOut.resize( nbNodes ); + } } } -void ElementsOnShape::SetShape (const TopoDS_Shape& theShape, - const SMDSAbs_ElementType theType) +bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut ) { - myType = theType; - myShape = theShape; - myIds.Clear(); + if ( n->GetID() >= (int) myNodeIsChecked.size() || + !myNodeIsChecked[ n->GetID() ]) + return false; - if (myMesh == 0) return; + isOut = myNodeIsOut[ n->GetID() ]; + return true; +} - switch (myType) +void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool isOut ) +{ + if ( n->GetID() < (int) myNodeIsChecked.size() ) { - case SMDSAbs_All: - myIds.ReSize(myMesh->NbEdges() + myMesh->NbFaces() + myMesh->NbVolumes()); - break; - case SMDSAbs_Node: - myIds.ReSize(myMesh->NbNodes()); - break; - case SMDSAbs_Edge: - myIds.ReSize(myMesh->NbEdges()); - break; - case SMDSAbs_Face: - myIds.ReSize(myMesh->NbFaces()); - break; - case SMDSAbs_Volume: - myIds.ReSize(myMesh->NbVolumes()); - break; - default: - break; + myNodeIsChecked[ n->GetID() ] = true; + myNodeIsOut [ n->GetID() ] = isOut; } +} - myShapesMap.Clear(); - addShape(myShape); +void ElementsOnShape::SetShape (const TopoDS_Shape& theShape, + const SMDSAbs_ElementType theType) +{ + bool shapeChanges = ( myShape != theShape ); + myType = theType; + myShape = theShape; + if ( myShape.IsNull() ) return; + + if ( shapeChanges ) + { + // find most complex shapes + TopTools_IndexedMapOfShape shapesMap; + TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }; + TopExp_Explorer sub; + for ( int i = 0; i < 4; ++i ) + { + if ( shapesMap.IsEmpty() ) + for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() ) + shapesMap.Add( sub.Current() ); + if ( i > 0 ) + for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() ) + shapesMap.Add( sub.Current() ); + } + + clearClassifiers(); + myClassifiers.resize( shapesMap.Extent() ); + for ( int i = 0; i < shapesMap.Extent(); ++i ) + myClassifiers[ i ].Init( shapesMap( i+1 ), myToler ); + } + + if ( theType == SMDSAbs_Node ) + { + SMESHUtils::FreeVector( myNodeIsChecked ); + SMESHUtils::FreeVector( myNodeIsOut ); + } + else + { + std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false ); + } } -void ElementsOnShape::addShape (const TopoDS_Shape& theShape) +void ElementsOnShape::clearClassifiers() { - if (theShape.IsNull() || myMesh == 0) - return; + // for ( size_t i = 0; i < myClassifiers.size(); ++i ) + // delete myClassifiers[ i ]; + myClassifiers.clear(); + + delete myOctree; + myOctree = 0; +} + +bool ElementsOnShape::IsSatisfy( long elemId ) +{ + if ( myClassifiers.empty() ) + return false; + + const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh(); + if ( myType == SMDSAbs_Node ) + return IsSatisfy( mesh->FindNode( elemId )); + return IsSatisfy( mesh->FindElement( elemId )); +} + +bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem) +{ + if ( !elem ) + return false; - if (!myShapesMap.Add(theShape)) return; + bool isSatisfy = myAllNodesFlag, isNodeOut; - myCurShapeType = theShape.ShapeType(); - switch (myCurShapeType) + 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)) { - case TopAbs_COMPOUND: - case TopAbs_COMPSOLID: - case TopAbs_SHELL: - case TopAbs_WIRE: + SMESH_TNodeXYZ aPnt( aNodeItr->next() ); + centerXYZ += aPnt; + + isNodeOut = true; + if ( !getNodeIsOut( aPnt._node, isNodeOut )) { - TopoDS_Iterator anIt (theShape, Standard_True, Standard_True); - for (; anIt.More(); anIt.Next()) addShape(anIt.Value()); + 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 ); } - break; + isSatisfy = !isNodeOut; + } + + // Check the center point for volumes MantisBug 0020168 + if ( isSatisfy && + myAllNodesFlag && + myClassifiers[0].ShapeType() == TopAbs_SOLID ) + { + centerXYZ /= elem->NbNodes(); + isSatisfy = false; + 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; +} + +bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node, + TopoDS_Shape* okShape) +{ + if ( !node ) + return false; + + if ( !myOctree && myClassifiers.size() > 5 ) + { + myWorkClassifiers.resize( myClassifiers.size() ); + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + myWorkClassifiers[ i ] = & myClassifiers[ i ]; + myOctree = new OctreeClassifier( myWorkClassifiers ); + } + + bool isNodeOut = true; + + if ( okShape || !getNodeIsOut( node, isNodeOut )) + { + SMESH_NodeXYZ aPnt = node; + if ( myOctree ) + { + myWorkClassifiers.clear(); + myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers ); + + for ( size_t i = 0; i < myWorkClassifiers.size(); ++i ) + myWorkClassifiers[i]->SetChecked( false ); + + for ( size_t i = 0; i < myWorkClassifiers.size(); ++i ) + if ( !myWorkClassifiers[i]->IsChecked() && + !myWorkClassifiers[i]->IsOut( aPnt )) + { + isNodeOut = false; + if ( okShape ) + *okShape = myWorkClassifiers[i]->Shape(); + break; + } + } + else + { + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + if ( !myClassifiers[i].IsOut( aPnt )) + { + isNodeOut = false; + if ( okShape ) + *okShape = myWorkClassifiers[i]->Shape(); + break; + } + } + setNodeIsOut( node, isNodeOut ); + } + + return !isNodeOut; +} + +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 (( isShapeBox = isBox( theShape ))) { - myCurSC.Load(theShape); - process(); + myIsOutFun = & ElementsOnShape::Classifier::isOutOfBox; } - break; - case TopAbs_FACE: + else { - TopoDS_Face aFace = TopoDS::Face(theShape); - BRepAdaptor_Surface SA (aFace, true); - Standard_Real - u1 = SA.FirstUParameter(), - u2 = SA.LastUParameter(), - v1 = SA.FirstVParameter(), - v2 = SA.LastVParameter(); - Handle(Geom_Surface) surf = BRep_Tool::Surface(aFace); - myCurProjFace.Init(surf, u1,u2, v1,v2); - myCurFace = aFace; - process(); + mySolidClfr = new BRepClass3d_SolidClassifier(theShape); + myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid; } break; + } + 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::Classifier::isOutOfFace; + break; + } case TopAbs_EDGE: - { - TopoDS_Edge anEdge = TopoDS::Edge(theShape); - Standard_Real u1, u2; - Handle(Geom_Curve) curve = BRep_Tool::Curve(anEdge, u1, u2); - myCurProjEdge.Init(curve, u1, u2); - process(); - } + { + 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; break; + } case TopAbs_VERTEX: - { - TopoDS_Vertex aV = TopoDS::Vertex(theShape); - myCurPnt = BRep_Tool::Pnt(aV); - process(); - } + { + myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) ); + myIsOutFun = & ElementsOnShape::Classifier::isOutOfVertex; break; + } default: - break; + 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 ); + } } } -void ElementsOnShape::process() +ElementsOnShape::Classifier::~Classifier() { - if (myShape.IsNull() || myMesh == 0) - return; + delete mySolidClfr; mySolidClfr = 0; +} - if (myType == SMDSAbs_Node) +bool ElementsOnShape::Classifier::isOutOfSolid (const gp_Pnt& p) +{ + mySolidClfr->Perform( p, myTol ); + return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON ); +} + +bool ElementsOnShape::Classifier::isOutOfBox (const gp_Pnt& p) +{ + return myBox.IsOut( p.XYZ() ); +} + +bool ElementsOnShape::Classifier::isOutOfFace (const gp_Pnt& p) +{ + myProjFace.Perform( p ); + if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol ) { - SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator(); - while (anIter->more()) - process(anIter->next()); + // check relatively to the face + Standard_Real u, v; + myProjFace.LowerDistanceParameters(u, v); + gp_Pnt2d aProjPnt (u, v); + BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol ); + if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON ) + return false; } - else + return true; +} + +bool ElementsOnShape::Classifier::isOutOfEdge (const gp_Pnt& p) +{ + myProjEdge.Perform( p ); + return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol ); +} + +bool ElementsOnShape::Classifier::isOutOfVertex(const gp_Pnt& p) +{ + return ( myVertexXYZ.Distance( p ) > myTol ); +} + +bool ElementsOnShape::Classifier::isBox (const TopoDS_Shape& theShape) +{ + TopTools_IndexedMapOfShape vMap; + TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap ); + if ( vMap.Extent() != 8 ) + return false; + + myBox.Clear(); + for ( int i = 1; i <= 8; ++i ) + myBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))).XYZ() ); + + gp_XYZ pMin = myBox.CornerMin(), pMax = myBox.CornerMax(); + for ( int i = 1; i <= 8; ++i ) { - if (myType == SMDSAbs_Edge || myType == SMDSAbs_All) + gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))); + for ( int iC = 1; iC <= 3; ++ iC ) { - SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator(); - while (anIter->more()) - process(anIter->next()); + double d1 = Abs( pMin.Coord( iC ) - p.Coord( iC )); + double d2 = Abs( pMax.Coord( iC ) - p.Coord( iC )); + if ( Min( d1, d2 ) > myTol ) + return false; } + } + myBox.Enlarge( myTol ); + return true; +} - if (myType == SMDSAbs_Face || myType == SMDSAbs_All) - { - SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); - while (anIter->more()) { - process(anIter->next()); - } - } +ElementsOnShape:: +OctreeClassifier::OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers ) + :SMESH_Octree( new SMESH_TreeLimit ) +{ + myClassifiers = classifiers; + compute(); +} - if (myType == SMDSAbs_Volume || myType == SMDSAbs_All) +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 ) { - SMDS_VolumeIteratorPtr anIter = myMesh->volumesIterator(); - while (anIter->more()) - process(anIter->next()); + 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::process (const SMDS_MeshElement* theElemPtr) +void ElementsOnShape:: +OctreeClassifier::GetClassifiersAtPoint( const gp_XYZ& point, + std::vector< ElementsOnShape::Classifier* >& result ) { - if (myShape.IsNull()) + if ( getBox()->IsOut( point )) return; - SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator(); - bool isSatisfy = myAllNodesFlag; - - while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) + if ( isLeaf() ) { - SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); - gp_Pnt aPnt (aNode->X(), aNode->Y(), aNode->Z()); + 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 }; - switch (myCurShapeType) + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + { + for ( int j = 0; j < nbChildren(); j++ ) { - case TopAbs_SOLID: + if ( !myClassifiers[i]->GetBndBox()->IsOut( *myChildren[j]->getBox() )) { - myCurSC.Perform(aPnt, myToler); - isSatisfy = (myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON); + myClassifiers[i]->SetFlag( childFlag[ j ]); + ++nbInChild[ j ]; } - break; - case TopAbs_FACE: + } + } + + 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 ])) { - myCurProjFace.Perform(aPnt); - isSatisfy = (myCurProjFace.IsDone() && myCurProjFace.LowerDistance() <= myToler); - if (isSatisfy) - { - // check relatively the face - Quantity_Parameter u, v; - myCurProjFace.LowerDistanceParameters(u, v); - gp_Pnt2d aProjPnt (u, v); - BRepClass_FaceClassifier aClsf (myCurFace, aProjPnt, myToler); - isSatisfy = (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON); - } + --nbInChild[ j ]; + child->myClassifiers[ nbInChild[ j ]] = myClassifiers[ i ]; + myClassifiers[ i ]->UnsetFlag( childFlag[ j ]); } - break; - case TopAbs_EDGE: - { - myCurProjEdge.Perform(aPnt); - isSatisfy = (myCurProjEdge.NbPoints() > 0 && myCurProjEdge.LowerDistance() <= myToler); + } + } + 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 + Description : Predicate for verifying whether entity belongs to + specified geometrical support +*/ + +BelongToGeom::BelongToGeom() + : myMeshDS(NULL), + 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 ) +{ + if ( myMeshDS != theMesh ) + { + myMeshDS = dynamic_cast(theMesh); + init(); + } +} + +void BelongToGeom::SetGeom( const TopoDS_Shape& theShape ) +{ + if ( myShape != theShape ) + { + myShape = theShape; + init(); + } +} + +static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap, + const TopoDS_Shape& theShape) +{ + if (theMap.Contains(theShape)) return true; + + if (theShape.ShapeType() == TopAbs_COMPOUND || + theShape.ShapeType() == TopAbs_COMPSOLID) + { + TopoDS_Iterator anIt (theShape, Standard_True, Standard_True); + for (; anIt.More(); anIt.Next()) + { + if (!IsSubShape(theMap, anIt.Value())) { + return false; } - break; - case TopAbs_VERTEX: + } + return true; + } + + return false; +} + +void BelongToGeom::init() +{ + if ( !myMeshDS || myShape.IsNull() ) return; + + // is sub-shape of main shape? + TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh(); + if (aMainShape.IsNull()) { + myIsSubshape = false; + } + else { + TopTools_IndexedMapOfShape aMap; + TopExp::MapShapes( aMainShape, aMap ); + myIsSubshape = IsSubShape( aMap, myShape ); + if ( myIsSubshape ) + { + aMap.Clear(); + TopExp::MapShapes( myShape, aMap ); + mySubShapesIDs.Clear(); + for ( int i = 1; i <= aMap.Extent(); ++i ) { - isSatisfy = (aPnt.Distance(myCurPnt) <= myToler); + int subID = myMeshDS->ShapeToIndex( aMap( i )); + if ( subID > 0 ) + mySubShapesIDs.Add( subID ); } - break; - default: + } + } + + //if (!myIsSubshape) // to be always ready to check an element not bound to geometry + { + 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 ); + } +} + +bool BelongToGeom::IsSatisfy (long theId) +{ + if (myMeshDS == 0 || myShape.IsNull()) + return false; + + if (!myIsSubshape) + { + return myElementsOnShapePtr->IsSatisfy(theId); + } + + // Case of sub-mesh + + if (myType == SMDSAbs_Node) + { + if ( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId )) + { + if ( aNode->getshapeId() < 1 ) + return myElementsOnShapePtr->IsSatisfy(theId); + else + return mySubShapesIDs.Contains( aNode->getshapeId() ); + } + } + else + { + if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId )) + { + if ( anElem->GetType() == myType ) { - isSatisfy = false; + if ( anElem->getshapeId() < 1 ) + return myElementsOnShapePtr->IsSatisfy(theId); + else + return mySubShapesIDs.Contains( anElem->getshapeId() ); } } } - if (isSatisfy) - myIds.Add(theElemPtr->GetID()); + return false; +} + +void BelongToGeom::SetType (SMDSAbs_ElementType theType) +{ + if ( myType != theType ) + { + myType = theType; + init(); + } +} + +SMDSAbs_ElementType BelongToGeom::GetType() const +{ + return myType; +} + +TopoDS_Shape BelongToGeom::GetShape() +{ + return myShape; +} + +const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const +{ + return myMeshDS; +} + +void BelongToGeom::SetTolerance (double theTolerance) +{ + myTolerance = theTolerance; + init(); +} + +double BelongToGeom::GetTolerance() +{ + return myTolerance; +} + +/* + Class : LyingOnGeom + Description : Predicate for verifying whether entiy lying or partially lying on + specified geometrical support +*/ + +LyingOnGeom::LyingOnGeom() + : myMeshDS(NULL), + myType(SMDSAbs_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 ) +{ + if ( myMeshDS != theMesh ) + { + myMeshDS = dynamic_cast(theMesh); + init(); + } +} + +void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape ) +{ + if ( myShape != theShape ) + { + myShape = theShape; + init(); + } +} + +void LyingOnGeom::init() +{ + if (!myMeshDS || myShape.IsNull()) return; + + // is sub-shape of main shape? + TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh(); + if (aMainShape.IsNull()) { + myIsSubshape = false; + } + else { + myIsSubshape = myMeshDS->IsGroupOfSubShapes( myShape ); + } + + if (myIsSubshape) + { + TopTools_IndexedMapOfShape shapes; + TopExp::MapShapes( myShape, shapes ); + mySubShapesIDs.Clear(); + for ( int i = 1; i <= shapes.Extent(); ++i ) + { + int subID = myMeshDS->ShapeToIndex( shapes( i )); + if ( subID > 0 ) + mySubShapesIDs.Add( subID ); + } + } + // else // to be always ready to check an element not bound to geometry + { + 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 ); + } +} + +bool LyingOnGeom::IsSatisfy( long theId ) +{ + if ( myMeshDS == 0 || myShape.IsNull() ) + return false; + + if (!myIsSubshape) + { + return myElementsOnShapePtr->IsSatisfy(theId); + } + + // Case of sub-mesh + + const SMDS_MeshElement* elem = + ( myType == SMDSAbs_Node ) ? myMeshDS->FindNode( theId ) : myMeshDS->FindElement( theId ); + + if ( mySubShapesIDs.Contains( elem->getshapeId() )) + return true; + + if ( elem->GetType() != SMDSAbs_Node && elem->GetType() == myType ) + { + SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator(); + while ( nodeItr->more() ) + { + const SMDS_MeshElement* aNode = nodeItr->next(); + if ( mySubShapesIDs.Contains( aNode->getshapeId() )) + return true; + } + } + + return false; +} + +void LyingOnGeom::SetType( SMDSAbs_ElementType theType ) +{ + if ( myType != theType ) + { + myType = theType; + init(); + } +} + +SMDSAbs_ElementType LyingOnGeom::GetType() const +{ + return myType; +} + +TopoDS_Shape LyingOnGeom::GetShape() +{ + return myShape; +} + +const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const +{ + return myMeshDS; +} + +void LyingOnGeom::SetTolerance (double theTolerance) +{ + myTolerance = theTolerance; + init(); +} + +double LyingOnGeom::GetTolerance() +{ + return myTolerance; +} + +TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0) +{} + +TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n), myElem(0) +{} + +TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t), myElem(0) +{} + +TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray), myElem(theSequenceOfXYZ.myElem) +{} + +template +TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd), myElem(0) +{} + +TSequenceOfXYZ::~TSequenceOfXYZ() +{} + +TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ) +{ + myArray = theSequenceOfXYZ.myArray; + myElem = theSequenceOfXYZ.myElem; + return *this; +} + +gp_XYZ& TSequenceOfXYZ::operator()(size_type n) +{ + return myArray[n-1]; +} + +const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const +{ + return myArray[n-1]; +} + +void TSequenceOfXYZ::clear() +{ + myArray.clear(); +} + +void TSequenceOfXYZ::reserve(size_type n) +{ + myArray.reserve(n); +} + +void TSequenceOfXYZ::push_back(const gp_XYZ& v) +{ + myArray.push_back(v); +} + +TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const +{ + return myArray.size(); +} + +SMDSAbs_EntityType TSequenceOfXYZ::getElementEntity() const +{ + return myElem ? myElem->GetEntityType() : SMDSEntity_Last; +} + +TMeshModifTracer::TMeshModifTracer(): + myMeshModifTime(0), myMesh(0) +{ +} +void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh ) +{ + if ( theMesh != myMesh ) + myMeshModifTime = 0; + myMesh = theMesh; +} +bool TMeshModifTracer::IsMeshModified() +{ + bool modified = false; + if ( myMesh ) + { + modified = ( myMeshModifTime != myMesh->GetMTime() ); + myMeshModifTime = myMesh->GetMTime(); + } + return modified; }