X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;ds=sidebyside;f=src%2FControls%2FSMESH_Controls.cxx;h=f6d48030d74f998beb4000b808c71f62e77b3641;hb=ed90d854b7d3444f34820631d88f62dc34f60ea5;hp=1aef723e89b63098ee06af41b623a08ef660286e;hpb=42c7eb97f9ec27537b638c98ef69dc55f19fa1cd;p=modules%2Fsmesh.git diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index 1aef723e8..f6d48030d 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.cxx @@ -1,4 +1,6 @@ -// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// Copyright (C) 2007-2008 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 @@ -15,37 +17,46 @@ // 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 - +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// #include "SMESH_ControlsDef.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 "SMDS_Mesh.hxx" #include "SMDS_Iterator.hxx" #include "SMDS_MeshElement.hxx" #include "SMDS_MeshNode.hxx" #include "SMDS_VolumeTool.hxx" +#include "SMDS_QuadraticFaceOfNodes.hxx" +#include "SMDS_QuadraticEdge.hxx" + +using namespace std; /* @@ -87,32 +98,69 @@ 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 = 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; } @@ -154,23 +202,41 @@ bool NumericalFunctor::GetPoints(const int theId, } bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem, - TSequenceOfXYZ& theRes ) + TSequenceOfXYZ& theRes ) { theRes.clear(); if ( anElem == 0) return false; + theRes.reserve( anElem->NbNodes() ); + // 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 ){ + SMDS_ElemIteratorPtr anIter; + + if ( anElem->IsQuadratic() ) { + switch ( anElem->GetType() ) { + case SMDSAbs_Edge: + anIter = static_cast + (anElem)->interlacedNodesElemIterator(); + break; + case SMDSAbs_Face: + anIter = static_cast + (anElem)->interlacedNodesElemIterator(); + break; + default: + anIter = anElem->nodesIterator(); + //return false; + } + } + else { + anIter = anElem->nodesIterator(); + } + + 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() ) ); - } } } @@ -189,6 +255,7 @@ void NumericalFunctor::SetPrecision( const long thePrecision ) double NumericalFunctor::GetValue( long theId ) { + myCurrElement = myMesh->FindElement( theId ); TSequenceOfXYZ P; if ( GetPoints( theId, P )) { @@ -204,6 +271,42 @@ double NumericalFunctor::GetValue( long theId ) return 0.; } +//======================================================================= +//function : GetValue +//purpose : +//======================================================================= + +double Volume::GetValue( long theElementId ) +{ + if ( theElementId && myMesh ) { + SMDS_VolumeTool aVolumeTool; + if ( aVolumeTool.Set( myMesh->FindElement( theElementId ))) + return aVolumeTool.GetSize(); + } + return 0; +} + +//======================================================================= +//function : GetBadRate +//purpose : meaningless as it is not quality control functor +//======================================================================= + +double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const +{ + return Value; +} + +//======================================================================= +//function : GetType +//purpose : +//======================================================================= + +SMDSAbs_ElementType Volume::GetType() const +{ + return SMDSAbs_Volume; +} + + /* Class : MinimumAngle Description : Functor for calculation of minimum angle @@ -246,52 +349,99 @@ SMDSAbs_ElementType MinimumAngle::GetType() const */ 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 - - //double aLen[ nbNodes ]; -#ifndef WNT - double aLen [nbNodes]; -#else - double* aLen = (double *)new double[nbNodes]; -#endif - - 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 + 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.; - double aMaxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) ); - static double aCoef = sqrt( 3. ) / 4; - - return aCoef * aMaxLen * aMaxLen / anArea; + return alfa * maxLen * half_perimeter / anArea; } - else - { - double aMinLen = aLen[ 0 ]; - double aMaxLen = aLen[ 0 ]; - - for(int i = 1; i < nbNodes ; i++ ){ - aMinLen = Min( aMinLen, aLen[ i ] ); - aMaxLen = Max( aMaxLen, aLen[ i ] ); - } -#ifdef WNT - delete [] aLen; -#endif - if ( aMinLen <= Precision::Confusion() ) + else if ( nbNodes == 6 ) { // quadratic triangles + // Compute lengths of the sides + 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 aMaxLen / aMinLen; + return alfa * maxLen * half_perimeter / anArea; + } + else if( nbNodes == 4 ) { // quadrangle + // return aspect ratio of the worst triange which can be built + // taking three nodes of the quadrangle + TSequenceOfXYZ triaPnts(3); + // triangle on nodes 1 3 2 + triaPnts(1) = P(1); + triaPnts(2) = P(3); + triaPnts(3) = P(2); + double ar = GetValue( triaPnts ); + // triangle on nodes 1 3 4 + triaPnts(3) = P(4); + ar = Max ( ar, GetValue( triaPnts )); + // triangle on nodes 1 2 4 + triaPnts(2) = P(2); + ar = Max ( ar, GetValue( triaPnts )); + // triangle on nodes 3 2 4 + triaPnts(1) = P(3); + ar = Max ( ar, GetValue( triaPnts )); + + return ar; + } + else { // nbNodes==8 - quadratic quadrangle + // return aspect ratio of the worst triange which can be built + // taking three nodes of the quadrangle + TSequenceOfXYZ triaPnts(3); + // triangle on nodes 1 3 2 + triaPnts(1) = P(1); + triaPnts(2) = P(5); + triaPnts(3) = P(3); + double ar = GetValue( triaPnts ); + // triangle on nodes 1 3 4 + triaPnts(3) = P(7); + ar = Max ( ar, GetValue( triaPnts )); + // triangle on nodes 1 2 4 + triaPnts(2) = P(3); + ar = Max ( ar, GetValue( triaPnts )); + // triangle on nodes 3 2 4 + triaPnts(1) = P(5); + ar = Max ( ar, GetValue( triaPnts )); + + return ar; } } @@ -378,7 +528,18 @@ namespace{ 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 return aQuality; + } + switch(nbNodes){ case 4:{ double aLen[6] = { @@ -411,8 +572,9 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) double aVolume = getVolume(P); //double aVolume = getVolume(aLen); double aHeight = getMaxHeight(aLen); - static double aCoeff = sqrt(6.0)/36.0; - aQuality = aCoeff*aHeight*aSumArea/aVolume; + static double aCoeff = sqrt(2.0)/12.0; + if ( aVolume > DBL_MIN ) + aQuality = aCoeff*aHeight*aSumArea/aVolume; break; } case 5:{ @@ -597,6 +759,21 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) break; } } + 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 = max( aQuality, aspect2D.GetValue( points )); + } + } return aQuality; } @@ -623,7 +800,7 @@ 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 ); @@ -654,7 +831,7 @@ double Warping::ComputeA( const gp_XYZ& thePnt1, N.Normalize(); double H = ( thePnt2 - theG ).Dot( N ); - return asin( fabs( H / L ) ) * 180 / PI; + return asin( fabs( H / L ) ) * 180. / PI; } double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const @@ -678,13 +855,13 @@ SMDSAbs_ElementType Warping::GetType() const 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 ) ) / 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() ) @@ -718,42 +895,46 @@ SMDSAbs_ElementType Taper::GetType() const */ 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; + static double PI2 = 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. / 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 ) ); - return A * 180 / PI; + //BUG SWP12743 + if ( A < Precision::Angular() ) + return 0.; + + return A * 180. / PI; } } @@ -777,21 +958,21 @@ SMDSAbs_ElementType Skew::GetType() const */ 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; + 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); + } + return SumVec.Magnitude() * 0.5; } double Area::GetBadRate( double Value, int /*nbNodes*/ ) const { + // meaningless as it is not a quality control functor return Value; } @@ -807,11 +988,16 @@ SMDSAbs_ElementType Area::GetType() const */ 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 { + // meaningless as it is not quality control functor return Value; } @@ -829,7 +1015,10 @@ double Length2D::GetValue( long theElementId) { TSequenceOfXYZ P; + //cout<<"Length2D::GetValue"<FindElement( theElementId ); @@ -843,7 +1032,11 @@ double Length2D::GetValue( long theElementId) case SMDSAbs_Edge: if (len == 2){ aVal = getDistance( P( 1 ), P( 2 ) ); - break; + 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 @@ -861,6 +1054,22 @@ double Length2D::GetValue( long theElementId) 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="<