diff --git a/src/ChFi3d/ChFi3d_Builder_CnCrn.cxx b/src/ChFi3d/ChFi3d_Builder_CnCrn.cxx index 36c4a7ca73..e6bbe9772b 100644 --- a/src/ChFi3d/ChFi3d_Builder_CnCrn.cxx +++ b/src/ChFi3d/ChFi3d_Builder_CnCrn.cxx @@ -419,7 +419,17 @@ static void CurveHermite (const TopOpeBRepDS_DataStructure& DStr, if (ext.NbExt()!=0){ Extrema_POnCurv POnC, POnL; ext.Points(1, POnC, POnL); - param.ChangeValue(nb) =POnC.Parameter(); + if (POnC.Value().Distance(POnL.Value()) < Precision::Confusion()) + param.ChangeValue(nb) =POnC.Parameter(); + else + { + if (!cproj.Value(nb).IsNull()) { + cproj.Value(nb)->D0(cproj.Value(nb)->LastParameter(),p01); + } + else if (!cproj.Value(nb+1).IsNull()) { + cproj.Value(nb+1)->D0(cproj.Value(nb+1)->FirstParameter(),p01); + } + } } } if (!ext.IsDone()||ext.NbExt()==0) { diff --git a/src/Extrema/Extrema.cdl b/src/Extrema/Extrema.cdl index 16cc793a40..b7b6942e90 100644 --- a/src/Extrema/Extrema.cdl +++ b/src/Extrema/Extrema.cdl @@ -92,9 +92,8 @@ is -- generic classes for CURVE-CURVE extremas: ---------------------------------------------- private generic class FuncExtCC, SeqPOnC; - generic class GenExtCC, CCF; + generic class GenExtCC; generic class GenLocateExtCC, CCLocF; - generic class CurveCache; ---------------------------------------------- @@ -126,36 +125,24 @@ is -- Curve-Curve: -- 3d: class ExtCC; - - class CCache instantiates CurveCache from Extrema - (Curve from Adaptor3d, - Pnt from gp, - HArray1OfPnt from TColgp); class ECC instantiates GenExtCC from Extrema (Curve from Adaptor3d, CurveTool from Extrema, Curve from Adaptor3d, CurveTool from Extrema, - CCache from Extrema, HArray1OfPnt from TColgp, POnCurv from Extrema, Pnt from gp, Vec from gp); class LocateExtCC; - - class LCCache instantiates CurveCache from Extrema - (Curve from Adaptor3d, - Pnt from gp, - HArray1OfPnt from TColgp); class ELCC instantiates GenExtCC from Extrema (Curve from Adaptor3d, CurveTool from Extrema, Curve from Adaptor3d, CurveTool from Extrema, - LCCache from Extrema, HArray1OfPnt from TColgp, POnCurv from Extrema, Pnt from gp, @@ -173,17 +160,11 @@ is -- 2d: class ExtCC2d; - class CCache2d instantiates CurveCache from Extrema - (Curve2d from Adaptor2d, - Pnt2d from gp, - HArray1OfPnt2d from TColgp); - class ECC2d instantiates GenExtCC from Extrema (Curve2d from Adaptor2d, Curve2dTool from Extrema, Curve2d from Adaptor2d, Curve2dTool from Extrema, - CCache2d from Extrema, HArray1OfPnt2d from TColgp, POnCurv2d from Extrema, Pnt2d from gp, @@ -191,17 +172,12 @@ is class LocateExtCC2d; - class LCCache2d instantiates CurveCache from Extrema - (Curve2d from Adaptor2d, - Pnt2d from gp, - HArray1OfPnt2d from TColgp); class ELCC2d instantiates GenExtCC from Extrema (Curve2d from Adaptor2d, Curve2dTool from Extrema, Curve2d from Adaptor2d, Curve2dTool from Extrema, - LCCache2d from Extrema, HArray1OfPnt2d from TColgp, POnCurv2d from Extrema, Pnt2d from gp, diff --git a/src/Extrema/Extrema_CurveCache.cdl b/src/Extrema/Extrema_CurveCache.cdl deleted file mode 100644 index 55a5a80a5d..0000000000 --- a/src/Extrema/Extrema_CurveCache.cdl +++ /dev/null @@ -1,114 +0,0 @@ --- Created on: 2008-12-28 --- Created by: Roman Lygin --- Copyright (c) 2008-2014 OPEN CASCADE SAS --- --- This file is part of Open CASCADE Technology software library. --- --- This library is free software; you can redistribute it and/or modify it under --- the terms of the GNU Lesser General Public License version 2.1 as published --- by the Free Software Foundation, with special exception defined in the file --- OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT --- distribution for complete text of the license and disclaimer of any warranty. --- --- Alternatively, this file may be used under the terms of Open CASCADE --- commercial license or contractual agreement. - --- roman.lygin@gmail.com - -generic class CurveCache from Extrema - (Curve as any; -- Adaptor3d_Curve or Adaptor2d_Curve2d - Pnt as any; -- gp_Pnt or gp_Pnt2d - ArrayOfPnt as Transient from Standard) -- TColgp_HArray1OfPnt or TColgp_HArray1OfPnt2d -inherits Transient from Standard - - ---Purpose: - -uses - Array1OfReal from TColStd, - HArray1OfReal from TColStd - -raises NotDone from StdFail - -is - - Create returns mutable CurveCache from Extrema; - - Create (theC: Curve; theUFirst, theULast: Real; theNbSamples: Integer; - theToCalculate: Boolean)returns mutable CurveCache from Extrema; - - Create (theC: Curve; theUFirst, theULast: Real; - IntervalsCN: Array1OfReal from TColStd; - StartIndex, EndIndex: Integer; - Coeff: Real) - returns mutable CurveCache from Extrema; - - SetCurve (me: mutable; theC: Curve; theNbSamples: Integer; - theToCalculate: Boolean); - ---Purpose: Sets the \a theRank-th curve (1 or 2) and its parameters. - -- Caches sample points for further reuse in Perform(). - - SetCurve (me: mutable; theC: Curve; - theUFirst, theULast: Real; theNbSamples: Integer; theToCalculate: Boolean); - ---Purpose: - - SetRange (me: mutable; Uinf, Usup: Real; theToCalculate: Boolean); - ---Purpose: - - CalculatePoints (me: mutable); - ---Purpose: Calculates points along the curve and stores - -- them in internal array for further reuse in Perform(). - - IsValid (me) returns Boolean; - ---C++: inline - ---Purpose: Returns True if the points array has already been calculated - -- for specified curve and range - - Parameters (me) returns HArray1OfReal from TColStd - raises NotDone from StdFail; - ---C++: inline - ---C++: return const & - ---Purpose: Raises StdFail_NotDone if points have not yet been calculated. - - Points (me) returns ArrayOfPnt - raises NotDone from StdFail; - ---C++: inline - ---C++: return const & - ---Purpose: Raises StdFail_NotDone if points have not yet been calculated. - - CurvePtr (me) returns Address from Standard; - ---C++: inline - ---Purpose: - - NbSamples (me) returns Integer; - ---C++: inline - ---Purpose: - - FirstParameter (me) returns Real; - ---C++: inline - ---Purpose: - - LastParameter (me) returns Real; - ---C++: inline - ---Purpose: - - TrimFirstParameter (me) returns Real; - ---C++: inline - ---Purpose: For bounded curves returns FirstParameter(), for unbounded - -1e-10. - - TrimLastParameter (me) returns Real; - ---C++: inline - ---Purpose: For bounded curves returns LastParameter(), for unbounded - 1e-10. - -fields - - myC : Address from Standard; - myFirst : Real; - myLast : Real; - myTrimFirst : Real; - myTrimLast : Real; - myNbSamples : Integer; - myParamArray : HArray1OfReal from TColStd; - myPntArray : ArrayOfPnt; - myIsArrayValid: Boolean; - -end; diff --git a/src/Extrema/Extrema_CurveCache.gxx b/src/Extrema/Extrema_CurveCache.gxx deleted file mode 100644 index 3d4d82792b..0000000000 --- a/src/Extrema/Extrema_CurveCache.gxx +++ /dev/null @@ -1,196 +0,0 @@ -// Created on: 2008-12-28 -// Created by: Roman Lygin -// Copyright (c) 2008-2014 OPEN CASCADE SAS -// -// This file is part of Open CASCADE Technology software library. -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License version 2.1 as published -// by the Free Software Foundation, with special exception defined in the file -// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT -// distribution for complete text of the license and disclaimer of any warranty. -// -// Alternatively, this file may be used under the terms of Open CASCADE -// commercial license or contractual agreement. - -// roman.lygin@gmail.com - -#include - -//======================================================================= -//function : Extrema_CurveCache -//purpose : -//======================================================================= - -Extrema_CurveCache::Extrema_CurveCache(const Curve& theC, - const Standard_Real theUFirst, - const Standard_Real theULast, - const Standard_Integer theNbSamples, - const Standard_Boolean theToCalculate) : - myC (0), myNbSamples (-1), myIsArrayValid (Standard_False) -{ - SetCurve (theC, theUFirst, theULast, theNbSamples, theToCalculate); -} - -//======================================================================= -//function : Extrema_CurveCache -//purpose : with sampling by knots and between them -//======================================================================= - -Extrema_CurveCache::Extrema_CurveCache(const Curve& theC, - const Standard_Real theUFirst, - const Standard_Real theULast, - const TColStd_Array1OfReal& IntervalsCN, - const Standard_Integer StartIndex, - const Standard_Integer EndIndex, - const Standard_Real Coeff) -{ - myC = (Standard_Address)&theC; - myIsArrayValid = Standard_False; - myParamArray.Nullify(); - myPntArray.Nullify(); - - myTrimFirst = myFirst = theUFirst; - myTrimLast = myLast = theULast; - - Standard_Integer Nbp = 3; - if (2 * Coeff < 10000.0) - Nbp = Max((Standard_Integer) (2 * Coeff), Nbp); - myNbSamples = (EndIndex - StartIndex)*Nbp + 1; - - const Standard_Integer aNbSTresh = 10000; - if (myNbSamples > aNbSTresh) - { - Nbp = aNbSTresh / (EndIndex - StartIndex); - myNbSamples = (EndIndex - StartIndex)*Nbp + 1; - } - - //Cache points - myParamArray = new TColStd_HArray1OfReal(1, myNbSamples); - myPntArray = new ArrayOfPnt (1, myNbSamples); - - const Curve& aCurve = *((Curve*)myC); - - Standard_Integer i, j, k = 1; - Standard_Real aPar; - for (i = StartIndex; i < EndIndex; i++) - { - Standard_Real delta = (IntervalsCN(i+1) - IntervalsCN(i)) / Nbp; - for (j = 0; j < Nbp; j++) - { - aPar = IntervalsCN(i) + j*delta; - myParamArray->SetValue(k, aPar); - myPntArray->SetValue(k++, aCurve.Value(aPar)); - } - } - Standard_Real aDelta = (myTrimLast - myTrimFirst) / myNbSamples / 200.; - myParamArray->SetValue(1, myTrimFirst + aDelta); - myPntArray->SetValue(1, aCurve.Value(myTrimFirst + aDelta)); - myParamArray->SetValue(myNbSamples, myTrimLast - aDelta); - myPntArray->SetValue(myNbSamples, aCurve.Value(myTrimLast - aDelta)); - - myIsArrayValid = Standard_True; //cache is available now -} - -//======================================================================= -//function : Extrema_CurveCache -//purpose : -//======================================================================= - -Extrema_CurveCache::Extrema_CurveCache() : myC (0), myNbSamples (-1), - myIsArrayValid (Standard_False) -{ -} - -//======================================================================= -//function : SetCurve -//purpose : -//======================================================================= - -void Extrema_CurveCache::SetCurve (const Curve& theC, - const Standard_Integer theNbSamples, - const Standard_Boolean theToCalculate) -{ - myC = (Standard_Address)&theC; - myNbSamples = theNbSamples; - myIsArrayValid = Standard_False; - myParamArray.Nullify(); - myPntArray.Nullify(); - if (theToCalculate) { - CalculatePoints(); - } -} - -//======================================================================= -//function : SetCurve -//purpose : -//======================================================================= - -void Extrema_CurveCache::SetCurve (const Curve& theC, - const Standard_Real theUFirst, - const Standard_Real theULast, - const Standard_Integer theNbSamples, - const Standard_Boolean theToCalculate) -{ - SetCurve (theC, theNbSamples, Standard_False); //no calculation - SetRange (theUFirst, theULast, theToCalculate); -} - -//======================================================================= -//function : SetRange -//purpose : -//======================================================================= - -void Extrema_CurveCache::SetRange (const Standard_Real theUFirst, - const Standard_Real theULast, - const Standard_Boolean theToCalculate) -{ - //myTrimFirst and myTrimLast are used to compute values on unlimited curves - myTrimFirst = myFirst = theUFirst; - if (Precision::IsInfinite(myTrimFirst)){ - myTrimFirst = -1.0e+10; - } - myTrimLast = myLast = theULast; - if (Precision::IsInfinite(myTrimLast)){ - myTrimLast = 1.0e+10; - } - - myIsArrayValid = Standard_False; - myParamArray.Nullify(); - myPntArray.Nullify(); - if (theToCalculate) { - CalculatePoints(); - } -} - -//======================================================================= -//function : CalculatePoints -//purpose : -//======================================================================= - -void Extrema_CurveCache::CalculatePoints() -{ - if (myIsArrayValid) return; //no need to recalculate if nothing has changed - const Curve& aCurve = *((Curve*)myC); - - // compute myNbSamples point along the [myTrimFirst, myTrimLast] range - - Standard_Real aDelta = myTrimLast - myTrimFirst; - Standard_Real aPar0 = aDelta / myNbSamples / 100.; - aDelta = (aDelta - aPar0) / (myNbSamples - 1); - aPar0 = myTrimFirst + (aPar0/2.); - - //Cache points - - myParamArray = new TColStd_HArray1OfReal(1, myNbSamples); - myPntArray = new ArrayOfPnt (1, myNbSamples); - - Standard_Integer i; - Standard_Real aPar; - for (i = 1, aPar = aPar0; i <= myNbSamples; i++, aPar += aDelta) { - myParamArray->SetValue(i, aPar); - myPntArray->SetValue (i, aCurve.Value (aPar)); - } - - myIsArrayValid = Standard_True; //cache is available now -} diff --git a/src/Extrema/Extrema_CurveCache.lxx b/src/Extrema/Extrema_CurveCache.lxx deleted file mode 100644 index 5672ef0f5f..0000000000 --- a/src/Extrema/Extrema_CurveCache.lxx +++ /dev/null @@ -1,111 +0,0 @@ -// Created on: 2008-12-28 -// Created by: Roman Lygin -// Copyright (c) 2008-2014 OPEN CASCADE SAS -// -// This file is part of Open CASCADE Technology software library. -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License version 2.1 as published -// by the Free Software Foundation, with special exception defined in the file -// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT -// distribution for complete text of the license and disclaimer of any warranty. -// -// Alternatively, this file may be used under the terms of Open CASCADE -// commercial license or contractual agreement. - -// roman.lygin@gmail.com - -#include -#include - -//======================================================================= -//function : IsValid -//purpose : -//======================================================================= - -inline Standard_Boolean Extrema_CurveCache::IsValid() const -{ - return myIsArrayValid; -} - -//======================================================================= -//function : Parameters -//purpose : -//======================================================================= - -inline const Handle(TColStd_HArray1OfReal)& Extrema_CurveCache::Parameters() const -{ - StdFail_NotDone_Raise_if (!myIsArrayValid, "Extrema_CurveCache::Parameters()") - return myParamArray; -} - -//======================================================================= -//function : Points -//purpose : -//======================================================================= - -inline const Handle(ArrayOfPnt)& Extrema_CurveCache::Points() const -{ - StdFail_NotDone_Raise_if (!myIsArrayValid, "Extrema_CurveCache::Points()") - return myPntArray; -} - -//======================================================================= -//function : CurvePtr -//purpose : -//======================================================================= - -inline Standard_Address Extrema_CurveCache::CurvePtr() const -{ - return myC; -} - -//======================================================================= -//function : NbSamples -//purpose : -//======================================================================= - -inline Standard_Integer Extrema_CurveCache::NbSamples() const -{ - return myNbSamples; -} - -//======================================================================= -//function : FirstParameter -//purpose : -//======================================================================= - -inline Standard_Real Extrema_CurveCache::FirstParameter() const -{ - return myFirst; -} - -//======================================================================= -//function : LastParameter -//purpose : -//======================================================================= - -inline Standard_Real Extrema_CurveCache::LastParameter() const -{ - return myLast; -} - -//======================================================================= -//function : TrimFirstParameter -//purpose : -//======================================================================= - -inline Standard_Real Extrema_CurveCache::TrimFirstParameter() const -{ - return myTrimFirst; -} - -//======================================================================= -//function : TrimLastParameter -//purpose : -//======================================================================= - -inline Standard_Real Extrema_CurveCache::TrimLastParameter() const -{ - return myTrimLast; -} diff --git a/src/Extrema/Extrema_ExtCC.cdl b/src/Extrema/Extrema_ExtCC.cdl index 889e62d9aa..2990897d15 100644 --- a/src/Extrema/Extrema_ExtCC.cdl +++ b/src/Extrema/Extrema_ExtCC.cdl @@ -147,7 +147,6 @@ fields myInf: Real [2]; mySup: Real [2]; myTol: Real [2]; - myCacheLists: ListOfTransient from TColStd [2]; -- lists of Handle(Extrema_CCache) P1f: Pnt; P1l: Pnt; P2f: Pnt; diff --git a/src/Extrema/Extrema_ExtCC.cxx b/src/Extrema/Extrema_ExtCC.cxx index cf0cae0d37..5b78646569 100644 --- a/src/Extrema/Extrema_ExtCC.cxx +++ b/src/Extrema/Extrema_ExtCC.cxx @@ -45,7 +45,6 @@ #include #include -#include //======================================================================= //function : Extrema_ExtCC @@ -72,8 +71,9 @@ Extrema_ExtCC::Extrema_ExtCC(const Adaptor3d_Curve& C1, const Standard_Real V1, const Standard_Real V2, const Standard_Real TolC1, - const Standard_Real TolC2) : - myDone (Standard_False) + const Standard_Real TolC2) +: myECC(C1, C2, U1, U2, V1, V2), + myDone (Standard_False) { SetCurve (1, C1, U1, U2); SetCurve (2, C2, V1, V2); @@ -91,8 +91,9 @@ Extrema_ExtCC::Extrema_ExtCC(const Adaptor3d_Curve& C1, Extrema_ExtCC::Extrema_ExtCC(const Adaptor3d_Curve& C1, const Adaptor3d_Curve& C2, const Standard_Real TolC1, - const Standard_Real TolC2) : - myDone (Standard_False) + const Standard_Real TolC2) +: myECC(C1, C2), + myDone (Standard_False) { SetCurve (1, C1, C1.FirstParameter(), C1.LastParameter()); SetCurve (2, C2, C2.FirstParameter(), C2.LastParameter()); @@ -111,9 +112,6 @@ void Extrema_ExtCC::SetCurve (const Standard_Integer theRank, const Adaptor3d_Cu Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_ExtCC::SetCurve()") Standard_Integer anInd = theRank - 1; myC[anInd] = (Standard_Address)&C; - - //clear the previous cache to rebuild it later in Perform() - myCacheLists[anInd].Clear(); } //======================================================================= @@ -163,6 +161,8 @@ void Extrema_ExtCC::SetTolerance (const Standard_Integer theRank, const Standard void Extrema_ExtCC::Perform() { Standard_NullObject_Raise_if (!myC[0] || !myC[1], "Extrema_ExtCC::Perform()") + myECC.SetParams(*((Adaptor3d_Curve*)myC[0]), + *((Adaptor3d_Curve*)myC[1]), myInf[0], mySup[0], myInf[1], mySup[1]); myDone = Standard_False; mypoints.Clear(); mySqDist.Clear(); @@ -194,8 +194,6 @@ void Extrema_ExtCC::Perform() if (Precision::IsInfinite(U12) || Precision::IsInfinite(U22)) mydist22 = RealLast(); else mydist22 = P1l.SquareDistance(P2l); - myECC.SetTolerance (Tol); - //Depending on the types of curves, the algorithm is chosen: //- _ExtElC, when one of the curves is a line and the other is elementary, // or there are two circles; @@ -249,169 +247,12 @@ void Extrema_ExtCC::Perform() Results(CCXtrem, U11, U12, U21, U22); } else { - Standard_Integer i; - Standard_Integer aNbS = 32; //default number of sample points per interval (why 32?) - for (i = 0; i < 2; i++) { - TColStd_ListOfTransient& aCacheList = myCacheLists[i]; - if (aCacheList.IsEmpty()) { - //no caches, build them - Adaptor3d_Curve& aC = *(Adaptor3d_Curve*)myC[i]; - //single interval from myInf[i] to mySup[i] - Handle(Extrema_CCache) aCache = new Extrema_CCache(aC, myInf[i], mySup[i], aNbS, Standard_True); - aCacheList.Append (aCache); - } - } - Handle(Extrema_CCache) aCache1 = Handle(Extrema_CCache)::DownCast (myCacheLists[0].First()); - myECC.SetCurveCache (1, aCache1); - Handle(Extrema_CCache) aCache2 = Handle(Extrema_CCache)::DownCast (myCacheLists[1].First()); - myECC.SetCurveCache (2, aCache2); myECC.Perform(); - Results (myECC, U11, U12, U21, U22); + Results(myECC, U11, U12, U21, U22); } } else { - //common case - use _GenExtCC - //1. check and prepare caches - - Standard_Integer i; - Standard_Integer aNbS[2] = {32, 32}, aNbInter[2] = {1, 1}; - Standard_Real Coeff[2] = {1., 1.}; - Standard_Real rf = 0., rl = 0., LL[2] = {-1., -1.}; - Standard_Boolean KnotSampling[2] = {Standard_False, Standard_False}; - for (i = 0; i < 2; i++) { - TColStd_ListOfTransient& aCacheList = myCacheLists[i]; - if (aCacheList.IsEmpty()) { - //no caches, build them - Adaptor3d_Curve& aC = *(Adaptor3d_Curve*)myC[i]; - - Standard_Real du1 = 0., t = 0.; - gp_Pnt P1, P2; - GeomAbs_CurveType aType = aC.GetType(); - switch (aType) { - case GeomAbs_BezierCurve: - aNbS[i] = aC.NbPoles() * 2; - break; - case GeomAbs_BSplineCurve: - if (aC.Degree() <= 3 && - aC.NbKnots() > 2) - KnotSampling[i] = Standard_True; - - aNbS[i] = aC.NbPoles() * 2; - rf = (Extrema_CurveTool::BSpline(*((Adaptor3d_Curve*)myC[i])))->FirstParameter(); - rl = (Extrema_CurveTool::BSpline(*((Adaptor3d_Curve*)myC[i])))->LastParameter(); - aNbS[i] = (Standard_Integer) ( aNbS[i] * ((mySup[i] - myInf[i]) / (rl - rf)) + 1 ); - case GeomAbs_OtherCurve: - case GeomAbs_Line: - //adjust number of sample points for Lines, B-Splines and Other curves - aNbInter[i] = aC.NbIntervals (GeomAbs_C2); - aNbS[i] = Max(aNbS[i] / aNbInter[i], 3); - LL[i] = 0.; - du1 = (mySup[i] - myInf[i]) / ((aNbS[i]-1)*aNbInter[i]); - P1 = Extrema_CurveTool::Value(*((Adaptor3d_Curve*)myC[i]), myInf[i]); - for(t = myInf[i] + du1; t <= mySup[i]; t += du1) { - P2 = Extrema_CurveTool::Value(*((Adaptor3d_Curve*)myC[i]), t); - LL[i] += P1.Distance(P2); - P1 = P2; - } - LL[i] /= ((aNbS[i]-1)*aNbInter[i]); - break; - default: - break; - } - } - } - - if(LL[0] > 0. && LL[1] > 0.) { - if(LL[0] > 4.*LL[1]) { - Coeff[0] = LL[0]/LL[1]/2.; - if (aNbS[0] * Coeff[0] <= 10000.0) - aNbS[0] = (Standard_Integer) ( aNbS[0] * Coeff[0] ); - } - else if(LL[1] > 4.*LL[0]) { - Coeff[1] = LL[1]/LL[0]/2.; - if (aNbS[1] * Coeff[1] <= 10000.0) - aNbS[1] = (Standard_Integer) (aNbS[1] * Coeff[1] ); - } - } - //modified by NIZNHY-PKV Tue Apr 17 10:01:32 2012f - Standard_Integer aNbSTresh; - // - aNbSTresh=10000; - // - for (i = 0; i < 2; ++i) { - if (aNbS[i]>aNbSTresh) { - aNbS[i]=aNbSTresh; - } - } - //modified by NIZNHY-PKV Tue Apr 17 10:01:34 2012t - for (i = 0; i < 2; i++) { - TColStd_ListOfTransient& aCacheList = myCacheLists[i]; - if (aCacheList.IsEmpty()) { - //no caches, build them - Adaptor3d_Curve& aC = *(Adaptor3d_Curve*)myC[i]; - - if (aC.GetType() >= GeomAbs_BSplineCurve) - { - //can be multiple intervals, one cache per one C2 interval - TColStd_Array1OfReal anArr (1, aNbInter[i] + 1); - aC.Intervals (anArr, GeomAbs_C2); - - if (KnotSampling[i]) - { - Standard_Integer NbIntervCN = aC.NbIntervals(GeomAbs_CN); - TColStd_Array1OfReal IntervalsCN(1, NbIntervCN + 1); - aC.Intervals(IntervalsCN, GeomAbs_CN); - - Standard_Integer j = 1, start_j = 1, k; - Standard_Real NextKnot; - for (k = 1; k <= aNbInter[i]; k++) - { - do - { - NextKnot = IntervalsCN(j+1); - j++; - } - while (NextKnot != anArr(k+1)); - - Handle(Extrema_CCache) aCache = - new Extrema_CCache(aC, anArr(k), anArr(k+1), - IntervalsCN, start_j, j, Coeff[i]); - aCacheList.Append (aCache); - - start_j = j; - } - } - else - { - for (Standard_Integer k = 1; k <= aNbInter[i]; k++) { - Handle(Extrema_CCache) aCache = - new Extrema_CCache(aC, anArr(k), anArr(k+1), aNbS[i], Standard_True); - aCacheList.Append (aCache); - } - } - } - else - { - //single interval from myInf[i] to mySup[i] - Handle(Extrema_CCache) aCache = new Extrema_CCache (aC, - myInf[i], mySup[i], aNbS[i], Standard_True); - aCacheList.Append (aCache); - } - } - } - - //2. process each cache from one list with each cache from the other - TColStd_ListIteratorOfListOfTransient anIt1 (myCacheLists[0]); - for (; anIt1.More(); anIt1.Next()) { - Handle(Extrema_CCache) aCache1 = Handle(Extrema_CCache)::DownCast (anIt1.Value()); - myECC.SetCurveCache (1, aCache1); - TColStd_ListIteratorOfListOfTransient anIt2 (myCacheLists[1]); - for (; anIt2.More(); anIt2.Next()) { - Handle(Extrema_CCache) aCache2 = Handle(Extrema_CCache)::DownCast (anIt2.Value()); - myECC.SetCurveCache (2, aCache2); - myECC.Perform(); - Results(myECC, U11, U12, U21, U22); - } - } + myECC.Perform(); + Results(myECC, U11, U12, U21, U22); } } diff --git a/src/Extrema/Extrema_ExtCC2d.cdl b/src/Extrema/Extrema_ExtCC2d.cdl index f3ffa2e7d4..b156e6bd42 100644 --- a/src/Extrema/Extrema_ExtCC2d.cdl +++ b/src/Extrema/Extrema_ExtCC2d.cdl @@ -27,7 +27,6 @@ uses POnCurv2d from Extrema, SequenceOfReal from TColStd, Curve2d from Adaptor2d, Curve2dTool from Extrema, - CCache2d from Extrema, ECC2d from Extrema @@ -126,9 +125,7 @@ is is static protected; --- modified by NIZHNY-EAP Thu Jan 27 16:53:25 2000 ___BEGIN___ - Results(me: in out;AlgExt: ECC2d from Extrema; C : Curve2d from Adaptor2d; --- modified by NIZHNY-EAP Thu Jan 27 16:53:26 2000 ___END___ + Results(me: in out;AlgExt: ECC2d from Extrema; Ut11, Ut12, Ut21, Ut22: Real; Period1 : Real from Standard = 0.0; Period2 : Real from Standard = 0.0) diff --git a/src/Extrema/Extrema_ExtCC2d.cxx b/src/Extrema/Extrema_ExtCC2d.cxx index c754578a2e..8448867d9e 100644 --- a/src/Extrema/Extrema_ExtCC2d.cxx +++ b/src/Extrema/Extrema_ExtCC2d.cxx @@ -38,21 +38,9 @@ #include -//======================================================================= -//function : IsParallelDot -//purpose : This function returns True if angle between and -// vectors or between and <-theV2> vectors is less than -//AngTol. -//======================================================================= -static Standard_Boolean IsParallelDot( gp_Vec2d theV1, - gp_Vec2d theV2, - Standard_Real AngTol) - { - return Abs(theV1.Dot(theV2)) >= - theV1.Magnitude()*theV2.Magnitude()*cos(AngTol); - } - -Extrema_ExtCC2d::Extrema_ExtCC2d() {} +Extrema_ExtCC2d::Extrema_ExtCC2d() +{ +} Extrema_ExtCC2d::Extrema_ExtCC2d(const Adaptor2d_Curve2d& C1, @@ -100,7 +88,6 @@ void Extrema_ExtCC2d::Perform (const Adaptor2d_Curve2d& C1, { mypoints.Clear(); mySqDist.Clear(); - Standard_Integer NbU = 32, NbV = 32; GeomAbs_CurveType type1 = Extrema_Curve2dTool::GetType(C1), type2 = Extrema_Curve2dTool::GetType(*((Adaptor2d_Curve2d*)myC)); Standard_Real U11, U12, U21, U22, Tol = Min(mytolc1, mytolc2); // Extrema_POnCurv2d P1, P2; @@ -148,11 +135,11 @@ void Extrema_ExtCC2d::Perform (const Adaptor2d_Curve2d& C1, case GeomAbs_BezierCurve: case GeomAbs_OtherCurve: case GeomAbs_BSplineCurve: { - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); Standard_Real Period2 = 0.; if (Extrema_Curve2dTool::IsPeriodic(*((Adaptor2d_Curve2d*)myC))) Period2 = Extrema_Curve2dTool::Period(*((Adaptor2d_Curve2d*)myC)); - Results(Xtrem, C1, U11, U12, U21, U22, 2*M_PI,Period2); + Results(Xtrem, U11, U12, U21, U22, 2*M_PI,Period2); } break; case GeomAbs_Line: { @@ -179,33 +166,33 @@ void Extrema_ExtCC2d::Perform (const Adaptor2d_Curve2d& C1, break; case GeomAbs_Ellipse: { //Extrema_ExtElC2d Xtrem(Extrema_Curve2dTool::Ellipse(C1), Extrema_Curve2dTool::Ellipse(*((Adaptor2d_Curve2d*)myC))); - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); - Results(Xtrem, C1, U11, U12, U21, U22,2*M_PI, 2*M_PI); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); + Results(Xtrem, U11, U12, U21, U22,2*M_PI, 2*M_PI); } break; case GeomAbs_Parabola: { //Extrema_ExtElC2d Xtrem(Extrema_Curve2dTool::Ellipse(C1), Extrema_Curve2dTool::Parabola(*((Adaptor2d_Curve2d*)myC))); - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); - Results(Xtrem, C1, U11, U12, U21, U22, 2*M_PI, 0.); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); + Results(Xtrem, U11, U12, U21, U22, 2*M_PI, 0.); } break; case GeomAbs_Hyperbola: { //Extrema_ExtElC2d Xtrem(Extrema_Curve2dTool::Ellipse(C1), Extrema_Curve2dTool::Hyperbola(*((Adaptor2d_Curve2d*)myC))); - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); - Results(Xtrem, C1, U11, U12, U21, U22, 2*M_PI, 0.); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); + Results(Xtrem, U11, U12, U21, U22, 2*M_PI, 0.); } break; case GeomAbs_BezierCurve: case GeomAbs_OtherCurve: case GeomAbs_BSplineCurve: { - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); Standard_Real Period2 = 0.; if (Extrema_Curve2dTool::IsPeriodic(*((Adaptor2d_Curve2d*)myC))) Period2 = Extrema_Curve2dTool::Period(*((Adaptor2d_Curve2d*)myC)); - Results(Xtrem, C1, U11, U12, U21, U22, 2*M_PI,Period2); + Results(Xtrem, U11, U12, U21, U22, 2*M_PI,Period2); } break; case GeomAbs_Line: { @@ -233,34 +220,34 @@ void Extrema_ExtCC2d::Perform (const Adaptor2d_Curve2d& C1, case GeomAbs_Ellipse: { //inverse = Standard_True; //Extrema_ExtElC2d Xtrem(Extrema_Curve2dTool::Ellipse(*((Adaptor2d_Curve2d*)myC)), Extrema_Curve2dTool::Parabola(C1)); - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); - Results(Xtrem, C1, U11, U12, U21, U22, 0., 2*M_PI); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); + Results(Xtrem, U11, U12, U21, U22, 0., 2*M_PI); } break; case GeomAbs_Parabola: { //Extrema_ExtElC2d Xtrem(Extrema_Curve2dTool::Parabola(C1), Extrema_Curve2dTool::Parabola(*((Adaptor2d_Curve2d*)myC))); - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); - Results(Xtrem, C1, U11, U12, U21, U22, 0., 0.); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); + Results(Xtrem, U11, U12, U21, U22, 0., 0.); } break; case GeomAbs_Hyperbola: { //inverse = Standard_True; //Extrema_ExtElC2d Xtrem(Extrema_Curve2dTool::Hyperbola(*((Adaptor2d_Curve2d*)myC)), Extrema_Curve2dTool::Parabola(C1)); - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); - Results(Xtrem, C1, U11, U12, U21, U22, 0., 0.); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); + Results(Xtrem, U11, U12, U21, U22, 0., 0.); } break; case GeomAbs_BezierCurve: case GeomAbs_OtherCurve: case GeomAbs_BSplineCurve: { - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); Standard_Real Period2 = 0.; if (Extrema_Curve2dTool::IsPeriodic(*((Adaptor2d_Curve2d*)myC))) Period2 = Extrema_Curve2dTool::Period(*((Adaptor2d_Curve2d*)myC)); - Results(Xtrem, C1, U11, U12, U21, U22, 0., Period2); + Results(Xtrem, U11, U12, U21, U22, 0., Period2); } break; case GeomAbs_Line: { @@ -288,33 +275,33 @@ void Extrema_ExtCC2d::Perform (const Adaptor2d_Curve2d& C1, case GeomAbs_Ellipse: { //inverse = Standard_True; //Extrema_ExtElC2d Xtrem(Extrema_Curve2dTool::Ellipse(*((Adaptor2d_Curve2d*)myC)), Extrema_Curve2dTool::Hyperbola(C1)); - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); - Results(Xtrem, C1, U11, U12, U21, U22, 0., 2*M_PI ); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); + Results(Xtrem, U11, U12, U21, U22, 0., 2*M_PI ); } break; case GeomAbs_Parabola: { //Extrema_ExtElC2d Xtrem(Extrema_Curve2dTool::Hyperbola(C1), Extrema_Curve2dTool::Parabola(*((Adaptor2d_Curve2d*)myC))); - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); - Results(Xtrem, C1, U11, U12, U21, U22, 0., 0.); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); + Results(Xtrem, U11, U12, U21, U22, 0., 0.); } break; case GeomAbs_Hyperbola: { //Extrema_ExtElC2d Xtrem(Extrema_Curve2dTool::Hyperbola(C1), Extrema_Curve2dTool::Hyperbola(*((Adaptor2d_Curve2d*)myC))); - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); - Results(Xtrem, C1, U11, U12, U21, U22, 0., 0.); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); + Results(Xtrem, U11, U12, U21, U22, 0., 0.); } break; case GeomAbs_OtherCurve: case GeomAbs_BezierCurve: case GeomAbs_BSplineCurve: { - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC) - , NbU, NbV, mytolc1, mytolc2); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); Standard_Real Period2 = 0.; if (Extrema_Curve2dTool::IsPeriodic(*((Adaptor2d_Curve2d*)myC))) Period2 = Extrema_Curve2dTool::Period(*((Adaptor2d_Curve2d*)myC)); - Results(Xtrem, C1, U11, U12, U21, U22, 0., Period2); + Results(Xtrem, U11, U12, U21, U22, 0., Period2); } break; case GeomAbs_Line: { @@ -333,13 +320,13 @@ void Extrema_ExtCC2d::Perform (const Adaptor2d_Curve2d& C1, case GeomAbs_BezierCurve: case GeomAbs_OtherCurve: case GeomAbs_BSplineCurve: { - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); Standard_Real Period1 = 0.; if (Extrema_Curve2dTool::IsPeriodic(C1)) Period1 = Extrema_Curve2dTool::Period(C1); Standard_Real Period2 = 0.; if (Extrema_Curve2dTool::IsPeriodic(*((Adaptor2d_Curve2d*)myC))) Period2 = Extrema_Curve2dTool::Period(*((Adaptor2d_Curve2d*)myC)); - Results(Xtrem, C1, U11, U12, U21, U22, Period1, Period2); + Results(Xtrem, U11, U12, U21, U22, Period1, Period2); } break; @@ -372,11 +359,11 @@ void Extrema_ExtCC2d::Perform (const Adaptor2d_Curve2d& C1, case GeomAbs_BezierCurve: case GeomAbs_OtherCurve: case GeomAbs_BSplineCurve: { - Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC), - NbU, NbV, mytolc1, mytolc2); + Extrema_ECC2d Xtrem(C1, *((Adaptor2d_Curve2d*)myC)); + Xtrem.Perform(); Standard_Real Period2 = 0.; if (Extrema_Curve2dTool::IsPeriodic(*((Adaptor2d_Curve2d*)myC))) Period2 = Extrema_Curve2dTool::Period(*((Adaptor2d_Curve2d*)myC)); - Results(Xtrem, C1, U11, U12, U21, U22, 0., Period2); + Results(Xtrem, U11, U12, U21, U22, 0., Period2); } break; case GeomAbs_Line: { @@ -511,54 +498,47 @@ void Extrema_ExtCC2d::Results(const Extrema_ExtElC2d& AlgExt, void Extrema_ExtCC2d::Results(const Extrema_ECC2d& AlgExt, -// modified by NIZHNY-EAP Wed Feb 23 14:51:24 2000 ___BEGIN___ - const Adaptor2d_Curve2d& C1, -// modified by NIZHNY-EAP Wed Feb 23 14:51:26 2000 ___END___ - const Standard_Real Ut11, - const Standard_Real Ut12, - const Standard_Real Ut21, - const Standard_Real Ut22, - const Standard_Real Period1, - const Standard_Real Period2) + const Standard_Real Ut11, + const Standard_Real Ut12, + const Standard_Real Ut21, + const Standard_Real Ut22, + const Standard_Real Period1, + const Standard_Real Period2) { Standard_Integer i, NbExt; Standard_Real Val, U, U2; Extrema_POnCurv2d P1, P2; - - myDone = AlgExt.IsDone(); - if (myDone) { - if (!myIsPar) { - NbExt = AlgExt.NbExt(); - for (i = 1; i <= NbExt; i++) { - // Verification de la validite des parametres pour le cas trimme: - AlgExt.Points(i, P1, P2); - U = P1.Parameter(); - if (Period1 != 0.0) U = ElCLib::InPeriod(U,Ut11,Ut11+Period1); - U2 = P2.Parameter(); - if (Period2 != 0.0) U2 = ElCLib::InPeriod(U2,Ut21,Ut21+Period2); - if ((U >= Ut11 - Precision::PConfusion()) && - (U <= Ut12 + Precision::PConfusion()) && - (U2 >= Ut21 - Precision::PConfusion()) && - (U2 <= Ut22 + Precision::PConfusion())) { -// modified by NIZHNY-EAP Thu Jan 27 16:40:55 2000 ___BEGIN___ - // to be sure that it's a real extrema - gp_Pnt2d p; - gp_Vec2d v1, v2; - Extrema_Curve2dTool::D1(C1,U,p, v1); - Extrema_Curve2dTool::D1(*((Adaptor2d_Curve2d*)myC),U2,p, v2); - if (IsParallelDot(v1, v2, Precision::Angular())) - { - mynbext++; - Val = AlgExt.SquareDistance(i); - P1.SetValues(U, P1.Value()); - P2.SetValues(U2, P2.Value()); - mySqDist.Append(Val); - mypoints.Append(P1); - mypoints.Append(P2); - } -// modified by NIZHNY-EAP Thu Jan 27 16:41:00 2000 ___END___ - } + myDone = AlgExt.IsDone(); + if (myDone) + { + if (!myIsPar) + { + NbExt = AlgExt.NbExt(); + for (i = 1; i <= NbExt; i++) + { + // Verification de la validite des parametres pour le cas trimme: + AlgExt.Points(i, P1, P2); + U = P1.Parameter(); + if (Period1 != 0.0) + U = ElCLib::InPeriod(U,Ut11,Ut11+Period1); + U2 = P2.Parameter(); + if (Period2 != 0.0) + U2 = ElCLib::InPeriod(U2,Ut21,Ut21+Period2); + + if ((U >= Ut11 - Precision::PConfusion()) && + (U <= Ut12 + Precision::PConfusion()) && + (U2 >= Ut21 - Precision::PConfusion()) && + (U2 <= Ut22 + Precision::PConfusion())) + { + mynbext++; + Val = AlgExt.SquareDistance(i); + P1.SetValues(U, P1.Value()); + P2.SetValues(U2, P2.Value()); + mySqDist.Append(Val); + mypoints.Append(P1); + mypoints.Append(P2); + } } } diff --git a/src/Extrema/Extrema_GenExtCC.cdl b/src/Extrema/Extrema_GenExtCC.cdl index 0c724a52cb..64b9fb52e8 100644 --- a/src/Extrema/Extrema_GenExtCC.cdl +++ b/src/Extrema/Extrema_GenExtCC.cdl @@ -19,7 +19,6 @@ generic class GenExtCC from Extrema Tool1 as any; -- as ToolCurve(Curve1) Curve2 as any; Tool2 as any; -- as ToolCurve(Curve2) - Cache as Transient from Standard; -- CurveCache from Extrema ArrayOfPnt as Transient from Standard; -- as returned by Extrema_CurveCache::Points() POnC as any; Pnt as any; @@ -27,14 +26,12 @@ generic class GenExtCC from Extrema ---Purpose: It calculates all the distance between two curves. -- These distances can be maximum or minimum. - +uses SequenceOfReal from TColStd, + Vector from math + raises InfiniteSolutions from StdFail, NotDone from StdFail, OutOfRange from Standard - -private class CCF instantiates FuncExtCC from Extrema (Curve1, Tool1, - Curve2, Tool2, - POnC, Pnt, Vec); is @@ -43,30 +40,20 @@ is -- between Uinf and Usup for C1 and between Vinf and Vsup -- for C2. - Create (C1: Curve1; C2: Curve2; NbU,NbV: Integer; TolU,TolV: Real) returns GenExtCC; + Create (C1: Curve1; C2: Curve2) returns GenExtCC; ---Purpose: It calculates all the distances. -- The function F(u,v)=distance(C1(u),C2(v)) has an - -- extremum when gradient(f)=0. The algorithm searchs - -- all the zeros. - -- NbU,NbV are used to locate the close points to - -- find the zeros. They must be great enough such that - -- if there is N extrema, there will be N extrema - -- between the segment and the grid. - -- TolU and TolV are used to determinethe conditions - -- to stop the iterations; at the iteration number n: - -- (Un - Un-1) < TolU and (Vn - Vn-1) < TolV . + -- extremum when gradient(f)=0. The algorithm uses + -- Evtushenko's global optimization solver.. - Create (C1: Curve1; C2: Curve2; Uinf, Usup, Vinf, Vsup: Real; - NbU,NbV: Integer; TolU,TolV: Real) returns GenExtCC; + Create (C1: Curve1; C2: Curve2; Uinf, Usup, Vinf, Vsup: Real) returns GenExtCC; ---Purpose: Calculates all the distances as above -- between Uinf and Usup for C1 and between Vinf and Vsup -- for C2. - SetCurveCache (me: in out; theRank: Integer; theCache: Cache); - ---Purpose: - - SetTolerance (me: in out; Tol: Real); - ---Purpose: + SetParams(me: in out; C1: Curve1; C2: Curve2; Uinf, Usup, Vinf, Vsup: Real) + ---Purpose: Set params in case of empty constructor is usage. + is static; Perform(me: in out) is static; ---Purpose: Performs calculations. @@ -106,8 +93,12 @@ is is static; fields - myF : CCF from Extrema; + myLowBorder : Vector from math; + myUppBorder : Vector from math; + mySolCount : Integer from Standard; + myPoints1 : SequenceOfReal from TColStd; + myPoints2 : SequenceOfReal from TColStd; + myC : Address from Standard [2]; myDone : Boolean; - myCache: Cache [2]; end GenExtCC; diff --git a/src/Extrema/Extrema_GenExtCC.gxx b/src/Extrema/Extrema_GenExtCC.gxx index 15932d9ccc..90e4e8a7d8 100644 --- a/src/Extrema/Extrema_GenExtCC.gxx +++ b/src/Extrema/Extrema_GenExtCC.gxx @@ -14,282 +14,171 @@ // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. -#include -#include -#include -#include -#include +#include +#include #include #include +#include +#include #include -Extrema_GenExtCC::Extrema_GenExtCC () : myDone (Standard_False) +Extrema_GenExtCC::Extrema_GenExtCC() +: myLowBorder(1,2), + myUppBorder(1,2), + myDone(Standard_False) { } Extrema_GenExtCC::Extrema_GenExtCC (const Curve1& C1, - const Curve2& C2, - const Standard_Integer NbU, - const Standard_Integer NbV, - const Standard_Real TolU, - const Standard_Real TolV) : myF (C1,C2, Min(TolU, TolV)), myDone (Standard_False) + const Curve2& C2) +: myLowBorder(1,2), + myUppBorder(1,2), + myDone(Standard_False) { - SetCurveCache (1, new Cache (C1, C1.FirstParameter(), C1.LastParameter(), NbU, Standard_True)); - SetCurveCache (2, new Cache (C2, C2.FirstParameter(), C2.LastParameter(), NbV, Standard_True)); - Perform(); + myC[0] = (Standard_Address)&C1; + myC[1] = (Standard_Address)&C2; + myLowBorder(1) = C1.FirstParameter(); + myLowBorder(2) = C2.FirstParameter(); + myUppBorder(1) = C1.LastParameter(); + myUppBorder(2) = C2.LastParameter(); } Extrema_GenExtCC::Extrema_GenExtCC (const Curve1& C1, - const Curve2& C2, - const Standard_Real Uinf, - const Standard_Real Usup, - const Standard_Real Vinf, - const Standard_Real Vsup, - const Standard_Integer NbU, - const Standard_Integer NbV, - const Standard_Real /*TolU*/, - const Standard_Real /*TolV*/) : myF (C1,C2), myDone (Standard_False) + const Curve2& C2, + const Standard_Real Uinf, + const Standard_Real Usup, + const Standard_Real Vinf, + const Standard_Real Vsup) +: myLowBorder(1,2), + myUppBorder(1,2), + myDone(Standard_False) { - SetCurveCache (1, new Cache (C1, Uinf, Usup, NbU, Standard_True)); - SetCurveCache (2, new Cache (C2, Vinf, Vsup, NbV, Standard_True)); - Perform(); + myC[0] = (Standard_Address)&C1; + myC[1] = (Standard_Address)&C2; + myLowBorder(1) = Uinf; + myLowBorder(2) = Vinf; + myUppBorder(1) = Usup; + myUppBorder(2) = Vsup; } -void Extrema_GenExtCC::SetCurveCache (const Standard_Integer theRank, - const Handle(Cache)& theCache) +void Extrema_GenExtCC::SetParams (const Curve1& C1, + const Curve2& C2, + const Standard_Real Uinf, + const Standard_Real Usup, + const Standard_Real Vinf, + const Standard_Real Vsup) { - Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_GenExtCC::SetCurveCache()") - myF.SetCurve (theRank, *(Curve1*)theCache->CurvePtr()); - Standard_Integer anInd = theRank - 1; - myCache[anInd] = theCache; + myC[0] = (Standard_Address)&C1; + myC[1] = (Standard_Address)&C2; + myLowBorder(1) = Uinf; + myLowBorder(2) = Vinf; + myUppBorder(1) = Usup; + myUppBorder(2) = Vsup; } -void Extrema_GenExtCC::SetTolerance (const Standard_Real Tol) -{ - myF.SetTolerance (Tol); -} - - //============================================================================= void Extrema_GenExtCC::Perform () -/*----------------------------------------------------------------------------- -Fonction: - Recherche de toutes les distances extremales entre les courbes C1 et C2. - a partir de 2 echantillonnages (NbU,NbV). - -Methode: - L'algorithme part de l'hypothese que les echantillonnages sont suffisamment - fins pour que, s'il existe N distances extremales entre les 2 courbes, - alors il existe aussi N extrema entre les 2 ensembles de points. - Ainsi, l'algorithme consiste a partir des extrema des echantillons - pour trouver les extrema des courbes. - Les extrema sont calcules par l'algorithme math_FunctionSetRoot avec les - arguments suivants: - - myF: Extrema_FuncExtCC cree a partir de C1 et C2, - - UV: math_Vector dont les composantes sont les parametres des points de - l'extremum sur les ensembles de points, - - Tol: Min(TolU,TolV), (Prov.:math_FunctionSetRoot n'autorise pas un vecteur) - - UVinf: math_Vector dont les composantes sont les bornes inferieures de u et - v, - - UVsup: math_Vector dont les composantes sont les bornes superieures de u et - v. - -Traitement: - a- Constitution du tableau des square distances (TbDist2(0,NbU+1,0,NbV+1)): - Le tableau est volontairement etendu; les lignes 0 et NbU+1 et les - colonnes 0 et NbV+1 seront initialisees a RealFirst() ou RealLast() - pour simplifier les tests effectues dans l'etape b - (on n'a pas besoin de tester si le point est sur une extremite). - b- Calcul des extrema: - On recherche d'abord les minima et ensuite les maxima. Ces 2 traitements - se passent de facon similaire: - b.a- Initialisations: - - des 'bords' du tableau TbDist2 (a RealLast() dans le cas des minima - et a RealLast() dans le cas des maxima), - - du tableau TbSel(0,NbU+1,0,NbV+1) de selection des points pour un - calcul d'extremum local (a 0). Lorsqu'un couple de points sera - selectionne, il ne sera plus selectionnable, ainsi que les couples - adjacents (8 au maximum). - Les adresses correspondantes seront mises a 1. - b.b- Calcul des minima (ou maxima): - On boucle sur toutes les square distances du tableau TbDist2: - - recherche d'un minimum (ou maximum) sur les echantillons, - - calcul de l'extremum sur les courbes, - - mise a jour du tableau TbSel. ------------------------------------------------------------------------------*/ { myDone = Standard_False; - const Handle(Cache)& aCache1 = myCache[0]; - const Handle(Cache)& aCache2 = myCache[1]; - Standard_NullObject_Raise_if ((aCache1.IsNull() || aCache2.IsNull()), - "Extrema_GenExtCC::Perform()") + Curve1 &C1 = *(Curve1*)myC[0]; + Curve2 &C2 = *(Curve2*)myC[1]; - Standard_Integer aNbU = aCache1->NbSamples(), aNbV = aCache2->NbSamples(); - Standard_OutOfRange_Raise_if ((aNbU < 2 ||aNbV < 2), "Extrema_GenExtCC::Perform()") + Standard_Integer aNbInter[2]; + aNbInter[0] = C1.NbIntervals(GeomAbs_C2); + aNbInter[1] = C2.NbIntervals(GeomAbs_C2); + TColStd_Array1OfReal anIntervals1(1, aNbInter[0] + 1); + TColStd_Array1OfReal anIntervals2(1, aNbInter[1] + 1); + C1.Intervals(anIntervals1, GeomAbs_C2); + C2.Intervals(anIntervals2, GeomAbs_C2); -/* -a- Constitution du tableau des distances (TbDist2(0,NbU+1,0,NbV+1)): - --------------------------------------------------------------- -*/ + math_MultipleVarFunction *aFunc = new Extrema_GlobOptFuncCCC2(C1, C2); + math_GlobOptMin aFinder(aFunc, myLowBorder, myUppBorder); - //ensure that caches have been calculated - if (!aCache1->IsValid()) - aCache1->CalculatePoints(); - if (!aCache2->IsValid()) - aCache2->CalculatePoints(); + Standard_Integer i,j,k; + math_Vector aFirstBorderInterval(1,2); + math_Vector aSecondBorderInterval(1,2); + Standard_Real aF = RealLast(); // best functional value + Standard_Real aCurrF = RealLast(); // current functional value computed on current interval + for(i = 1; i <= aNbInter[0]; i++) + { + 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); -// Calcul des distances - const Handle(ArrayOfPnt)& aPntArray1 = aCache1->Points(); - const Handle(ArrayOfPnt)& aPntArray2 = aCache2->Points(); - Standard_Integer NoU, NoV; - TColStd_Array2OfReal TbDist2(0, aNbU + 1, 0, aNbV + 1); - for (NoU = 1; NoU <= aNbU; NoU++) { - const Pnt& P1 = aPntArray1->Value (NoU); - for (NoV = 1; NoV <= aNbV; NoV++) { - const Pnt& P2 = aPntArray2->Value (NoV); - TbDist2(NoU,NoV) = P1.SquareDistance(P2); + aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval); + aFinder.Perform(); + + aCurrF = aFinder.GetF(); + if (aCurrF < aF + Precision::Confusion()) + { + if (aCurrF > Abs(aF - Precision::Confusion()) || (aCurrF < 1.0e-15 && aF < 1.0e-15)) + { + Standard_Integer myTmpSolCount = aFinder.NbExtrema(); + math_Vector sol(1,2); + for(k = 1; k <= myTmpSolCount; k++) + { + aFinder.Points(k, sol); + myPoints1.Append(sol(1)); + myPoints2.Append(sol(2)); + } + mySolCount += myTmpSolCount; + } // if (aCurrF > aF - Precision::Confusion()) + else + { + aF = aCurrF; + mySolCount = aFinder.NbExtrema(); + myPoints1.Clear(); + myPoints2.Clear(); + math_Vector sol(1,2); + for(k = 1; k <= mySolCount; k++) + { + aFinder.Points(k, sol); + myPoints1.Append(sol(1)); + myPoints2.Append(sol(2)); + } + } // else + } //if (aCurrF < aF + Precision::Confusion()) } } -/* -b- Calcul des minima: - ----------------- - b.a) Initialisations: -*/ -// - generales - math_Vector Tol(1, 2); - Tol(1) = myF.Tolerance(); - Tol(2) = myF.Tolerance(); - math_Vector UV(1,2), UVinf(1,2), UVsup(1,2); - UVinf(1) = aCache1->TrimFirstParameter(); - UVinf(2) = aCache2->TrimFirstParameter(); - UVsup(1) = aCache1->TrimLastParameter(); - UVsup(2) = aCache2->TrimLastParameter(); - - myF.SubIntervalInitialize(UVinf,UVsup); - -// - des 'bords' du tableau TbDist2 - for (NoV = 0; NoV <= aNbV+1; NoV++) { - TbDist2(0,NoV) = RealLast(); - TbDist2(aNbU+1,NoV) = RealLast(); - } - for (NoU = 1; NoU <= aNbU; NoU++) { - TbDist2(NoU,0) = RealLast(); - TbDist2(NoU,aNbV+1) = RealLast(); - } - -// - du tableau TbSel(0,aNbU+1,0,aNbV+1) de selection des points - TColStd_Array2OfInteger TbSel(0,aNbU+1,0,aNbV+1); - TbSel.Init(0); - -/* - b.b) Calcul des minima: -*/ -// - recherche d un minimum sur la grille - Standard_Integer Nu, Nv; - Standard_Real Dist2; - const Handle(TColStd_HArray1OfReal)& aParamArray1 = aCache1->Parameters(); - const Handle(TColStd_HArray1OfReal)& aParamArray2 = aCache2->Parameters(); - for (NoU = 1; NoU <= aNbU; NoU++) { - for (NoV = 1; NoV <= aNbV; NoV++) { - if (TbSel(NoU,NoV) == 0) { - Dist2 = TbDist2(NoU,NoV); - if ((TbDist2(NoU-1,NoV-1) >= Dist2) && - (TbDist2(NoU-1,NoV ) >= Dist2) && - (TbDist2(NoU-1,NoV+1) >= Dist2) && - (TbDist2(NoU ,NoV-1) >= Dist2) && - (TbDist2(NoU ,NoV+1) >= Dist2) && - (TbDist2(NoU+1,NoV-1) >= Dist2) && - (TbDist2(NoU+1,NoV ) >= Dist2) && - (TbDist2(NoU+1,NoV+1) >= Dist2)) { - -// - calcul de l extremum sur la surface: - - UV(1) = aParamArray1->Value(NoU); - UV(2) = aParamArray2->Value(NoV); - - math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup); - - -// - mise a jour du tableau TbSel - for (Nu = NoU-1; Nu <= NoU+1; Nu++) { - for (Nv = NoV-1; Nv <= NoV+1; Nv++) { - TbSel(Nu,Nv) = 1; - } - } - } - } // if (TbSel(NoU,NoV) - } // for (NoV = 1; ... - } // for (NoU = 1; ... -/* -c- Calcul des maxima: - ----------------- - c.a) Initialisations: -*/ -// - des 'bords' du tableau TbDist2 - for (NoV = 0; NoV <= aNbV+1; NoV++) { - TbDist2(0,NoV) = RealFirst(); - TbDist2(aNbU+1,NoV) = RealFirst(); - } - for (NoU = 1; NoU <= aNbU; NoU++) { - TbDist2(NoU,0) = RealFirst(); - TbDist2(NoU,aNbV+1) = RealFirst(); - } - -// - du tableau TbSel(0,aNbU+1,0,aNbV+1) de selection des points - TbSel.Init(0); - /*for (NoU = 0; NoU <= aNbU+1; NoU++) { - for (NoV = 0; NoV <= aNbV+1; NoV++) { - TbSel(NoU,NoV) = 0; + // Clear solutions clusters if it is necessary. + for(i = 1; i <= mySolCount - 1; i++) + { + for(j = i + 1; j <= mySolCount; j++) + { + if ((myPoints1(i) - myPoints1(j)) < (myUppBorder(1) - myLowBorder(1)) * Precision::Confusion() && + (myPoints2(i) - myPoints2(j)) < (myUppBorder(2) - myLowBorder(2)) * Precision::Confusion()) + { + // Points with indexes i and j is in same cluster, delete j point from it. + myPoints1.Remove(j); + myPoints2.Remove(j); + j--; + mySolCount--; + } } - }*/ -/* - c.b) Calcul des maxima: -*/ -// - recherche d un maximum sur la grille - for (NoU = 1; NoU <= aNbU; NoU++) { - for (NoV = 1; NoV <= aNbV; NoV++) { - if (TbSel(NoU,NoV) == 0) { - Dist2 = TbDist2(NoU,NoV); - if ((TbDist2(NoU-1,NoV-1) <= Dist2) && - (TbDist2(NoU-1,NoV ) <= Dist2) && - (TbDist2(NoU-1,NoV+1) <= Dist2) && - (TbDist2(NoU ,NoV-1) <= Dist2) && - (TbDist2(NoU ,NoV+1) <= Dist2) && - (TbDist2(NoU+1,NoV-1) <= Dist2) && - (TbDist2(NoU+1,NoV ) <= Dist2) && - (TbDist2(NoU+1,NoV+1) <= Dist2)) { + } -// - calcul de l extremum sur la surface: - - UV(1) = aParamArray1->Value(NoU); - UV(2) = aParamArray2->Value(NoV); - - math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup); - -// - mise a jour du tableau TbSel - for (Nu = NoU-1; Nu <= NoU+1; Nu++) { - for (Nv = NoV-1; Nv <= NoV+1; Nv++) { - TbSel(Nu,Nv) = 1; - } - } - } - } // if (TbSel(NoU,NoV)) - } // for (NoV = 1; ...) - } // for (NoU = 1; ...) + delete aFunc; myDone = Standard_True; } //============================================================================= -Standard_Boolean Extrema_GenExtCC::IsDone () const { return myDone; } +Standard_Boolean Extrema_GenExtCC::IsDone () const +{ + return myDone; +} //============================================================================= Standard_Integer Extrema_GenExtCC::NbExt () const { StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()") - return myF.NbExt(); + + return mySolCount; } //============================================================================= @@ -297,14 +186,18 @@ 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()") - return myF.SquareDistance(N); + + return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N))); } //============================================================================= void Extrema_GenExtCC::Points (const Standard_Integer N, - POnC& P1, POnC& P2) const + POnC& P1, + POnC& P2) const { - StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()") - Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()") - myF.Points(N,P1,P2); -} + StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::Points()") + Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::Points()") + + P1.SetValues(myPoints1(N), Tool1::Value(*((Curve1*)myC[0]), myPoints1(N))); + P2.SetValues(myPoints2(N), Tool2::Value(*((Curve2*)myC[1]), myPoints2(N))); +} \ No newline at end of file diff --git a/src/Extrema/Extrema_GlobOptFuncCC.cxx b/src/Extrema/Extrema_GlobOptFuncCC.cxx new file mode 100644 index 0000000000..5246d3c05a --- /dev/null +++ b/src/Extrema/Extrema_GlobOptFuncCC.cxx @@ -0,0 +1,410 @@ +// Created on: 2014-01-20 +// Created by: Alexaner Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement + +#include + +#include +#include +#include +#include +#include +#include +#include + +static Standard_Integer _NbVariables() +{ + return 2; +} + +// 3d _Value +static Standard_Boolean _Value(const Adaptor3d_Curve& C1, + const Adaptor3d_Curve& C2, + const math_Vector& X, + Standard_Real& F) +{ + Standard_Real u = X(1); + Standard_Real v = X(2); + + if (u < C1.FirstParameter() || + u > C1.LastParameter() || + v < C2.FirstParameter() || + v > C2.LastParameter()) + { + return Standard_False; + } + + F = C2.Value(v).Distance(C1.Value(u)); + return Standard_True; +} + +// 2d _Value +static Standard_Boolean _Value(const Adaptor2d_Curve2d& C1, + const Adaptor2d_Curve2d& C2, + const math_Vector& X, + Standard_Real& F) +{ + Standard_Real u = X(1); + Standard_Real v = X(2); + + if (u < C1.FirstParameter() || + u > C1.LastParameter() || + v < C2.FirstParameter() || + v > C2.LastParameter()) + { + return Standard_False; + } + + F = C2.Value(v).Distance(C1.Value(u)); + return Standard_True; +} + +//! F = (x2(v) - x1(u))^2 + (y2(v) - y1(u))^2 + (z2(v) - z1(u))^2 + +// 3d _Gradient +static Standard_Boolean _Gradient(const Adaptor3d_Curve& C1, + const Adaptor3d_Curve& C2, + const math_Vector& X, + math_Vector& G) +{ + gp_Pnt C1D0, C2D0; + gp_Vec C1D1, C2D1; + + if(X(1) < C1.FirstParameter() || + X(1) > C1.LastParameter() || + X(2) < C2.FirstParameter() || + X(2) > C2.LastParameter()) + { + return Standard_False; + } + + 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(); + return Standard_True; +} + +// 2d _Graient +static Standard_Boolean _Gradient(const Adaptor2d_Curve2d& C1, + const Adaptor2d_Curve2d& C2, + const math_Vector& X, + math_Vector& G) +{ + gp_Pnt2d C1D0, C2D0; + gp_Vec2d C1D1, C2D1; + + if(X(1) < C1.FirstParameter() || + X(1) > C1.LastParameter() || + X(2) < C2.FirstParameter() || + X(2) > C2.LastParameter()) + { + return Standard_False; + } + + 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(); + G(2) = (C2D0.X() - C1D0.X()) * C2D1.X() + + (C2D0.Y() - C1D0.Y()) * C2D1.Y(); + return Standard_True; +} + +// 3d _Hessian +static Standard_Boolean _Hessian (const Adaptor3d_Curve& C1, + const Adaptor3d_Curve& C2, + const math_Vector& X, + math_Matrix & H) +{ + gp_Pnt C1D0, C2D0; + gp_Vec C1D1, C2D1; + gp_Vec C1D2, C2D2; + + if(X(1) < C1.FirstParameter() || + X(1) > C1.LastParameter() || + X(2) < C2.FirstParameter() || + X(2) > C2.LastParameter()) + { + return Standard_False; + } + + C1.D2(X(1), C1D0, C1D1, C1D2); + C2.D2(X(2), C2D0, C2D1, C2D2); + + H(1, 1) = C1D1.X() * C1D1.X() + + C1D1.Y() * C1D1.Y() + + C1D1.Z() * C1D1.Z() + - (C2D0.X() - C1D0.X()) * C1D2.X() + - (C2D0.Y() - C1D0.Y()) * C1D2.Y() + - (C2D0.Z() - C1D0.Z()) * C1D2.Z(); + + H(1, 2) = - C2D1.X() * C1D1.X() + - C2D1.Y() * C1D1.Y() + - C2D1.Z() * C1D1.Z(); + + H(2,1) = H(1,2); + + H(2,2) = C2D1.X() * C2D1.X() + + C2D1.Y() * C2D1.Y() + + C2D1.Z() * C2D1.Z() + + (C2D0.X() - C1D0.X()) * C2D2.X() + + (C2D0.Y() - C1D0.Y()) * C2D2.Y() + + (C2D0.Z() - C1D0.Z()) * C2D2.Z(); + return Standard_True; +} + +// 2d _Hessian +static Standard_Boolean _Hessian (const Adaptor2d_Curve2d& C1, + const Adaptor2d_Curve2d& C2, + const math_Vector& X, + math_Matrix & H) +{ + gp_Pnt2d C1D0, C2D0; + gp_Vec2d C1D1, C2D1; + gp_Vec2d C1D2, C2D2; + + if(X(1) < C1.FirstParameter() || + X(1) > C1.LastParameter() || + X(2) < C2.FirstParameter() || + X(2) > C2.LastParameter()) + { + return Standard_False; + } + + C1.D2(X(1), C1D0, C1D1, C1D2); + C2.D2(X(2), C2D0, C2D1, C2D2); + + H(1, 1) = C1D1.X() * C1D1.X() + + C1D1.Y() * C1D1.Y() + - (C2D0.X() - C1D0.X()) * C1D2.X() + - (C2D0.Y() - C1D0.Y()) * C1D2.Y(); + + H(1, 2) = - C2D1.X() * C1D1.X() + - C2D1.Y() * C1D1.Y(); + + H(2,1) = H(1,2); + + H(2,2) = C2D1.X() * C2D1.X() + + C2D1.Y() * C2D1.Y() + + (C2D0.X() - C1D0.X()) * C2D2.X() + + (C2D0.Y() - C1D0.Y()) * C2D2.Y(); + return Standard_True; +} + +// C0 + +//======================================================================= +//function : Extrema_GlobOptFuncCCC0 +//purpose : Constructor +//======================================================================= +Extrema_GlobOptFuncCCC0::Extrema_GlobOptFuncCCC0(const Adaptor3d_Curve& C1, + const Adaptor3d_Curve& C2) +: myC1_3d(&C1), + myC2_3d(&C2) +{ + myType = 1; +} + +//======================================================================= +//function : Extrema_GlobOptFuncCCC0 +//purpose : Constructor +//======================================================================= +Extrema_GlobOptFuncCCC0::Extrema_GlobOptFuncCCC0(const Adaptor2d_Curve2d& C1, + const Adaptor2d_Curve2d& C2) +: myC1_2d(&C1), + myC2_2d(&C2) +{ + myType = 2; +} + + +//======================================================================= +//function : NbVariables +//purpose : +//======================================================================= +Standard_Integer Extrema_GlobOptFuncCCC0::NbVariables() const +{ + return _NbVariables(); +} + +//======================================================================= +//function : Value +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCCC0::Value(const math_Vector& X,Standard_Real& F) +{ + if (myType == 1) + return _Value(*myC1_3d, *myC2_3d, X, F); + else + return _Value(*myC1_2d, *myC2_2d, X, F); +} + +// C1 + +//======================================================================= +//function : Extrema_GlobOptFuncCCC1 +//purpose : Constructor +//======================================================================= +Extrema_GlobOptFuncCCC1::Extrema_GlobOptFuncCCC1(const Adaptor3d_Curve& C1, + const Adaptor3d_Curve& C2) +: myC1_3d(&C1), + myC2_3d(&C2) +{ + myType = 1; +} + +//======================================================================= +//function : Extrema_GlobOptFuncCCC1 +//purpose : Constructor +//======================================================================= +Extrema_GlobOptFuncCCC1::Extrema_GlobOptFuncCCC1(const Adaptor2d_Curve2d& C1, + const Adaptor2d_Curve2d& C2) +: myC1_2d(&C1), + myC2_2d(&C2) +{ + myType = 2; +} + +//======================================================================= +//function : NbVariables +//purpose : +//======================================================================= +Standard_Integer Extrema_GlobOptFuncCCC1::NbVariables() const +{ + return _NbVariables(); +} + +//======================================================================= +//function : Value +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCCC1::Value(const math_Vector& X,Standard_Real& F) +{ + if (myType == 1) + return _Value(*myC1_3d, *myC2_3d, X, F); + else + return _Value(*myC1_2d, *myC2_2d, X, F); +} + +//======================================================================= +//function : Gradient +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCCC1::Gradient(const math_Vector& X,math_Vector& G) +{ + if (myType == 1) + return _Gradient(*myC1_3d, *myC2_3d, X, G); + else + return _Gradient(*myC1_2d, *myC2_2d, X, G); +} + +//======================================================================= +//function : Values +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCCC1::Values(const math_Vector& X,Standard_Real& F,math_Vector& G) +{ + return (Value(X, F) && Gradient(X, G)); +} + +// C2 + +//======================================================================= +//function : Extrema_GlobOptFuncCCC2 +//purpose : Constructor +//======================================================================= +Extrema_GlobOptFuncCCC2::Extrema_GlobOptFuncCCC2(const Adaptor3d_Curve& C1, + const Adaptor3d_Curve& C2) +: myC1_3d(&C1), + myC2_3d(&C2) +{ + myType = 1; +} + +//======================================================================= +//function : Extrema_GlobOptFuncCCC2 +//purpose : Constructor +//======================================================================= +Extrema_GlobOptFuncCCC2::Extrema_GlobOptFuncCCC2(const Adaptor2d_Curve2d& C1, + const Adaptor2d_Curve2d& C2) +: myC1_2d(&C1), + myC2_2d(&C2) +{ + myType = 2; +} + +//======================================================================= +//function : NbVariables +//purpose : +//======================================================================= +Standard_Integer Extrema_GlobOptFuncCCC2::NbVariables() const +{ + return _NbVariables(); +} + +//======================================================================= +//function : Value +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCCC2::Value(const math_Vector& X,Standard_Real& F) +{ + if (myType == 1) + return _Value(*myC1_3d, *myC2_3d, X, F); + else + return _Value(*myC1_2d, *myC2_2d, X, F); +} + +//======================================================================= +//function : Gradient +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCCC2::Gradient(const math_Vector& X,math_Vector& G) +{ + if (myType == 1) + return _Gradient(*myC1_3d, *myC2_3d, X, G); + else + return _Gradient(*myC1_2d, *myC2_2d, X, G); +} + +//======================================================================= +//function : Values +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCCC2::Values(const math_Vector& X,Standard_Real& F,math_Vector& G) +{ + return (Value(X, F) && Gradient(X, G)); +} + +//======================================================================= +//function : Values +//purpose : +//======================================================================= +Standard_Boolean Extrema_GlobOptFuncCCC2::Values(const math_Vector& X,Standard_Real& F,math_Vector& G,math_Matrix& H) +{ + Standard_Boolean isHessianComputed = Standard_False; + if (myType == 1) + isHessianComputed = _Hessian(*myC1_3d, *myC2_3d, X, H); + else + isHessianComputed = _Hessian(*myC1_2d, *myC2_2d, X, H); + + + return (Value(X, F) && Gradient(X, G) && isHessianComputed); +} diff --git a/src/Extrema/Extrema_GlobOptFuncCC.hxx b/src/Extrema/Extrema_GlobOptFuncCC.hxx new file mode 100644 index 0000000000..0369d956b8 --- /dev/null +++ b/src/Extrema/Extrema_GlobOptFuncCC.hxx @@ -0,0 +1,119 @@ +// Created on: 2014-01-20 +// Created by: Alexaner Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement + +#ifndef _Extrema_GlobOptFuncCC_HeaderFile +#define _Extrema_GlobOptFuncCC_HeaderFile + +#include +#include +#include +#include +#include +#include +#include + +//! This class implements function which operates with Eucluidean distance
+//! between 2 curves in case of C1 and C2 continuity is C0. +class Extrema_GlobOptFuncCCC0 : public math_MultipleVarFunction +{ +public: + + Standard_EXPORT Extrema_GlobOptFuncCCC0(const Adaptor3d_Curve& C1, + const Adaptor3d_Curve& C2); + + Standard_EXPORT Extrema_GlobOptFuncCCC0(const Adaptor2d_Curve2d& C1, + const Adaptor2d_Curve2d& C2); + + + + Standard_EXPORT virtual Standard_Integer NbVariables() const; + + Standard_EXPORT virtual Standard_Boolean Value(const math_Vector& X,Standard_Real& F); + + +private: + + Extrema_GlobOptFuncCCC0 & operator = (const Extrema_GlobOptFuncCCC0 & theOther); + + const Adaptor3d_Curve *myC1_3d, *myC2_3d; + const Adaptor2d_Curve2d *myC1_2d, *myC2_2d; + Standard_Integer myType; +}; + + +//! This class implements function which operates with Eucluidean distance
+//! between 2 curves in case of C1 and C2 continuity is C1. +class Extrema_GlobOptFuncCCC1 : public math_MultipleVarFunctionWithGradient +{ +public: + + Standard_EXPORT Extrema_GlobOptFuncCCC1(const Adaptor3d_Curve& C1, + const Adaptor3d_Curve& C2); + + Standard_EXPORT Extrema_GlobOptFuncCCC1(const Adaptor2d_Curve2d& C1, + const Adaptor2d_Curve2d& C2); + + Standard_EXPORT virtual Standard_Integer NbVariables() const; + + Standard_EXPORT virtual Standard_Boolean Value(const math_Vector& X,Standard_Real& F); + + Standard_EXPORT virtual Standard_Boolean Gradient(const math_Vector& X,math_Vector& G); + + Standard_EXPORT virtual Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G); + + +private: + + Extrema_GlobOptFuncCCC1 & operator = (const Extrema_GlobOptFuncCCC1 & theOther); + + const Adaptor3d_Curve *myC1_3d, *myC2_3d; + const Adaptor2d_Curve2d *myC1_2d, *myC2_2d; + Standard_Integer myType; +}; + + +//! This class implements function which operates with Eucluidean distance
+//! between 2 curves in case of C1 and C2 continuity is C2 or more. +class Extrema_GlobOptFuncCCC2 : public math_MultipleVarFunctionWithHessian +{ +public: + + Standard_EXPORT Extrema_GlobOptFuncCCC2(const Adaptor3d_Curve& C1, + const Adaptor3d_Curve& C2); + + Standard_EXPORT Extrema_GlobOptFuncCCC2(const Adaptor2d_Curve2d& C1, + const Adaptor2d_Curve2d& C2); + + Standard_EXPORT virtual Standard_Integer NbVariables() const; + + Standard_EXPORT virtual Standard_Boolean Value(const math_Vector& X,Standard_Real& F); + + Standard_EXPORT virtual Standard_Boolean Gradient(const math_Vector& X,math_Vector& G); + + Standard_EXPORT virtual Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G); + + Standard_EXPORT virtual Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G,math_Matrix& H); + + +private: + + Extrema_GlobOptFuncCCC2 & operator = (const Extrema_GlobOptFuncCCC2 & theOther); + + const Adaptor3d_Curve *myC1_3d, *myC2_3d; + const Adaptor2d_Curve2d *myC1_2d, *myC2_2d; + Standard_Integer myType; +}; + +#endif diff --git a/src/Extrema/Extrema_LocateExtCC2d.cdl b/src/Extrema/Extrema_LocateExtCC2d.cdl index 2a4dbd4c6d..330ec5af29 100644 --- a/src/Extrema/Extrema_LocateExtCC2d.cdl +++ b/src/Extrema/Extrema_LocateExtCC2d.cdl @@ -25,8 +25,7 @@ uses POnCurv2d from Extrema, Vec2d from gp, HArray1OfPnt2d from TColgp, Curve2d from Adaptor2d, - Curve2dTool from Extrema, - LCCache2d from Extrema + Curve2dTool from Extrema raises DomainError from Standard, NotDone from StdFail diff --git a/src/Extrema/FILES b/src/Extrema/FILES index 9aa286ee6d..f091b03b59 100644 --- a/src/Extrema/FILES +++ b/src/Extrema/FILES @@ -1 +1,3 @@ Extrema_HUBTreeOfSphere.hxx +Extrema_GlobOptFuncCC.hxx +Extrema_GlobOptFuncCC.cxx diff --git a/src/LocOpe/LocOpe_WiresOnShape.cxx b/src/LocOpe/LocOpe_WiresOnShape.cxx index 2bd2673df4..ed42ae8fd8 100644 --- a/src/LocOpe/LocOpe_WiresOnShape.cxx +++ b/src/LocOpe/LocOpe_WiresOnShape.cxx @@ -1158,7 +1158,7 @@ void FindInternalIntersections(const TopoDS_Edge& theEdge, } } } - if (!IntersFound) //intersection is inside "theEdge" => split + if (!IntersFound && aDist <= Precision::Confusion()) //intersection is inside "theEdge" => split { gp_Pnt aPoint = aCurve->Value(anIntPar); if (aPoint.Distance(thePnt[0]) > BRep_Tool::Tolerance(theVertices[0]) && diff --git a/src/math/FILES b/src/math/FILES index 7b4ee71660..60b777e2d9 100755 --- a/src/math/FILES +++ b/src/math/FILES @@ -8,4 +8,6 @@ math_Vector.hxx math_Vector.cxx math_IntegerVector.hxx math_IntegerVector.cxx -math_SingleTab.hxx \ No newline at end of file +math_SingleTab.hxx +math_GlobOptMin.hxx +math_GlobOptMin.cxx \ No newline at end of file diff --git a/src/math/math.cdl b/src/math/math.cdl index 2e42680fe4..ea190ffe9d 100644 --- a/src/math/math.cdl +++ b/src/math/math.cdl @@ -51,6 +51,7 @@ is deferred class MultipleVarFunctionWithHessian; deferred class FunctionSet; deferred class FunctionSetWithDerivatives; + imported GlobOptMin; class IntegerRandom; class Gauss; class GaussLeastSquare; diff --git a/src/math/math_GlobOptMin.cxx b/src/math/math_GlobOptMin.cxx new file mode 100644 index 0000000000..cd52c606f0 --- /dev/null +++ b/src/math/math_GlobOptMin.cxx @@ -0,0 +1,386 @@ +// Created on: 2014-01-20 +// Created by: Alexaner Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +const Handle(Standard_Type)& STANDARD_TYPE(math_GlobOptMin) +{ + static Handle(Standard_Type) _atype = new Standard_Type ("math_GlobOptMin", sizeof (math_GlobOptMin)); + return _atype; +} + +//======================================================================= +//function : math_GlobOptMin +//purpose : Constructor +//======================================================================= +math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc, + const math_Vector& theA, + const math_Vector& theB, + Standard_Real theC) +: myN(theFunc->NbVariables()), + myA(1, myN), + myB(1, myN), + myGlobA(1, myN), + myGlobB(1, myN), + myX(1, myN), + myTmp(1, myN), + myV(1, myN) +{ + Standard_Integer i; + + myFunc = theFunc; + myC = theC; + myZ = -1; + mySolCount = 0; + + for(i = 1; i <= myN; i++) + { + myGlobA(i) = theA(i); + myGlobB(i) = theB(i); + + myA(i) = theA(i); + myB(i) = theB(i); + } + + myDone = Standard_False; +} + +//======================================================================= +//function : SetGlobalParams +//purpose : Set params without memory allocation. +//======================================================================= +void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc, + const math_Vector& theA, + const math_Vector& theB, + Standard_Real theC) +{ + Standard_Integer i; + + myFunc = theFunc; + myC = theC; + myZ = -1; + mySolCount = 0; + + for(i = 1; i <= myN; i++) + { + myGlobA(i) = theA(i); + myGlobB(i) = theB(i); + + myA(i) = theA(i); + myB(i) = theB(i); + } + + myDone = Standard_False; +} + +//======================================================================= +//function : SetLocalParams +//purpose : Set params without memory allocation. +//======================================================================= +void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA, + const math_Vector& theLocalB) +{ + Standard_Integer i; + + myZ = -1; + mySolCount = 0; + + for(i = 1; i <= myN; i++) + { + myA(i) = theLocalA(i); + myB(i) = theLocalB(i); + } + + myDone = Standard_False; +} + +//======================================================================= +//function : ~math_GlobOptMin +//purpose : +//======================================================================= +math_GlobOptMin::~math_GlobOptMin() +{ +} + +//======================================================================= +//function : Perform +//purpose : Compute Global extremum point +//======================================================================= +// In this algo indexes started from 1, not from 0. +void math_GlobOptMin::Perform() +{ + Standard_Integer i; + + // Compute parameters range + Standard_Real minLength = RealLast(); + Standard_Real maxLength = RealFirst(); + for(i = 1; i <= myN; i++) + { + Standard_Real currentLength = myB(i) - myA(i); + if (currentLength < minLength) + minLength = currentLength; + if (currentLength > maxLength) + maxLength = currentLength; + } + + Standard_Real myTol = 1e-2; + + myE1 = minLength * myTol / myC; + myE2 = 2.0 * myTol / myC * maxLength; + myE3 = - maxLength * myTol / 4; + + // Compure start point. + math_Vector aPnt(1,2); + for(i = 1; i <= myN; i++) + { + Standard_Real currCentral = (myA(i) + myB(i)) / 2.0; + aPnt(i) = currCentral; + myY.Append(currCentral); + } + + myFunc->Value(aPnt, myF); + mySolCount++; + + computeGlobalExtremum(myN); + + myDone = Standard_True; + +} + +//======================================================================= +//function : computeLocalExtremum +//purpose : +//======================================================================= +Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt, + Standard_Real& theVal, + math_Vector& theOutPnt) +{ + Standard_Integer i; + + //Newton method + if (dynamic_cast(myFunc)) + { + math_MultipleVarFunctionWithHessian* myTmp = + dynamic_cast (myFunc); + + math_NewtonMinimum newtonMinimum(*myTmp, thePnt); + if (newtonMinimum.IsDone()) + { + newtonMinimum.Location(theOutPnt); + theVal = newtonMinimum.Minimum(); + } + else return Standard_False; + } else + + // BFGS method used. + if (dynamic_cast(myFunc)) + { + math_MultipleVarFunctionWithGradient* myTmp = + dynamic_cast (myFunc); + math_BFGS bfgs(*myTmp, thePnt); + if (bfgs.IsDone()) + { + bfgs.Location(theOutPnt); + theVal = bfgs.Minimum(); + } + else return Standard_False; + } else + + // Powell method used. + if (dynamic_cast(myFunc)) + { + math_Matrix m(1, myN, 1, myN, 0.0); + for(i = 1; i <= myN; i++) + m(1, 1) = 1.0; + + math_Powell powell(*myFunc, thePnt, m, 1e-10); + + if (powell.IsDone()) + { + powell.Location(theOutPnt); + theVal = powell.Minimum(); + } + else return Standard_False; + } + + if (isInside(theOutPnt)) + return Standard_True; + else + return Standard_False; +} + +//======================================================================= +//function : ComputeGlobalExtremum +//purpose : +//======================================================================= +void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j) +{ + Standard_Integer i; + Standard_Real d; // Functional in moved point. + Standard_Real val = RealLast(); // Local extrema computed in moved point. + Standard_Real stepBestValue = RealLast(); + math_Vector stepBestPoint(1,2); + Standard_Boolean isInside = Standard_False; + Standard_Real r; + + for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j)) + { + if (myX(j) > myB(j)) + myX(j) = myB(j); + + if (j == 1) + { + isInside = Standard_False; + myFunc->Value(myX, d); + r = (d - myF) * myZ; + if(r > myE3) + { + isInside = computeLocalExtremum(myX, val, myTmp); + } + stepBestValue = (isInside && (val < d))? val : d; + stepBestPoint = (isInside && (val < d))? myTmp : myX; + + // Solutions are close to each other. + if (Abs(stepBestValue - myF) < Precision::Confusion() * 0.01) + { + if (!isStored(stepBestPoint)) + { + if ((stepBestValue - myF) * myZ > 0.0) + myF = stepBestValue; + for(i = 1; i <= myN; i++) + myY.Append(stepBestPoint(i)); + mySolCount++; + } + } + + // New best solution. + if ((stepBestValue - myF) * myZ > Precision::Confusion() * 0.01) + { + mySolCount = 0; + myF = stepBestValue; + myY.Clear(); + for(i = 1; i <= myN; i++) + myY.Append(stepBestPoint(i)); + mySolCount++; + } + + myV(1) = myE2 + Abs(myF - d) / myC; + } + else + { + myV(j) = RealLast() / 2.0; + computeGlobalExtremum(j - 1); + } + if ((j < myN) && (myV(j + 1) > myV(j))) + { + if (myV(j) > (myB(j + 1) - myA(j + 1)) / 3.0) // Case of too big step. + myV(j + 1) = (myB(j + 1) - myA(j + 1)) / 3.0; + else + myV(j + 1) = myV(j); + } + } +} + +//======================================================================= +//function : IsInside +//purpose : +//======================================================================= +Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt) +{ + Standard_Integer i; + + for(i = 1; i <= myN; i++) + { + if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i)) + return Standard_False; + } + + return Standard_True; +} +//======================================================================= +//function : IsStored +//purpose : +//======================================================================= +Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt) +{ + Standard_Integer i,j; + Standard_Boolean isSame = Standard_True; + + for(i = 0; i < mySolCount; i++) + { + isSame = Standard_True; + for(j = 1; j <= myN; j++) + { + if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) - myA(j)) * Precision::Confusion()) + { + isSame = Standard_False; + break; + } + } + if (isSame == Standard_True) + return Standard_True; + + } + return Standard_False; +} + +//======================================================================= +//function : NbExtrema +//purpose : +//======================================================================= +Standard_Integer math_GlobOptMin::NbExtrema() +{ + return mySolCount; +} + +//======================================================================= +//function : GetF +//purpose : +//======================================================================= +Standard_Real math_GlobOptMin::GetF() +{ + return myF; +} + +//======================================================================= +//function : IsDone +//purpose : +//======================================================================= +Standard_Boolean math_GlobOptMin::isDone() +{ + return myDone; +} + +//======================================================================= +//function : Points +//purpose : +//======================================================================= +void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol) +{ + Standard_Integer j; + + for(j = 1; j <= myN; j++) + theSol(j) = myY((theIndex - 1) * myN + j); +} diff --git a/src/math/math_GlobOptMin.hxx b/src/math/math_GlobOptMin.hxx new file mode 100644 index 0000000000..55f0cca7d3 --- /dev/null +++ b/src/math/math_GlobOptMin.hxx @@ -0,0 +1,101 @@ +// Created on: 2014-01-20 +// Created by: Alexaner Malyshev +// Copyright (c) 2014-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _math_GlobOptMin_HeaderFile +#define _math_GlobOptMin_HeaderFile + +#include +#include +#include + +//! This class represents Evtushenko's algorithm of global optimization based on nonuniform mesh.
+//! Article: Yu. Evtushenko. Numerical methods for finding global extreme (case of a non-uniform mesh).
+//! U.S.S.R. Comput. Maths. Math. Phys., Vol. 11, N 6, pp. 38-54. + +class math_GlobOptMin +{ +public: + + Standard_EXPORT math_GlobOptMin(math_MultipleVarFunction* theFunc, + const math_Vector& theA, + const math_Vector& theB, + Standard_Real theC = 9); + + Standard_EXPORT void SetGlobalParams(math_MultipleVarFunction* theFunc, + const math_Vector& theA, + const math_Vector& theB, + Standard_Real theC = 9); + + Standard_EXPORT void SetLocalParams(const math_Vector& theLocalA, + const math_Vector& theLocalB); + + Standard_EXPORT ~math_GlobOptMin(); + + Standard_EXPORT void Perform(); + + //! Get best functional value. + Standard_EXPORT Standard_Real GetF(); + + //! Return count of global extremas. NbExtrema <= MAX_SOLUTIONS. + Standard_EXPORT Standard_Integer NbExtrema(); + + //! Return solution i, 1 <= i <= NbExtrema. + Standard_EXPORT void Points(const Standard_Integer theIndex, math_Vector& theSol); + +private: + + math_GlobOptMin & operator = (const math_GlobOptMin & theOther); + + Standard_Boolean computeLocalExtremum(const math_Vector& thePnt, Standard_Real& theVal, math_Vector& theOutPnt); + + void computeGlobalExtremum(Standard_Integer theIndex); + + //! Check that myA <= pnt <= myB + Standard_Boolean isInside(const math_Vector& thePnt); + + Standard_Boolean isStored(const math_Vector &thePnt); + + Standard_Boolean isDone(); + + // Input. + math_MultipleVarFunction* myFunc; + Standard_Integer myN; + math_Vector myA; // Left border on current C2 interval. + math_Vector myB; // Right border on current C2 interval. + math_Vector myGlobA; // Global left border. + math_Vector myGlobB; // Global right border. + + // Output. + Standard_Boolean myDone; + NCollection_Sequence myY;// Current solutions. + Standard_Integer mySolCount; // Count of solutions. + + // Algorithm data. + Standard_Real myZ; + Standard_Real myC; //Lipschitz constant + Standard_Real myE1; // Border coeff. + Standard_Real myE2; // Minimum step size. + Standard_Real myE3; // Local extrema starting parameter. + + math_Vector myX; // Current modified solution + math_Vector myTmp; // Current modified solution + math_Vector myV; // Steps array. + + Standard_Real myF; // Current value of Global optimum. +}; + +const Handle(Standard_Type)& TYPE(math_GlobOptMin); + +#endif diff --git a/src/math/math_NewtonMinimum.cxx b/src/math/math_NewtonMinimum.cxx index f337945220..7788c6915f 100644 --- a/src/math/math_NewtonMinimum.cxx +++ b/src/math/math_NewtonMinimum.cxx @@ -143,12 +143,19 @@ void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F, } LU.Solve(TheGradient, TheStep); - *suivant = *precedent - TheStep; + Standard_Boolean hasProblem = Standard_False; + do + { + *suivant = *precedent - TheStep; + // Gestion de la convergence + hasProblem = !(F.Value(*suivant, TheMinimum)); - // Gestion de la convergence - - F.Value(*suivant, TheMinimum); + if (hasProblem) + { + TheStep /= 2.0; + } + } while (hasProblem); if (IsConverged()) { NbConv++; } else { NbConv=0; } diff --git a/tests/bugs/modalg_5/bug23706_10 b/tests/bugs/modalg_5/bug23706_10 index 491b9f083e..3f1e0ecd1e 100755 --- a/tests/bugs/modalg_5/bug23706_10 +++ b/tests/bugs/modalg_5/bug23706_10 @@ -8,9 +8,19 @@ puts "" bsplinecurve r3 2 6 1 3 2 1 3 1 4 1 5 1 6 3 2 5 3 1 3 7 3 1 4 8 3 1 4 8 3 1 4 8 3 1 5 9 3 1 9 7 3 1 bsplinecurve r4 2 6 2 3 2.5 1 3 1 3.5 1 4 1 4.5 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 r3 r4] +extrema r3 r4 -if { [regexp "Extrema 3 is point : 4 8 3" $info] != 1 || [regexp "Extrema 8 is point : 4 8 3" $info] != 1 } { +cvalue ext_1 0 x y z +set info [dump x] +regexp "(\[-0-9.+eE\])" $info full xx +set info [dump y] +regexp "(\[-0-9.+eE\])" $info full yy +set info [dump z] +regexp "(\[-0-9.+eE\])" $info full zz + +if { sqrt(($xx - 4.0) * ($xx - 4.0) + + ($yy - 8.0) * ($yy - 8.0) + + ($zz - 3.0) * ($zz - 3.0)) > 1.0e-7} { puts "Error : Point of extrema is wrong" } else { puts "OK: Point of extrema is valid" diff --git a/tests/bugs/modalg_5/bug23706_11 b/tests/bugs/modalg_5/bug23706_11 index 207023b2a8..b1afdfe5a0 100755 --- a/tests/bugs/modalg_5/bug23706_11 +++ b/tests/bugs/modalg_5/bug23706_11 @@ -10,7 +10,7 @@ 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] != 3 } { +if { [llength $info] != 1 } { puts "Error : Extrema is wrong" } else { puts "OK: Extrema is valid" diff --git a/tests/bugs/modalg_5/bug23706_12 b/tests/bugs/modalg_5/bug23706_12 index 743e4e01ac..e9aea500f4 100755 --- a/tests/bugs/modalg_5/bug23706_12 +++ b/tests/bugs/modalg_5/bug23706_12 @@ -11,7 +11,7 @@ bsplinecurve r10 2 6 2 3 2.5 1 3 1 3.5 1 4 1 4.5 3 5 20 3 1 8 15 3 1 12 18 3 1 1 set info [extrema r9 r10] -if { [llength $info] != 6 } { +if { [llength $info] != 1 } { puts "Error : Extrema is wrong" } else { puts "OK: Extrema is valid" diff --git a/tests/bugs/modalg_5/bug23706_13 b/tests/bugs/modalg_5/bug23706_13 index 6a1f221f34..0ba9d96038 100755 --- a/tests/bugs/modalg_5/bug23706_13 +++ b/tests/bugs/modalg_5/bug23706_13 @@ -11,13 +11,11 @@ puts "" set info [2dextrema b3 b4] set status 0 -for { set i 1 } { $i <= 15 } { incr i 1 } { - regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp - if { $pp != 1.4142135623730951 } { - puts "Error : Extrema is wrong on dist $i" + regexp "dist 1: +(\[-0-9.+eE\]+)" $info full pp + if { $pp > 1.0e-7 } { + puts "Error : Extrema is wrong" set status 1 } -} if { $status != 0 } { puts "Error : Extrema is wrong" diff --git a/tests/bugs/modalg_5/bug23706_14 b/tests/bugs/modalg_5/bug23706_14 index f15381d2b1..f27011589b 100755 --- a/tests/bugs/modalg_5/bug23706_14 +++ b/tests/bugs/modalg_5/bug23706_14 @@ -11,37 +11,14 @@ puts "" set info [2dextrema b9 b10] set status 0 -for { set i 1 } { $i <= 9 } { incr i } { - regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp1 - if { $pp1 != 7.2801098892805181 } { +for { set i 1 } { $i <= 6 } { incr i 1 } { + regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp + if { abs($pp - 4.316921907096100) > 1.0e-7 } { puts "Error : Extrema is wrong on dist $i" set status 1 } } -for { set j 11 } { $j <= 19 } { incr j 1 } { - regexp "dist $j: +(\[-0-9.+eE\]+)" $info full pp2 - if { $pp2 != 7.2801098892805181 } { - puts "Error : Extrema is wrong on dist $j" - set status 1 - } -} - -regexp {dist 10: +([-0-9.+eE]+)} $info full pp3 -regexp {dist 20: +([-0-9.+eE]+)} $info full pp4 -regexp {dist 21: +([-0-9.+eE]+)} $info full pp5 -set pp_c 4.316921907096102 -set pp_ch 4.3169219070961038 - -if { $pp3 != $pp_c } { - puts "Error : Extrema is wrong on dist 10" - set status 1 -} -if { $pp4 != $pp_ch || $pp5 != $pp_ch} { - puts "Error : Extrema is wrong on dist 20 or 21" - set status 1 -} - if { $status != 0 } { puts "Error : Extrema is wrong" } else { diff --git a/tests/bugs/modalg_5/bug23706_15 b/tests/bugs/modalg_5/bug23706_15 index c5a7f0672e..0148dde86e 100755 --- a/tests/bugs/modalg_5/bug23706_15 +++ b/tests/bugs/modalg_5/bug23706_15 @@ -11,21 +11,12 @@ puts "" set info [2dextrema b7 b8] set status 0 -for { set i 2 } { $i <= 5 } { incr i } { - regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp1 - if { $pp1 !=4.3624023150195192 } { - puts "Error : Extrema is wrong on dist $i" - set status 1 - } -} -regexp {dist 1: +([-0-9.+eE]+)} $info full pp2 -set pp_ch 4.3624023150195184 - -if { $pp2 != $pp_ch } { - puts "Error : Extrema is wrong on dist 1" +regexp "dist 1: +(\[-0-9.+eE\]+)" $info full pp +if { abs($pp - 2.3246777409727225) < 1.0e-7 } { + puts "Error : Extrema is wrong" set status 1 -} + } if { $status != 0 } { puts "Error : Extrema is wrong" diff --git a/tests/bugs/modalg_5/bug24200 b/tests/bugs/modalg_5/bug24200 index 012de9bd26..e1e64efb90 100644 --- a/tests/bugs/modalg_5/bug24200 +++ b/tests/bugs/modalg_5/bug24200 @@ -10,12 +10,6 @@ restore [locate_data_file bug24200_c1] c1 restore [locate_data_file bug24200_c2] c2 set info_1 [extrema c1 c2] -if { [regexp "ext_15" $info_1] != 1 } { - puts "Error : Extrema is wrong" -} else { - puts "OK : Extrema is correct" -} - trim c1t c1 677.8 678.8 trim c2t c2 2477 2479 extrema c1t c2t @@ -34,4 +28,3 @@ if { [expr 1.*abs($checkdist - $dist)/$checkdist] > 0.1 } { } else { puts "OK: Distance is correct" } - diff --git a/tests/bugs/moddata_1/buc60825 b/tests/bugs/moddata_1/buc60825 index 5bcf9f023c..4c9920419c 100755 --- a/tests/bugs/moddata_1/buc60825 +++ b/tests/bugs/moddata_1/buc60825 @@ -11,7 +11,7 @@ checkshape aEdge2 distmini d aEdge1 aEdge2 regexp {NB RESULTS +: +([-0-9.+eE]+)} [BUC60825 aEdge1 aEdge2] full ext -if { $ext != 0 } { +if { $ext != 1 } { puts "Error : The extrema has not been calculated." } diff --git a/tests/bugs/moddata_1/buc60890 b/tests/bugs/moddata_1/buc60890 index ea2a5c8207..e2cb21f72c 100755 --- a/tests/bugs/moddata_1/buc60890 +++ b/tests/bugs/moddata_1/buc60890 @@ -15,7 +15,7 @@ circle curve_2 5.3 -0.5 2 set err [llength [2dextrema curve_1 curve_2]] -if {$err == 16} { +if {$err == 6} { puts "BUC60890 OK: Function 2dextrema works properly." } else { puts "Faulty BUC60890 : Function 2dextrema works wrongly." diff --git a/tests/bugs/moddata_2/bug862 b/tests/bugs/moddata_2/bug862 index a1a7fad1e8..c8cdbc1b51 100755 --- a/tests/bugs/moddata_2/bug862 +++ b/tests/bugs/moddata_2/bug862 @@ -11,7 +11,7 @@ restore [locate_data_file OCC862_2.draw] c2 set result [extrema c1 c2] set err [llength $result] -if { $err <= 1} { +if { $err != 1} { puts "Faulty OCC862 (amount of solution): command extrema does NOT work properly" } else { puts "OCC862 OK (amount of solution): command extrema works properly" diff --git a/tests/bugs/moddata_3/bug23994 b/tests/bugs/moddata_3/bug23994 index 9b1dd91544..be1ab53b56 100755 --- a/tests/bugs/moddata_3/bug23994 +++ b/tests/bugs/moddata_3/bug23994 @@ -27,13 +27,6 @@ extrema airfl rhomb if { [isdraw ext_1] } { mkedge result ext_1 - set length 9.09515 -} else { - puts "${BugNumber}: invalid result" -} - -if { [isdraw ext_2] } { - mkedge result ext_2 set length 5.14563 } else { puts "${BugNumber}: invalid result"