PAL10872. Fix aspect ratio for quadrangles

This commit is contained in:
eap 2005-12-29 15:45:10 +00:00
parent e69e13adbf
commit fa1f85dfe5

View File

@ -254,38 +254,65 @@ 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
// See PAL10872
int nbNodes = P.size();
if ( nbNodes != 3 && nbNodes != 4 )
if ( nbNodes < 3 )
return 0;
// Compute lengths of the sides
double aLen[ nbNodes ];
vector< double > aLen (nbNodes);
for ( int i = 0; i < nbNodes - 1; i++ )
aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
// Compute aspect ratio
if ( nbNodes == 3 )
if ( nbNodes == 3 )
{
// 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 = Min( Min( aLen[ 0 ], aLen[ 1 ] ), Min( aLen[ 2 ], aLen[ 3 ] ) );
if ( aMinLen <= Precision::Confusion() )
return 0.;
double aMaxLen = Max( Max( aLen[ 0 ], aLen[ 1 ] ), Max( aLen[ 2 ], aLen[ 3 ] ) );
return aMaxLen / aMinLen;
// 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;
}
}
@ -1116,7 +1143,7 @@ void MultiConnection2D::GetValues(MValues& theValues){
SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
for(; anIter->more(); ){
const SMDS_MeshFace* anElem = anIter->next();
long anElemId = anElem->GetID();
//long anElemId = anElem->GetID();
SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
long aNodeId[3];