diff --git a/src/Extrema/Extrema_GenExtCC.gxx b/src/Extrema/Extrema_GenExtCC.gxx index 10e497d2e2..86719def54 100644 --- a/src/Extrema/Extrema_GenExtCC.gxx +++ b/src/Extrema/Extrema_GenExtCC.gxx @@ -26,6 +26,7 @@ #include #include #include +#include // Comparator, used in std::sort. static Standard_Boolean comp(const gp_XY& theA, @@ -47,6 +48,57 @@ static Standard_Boolean comp(const gp_XY& theA, return Standard_False; } +static void ChangeIntervals(Handle(TColStd_HArray1OfReal)& theInts, const Standard_Integer theNbInts) +{ + Standard_Integer aNbInts = theInts->Length() - 1; + Standard_Integer aNbAdd = theNbInts - aNbInts; + Handle(TColStd_HArray1OfReal) aNewInts = new TColStd_HArray1OfReal(1, theNbInts + 1); + Standard_Integer aNbLast = theInts->Length(); + Standard_Integer i; + if (aNbInts == 1) + { + aNewInts->SetValue(1, theInts->First()); + aNewInts->SetValue(theNbInts + 1, theInts->Last()); + Standard_Real dt = (theInts->Last() - theInts->First()) / theNbInts; + Standard_Real t = theInts->First() + dt; + for (i = 2; i <= theNbInts; ++i, t += dt) + { + aNewInts->SetValue(i, t); + } + theInts = aNewInts; + return; + } + // + for (i = 1; i <= aNbLast; ++i) + { + aNewInts->SetValue(i, theInts->Value(i)); + } + // + while (aNbAdd > 0) + { + Standard_Real anLIntMax = -1.; + Standard_Integer aMaxInd = -1; + for (i = 1; i < aNbLast; ++i) + { + Standard_Real anL = aNewInts->Value(i + 1) - aNewInts->Value(i); + if (anL > anLIntMax) + { + anLIntMax = anL; + aMaxInd = i; + } + } + + Standard_Real t = (aNewInts->Value(aMaxInd + 1) + aNewInts->Value(aMaxInd)) / 2.; + for (i = aNbLast; i > aMaxInd; --i) + { + aNewInts->SetValue(i + 1, aNewInts->Value(i)); + } + aNbLast++; + aNbAdd--; + aNewInts->SetValue(aMaxInd + 1, t); + } + theInts = aNewInts; +} class Extrema_CCPointsInspector : public NCollection_CellFilter_InspectorXY { public: @@ -111,7 +163,6 @@ static Standard_Real ProjPOnC(const Pnt& theP, if (aD < aDist) aDist = aD; } - aDist = sqrt(aDist); } return aDist; } @@ -229,14 +280,79 @@ void Extrema_GenExtCC::Perform() aNbInter[1] = C2.NbIntervals(aContinuity); } - TColStd_Array1OfReal anIntervals1(1, aNbInter[0] + 1); - TColStd_Array1OfReal anIntervals2(1, aNbInter[1] + 1); - C1.Intervals(anIntervals1, aContinuity); - C2.Intervals(anIntervals2, aContinuity); + Standard_Real anL[2]; + Standard_Integer indmax = -1, indmin = -1; + const Standard_Real mult = 20.; + if (!(Precision::IsInfinite(C1.FirstParameter()) || Precision::IsInfinite(C1.LastParameter()) || + Precision::IsInfinite(C2.FirstParameter()) || Precision::IsInfinite(C2.LastParameter()))) + { + anL[0] = GCPnts_AbscissaPoint::Length(C1); + anL[1] = GCPnts_AbscissaPoint::Length(C2); + if (anL[0] / aNbInter[0] > mult * anL[1] / aNbInter[1]) + { + indmax = 0; + indmin = 1; + } + else if (anL[1] / aNbInter[1] > mult * anL[0] / aNbInter[0]) + { + indmax = 1; + indmin = 0; + } + } + Standard_Integer aNbIntOpt = 0; + if (indmax >= 0) + { + aNbIntOpt = RealToInt(anL[indmax] * aNbInter[indmin] / anL[indmin] / (mult / 4.)) + 1; + if (aNbIntOpt > 100 || aNbIntOpt < aNbInter[indmax]) + { + indmax = -1; + } + else + { + if (aNbIntOpt * aNbInter[indmin] > 100) + { + aNbIntOpt = 100 / aNbInter[indmin]; + if (aNbIntOpt < aNbInter[indmax]) + { + indmax = -1; + } + } + } + } + + Handle(TColStd_HArray1OfReal) anIntervals1 = new TColStd_HArray1OfReal(1, aNbInter[0] + 1); + Handle(TColStd_HArray1OfReal) anIntervals2 = new TColStd_HArray1OfReal(1, aNbInter[1] + 1); + C1.Intervals(anIntervals1->ChangeArray1(), aContinuity); + C2.Intervals(anIntervals2->ChangeArray1(), aContinuity); + if (indmax >= 0) + { + if (indmax == 0) + { + //Change anIntervals1 + ChangeIntervals(anIntervals1, aNbIntOpt); + aNbInter[0] = anIntervals1->Length() - 1; + } + else + { + //Change anIntervals2; + ChangeIntervals(anIntervals2, aNbIntOpt); + aNbInter[1] = anIntervals2->Length() - 1; + } + } + if (C1.IsClosed() && aNbInter[0] == 1) + { + ChangeIntervals(anIntervals1, 3); + aNbInter[0] = anIntervals1->Length() - 1; + } + if (C2.IsClosed() && aNbInter[1] == 1) + { + ChangeIntervals(anIntervals2, 3); + aNbInter[1] = anIntervals2->Length() - 1; + } // Lipchitz constant computation. const Standard_Real aMaxLC = 10000.; - Standard_Real aLC = 9.0; // Default value. + Standard_Real aLC = 100.0; // Default value. const Standard_Real aMaxDer1 = 1.0 / C1.Resolution(1.0); const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0); Standard_Real aMaxDer = Max(aMaxDer1, aMaxDer2) * Sqrt(2.0); @@ -276,6 +392,43 @@ void Extrema_GenExtCC::Perform() } Extrema_GlobOptFuncCCC2 aFunc (C1, C2); + if (aLC < aMaxLC || aMaxDer > aMaxLC) + { + //Estimation of Lipschitz constant by gradient of optimization function + //using sampling in parameter space. + math_Vector aT(1, 2), aG(1, 2); + Standard_Real aF, aMaxG = 0.; + Standard_Real t1, t2, dt1, dt2; + Standard_Integer n1 = 21, n2 = 21, i1, i2; + dt1 = (C1.LastParameter() - C1.FirstParameter()) / (n1 - 1); + dt2 = (C2.LastParameter() - C2.FirstParameter()) / (n2 - 1); + for (i1 = 1, t1 = C1.FirstParameter(); i1 <= n1; ++i1, t1 += dt1) + { + aT(1) = t1; + for (i2 = 1, t2 = C2.FirstParameter(); i2 <= n2; ++i2, t2 += dt2) + { + aT(2) = t2; + aFunc.Values(aT, aF, aG); + Standard_Real aMod = aG(1)*aG(1) + aG(2)*aG(2); + aMaxG = Max(aMaxG, aMod); + } + } + aMaxG = Sqrt(aMaxG); + if (aMaxG > aMaxDer) + { + aLC = Min(aMaxG, aMaxLC); + isConstLockedFlag = Standard_True; + } + if (aMaxG > 100. * aMaxLC) + { + aLC = 100. * aMaxLC; + isConstLockedFlag = Standard_True; + } + else if (aMaxG < 0.1 * aMaxDer) + { + isConstLockedFlag = Standard_True; + } + } math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder, aLC); aFinder.SetLipConstState(isConstLockedFlag); aFinder.SetContinuity(aContinuity == GeomAbs_C2 ? 2 : 1); @@ -286,8 +439,8 @@ void Extrema_GenExtCC::Perform() aFinder.SetFunctionalMinimalValue(0.0); // Best distance cannot be lower than 0.0. // Size computed to have cell index inside of int32 value. - const Standard_Real aCellSize = Max(Max(anIntervals1.Last() - anIntervals1.First(), - anIntervals2.Last() - anIntervals2.First()) + const Standard_Real aCellSize = Max(Max(anIntervals1->Last() - anIntervals1->First(), + anIntervals2->Last() - anIntervals2->First()) * Precision::PConfusion() / (2.0 * Sqrt(2.0)), Precision::PConfusion()); Extrema_CCPointsInspector anInspector(aCellSize); @@ -303,10 +456,10 @@ void Extrema_GenExtCC::Perform() { for(j = 1; j <= aNbInter[1]; j++) { - aFirstBorderInterval(1) = anIntervals1(i); - aFirstBorderInterval(2) = anIntervals2(j); - aSecondBorderInterval(1) = anIntervals1(i + 1); - aSecondBorderInterval(2) = anIntervals2(j + 1); + aFirstBorderInterval(1) = anIntervals1->Value(i); + aFirstBorderInterval(2) = anIntervals2->Value(j); + aSecondBorderInterval(1) = anIntervals1->Value(i + 1); + aSecondBorderInterval(2) = anIntervals2->Value(j + 1); aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval); aFinder.Perform(GetSingleSolutionFlag()); @@ -515,6 +668,7 @@ Standard_Boolean Extrema_GenExtCC::IsDone() const //======================================================================= Standard_Boolean Extrema_GenExtCC::IsParallel() const { + if (!IsDone()) throw StdFail_NotDone(); return myParallel; } @@ -524,8 +678,7 @@ Standard_Boolean Extrema_GenExtCC::IsParallel() const //======================================================================= Standard_Integer Extrema_GenExtCC::NbExt() const { - StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()") - + if (!IsDone()) throw StdFail_NotDone(); return myPoints1.Length(); } @@ -535,8 +688,10 @@ Standard_Integer Extrema_GenExtCC::NbExt() const //======================================================================= Standard_Real Extrema_GenExtCC::SquareDistance(const Standard_Integer N) const { - StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()") - Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()") + if (N < 1 || N > NbExt()) + { + throw Standard_OutOfRange(); + } return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N))); } @@ -549,8 +704,10 @@ void Extrema_GenExtCC::Points(const Standard_Integer N, POnC& P1, POnC& P2) const { - StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::Points()") - Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::Points()") + if (N < 1 || N > NbExt()) + { + throw Standard_OutOfRange(); + } P1.SetValues(myPoints1(N), Tool1::Value(*((Curve1*)myC[0]), myPoints1(N))); P2.SetValues(myPoints2(N), Tool2::Value(*((Curve2*)myC[1]), myPoints2(N))); diff --git a/src/Extrema/Extrema_GlobOptFuncCC.cxx b/src/Extrema/Extrema_GlobOptFuncCC.cxx index 5246d3c05a..1c0f1f07cc 100644 --- a/src/Extrema/Extrema_GlobOptFuncCC.cxx +++ b/src/Extrema/Extrema_GlobOptFuncCC.cxx @@ -45,7 +45,7 @@ static Standard_Boolean _Value(const Adaptor3d_Curve& C1, return Standard_False; } - F = C2.Value(v).Distance(C1.Value(u)); + F = C2.Value(v).SquareDistance(C1.Value(u)); return Standard_True; } @@ -66,7 +66,7 @@ static Standard_Boolean _Value(const Adaptor2d_Curve2d& C1, return Standard_False; } - F = C2.Value(v).Distance(C1.Value(u)); + F = C2.Value(v).SquareDistance(C1.Value(u)); return Standard_True; } @@ -91,13 +91,14 @@ static Standard_Boolean _Gradient(const Adaptor3d_Curve& C1, C1.D1(X(1), C1D0, C1D1); C2.D1(X(2), C2D0, C2D1); - + G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X() - (C2D0.Y() - C1D0.Y()) * C1D1.Y() - (C2D0.Z() - C1D0.Z()) * C1D1.Z(); G(2) = (C2D0.X() - C1D0.X()) * C2D1.X() + (C2D0.Y() - C1D0.Y()) * C2D1.Y() + (C2D0.Z() - C1D0.Z()) * C2D1.Z(); + G *= 2.; return Standard_True; } @@ -123,8 +124,11 @@ static Standard_Boolean _Gradient(const Adaptor2d_Curve2d& C1, G(1) = - (C2D0.X() - C1D0.X()) * C1D1.X() - (C2D0.Y() - C1D0.Y()) * C1D1.Y(); + G(2) = (C2D0.X() - C1D0.X()) * C2D1.X() + (C2D0.Y() - C1D0.Y()) * C2D1.Y(); + G *= 2.; + return Standard_True; } @@ -168,6 +172,7 @@ static Standard_Boolean _Hessian (const Adaptor3d_Curve& C1, + (C2D0.X() - C1D0.X()) * C2D2.X() + (C2D0.Y() - C1D0.Y()) * C2D2.Y() + (C2D0.Z() - C1D0.Z()) * C2D2.Z(); + H *= 2.; return Standard_True; } @@ -206,10 +211,11 @@ static Standard_Boolean _Hessian (const Adaptor2d_Curve2d& C1, + C2D1.Y() * C2D1.Y() + (C2D0.X() - C1D0.X()) * C2D2.X() + (C2D0.Y() - C1D0.Y()) * C2D2.Y(); + H *= 2.; return Standard_True; } -// C0 +//C0 //======================================================================= //function : Extrema_GlobOptFuncCCC0 @@ -405,6 +411,5 @@ Standard_Boolean Extrema_GlobOptFuncCCC2::Values(const math_Vector& X,Standard_R else isHessianComputed = _Hessian(*myC1_2d, *myC2_2d, X, H); - return (Value(X, F) && Gradient(X, G) && isHessianComputed); } diff --git a/src/IntTools/IntTools_FaceFace.cxx b/src/IntTools/IntTools_FaceFace.cxx index 8726f0ad56..6527e5715f 100644 --- a/src/IntTools/IntTools_FaceFace.cxx +++ b/src/IntTools/IntTools_FaceFace.cxx @@ -1188,7 +1188,7 @@ void IntTools_FaceFace::MakeCurve(const Standard_Integer Index, bAvoidLineConstructor, myTol, aSeqOfL, - aReachedTol, + aReachedTol, // obsolete parameter myContext); // aNbSeqOfL=aSeqOfL.Length(); @@ -1196,7 +1196,7 @@ void IntTools_FaceFace::MakeCurve(const Standard_Integer Index, Standard_Real aTolC = 0.; if (bIsDecomposited) { nbiter=aNbSeqOfL; - aTolC = aReachedTol; + aTolC = Precision::Confusion(); } else { nbiter=1; diff --git a/src/IntTools/IntTools_WLineTool.cxx b/src/IntTools/IntTools_WLineTool.cxx index f36c78833b..b9a964f220 100644 --- a/src/IntTools/IntTools_WLineTool.cxx +++ b/src/IntTools/IntTools_WLineTool.cxx @@ -229,203 +229,6 @@ Standard_Boolean IntTools_WLineTool::NotUseSurfacesForApprox(const TopoDS_Face& /////////////////////// DecompositionOfWLine //////////////////////////// -//======================================================================= -//function : CheckTangentZonesExist -//purpose : static subfunction in ComputeTangentZones -//======================================================================= -static -Standard_Boolean CheckTangentZonesExist(const Handle(GeomAdaptor_HSurface)& theSurface1, - const Handle(GeomAdaptor_HSurface)& theSurface2) -{ - if ( ( theSurface1->GetType() != GeomAbs_Torus ) || - ( theSurface2->GetType() != GeomAbs_Torus ) ) - return Standard_False; - - gp_Torus aTor1 = theSurface1->Torus(); - gp_Torus aTor2 = theSurface2->Torus(); - - if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() ) - return Standard_False; - - if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) || - ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) ) - return Standard_False; - - if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) || - ( aTor2.MajorRadius() < aTor2.MinorRadius() ) ) - return Standard_False; - - return Standard_True; -} - - -//======================================================================= -//function : ComputeTangentZones -//purpose : static subfunction in DecompositionOfWLine -//======================================================================= -static -Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1, - const Handle(GeomAdaptor_HSurface)& theSurface2, - const TopoDS_Face& theFace1, - const TopoDS_Face& theFace2, - Handle(TColgp_HArray1OfPnt2d)& theResultOnS1, - Handle(TColgp_HArray1OfPnt2d)& theResultOnS2, - Handle(TColStd_HArray1OfReal)& theResultRadius, - const Handle(IntTools_Context)& aContext) -{ - Standard_Integer aResult = 0; - if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) ) - return aResult; - - - TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2; - TColStd_SequenceOfReal aSeqResultRad; - - gp_Torus aTor1 = theSurface1->Torus(); - gp_Torus aTor2 = theSurface2->Torus(); - - gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() ); - gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() ); - Standard_Integer j = 0; - - for ( j = 0; j < 2; j++ ) { - Standard_Real aCoef = ( j == 0 ) ? -1 : 1; - Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius()); - Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius()); - - gp_Circ aCircle1( anax1, aRadius1 ); - gp_Circ aCircle2( anax2, aRadius2 ); - - // roughly compute radius of tangent zone for perpendicular case - Standard_Real aCriteria = Precision::Confusion() * 0.5; - - Standard_Real aT1 = aCriteria; - Standard_Real aT2 = aCriteria; - if ( j == 0 ) { - // internal tangency - Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius(); - //aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria ); - aT1 = 2. * aR * aCriteria; - aT2 = aT1; - } - else { - // external tangency - Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius(); - Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius(); - Standard_Real aDelta = aRb - aCriteria; - aDelta *= aDelta; - aDelta -= aRm * aRm; - aDelta /= 2. * (aRb - aRm); - aDelta -= 0.5 * (aRb - aRm); - - aT1 = 2. * aRm * (aRm - aDelta); - aT2 = aT1; - } - aCriteria = ( aT1 > aT2) ? aT1 : aT2; - if ( aCriteria > 0 ) - aCriteria = sqrt( aCriteria ); - - if ( aCriteria > 0.5 * aTor1.MinorRadius() ) { - // too big zone -> drop to minimum - aCriteria = Precision::Confusion(); - } - - GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) ); - GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) ); - Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI, - Precision::PConfusion(), Precision::PConfusion()); - - if ( anExtrema.IsDone() ) { - - Standard_Integer i = 0; - for ( i = 1; i <= anExtrema.NbExt(); i++ ) { - if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria ) - continue; - - Extrema_POnCurv P1, P2; - anExtrema.Points( i, P1, P2 ); - - Standard_Boolean bFoundResult = Standard_True; - gp_Pnt2d pr1, pr2; - - Standard_Integer surfit = 0; - for ( surfit = 0; surfit < 2; surfit++ ) { - GeomAPI_ProjectPointOnSurf& aProjector = - (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2); - - gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value(); - aProjector.Perform(aP3d); - - if(!aProjector.IsDone()) - bFoundResult = Standard_False; - else { - if(aProjector.LowerDistance() > aCriteria) { - bFoundResult = Standard_False; - } - else { - Standard_Real foundU = 0, foundV = 0; - aProjector.LowerDistanceParameters(foundU, foundV); - if ( surfit == 0 ) - pr1 = gp_Pnt2d( foundU, foundV ); - else - pr2 = gp_Pnt2d( foundU, foundV ); - } - } - } - if ( bFoundResult ) { - aSeqResultS1.Append( pr1 ); - aSeqResultS2.Append( pr2 ); - aSeqResultRad.Append( aCriteria ); - - // torus is u and v periodic - const Standard_Real twoPI = M_PI + M_PI; - Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()}; - Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()}; - - // iteration on period bounds - for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) { - Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI; - Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI; - - // iteration on surfaces - for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) { - Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp; - Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp; - TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2; - TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2; - - if (fabs(arr1[0] - aBound) < Precision::PConfusion()) { - aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) ); - aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) ); - aSeqResultRad.Append( aCriteria ); - } - if (fabs(arr1[1] - aBound) < Precision::PConfusion()) { - aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) ); - aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) ); - aSeqResultRad.Append( aCriteria ); - } - } - } // - } - } - } - } - aResult = aSeqResultRad.Length(); - - if ( aResult > 0 ) { - theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult ); - theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult ); - theResultRadius = new TColStd_HArray1OfReal( 1, aResult ); - - for ( Standard_Integer i = 1 ; i <= aResult; i++ ) { - theResultOnS1->SetValue( i, aSeqResultS1.Value(i) ); - theResultOnS2->SetValue( i, aSeqResultS2.Value(i) ); - theResultRadius->SetValue( i, aSeqResultRad.Value(i) ); - } - } - return aResult; -} - //======================================================================= //function : IsPointOnBoundary //purpose : static subfunction in DecompositionOfWLine @@ -457,27 +260,6 @@ Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter, return bRet; } -//======================================================================= -//function : IsInsideTanZone -//purpose : Check if point is inside a radial tangent zone. -// static subfunction in DecompositionOfWLine and FindPoint -//======================================================================= -static -Standard_Boolean IsInsideTanZone(const gp_Pnt2d& thePoint, - const gp_Pnt2d& theTanZoneCenter, - const Standard_Real theZoneRadius, - Handle(GeomAdaptor_HSurface) theGASurface) -{ - Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius ); - Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius ); - Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution; - aRadiusSQR *= aRadiusSQR; - if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR ) - return Standard_True; - - return Standard_False; -} - //======================================================================= //function : AdjustByNeighbour //purpose : static subfunction in DecompositionOfWLine @@ -652,72 +434,6 @@ Standard_Boolean FindPoint(const gp_Pnt2d& theFirstPoint, return Standard_False; } -//======================================================================= -//function : FindPoint -//purpose : Find point on the boundary of radial tangent zone -// static subfunction in DecompositionOfWLine -//======================================================================= -static -Standard_Boolean FindPoint(const gp_Pnt2d& theFirstPoint, - const gp_Pnt2d& theLastPoint, - const Standard_Real theUmin, - const Standard_Real theUmax, - const Standard_Real theVmin, - const Standard_Real theVmax, - const gp_Pnt2d& theTanZoneCenter, - const Standard_Real theZoneRadius, - Handle(GeomAdaptor_HSurface) theGASurface, - gp_Pnt2d& theNewPoint) { - theNewPoint = theLastPoint; - - if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) ) - return Standard_False; - - Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius ); - Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius ); - - Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution; - gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) ); - gp_Circ2d aCircle( anAxis, aRadius ); - - // - gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() ); - Standard_Real aLength = aDir.Magnitude(); - if ( aLength <= gp::Resolution() ) - return Standard_False; - gp_Lin2d aLine( theFirstPoint, aDir ); - - // - Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine ); - Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength ); - Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle ); - - Standard_Real aTol = aRadius * 0.001; - aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol; - - Geom2dAPI_InterCurveCurve anIntersector; - anIntersector.Init( aC1, aC2, aTol ); - - if ( anIntersector.NbPoints() == 0 ) - return Standard_False; - - Standard_Boolean aFound = Standard_False; - Standard_Real aMinDist = aLength * aLength; - Standard_Integer i = 0; - for ( i = 1; i <= anIntersector.NbPoints(); i++ ) { - gp_Pnt2d aPInt = anIntersector.Point( i ); - if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) { - if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) && - ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) { - theNewPoint = aPInt; - aFound = Standard_True; - } - } - } - - return aFound; -} - //======================================================================= //function : DecompositionOfWLine //purpose : @@ -732,7 +448,7 @@ Standard_Boolean IntTools_WLineTool:: const Standard_Boolean theAvoidLConstructor, const Standard_Real theTol, IntPatch_SequenceOfLine& theNewLines, - Standard_Real& theReachedTol3d, + Standard_Real& /*theReachedTol3d*/, const Handle(IntTools_Context)& aContext) { Standard_Boolean bRet, bAvoidLineConstructor; @@ -758,13 +474,7 @@ Standard_Boolean IntTools_WLineTool:: TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts); TColStd_Array1OfInteger anArrayOfLineType(1, aNbPnts); TColStd_ListOfInteger aListOfPointIndex; - - Handle(TColgp_HArray1OfPnt2d) aTanZoneS1; - Handle(TColgp_HArray1OfPnt2d) aTanZoneS2; - Handle(TColStd_HArray1OfReal) aTanZoneRadius; - Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2, - aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext); - + // nblines=0; aTol=Precision::Confusion(); @@ -835,24 +545,6 @@ Standard_Boolean IntTools_WLineTool:: bIsCurrentPointOnBoundary = Standard_True; break; } - else { - // check if a point belong to a tangent zone. Begin - Standard_Integer zIt = 0; - for ( zIt = 1; zIt <= aNbZone; zIt++ ) { - gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt); - Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt); - - if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) { - // set boundary flag to split the curve by a tangent zone - bIsPointOnBoundary = Standard_True; - bIsCurrentPointOnBoundary = Standard_True; - if ( theReachedTol3d < aZoneRadius ) { - theReachedTol3d = aZoneRadius; - } - break; - } - } - } }//for(j = 0; j < 2; j++) { if(bIsCurrentPointOnBoundary){ @@ -931,7 +623,7 @@ Standard_Boolean IntTools_WLineTool:: Standard_Integer nbboundaries = 0; Standard_Boolean bIsNearBoundary = Standard_False; - Standard_Integer aZoneIndex = 0; + //Standard_Integer aZoneIndex = 0; Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1 Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1 @@ -981,32 +673,6 @@ Standard_Boolean IntTools_WLineTool:: } } - // check if a point belong to a tangent zone. Begin - for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) { - gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt); - Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt); - - Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast; - const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1); - Standard_Real nU1, nV1; - - if(surfit == 0) - aNeighbourPoint.ParametersOnS1(nU1, nV1); - else - aNeighbourPoint.ParametersOnS2(nU1, nV1); - gp_Pnt2d ap1(nU1, nV1); - gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface ); - - - if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) { - aZoneIndex = zIt; - bIsNearBoundary = Standard_True; - if ( theReachedTol3d < aZoneRadius ) { - theReachedTol3d = aZoneRadius; - } - } - } - // check if a point belong to a tangent zone. End Standard_Boolean bComputeLineEnd = Standard_False; if(nbboundaries == 2) { @@ -1145,20 +811,7 @@ Standard_Boolean IntTools_WLineTool:: gp_Pnt2d ap1(nU1, nV1); gp_Pnt2d ap2; - - if ( aZoneIndex ) { - // exclude point from a tangent zone - anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface ); - gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex); - Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex); - - if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax, - aPZone, aZoneRadius, aGASurface, ap2) ) { - anewpoint = ap2; - found = Standard_True; - } - } - else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) { + if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) { // re-compute point near boundary if shifted on a period ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface ); diff --git a/src/LocOpe/LocOpe_WiresOnShape.cxx b/src/LocOpe/LocOpe_WiresOnShape.cxx index 76fcc8191f..450ef7dea9 100644 --- a/src/LocOpe/LocOpe_WiresOnShape.cxx +++ b/src/LocOpe/LocOpe_WiresOnShape.cxx @@ -1300,6 +1300,10 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge, Standard_Real aTolV[2]; aTolV[0] =BRep_Tool::Tolerance(theVertices[0]); aTolV[1] =BRep_Tool::Tolerance(theVertices[1]); + Standard_Real ext = 16.; // = 4 * 4 - to avoid creating microedges, area around vertices is increased + // up to 4 vertex tolerance. Such approach is usual for other topological + // algorithms, for example, Boolean Operations. + Standard_Real aTolVExt[2] = { ext * aTolV[0] * aTolV[0], ext * aTolV[1] * aTolV[1] }; BRepAdaptor_Curve2d thePCurve(theEdge, theFace); Bnd_Box2d theBox; @@ -1309,6 +1313,9 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge, Standard_Real aFpar, aLpar; const Handle(Geom_Curve)& theCurve = BRep_Tool::Curve(theEdge, thePar[0], thePar[1]); GeomAdaptor_Curve theGAcurve(theCurve, thePar[0], thePar[1]); + Standard_Real aTolV2d[2] = { theGAcurve.Resolution(aTolV[0]), theGAcurve.Resolution(aTolV[1]) }; + aTolV2d[0] = Max(aTolV2d[0], Precision::PConfusion()); + aTolV2d[1] = Max(aTolV2d[1], Precision::PConfusion()); Standard_Real aDistMax = Precision::Confusion() * Precision::Confusion(); TopExp_Explorer Explo(theFace, TopAbs_EDGE); for (; Explo.More(); Explo.Next()) @@ -1335,6 +1342,23 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge, isOverlapped = Standard_True; return; } + // Check extremity distances + Standard_Real dists[4]; + gp_Pnt aP11, aP12, aP21, aP22; + anExtrema.TrimmedSquareDistances(dists[0], dists[1], dists[2], dists[3], + aP11, aP12, aP21, aP22); + for (i = 0; i < 4; ++i) + { + if (i < 2) + j = 0; + else + j = 1; + if (dists[i] < aTolVExt[j] / ext) + { + return; + } + } + for (i = 1; i <= aNbExt; i++) { Standard_Real aDist = anExtrema.SquareDistance(i); @@ -1347,7 +1371,7 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge, Standard_Real anIntPar = aPOnC2.Parameter(); for (j = 0; j < 2; j++) //try to find intersection on an extremity of "theEdge" { - if (Abs(theIntPar - thePar[j]) <= Precision::PConfusion()) + if (Abs(theIntPar - thePar[j]) <= aTolV2d[j]) break; } //intersection found in the middle of the edge @@ -1356,10 +1380,10 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge, gp_Pnt aPoint = aCurve->Value(anIntPar); gp_Pnt aPointInt = theCurve->Value(theIntPar); - if (aPointInt.SquareDistance(thePnt[0]) > aTolV[0] * aTolV[0] && - aPointInt.SquareDistance(thePnt[1]) > aTolV[1] * aTolV[1] && - aPoint.SquareDistance(thePnt[0]) > aTolV[0] * aTolV[0] && - aPoint.SquareDistance(thePnt[1]) > aTolV[1] * aTolV[1]) + if (aPointInt.SquareDistance(thePnt[0]) > aTolVExt[0] && + aPointInt.SquareDistance(thePnt[1]) > aTolVExt[1] && + aPoint.SquareDistance(thePnt[0]) > aTolVExt[0] && + aPoint.SquareDistance(thePnt[1]) > aTolVExt[1]) { SplitPars.Append(theIntPar); if( aDist > aDistMax) diff --git a/tests/bugs/modalg_5/bug23706_10 b/tests/bugs/modalg_5/bug23706_10 index d80c11f0e4..6a5d69e19d 100755 --- a/tests/bugs/modalg_5/bug23706_10 +++ b/tests/bugs/modalg_5/bug23706_10 @@ -13,7 +13,7 @@ set info [extrema r3 r4] if {[regexp "ext_1" $info]} { set dist [lindex [length ext_1] end] - if { $dist > 4.0e-13 } { + if { $dist > 5.0e-11 } { puts "Error: Extrema distance is too big" } } else { diff --git a/tests/bugs/modalg_5/bug23706_11 b/tests/bugs/modalg_5/bug23706_11 index b1afdfe5a0..9df8d5c5f4 100755 --- a/tests/bugs/modalg_5/bug23706_11 +++ b/tests/bugs/modalg_5/bug23706_11 @@ -10,8 +10,11 @@ bsplinecurve r1 2 5 1 3 2 1 3 1 4 1 5 3 2 5 3 1 3 7 3 1 4 8 3 1 4 8 3 1 5 9 3 1 bsplinecurve r2 2 5 2 3 2.5 1 3 1 3.5 1 4 3 -1 2 3 1 1 11 3 1 3 9 3 1 3 9 3 1 3 9 3 1 5 7 3 1 7 4 3 1 set info [extrema r1 r2] -if { [llength $info] != 1 } { - puts "Error : Extrema is wrong" +if {[regexp "ext_1" $info]} { + set dist [lindex [length ext_1] end] + if { $dist > 1.0e-8 } { + puts "Error: Extrema distance is too big" + } } else { - puts "OK: Extrema is valid" + puts "Error: Extrema is not found" } diff --git a/tests/bugs/modalg_6/bug27665 b/tests/bugs/modalg_6/bug27665 index 85fe8294f2..7fa5fc9be7 100644 --- a/tests/bugs/modalg_6/bug27665 +++ b/tests/bugs/modalg_6/bug27665 @@ -5,7 +5,7 @@ puts "" ################################################# # BrepExrtrema_DistShapeShape bad performance on OCCT 6.7.0 ################################################# -cpulimit 100 +cpulimit 500 restore [locate_data_file bug27665_wircmpd.brep] w explode w diff --git a/tests/lowalgos/extcc/bug29858_03 b/tests/lowalgos/extcc/bug29858_03 new file mode 100644 index 0000000000..51702e229d --- /dev/null +++ b/tests/lowalgos/extcc/bug29858_03 @@ -0,0 +1,29 @@ +puts "========" +puts "OCC29858" +puts "========" +puts "" +################################################# +# Regression in GeomAPI_ExtremaCurveCurve +################################################# + +# Read input +restore [locate_data_file bug29858_03_e1.brep] e1 +restore [locate_data_file bug29858_03_e2.brep] e2 + +# Extract geometry from topology +mkcurve c1 e1 +mkcurve c2 e2 + +# Run extrema +set info [extrema c1 c2] + +# Check result +if {[regexp "ext_1" $info]} { + set dist [lindex [length ext_1] end] + if { $dist > 1.0e-10 } { + puts "Error: Extrema distance is too big" + } +} else { + puts "Error: Extrema is not found" +} + diff --git a/tests/lowalgos/extcc/bug32882 b/tests/lowalgos/extcc/bug32882 new file mode 100644 index 0000000000..9e91bbedf6 --- /dev/null +++ b/tests/lowalgos/extcc/bug32882 @@ -0,0 +1,36 @@ +puts "========" +puts "OCC32882: Modeling Data - Extrema curve/curve cannot find all solutions" +puts "========" +puts "" + + +# Read input +restore [locate_data_file bug32882.brep] cc +explode cc + +# Extract geometry from topology +mkcurve c1 cc_1 +mkcurve c2 cc_2 +mkcurve c3 cc_3 + +# Run extrema c1/c3 +set info [extrema c1 c3] + +# Check number of solution +if { [llength $info] != 3 } { + puts "Error: Invalid extrema number in extrema c1-c3 output" +} + +# Check result +checklength ext_1 -l 2.929642751e-14 -eps .01 +checklength ext_2 -l 3.480934286e-14 -eps .01 +checklength ext_3 -l 3.177643716e-14 -eps .01 + +# Run extrema c3/c2 +set info [extrema c3 c2] +# Check number of solutions +if { [llength $info] != 2 } { + puts "Error: Invalid extrema number in extrema c3-c2 output" +} +checklength ext_1 -l 5.684341886e-14 -eps .01 +checklength ext_2 -l 2.929642751e-14 -eps .01