IPAL0052448: Too thin viscous layers are constructed

This commit is contained in:
eap 2014-07-08 20:07:00 +04:00
parent bd39d7c3ed
commit 2e43961579
5 changed files with 32 additions and 27 deletions

View File

@ -2735,21 +2735,19 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1, double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
const TopoDS_Edge & theE2, const TopoDS_Edge & theE2,
const TopoDS_Face & theFace, const TopoDS_Face & theFace,
const TopoDS_Vertex & theCommonV,
gp_Vec* theFaceNormal) gp_Vec* theFaceNormal)
{ {
double angle = 1e100; double angle = 1e100;
try try
{ {
TopoDS_Vertex vCommon;
if ( !TopExp::CommonVertex( theE1, theE2, vCommon ))
return angle;
double f,l; double f,l;
Handle(Geom_Curve) c1 = BRep_Tool::Curve( theE1, f,l ); Handle(Geom_Curve) c1 = BRep_Tool::Curve( theE1, f,l );
Handle(Geom_Curve) c2 = BRep_Tool::Curve( theE2, f,l ); Handle(Geom_Curve) c2 = BRep_Tool::Curve( theE2, f,l );
Handle(Geom2d_Curve) c2d1 = BRep_Tool::CurveOnSurface( theE1, theFace, f,l ); Handle(Geom2d_Curve) c2d1 = BRep_Tool::CurveOnSurface( theE1, theFace, f,l );
Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace ); Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace );
double p1 = BRep_Tool::Parameter( vCommon, theE1 ); double p1 = BRep_Tool::Parameter( theCommonV, theE1 );
double p2 = BRep_Tool::Parameter( vCommon, theE2 ); double p2 = BRep_Tool::Parameter( theCommonV, theE2 );
if ( c1.IsNull() || c2.IsNull() ) if ( c1.IsNull() || c2.IsNull() )
return angle; return angle;
gp_XY uv = c2d1->Value( p1 ).XY(); gp_XY uv = c2d1->Value( p1 ).XY();
@ -2758,10 +2756,10 @@ double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
gp_Vec vec1, vec2, vecRef = du ^ dv; gp_Vec vec1, vec2, vecRef = du ^ dv;
int nbLoops = 0; int nbLoops = 0;
double p1tmp = p1; double p1tmp = p1;
while ( vecRef.SquareMagnitude() < std::numeric_limits<double>::min() ) while ( vecRef.SquareMagnitude() < 1e-25 )
{ {
double dp = ( l - f ) / 1000.; double dp = ( l - f ) / 1000.;
p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? +1. : -1.); p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? -1. : +1.);
uv = c2d1->Value( p1tmp ).XY(); uv = c2d1->Value( p1tmp ).XY();
surf->D1( uv.X(), uv.Y(), p, du, dv ); surf->D1( uv.X(), uv.Y(), p, du, dv );
vecRef = du ^ dv; vecRef = du ^ dv;

View File

@ -218,7 +218,8 @@ class SMESH_EXPORT SMESH_MesherHelper
static double MaxTolerance( const TopoDS_Shape& shape ); static double MaxTolerance( const TopoDS_Shape& shape );
static double GetAngle( const TopoDS_Edge & E1, const TopoDS_Edge & E2, static double GetAngle( const TopoDS_Edge & E1, const TopoDS_Edge & E2,
const TopoDS_Face & F, gp_Vec* faceNormal=0); const TopoDS_Face & F, const TopoDS_Vertex & V,
gp_Vec* faceNormal=0);
static bool IsClosedEdge( const TopoDS_Edge& anEdge ); static bool IsClosedEdge( const TopoDS_Edge& anEdge );

View File

@ -427,10 +427,10 @@ namespace {
*/ */
//================================================================================ //================================================================================
double normAngle(const TopoDS_Edge & E1, const TopoDS_Edge & E2, const TopoDS_Face & F) // double normAngle(const TopoDS_Edge & E1, const TopoDS_Edge & E2, const TopoDS_Face & F)
{ // {
return SMESH_MesherHelper::GetAngle( E1, E2, F ) / ( 0.5 * M_PI ); // return SMESH_MesherHelper::GetAngle( E1, E2, F ) / ( 0.5 * M_PI );
} // }
//================================================================================ //================================================================================
/*! /*!

View File

@ -4070,7 +4070,7 @@ bool StdMeshers_Quadrangle_2D::check()
int iPrev = myHelper->WrapIndex( i-1, wire->NbEdges() ); int iPrev = myHelper->WrapIndex( i-1, wire->NbEdges() );
const TopoDS_Edge& e1 = wire->Edge( iPrev ); const TopoDS_Edge& e1 = wire->Edge( iPrev );
const TopoDS_Edge& e2 = wire->Edge( i ); const TopoDS_Edge& e2 = wire->Edge( i );
double angle = myHelper->GetAngle( e1, e2, geomFace ); double angle = myHelper->GetAngle( e1, e2, geomFace, wire->FirstVertex( i ));
if (( maxAngle < angle ) && if (( maxAngle < angle ) &&
( 5.* M_PI/180 < angle && angle < 175.* M_PI/180 )) ( 5.* M_PI/180 < angle && angle < 175.* M_PI/180 ))
{ {
@ -4224,7 +4224,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace,
TopoDS_Vertex v = helper.IthVertex( 0, *edge ); TopoDS_Vertex v = helper.IthVertex( 0, *edge );
if ( !theConsiderMesh || SMESH_Algo::VertexNode( v, helper.GetMeshDS() )) if ( !theConsiderMesh || SMESH_Algo::VertexNode( v, helper.GetMeshDS() ))
{ {
double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace ); double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace, v );
vertexByAngle.insert( make_pair( angle, v )); vertexByAngle.insert( make_pair( angle, v ));
angleByVertex.Bind( v, angle ); angleByVertex.Bind( v, angle );
} }

View File

@ -922,7 +922,7 @@ namespace
// get angle between the 2 edges // get angle between the 2 edges
gp_Vec faceNormal; gp_Vec faceNormal;
double angle = helper.GetAngle( edges[0], edges[1], faceFrw, &faceNormal ); double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
if ( Abs( angle ) < 5 * M_PI/180 ) if ( Abs( angle ) < 5 * M_PI/180 )
{ {
dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] ); dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
@ -1004,7 +1004,8 @@ namespace
while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 ))) while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
iE2 = ( iE2 + 1 ) % nbEdges; iE2 = ( iE2 + 1 ) % nbEdges;
double angle = helper.GetAngle( wires[iW]->Edge( iE1 ), double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
wires[iW]->Edge( iE2 ), F ); wires[iW]->Edge( iE2 ), F,
wires[iW]->FirstVertex( iE2 ));
if ( angle < -5. * M_PI / 180. ) if ( angle < -5. * M_PI / 180. )
return true; return true;
} }
@ -2283,12 +2284,19 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
isOK = false; isOK = false;
Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 ) int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
if ( pointKind < IMPOSSIBLE )
{ {
normal; if ( pointKind != REGULAR && !shiftInside )
{
gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
if ( normShift * normal.XYZ() < 0. )
normal = normShift;
}
isOK = true; isOK = true;
} }
else // hard singularity else // hard singularity, to call with shiftInside=true ?
{ {
const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face ); const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
@ -2371,12 +2379,10 @@ gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n,
} }
else else
{ {
TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]); if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]);
if ( !v10.IsSame( v01 ))
std::swap( ee[0], ee[1] ); std::swap( ee[0], ee[1] );
} }
angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F ); angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
} }
// compute a weighted normal // compute a weighted normal