X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=6a97d47246cc5928909c9a7037536ebf849d245b;hp=314a0f5c7bcfd13f9e5f79811eef5d769783ff20;hb=1821a9c35f90be93d3d55b1e1db9114902eaf6f8;hpb=8506f58d1068f8c0e4d535851a26ddd7dc85c2ae diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index 314a0f5c7..6a97d4724 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2014 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 @@ -6,7 +6,7 @@ // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -245,8 +245,8 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) TopLoc_Location loc; Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc ); - if ( surface->IsUPeriodic() || surface->IsVPeriodic() || - surface->IsUClosed() || surface->IsVClosed() ) + // if ( surface->IsUPeriodic() || surface->IsVPeriodic() || + // surface->IsUClosed() || surface->IsVClosed() ) { //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface ))) //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface(); @@ -281,7 +281,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) } // look for a degenerated edge - if ( BRep_Tool::Degenerated( edge )) { + if ( SMESH_Algo::isDegenerated( edge )) { myDegenShapeIds.insert( meshDS->ShapeToIndex( edge )); for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() ) myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() )); @@ -478,7 +478,10 @@ bool SMESH_MesherHelper::toCheckPosOnShape(int shapeID ) const void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const { - ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok)); + std::map< int,bool >::iterator sh_ok = + ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok)).first; + if ( !ok ) + sh_ok->second = ok; } //======================================================================= @@ -600,10 +603,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() ) uvOK = ( V == vert.Current() ); if ( !uvOK ) { -#ifdef _DEBUG_ MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID - << " not in face " << GetMeshDS()->ShapeToIndex( F ) ); -#endif + << " not in face " << GetMeshDS()->ShapeToIndex( F ) ); // get UV of a vertex closest to the node double dist = 1e100; gp_Pnt pn = XYZ( n ); @@ -663,9 +664,10 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, const bool force, double distXYZ[4]) const { - int shapeID = n->getshapeId(); + int shapeID = n->getshapeId(); bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() )); - if ( force || toCheckPosOnShape( shapeID ) || infinit ) + bool zero = ( uv.X() == 0. && uv.Y() == 0. ); + if ( force || toCheckPosOnShape( shapeID ) || infinit || zero ) { // check that uv is correct TopLoc_Location loc; @@ -845,7 +847,7 @@ gp_XY SMESH_MesherHelper::GetCenterUV(const gp_XY& uv1, double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E, const SMDS_MeshNode* n, const SMDS_MeshNode* inEdgeNode, - bool* check) + bool* check) const { double param = Precision::Infinite(); @@ -898,8 +900,10 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, const bool force, double distXYZ[4]) const { - int shapeID = n->getshapeId(); - if ( force || toCheckPosOnShape( shapeID )) + int shapeID = n->getshapeId(); + bool infinit = Precision::IsInfinite( u ); + bool zero = ( u == 0. ); + if ( force || toCheckPosOnShape( shapeID ) || infinit || zero ) { TopLoc_Location loc; double f,l; Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l ); @@ -915,12 +919,17 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, { gp_Pnt nodePnt = SMESH_TNodeXYZ( n ); if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() ); - gp_Pnt curvPnt = curve->Value( u ); - double dist = nodePnt.Distance( curvPnt ); - if ( distXYZ ) { - curvPnt.Transform( loc ); - distXYZ[0] = dist; - distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z(); + gp_Pnt curvPnt; + double dist = u; + if ( !infinit ) + { + curvPnt = curve->Value( u ); + dist = nodePnt.Distance( curvPnt ); + if ( distXYZ ) { + curvPnt.Transform( loc ); + distXYZ[0] = dist; + distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z(); + } } if ( dist > tol ) { @@ -944,6 +953,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, } Quantity_Parameter U = projector->LowerDistanceParameter(); u = double( U ); + MESSAGE(" f " << f << " l " << l << " u " << u); curvPnt = curve->Value( u ); dist = nodePnt.Distance( curvPnt ); if ( distXYZ ) { @@ -968,6 +978,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, } if (( u < f-tol || u > l+tol ) && force ) { + MESSAGE("u < f-tol || u > l+tol ; u " << u << " f " << f << " l " << l); // node is on vertex but is set on periodic but trimmed edge (issue 0020890) try { @@ -1172,11 +1183,12 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, gp_XY uvAvg; gp_Pnt P; + bool toCheck = true; if ( !F.IsNull() && !force3d ) { uvAvg = calcTFI (0.5, 0.5, - GetNodeUV(F,n1,n3), GetNodeUV(F,n2,n4), - GetNodeUV(F,n3,n1), GetNodeUV(F,n4,n2), + GetNodeUV(F,n1,n3,&toCheck), GetNodeUV(F,n2,n4,&toCheck), + GetNodeUV(F,n3,n1,&toCheck), GetNodeUV(F,n4,n2,&toCheck), GetNodeUV(F,n12,n3), GetNodeUV(F,n23,n4), GetNodeUV(F,n34,n2), GetNodeUV(F,n41,n2)); TopLoc_Location loc; @@ -1196,10 +1208,10 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, if ( !F.IsNull() ) // force3d { - uvAvg = (GetNodeUV(F,n1,n3) + - GetNodeUV(F,n2,n4) + - GetNodeUV(F,n3,n1) + - GetNodeUV(F,n4,n2)) / 4; + uvAvg = (GetNodeUV(F,n1,n3,&toCheck) + + GetNodeUV(F,n2,n4,&toCheck) + + GetNodeUV(F,n3,n1,&toCheck) + + GetNodeUV(F,n4,n2,&toCheck)) / 4; //CheckNodeUV( F, centralNode, uvAvg, 2*BRep_Tool::Tolerance( F ), /*force=*/true); meshDS->SetNodeOnFace( centralNode, faceID, uvAvg.X(), uvAvg.Y() ); } @@ -1416,8 +1428,15 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, return getMediumNodeOnComposedWire(n1,n2,force3d); } E = TopoDS::Edge(meshDS->IndexToShape( edgeID = pos.first )); - u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]); - u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]); + try { + u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]); + u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]); + } + catch ( Standard_Failure& f ) + { + // issue 22502 / a node is on VERTEX not belonging to E + return getMediumNodeOnComposedWire(n1,n2,force3d); + } } if ( !force3d & uvOK[0] && uvOK[1] ) @@ -2489,9 +2508,16 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) { SMESH_TNodeXYZ nPnt[3]; SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator(); + int iNodeOnFace = 0, iPosDim = SMDS_TOP_VERTEX; for ( int iN = 0; nodesIt->more() && iN < 3; ++iN) // loop on nodes + { nPnt[ iN ] = nodesIt->next(); - + if ( nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition() > iPosDim ) + { + iNodeOnFace = iN; + iPosDim = nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition(); + } + } // compute normal gp_Vec v01( nPnt[0], nPnt[1] ), v02( nPnt[0], nPnt[2] ); if ( v01.SquareMagnitude() > RealSmall() && @@ -2499,7 +2525,7 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) { Ne = v01 ^ v02; if (( normalOK = ( Ne.SquareMagnitude() > RealSmall() ))) - uv = GetNodeUV( theFace, nPnt[0]._node, nPnt[2]._node, &normalOK ); + uv = GetNodeUV( theFace, nPnt[iNodeOnFace]._node, 0, &normalOK ); } } } @@ -2646,6 +2672,76 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape ) return tol; } +//================================================================================ +/*! + * \brief Return an angle between two EDGEs sharing a common VERTEX with reference + * of the FACE normal + * \return double - the angle (between -Pi and Pi), negative if the angle is concave, + * 1e100 in case of failure + * \waring Care about order of the EDGEs and their orientation to be as they are + * within the FACE! Don't pass degenerated EDGEs neither! + */ +//================================================================================ + +double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1, + const TopoDS_Edge & theE2, + const TopoDS_Face & theFace) +{ + 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 ); + if ( c1.IsNull() || c2.IsNull() ) + return angle; + gp_XY uv = c2d1->Value( p1 ).XY(); + gp_Vec du, dv; gp_Pnt p; + surf->D1( uv.X(), uv.Y(), p, du, dv ); + gp_Vec vec1, vec2, vecRef = du ^ dv; + int nbLoops = 0; + double p1tmp = p1; + while ( vecRef.SquareMagnitude() < std::numeric_limits::min() ) + { + double dp = ( l - f ) / 1000.; + 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; + if ( ++nbLoops > 10 ) + { +#ifdef _DEBUG_ + cout << "SMESH_MesherHelper::GetAngle(): Captured in a sigularity" << endl; +#endif + return angle; + } + } + if ( theFace.Orientation() == TopAbs_REVERSED ) + vecRef.Reverse(); + c1->D1( p1, p, vec1 ); + c2->D1( p2, p, vec2 ); + 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 ); + } + catch (...) + { + } + return angle; +} + //================================================================================ /*! * \brief Check if the first and last vertices of an edge are the same @@ -3818,7 +3914,7 @@ namespace { // Structures used by FixQuadraticElements() { // check if the EDGE needs checking const TopoDS_Edge& edge = TopoDS::Edge( edgeIt.Current() ); - if ( BRep_Tool::Degenerated( edge ) ) + if ( SMESH_Algo::isDegenerated( edge ) ) continue; if ( theHelper.IsRealSeam( edge ) && edge.Orientation() == TopAbs_REVERSED ) @@ -4666,4 +4762,3 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, } } } -