lili/OCCT
1
0
Fork 0
mirror of https://github.com/Open-Cascade-SAS/OCCT.git synced 2026-05-23 05:49:13 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
ae33c4ef42e6f49b9e222ed68fda541013bfa85a
OCCT/src/ModelingData/TKGeomBase/Extrema/Extrema_FuncExtCC.gxx
Pasukhin Dmitry 08c8482271 Coding - Applying formatting to gxx files (#730) 
- Updated GitHub Action for ASCII Check and Clang-Format
- Formatted all gxx files in the src directory
2025-09-21 09:53:01 +01:00

727 lines
21 KiB
Plaintext
Raw Blame History

// Copyright (c) 1995-1999 Matra Datavision
// Copyright (c) 1999-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.
// Modified by MPS : 21-05-97 PRO 7598
// Le nombre de solutions rendu est mySqDist.Length() et non
// mySqDist.Length()/2
// ajout des valeurs absolues dans le test d'orthogonalite de
// GetStateNumber()
/*-----------------------------------------------------------------------------
Fonctions permettant de rechercher une distance extremale entre 2 courbes
C1 et C2 (en partant de points approches C1(u0) et C2(v0)).
Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par
l'algorithme math_FunctionSetRoot.
Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont:
{ F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| }
{ F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| }
Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1
et F2 sont egales a:
{ Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2
{ Dvf1(u,v) = Dv.Du/||Du||
{ Duf2(u,v) = -Du.Dv/||Dv||
{ Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2
----------------------------------------------------------------------------*/
#include <Precision.hxx>
static const Standard_Real MinTol = 1.e-20;
static const Standard_Real TolFactor = 1.e-12;
static const Standard_Real MinStep = 1e-7;
static const Standard_Integer MaxOrder = 3;
//=================================================================================================
Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C)
{
const Standard_Integer NPoint = 10;
Standard_Real aStartParam, anEndParam;
if (C == myC1)
{
aStartParam = myUinfium;
anEndParam = myUsupremum;
}
else if (C == myC2)
{
aStartParam = myVinfium;
anEndParam = myVsupremum;
}
else
{
// Warning: No curve for tolerance computing!
return MinTol;
}
const Standard_Real aStep = (anEndParam - aStartParam) / (Standard_Real)NPoint;
Standard_Integer aNum = 0;
Standard_Real aMax = -Precision::Infinite(); // Maximum value of 1st derivative
//(it is computed with using NPoint point)
do
{
Standard_Real u = aStartParam + aNum * aStep; // parameter for every point
if (u > anEndParam)
u = anEndParam;
Pnt Ptemp; // empty point (is not used below)
Vec VDer; // 1st derivative vector
Tool1::D1(*((Curve1*)C), u, Ptemp, VDer);
Standard_Real vm = VDer.Magnitude();
if (vm > aMax)
aMax = vm;
} while (++aNum < NPoint + 1);
return Max(aMax * TolFactor, MinTol);
}
//=================================================================================================
Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol)
: myC1(0),
myC2(0),
myTol(thetol)
{
math_Vector V1(1, 2), V2(1, 2);
V1(1) = 0.0;
V2(1) = 0.0;
V1(2) = 0.0;
V2(2) = 0.0;
SubIntervalInitialize(V1, V2);
myMaxDerivOrderC1 = 0;
myTolC1 = MinTol;
myMaxDerivOrderC2 = 0;
myTolC2 = MinTol;
}
//=================================================================================================
Extrema_FuncExtCC::Extrema_FuncExtCC(const Curve1& C1, const Curve2& C2, const Standard_Real thetol)
: myC1((Standard_Address)&C1),
myC2((Standard_Address)&C2),
myTol(thetol)
{
math_Vector V1(1, 2), V2(1, 2);
V1(1) = Tool1::FirstParameter(*((Curve1*)myC1));
V2(1) = Tool1::LastParameter(*((Curve1*)myC1));
V1(2) = Tool2::FirstParameter(*((Curve2*)myC2));
V2(2) = Tool2::LastParameter(*((Curve2*)myC2));
SubIntervalInitialize(V1, V2);
switch (Tool1::GetType(*((Curve1*)myC1)))
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OffsetCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrderC1 = MaxOrder;
myTolC1 = SearchOfTolerance((Standard_Address)&C1);
break;
default:
myMaxDerivOrderC1 = 0;
myTolC1 = MinTol;
break;
}
switch (Tool2::GetType(*((Curve2*)myC2)))
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OffsetCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrderC2 = MaxOrder;
myTolC2 = SearchOfTolerance((Standard_Address)&C2);
break;
default:
myMaxDerivOrderC2 = 0;
myTolC2 = MinTol;
break;
}
}
//=================================================================================================
void Extrema_FuncExtCC::SetCurve(const Standard_Integer theRank, const Curve1& C)
{
Standard_OutOfRange_Raise_if(theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
if (theRank == 1)
{
myC1 = (Standard_Address)&C;
switch (/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType())
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OffsetCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrderC1 = MaxOrder;
myTolC1 = SearchOfTolerance((Standard_Address)&C);
break;
default:
myMaxDerivOrderC1 = 0;
myTolC1 = MinTol;
break;
}
}
else if (theRank == 2)
{
myC2 = (Standard_Address)&C;
switch (/*Tool2::GetType(*((Curve2*)myC2))*/ C.GetType())
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OffsetCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrderC2 = MaxOrder;
myTolC2 = SearchOfTolerance((Standard_Address)&C);
break;
default:
myMaxDerivOrderC2 = 0;
myTolC2 = MinTol;
break;
}
}
}
//=================================================================================================
Standard_Boolean Extrema_FuncExtCC::Value(const math_Vector& UV, math_Vector& F)
{
myU = UV(1);
myV = UV(2);
Tool1::D1(*((Curve1*)myC1), myU, myP1, myDu);
Tool2::D1(*((Curve2*)myC2), myV, myP2, myDv);
Vec P1P2(myP1, myP2);
Standard_Real Ndu = myDu.Magnitude();
if (myMaxDerivOrderC1 != 0)
{
if (Ndu <= myTolC1)
{
// Derivative is approximated by Taylor-series
const Standard_Real DivisionFactor = 1.e-3;
Standard_Real du;
if ((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
du = 0.0;
else
du = myUsupremum - myUinfium;
const Standard_Real aDelta = Max(du * DivisionFactor, MinStep);
Standard_Integer n = 1; // Derivative order
Vec V;
Standard_Boolean IsDeriveFound;
do
{
V = Tool1::DN(*((Curve1*)myC1), myU, ++n);
Ndu = V.Magnitude();
IsDeriveFound = (Ndu > myTolC1);
} while (!IsDeriveFound && n < myMaxDerivOrderC1);
if (IsDeriveFound)
{
Standard_Real u;
if (myU - myUinfium < aDelta)
u = myU + aDelta;
else
u = myU - aDelta;
Pnt P1, P2;
Tool1::D0(*((Curve1*)myC1), Min(myU, u), P1);
Tool1::D0(*((Curve1*)myC1), Max(myU, u), P2);
Vec V1(P1, P2);
Standard_Real aDirFactor = V.Dot(V1);
if (aDirFactor < 0.0)
myDu = -V;
else
myDu = V;
} // if(IsDeriveFound)
else
{
// Derivative is approximated by three points
Pnt Ptemp; //(0,0,0)-coordinate
Pnt P1, P2, P3;
Standard_Boolean IsParameterGrown;
if (myU - myUinfium < 2 * aDelta)
{
Tool1::D0(*((Curve1*)myC1), myU, P1);
Tool1::D0(*((Curve1*)myC1), myU + aDelta, P2);
Tool1::D0(*((Curve1*)myC1), myU + 2 * aDelta, P3);
IsParameterGrown = Standard_True;
}
else
{
Tool1::D0(*((Curve1*)myC1), myU - 2 * aDelta, P1);
Tool1::D0(*((Curve1*)myC1), myU - aDelta, P2);
Tool1::D0(*((Curve1*)myC1), myU, P3);
IsParameterGrown = Standard_False;
}
Vec V1(Ptemp, P1), V2(Ptemp, P2), V3(Ptemp, P3);
if (IsParameterGrown)
myDu = -3 * V1 + 4 * V2 - V3;
else
myDu = V1 - 4 * V2 + 3 * V3;
} // else of if(IsDeriveFound)
Ndu = myDu.Magnitude();
} // if (Ndu <= myTolC1) condition
} // if(myMaxDerivOrder != 0)
if (Ndu <= MinTol)
{
// Warning: 1st derivative of C1 is equal to zero!
return Standard_False;
}
Standard_Real Ndv = myDv.Magnitude();
if (myMaxDerivOrderC2 != 0)
{
if (Ndv <= myTolC2)
{
const Standard_Real DivisionFactor = 1.e-3;
Standard_Real dv;
if ((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
dv = 0.0;
else
dv = myVsupremum - myVinfium;
const Standard_Real aDelta = Max(dv * DivisionFactor, MinStep);
// Derivative is approximated by Taylor-series
Standard_Integer n = 1; // Derivative order
Vec V;
Standard_Boolean IsDeriveFound;
do
{
V = Tool2::DN(*((Curve2*)myC2), myV, ++n);
Ndv = V.Magnitude();
IsDeriveFound = (Ndv > myTolC2);
} while (!IsDeriveFound && n < myMaxDerivOrderC2);
if (IsDeriveFound)
{
Standard_Real v;
if (myV - myVinfium < aDelta)
v = myV + aDelta;
else
v = myV - aDelta;
Pnt P1, P2;
Tool2::D0(*((Curve2*)myC2), Min(myV, v), P1);
Tool2::D0(*((Curve2*)myC2), Max(myV, v), P2);
Vec V1(P1, P2);
Standard_Real aDirFactor = V.Dot(V1);
if (aDirFactor < 0.0)
myDv = -V;
else
myDv = V;
} // if(IsDeriveFound)
else
{
// Derivative is approximated by three points
Pnt Ptemp; //(0,0,0)-coordinate
Pnt P1, P2, P3;
Standard_Boolean IsParameterGrown;
if (myV - myVinfium < 2 * aDelta)
{
Tool2::D0(*((Curve2*)myC2), myV, P1);
Tool2::D0(*((Curve2*)myC2), myV + aDelta, P2);
Tool2::D0(*((Curve2*)myC2), myV + 2 * aDelta, P3);
IsParameterGrown = Standard_True;
}
else
{
Tool2::D0(*((Curve2*)myC2), myV - 2 * aDelta, P1);
Tool2::D0(*((Curve2*)myC2), myV - aDelta, P2);
Tool2::D0(*((Curve2*)myC2), myV, P3);
IsParameterGrown = Standard_False;
}
Vec V1(Ptemp, P1), V2(Ptemp, P2), V3(Ptemp, P3);
if (IsParameterGrown)
myDv = -3 * V1 + 4 * V2 - V3;
else
myDv = V1 - 4 * V2 + 3 * V3;
} // else of if(IsDeriveFound)
Ndv = myDv.Magnitude();
} // if (Ndv <= myTolC2)
} // if(myMaxDerivOrder != 0)
if (Ndv <= MinTol)
{
// 1st derivative of C2 is equal to zero!
return Standard_False;
}
F(1) = P1P2.Dot(myDu) / Ndu;
F(2) = P1P2.Dot(myDv) / Ndv;
return Standard_True;
}
//=============================================================================
Standard_Boolean Extrema_FuncExtCC::Derivatives(const math_Vector& UV, math_Matrix& Df)
{
math_Vector F(1, 2);
return Values(UV, F, Df);
}
//=============================================================================
Standard_Boolean Extrema_FuncExtCC::Values(const math_Vector& UV, math_Vector& F, math_Matrix& Df)
{
myU = UV(1);
myV = UV(2);
if (Value(UV, F) == Standard_False) // Computes F, myDu, myDv
{
// Warning: No function value found!
return Standard_False;
} // if(Value(UV, F) == Standard_False)
Vec Du, Dv, Duu, Dvv;
Tool1::D2(*((Curve1*)myC1), myU, myP1, Du, Duu);
Tool2::D2(*((Curve2*)myC2), myV, myP2, Dv, Dvv);
// Calling of "Value(...)" function change class member values.
// After running it is necessary to return to previous values.
const Standard_Real myU_old = myU, myV_old = myV;
const Pnt myP1_old = myP1, myP2_old = myP2;
const Vec myDu_old = myDu, myDv_old = myDv;
// Attention: aDelta value must be greater than same value for "Value(...)"
// function to avoid of points' collisions.
const Standard_Real DivisionFactor = 0.01;
Standard_Real du;
if ((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
du = 0.0;
else
du = myUsupremum - myUinfium;
const Standard_Real aDeltaU = Max(du * DivisionFactor, MinStep);
Standard_Real dv;
if ((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
dv = 0.0;
else
dv = myVsupremum - myVinfium;
const Standard_Real aDeltaV = Max(dv * DivisionFactor, MinStep);
Vec P1P2(myP1, myP2);
if ((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
{
// Derivative is approximated by three points
math_Vector FF1(1, 2), FF2(1, 2), FF3(1, 2);
Standard_Real F1, F2, F3;
/////////////////////////// Search of DF1_u derivative (begin) ///////////////////
if (myU - myUinfium < 2 * aDeltaU)
{
F1 = F(1);
math_Vector UV2(1, 2), UV3(1, 2);
UV2(1) = myU + aDeltaU;
UV2(2) = myV;
UV3(1) = myU + 2 * aDeltaU;
UV3(2) = myV;
if (!((Value(UV2, FF2)) && (Value(UV3, FF3))))
{
// There are many points close to singularity points and which have zero-derivative.
// Try to decrease aDelta variable's value.
return Standard_False;
}
F2 = FF2(1);
F3 = FF3(1);
Df(1, 1) = (-3 * F1 + 4 * F2 - F3) / (2.0 * aDeltaU);
} // if(myU-myUinfium < 2*aDeltaU)
else
{
F3 = F(1);
math_Vector UV2(1, 2), UV1(1, 2);
UV2(1) = myU - aDeltaU;
UV2(2) = myV;
UV1(1) = myU - 2 * aDeltaU;
UV1(2) = myV;
if (!((Value(UV2, FF2)) && (Value(UV1, FF1))))
{
// There are many points close to singularity points and which have zero-derivative.
// Try to decrease aDelta variable's value.
return Standard_False;
}
F1 = FF1(1);
F2 = FF2(1);
Df(1, 1) = (F1 - 4 * F2 + 3 * F3) / (2.0 * aDeltaU);
} // else of if(myU-myUinfium < 2*aDeltaU) condition
/////////////////////////// Search of DF1_u derivative (end) ///////////////////
// Return to previous values
myU = myU_old;
myV = myV_old;
/////////////////////////// Search of DF1_v derivative (begin) ///////////////////
if (myV - myVinfium < 2 * aDeltaV)
{
F1 = F(1);
math_Vector UV2(1, 2), UV3(1, 2);
UV2(1) = myU;
UV2(2) = myV + aDeltaV;
UV3(1) = myU;
UV3(2) = myV + 2 * aDeltaV;
if (!((Value(UV2, FF2)) && (Value(UV3, FF3))))
{
// There are many points close to singularity points and which have zero-derivative.
// Try to decrease aDelta variable's value.
return Standard_False;
}
F2 = FF2(1);
F3 = FF3(1);
Df(1, 2) = (-3 * F1 + 4 * F2 - F3) / (2.0 * aDeltaV);
} // if(myV-myVinfium < 2*aDeltaV)
else
{
F3 = F(1);
math_Vector UV2(1, 2), UV1(1, 2);
UV2(1) = myU;
UV2(2) = myV - aDeltaV;
UV1(1) = myU;
UV1(2) = myV - 2 * aDeltaV;
if (!((Value(UV2, FF2)) && (Value(UV1, FF1))))
{
// There are many points close to singularity points and which have zero-derivative.
// Try to decrease aDelta variable's value.
return Standard_False;
}
F1 = FF1(1);
F2 = FF2(1);
Df(1, 2) = (F1 - 4 * F2 + 3 * F3) / (2.0 * aDeltaV);
} // else of if(myV-myVinfium < 2*aDeltaV)
/////////////////////////// Search of DF1_v derivative (end) ///////////////////
// Return to previous values
myU = myU_old;
myV = myV_old;
myP1 = myP1_old, myP2 = myP2_old;
myDu = myDu_old, myDv = myDv_old;
} // if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
else
{
const Standard_Real Ndu = myDu.Magnitude();
Df(1, 1) = -Ndu + (P1P2.Dot(Duu) / Ndu) - F(1) * (myDu.Dot(Duu) / (Ndu * Ndu));
Df(1, 2) = myDv.Dot(myDu) / Ndu;
} // else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
if ((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
{
// Derivative is approximated by three points
math_Vector FF1(1, 2), FF2(1, 2), FF3(1, 2);
Standard_Real F1, F2, F3;
/////////////////////////// Search of DF2_v derivative (begin) ///////////////////
if (myV - myVinfium < 2 * aDeltaV)
{
F1 = F(2);
math_Vector UV2(1, 2), UV3(1, 2);
UV2(1) = myU;
UV2(2) = myV + aDeltaV;
UV3(1) = myU;
UV3(2) = myV + 2 * aDeltaV;
if (!((Value(UV2, FF2)) && (Value(UV3, FF3))))
{
// There are many points close to singularity points and which have zero-derivative.
// Try to decrease aDelta variable's value.
return Standard_False;
}
F2 = FF2(2);
F3 = FF3(2);
Df(2, 2) = (-3 * F1 + 4 * F2 - F3) / (2.0 * aDeltaV);
} // if(myV-myVinfium < 2*aDeltaV)
else
{
F3 = F(2);
math_Vector UV2(1, 2), UV1(1, 2);
UV2(1) = myU;
UV2(2) = myV - aDeltaV;
UV1(1) = myU;
UV1(2) = myV - 2 * aDeltaV;
if (!((Value(UV2, FF2)) && (Value(UV1, FF1))))
{
// There are many points close to singularity points and which have zero-derivative.
// Try to decrease aDelta variable's value.
return Standard_False;
}
F1 = FF1(2);
F2 = FF2(2);
Df(2, 2) = (F1 - 4 * F2 + 3 * F3) / (2.0 * aDeltaV);
} // else of if(myV-myVinfium < 2*aDeltaV)
/////////////////////////// Search of DF2_v derivative (end) ///////////////////
// Return to previous values
myU = myU_old;
myV = myV_old;
/////////////////////////// Search of DF2_u derivative (begin) ///////////////////
if (myU - myUinfium < 2 * aDeltaU)
{
F1 = F(2);
math_Vector UV2(1, 2), UV3(1, 2);
UV2(1) = myU + aDeltaU;
UV2(2) = myV;
UV3(1) = myU + 2 * aDeltaU;
UV3(2) = myV;
if (!((Value(UV2, FF2)) && (Value(UV3, FF3))))
{
// There are many points close to singularity points and which have zero-derivative.
// Try to decrease aDelta variable's value.
return Standard_False;
}
F2 = FF2(2);
F3 = FF3(2);
Df(2, 1) = (-3 * F1 + 4 * F2 - F3) / (2.0 * aDeltaU);
} // if(myU-myUinfium < 2*aDelta)
else
{
F3 = F(2);
math_Vector UV2(1, 2), UV1(1, 2);
UV2(1) = myU - aDeltaU;
UV2(2) = myV;
UV1(1) = myU - 2 * aDeltaU;
UV1(2) = myV;
if (!((Value(UV2, FF2)) && (Value(UV1, FF1))))
{
// There are many points close to singularity points
// and which have zero-derivative.
// Try to decrease aDelta variable's value.
return Standard_False;
}
F1 = FF1(2);
F2 = FF2(2);
Df(2, 1) = (F1 - 4 * F2 + 3 * F3) / (2.0 * aDeltaU);
} // else of if(myU-myUinfium < 2*aDeltaU)
/////////////////////////// Search of DF2_u derivative (end) ///////////////////
// Return to previous values
myU = myU_old;
myV = myV_old;
myP1 = myP1_old;
myP2 = myP2_old;
myDu = myDu_old;
myDv = myDv_old;
} // if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
else
{
Standard_Real Ndv = myDv.Magnitude();
Df(2, 2) = Ndv + (P1P2.Dot(Dvv) / Ndv) - F(2) * (myDv.Dot(Dvv) / (Ndv * Ndv));
Df(2, 1) = -myDu.Dot(myDv) / Ndv;
} // else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
return Standard_True;
} // end of function
//=============================================================================
Standard_Integer Extrema_FuncExtCC::GetStateNumber()
{
Vec Du(myDu), Dv(myDv);
Vec P1P2(myP1, myP2);
Standard_Real mod = Du.Magnitude();
if (mod > myTolC1)
{
Du /= mod;
}
mod = Dv.Magnitude();
if (mod > myTolC2)
{
Dv /= mod;
}
if (Abs(P1P2.Dot(Du)) <= myTol && Abs(P1P2.Dot(Dv)) <= myTol)
{
mySqDist.Append(myP1.SquareDistance(myP2));
myPoints.Append(POnC(myU, myP1));
myPoints.Append(POnC(myV, myP2));
}
return 0;
}
//=============================================================================
void Extrema_FuncExtCC::Points(const Standard_Integer N, POnC& P1, POnC& P2) const
{
P1 = myPoints.Value(2 * N - 1);
P2 = myPoints.Value(2 * N);
}
//=============================================================================
void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound,
const math_Vector& theSupBound)
{
myUinfium = theInfBound(1);
myUsupremum = theSupBound(1);
myVinfium = theInfBound(2);
myVsupremum = theSupBound(2);
}
//=============================================================================
Reference in New Issue View Git Blame Copy Permalink