From 2e439615792167de7907f09cc8c897c8a3f7e211 Mon Sep 17 00:00:00 2001 From: eap Date: Tue, 8 Jul 2014 20:07:00 +0400 Subject: [PATCH] IPAL0052448: Too thin viscous layers are constructed --- src/SMESH/SMESH_MesherHelper.cxx | 20 ++++++++--------- src/SMESH/SMESH_MesherHelper.hxx | 3 ++- src/StdMeshers/StdMeshers_Prism_3D.cxx | 8 +++---- src/StdMeshers/StdMeshers_Quadrangle_2D.cxx | 4 ++-- src/StdMeshers/StdMeshers_ViscousLayers.cxx | 24 +++++++++++++-------- 5 files changed, 32 insertions(+), 27 deletions(-) diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index d41ff0578..819b13c4c 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -2732,24 +2732,22 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape ) */ //================================================================================ -double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1, - const TopoDS_Edge & theE2, - const TopoDS_Face & theFace, - gp_Vec* theFaceNormal) +double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1, + const TopoDS_Edge & theE2, + const TopoDS_Face & theFace, + const TopoDS_Vertex & theCommonV, + gp_Vec* theFaceNormal) { double angle = 1e100; try { - TopoDS_Vertex vCommon; - if ( !TopExp::CommonVertex( theE1, theE2, vCommon )) - return angle; double f,l; Handle(Geom_Curve) c1 = BRep_Tool::Curve( theE1, f,l ); Handle(Geom_Curve) c2 = BRep_Tool::Curve( theE2, f,l ); Handle(Geom2d_Curve) c2d1 = BRep_Tool::CurveOnSurface( theE1, theFace, f,l ); Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace ); - double p1 = BRep_Tool::Parameter( vCommon, theE1 ); - double p2 = BRep_Tool::Parameter( vCommon, theE2 ); + double p1 = BRep_Tool::Parameter( theCommonV, theE1 ); + double p2 = BRep_Tool::Parameter( theCommonV, theE2 ); if ( c1.IsNull() || c2.IsNull() ) return angle; 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; int nbLoops = 0; double p1tmp = p1; - while ( vecRef.SquareMagnitude() < std::numeric_limits::min() ) + while ( vecRef.SquareMagnitude() < 1e-25 ) { 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(); surf->D1( uv.X(), uv.Y(), p, du, dv ); vecRef = du ^ dv; diff --git a/src/SMESH/SMESH_MesherHelper.hxx b/src/SMESH/SMESH_MesherHelper.hxx index 91d7d1041..7c8ddbdcf 100644 --- a/src/SMESH/SMESH_MesherHelper.hxx +++ b/src/SMESH/SMESH_MesherHelper.hxx @@ -218,7 +218,8 @@ class SMESH_EXPORT SMESH_MesherHelper static double MaxTolerance( const TopoDS_Shape& shape ); 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 ); diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx index f06261e78..99fdbe90f 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.cxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx @@ -427,10 +427,10 @@ namespace { */ //================================================================================ - 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 ); - } + // 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 ); + // } //================================================================================ /*! diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 6ff221049..f720468f7 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -4070,7 +4070,7 @@ bool StdMeshers_Quadrangle_2D::check() int iPrev = myHelper->WrapIndex( i-1, wire->NbEdges() ); const TopoDS_Edge& e1 = wire->Edge( iPrev ); 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 ) && ( 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 ); 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 )); angleByVertex.Bind( v, angle ); } diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 10006f481..241a21636 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -922,7 +922,7 @@ namespace // get angle between the 2 edges 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 ) { dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] ); @@ -1004,7 +1004,8 @@ namespace while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 ))) iE2 = ( iE2 + 1 ) % nbEdges; 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. ) return true; } @@ -2283,12 +2284,19 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, isOK = false; 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; } - else // hard singularity + else // hard singularity, to call with shiftInside=true ? { const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face ); @@ -2371,12 +2379,10 @@ gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n, } else { - TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]); - TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]); - if ( !v10.IsSame( v01 )) + if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 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 -- 2.39.2