X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=6b400350d09daae940b87994fc36ecec6be6ac76;hb=f0bf265fc0ea849ad1a359e10ad0c9cf6b0b9bb8;hp=06bfd5304a9fb4422a87b325dee2880818a9cc41;hpb=c3bf92bd87b770fd81631a3853f7f5bb1ac6a4e8;p=modules%2Fsmesh.git diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 06bfd5304..6b400350d 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -1,111 +1,212 @@ -// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, -// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS -// -// This library is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 2.1 of the License. -// -// This library is 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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org - -#include - +// Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is 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 +// + +#include "SMESH_ControlsDef.hxx" + +#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_Mesh.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 "SMDS_Mesh.hxx" -#include "SMDS_Iterator.hxx" -#include "SMDS_MeshElement.hxx" -#include "SMDS_MeshNode.hxx" +#include + +#include +#include -#include "SMESH_Controls.hxx" /* - AUXILIARY METHODS + AUXILIARY METHODS */ -static inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) -{ - return gp_Vec( P1 - P2 ).Angle( gp_Vec( P3 - P2 ) ); -} +namespace{ -static inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) -{ - gp_Vec aVec1( P2 - P1 ); - gp_Vec aVec2( P3 - P1 ); - return ( aVec1 ^ aVec2 ).Magnitude() * 0.5; -} + inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode ) + { + return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() ); + } -static inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 ) -{ - return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() ); -} + inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) + { + gp_Vec v1( P1 - P2 ), v2( P3 - P2 ); -static inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 ) -{ - double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) ); - return aDist; -} + return v1.Magnitude() < gp::Resolution() || + v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 ); + } -static int getNbMultiConnection( SMDS_Mesh* theMesh, const int theId ) -{ - if ( theMesh == 0 ) - return 0; + inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 ) + { + gp_Vec aVec1( P2 - P1 ); + gp_Vec aVec2( P3 - P1 ); + return ( aVec1 ^ aVec2 ).Magnitude() * 0.5; + } - const SMDS_MeshElement* anEdge = theMesh->FindElement( theId ); - if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 ) - return 0; + inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 ) + { + return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() ); + } - TColStd_MapOfInteger aMap; - int aResult = 0; - SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator(); - if ( anIter != 0 ) + + inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 ) { - 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(); + double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) ); + return aDist; + } - if ( anIter->more() ) // i.e. first node - aMap.Add( anId ); - else if ( aMap.Contains( anId ) ) - aResult++; + int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId ) + { + if ( theMesh == 0 ) + return 0; + + const SMDS_MeshElement* anEdge = theMesh->FindElement( theId ); + if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */) + return 0; + + // 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 ); + +// TColStd_MapOfInteger aMap; + +// SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator(); +// if ( anIter != 0 ) { +// while( anIter->more() ) { +// const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); +// if ( aNode == 0 ) +// return 0; +// SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator(); +// while( anElemIter->more() ) { +// const SMDS_MeshElement* anElem = anElemIter->next(); +// if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) { +// int anId = anElem->GetID(); + +// if ( anIter->more() ) // i.e. first node +// aMap.Add( anId ); +// else if ( aMap.Contains( anId ) ) +// aResult++; +// } +// } +// } +// } + + return aResult; } - 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; + } } + using namespace SMESH::Controls; /* - FUNCTORS -*/ + * FUNCTORS + */ /* Class : NumericalFunctor @@ -113,662 +214,3732 @@ using namespace SMESH::Controls; */ NumericalFunctor::NumericalFunctor(): myMesh(NULL) -{} +{ + myPrecision = -1; +} -void NumericalFunctor::SetMesh( SMDS_Mesh* theMesh ) +void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh ) { myMesh = theMesh; } -bool NumericalFunctor::getPoints( const int theId, - TColgp_SequenceOfXYZ& theRes ) const +bool NumericalFunctor::GetPoints(const int theId, + TSequenceOfXYZ& theRes ) const { - theRes.Clear(); + theRes.clear(); if ( myMesh == 0 ) return false; - // Get nodes of the face - const SMDS_MeshElement* anElem = myMesh->FindElement( theId ); - if ( anElem == 0 || anElem->GetType() != GetType() ) + return GetPoints( myMesh->FindElement( theId ), theRes ); +} + +bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, + TSequenceOfXYZ& theRes ) +{ + theRes.clear(); + + if ( anElem == 0) return false; - int nbNodes = anElem->NbNodes(); + theRes.reserve( anElem->NbNodes() ); + + // Get nodes of the element + SMDS_ElemIteratorPtr anIter; + + if ( anElem->IsQuadratic() ) { + switch ( anElem->GetType() ) { + case SMDSAbs_Edge: + anIter = dynamic_cast + (anElem)->interlacedNodesElemIterator(); + break; + case SMDSAbs_Face: + anIter = dynamic_cast + (anElem)->interlacedNodesElemIterator(); + break; + default: + anIter = anElem->nodesIterator(); + //return false; + } + } + else { + anIter = anElem->nodesIterator(); + } - SMDS_ElemIteratorPtr anIter = anElem->nodesIterator(); - if ( anIter != 0 ) - { - while( anIter->more() ) - { - const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next(); - if ( aNode != 0 ) - theRes.Append( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); + if ( anIter ) { + while( anIter->more() ) { + if ( const SMDS_MeshNode* aNode = static_cast( anIter->next() )) + theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) ); } } return true; } - -/* - Class : MinimumAngle - Description : Functor for calculation of minimum angle -*/ -double MinimumAngle::GetValue( long theId ) +long NumericalFunctor::GetPrecision() const { - TColgp_SequenceOfXYZ P; - if ( !getPoints( theId, P ) || P.Length() != 3 && P.Length() != 4 ) - return 0; - - double aMin; - - if ( P.Length() == 3 ) - { - double A0 = getAngle( P( 3 ), P( 1 ), P( 2 ) ); - double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) ); - double A2 = getAngle( P( 2 ), P( 3 ), P( 1 ) ); - - aMin = Min( A0, Min( A1, A2 ) ); - } - else - { - double A0 = getAngle( P( 4 ), P( 1 ), P( 2 ) ); - double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) ); - double A2 = getAngle( P( 2 ), P( 3 ), P( 4 ) ); - double A3 = getAngle( P( 3 ), P( 4 ), P( 1 ) ); - - aMin = Min( Min( A0, A1 ), Min( A2, A3 ) ); - } - - return aMin * 180 / PI; + return myPrecision; } -SMDSAbs_ElementType MinimumAngle::GetType() const +void NumericalFunctor::SetPrecision( const long thePrecision ) { - return SMDSAbs_Face; + myPrecision = thePrecision; + myPrecisionValue = pow( 10., (double)( myPrecision ) ); } - -/* - Class : AspectRatio - Description : Functor for calculating aspect ratio -*/ -double AspectRatio::GetValue( long theId ) +double NumericalFunctor::GetValue( long theId ) { - TColgp_SequenceOfXYZ P; - if ( !getPoints( theId, P ) || P.Length() != 3 && P.Length() != 4 ) - return 0; + double aVal = 0; - int nbNodes = P.Length(); + myCurrElement = myMesh->FindElement( theId ); - // Compute lengths of the sides + TSequenceOfXYZ P; + if ( GetPoints( theId, P )) + aVal = Round( GetValue( P )); - 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 ) ); + return aVal; +} - // Compute aspect ratio +double NumericalFunctor::Round( const double & aVal ) +{ + return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal; +} - if ( nbNodes == 3 ) +//================================================================================ +/*! + * \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) +{ + 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 aMaxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) ); - double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) ); - static double aCoef = sqrt( 3. ) / 4; + 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 )); + } - return anArea != 0 ? aCoef * aMaxLen * aMaxLen / anArea : 0; + if ( minmax ) + { + funValues[0] = minmax[0]; + funValues[nbIntervals] = minmax[1]; } else { - double aMaxLen = Max( Max( aLen[ 0 ], aLen[ 1 ] ), Max( aLen[ 2 ], aLen[ 3 ] ) ); - double aMinLen = Min( Min( aLen[ 0 ], aLen[ 1 ] ), Min( aLen[ 2 ], aLen[ 3 ] ) ); - - return aMinLen != 0 ? aMaxLen / aMinLen : 0; + 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 ); + funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r; -SMDSAbs_ElementType AspectRatio::GetType() const -{ - return SMDSAbs_Face; + // count values in the i-th interval if there are any + if ( min != values.end() && *min <= funValues[i+1] ) + { + // 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; + } + } + // add values larger than minmax[1] + nbEvents.back() += std::distance( min, values.end() ); } +//======================================================================= +//function : GetValue +//purpose : +//======================================================================= -/* - Class : Warping - Description : Functor for calculating warping -*/ -double Warping::GetValue( long theId ) +double Volume::GetValue( long theElementId ) { - TColgp_SequenceOfXYZ P; - if ( !getPoints( theId, P ) || P.Length() != 4 ) - return 0; - - gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4; - - double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G ); - double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G ); - double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G ); - double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G ); - - return Max( Max( A1, A2 ), Max( A3, A4 ) ); + if ( theElementId && myMesh ) { + SMDS_VolumeTool aVolumeTool; + if ( aVolumeTool.Set( myMesh->FindElement( theElementId ))) + return aVolumeTool.GetSize(); + } + return 0; } -double Warping::ComputeA( const gp_XYZ& thePnt1, - const gp_XYZ& thePnt2, - const gp_XYZ& thePnt3, - const gp_XYZ& theG ) const -{ - 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; - - gp_XYZ GI = ( thePnt2 - thePnt1 ) / 2. - theG; - gp_XYZ GJ = ( thePnt3 - thePnt2 ) / 2. - theG; - gp_XYZ N = GI.Crossed( GJ ); - N.Normalize(); +//======================================================================= +//function : GetBadRate +//purpose : meaningless as it is not quality control functor +//======================================================================= - double H = gp_Vec( thePnt2 - theG ).Dot( gp_Vec( N ) ); - return asin( fabs( H / L ) ) * 180 / PI; +double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; } -SMDSAbs_ElementType Warping::GetType() const +//======================================================================= +//function : GetType +//purpose : +//======================================================================= + +SMDSAbs_ElementType Volume::GetType() const { - return SMDSAbs_Face; + return SMDSAbs_Volume; } /* - Class : Taper - Description : Functor for calculating taper + Class : MaxElementLength2D + Description : Functor calculating maximum length of 2D element */ -double Taper::GetValue( long theId ) +double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P ) { - TColgp_SequenceOfXYZ P; - if ( !getPoints( theId, P ) || P.Length() != 4 ) - return 0; - - // Compute taper - double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2; - double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2; - double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2; - double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2; + 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)); + } - double JA = 0.25 * ( J1 + J2 + J3 + J4 ); + if( myPrecision >= 0 ) + { + double prec = pow( 10., (double)myPrecision ); + aVal = floor( aVal * prec + 0.5 ) / prec; + } + return aVal; +} - double T1 = fabs( ( J1 - JA ) / JA ); - double T2 = fabs( ( J2 - JA ) / JA ); - double T3 = fabs( ( J3 - JA ) / JA ); - double T4 = fabs( ( J4 - JA ) / JA ); +double MaxElementLength2D::GetValue( long theElementId ) +{ + TSequenceOfXYZ P; + if( GetPoints( theElementId, P ) && myMesh->FindElement( theElementId )->GetType() == SMDSAbs_Face) { + GetValue(P); + } + return 0.; +} - return Max( Max( T1, T2 ), Max( T3, T4 ) ); +double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; } -SMDSAbs_ElementType Taper::GetType() const +SMDSAbs_ElementType MaxElementLength2D::GetType() const { return SMDSAbs_Face; } - /* - Class : Skew - Description : Functor for calculating skew in degrees + Class : MaxElementLength3D + Description : Functor calculating maximum length of 3D element */ -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; - - return gp_Vec( p31 - p2 ).Angle( p12 - p23 ); -} -double Skew::GetValue( long theId ) +double MaxElementLength3D::GetValue( long theElementId ) { - TColgp_SequenceOfXYZ P; - if ( !getPoints( theId, P ) || P.Length() != 3 && P.Length() != 4 ) - return 0; - - // Compute skew - static double PI2 = PI / 2; - if ( P.Length() == 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 ) ) ); + 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 ); + } + } + } + } + } - return Max( A0, Max( A1, A2 ) ) * 180 / PI; + if( myPrecision >= 0 ) + { + double prec = pow( 10., (double)myPrecision ); + aVal = floor( aVal * prec + 0.5 ) / prec; + } + return aVal; } - 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; - - double A = fabs( PI2 - gp_Vec( p34 - p12 ).Angle( p23 - p41 ) ); + return 0.; +} - return A * 180 / PI; - } +double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; } -SMDSAbs_ElementType Skew::GetType() const +SMDSAbs_ElementType MaxElementLength3D::GetType() const { - return SMDSAbs_Face; + return SMDSAbs_Volume; } /* - Class : Area - Description : Functor for calculating area + Class : MinimumAngle + Description : Functor for calculation of minimum angle */ -double Area::GetValue( long theId ) + +double MinimumAngle::GetValue( const TSequenceOfXYZ& P ) { - TColgp_SequenceOfXYZ P; - if ( !getPoints( theId, P ) || P.Length() != 3 && P.Length() != 4 ) - return 0; + double aMin; - if ( P.Length() == 3 ) - return getArea( P( 1 ), P( 2 ), P( 3 ) ); - else - return getArea( P( 1 ), P( 2 ), P( 3 ) ) + getArea( P( 1 ), P( 3 ), P( 4 ) ); -} + if (P.size() <3) + return 0.; -SMDSAbs_ElementType Area::GetType() const -{ - return SMDSAbs_Face; -} + 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; } -SMDSAbs_ElementType MultiConnection::GetType() const +double AspectRatio::GetValue( const TSequenceOfXYZ& P ) { - return SMDSAbs_Edge; -} + // 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(); -/* - PREDICATES -*/ + if ( nbNodes < 3 ) + return 0; -/* - Class : FreeBorders - Description : Predicate for free borders -*/ + // Compute aspect ratio -FreeBorders::FreeBorders() -{ - myMesh = 0; + 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.; + return alfa * maxLen * half_perimeter / anArea; + } + else if ( nbNodes == 6 ) { // quadratic triangles + // Compute lengths of the sides + std::vector< double > aLen (3); + aLen[0] = getDistance( P(1), P(3) ); + aLen[1] = getDistance( P(3), P(5) ); + aLen[2] = getDistance( P(5), P(1) ); + // Q = alfa * h * p / S, where + // + // alfa = sqrt( 3 ) / 6 + // h - length of the longest edge + // p - half perimeter + // S - triangle surface + const double alfa = sqrt( 3. ) / 6.; + double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) ); + double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.; + double anArea = getArea( P(1), P(3), P(5) ); + if ( anArea <= Precision::Confusion() ) + return 0.; + 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 <= Precision::Confusion() ) + return 0.; + 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 + 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 <= Precision::Confusion() ) + return 0.; + return alpha * L * C1 / C2; + } + return 0; } -void FreeBorders::SetMesh( SMDS_Mesh* theMesh ) +double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const { - myMesh = theMesh; + // the aspect ratio is in the range [1.0,infinity] + // < 1.0 = very bad, zero area + // 1.0 = good + // infinity = bad + return ( Value < 0.9 ) ? 1000 : Value / 1000.; } -bool FreeBorders::IsSatisfy( long theId ) +SMDSAbs_ElementType AspectRatio::GetType() const { - return getNbMultiConnection( myMesh, theId ) == 1; + return SMDSAbs_Face; +} + + +/* + Class : AspectRatio3D + Description : Functor for calculating aspect ratio +*/ +namespace{ + + inline double getHalfPerimeter(double theTria[3]){ + return (theTria[0] + theTria[1] + theTria[2])/2.0; + } + + inline double getArea(double theHalfPerim, double theTria[3]){ + return sqrt(theHalfPerim* + (theHalfPerim-theTria[0])* + (theHalfPerim-theTria[1])* + (theHalfPerim-theTria[2])); + } + + inline double getVolume(double theLen[6]){ + double a2 = theLen[0]*theLen[0]; + double b2 = theLen[1]*theLen[1]; + double c2 = theLen[2]*theLen[2]; + double d2 = theLen[3]*theLen[3]; + double e2 = theLen[4]*theLen[4]; + double f2 = theLen[5]*theLen[5]; + double P = 4.0*a2*b2*d2; + double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2); + double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2); + return sqrt(P-Q+R)/12.0; + } + + inline double getVolume2(double theLen[6]){ + double a2 = theLen[0]*theLen[0]; + double b2 = theLen[1]*theLen[1]; + double c2 = theLen[2]*theLen[2]; + double d2 = theLen[3]*theLen[3]; + double e2 = theLen[4]*theLen[4]; + double f2 = theLen[5]*theLen[5]; + + double P = a2*e2*(b2+c2+d2+f2-a2-e2); + double Q = b2*f2*(a2+c2+d2+e2-b2-f2); + double R = c2*d2*(a2+b2+e2+f2-c2-d2); + double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2; + + return sqrt(P+Q+R-S)/12.0; + } + + inline double getVolume(const TSequenceOfXYZ& P){ + gp_Vec aVec1( P( 2 ) - P( 1 ) ); + gp_Vec aVec2( P( 3 ) - P( 1 ) ); + gp_Vec aVec3( P( 4 ) - P( 1 ) ); + gp_Vec anAreaVec( aVec1 ^ aVec2 ); + return fabs(aVec3 * anAreaVec) / 6.0; + } + + inline double getMaxHeight(double theLen[6]) + { + 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(); + + 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 + getDistance(P( 2 ),P( 3 )), // b + getDistance(P( 3 ),P( 1 )), // c + getDistance(P( 2 ),P( 4 )), // d + getDistance(P( 3 ),P( 4 )), // e + getDistance(P( 1 ),P( 4 )) // f + }; + double aTria[4][3] = { + {aLen[0],aLen[1],aLen[2]}, // abc + {aLen[0],aLen[3],aLen[5]}, // adf + {aLen[1],aLen[3],aLen[4]}, // bde + {aLen[2],aLen[4],aLen[5]} // cef + }; + double aSumArea = 0.0; + double aHalfPerimeter = getHalfPerimeter(aTria[0]); + double anArea = getArea(aHalfPerimeter,aTria[0]); + aSumArea += anArea; + aHalfPerimeter = getHalfPerimeter(aTria[1]); + anArea = getArea(aHalfPerimeter,aTria[1]); + aSumArea += anArea; + aHalfPerimeter = getHalfPerimeter(aTria[2]); + anArea = getArea(aHalfPerimeter,aTria[2]); + aSumArea += anArea; + aHalfPerimeter = getHalfPerimeter(aTria[3]); + anArea = getArea(aHalfPerimeter,aTria[3]); + aSumArea += anArea; + double aVolume = getVolume(P); + //double aVolume = getVolume(aLen); + double aHeight = getMaxHeight(aLen); + static double aCoeff = sqrt(2.0)/12.0; + 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 = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )}; + 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 = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )}; + 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 = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )}; + aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + } + { + gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )}; + 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; + SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes ); + int nbFaces = SMDS_VolumeTool::NbFaces( type ); + TSequenceOfXYZ points(4); + for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume + if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 ) + continue; + const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true ); + for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face + points( p + 1 ) = P( pInd[ p ] + 1 ); + aQuality = std::max( aQuality, aspect2D.GetValue( points )); + } + } + return aQuality; +} + +double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + // the aspect ratio is in the range [1.0,infinity] + // 1.0 = good + // infinity = bad + return Value / 1000.; +} + +SMDSAbs_ElementType AspectRatio3D::GetType() const +{ + return SMDSAbs_Volume; +} + + +/* + 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.; + + 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 Warping::ComputeA( const gp_XYZ& thePnt1, + const gp_XYZ& thePnt2, + const gp_XYZ& thePnt3, + const gp_XYZ& theG ) const +{ + 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.; + + 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 M_PI / 2; + + N.Normalize(); + + double H = ( thePnt2 - theG ).Dot( N ); + return asin( fabs( H / L ) ) * 180. / M_PI; +} + +double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + // the warp is in the range [0.0,PI/2] + // 0.0 = good (no warp) + // PI/2 = bad (face pliee) + return Value; +} + +SMDSAbs_ElementType Warping::GetType() const +{ + return SMDSAbs_Face; +} + + +/* + Class : Taper + Description : Functor for calculating taper +*/ +double Taper::GetValue( const TSequenceOfXYZ& P ) +{ + if ( P.size() != 4 ) + return 0.; + + // Compute taper + double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2.; + double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2.; + double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2.; + double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2.; + + double JA = 0.25 * ( J1 + J2 + J3 + J4 ); + if ( JA <= Precision::Confusion() ) + return 0.; + + 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 Taper::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + // the taper is in the range [0.0,1.0] + // 0.0 = good (no taper) + // 1.0 = bad (les cotes opposes sont allignes) + return Value; +} + +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_Vec v1( p31 - p2 ), v2( p12 - p23 ); + + 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.; + + // Compute skew + static 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. / 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_Vec v1( p34 - p12 ), v2( p23 - p41 ); + double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution() + ? 0. : fabs( PI2 - v1.Angle( v2 ) ); + + //BUG SWP12743 + if ( A < Precision::Angular() ) + return 0.; + + return A * 180. / M_PI; + } +} + +double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + // the skew is in the range [0.0,PI/2]. + // 0.0 = good + // PI/2 = bad + return Value; +} + +SMDSAbs_ElementType Skew::GetType() const +{ + return SMDSAbs_Face; +} + + +/* + Class : Area + Description : Functor for calculating area +*/ +double Area::GetValue( const TSequenceOfXYZ& P ) +{ + 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 a quality control functor + return Value; +} + +SMDSAbs_ElementType Area::GetType() const +{ + return SMDSAbs_Face; +} + + +/* + Class : Length + Description : Functor for calculating length of edge +*/ +double Length::GetValue( const TSequenceOfXYZ& P ) +{ + 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 +{ + // meaningless as it is not quality control functor + return Value; +} + +SMDSAbs_ElementType Length::GetType() const +{ + return SMDSAbs_Edge; +} + +/* + Class : Length2D + Description : Functor for calculating length of edge +*/ + +double Length2D::GetValue( long theElementId) +{ + TSequenceOfXYZ P; + + //cout<<"Length2D::GetValue"<FindElement( theElementId ); + SMDSAbs_ElementType aType = aElem->GetType(); + + int len = P.size(); + + switch (aType){ + case SMDSAbs_All: + case SMDSAbs_Node: + case SMDSAbs_Edge: + if (len == 2){ + aVal = getDistance( P( 1 ), P( 2 ) ); + break; + } + else if (len == 3){ // quadratic edge + aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 )); + break; + } + case SMDSAbs_Face: + if (len == 3){ // triangles + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 1 )); + aVal = Max(L1,Max(L2,L3)); + break; + } + else if (len == 4){ // quadrangles + double L1 = getDistance(P( 1 ),P( 2 )); + double L2 = getDistance(P( 2 ),P( 3 )); + double L3 = getDistance(P( 3 ),P( 4 )); + double L4 = getDistance(P( 4 ),P( 1 )); + aVal = Max(Max(L1,L2),Max(L3,L4)); + break; + } + if (len == 6){ // quadratic triangles + double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 )); + double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 )); + double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 )); + aVal = Max(L1,Max(L2,L3)); + //cout<<"L1="<FindElement( theId ); + if ( !anElem ) + return false; + const SMDSAbs_ElementType anElemType = anElem->GetType(); + if ( myType != SMDSAbs_All && anElemType != myType ) + return false; + bool isOk = ( anElem->GetGeomType() == myGeomType ); + return isOk; +} + +void ElemGeomType::SetType( SMDSAbs_ElementType theType ) +{ + myType = theType; +} + +SMDSAbs_ElementType ElemGeomType::GetType() const +{ + return myType; +} + +void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType ) +{ + myGeomType = theType; +} + +SMDSAbs_GeometryType ElemGeomType::GetGeomType() const +{ + return myGeomType; +} + +//================================================================================ +/*! + * \brief Class CoplanarFaces + */ +//================================================================================ + +CoplanarFaces::CoplanarFaces() + : myFaceID(0), myToler(0) +{ +} +void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified() ) + { + // Build a set of coplanar face ids + + myCoplanarIDs.clear(); + + if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler ) + return; + + const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID ); + if ( !face || face->GetType() != SMDSAbs_Face ) + return; + + bool normOK; + gp_Vec myNorm = getNormale( static_cast(face), &normOK ); + if (!normOK) + return; + + const double radianTol = myToler * M_PI / 180.; + std::set< SMESH_TLink > checkedLinks; + + 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(); + + 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 CoplanarFaces::IsSatisfy( long theElementId ) +{ + return myCoplanarIDs.count( theElementId ); +} + +/* + *Class : RangeOfIds + *Description : Predicate for Range of Ids. + * Range may be specified with two ways. + * 1. Using AddToRange method + * 2. With SetRangeStr method. Parameter of this method is a string + * like as "1,2,3,50-60,63,67,70-" +*/ + +//======================================================================= +// name : RangeOfIds +// Purpose : Constructor +//======================================================================= +RangeOfIds::RangeOfIds() +{ + myMesh = 0; + myType = SMDSAbs_All; +} + +//======================================================================= +// name : SetMesh +// Purpose : Set mesh +//======================================================================= +void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh ) +{ + myMesh = theMesh; +} + +//======================================================================= +// name : AddToRange +// Purpose : Add ID to the range +//======================================================================= +bool RangeOfIds::AddToRange( long theEntityId ) +{ + myIds.Add( theEntityId ); + return true; +} + +//======================================================================= +// name : GetRangeStr +// Purpose : Get range as a string. +// Example: "1,2,3,50-60,63,67,70-" +//======================================================================= +void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr ) +{ + theResStr.Clear(); + + 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 ); + + TCollection_AsciiString aStr; + if ( aMinId != IntegerFirst() ) + aStr += aMinId; + + aStr += "-"; + + 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 aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) ) + aStr.Remove( aPos, 2 ); + + 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 : FreeEdges - Description : Predicate for free Edges + Class : LessThan + Description : Comparator "<" */ -FreeEdges::FreeEdges() +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 +*/ + +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 FreeEdges::SetMesh( SMDS_Mesh* theMesh ) -{ - myMesh = theMesh; +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() ) + { + 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(); } -bool FreeEdges::IsSatisfy( long theId ) +static void getLinks( const SMDS_MeshFace* theFace, + ManifoldPart::TVectorOfLink& theLinks ) { - return getNbMultiConnection( myMesh, theId ) == 1; + 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 ); + } } -SMDSAbs_ElementType FreeEdges::GetType() const +bool ManifoldPart::findConnected + ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt, + SMDS_MeshFace* theStartFace, + ManifoldPart::TMapOfLink& theNonManifold, + TColStd_MapOfInteger& theResFaces ) { - return SMDSAbs_Face; -} + theResFaces.Clear(); + if ( !theAllFacePtrInt.size() ) + return false; -FreeEdges::Border::Border(long thePntId1, long thePntId2){ - PntId[0] = thePntId1; PntId[1] = thePntId2; - if(thePntId1 > thePntId2){ - PntId[1] = thePntId1; PntId[0] = thePntId2; + if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() ) + { + myMapBadGeomIds.Add( theStartFace->GetID() ); + return false; } -} -//bool operator<(const FreeEdges::Border& x, const FreeEdges::Border& y){ -// if(x.PntId[0] < y.PntId[0]) return true; -// if(x.PntId[0] == y.PntId[0]) -// if(x.PntId[1] < y.PntId[1]) return true; -// return false; -//} - -namespace SMESH{ - namespace Controls{ - struct EdgeBorder: public FreeEdges::Border{ - long ElemId; - EdgeBorder(long theElemId, long thePntId1, long thePntId2): - FreeEdges::Border(thePntId1,thePntId2), - ElemId(theElemId) - {} - }; - - bool operator<(const FreeEdges::Border& x, const FreeEdges::Border& y){ - if(x.PntId[0] < y.PntId[0]) return true; - if(x.PntId[0] == y.PntId[0]) - if(x.PntId[1] < y.PntId[1]) return true; - return false; - } + ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip; + ManifoldPart::TVectorOfLink aSeqOfBoundary; + theResFaces.Add( theStartFace->GetID() ); + ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace; - typedef std::set EdgeBorderS; + expandBoundary( aMapOfBoundary, aSeqOfBoundary, + aDMapLinkFace, theNonManifold, theStartFace ); - inline void UpdateBorders(const EdgeBorder& theBorder, - EdgeBorderS& theRegistry, - EdgeBorderS& theContainer) + bool isDone = false; + while ( !isDone && aMapOfBoundary.size() != 0 ) + { + bool isToReset = false; + ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin(); + for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink ) { - if(theRegistry.find(theBorder) == theRegistry.end()){ - theRegistry.insert(theBorder); - theContainer.insert(theBorder); - }else{ - theContainer.erase(theBorder); + 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 FreeEdges::GetBoreders(Borders& theBorders) +void ManifoldPart::expandBoundary + ( ManifoldPart::TMapOfLink& theMapOfBoundary, + ManifoldPart::TVectorOfLink& theSeqOfBoundary, + ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr, + ManifoldPart::TMapOfLink& theNonManifold, + SMDS_MeshFace* theNextFace ) const { - EdgeBorderS aRegistry; - EdgeBorderS aContainer; - SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); - for(; anIter->more(); ){ - const SMDS_MeshFace* anElem = anIter->next(); - long anElemId = anElem->GetID(); - SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator(); - long aNodeId[2]; - const SMDS_MeshElement* aNode; - if(aNodesIter->more()){ - aNode = aNodesIter->next(); - aNodeId[0] = aNodeId[1] = aNode->GetID(); - } - for(; aNodesIter->more(); ){ - aNode = aNodesIter->next(); - EdgeBorder aBorder(anElemId,aNodeId[1],aNode->GetID()); - aNodeId[1] = aNode->GetID(); - //std::cout< 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 : Comparator - Description : Base class for comparators + ElementsOnSurface */ -Comparator::Comparator(): - myMargin(0) -{} - -Comparator::~Comparator() -{} -void Comparator::SetMesh( SMDS_Mesh* theMesh ) +ElementsOnSurface::ElementsOnSurface() { - if ( myFunctor ) - myFunctor->SetMesh( theMesh ); + myIds.Clear(); + myType = SMDSAbs_All; + mySurf.Nullify(); + myToler = Precision::Confusion(); + myUseBoundaries = false; } -void Comparator::SetMargin( double theValue ) +ElementsOnSurface::~ElementsOnSurface() { - myMargin = theValue; } -void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct ) +void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh ) { - myFunctor = theFunct; + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified()) + process(); } -SMDSAbs_ElementType Comparator::GetType() const +bool ElementsOnSurface::IsSatisfy( long theElementId ) { - return myFunctor ? myFunctor->GetType() : SMDSAbs_All; + return myIds.Contains( theElementId ); } +SMDSAbs_ElementType ElementsOnSurface::GetType() const +{ return myType; } -/* - Class : LessThan - Description : Comparator "<" -*/ -bool LessThan::IsSatisfy( long theId ) +void ElementsOnSurface::SetTolerance( const double theToler ) { - return myFunctor && myFunctor->GetValue( theId ) < myMargin; + if ( myToler != theToler ) + myIds.Clear(); + myToler = theToler; } +double ElementsOnSurface::GetTolerance() const +{ return myToler; } -/* - Class : MoreThan - Description : Comparator ">" -*/ -bool MoreThan::IsSatisfy( long theId ) +void ElementsOnSurface::SetUseBoundaries( bool theUse ) { - return myFunctor && myFunctor->GetValue( theId ) > myMargin; + 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; -/* - Class : EqualTo - Description : Comparator "=" -*/ -EqualTo::EqualTo(): - myToler(Precision::Confusion()) -{} + if ( !myMeshModifTracer.GetMesh() ) + return; -bool EqualTo::IsSatisfy( long theId ) + 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 ) { - return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler; + SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator(); + bool isSatisfy = true; + for ( ; aNodeItr->more(); ) + { + SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next(); + if ( !isOnSurface( aNode ) ) + { + isSatisfy = false; + break; + } + } + if ( isSatisfy ) + myIds.Add( theElemPtr->GetID() ); } -void EqualTo::SetTolerance( double theToler ) +bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) { - myToler = theToler; + 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; } /* - Class : LogicalNOT - Description : Logical NOT predicate + ElementsOnShape */ -LogicalNOT::LogicalNOT() -{} -LogicalNOT::~LogicalNOT() -{} +ElementsOnShape::ElementsOnShape() + : //myMesh(0), + myType(SMDSAbs_All), + myToler(Precision::Confusion()), + myAllNodesFlag(false) +{ + myCurShapeType = TopAbs_SHAPE; +} -bool LogicalNOT::IsSatisfy( long theId ) +ElementsOnShape::~ElementsOnShape() { - return myPredicate && !myPredicate->IsSatisfy( theId ); } -void LogicalNOT::SetMesh( SMDS_Mesh* theMesh ) +void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh) { - if ( myPredicate ) - myPredicate->SetMesh( theMesh ); + myMeshModifTracer.SetMesh( theMesh ); + if ( myMeshModifTracer.IsMeshModified()) + SetShape(myShape, myType); } -void LogicalNOT::SetPredicate( PredicatePtr thePred ) +bool ElementsOnShape::IsSatisfy (long theElementId) { - myPredicate = thePred; + return myIds.Contains(theElementId); } -SMDSAbs_ElementType LogicalNOT::GetType() const +SMDSAbs_ElementType ElementsOnShape::GetType() const { - return myPredicate ? myPredicate->GetType() : SMDSAbs_All; + return myType; } +void ElementsOnShape::SetTolerance (const double theToler) +{ + if (myToler != theToler) { + myToler = theToler; + SetShape(myShape, myType); + } +} -/* - Class : LogicalBinary - Description : Base class for binary logical predicate -*/ -LogicalBinary::LogicalBinary() -{} +double ElementsOnShape::GetTolerance() const +{ + return myToler; +} -LogicalBinary::~LogicalBinary() -{} +void ElementsOnShape::SetAllNodes (bool theAllNodes) +{ + if (myAllNodesFlag != theAllNodes) { + myAllNodesFlag = theAllNodes; + SetShape(myShape, myType); + } +} -void LogicalBinary::SetMesh( SMDS_Mesh* theMesh ) +void ElementsOnShape::SetShape (const TopoDS_Shape& theShape, + const SMDSAbs_ElementType theType) { - if ( myPredicate1 ) - myPredicate1->SetMesh( theMesh ); + myType = theType; + myShape = theShape; + myIds.Clear(); - if ( myPredicate2 ) - myPredicate2->SetMesh( theMesh ); + const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh(); + + if ( !myMesh ) return; + + myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType )); + + myShapesMap.Clear(); + addShape(myShape); } -void LogicalBinary::SetPredicate1( PredicatePtr thePredicate ) +void ElementsOnShape::addShape (const TopoDS_Shape& theShape) { - myPredicate1 = thePredicate; + if (theShape.IsNull() || myMeshModifTracer.GetMesh() == 0) + return; + + if (!myShapesMap.Add(theShape)) return; + + myCurShapeType = theShape.ShapeType(); + switch (myCurShapeType) + { + case TopAbs_COMPOUND: + case TopAbs_COMPSOLID: + case TopAbs_SHELL: + case TopAbs_WIRE: + { + TopoDS_Iterator anIt (theShape, Standard_True, Standard_True); + for (; anIt.More(); anIt.Next()) addShape(anIt.Value()); + } + break; + case TopAbs_SOLID: + { + myCurSC.Load(theShape); + process(); + } + break; + case TopAbs_FACE: + { + TopoDS_Face aFace = TopoDS::Face(theShape); + BRepAdaptor_Surface SA (aFace, true); + Standard_Real + u1 = SA.FirstUParameter(), + u2 = SA.LastUParameter(), + v1 = SA.FirstVParameter(), + v2 = SA.LastVParameter(); + Handle(Geom_Surface) surf = BRep_Tool::Surface(aFace); + myCurProjFace.Init(surf, u1,u2, v1,v2); + myCurFace = aFace; + process(); + } + break; + case TopAbs_EDGE: + { + TopoDS_Edge anEdge = TopoDS::Edge(theShape); + Standard_Real u1, u2; + Handle(Geom_Curve) curve = BRep_Tool::Curve(anEdge, u1, u2); + myCurProjEdge.Init(curve, u1, u2); + process(); + } + break; + case TopAbs_VERTEX: + { + TopoDS_Vertex aV = TopoDS::Vertex(theShape); + myCurPnt = BRep_Tool::Pnt(aV); + process(); + } + break; + default: + break; + } } -void LogicalBinary::SetPredicate2( PredicatePtr thePredicate ) +void ElementsOnShape::process() { - myPredicate2 = thePredicate; + const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh(); + if (myShape.IsNull() || myMesh == 0) + return; + + SMDS_ElemIteratorPtr anIter = myMesh->elementsIterator(myType); + while (anIter->more()) + process(anIter->next()); } -SMDSAbs_ElementType LogicalBinary::GetType() const +void ElementsOnShape::process (const SMDS_MeshElement* theElemPtr) { - if ( !myPredicate1 || !myPredicate2 ) - return SMDSAbs_All; + if (myShape.IsNull()) + return; - SMDSAbs_ElementType aType1 = myPredicate1->GetType(); - SMDSAbs_ElementType aType2 = myPredicate2->GetType(); + SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator(); + bool isSatisfy = myAllNodesFlag; - return aType1 == aType2 ? aType1 : SMDSAbs_All; -} + gp_XYZ centerXYZ (0, 0, 0); + while (aNodeItr->more() && (isSatisfy == myAllNodesFlag)) + { + SMESH_TNodeXYZ aPnt ( aNodeItr->next() ); + centerXYZ += aPnt; -/* - Class : LogicalAND - Description : Logical AND -*/ -bool LogicalAND::IsSatisfy( long theId ) -{ - return - myPredicate1 && - myPredicate2 && - myPredicate1->IsSatisfy( theId ) && - myPredicate2->IsSatisfy( theId ); -} + switch (myCurShapeType) + { + case TopAbs_SOLID: + { + myCurSC.Perform(aPnt, myToler); + isSatisfy = (myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON); + } + break; + case TopAbs_FACE: + { + myCurProjFace.Perform(aPnt); + isSatisfy = (myCurProjFace.IsDone() && myCurProjFace.LowerDistance() <= myToler); + if (isSatisfy) + { + // check relatively the face + Quantity_Parameter u, v; + myCurProjFace.LowerDistanceParameters(u, v); + gp_Pnt2d aProjPnt (u, v); + BRepClass_FaceClassifier aClsf (myCurFace, aProjPnt, myToler); + isSatisfy = (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON); + } + } + break; + case TopAbs_EDGE: + { + myCurProjEdge.Perform(aPnt); + isSatisfy = (myCurProjEdge.NbPoints() > 0 && myCurProjEdge.LowerDistance() <= myToler); + } + break; + case TopAbs_VERTEX: + { + isSatisfy = (myCurPnt.Distance(aPnt) <= myToler); + } + break; + default: + { + isSatisfy = false; + } + } + } + if (isSatisfy && myCurShapeType == TopAbs_SOLID) { // Check the center point for volumes MantisBug 0020168 + centerXYZ /= theElemPtr->NbNodes(); + gp_Pnt aCenterPnt (centerXYZ); + myCurSC.Perform(aCenterPnt, myToler); + if ( !(myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON)) + isSatisfy = false; + } -/* - Class : LogicalOR - Description : Logical OR -*/ -bool LogicalOR::IsSatisfy( long theId ) -{ - return - myPredicate1 && - myPredicate2 && - myPredicate1->IsSatisfy( theId ) || - myPredicate2->IsSatisfy( theId ); + if (isSatisfy) + myIds.Add(theElemPtr->GetID()); } +TSequenceOfXYZ::TSequenceOfXYZ() +{} -/* - FILTER -*/ +TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n) +{} -Filter::Filter() +TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t) {} -Filter::~Filter() +TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray) {} -void Filter::SetPredicate( PredicatePtr thePredicate ) +template +TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd) +{} + +TSequenceOfXYZ::~TSequenceOfXYZ() +{} + +TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ) { - myPredicate = thePredicate; + myArray = theSequenceOfXYZ.myArray; + return *this; } -Filter::TIdSequence -Filter::GetElementsId( SMDS_Mesh* theMesh ) +gp_XYZ& TSequenceOfXYZ::operator()(size_type n) { - TIdSequence aSequence; - if ( !theMesh || !myPredicate ) return aSequence; + return myArray[n-1]; +} - myPredicate->SetMesh( theMesh ); +const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const +{ + return myArray[n-1]; +} - SMDSAbs_ElementType aType = myPredicate->GetType(); - switch(aType){ - case SMDSAbs_Edge:{ - SMDS_EdgeIteratorPtr anIter = theMesh->edgesIterator(); - if ( anIter != 0 ) { - while( anIter->more() ) { - const SMDS_MeshElement* anElem = anIter->next(); - long anId = anElem->GetID(); - if ( myPredicate->IsSatisfy( anId ) ) - aSequence.push_back( anId ); - } - } - } - case SMDSAbs_Face:{ - SMDS_FaceIteratorPtr anIter = theMesh->facesIterator(); - if ( anIter != 0 ) { - while( anIter->more() ) { - const SMDS_MeshElement* anElem = anIter->next(); - long anId = anElem->GetID(); - if ( myPredicate->IsSatisfy( anId ) ) - aSequence.push_back( anId ); - } - } - } +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(); +} + +TMeshModifTracer::TMeshModifTracer(): + myMeshModifTime(0), myMesh(0) +{ +} +void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh ) +{ + if ( theMesh != myMesh ) + myMeshModifTime = 0; + myMesh = theMesh; +} +bool TMeshModifTracer::IsMeshModified() +{ + bool modified = false; + if ( myMesh ) + { + modified = ( myMeshModifTime != myMesh->GetMTime() ); + myMeshModifTime = myMesh->GetMTime(); } - return aSequence; + return modified; }