X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=769096df4344b7e69a9d61ac59e8d7d6d903b5f9;hp=717f8e66ae41e82d0d06fdeef300d432aca0ff2d;hb=07589de92eb5d8e7d4e3bfde8b26f0d69251acf0;hpb=6b471bcc54cbeb90f0d977323db8c76d3d2cce09 diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 717f8e66a..769096df4 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -1,58 +1,89 @@ -// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, -// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE // -// 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. +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS // -// 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 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. // -// 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 +// 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 +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -// See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org #include "SMESH_ControlsDef.hxx" -#include +#include "SMDS_BallElement.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 +#include -#include "SMDS_Mesh.hxx" -#include "SMDS_Iterator.hxx" -#include "SMDS_MeshElement.hxx" -#include "SMDS_MeshNode.hxx" -#include "SMDS_VolumeTool.hxx" +#include +#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 ); @@ -87,36 +118,70 @@ namespace{ return 0; const SMDS_MeshElement* anEdge = theMesh->FindElement( theId ); - if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 ) + if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */) return 0; - TColStd_MapOfInteger aMap; - - int aResult = 0; - SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator(); - if ( anIter != 0 ) { - while( anIter->more() ) { - const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); - if ( aNode == 0 ) - return 0; - SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); - while( anElemIter->more() ) { - const SMDS_MeshElement* anElem = anElemIter->next(); - if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) { - int anId = anElem->GetID(); - - if ( anIter->more() ) // i.e. first node - aMap.Add( anId ); - else if ( aMap.Contains( anId ) ) - aResult++; - } - } + // for each pair of nodes in anEdge (there are 2 pairs in a quadratic edge) + // count elements containing both nodes of the pair. + // Note that there may be such cases for a quadratic edge (a horizontal line): + // + // Case 1 Case 2 + // | | | | | + // | | | | | + // +-----+------+ +-----+------+ + // | | | | + // | | | | + // result sould 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 ); + if ( aNode1 == aLastNode ) aNode1 = 0; + + SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator(); + while( anElemIter->more() ) { + const SMDS_MeshElement* anElem = anElemIter->next(); + if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) { + SMDS_ElemIteratorPtr anIter = anElem->nodesIterator(); + while ( anIter->more() ) { + if ( const SMDS_MeshElement* anElemNode = anIter->next() ) { + if ( anElemNode == aNode0 ) { + aResult0++; + if ( !aNode1 ) break; // not a quadratic edge + } + else if ( anElemNode == aNode1 ) + aResult1++; + } + } } } + int aResult = std::max ( aResult0, aResult1 ); 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 <= numeric_limits::min()); + if ( !zeroLen ) + n /= len; + + if (ok) *ok = !zeroLen; + + return n; + } } @@ -124,13 +189,16 @@ namespace{ using namespace SMESH::Controls; /* - FUNCTORS -*/ + * FUNCTORS + */ +//================================================================================ /* Class : NumericalFunctor Description : Base class for numerical functors */ +//================================================================================ + NumericalFunctor::NumericalFunctor(): myMesh(NULL) { @@ -142,7 +210,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(); @@ -150,26 +218,52 @@ 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, - TSequenceOfXYZ& theRes ) + TSequenceOfXYZ& theRes ) { 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 = anElem->nodesIterator(); - if ( anIter != 0 ) - { - while( anIter->more() ) - { - const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); - if ( aNode != 0 ){ - theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); + 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(); + } + + if ( anIter ) { + double xyz[3]; + while( anIter->more() ) { + if ( const SMDS_MeshNode* aNode = static_cast( anIter->next() )) + { + aNode->GetXYZ( xyz ); + theRes.push_back( gp_XYZ( xyz[0], xyz[1], xyz[2] )); } } } @@ -185,29 +279,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 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() ) { - double aVal = GetValue( P ); - if ( myPrecision >= 0 ) + SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator( GetType() ); + while ( elemIt->more() ) + values.insert( GetValue( elemIt->next()->GetID() )); + } + else + { + 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 ) + { + // 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 ) { @@ -219,31 +409,287 @@ 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_ElementType aType = aElem->GetType(); + int len = P.size(); + switch( aType ) { + case SMDSAbs_Volume: + if( len == 4 ) { // tetras + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + double L4 = getDistance(P( 1 ),P( 4 )); + double L5 = getDistance(P( 2 ),P( 4 )); + double L6 = getDistance(P( 3 ),P( 4 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + break; + } + else if( len == 5 ) { // pyramids + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + double L5 = getDistance(P( 1 ),P( 5 )); + double L6 = getDistance(P( 2 ),P( 5 )); + double L7 = getDistance(P( 3 ),P( 5 )); + double L8 = getDistance(P( 4 ),P( 5 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(L7,L8)); + break; + } + else if( len == 6 ) { // pentas + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 4 )); + double L7 = getDistance(P( 1 ),P( 4 )); + double L8 = getDistance(P( 2 ),P( 5 )); + double L9 = getDistance(P( 3 ),P( 6 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),L9)); + break; + } + else if( len == 8 ) { // hexas + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 7 )); + double L7 = getDistance(P( 7 ),P( 8 )); + double L8 = getDistance(P( 8 ),P( 5 )); + double L9 = getDistance(P( 1 ),P( 5 )); + double L10= getDistance(P( 2 ),P( 6 )); + double L11= getDistance(P( 3 ),P( 7 )); + double L12= getDistance(P( 4 ),P( 8 )); + double D1 = getDistance(P( 1 ),P( 7 )); + double D2 = getDistance(P( 2 ),P( 8 )); + double D3 = getDistance(P( 3 ),P( 5 )); + double D4 = getDistance(P( 4 ),P( 6 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10))); + aVal = Max(aVal,Max(L11,L12)); + aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); + break; + } + else if( len == 12 ) { // hexagonal prism + for ( int i1 = 1; i1 < 12; ++i1 ) + for ( int i2 = i1+1; i1 <= 12; ++i1 ) + aVal = Max( aVal, getDistance(P( i1 ),P( i2 ))); + break; + } + else if( len == 10 ) { // quadratic tetras + double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 )); + double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 )); + double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + break; + } + else if( len == 13 ) { // quadratic pyramids + double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 )); + double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 )); + double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(L7,L8)); + break; + } + else if( len == 15 ) { // quadratic pentas + double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 )); + double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 )); + double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 )); + double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),L9)); + break; + } + else if( len == 20 || len == 27 ) { // quadratic hexas + double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 )); + double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 )); + double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 )); + double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 )); + double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 )); + double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 )); + double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 )); + double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 )); + double D1 = getDistance(P( 1 ),P( 7 )); + double D2 = getDistance(P( 2 ),P( 8 )); + double D3 = getDistance(P( 3 ),P( 5 )); + double D4 = getDistance(P( 4 ),P( 6 )); + aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); + aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10))); + aVal = Max(aVal,Max(L11,L12)); + aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4))); + break; + } + else if( len > 1 && aElem->IsPoly() ) { // polys + // get the maximum distance between all pairs of nodes + for( int i = 1; i <= len; i++ ) { + for( int j = 1; j <= len; j++ ) { + if( j > i ) { // optimization of the loop + double D = getDistance( P(i), P(j) ); + aVal = Max( aVal, D ); + } + } + } + } + } + + 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 ) { @@ -255,12 +701,13 @@ double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) aMin = getAngle(P( P.size() ), P( 1 ), P( 2 )); aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 ))); - for (int i=2; iFindElement( 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 // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis. // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4 + // PAL10872 int nbNodes = P.size(); if ( nbNodes < 3 ) return 0; - // Compute lengths of the sides - - 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 ) ); - // Compute aspect ratio - if ( nbNodes == 3 ) - { + 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 ) ); // 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.; - + if ( anArea <= theEps ) + return theInf; return alfa * maxLen * half_perimeter / anArea; } - else - { - // return aspect ratio of the worst triange which can be built + 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 <= theEps ) + return theInf; + return alfa * maxLen * half_perimeter / anArea; + } + else if( nbNodes == 4 ) { // quadrangle + // Compute lengths of the sides + std::vector< double > aLen (4); + aLen[0] = getDistance( P(1), P(2) ); + aLen[1] = getDistance( P(2), P(3) ); + aLen[2] = getDistance( P(3), P(4) ); + aLen[3] = getDistance( P(4), P(1) ); + // Compute lengths of the diagonals + std::vector< double > aDia (2); + aDia[0] = getDistance( P(1), P(3) ); + aDia[1] = getDistance( P(2), P(4) ); + // Compute areas of all triangles which can be built + // taking three nodes of the quadrangle + std::vector< double > anArea (4); + 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 + std::vector< double > aLen (4); + aLen[0] = getDistance( P(1), P(3) ); + aLen[1] = getDistance( P(3), P(5) ); + aLen[2] = getDistance( P(5), P(7) ); + aLen[3] = getDistance( P(7), P(1) ); + // Compute lengths of the diagonals + std::vector< double > aDia (2); + 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(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; + std::vector< 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 @@ -357,10 +910,13 @@ SMDSAbs_ElementType AspectRatio::GetType() const } +//================================================================================ /* Class : AspectRatio3D Description : Functor for calculating aspect ratio */ +//================================================================================ + namespace{ inline double getHalfPerimeter(double theTria[3]){ @@ -369,9 +925,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]){ @@ -413,21 +969,55 @@ namespace{ inline double getMaxHeight(double theLen[6]) { - double aHeight = max(theLen[0],theLen[1]); - aHeight = max(aHeight,theLen[2]); - aHeight = max(aHeight,theLen[3]); - aHeight = max(aHeight,theLen[4]); - aHeight = max(aHeight,theLen[5]); + double aHeight = std::max(theLen[0],theLen[1]); + aHeight = std::max(aHeight,theLen[2]); + aHeight = std::max(aHeight,theLen[3]); + aHeight = std::max(aHeight,theLen[4]); + aHeight = std::max(aHeight,theLen[5]); return aHeight; } } +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; + if(myCurrElement->IsPoly()) return aQuality; + int nbNodes = P.size(); - switch(nbNodes){ + + 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 return aQuality; + } + + switch(nbNodes) { case 4:{ double aLen[6] = { getDistance(P( 1 ),P( 2 )), // a @@ -460,191 +1050,207 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) //double aVolume = getVolume(aLen); double aHeight = getMaxHeight(aLen); static double aCoeff = sqrt(2.0)/12.0; - aQuality = aCoeff*aHeight*aSumArea/aVolume; + if ( aVolume > DBL_MIN ) + aQuality = aCoeff*aHeight*aSumArea/aVolume; break; } case 5:{ { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } break; } case 6:{ { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } break; } case 8:{ { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } { gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )}; - aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); } 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; @@ -657,7 +1263,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true ); for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face points( p + 1 ) = P( pInd[ p ] + 1 ); - aQuality = max( aQuality, aspect2D.GetValue( points )); + aQuality = std::max( aQuality, aspect2D.GetValue( points )); } } return aQuality; @@ -677,23 +1283,30 @@ SMDSAbs_ElementType AspectRatio3D::GetType() const } +//================================================================================ /* Class : Warping Description : Functor for calculating warping */ +//================================================================================ + double Warping::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 4 ) return 0; - gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4; + 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 ); - 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, @@ -704,20 +1317,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 @@ -734,37 +1347,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; + 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; } @@ -774,49 +1394,55 @@ 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; - gp_XYZ p23 = ( p3 + p2 ) / 2; - gp_XYZ p31 = ( p3 + p1 ) / 2; + gp_XYZ p12 = ( p2 + p1 ) / 2.; + gp_XYZ p23 = ( p3 + p2 ) / 2.; + gp_XYZ p31 = ( p3 + p1 ) / 2.; gp_Vec v1( p31 - p2 ), v2( p12 - p23 ); - return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 ); + return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 ); } double Skew::GetValue( const TSequenceOfXYZ& P ) { if ( P.size() != 3 && P.size() != 4 ) - return 0; + 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 { - gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2; - gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2; - gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2; - gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2; + gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.; + gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.; + gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.; + gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.; gp_Vec v1( p34 - p12 ), v2( p23 - p41 ); double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution() - ? 0 : fabs( PI2 - v1.Angle( v2 ) ); + ? 0. : fabs( PI2 - v1.Angle( v2 ) ); + + double val = A * 180. / M_PI; - return A * 180 / PI; + const double eps = 0.1; // val is in degrees + + return val < eps ? 0. : val; } } @@ -834,28 +1460,37 @@ SMDSAbs_ElementType Skew::GetType() const } +//================================================================================ /* Class : Area Description : Functor for calculating area */ +//================================================================================ + double Area::GetValue( const TSequenceOfXYZ& P ) { - double aArea = 0; - if ( P.size() == 3 ) - return getArea( P( 1 ), P( 2 ), P( 3 ) ); - else if (P.size() > 3) - aArea = getArea( P( 1 ), P( 2 ), P( 3 ) ); - else - return 0; - - for (int i=4; i<=P.size(); i++) - aArea += getArea(P(1),P(i-1),P(i)); - return aArea; -} + 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 (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); + } + val = SumVec.Magnitude() * 0.5; + } + return val; +} double Area::GetBadRate( double Value, int /*nbNodes*/ ) const { - // meaningless as it is not quality control functor + // meaningless as it is not a quality control functor return Value; } @@ -864,14 +1499,20 @@ 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 ) { - return ( P.size() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0 ); + switch ( P.size() ) { + case 2: return getDistance( P( 1 ), P( 2 ) ); + case 3: return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) ); + default: return 0.; + } } double Length::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -885,112 +1526,240 @@ 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( long theElementId ) { TSequenceOfXYZ P; - if (GetPoints(theElementId,P)){ - - double aVal;// = GetValue( P ); - const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId ); - SMDSAbs_ElementType aType = aElem->GetType(); - + if ( GetPoints( theElementId, P )) + { + double aVal = 0; int len = P.size(); + SMDSAbs_EntityType aType = P.getElementEntity(); - switch (aType){ - case SMDSAbs_All: - case SMDSAbs_Node: - case SMDSAbs_Edge: - if (len == 2){ - aVal = getDistance( P( 1 ), P( 2 ) ); - break; - } - case SMDSAbs_Face: + switch (aType) { + case SMDSEntity_Edge: + if (len == 2) + aVal = getDistance( P( 1 ), P( 2 ) ); + break; + case SMDSEntity_Quad_Edge: + if (len == 3) // quadratic edge + aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 )); + break; + case SMDSEntity_Triangle: if (len == 3){ // triangles - double L1 = getDistance(P( 1 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 1 )); - aVal = Max(L1,Max(L2,L3)); - break; + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + aVal = Min(L1,Min(L2,L3)); } - 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; + break; + case SMDSEntity_Quadrangle: + if (len == 4){ // quadrangles + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + aVal = Min(Min(L1,L2),Min(L3,L4)); } - case SMDSAbs_Volume: - if (len == 4){ // tetraidrs - double L1 = getDistance(P( 1 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 1 )); - double L4 = getDistance(P( 1 ),P( 4 )); - double L5 = getDistance(P( 2 ),P( 4 )); - double L6 = getDistance(P( 3 ),P( 4 )); - aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); - break; + break; + case SMDSEntity_Quad_Triangle: + case SMDSEntity_BiQuad_Triangle: + if (len >= 6){ // quadratic triangles + double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 )); + double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 )); + double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 )); + aVal = Min(L1,Min(L2,L3)); } - else if (len == 5){ // piramids - double L1 = getDistance(P( 1 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 1 )); - double L4 = getDistance(P( 4 ),P( 1 )); - double L5 = getDistance(P( 1 ),P( 5 )); - double L6 = getDistance(P( 2 ),P( 5 )); - double L7 = getDistance(P( 3 ),P( 5 )); - double L8 = getDistance(P( 4 ),P( 5 )); - - aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); - aVal = Max(aVal,Max(L7,L8)); - break; + break; + case SMDSEntity_Quad_Quadrangle: + case SMDSEntity_BiQuad_Quadrangle: + if (len >= 8){ // quadratic quadrangles + double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 )); + double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 )); + double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 )); + double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 )); + aVal = Min(Min(L1,L2),Min(L3,L4)); } - else if (len == 6){ // pentaidres - double L1 = getDistance(P( 1 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 1 )); - double L4 = getDistance(P( 4 ),P( 5 )); - double L5 = getDistance(P( 5 ),P( 6 )); - double L6 = getDistance(P( 6 ),P( 4 )); - double L7 = getDistance(P( 1 ),P( 4 )); - double L8 = getDistance(P( 2 ),P( 5 )); - double L9 = getDistance(P( 3 ),P( 6 )); - - aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); - aVal = Max(aVal,Max(Max(L7,L8),L9)); - break; + break; + case SMDSEntity_Tetra: + if (len == 4){ // tetrahedra + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + double L4 = getDistance(P( 1 ),P( 4 )); + double L5 = getDistance(P( 2 ),P( 4 )); + double L6 = getDistance(P( 3 ),P( 4 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); } - else if (len == 8){ // hexaider - double L1 = getDistance(P( 1 ),P( 2 )); - double L2 = getDistance(P( 2 ),P( 3 )); - double L3 = getDistance(P( 3 ),P( 4 )); - double L4 = getDistance(P( 4 ),P( 1 )); - double L5 = getDistance(P( 5 ),P( 6 )); - double L6 = getDistance(P( 6 ),P( 7 )); - double L7 = getDistance(P( 7 ),P( 8 )); - double L8 = getDistance(P( 8 ),P( 5 )); - double L9 = getDistance(P( 1 ),P( 5 )); - double L10= getDistance(P( 2 ),P( 6 )); - double L11= getDistance(P( 3 ),P( 7 )); - double L12= getDistance(P( 4 ),P( 8 )); - - aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6)); - aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10))); - aVal = Max(aVal,Max(L11,L12)); - break; - + break; + case SMDSEntity_Pyramid: + if (len == 5){ // piramids + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + double L5 = getDistance(P( 1 ),P( 5 )); + double L6 = getDistance(P( 2 ),P( 5 )); + double L7 = getDistance(P( 3 ),P( 5 )); + double L8 = getDistance(P( 4 ),P( 5 )); + + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(L7,L8)); } - - default: aVal=-1; + break; + case SMDSEntity_Penta: + if (len == 6) { // pentaidres + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 4 )); + double L7 = getDistance(P( 1 ),P( 4 )); + double L8 = getDistance(P( 2 ),P( 5 )); + double L9 = getDistance(P( 3 ),P( 6 )); + + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(Min(L7,L8),L9)); + } + break; + case SMDSEntity_Hexa: + if (len == 8){ // hexahedron + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 7 )); + double L7 = getDistance(P( 7 ),P( 8 )); + double L8 = getDistance(P( 8 ),P( 5 )); + double L9 = getDistance(P( 1 ),P( 5 )); + double L10= getDistance(P( 2 ),P( 6 )); + double L11= getDistance(P( 3 ),P( 7 )); + double L12= getDistance(P( 4 ),P( 8 )); + + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10))); + aVal = Min(aVal,Min(L11,L12)); + } + break; + case SMDSEntity_Quad_Tetra: + if (len == 10){ // quadratic tetraidrs + double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 )); + double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 )); + double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + } + break; + case SMDSEntity_Quad_Pyramid: + if (len == 13){ // quadratic piramids + double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 )); + double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 )); + double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(L7,L8)); + } + break; + case SMDSEntity_Quad_Penta: + if (len == 15){ // quadratic pentaidres + double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 )); + double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 )); + double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 )); + double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 )); + double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(Min(L7,L8),L9)); + } + break; + case SMDSEntity_Quad_Hexa: + case SMDSEntity_TriQuad_Hexa: + if (len >= 20) { // quadratic hexaider + double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 )); + double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 )); + double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 )); + double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 )); + double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 )); + double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 )); + double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 )); + double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10))); + aVal = Min(aVal,Min(L11,L12)); + } + break; + case SMDSEntity_Polygon: + if ( len > 1 ) { + aVal = getDistance( P(1), P( P.size() )); + for ( size_t i = 1; i < P.size(); ++i ) + aVal = Min( aVal, getDistance( P( i ), P( i+1 ))); + } + break; + case SMDSEntity_Quad_Polygon: + if ( len > 2 ) { + aVal = getDistance( P(1), P( P.size() )) + getDistance( P(P.size()), P( P.size()-1 )); + for ( size_t i = 1; i < P.size()-1; i += 2 ) + aVal = Min( aVal, getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 ))); + } + break; + case SMDSEntity_Hexagonal_Prism: + if (len == 12) { // hexagonal prism + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 5 )); + double L5 = getDistance(P( 5 ),P( 6 )); + double L6 = getDistance(P( 6 ),P( 1 )); + + double L7 = getDistance(P( 7 ), P( 8 )); + double L8 = getDistance(P( 8 ), P( 9 )); + double L9 = getDistance(P( 9 ), P( 10 )); + double L10= getDistance(P( 10 ),P( 11 )); + double L11= getDistance(P( 11 ),P( 12 )); + double L12= getDistance(P( 12 ),P( 7 )); + + double L13 = getDistance(P( 1 ),P( 7 )); + double L14 = getDistance(P( 2 ),P( 8 )); + double L15 = getDistance(P( 3 ),P( 9 )); + double L16 = getDistance(P( 4 ),P( 10 )); + double L17 = getDistance(P( 5 ),P( 11 )); + double L18 = getDistance(P( 6 ),P( 12 )); + aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6)); + aVal = Min(aVal, Min(Min(Min(L7,L8),Min(L9,L10)),Min(L11,L12))); + aVal = Min(aVal, Min(Min(Min(L13,L14),Min(L15,L16)),Min(L17,L18))); + } + break; + case SMDSEntity_Polyhedra: + { + } + break; + default: + return 0; } - if (aVal <0){ + if (aVal < 0 ) { return 0.; } @@ -1008,7 +1777,7 @@ double Length2D::GetValue( long theElementId) double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const { - // meaningless as it is not quality control functor + // meaningless as it is not a quality control functor return Value; } @@ -1026,57 +1795,105 @@ Length2D::Value::Value(double theLength,long thePntId1, long thePntId2): } } -bool Length2D::Value::operator<(const Length2D::Value& x) const{ +bool Length2D::Value::operator<(const Length2D::Value& x) const +{ if(myPntId[0] < x.myPntId[0]) return true; if(myPntId[0] == x.myPntId[0]) if(myPntId[1] < x.myPntId[1]) return true; return false; } -void Length2D::GetValues(TValues& theValues){ +void Length2D::GetValues(TValues& theValues) +{ TValues aValues; SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); for(; anIter->more(); ){ const SMDS_MeshFace* anElem = anIter->next(); - SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator(); - long aNodeId[2]; - gp_Pnt P[3]; - double aLength; - const SMDS_MeshElement* aNode; - if(aNodesIter->more()){ - aNode = aNodesIter->next(); - const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode; - P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z()); - aNodeId[0] = aNodeId[1] = aNode->GetID(); - aLength = 0; + if(anElem->IsQuadratic()) { + const SMDS_VtkFace* F = + dynamic_cast(anElem); + // use special nodes iterator + SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator(); + long aNodeId[4]; + gp_Pnt P[4]; + + double aLength; + const SMDS_MeshElement* aNode; + if(anIter->more()){ + aNode = anIter->next(); + const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode; + P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z()); + aNodeId[0] = aNodeId[1] = aNode->GetID(); + aLength = 0; + } + for(; anIter->more(); ){ + const SMDS_MeshNode* N1 = static_cast (anIter->next()); + P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z()); + aNodeId[2] = N1->GetID(); + aLength = P[1].Distance(P[2]); + if(!anIter->more()) break; + const SMDS_MeshNode* N2 = static_cast (anIter->next()); + P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z()); + aNodeId[3] = N2->GetID(); + aLength += P[2].Distance(P[3]); + 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[0]); + theValues.insert(aValue1); + theValues.insert(aValue2); } - 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()); + else { + SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator(); + long aNodeId[2]; + gp_Pnt P[3]; + + double aLength; + const SMDS_MeshElement* aNode; + if(aNodesIter->more()){ + aNode = aNodesIter->next(); + const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode; + P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z()); + aNodeId[0] = aNodeId[1] = aNode->GetID(); + aLength = 0; + } + for(; aNodesIter->more(); ){ + aNode = aNodesIter->next(); + const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode; + long anId = aNode->GetID(); + + P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z()); + + aLength = P[1].Distance(P[2]); + + Value aValue(aLength,aNodeId[1],anId); + aNodeId[1] = anId; + P[1] = P[2]; + theValues.insert(aValue); + } - aLength = P[1].Distance(P[2]); + aLength = P[0].Distance(P[1]); - Value aValue(aLength,aNodeId[1],anId); - aNodeId[1] = anId; - P[1] = P[2]; + Value aValue(aLength,aNodeId[0],aNodeId[1]); theValues.insert(aValue); } - - aLength = P[0].Distance(P[1]); - - Value aValue(aLength,aNodeId[0],aNodeId[1]); - theValues.insert(aValue); } } +//================================================================================ /* Class : MultiConnection Description : Functor for calculating number of faces conneted to the edge */ +//================================================================================ + double MultiConnection::GetValue( const TSequenceOfXYZ& P ) { return 0; @@ -1097,10 +1914,13 @@ SMDSAbs_ElementType MultiConnection::GetType() const return SMDSAbs_Edge; } +//================================================================================ /* Class : MultiConnection2D Description : Functor for calculating number of faces conneted to the edge */ +//================================================================================ + double MultiConnection2D::GetValue( const TSequenceOfXYZ& P ) { return 0; @@ -1108,60 +1928,58 @@ double MultiConnection2D::GetValue( const TSequenceOfXYZ& P ) double MultiConnection2D::GetValue( long theElementId ) { - TSequenceOfXYZ P; int aResult = 0; - if (GetPoints(theElementId,P)){ - const SMDS_MeshElement* anFaceElem = myMesh->FindElement( theElementId ); - SMDSAbs_ElementType aType = anFaceElem->GetType(); - - int len = P.size(); - - TColStd_MapOfInteger aMap; - int aResult = 0; + const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId); + SMDSAbs_ElementType aType = aFaceElem->GetType(); - switch (aType){ - case SMDSAbs_All: - case SMDSAbs_Node: - case SMDSAbs_Edge: - case SMDSAbs_Face: - if (len == 3){ // triangles - int Nb[3] = {0,0,0}; - - int i=0; - SMDS_ElemIteratorPtr anIter = anFaceElem->nodesIterator(); - if ( anIter != 0 ) { - while( anIter->more() ) { - const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); - if ( aNode == 0 ){ - break; - } - SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); - while( anElemIter->more() ) { - const SMDS_MeshElement* anElem = anElemIter->next(); - if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) { - int anId = anElem->GetID(); - - if ( anIter->more() ) // i.e. first node - aMap.Add( anId ); - else if ( aMap.Contains( anId ) ){ - Nb[i]++; - } - } - else if ( anElem != 0 && anElem->GetType() == SMDSAbs_Edge ) i++; - } - } - } - - aResult = Max(Max(Nb[0],Nb[1]),Nb[2]); + switch (aType) { + case SMDSAbs_Face: + { + int i = 0, len = aFaceElem->NbNodes(); + SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator(); + if (!anIter) break; + + const SMDS_MeshNode *aNode, *aNode0; + TColStd_MapOfInteger aMap, aMapPrev; + + for (i = 0; i <= len; i++) { + aMapPrev = aMap; + aMap.Clear(); + + int aNb = 0; + if (anIter->more()) { + aNode = (SMDS_MeshNode*)anIter->next(); + } else { + if (i == len) + aNode = aNode0; + else + break; + } + if (!aNode) break; + if (i == 0) aNode0 = aNode; + + SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); + while (anElemIter->more()) { + const SMDS_MeshElement* anElem = anElemIter->next(); + if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) { + int anId = anElem->GetID(); + + aMap.Add(anId); + if (aMapPrev.Contains(anId)) { + aNb++; + } + } + } + aResult = Max(aResult, aNb); } - break; - case SMDSAbs_Volume: - default: aResult=0; } - + break; + default: + aResult = 0; } - return aResult;//getNbMultiConnection( myMesh, theId ); + + return aResult; } double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -1183,18 +2001,26 @@ MultiConnection2D::Value::Value(long thePntId1, long thePntId2) } } -bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{ +bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const +{ if(myPntId[0] < x.myPntId[0]) return true; if(myPntId[0] == x.myPntId[0]) if(myPntId[1] < x.myPntId[1]) return true; return false; } -void MultiConnection2D::GetValues(MValues& theValues){ +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 = anElem->nodesIterator(); + SMDS_ElemIteratorPtr aNodesIter; + if ( anElem->IsQuadratic() ) + aNodesIter = dynamic_cast + (anElem)->interlacedNodesElemIterator(); + else + aNodesIter = anElem->nodesIterator(); long aNodeId[3]; //int aNbConnects=0; @@ -1207,7 +2033,7 @@ void MultiConnection2D::GetValues(MValues& theValues){ const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1; aNodeId[0] = aNodeId[1] = aNodes->GetID(); } - for(; aNodesIter->more(); ){ + for(; aNodesIter->more(); ) { aNode2 = (SMDS_MeshNode*) aNodesIter->next(); long anId = aNode2->GetID(); aNodeId[2] = anId; @@ -1215,11 +2041,12 @@ void MultiConnection2D::GetValues(MValues& theValues){ 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; + aItr->second += 1; + //aNbConnects = nb; + } + else { + theValues[aValue] = 1; + //aNbConnects = 1; } //cout << "NodeIds: "<FindNode( theNodeId ); + if (!aNode) + return false; + + return (aNode->NbInverseElements() < 1); } -//======================================================================= -// name : GetRangeStr -// Purpose : Get range as a string. -// Example: "1,2,3,50-60,63,67,70-" -//======================================================================= -void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr ) +SMDSAbs_ElementType FreeNodes::GetType() const { - theResStr.Clear(); + return SMDSAbs_Node; +} - TColStd_SequenceOfInteger anIntSeq; - TColStd_SequenceOfAsciiString aStrSeq; - TColStd_MapIteratorOfMapOfInteger anIter( myIds ); - for ( ; anIter.More(); anIter.Next() ) - { - int anId = anIter.Key(); - TCollection_AsciiString aStr( anId ); - anIntSeq.Append( anId ); - aStrSeq.Append( aStr ); - } - - for ( int i = 1, n = myMin.Length(); i <= n; i++ ) - { - int aMinId = myMin( i ); - int aMaxId = myMax( i ); +//================================================================================ +/* + Class : FreeFaces + Description : Predicate for free faces +*/ +//================================================================================ - TCollection_AsciiString aStr; - if ( aMinId != IntegerFirst() ) - aStr += aMinId; +FreeFaces::FreeFaces() +{ + myMesh = 0; +} - aStr += "-"; +void FreeFaces::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; +} - if ( aMaxId != IntegerLast() ) - aStr += aMaxId; +bool FreeFaces::IsSatisfy( long theId ) +{ + if (!myMesh) return false; + // check that faces nodes refers to less than two common volumes + const SMDS_MeshElement* aFace = myMesh->FindElement( theId ); + if ( !aFace || aFace->GetType() != SMDSAbs_Face ) + return false; - // find position of the string in result sequence and insert string in it - if ( anIntSeq.Length() == 0 ) - { - anIntSeq.Append( aMinId ); - aStrSeq.Append( aStr ); - } - else - { - if ( aMinId < anIntSeq.First() ) - { - anIntSeq.Prepend( aMinId ); - aStrSeq.Prepend( aStr ); - } - else if ( aMinId > anIntSeq.Last() ) - { - anIntSeq.Append( aMinId ); - aStrSeq.Append( aStr ); - } - else - for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ ) - if ( aMinId < anIntSeq( j ) ) - { - anIntSeq.InsertBefore( j, aMinId ); - aStrSeq.InsertBefore( j, aStr ); - break; - } - } + int nbNode = aFace->NbNodes(); + + // collect volumes to check that number of volumes with count equal nbNode not less than 2 + typedef map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters + typedef map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator + TMapOfVolume mapOfVol; + + SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator(); + while ( nodeItr->more() ) { + const SMDS_MeshNode* aNode = static_cast(nodeItr->next()); + if ( !aNode ) continue; + SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume); + while ( volItr->more() ) { + SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next(); + TItrMapOfVolume itr = mapOfVol.insert(make_pair(aVol, 0)).first; + (*itr).second++; + } } + int nbVol = 0; + TItrMapOfVolume volItr = mapOfVol.begin(); + TItrMapOfVolume volEnd = mapOfVol.end(); + for ( ; volItr != volEnd; ++volItr ) + if ( (*volItr).second >= nbNode ) + nbVol++; + // face is not free if number of volumes constructed on thier nodes more than one + return (nbVol < 2); +} - if ( aStrSeq.Length() == 0 ) - return; - - theResStr = aStrSeq( 1 ); - for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ ) - { - theResStr += ","; - theResStr += aStrSeq( j ); - } +SMDSAbs_ElementType FreeFaces::GetType() const +{ + return SMDSAbs_Face; } -//======================================================================= -// name : SetRangeStr -// Purpose : Define range with string -// Example of entry string: "1,2,3,50-60,63,67,70-" -//======================================================================= -bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr ) +//================================================================================ +/* + Class : LinearOrQuadratic + Description : Predicate to verify whether a mesh element is linear +*/ +//================================================================================ + +LinearOrQuadratic::LinearOrQuadratic() { - myMin.Clear(); - myMax.Clear(); - myIds.Clear(); + myMesh = 0; +} - TCollection_AsciiString aStr = theStr; - aStr.RemoveAll( ' ' ); - aStr.RemoveAll( '\t' ); +void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; +} - for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) ) - aStr.Remove( aPos, 2 ); +bool LinearOrQuadratic::IsSatisfy( long theId ) +{ + if (!myMesh) return false; + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) ) + return false; + return (!anElem->IsQuadratic()); +} - TCollection_AsciiString tmpStr = aStr.Token( ",", 1 ); - int i = 1; - while ( tmpStr != "" ) - { - tmpStr = aStr.Token( ",", i++ ); - int aPos = tmpStr.Search( '-' ); +void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType ) +{ + myType = theType; +} - if ( aPos == -1 ) - { - if ( tmpStr.IsIntegerValue() ) - myIds.Add( tmpStr.IntegerValue() ); - else - return false; - } - else - { - TCollection_AsciiString aMaxStr = tmpStr.Split( aPos ); - TCollection_AsciiString aMinStr = tmpStr; +SMDSAbs_ElementType LinearOrQuadratic::GetType() const +{ + return myType; +} - while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' ); - while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' ); +//================================================================================ +/* + Class : GroupColor + Description : Functor for check color of group to whic mesh element belongs to +*/ +//================================================================================ - if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() || - !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() ) - return false; +GroupColor::GroupColor() +{ +} - myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() ); - myMax.Append( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() ); - } - } +bool GroupColor::IsSatisfy( long theId ) +{ + return myIDs.count( theId ); +} - return true; +void GroupColor::SetType( SMDSAbs_ElementType theType ) +{ + myType = theType; } -//======================================================================= -// name : GetType -// Purpose : Get type of supported entities -//======================================================================= -SMDSAbs_ElementType RangeOfIds::GetType() const +SMDSAbs_ElementType GroupColor::GetType() const { return myType; } -//======================================================================= -// name : SetType -// Purpose : Set type of supported entities -//======================================================================= -void RangeOfIds::SetType( SMDSAbs_ElementType theType ) +static bool isEqual( const Quantity_Color& theColor1, + const Quantity_Color& theColor2 ) { - myType = theType; + // tolerance to compare colors + const double tol = 5*1e-3; + return ( fabs( theColor1.Red() - theColor2.Red() ) < tol && + fabs( theColor1.Green() - theColor2.Green() ) < tol && + fabs( theColor1.Blue() - theColor2.Blue() ) < tol ); } -//======================================================================= -// name : IsSatisfy -// Purpose : Verify whether entity satisfies to this rpedicate -//======================================================================= -bool RangeOfIds::IsSatisfy( long theId ) +void GroupColor::SetMesh( const SMDS_Mesh* theMesh ) { - if ( !myMesh ) - return false; + myIDs.clear(); - if ( myType == SMDSAbs_Node ) - { - if ( myMesh->FindNode( theId ) == 0 ) - return false; - } - else + const SMESHDS_Mesh* aMesh = dynamic_cast(theMesh); + if ( !aMesh ) + return; + + int nbGrp = aMesh->GetNbGroups(); + if ( !nbGrp ) + return; + + // iterates on groups and find necessary elements ids + const std::set& aGroups = aMesh->GetGroups(); + set::const_iterator GrIt = aGroups.begin(); + for (; GrIt != aGroups.end(); GrIt++) { - const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); - if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All ) - return false; + SMESHDS_GroupBase* aGrp = (*GrIt); + if ( !aGrp ) + continue; + // check type and color of group + if ( !isEqual( myColor, aGrp->GetColor() )) + continue; + + // IPAL52867 (prevent infinite recursion via GroupOnFilter) + if ( SMESHDS_GroupOnFilter * gof = dynamic_cast< SMESHDS_GroupOnFilter* >( aGrp )) + if ( gof->GetPredicate().get() == this ) + continue; + + SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType(); + if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) { + // add elements IDS into control + int aSize = aGrp->Extent(); + for (int i = 0; i < aSize; i++) + myIDs.insert( aGrp->GetID(i+1) ); + } } +} - if ( myIds.Contains( theId ) ) - return true; +void GroupColor::SetColorStr( const TCollection_AsciiString& theStr ) +{ + Kernel_Utils::Localizer loc; + TCollection_AsciiString aStr = theStr; + aStr.RemoveAll( ' ' ); + aStr.RemoveAll( '\t' ); + for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) ) + aStr.Remove( aPos, 2 ); + Standard_Real clr[3]; + clr[0] = clr[1] = clr[2] = 0.; + for ( int i = 0; i < 3; i++ ) { + TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 ); + if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() ) + clr[i] = tmpStr.RealValue(); + } + myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB ); +} - for ( int i = 1, n = myMin.Length(); i <= n; i++ ) - if ( theId >= myMin( i ) && theId <= myMax( i ) ) - return true; +//======================================================================= +// name : GetRangeStr +// Purpose : Get range as a string. +// Example: "1,2,3,50-60,63,67,70-" +//======================================================================= - return false; +void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const +{ + theResStr.Clear(); + theResStr += TCollection_AsciiString( myColor.Red() ); + theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() ); + theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() ); } +//================================================================================ /* - Class : Comparator - Description : Base class for comparators + Class : ElemGeomType + Description : Predicate to check element geometry type */ -Comparator::Comparator(): - myMargin(0) -{} - -Comparator::~Comparator() -{} +//================================================================================ -void Comparator::SetMesh( const SMDS_Mesh* theMesh ) +ElemGeomType::ElemGeomType() { - if ( myFunctor ) - myFunctor->SetMesh( theMesh ); + myMesh = 0; + myType = SMDSAbs_All; + myGeomType = SMDSGeom_TRIANGLE; } -void Comparator::SetMargin( double theValue ) +void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh ) { - myMargin = theValue; + myMesh = theMesh; } -void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct ) +bool ElemGeomType::IsSatisfy( long theId ) { - myFunctor = theFunct; + if (!myMesh) return false; + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + if ( !anElem ) + return false; + const SMDSAbs_ElementType anElemType = anElem->GetType(); + if ( myType != SMDSAbs_All && anElemType != myType ) + return false; + bool isOk = ( anElem->GetGeomType() == myGeomType ); + return isOk; } -SMDSAbs_ElementType Comparator::GetType() const +void ElemGeomType::SetType( SMDSAbs_ElementType theType ) { - return myFunctor ? myFunctor->GetType() : SMDSAbs_All; + myType = theType; } -double Comparator::GetMargin() +SMDSAbs_ElementType ElemGeomType::GetType() const { - return myMargin; + return myType; } - -/* - Class : LessThan - Description : Comparator "<" -*/ -bool LessThan::IsSatisfy( long theId ) +void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType ) { - return myFunctor && myFunctor->GetValue( theId ) < myMargin; + myGeomType = theType; } - -/* - Class : MoreThan - Description : Comparator ">" -*/ -bool MoreThan::IsSatisfy( long theId ) +SMDSAbs_GeometryType ElemGeomType::GetGeomType() const { - return myFunctor && myFunctor->GetValue( theId ) > myMargin; + return myGeomType; } - +//================================================================================ /* - Class : EqualTo - Description : Comparator "=" + Class : ElemEntityType + Description : Predicate to check element entity type */ -EqualTo::EqualTo(): - myToler(Precision::Confusion()) -{} +//================================================================================ -bool EqualTo::IsSatisfy( long theId ) +ElemEntityType::ElemEntityType(): + myMesh( 0 ), + myType( SMDSAbs_All ), + myEntityType( SMDSEntity_0D ) { - return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler; } -void EqualTo::SetTolerance( double theToler ) +void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh ) { - myToler = theToler; + myMesh = theMesh; } -double EqualTo::GetTolerance() +bool ElemEntityType::IsSatisfy( long theId ) { - return myToler; + if ( !myMesh ) return false; + if ( myType == SMDSAbs_Node ) + return myMesh->FindNode( theId ); + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + return ( anElem && + myEntityType == anElem->GetEntityType() ); } -/* - Class : LogicalNOT - Description : Logical NOT predicate -*/ -LogicalNOT::LogicalNOT() -{} - -LogicalNOT::~LogicalNOT() -{} - -bool LogicalNOT::IsSatisfy( long theId ) +void ElemEntityType::SetType( SMDSAbs_ElementType theType ) { - return myPredicate && !myPredicate->IsSatisfy( theId ); + myType = theType; } -void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh ) +SMDSAbs_ElementType ElemEntityType::GetType() const { - if ( myPredicate ) - myPredicate->SetMesh( theMesh ); + return myType; } -void LogicalNOT::SetPredicate( PredicatePtr thePred ) +void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType ) { - myPredicate = thePred; + myEntityType = theEntityType; } -SMDSAbs_ElementType LogicalNOT::GetType() const +SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const { - return myPredicate ? myPredicate->GetType() : SMDSAbs_All; + return myEntityType; } +//================================================================================ +/*! + * \brief Class ConnectedElements + */ +//================================================================================ -/* - Class : LogicalBinary - Description : Base class for binary logical predicate -*/ -LogicalBinary::LogicalBinary() -{} +ConnectedElements::ConnectedElements(): + myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {} -LogicalBinary::~LogicalBinary() -{} +SMDSAbs_ElementType ConnectedElements::GetType() const +{ return myType; } -void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh ) -{ - if ( myPredicate1 ) - myPredicate1->SetMesh( theMesh ); +int ConnectedElements::GetNode() const +{ return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ - if ( myPredicate2 ) - myPredicate2->SetMesh( theMesh ); -} +std::vector ConnectedElements::GetPoint() const +{ return myXYZ; } -void LogicalBinary::SetPredicate1( PredicatePtr thePredicate ) -{ - myPredicate1 = thePredicate; -} +void ConnectedElements::clearOkIDs() +{ myOkIDsReady = false; myOkIDs.clear(); } -void LogicalBinary::SetPredicate2( PredicatePtr thePredicate ) +void ConnectedElements::SetType( SMDSAbs_ElementType theType ) { - myPredicate2 = thePredicate; + if ( myType != theType || myMeshModifTracer.IsMeshModified() ) + clearOkIDs(); + myType = theType; } -SMDSAbs_ElementType LogicalBinary::GetType() const +void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh ) { - if ( !myPredicate1 || !myPredicate2 ) - return SMDSAbs_All; - - SMDSAbs_ElementType aType1 = myPredicate1->GetType(); - SMDSAbs_ElementType aType2 = myPredicate2->GetType(); - - return aType1 == aType2 ? aType1 : SMDSAbs_All; + 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 + } } - -/* - Class : LogicalAND - Description : Logical AND -*/ -bool LogicalAND::IsSatisfy( long theId ) +void ConnectedElements::SetNode( int nodeID ) { - return - myPredicate1 && - myPredicate2 && - myPredicate1->IsSatisfy( theId ) && - myPredicate2->IsSatisfy( theId ); -} + myNodeID = nodeID; + myXYZ.clear(); - -/* - Class : LogicalOR - Description : Logical OR -*/ -bool LogicalOR::IsSatisfy( long theId ) -{ - return - myPredicate1 && - myPredicate2 && - myPredicate1->IsSatisfy( theId ) || - myPredicate2->IsSatisfy( theId ); + 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; -/* - FILTER -*/ - -Filter::Filter() -{} + bool isSameDomain = false; -Filter::~Filter() -{} + // find myNodeID by myXYZ if possible + if ( myMeshModifTracer.GetMesh() ) + { + auto_ptr searcher + ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() )); -void Filter::SetPredicate( PredicatePtr thePredicate ) -{ - myPredicate = thePredicate; -} + vector< const SMDS_MeshElement* > foundElems; + searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems ); -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 ); + if ( !foundElems.empty() ) + { + myNodeID = foundElems[0]->GetNode(0)->GetID(); + if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() ) + isSameDomain = IsSatisfy( foundElems[0]->GetID() ); } } + if ( !isSameDomain ) + clearOkIDs(); } -void -Filter:: -GetElementsId( const SMDS_Mesh* theMesh, - PredicatePtr thePredicate, - TIdSequence& theSequence ) +bool ConnectedElements::IsSatisfy( long theElementId ) { - theSequence.clear(); + // Here we do NOT check if the mesh has changed, we do it in Set...() only!!! - if ( !theMesh || !thePredicate ) - return; + if ( !myOkIDsReady ) + { + if ( !myMeshModifTracer.GetMesh() ) + return false; + const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID ); + if ( !node0 ) + return false; - thePredicate->SetMesh( theMesh ); + 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(); - 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; + // 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 ); } -void -Filter::GetElementsId( const SMDS_Mesh* theMesh, - Filter::TIdSequence& theSequence ) +//================================================================================ +/*! + * \brief Class CoplanarFaces + */ +//================================================================================ + +CoplanarFaces::CoplanarFaces() + : myFaceID(0), myToler(0) { - GetElementsId(theMesh,myPredicate,theSequence); } +void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified() ) + { + // Build a set of coplanar face ids -/* - ManifoldPart -*/ + myCoplanarIDs.clear(); -typedef std::set TMapOfFacePtr; + if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler ) + return; -/* - Internal class Link -*/ + const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID ); + if ( !face || face->GetType() != SMDSAbs_Face ) + return; -ManifoldPart::Link::Link( SMDS_MeshNode* theNode1, - SMDS_MeshNode* theNode2 ) -{ - myNode1 = theNode1; - myNode2 = theNode2; -} + bool normOK; + gp_Vec myNorm = getNormale( static_cast(face), &normOK ); + if (!normOK) + return; -ManifoldPart::Link::~Link() -{ - myNode1 = 0; - myNode2 = 0; -} + const double radianTol = myToler * M_PI / 180.; + std::set< SMESH_TLink > checkedLinks; -bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const -{ - if ( myNode1 == theLink.myNode1 && - myNode2 == theLink.myNode2 ) - return true; - else if ( myNode1 == theLink.myNode2 && - myNode2 == theLink.myNode1 ) - return true; - else - return false; -} + std::list< pair< const SMDS_MeshElement*, gp_Vec > > faceQueue; + faceQueue.push_back( make_pair( face, myNorm )); + while ( !faceQueue.empty() ) + { + face = faceQueue.front().first; + myNorm = faceQueue.front().second; + faceQueue.pop_front(); -bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const -{ - if(myNode1 < x.myNode1) return true; - if(myNode1 == x.myNode1) - if(myNode2 < x.myNode2) return true; - return false; + 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.insert( SMESH_TLink( n1, n2 )).second ) + 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 || myNorm.Angle( norm ) <= radianTol) + { + myCoplanarIDs.insert( f->GetID() ); + faceQueue.push_back( make_pair( f, norm )); + } + } + } + } + } + } } - -bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1, - const ManifoldPart::Link& theLink2 ) +bool CoplanarFaces::IsSatisfy( long theElementId ) { - return theLink1.IsEqual( theLink2 ); + return myCoplanarIDs.count( theElementId ); } -ManifoldPart::ManifoldPart() -{ - myMesh = 0; - myAngToler = Precision::Angular(); - myIsOnlyManifold = true; -} +/* + *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-" +*/ -ManifoldPart::~ManifoldPart() +//======================================================================= +// name : RangeOfIds +// Purpose : Constructor +//======================================================================= +RangeOfIds::RangeOfIds() { myMesh = 0; + myType = SMDSAbs_All; } -void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh ) +//======================================================================= +// name : SetMesh +// Purpose : Set mesh +//======================================================================= +void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh ) { myMesh = theMesh; - process(); } -SMDSAbs_ElementType ManifoldPart::GetType() const -{ return SMDSAbs_Face; } - -bool ManifoldPart::IsSatisfy( long theElementId ) +//======================================================================= +// name : AddToRange +// Purpose : Add ID to the range +//======================================================================= +bool RangeOfIds::AddToRange( long theEntityId ) { - return myMapIds.Contains( theElementId ); + myIds.Add( theEntityId ); + return true; } -void ManifoldPart::SetAngleTolerance( const double theAngToler ) -{ myAngToler = theAngToler; } - -double ManifoldPart::GetAngleTolerance() const -{ return myAngToler; } - -void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly ) -{ myIsOnlyManifold = theIsOnly; } - -void ManifoldPart::SetStartElem( const long theStartId ) -{ myStartElemId = theStartId; } - -bool ManifoldPart::process() +//======================================================================= +// name : GetRangeStr +// Purpose : Get range as a string. +// Example: "1,2,3,50-60,63,67,70-" +//======================================================================= +void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr ) { - myMapIds.Clear(); - myMapBadGeomIds.Clear(); + theResStr.Clear(); - myAllFacePtr.clear(); - myAllFacePtrIntDMap.clear(); - if ( !myMesh ) - return false; + TColStd_SequenceOfInteger anIntSeq; + TColStd_SequenceOfAsciiString aStrSeq; - // collect all faces into own map - SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator(); - for (; anFaceItr->more(); ) + TColStd_MapIteratorOfMapOfInteger anIter( myIds ); + for ( ; anIter.More(); anIter.Next() ) { - SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next(); - myAllFacePtr.push_back( aFacePtr ); - myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1; + int anId = anIter.Key(); + TCollection_AsciiString aStr( anId ); + anIntSeq.Append( anId ); + aStrSeq.Append( aStr ); } - SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId ); - if ( !aStartFace ) - return false; + for ( int i = 1, n = myMin.Length(); i <= n; i++ ) + { + int aMinId = myMin( i ); + int aMaxId = myMax( i ); - // the map of non manifold links and bad geometry - TMapOfLink aMapOfNonManifold; - TColStd_MapOfInteger aMapOfTreated; + TCollection_AsciiString aStr; + if ( aMinId != IntegerFirst() ) + aStr += aMinId; - // begin cycle on faces from start index and run on vector till the end - // and from begin to start index to cover whole vector - const int aStartIndx = myAllFacePtrIntDMap[aStartFace]; - bool isStartTreat = false; - for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ ) - { - if ( fi == aStartIndx ) - isStartTreat = true; - // as result next time when fi will be equal to aStartIndx + aStr += "-"; - SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ]; + if ( aMaxId != IntegerLast() ) + 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 ); + } + else + { + if ( aMinId < anIntSeq.First() ) + { + anIntSeq.Prepend( aMinId ); + aStrSeq.Prepend( aStr ); + } + else if ( aMinId > anIntSeq.Last() ) + { + anIntSeq.Append( aMinId ); + aStrSeq.Append( aStr ); + } + else + for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ ) + if ( aMinId < anIntSeq( j ) ) + { + anIntSeq.InsertBefore( j, aMinId ); + aStrSeq.InsertBefore( j, aStr ); + break; + } + } + } + + if ( aStrSeq.Length() == 0 ) + return; + + theResStr = aStrSeq( 1 ); + for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ ) + { + theResStr += ","; + theResStr += aStrSeq( j ); + } +} + +//======================================================================= +// name : SetRangeStr +// Purpose : Define range with string +// Example of entry string: "1,2,3,50-60,63,67,70-" +//======================================================================= +bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr ) +{ + myMin.Clear(); + myMax.Clear(); + myIds.Clear(); + + TCollection_AsciiString aStr = theStr; + //aStr.RemoveAll( ' ' ); + //aStr.RemoveAll( '\t' ); + for ( int i = 1; i <= aStr.Length(); ++i ) + if ( isspace( aStr.Value( i ))) + aStr.SetValue( i, ','); + + for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) ) + aStr.Remove( aPos, 1 ); + + TCollection_AsciiString tmpStr = aStr.Token( ",", 1 ); + int i = 1; + while ( tmpStr != "" ) + { + tmpStr = aStr.Token( ",", i++ ); + int aPos = tmpStr.Search( '-' ); + + if ( aPos == -1 ) + { + if ( tmpStr.IsIntegerValue() ) + myIds.Add( tmpStr.IntegerValue() ); + else + return false; + } + else + { + TCollection_AsciiString aMaxStr = tmpStr.Split( aPos ); + TCollection_AsciiString aMinStr = tmpStr; + + while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' ); + while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' ); + + if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) || + (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) ) + return false; + + myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() ); + myMax.Append( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() ); + } + } + + return true; +} + +//======================================================================= +// name : GetType +// Purpose : Get type of supported entities +//======================================================================= +SMDSAbs_ElementType RangeOfIds::GetType() const +{ + return myType; +} + +//======================================================================= +// name : SetType +// Purpose : Set type of supported entities +//======================================================================= +void RangeOfIds::SetType( SMDSAbs_ElementType theType ) +{ + myType = theType; +} + +//======================================================================= +// name : IsSatisfy +// Purpose : Verify whether entity satisfies to this rpedicate +//======================================================================= +bool RangeOfIds::IsSatisfy( long theId ) +{ + if ( !myMesh ) + return false; + + if ( myType == SMDSAbs_Node ) + { + if ( myMesh->FindNode( theId ) == 0 ) + return false; + } + else + { + const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); + if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All )) + return false; + } + + if ( myIds.Contains( theId ) ) + return true; + + for ( int i = 1, n = myMin.Length(); i <= n; i++ ) + if ( theId >= myMin( i ) && theId <= myMax( i ) ) + return true; + + return false; +} + +/* + Class : Comparator + Description : Base class for comparators +*/ +Comparator::Comparator(): + myMargin(0) +{} + +Comparator::~Comparator() +{} + +void Comparator::SetMesh( const SMDS_Mesh* theMesh ) +{ + if ( myFunctor ) + myFunctor->SetMesh( theMesh ); +} + +void Comparator::SetMargin( double theValue ) +{ + myMargin = theValue; +} + +void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct ) +{ + myFunctor = theFunct; +} + +SMDSAbs_ElementType Comparator::GetType() const +{ + return myFunctor ? myFunctor->GetType() : SMDSAbs_All; +} + +double Comparator::GetMargin() +{ + return myMargin; +} + + +/* + Class : LessThan + Description : Comparator "<" +*/ +bool LessThan::IsSatisfy( long theId ) +{ + return myFunctor && myFunctor->GetValue( theId ) < myMargin; +} + + +/* + Class : MoreThan + Description : Comparator ">" +*/ +bool MoreThan::IsSatisfy( long theId ) +{ + return myFunctor && myFunctor->GetValue( theId ) > myMargin; +} + + +/* + Class : EqualTo + Description : Comparator "=" +*/ +EqualTo::EqualTo(): + myToler(Precision::Confusion()) +{} + +bool EqualTo::IsSatisfy( long theId ) +{ + return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler; +} + +void EqualTo::SetTolerance( double theToler ) +{ + myToler = theToler; +} + +double EqualTo::GetTolerance() +{ + return myToler; +} + +/* + Class : LogicalNOT + Description : Logical NOT predicate +*/ +LogicalNOT::LogicalNOT() +{} + +LogicalNOT::~LogicalNOT() +{} + +bool LogicalNOT::IsSatisfy( long theId ) +{ + return myPredicate && !myPredicate->IsSatisfy( theId ); +} + +void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh ) +{ + if ( myPredicate ) + myPredicate->SetMesh( theMesh ); +} + +void LogicalNOT::SetPredicate( PredicatePtr thePred ) +{ + myPredicate = thePred; +} + +SMDSAbs_ElementType LogicalNOT::GetType() const +{ + return myPredicate ? myPredicate->GetType() : SMDSAbs_All; +} + + +/* + Class : LogicalBinary + Description : Base class for binary logical predicate +*/ +LogicalBinary::LogicalBinary() +{} + +LogicalBinary::~LogicalBinary() +{} + +void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh ) +{ + if ( myPredicate1 ) + myPredicate1->SetMesh( theMesh ); + + if ( myPredicate2 ) + myPredicate2->SetMesh( theMesh ); +} + +void LogicalBinary::SetPredicate1( PredicatePtr thePredicate ) +{ + myPredicate1 = thePredicate; +} + +void LogicalBinary::SetPredicate2( PredicatePtr thePredicate ) +{ + myPredicate2 = thePredicate; +} + +SMDSAbs_ElementType LogicalBinary::GetType() const +{ + if ( !myPredicate1 || !myPredicate2 ) + return SMDSAbs_All; + + SMDSAbs_ElementType aType1 = myPredicate1->GetType(); + SMDSAbs_ElementType aType2 = myPredicate2->GetType(); + + return aType1 == aType2 ? aType1 : SMDSAbs_All; +} + + +/* + Class : LogicalAND + Description : Logical AND +*/ +bool LogicalAND::IsSatisfy( long theId ) +{ + return + myPredicate1 && + myPredicate2 && + myPredicate1->IsSatisfy( theId ) && + myPredicate2->IsSatisfy( theId ); +} + + +/* + Class : LogicalOR + Description : Logical OR +*/ +bool LogicalOR::IsSatisfy( long theId ) +{ + return + myPredicate1 && + myPredicate2 && + (myPredicate1->IsSatisfy( theId ) || + myPredicate2->IsSatisfy( 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() +{} + +Filter::~Filter() +{} + +void Filter::SetPredicate( PredicatePtr thePredicate ) +{ + myPredicate = thePredicate; +} + +void Filter::GetElementsId( const SMDS_Mesh* theMesh, + PredicatePtr thePredicate, + TIdSequence& theSequence ) +{ + theSequence.clear(); + + if ( !theMesh || !thePredicate ) + return; + + 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 ); + } + } +} + +void Filter::GetElementsId( const SMDS_Mesh* theMesh, + Filter::TIdSequence& theSequence ) +{ + GetElementsId(theMesh,myPredicate,theSequence); +} + +/* + ManifoldPart +*/ + +typedef std::set TMapOfFacePtr; + +/* + Internal class Link +*/ + +ManifoldPart::Link::Link( SMDS_MeshNode* theNode1, + SMDS_MeshNode* theNode2 ) +{ + myNode1 = theNode1; + myNode2 = theNode2; +} + +ManifoldPart::Link::~Link() +{ + myNode1 = 0; + myNode2 = 0; +} + +bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const +{ + if ( myNode1 == theLink.myNode1 && + myNode2 == theLink.myNode2 ) + return true; + else if ( myNode1 == theLink.myNode2 && + myNode2 == theLink.myNode1 ) + return true; + else + return false; +} + +bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const +{ + if(myNode1 < x.myNode1) return true; + if(myNode1 == x.myNode1) + if(myNode2 < x.myNode2) return true; + return false; +} + +bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1, + const ManifoldPart::Link& theLink2 ) +{ + return theLink1.IsEqual( theLink2 ); +} + +ManifoldPart::ManifoldPart() +{ + myMesh = 0; + myAngToler = Precision::Angular(); + myIsOnlyManifold = true; +} + +ManifoldPart::~ManifoldPart() +{ + myMesh = 0; +} + +void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; + process(); +} + +SMDSAbs_ElementType ManifoldPart::GetType() const +{ return SMDSAbs_Face; } + +bool ManifoldPart::IsSatisfy( long theElementId ) +{ + return myMapIds.Contains( theElementId ); +} + +void ManifoldPart::SetAngleTolerance( const double theAngToler ) +{ myAngToler = theAngToler; } + +double ManifoldPart::GetAngleTolerance() const +{ return myAngToler; } + +void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly ) +{ myIsOnlyManifold = theIsOnly; } + +void ManifoldPart::SetStartElem( const long theStartId ) +{ myStartElemId = theStartId; } + +bool ManifoldPart::process() +{ + myMapIds.Clear(); + myMapBadGeomIds.Clear(); + + myAllFacePtr.clear(); + myAllFacePtrIntDMap.clear(); + if ( !myMesh ) + return false; + + // collect all faces into own map + SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator(); + for (; anFaceItr->more(); ) + { + SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next(); + myAllFacePtr.push_back( aFacePtr ); + myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1; + } + + SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId ); + if ( !aStartFace ) + return false; + + // the map of non manifold links and bad geometry + TMapOfLink aMapOfNonManifold; + TColStd_MapOfInteger aMapOfTreated; + + // begin cycle on faces from start index and run on vector till the end + // and from begin to start index to cover whole vector + const int aStartIndx = myAllFacePtrIntDMap[aStartFace]; + bool isStartTreat = false; + for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ ) + { + if ( fi == aStartIndx ) + isStartTreat = true; + // as result next time when fi will be equal to aStartIndx + + SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ]; if ( aMapOfTreated.Contains( aFacePtr->GetID() ) ) continue; - aMapOfTreated.Add( aFacePtr->GetID() ); - TColStd_MapOfInteger aResFaces; - if ( !findConnected( myAllFacePtrIntDMap, aFacePtr, - aMapOfNonManifold, aResFaces ) ) - continue; - TColStd_MapIteratorOfMapOfInteger anItr( aResFaces ); - for ( ; anItr.More(); anItr.Next() ) + aMapOfTreated.Add( aFacePtr->GetID() ); + TColStd_MapOfInteger aResFaces; + if ( !findConnected( myAllFacePtrIntDMap, aFacePtr, + aMapOfNonManifold, aResFaces ) ) + continue; + TColStd_MapIteratorOfMapOfInteger anItr( aResFaces ); + for ( ; anItr.More(); anItr.Next() ) + { + int aFaceId = anItr.Key(); + aMapOfTreated.Add( aFaceId ); + myMapIds.Add( aFaceId ); + } + + if ( fi == ( myAllFacePtr.size() - 1 ) ) + fi = 0; + } // end run on vector of faces + return !myMapIds.IsEmpty(); +} + +static void getLinks( const SMDS_MeshFace* theFace, + ManifoldPart::TVectorOfLink& theLinks ) +{ + int aNbNode = theFace->NbNodes(); + SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator(); + int i = 1; + SMDS_MeshNode* aNode = 0; + for ( ; aNodeItr->more() && i <= aNbNode; ) + { + + SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next(); + if ( i == 1 ) + aNode = aN1; + i++; + SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next(); + i++; + ManifoldPart::Link aLink( aN1, aN2 ); + theLinks.push_back( aLink ); + } +} + +bool ManifoldPart::findConnected + ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt, + SMDS_MeshFace* theStartFace, + ManifoldPart::TMapOfLink& theNonManifold, + TColStd_MapOfInteger& theResFaces ) +{ + theResFaces.Clear(); + if ( !theAllFacePtrInt.size() ) + return false; + + if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() ) + { + myMapBadGeomIds.Add( theStartFace->GetID() ); + return false; + } + + ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip; + ManifoldPart::TVectorOfLink aSeqOfBoundary; + theResFaces.Add( theStartFace->GetID() ); + ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace; + + expandBoundary( aMapOfBoundary, aSeqOfBoundary, + aDMapLinkFace, theNonManifold, theStartFace ); + + bool isDone = false; + while ( !isDone && aMapOfBoundary.size() != 0 ) + { + bool isToReset = false; + ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin(); + for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink ) + { + ManifoldPart::Link aLink = *pLink; + if ( aMapToSkip.find( aLink ) != aMapToSkip.end() ) + continue; + // each link could be treated only once + aMapToSkip.insert( aLink ); + + ManifoldPart::TVectorOfFacePtr aFaces; + // find next + if ( myIsOnlyManifold && + (theNonManifold.find( aLink ) != theNonManifold.end()) ) + continue; + else + { + getFacesByLink( aLink, aFaces ); + // filter the element to keep only indicated elements + ManifoldPart::TVectorOfFacePtr aFiltered; + ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin(); + for ( ; pFace != aFaces.end(); ++pFace ) + { + SMDS_MeshFace* aFace = *pFace; + if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() ) + aFiltered.push_back( aFace ); + } + aFaces = aFiltered; + if ( aFaces.size() < 2 ) // no neihgbour faces + continue; + else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case + { + theNonManifold.insert( aLink ); + continue; + } + } + + // compare normal with normals of neighbor element + SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ]; + ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin(); + for ( ; pFace != aFaces.end(); ++pFace ) + { + SMDS_MeshFace* aNextFace = *pFace; + if ( aPrevFace == aNextFace ) + continue; + int anNextFaceID = aNextFace->GetID(); + if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) ) + // should not be with non manifold restriction. probably bad topology + continue; + // check if face was treated and skipped + if ( myMapBadGeomIds.Contains( anNextFaceID ) || + !isInPlane( aPrevFace, aNextFace ) ) + continue; + // add new element to connected and extend the boundaries. + theResFaces.Add( anNextFaceID ); + expandBoundary( aMapOfBoundary, aSeqOfBoundary, + aDMapLinkFace, theNonManifold, aNextFace ); + isToReset = true; + } + } + isDone = !isToReset; + } + + return !theResFaces.IsEmpty(); +} + +bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1, + const SMDS_MeshFace* theFace2 ) +{ + gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) ); + gp_XYZ aNorm2XYZ = getNormale( theFace2 ); + if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() ) + { + myMapBadGeomIds.Add( theFace2->GetID() ); + return false; + } + if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) ) + return true; + + return false; +} + +void ManifoldPart::expandBoundary + ( ManifoldPart::TMapOfLink& theMapOfBoundary, + ManifoldPart::TVectorOfLink& theSeqOfBoundary, + ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr, + ManifoldPart::TMapOfLink& theNonManifold, + SMDS_MeshFace* theNextFace ) const +{ + ManifoldPart::TVectorOfLink aLinks; + getLinks( theNextFace, aLinks ); + int aNbLink = (int)aLinks.size(); + for ( int i = 0; i < aNbLink; i++ ) + { + ManifoldPart::Link aLink = aLinks[ i ]; + if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) ) + continue; + if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() ) + { + if ( myIsOnlyManifold ) + { + // remove from boundary + theMapOfBoundary.erase( aLink ); + ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin(); + for ( ; pLink != theSeqOfBoundary.end(); ++pLink ) + { + ManifoldPart::Link aBoundLink = *pLink; + if ( aBoundLink.IsEqual( aLink ) ) + { + theSeqOfBoundary.erase( pLink ); + break; + } + } + } + } + else + { + theMapOfBoundary.insert( aLink ); + theSeqOfBoundary.push_back( aLink ); + theDMapLinkFacePtr[ aLink ] = theNextFace; + } + } +} + +void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink, + ManifoldPart::TVectorOfFacePtr& theFaces ) const +{ + 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 ); + } + // take all faces that shared second node + anItr = theLink.myNode2->facesIterator(); + // find the common part of two sets + for ( ; anItr->more(); ) + { + SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next(); + if ( aSetOfFaces.count( aFace ) ) + theFaces.push_back( aFace ); + } +} + +/* + 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() +{ + myIds.Clear(); + myType = SMDSAbs_All; + mySurf.Nullify(); + myToler = Precision::Confusion(); + myUseBoundaries = false; +} + +ElementsOnSurface::~ElementsOnSurface() +{ +} + +void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified()) + process(); +} + +bool ElementsOnSurface::IsSatisfy( long theElementId ) +{ + return myIds.Contains( theElementId ); +} + +SMDSAbs_ElementType ElementsOnSurface::GetType() const +{ return myType; } + +void ElementsOnSurface::SetTolerance( const double theToler ) +{ + if ( myToler != theToler ) + myIds.Clear(); + myToler = theToler; +} + +double ElementsOnSurface::GetTolerance() const +{ return myToler; } + +void ElementsOnSurface::SetUseBoundaries( bool theUse ) +{ + if ( myUseBoundaries != theUse ) { + myUseBoundaries = theUse; + SetSurface( mySurf, myType ); + } +} + +void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape, + const SMDSAbs_ElementType theType ) +{ + myIds.Clear(); + myType = theType; + mySurf.Nullify(); + if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE ) + return; + mySurf = TopoDS::Face( theShape ); + BRepAdaptor_Surface SA( mySurf, myUseBoundaries ); + Standard_Real + u1 = SA.FirstUParameter(), + u2 = SA.LastUParameter(), + v1 = SA.FirstVParameter(), + v2 = SA.LastVParameter(); + Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf ); + myProjector.Init( surf, u1,u2, v1,v2 ); + process(); +} + +void ElementsOnSurface::process() +{ + myIds.Clear(); + if ( mySurf.IsNull() ) + return; + + if ( !myMeshModifTracer.GetMesh() ) + return; + + myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType )); + + SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType ); + for(; anIter->more(); ) + process( anIter->next() ); +} + +void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr ) +{ + SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator(); + bool isSatisfy = true; + for ( ; aNodeItr->more(); ) + { + SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); + if ( !isOnSurface( aNode ) ) { - int aFaceId = anItr.Key(); - aMapOfTreated.Add( aFaceId ); - myMapIds.Add( aFaceId ); + isSatisfy = false; + break; } + } + if ( isSatisfy ) + myIds.Add( theElemPtr->GetID() ); +} - if ( fi == ( myAllFacePtr.size() - 1 ) ) - fi = 0; - } // end run on vector of faces - return !myMapIds.IsEmpty(); +bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) +{ + if ( mySurf.IsNull() ) + return false; + + gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() ); + // double aToler2 = myToler * myToler; +// if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane))) +// { +// gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln(); +// if ( aPln.SquareDistance( aPnt ) > aToler2 ) +// return false; +// } +// else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface))) +// { +// gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder(); +// double aRad = aCyl.Radius(); +// gp_Ax3 anAxis = aCyl.Position(); +// gp_XYZ aLoc = aCyl.Location().XYZ(); +// double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc ); +// double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc ); +// if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 ) +// return false; +// } +// else +// return false; + myProjector.Perform( aPnt ); + bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler ); + + return isOn; } -static void getLinks( const SMDS_MeshFace* theFace, - ManifoldPart::TVectorOfLink& theLinks ) + +/* + ElementsOnShape +*/ + +ElementsOnShape::ElementsOnShape() + : //myMesh(0), + myType(SMDSAbs_All), + myToler(Precision::Confusion()), + myAllNodesFlag(false) { - int aNbNode = theFace->NbNodes(); - SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator(); - int i = 1; - SMDS_MeshNode* aNode = 0; - for ( ; aNodeItr->more() && i <= aNbNode; ) +} + +ElementsOnShape::~ElementsOnShape() +{ + clearClassifiers(); +} + +SMDSAbs_ElementType ElementsOnShape::GetType() const +{ + return myType; +} + +void ElementsOnShape::SetTolerance (const double theToler) +{ + if (myToler != theToler) { + myToler = theToler; + SetShape(myShape, myType); + } +} + +double ElementsOnShape::GetTolerance() const +{ + return myToler; +} + +void ElementsOnShape::SetAllNodes (bool theAllNodes) +{ + 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 ); + } + } +} - SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next(); - if ( i == 1 ) - aNode = aN1; - i++; - SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next(); - i++; - ManifoldPart::Link aLink( aN1, aN2 ); - theLinks.push_back( aLink ); +bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut ) +{ + if ( n->GetID() >= (int) myNodeIsChecked.size() || + !myNodeIsChecked[ n->GetID() ]) + return false; + + isOut = myNodeIsOut[ n->GetID() ]; + return true; +} + +void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool isOut ) +{ + if ( n->GetID() < (int) myNodeIsChecked.size() ) + { + myNodeIsChecked[ n->GetID() ] = true; + myNodeIsOut [ n->GetID() ] = isOut; } } -static gp_XYZ getNormale( const SMDS_MeshFace* theFace ) +void ElementsOnShape::SetShape (const TopoDS_Shape& theShape, + const SMDSAbs_ElementType theType) { - 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++ ) + myType = theType; + myShape = theShape; + if ( myShape.IsNull() ) return; + + TopTools_IndexedMapOfShape shapesMap; + TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }; + TopExp_Explorer sub; + for ( int i = 0; i < 4; ++i ) { - SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); - anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); + 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() ); } - gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1); - gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1); - n = q1 ^ q2; - if ( aNbNode > 3 ) + clearClassifiers(); + myClassifiers.resize( shapesMap.Extent() ); + for ( int i = 0; i < shapesMap.Extent(); ++i ) + myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler ); + + if ( theType == SMDSAbs_Node ) + { + SMESHUtils::FreeVector( myNodeIsChecked ); + SMESHUtils::FreeVector( myNodeIsOut ); + } + else { - gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1); - n += q2 ^ q3; + std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false ); } - double len = n.Modulus(); - if ( len > 0 ) - n /= len; +} - return n; +void ElementsOnShape::clearClassifiers() +{ + for ( size_t i = 0; i < myClassifiers.size(); ++i ) + delete myClassifiers[ i ]; + myClassifiers.clear(); } -bool ManifoldPart::findConnected - ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt, - SMDS_MeshFace* theStartFace, - ManifoldPart::TMapOfLink& theNonManifold, - TColStd_MapOfInteger& theResFaces ) +bool ElementsOnShape::IsSatisfy (long elemId) { - theResFaces.Clear(); - if ( !theAllFacePtrInt.size() ) + const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh(); + const SMDS_MeshElement* elem = + ( myType == SMDSAbs_Node ? mesh->FindNode( elemId ) : mesh->FindElement( elemId )); + if ( !elem || myClassifiers.empty() ) return false; - if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() ) + bool isSatisfy = myAllNodesFlag, isNodeOut; + + gp_XYZ centerXYZ (0, 0, 0); + + SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator(); + while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) { - myMapBadGeomIds.Add( theStartFace->GetID() ); - return false; + SMESH_TNodeXYZ aPnt( aNodeItr->next() ); + centerXYZ += aPnt; + + isNodeOut = true; + if ( !getNodeIsOut( aPnt._node, isNodeOut )) + { + for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i ) + isNodeOut = myClassifiers[i]->IsOut( aPnt ); + + setNodeIsOut( aPnt._node, isNodeOut ); + } + isSatisfy = !isNodeOut; } - ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip; - ManifoldPart::TVectorOfLink aSeqOfBoundary; - theResFaces.Add( theStartFace->GetID() ); - ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace; + // Check the center point for volumes MantisBug 0020168 + if (isSatisfy && + myAllNodesFlag && + myClassifiers[0]->ShapeType() == TopAbs_SOLID) + { + centerXYZ /= elem->NbNodes(); + isSatisfy = false; + for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i ) + isSatisfy = ! myClassifiers[i]->IsOut( centerXYZ ); + } - expandBoundary( aMapOfBoundary, aSeqOfBoundary, - aDMapLinkFace, theNonManifold, theStartFace ); + return isSatisfy; +} - bool isDone = false; - while ( !isDone && aMapOfBoundary.size() != 0 ) +TopAbs_ShapeEnum ElementsOnShape::TClassifier::ShapeType() const +{ + return myShape.ShapeType(); +} + +bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p) +{ + return (this->*myIsOutFun)( p ); +} + +void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double theTol) +{ + myShape = theShape; + myTol = theTol; + switch ( myShape.ShapeType() ) { - bool isToReset = false; - ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin(); - for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink ) + case TopAbs_SOLID: { + if ( isBox( theShape )) { - ManifoldPart::Link aLink = *pLink; - if ( aMapToSkip.find( aLink ) != aMapToSkip.end() ) - continue; - // each link could be treated only once - aMapToSkip.insert( aLink ); + myIsOutFun = & ElementsOnShape::TClassifier::isOutOfBox; + } + else + { + mySolidClfr.Load(theShape); + myIsOutFun = & ElementsOnShape::TClassifier::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::TClassifier::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::TClassifier::isOutOfEdge; + break; + } + case TopAbs_VERTEX:{ + myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) ); + myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex; + break; + } + default: + throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier"); + } +} - ManifoldPart::TVectorOfFacePtr aFaces; - // find next - if ( myIsOnlyManifold && - (theNonManifold.find( aLink ) != theNonManifold.end()) ) - continue; - else - { - getFacesByLink( aLink, aFaces ); - // filter the element to keep only indicated elements - ManifoldPart::TVectorOfFacePtr aFiltered; - ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin(); - for ( ; pFace != aFaces.end(); ++pFace ) - { - SMDS_MeshFace* aFace = *pFace; - if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() ) - aFiltered.push_back( aFace ); - } - aFaces = aFiltered; - if ( aFaces.size() < 2 ) // no neihgbour faces - continue; - else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case - { - theNonManifold.insert( aLink ); - continue; - } - } +bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p) +{ + mySolidClfr.Perform( p, myTol ); + return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON ); +} - // compare normal with normals of neighbor element - SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ]; - ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin(); - for ( ; pFace != aFaces.end(); ++pFace ) - { - SMDS_MeshFace* aNextFace = *pFace; - if ( aPrevFace == aNextFace ) - continue; - int anNextFaceID = aNextFace->GetID(); - if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) ) - // should not be with non manifold restriction. probably bad topology - continue; - // check if face was treated and skipped - if ( myMapBadGeomIds.Contains( anNextFaceID ) || - !isInPlane( aPrevFace, aNextFace ) ) - continue; - // add new element to connected and extend the boundaries. - theResFaces.Add( anNextFaceID ); - expandBoundary( aMapOfBoundary, aSeqOfBoundary, - aDMapLinkFace, theNonManifold, aNextFace ); - isToReset = true; +bool ElementsOnShape::TClassifier::isOutOfBox (const gp_Pnt& p) +{ + return myBox.IsOut( p.XYZ() ); +} + +bool ElementsOnShape::TClassifier::isOutOfFace (const gp_Pnt& p) +{ + myProjFace.Perform( p ); + if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol ) + { + // check relatively to the face + Quantity_Parameter 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; + } + return true; +} + +bool ElementsOnShape::TClassifier::isOutOfEdge (const gp_Pnt& p) +{ + myProjEdge.Perform( p ); + return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol ); +} + +bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p) +{ + return ( myVertexXYZ.Distance( p ) > myTol ); +} + +bool ElementsOnShape::TClassifier::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 ) + { + gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))); + for ( int iC = 1; iC <= 3; ++ iC ) + { + 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; +} + + +/* + Class : BelongToGeom + Description : Predicate for verifying whether entity belongs to + specified geometrical support +*/ + +BelongToGeom::BelongToGeom() + : myMeshDS(NULL), + myType(SMDSAbs_All), + myIsSubshape(false), + myTolerance(Precision::Confusion()) +{} + +void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMeshDS = dynamic_cast(theMesh); + init(); +} + +void BelongToGeom::SetGeom( const TopoDS_Shape& 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; } } - isDone = !isToReset; + return true; } - return !theResFaces.IsEmpty(); + return false; } -bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1, - const SMDS_MeshFace* theFace2 ) +void BelongToGeom::init() { - gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) ); - gp_XYZ aNorm2XYZ = getNormale( theFace2 ); - if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() ) + 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) // to be always ready to check an element not bound to geometry { - myMapBadGeomIds.Add( theFace2->GetID() ); - return false; + myElementsOnShapePtr.reset(new ElementsOnShape()); + myElementsOnShapePtr->SetTolerance(myTolerance); + myElementsOnShapePtr->SetAllNodes(true); // "belong", while false means "lays on" + myElementsOnShapePtr->SetMesh(myMeshDS); + myElementsOnShapePtr->SetShape(myShape, myType); } - if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) ) - return true; +} +static bool IsContains( const SMESHDS_Mesh* theMeshDS, + const TopoDS_Shape& theShape, + const SMDS_MeshElement* theElem, + TopAbs_ShapeEnum theFindShapeEnum, + TopAbs_ShapeEnum theAvoidShapeEnum = TopAbs_SHAPE ) +{ + TopExp_Explorer anExp( theShape,theFindShapeEnum,theAvoidShapeEnum ); + + while( anExp.More() ) + { + const TopoDS_Shape& aShape = anExp.Current(); + if( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ) ){ + if( aSubMesh->Contains( theElem ) ) + return true; + } + anExp.Next(); + } return false; } -void ManifoldPart::expandBoundary - ( ManifoldPart::TMapOfLink& theMapOfBoundary, - ManifoldPart::TVectorOfLink& theSeqOfBoundary, - ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr, - ManifoldPart::TMapOfLink& theNonManifold, - SMDS_MeshFace* theNextFace ) const +bool BelongToGeom::IsSatisfy (long theId) { - ManifoldPart::TVectorOfLink aLinks; - getLinks( theNextFace, aLinks ); - int aNbLink = (int)aLinks.size(); - for ( int i = 0; i < aNbLink; i++ ) + if (myMeshDS == 0 || myShape.IsNull()) + return false; + + if (!myIsSubshape) { - ManifoldPart::Link aLink = aLinks[ i ]; - if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) ) - continue; - if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() ) + return myElementsOnShapePtr->IsSatisfy(theId); + } + + // Case of submesh + if (myType == SMDSAbs_Node) + { + if( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ) ) { - if ( myIsOnlyManifold ) + if ( aNode->getshapeId() < 1 ) + return myElementsOnShapePtr->IsSatisfy(theId); + + const SMDS_PositionPtr& aPosition = aNode->GetPosition(); + SMDS_TypeOfPosition aTypeOfPosition = aPosition->GetTypeOfPosition(); + switch( aTypeOfPosition ) { - // remove from boundary - theMapOfBoundary.erase( aLink ); - ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin(); - for ( ; pLink != theSeqOfBoundary.end(); ++pLink ) - { - ManifoldPart::Link aBoundLink = *pLink; - if ( aBoundLink.IsEqual( aLink ) ) - { - theSeqOfBoundary.erase( pLink ); - break; - } - } + case SMDS_TOP_VERTEX : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_VERTEX )); + case SMDS_TOP_EDGE : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_EDGE )); + case SMDS_TOP_FACE : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_FACE )); + case SMDS_TOP_3DSPACE: return ( IsContains( myMeshDS,myShape,aNode,TopAbs_SOLID ) || + IsContains( myMeshDS,myShape,aNode,TopAbs_SHELL )); } } - else + } + else + { + if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId )) { - theMapOfBoundary.insert( aLink ); - theSeqOfBoundary.push_back( aLink ); - theDMapLinkFacePtr[ aLink ] = theNextFace; + if ( anElem->getshapeId() < 1 ) + return myElementsOnShapePtr->IsSatisfy(theId); + + if( myType == SMDSAbs_All ) + { + return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ) || + IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ) || + IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )|| + IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL )); + } + else if( myType == anElem->GetType() ) + { + switch( myType ) + { + case SMDSAbs_Edge : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE )); + case SMDSAbs_Face : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_FACE )); + case SMDSAbs_Volume: return ( IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )|| + IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL )); + } + } } } + + return false; } -void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink, - ManifoldPart::TVectorOfFacePtr& theFaces ) const +void BelongToGeom::SetType (SMDSAbs_ElementType theType) { - SMDS_Mesh::SetOfFaces aSetOfFaces; - // take all faces that shared first node - SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator(); - for ( ; anItr->more(); ) - { - SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next(); - if ( !aFace ) - continue; - aSetOfFaces.Add( aFace ); - } - // take all faces that shared second node - anItr = theLink.myNode2->facesIterator(); - // find the common part of two sets - for ( ; anItr->more(); ) - { - SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next(); - if ( aSetOfFaces.Contains( aFace ) ) - theFaces.push_back( aFace ); - } + myType = theType; + init(); } +SMDSAbs_ElementType BelongToGeom::GetType() const +{ + return myType; +} -/* - ElementsOnSurface -*/ - -ElementsOnSurface::ElementsOnSurface() +TopoDS_Shape BelongToGeom::GetShape() { - myMesh = 0; - myIds.Clear(); - myType = SMDSAbs_All; - mySurf.Nullify(); - myToler = Precision::Confusion(); + return myShape; } -ElementsOnSurface::~ElementsOnSurface() +const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const { - myMesh = 0; + return myMeshDS; } -void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh ) +void BelongToGeom::SetTolerance (double theTolerance) { - if ( myMesh == theMesh ) - return; - myMesh = theMesh; - myIds.Clear(); - process(); + myTolerance = theTolerance; + if (!myIsSubshape) + init(); } -bool ElementsOnSurface::IsSatisfy( long theElementId ) +double BelongToGeom::GetTolerance() { - return myIds.Contains( theElementId ); + return myTolerance; } -SMDSAbs_ElementType ElementsOnSurface::GetType() const -{ return myType; } +/* + Class : LyingOnGeom + Description : Predicate for verifying whether entiy lying or partially lying on + specified geometrical support +*/ -void ElementsOnSurface::SetTolerance( const double theToler ) -{ myToler = theToler; } +LyingOnGeom::LyingOnGeom() + : myMeshDS(NULL), + myType(SMDSAbs_All), + myIsSubshape(false), + myTolerance(Precision::Confusion()) +{} -double ElementsOnSurface::GetTolerance() const +void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh ) { - return myToler; + myMeshDS = dynamic_cast(theMesh); + init(); } -void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape, - const SMDSAbs_ElementType theType ) +void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape ) { - myType = theType; - mySurf.Nullify(); - if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE ) - { - mySurf.Nullify(); - return; - } - TopoDS_Face aFace = TopoDS::Face( theShape ); - mySurf = BRep_Tool::Surface( aFace ); + myShape = theShape; + init(); } -void ElementsOnSurface::process() +void LyingOnGeom::init() { - myIds.Clear(); - if ( mySurf.IsNull() ) - return; - - if ( myMesh == 0 ) - return; + if (!myMeshDS || myShape.IsNull()) return; - if ( myType == SMDSAbs_Face || myType == SMDSAbs_All ) - { - SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); - for(; anIter->more(); ) - process( anIter->next() ); + // is sub-shape of main shape? + TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh(); + if (aMainShape.IsNull()) { + myIsSubshape = false; + } + else { + myIsSubshape = myMeshDS->IsGroupOfSubShapes( myShape ); } - if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All ) + if (myIsSubshape) { - SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator(); - for(; anIter->more(); ) - process( anIter->next() ); + 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 ); + } } - - if ( myType == SMDSAbs_Node ) + else { - SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator(); - for(; anIter->more(); ) - process( anIter->next() ); + myElementsOnShapePtr.reset(new ElementsOnShape()); + myElementsOnShapePtr->SetTolerance(myTolerance); + myElementsOnShapePtr->SetAllNodes(false); // lays on, while true means "belong" + myElementsOnShapePtr->SetMesh(myMeshDS); + myElementsOnShapePtr->SetShape(myShape, myType); } } -void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr ) +bool LyingOnGeom::IsSatisfy( long theId ) { - SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator(); - bool isSatisfy = true; - for ( ; aNodeItr->more(); ) + if ( myMeshDS == 0 || myShape.IsNull() ) + return false; + + if (!myIsSubshape) { - SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); - if ( !isOnSurface( aNode ) ) + 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 ) + { + SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator(); + while ( nodeItr->more() ) { - isSatisfy = false; - break; + const SMDS_MeshElement* aNode = nodeItr->next(); + if ( mySubShapesIDs.Contains( aNode->getshapeId() )) + return true; } } - if ( isSatisfy ) - myIds.Add( theElemPtr->GetID() ); + + return false; } -bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) const +void LyingOnGeom::SetType( SMDSAbs_ElementType theType ) { - if ( mySurf.IsNull() ) - return false; + myType = theType; + init(); +} - gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() ); - double aToler2 = myToler * myToler; - if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane))) - { - gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln(); - if ( aPln.SquareDistance( aPnt ) > aToler2 ) - return false; - } - else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface))) +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; + if (!myIsSubshape) + init(); +} + +double LyingOnGeom::GetTolerance() +{ + return myTolerance; +} + +bool LyingOnGeom::Contains( const SMESHDS_Mesh* theMeshDS, + const TopoDS_Shape& theShape, + const SMDS_MeshElement* theElem, + TopAbs_ShapeEnum theFindShapeEnum, + TopAbs_ShapeEnum theAvoidShapeEnum ) +{ + // if (IsContains(theMeshDS, theShape, theElem, theFindShapeEnum, theAvoidShapeEnum)) + // return true; + + // TopTools_MapOfShape aSubShapes; + // TopExp_Explorer exp( theShape, theFindShapeEnum, theAvoidShapeEnum ); + // for ( ; exp.More(); exp.Next() ) + // { + // const TopoDS_Shape& aShape = exp.Current(); + // if ( !aSubShapes.Add( aShape )) continue; + + // if ( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape )) + // { + // if ( aSubMesh->Contains( theElem )) + // return true; + + // SMDS_ElemIteratorPtr nodeItr = theElem->nodesIterator(); + // while ( nodeItr->more() ) + // { + // const SMDS_MeshElement* aNode = nodeItr->next(); + // if ( aSubMesh->Contains( aNode )) + // return true; + // } + // } + // } + return false; +} + +TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0) +{} + +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 ) { - gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder(); - double aRad = aCyl.Radius(); - gp_Ax3 anAxis = aCyl.Position(); - gp_XYZ aLoc = aCyl.Location().XYZ(); - double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc ); - double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc ); - if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 ) - return false; + modified = ( myMeshModifTime != myMesh->GetMTime() ); + myMeshModifTime = myMesh->GetMTime(); } - else - return false; - - return true; + return modified; }