From c85da25c707edcbe163a551188f8490e2e08ddcd Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 21 May 2014 14:48:00 +0400 Subject: [PATCH] 52426, 22582: EDF 8036 SMESH: ConvertToQuadratic fails with theForce3d off For ConvertToQuadratic: fix a case with a seam edge located not at zero parameter For Viscous layers: fix getting normal at a corner classified as 'Reversed' (+-PI) --- src/SMDS/SMDS_MeshElement.hxx | 2 +- src/SMESH/SMESH_MeshEditor.cxx | 1 + src/SMESH/SMESH_MesherHelper.cxx | 81 +++-- src/SMESH/SMESH_MesherHelper.hxx | 3 +- src/StdMeshers/StdMeshers_ViscousLayers.cxx | 330 +++++++++++++++----- 5 files changed, 311 insertions(+), 106 deletions(-) diff --git a/src/SMDS/SMDS_MeshElement.hxx b/src/SMDS/SMDS_MeshElement.hxx index dc53eaf22..f28ff597a 100644 --- a/src/SMDS/SMDS_MeshElement.hxx +++ b/src/SMDS/SMDS_MeshElement.hxx @@ -154,7 +154,7 @@ public: struct Filter { virtual bool operator()(const SMDS_MeshElement* e) const = 0; - ~Filter() {} + virtual ~Filter() {} }; struct NonNullFilter: public Filter { diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index 0e4aa3de2..8b6afbfe5 100644 --- a/src/SMESH/SMESH_MeshEditor.cxx +++ b/src/SMESH/SMESH_MeshEditor.cxx @@ -8751,6 +8751,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT aHelper.SetIsQuadratic( true ); aHelper.SetIsBiQuadratic( theToBiQuad ); aHelper.SetElementsOnShape(true); + aHelper.ToFixNodeParameters( true ); // convert elements assigned to sub-meshes int nbCheckedElems = 0; diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index 784853234..16f555477 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -242,34 +242,41 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() ) { const TopoDS_Face& face = TopoDS::Face( eF.Current() ); - TopLoc_Location loc; - Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc ); + // TopLoc_Location loc; + // Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc ); // if ( surface->IsUPeriodic() || surface->IsVPeriodic() || // surface->IsUClosed() || surface->IsVClosed() ) { //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface ))) //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface(); - GeomAdaptor_Surface surf( surface ); for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next()) { // look for a seam edge - const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() ); + TopoDS_Edge edge = TopoDS::Edge( exp.Current() ); if ( BRep_Tool::IsClosed( edge, face )) { // initialize myPar1, myPar2 and myParIndex gp_Pnt2d uv1, uv2; BRep_Tool::UVPoints( edge, face, uv1, uv2 ); if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) )) { + double u1 = uv1.Coord(1); + edge.Reverse(); + BRep_Tool::UVPoints( edge, face, uv1, uv2 ); + double u2 = uv1.Coord(1); myParIndex |= U_periodic; - myPar1[0] = surf.FirstUParameter(); - myPar2[0] = surf.LastUParameter(); + myPar1[0] = Min( u1, u2 ); + myPar2[0] = Max( u1, u2 ); } else { + double v1 = uv1.Coord(2); + edge.Reverse(); + BRep_Tool::UVPoints( edge, face, uv1, uv2 ); + double v2 = uv1.Coord(2); myParIndex |= V_periodic; - myPar1[1] = surf.FirstVParameter(); - myPar2[1] = surf.LastVParameter(); + myPar1[1] = Min( v1, v2 ); + myPar2[1] = Max( v1, v2 ); } // store seam shape indices, negative if shape encounters twice int edgeID = meshDS->ShapeToIndex( edge ); @@ -287,13 +294,15 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() )); } } - if ( !myDegenShapeIds.empty() && !myParIndex ) { - if ( surface->IsUPeriodic() || surface->IsUClosed() ) { + if ( !myDegenShapeIds.empty() && !myParIndex ) + { + BRepAdaptor_Surface surf( face, false ); + if ( surf.IsUPeriodic() || surf.IsUClosed() ) { myParIndex |= U_periodic; myPar1[0] = surf.FirstUParameter(); myPar2[0] = surf.LastUParameter(); } - else if ( surface->IsVPeriodic() || surface->IsVClosed() ) { + else if ( surf.IsVPeriodic() || surf.IsVClosed() ) { myParIndex |= V_periodic; myPar1[1] = surf.FirstVParameter(); myPar2[1] = surf.LastVParameter(); @@ -513,8 +522,8 @@ gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& u double p1 = uv1.Coord( i ); double dp1 = Abs( p1-myPar1[i-1]), dp2 = Abs( p1-myPar2[i-1]); if ( myParIndex == i || - dp1 < ( myPar2[i-1] - myPar2[i-1] ) / 100. || - dp2 < ( myPar2[i-1] - myPar2[i-1] ) / 100. ) + dp1 < ( myPar2[i-1] - myPar1[i-1] ) / 100. || + dp2 < ( myPar2[i-1] - myPar1[i-1] ) / 100. ) { double p2 = uv2.Coord( i ); double p1Alt = ( dp1 < dp2 ) ? myPar2[i-1] : myPar1[i-1]; @@ -544,7 +553,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, { // node has position on face const SMDS_FacePosition* fpos = - static_cast(n->GetPosition()); + static_cast( Pos ); uv.SetCoord(fpos->GetUParameter(),fpos->GetVParameter()); if ( check ) uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F )); @@ -555,7 +564,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, // corresponding edge from face, get pcurve for this // edge and retrieve value from this pcurve const SMDS_EdgePosition* epos = - static_cast(n->GetPosition()); + static_cast( Pos ); int edgeID = n->getshapeId(); TopoDS_Edge E = TopoDS::Edge(GetMeshDS()->IndexToShape(edgeID)); double f, l, u = epos->GetUParameter(); @@ -1457,7 +1466,6 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 )); else uv[1].SetCoord( 2, uv[0].Coord( 2 )); } - TopLoc_Location loc; Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc); gp_XY UV = GetMiddleUV( S, uv[0], uv[1] ); @@ -2728,7 +2736,8 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape ) double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1, const TopoDS_Edge & theE2, - const TopoDS_Face & theFace) + const TopoDS_Face & theFace, + gp_Vec* theFaceNormal) { double angle = 1e100; try @@ -2768,16 +2777,33 @@ double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1, } if ( theFace.Orientation() == TopAbs_REVERSED ) vecRef.Reverse(); + if ( theFaceNormal ) *theFaceNormal = vecRef; + c1->D1( p1, p, vec1 ); c2->D1( p2, p, vec2 ); - TopoDS_Face F = theFace; - if ( F.Orientation() == TopAbs_INTERNAL ) - F.Orientation( TopAbs_FORWARD ); + // TopoDS_Face F = theFace; + // if ( F.Orientation() == TopAbs_INTERNAL ) + // F.Orientation( TopAbs_FORWARD ); if ( theE1.Orientation() /*GetSubShapeOri( F, theE1 )*/ == TopAbs_REVERSED ) vec1.Reverse(); if ( theE2.Orientation() /*GetSubShapeOri( F, theE2 )*/ == TopAbs_REVERSED ) vec2.Reverse(); angle = vec1.AngleWithRef( vec2, vecRef ); + + if ( Abs ( angle ) >= 0.99 * M_PI ) + { + BRep_Tool::Range( theE1, f, l ); + p1 += 1e-7 * ( p1-f < l-p1 ? +1. : -1. ); + c1->D1( p1, p, vec1 ); + if ( theE1.Orientation() == TopAbs_REVERSED ) + vec1.Reverse(); + BRep_Tool::Range( theE2, f, l ); + p2 += 1e-7 * ( p2-f < l-p2 ? +1. : -1. ); + c2->D1( p2, p, vec2 ); + if ( theE2.Orientation() == TopAbs_REVERSED ) + vec2.Reverse(); + angle = vec1.AngleWithRef( vec2, vecRef ); + } } catch (...) { @@ -4136,9 +4162,11 @@ namespace { // Structures used by FixQuadraticElements() //BRepClass3d_SolidClassifier solidClassifier( shape ); TIDSortedElemSet checkedVols, movedNodes; - for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID + //for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID + for ( size_t iF = 0; iF < concaveFaces.size(); ++iF ) // loop on concave FACEs { - const TopoDS_Shape& face = faceIt.Current(); + //const TopoDS_Shape& face = faceIt.Current(); + const TopoDS_Shape& face = concaveFaces[ iF ]; SMESHDS_SubMesh* faceSM = meshDS->MeshElements( face ); if ( !faceSM ) continue; @@ -4153,11 +4181,16 @@ namespace { // Structures used by FixQuadraticElements() if ( !vertexSM ) continue; nodeIt = vertexSM->GetNodes(); } + // get ids of sub-shapes of the FACE + set< int > subIDs; + SMESH_subMeshIteratorPtr smIt = + theHelper.GetMesh()->GetSubMesh( face )->getDependsOnIterator(/*includeSelf=*/true); + while ( smIt->more() ) + subIDs.insert( smIt->next()->GetId() ); // find suspicious volumes adjacent to the FACE vector< const SMDS_MeshNode* > nOnFace( 4 ); const SMDS_MeshNode* nInSolid; - //vector< const SMDS_MeshElement* > intersectedFaces; while ( nodeIt->more() ) { const SMDS_MeshNode* n = nodeIt->next(); @@ -4179,7 +4212,7 @@ namespace { // Structures used by FixQuadraticElements() n = *volNode; if ( n->GetPosition()->GetDim() == 3 ) nInSolid = n; - else + else if ( subIDs.count( n->getshapeId() )) nOnFace.push_back( n ); } if ( !nInSolid || nOnFace.size() != nbN - 1 ) diff --git a/src/SMESH/SMESH_MesherHelper.hxx b/src/SMESH/SMESH_MesherHelper.hxx index d67e335b5..0c05c4425 100644 --- a/src/SMESH/SMESH_MesherHelper.hxx +++ b/src/SMESH/SMESH_MesherHelper.hxx @@ -217,7 +217,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); + static double GetAngle( const TopoDS_Edge & E1, const TopoDS_Edge & E2, + const TopoDS_Face & F, gp_Vec* faceNormal=0); static bool IsClosedEdge( const TopoDS_Edge& anEdge ); diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index b2f6855b0..71de50af3 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -83,7 +83,7 @@ #include #include -//#define __myDEBUG +#define __myDEBUG using namespace std; @@ -486,6 +486,11 @@ namespace VISCOUS_3D bool makeLayer(_SolidData& data); bool setEdgeData(_LayerEdge& edge, const set& subIds, SMESH_MesherHelper& helper, _SolidData& data); + gp_XYZ getFaceNormal(const SMDS_MeshNode* n, + const TopoDS_Face& face, + SMESH_MesherHelper& helper, + bool& isOK, + bool shiftInside=false); gp_XYZ getWeigthedNormal( const SMDS_MeshNode* n, std::pair< TGeomID, gp_XYZ > fId2Normal[], const int nbFaces ); @@ -501,6 +506,8 @@ namespace VISCOUS_3D vector< vector<_LayerEdge*> >& edgesByGeom); bool sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); + void limitStepSizeByCurvature( _SolidData& data, + vector< vector<_LayerEdge*> >& edgesByGeom); void limitStepSize( _SolidData& data, const SMDS_MeshElement* face, const double cosin); @@ -778,8 +785,8 @@ namespace if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX && helper.IsClosedEdge( fromE )) { - if ( fabs(u-f) < fabs(u-l )) c->D1( l, p, dv ); - else c->D1( f, p, dv ); + if ( fabs(u-f) < fabs(u-l)) c->D1( l, p, dv ); + else c->D1( f, p, dv ); if ( o == TopAbs_REVERSED ) dv.Reverse(); gp_Vec dir2 = norm ^ dv; @@ -792,51 +799,71 @@ namespace const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok, double* cosin=0) { + TopoDS_Face faceFrw = F; + faceFrw.Orientation( TopAbs_FORWARD ); double f,l; TopLoc_Location loc; - vector< TopoDS_Edge > edges; // sharing a vertex - PShapeIteratorPtr eIt = helper.GetAncestors( fromV, *helper.GetMesh(), TopAbs_EDGE); - while ( eIt->more()) - { - const TopoDS_Edge* e = static_cast( eIt->next() ); - if ( helper.IsSubShape( *e, F ) && !BRep_Tool::Curve( *e, loc,f,l).IsNull() ) - edges.push_back( *e ); - } - gp_XYZ dir(0,0,0); - if ( !( ok = ( edges.size() > 0 ))) return dir; - // get average dir of edges going fromV - gp_XYZ edgeDir; - //if ( edges.size() > 1 ) - for ( size_t i = 0; i < edges.size(); ++i ) - { - edgeDir = getEdgeDir( edges[i], fromV ); - double size2 = edgeDir.SquareModulus(); - if ( size2 > numeric_limits::min() ) - edgeDir /= sqrt( size2 ); - else - ok = false; - dir += edgeDir; - } - gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok ); - try { - if ( edges.size() == 1 ) - dir = fromEdgeDir; - else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees - dir = fromEdgeDir.Normalized() + getFaceDir( F, edges[1], node, helper, ok ).Normalized(); - else if ( dir * fromEdgeDir < 0 ) - dir *= -1; - } - catch ( Standard_Failure ) + TopoDS_Edge edges[2]; // sharing a vertex + int nbEdges = 0; { - ok = false; + TopoDS_Vertex VV[2]; + TopExp_Explorer exp( faceFrw, TopAbs_EDGE ); + for ( ; exp.More() && nbEdges < 2; exp.Next() ) + { + const TopoDS_Edge& e = TopoDS::Edge( exp.Current() ); + if ( SMESH_Algo::isDegenerated( e )) continue; + TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true ); + if ( VV[1].IsSame( fromV )) { + edges[ 0 ] = e; + nbEdges++; + } + else if ( VV[0].IsSame( fromV )) { + edges[ 1 ] = e; + nbEdges++; + } + } } - if ( ok ) + gp_XYZ dir(0,0,0), edgeDir[2]; + if ( nbEdges == 2 ) { - //dir /= edges.size(); + // get dirs of edges going fromV + ok = true; + for ( size_t i = 0; i < nbEdges && ok; ++i ) + { + edgeDir[i] = getEdgeDir( edges[i], fromV ); + double size2 = edgeDir[i].SquareModulus(); + if (( ok = size2 > numeric_limits::min() )) + edgeDir[i] /= sqrt( size2 ); + } + if ( !ok ) return dir; + + // get angle between the 2 edges + gp_Vec faceNormal; + double angle = helper.GetAngle( edges[0], edges[1], faceFrw, &faceNormal ); + if ( Abs( angle ) < 5 * M_PI/180 ) + { + dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] ); + } + else + { + dir = edgeDir[0] + edgeDir[1]; + if ( angle < 0 ) + dir.Reverse(); + } if ( cosin ) { - double angle = gp_Vec( edgeDir ).Angle( dir ); - *cosin = cos( angle ); + double angle = gp_Vec( edgeDir[0] ).Angle( dir ); + *cosin = Cos( angle ); } } + else if ( nbEdges == 1 ) + { + dir = getFaceDir( faceFrw, edges[0], node, helper, ok ); + if ( cosin ) *cosin = 1.; + } + else + { + ok = false; + } + return dir; } //================================================================================ @@ -1557,6 +1584,9 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) // fill data._simplexTestEdges findSimplexTestEdges( data, edgesByGeom ); + // limit data._stepSize depending on surface curvature + limitStepSizeByCurvature( data, edgesByGeom ); + // Put _LayerEdge's into the vector data._edges if ( !sortEdges( data, edgesByGeom )) return false; @@ -1644,6 +1674,62 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize ) } } +//================================================================================ +/*! + * \brief Limit data._stepSize by evaluating curvature of shapes + */ +//================================================================================ + +void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data, + vector< vector<_LayerEdge*> >& edgesByGeom) +{ + const int nbTestPnt = 5; + const double minCurvature = 0.9 / data._hyp->GetTotalThickness(); + + BRepLProp_SLProps surfProp( 2, 1e-6 ); + SMESH_MesherHelper helper( *_mesh ); + + TopExp_Explorer face( data._solid, TopAbs_FACE ); + for ( ; face.More(); face.Next() ) + { + const TopoDS_Face& F = TopoDS::Face( face.Current() ); + BRepAdaptor_Surface surface( F, false ); + surfProp.SetSurface( surface ); + + SMESH_subMesh * sm = _mesh->GetSubMesh( F ); + SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true); + while ( smIt->more() ) + { + sm = smIt->next(); + const vector<_LayerEdge*>& ledges = edgesByGeom[ sm->GetId() ]; + int step = Max( 1, int( ledges.size()) / nbTestPnt ); + for ( size_t i = 0; i < ledges.size(); i += step ) + { + gp_XY uv = helper.GetNodeUV( F, ledges[i]->_nodes[0] ); + surfProp.SetParameters( uv.X(), uv.Y() ); + if ( !surfProp.IsCurvatureDefined() ) + continue; + double surfCurvature = Max( Abs( surfProp.MaxCurvature() ), + Abs( surfProp.MinCurvature() )); + if ( surfCurvature < minCurvature ) + continue; + + gp_Dir minDir, maxDir; + surfProp.CurvatureDirections( maxDir, minDir ); + if ( F.Orientation() == TopAbs_REVERSED ) { + maxDir.Reverse(); minDir.Reverse(); + } + const gp_XYZ& inDir = ledges[i]->_normal; + if ( inDir * maxDir.XYZ() < 0 && + inDir * minDir.XYZ() < 0 ) + continue; + + limitStepSize( data, 0.9 / surfCurvature ); + } + } + } +} + //================================================================================ /*! * Fill data._simplexTestEdges. These _LayerEdge's are used to stop inflation @@ -1666,7 +1752,7 @@ void _ViscousBuilder::findSimplexTestEdges( _SolidData& data, const double minCurvature = 1. / data._hyp->GetTotalThickness(); - for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) + for ( size_t iS = 1; iS < edgesByGeom.size(); ++iS ) { // look for a FACE with layers and w/o _LayerEdge's const vector<_LayerEdge*>& eS = edgesByGeom[iS]; @@ -1935,37 +2021,9 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id )) continue; F = TopoDS::Face( s ); + geomNorm = getFaceNormal( node, F, helper, normOK ); + if ( !normOK ) continue; - gp_XY uv = helper.GetNodeUV( F, node, 0, &normOK ); - Handle(Geom_Surface) surface = BRep_Tool::Surface( F ); - { - gp_Dir normal; - if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 ) - { - geomNorm = normal; - normOK = true; - } - else // hard singularity - { - SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face); - while ( fIt->more() ) - { - const SMDS_MeshElement* f = fIt->next(); - if ( editor.FindShape( f ) == *id ) - { - SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) geomNorm.XYZ(), /*normalized=*/false ); - if ( helper.IsReversedSubMesh( F )) - geomNorm.Reverse(); - break; - } - } - double size2 = geomNorm.SquareMagnitude(); - if ( size2 > numeric_limits::min() ) - geomNorm /= sqrt( size2 ); - else - normOK = false; - } - } if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED ) geomNorm.Reverse(); id2Norm[ totalNbFaces ].first = *id; @@ -1976,6 +2034,22 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, if ( totalNbFaces == 0 ) return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index); + if ( normOK && edge._normal.Modulus() < 1e-3 && totalNbFaces > 1 ) + { + // opposite normals, re-get normals at shifted positions (IPAL 52426) + edge._normal.SetCoord( 0,0,0 ); + for ( int i = 0; i < totalNbFaces; ++i ) + { + const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( id2Norm[i].first )); + geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true ); + if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED ) + geomNorm.Reverse(); + if ( normOK ) + id2Norm[ i ].second = geomNorm.XYZ(); + edge._normal += id2Norm[ i ].second; + } + } + if ( totalNbFaces < 3 ) { //edge._normal /= totalNbFaces; @@ -1992,7 +2066,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, case SMDS_TOP_EDGE: { TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS())); - gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK); + gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK ); double angle = inFaceDir.Angle( edge._normal ); // [0,PI] edge._cosin = cos( angle ); //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl; @@ -2000,8 +2074,8 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, } case SMDS_TOP_VERTEX: { TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS())); - gp_XYZ inFaceDir = getFaceDir( F, V, node, helper, normOK); - double angle = gp_Vec( inFaceDir).Angle( edge._normal ); // [0,PI] + gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK ); + double angle = inFaceDir.Angle( edge._normal ); // [0,PI] edge._cosin = cos( angle ); //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl; break; @@ -2009,7 +2083,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, default: return error(SMESH_Comment("Invalid shape position of node ")<::min() ) @@ -2087,7 +2161,102 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, //================================================================================ /*! - * \brief Return normal at a node weighted with angles taken by FACEs + * \brief Return normal to a FACE at a node + * \param [in] n - node + * \param [in] face - FACE + * \param [in] helper - helper + * \param [out] isOK - true or false + * \param [in] shiftInside - to find normal at a position shifted inside the face + * \return gp_XYZ - normal + */ +//================================================================================ + +gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, + const TopoDS_Face& face, + SMESH_MesherHelper& helper, + bool& isOK, + bool shiftInside) +{ + gp_XY uv; + if ( shiftInside ) + { + // get a shifted position + gp_Pnt p = SMESH_TNodeXYZ( node ); + gp_XYZ shift( 0,0,0 ); + TopoDS_Shape S = helper.GetSubShapeByNode( node, helper.GetMeshDS() ); + switch ( S.ShapeType() ) { + case TopAbs_VERTEX: + { + shift = getFaceDir( face, TopoDS::Vertex( S ), node, helper, isOK ); + break; + } + case TopAbs_EDGE: + { + shift = getFaceDir( face, TopoDS::Edge( S ), node, helper, isOK ); + break; + } + default: + isOK = false; + } + if ( isOK ) + shift.Normalize(); + p.Translate( shift * 1e-5 ); + + TopLoc_Location loc; + GeomAPI_ProjectPointOnSurf& projector = helper.GetProjector( face, loc, 1e-7 ); + + if ( !loc.IsIdentity() ) p.Transform( loc.Transformation().Inverted() ); + + projector.Perform( p ); + if ( !projector.IsDone() || projector.NbPoints() < 1 ) + { + isOK = false; + return p.XYZ(); + } + Quantity_Parameter U,V; + projector.LowerDistanceParameters(U,V); + uv.SetCoord( U,V ); + } + else + { + uv = helper.GetNodeUV( face, node, 0, &isOK ); + } + + gp_Dir normal; + isOK = false; + + Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); + if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 ) + { + normal; + isOK = true; + } + else // hard singularity + { + const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face ); + + SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() ) + { + const SMDS_MeshElement* f = fIt->next(); + if ( f->getshapeId() == faceID ) + { + isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true ); + if ( isOK ) + { + if ( helper.IsReversedSubMesh( face )) + normal.Reverse(); + break; + } + } + } + } + return normal.XYZ(); +} + +//================================================================================ +/*! + * \brief Return a normal at a node weighted with angles taken by FACEs * \param [in] n - the node * \param [in] fId2Normal - FACE ids and normals * \param [in] nbFaces - nb of FACEs meeting at the node @@ -2311,7 +2480,8 @@ gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper ) void _LayerEdge::SetCosin( double cosin ) { _cosin = cosin; - _lenFactor = ( Abs( _cosin ) > 0.1 ) ? 1./sqrt(1-_cosin*_cosin) : 1.0; + cosin = Abs( _cosin ); + _lenFactor = ( /*0.1 < cosin &&*/ cosin < 1-1e-12 ) ? 1./sqrt(1-cosin*cosin) : 1.0; } //================================================================================ @@ -2939,7 +3109,7 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, gp_Vec2d vec1( center, uv1 ); double uLast = vec0.Angle( vec1 ); // -PI - +PI double uMidl = vec0.Angle( vecM ); - if ( uLast * uMidl < 0. ) + if ( uLast * uMidl <= 0. ) uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI; const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() ); @@ -3145,7 +3315,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, if ( edge1->_cosin < 0 ) dirInFace = dir1; else - getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ); + dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ); double angle = dir1.Angle( edge1->_normal ); // [0,PI] edge1->SetCosin( cos( angle )); -- 2.30.2