diff --git a/src/ModelingData/TKG2d/Geom2d/FILES.cmake b/src/ModelingData/TKG2d/Geom2d/FILES.cmake index a207fb8251..c8767be16e 100644 --- a/src/ModelingData/TKG2d/Geom2d/FILES.cmake +++ b/src/ModelingData/TKG2d/Geom2d/FILES.cmake @@ -31,6 +31,7 @@ set(OCCT_Geom2d_FILES Geom2d_Line.hxx Geom2d_OffsetCurve.cxx Geom2d_OffsetCurve.hxx + Geom2d_OffsetCurveUtils.pxx Geom2d_Parabola.cxx Geom2d_Parabola.hxx Geom2d_Point.cxx diff --git a/src/ModelingData/TKG2d/Geom2d/Geom2d_OffsetCurveUtils.pxx b/src/ModelingData/TKG2d/Geom2d/Geom2d_OffsetCurveUtils.pxx new file mode 100644 index 0000000000..5db5d6b2aa --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2d/Geom2d_OffsetCurveUtils.pxx @@ -0,0 +1,507 @@ +// Copyright (c) 2015-2025 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 _Geom2d_OffsetCurveUtils_HeaderFile +#define _Geom2d_OffsetCurveUtils_HeaderFile + +#include +#include +#include +#include +#include + +#include + +//! Internal helper namespace for 2D offset curve calculations. +//! Provides static inline functions to compute offset curve point and derivatives +//! from basis curve derivatives. +//! +//! These functions are used by Geom2d_OffsetCurve and Geom2dAdaptor_Curve. +//! +//! Mathematical basis: +//! P(u) = p(u) + Offset * N / ||N|| +//! where N = (p'(u).Y, -p'(u).X) is the normal direction (rotated tangent by 90 degrees) +namespace Geom2d_OffsetCurveUtils +{ + +//! Calculates D0 (point) for 2D offset curve. +//! @param[in,out] theValue on input: basis curve point; on output: offset point +//! @param[in] theD1 first derivative of basis curve at the point +//! @param[in] theOffset offset distance value +//! @return true if successful, false if tangent vector has zero magnitude +inline bool CalculateD0(gp_Pnt2d& theValue, const gp_Vec2d& theD1, double theOffset) +{ + if (theD1.SquareMagnitude() <= gp::Resolution()) + { + return false; + } + + gp_Dir2d aNormal(theD1.Y(), -theD1.X()); + theValue.ChangeCoord().Add(aNormal.XY() * theOffset); + return true; +} + +//! Calculates D0 and D1 for 2D offset curve. +//! @param[in,out] theValue on input: basis curve point; on output: offset point +//! @param[in,out] theD1 on input: first derivative of basis; on output: offset curve D1 +//! @param[in] theD2 second derivative of basis curve +//! @param[in] theOffset offset distance value +//! @return true if successful, false if computation failed +inline bool CalculateD1(gp_Pnt2d& theValue, + gp_Vec2d& theD1, + const gp_Vec2d& theD2, + double theOffset) +{ + // P(u) = p(u) + Offset * Ndir / R + // with R = || p' ^ Z|| and Ndir = P' ^ Z + + // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) + + gp_XY Ndir(theD1.Y(), -theD1.X()); + gp_XY DNdir(theD2.Y(), -theD2.X()); + double R2 = Ndir.SquareModulus(); + double R = std::sqrt(R2); + double R3 = R * R2; + double Dr = Ndir.Dot(DNdir); + if (R3 <= gp::Resolution()) + { + if (R2 <= gp::Resolution()) + { + return false; + } + // We try another computation but the stability is not very good. + DNdir.Multiply(R); + DNdir.Subtract(Ndir.Multiplied(Dr / R)); + DNdir.Multiply(theOffset / R2); + } + else + { + // Same computation as IICURV in EUCLID-IS because the stability is better + DNdir.Multiply(theOffset / R); + DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3)); + } + + Ndir.Multiply(theOffset / R); + // P(u) + theValue.ChangeCoord().Add(Ndir); + // P'(u) + theD1.Add(gp_Vec2d(DNdir)); + return true; +} + +//! Calculates D0, D1, D2 for 2D offset curve. +//! @param[in,out] theValue on input: basis curve point; on output: offset point +//! @param[in,out] theD1 on input: first derivative of basis; on output: offset curve D1 +//! @param[in,out] theD2 on input: second derivative of basis; on output: offset curve D2 +//! @param[in] theD3 third derivative of basis curve +//! @param[in] theIsDirChange flag indicating direction change at singular point +//! @param[in] theOffset offset distance value +//! @return true if successful, false if computation failed +inline bool CalculateD2(gp_Pnt2d& theValue, + gp_Vec2d& theD1, + gp_Vec2d& theD2, + const gp_Vec2d& theD3, + bool theIsDirChange, + double theOffset) +{ + // P(u) = p(u) + Offset * Ndir / R + // with R = || p' ^ Z|| and Ndir = P' ^ Z + + // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) + + // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) + + // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2))) + + gp_XY Ndir(theD1.Y(), -theD1.X()); + gp_XY DNdir(theD2.Y(), -theD2.X()); + gp_XY D2Ndir(theD3.Y(), -theD3.X()); + double R2 = Ndir.SquareModulus(); + double R = std::sqrt(R2); + double R3 = R2 * R; + double R4 = R2 * R2; + double R5 = R3 * R2; + double Dr = Ndir.Dot(DNdir); + double D2r = Ndir.Dot(D2Ndir) + DNdir.Dot(DNdir); + if (R5 <= gp::Resolution()) + { + if (R4 <= gp::Resolution()) + { + return false; + } + // We try another computation but the stability is not very good dixit ISG. + // V2 = P" (U) : + D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R2)); + D2Ndir.Add(Ndir.Multiplied(((3.0 * Dr * Dr) / R4) - (D2r / R2))); + D2Ndir.Multiply(theOffset / R); + + // V1 = P' (U) : + DNdir.Multiply(R); + DNdir.Subtract(Ndir.Multiplied(Dr / R)); + DNdir.Multiply(theOffset / R2); + } + else + { + // Same computation as IICURV in EUCLID-IS because the stability is better. + // V2 = P" (U) : + D2Ndir.Multiply(theOffset / R); + D2Ndir.Subtract(DNdir.Multiplied(2.0 * theOffset * Dr / R3)); + D2Ndir.Add(Ndir.Multiplied(theOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3)))); + + // V1 = P' (U) + DNdir.Multiply(theOffset / R); + DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3)); + } + + Ndir.Multiply(theOffset / R); + // P(u) + theValue.ChangeCoord().Add(Ndir); + // P'(u) : + theD1.Add(gp_Vec2d(DNdir)); + // P"(u) : + if (theIsDirChange) + { + theD2.Reverse(); + } + theD2.Add(gp_Vec2d(D2Ndir)); + return true; +} + +//! Calculates D0, D1, D2, D3 for 2D offset curve. +//! @param[in,out] theValue on input: basis curve point; on output: offset point +//! @param[in,out] theD1 on input: first derivative of basis; on output: offset curve D1 +//! @param[in,out] theD2 on input: second derivative of basis; on output: offset curve D2 +//! @param[in,out] theD3 on input: third derivative of basis; on output: offset curve D3 +//! @param[in] theD4 fourth derivative of basis curve +//! @param[in] theIsDirChange flag indicating direction change at singular point +//! @param[in] theOffset offset distance value +//! @return true if successful, false if computation failed +inline bool CalculateD3(gp_Pnt2d& theValue, + gp_Vec2d& theD1, + gp_Vec2d& theD2, + gp_Vec2d& theD3, + const gp_Vec2d& theD4, + bool theIsDirChange, + double theOffset) +{ + // P(u) = p(u) + Offset * Ndir / R + // with R = || p' ^ Z|| and Ndir = P' ^ Z + + // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) + + // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) + + // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2))) + + // P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2 ) * D2Ndir - + // (3.0 * D2r / R2) * DNdir) + (3.0 * Dr * Dr / R4) * DNdir - + // (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir + + // (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir + + gp_XY Ndir(theD1.Y(), -theD1.X()); + gp_XY DNdir(theD2.Y(), -theD2.X()); + gp_XY D2Ndir(theD3.Y(), -theD3.X()); + gp_XY D3Ndir(theD4.Y(), -theD4.X()); + double R2 = Ndir.SquareModulus(); + double R = std::sqrt(R2); + double R3 = R2 * R; + double R4 = R2 * R2; + double R5 = R3 * R2; + double R6 = R3 * R3; + double R7 = R5 * R2; + double Dr = Ndir.Dot(DNdir); + double D2r = Ndir.Dot(D2Ndir) + DNdir.Dot(DNdir); + double D3r = Ndir.Dot(D3Ndir) + 3.0 * DNdir.Dot(D2Ndir); + + if (R7 <= gp::Resolution()) + { + if (R6 <= gp::Resolution()) + { + return false; + } + // We try another computation but the stability is not very good dixit ISG. + // V3 = P"' (U) : + D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * Dr / R2)); + D3Ndir.Subtract(DNdir.Multiplied(3.0 * ((D2r / R2) + (Dr * Dr / R4)))); + D3Ndir.Add( + Ndir.Multiplied(6.0 * Dr * Dr / R4 + 6.0 * Dr * D2r / R4 - 15.0 * Dr * Dr * Dr / R6 - D3r)); + D3Ndir.Multiply(theOffset / R); + // V2 = P" (U) : + R4 = R2 * R2; + D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R2)); + D2Ndir.Subtract(Ndir.Multiplied((3.0 * Dr * Dr / R4) - (D2r / R2))); + D2Ndir.Multiply(theOffset / R); + // V1 = P' (U) : + DNdir.Multiply(R); + DNdir.Subtract(Ndir.Multiplied(Dr / R)); + DNdir.Multiply(theOffset / R2); + } + else + { + // Same computation as IICURV in EUCLID-IS because the stability is better. + // V3 = P"' (U) : + D3Ndir.Multiply(theOffset / R); + D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * theOffset * Dr / R3)); + D3Ndir.Subtract(DNdir.Multiplied(((3.0 * theOffset) * ((D2r / R3) + (Dr * Dr) / R5)))); + D3Ndir.Add(Ndir.Multiplied( + (theOffset * (6.0 * Dr * Dr / R5 + 6.0 * Dr * D2r / R5 - 15.0 * Dr * Dr * Dr / R7 - D3r)))); + // V2 = P" (U) : + D2Ndir.Multiply(theOffset / R); + D2Ndir.Subtract(DNdir.Multiplied(2.0 * theOffset * Dr / R3)); + D2Ndir.Subtract(Ndir.Multiplied(theOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3)))); + // V1 = P' (U) : + DNdir.Multiply(theOffset / R); + DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3)); + } + + Ndir.Multiply(theOffset / R); + // P(u) + theValue.ChangeCoord().Add(Ndir); + // P'(u) : + theD1.Add(gp_Vec2d(DNdir)); + // P"(u) + theD2.Add(gp_Vec2d(D2Ndir)); + // P"'(u) + if (theIsDirChange) + { + theD3.Reverse(); + } + theD3.Add(gp_Vec2d(D3Ndir)); + return true; +} + +//! Adjusts derivatives at singular points where the first derivative is nearly zero. +//! Uses Taylor series approximation to find a valid tangent direction. +//! @tparam CurveType type supporting D0, DN methods (Geom2d_Curve or adaptor) +//! @param[in] theCurve basis curve for derivative evaluation +//! @param[in] theMaxDerivative maximum derivative order to compute (3 or 4) +//! @param[in] theU parameter value +//! @param[in,out] theD1 first derivative (will be adjusted) +//! @param[in,out] theD2 second derivative (will be adjusted) +//! @param[in,out] theD3 third derivative (will be adjusted) +//! @param[in,out] theD4 fourth derivative (will be adjusted if theMaxDerivative >= 4) +//! @return true if direction change detected +template +bool AdjustDerivative(const CurveType& theCurve, + int theMaxDerivative, + double theU, + gp_Vec2d& theD1, + gp_Vec2d& theD2, + gp_Vec2d& theD3, + gp_Vec2d& theD4) +{ + static const double aTol = gp::Resolution(); + static const double aMinStep = 1e-7; + static const int aMaxDerivOrder = 3; + + bool isDirectionChange = false; + const double anUinfium = theCurve.FirstParameter(); + const double anUsupremum = theCurve.LastParameter(); + + static const double DivisionFactor = 1.e-3; + double du; + if ((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + { + du = 0.0; + } + else + { + du = anUsupremum - anUinfium; + } + + const double aDelta = std::max(du * DivisionFactor, aMinStep); + + // Derivative is approximated by Taylor-series + int anIndex = 1; // Derivative order + gp_Vec2d V; + + do + { + V = theCurve.DN(theU, ++anIndex); + } while ((V.SquareMagnitude() <= aTol) && anIndex < aMaxDerivOrder); + + double u; + + if (theU - anUinfium < aDelta) + { + u = theU + aDelta; + } + else + { + u = theU - aDelta; + } + + gp_Pnt2d P1, P2; + theCurve.D0(std::min(theU, u), P1); + theCurve.D0(std::max(theU, u), P2); + + gp_Vec2d V1(P1, P2); + isDirectionChange = V.Dot(V1) < 0.0; + const double aSign = isDirectionChange ? -1.0 : 1.0; + + theD1 = V * aSign; + gp_Vec2d* aDeriv[3] = {&theD2, &theD3, &theD4}; + for (int i = 1; i < theMaxDerivative; i++) + { + *(aDeriv[i - 1]) = theCurve.DN(theU, anIndex + i) * aSign; + } + + return isDirectionChange; +} + +//! Template function for D0 evaluation of 2D offset curve. +//! Gets D1 from basis curve and computes offset point. +//! +//! @tparam BasisCurveType type of basis curve (must have D1 method) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theOffset offset distance +//! @param[out] theValue computed offset point +//! @return true if evaluation succeeded, false if normal is undefined +template +bool EvaluateD0(double theU, + const BasisCurveType& theBasisCurve, + double theOffset, + gp_Pnt2d& theValue) +{ + gp_Vec2d aD1; + theBasisCurve->D1(theU, theValue, aD1); + return CalculateD0(theValue, aD1, theOffset); +} + +//! Template function for D1 evaluation of 2D offset curve. +//! Gets D2 from basis curve and computes offset point and first derivative. +//! +//! @tparam BasisCurveType type of basis curve (must have D2 method) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theOffset offset distance +//! @param[out] theValue computed offset point +//! @param[out] theD1 computed first derivative +//! @return true if evaluation succeeded, false if computation failed +template +bool EvaluateD1(double theU, + const BasisCurveType& theBasisCurve, + double theOffset, + gp_Pnt2d& theValue, + gp_Vec2d& theD1) +{ + gp_Vec2d aD2; + theBasisCurve->D2(theU, theValue, theD1, aD2); + return CalculateD1(theValue, theD1, aD2, theOffset); +} + +//! Template function for D2 evaluation of 2D offset curve. +//! Gets D3 from basis curve and computes offset point and derivatives. +//! Handles singular points where first derivative is nearly zero. +//! +//! @tparam BasisCurveType type of basis curve (must have D3 and DN methods) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theOffset offset distance +//! @param[out] theValue computed offset point +//! @param[out] theD1 computed first derivative +//! @param[out] theD2 computed second derivative +//! @return true if evaluation succeeded, false if computation failed +template +bool EvaluateD2(double theU, + const BasisCurveType& theBasisCurve, + double theOffset, + gp_Pnt2d& theValue, + gp_Vec2d& theD1, + gp_Vec2d& theD2) +{ + gp_Vec2d aD3; + theBasisCurve->D3(theU, theValue, theD1, theD2, aD3); + + bool isDirectionChange = false; + if (theD1.SquareMagnitude() <= gp::Resolution()) + { + gp_Vec2d aDummyD4; + isDirectionChange = AdjustDerivative(*theBasisCurve, 3, theU, theD1, theD2, aD3, aDummyD4); + } + + return CalculateD2(theValue, theD1, theD2, aD3, isDirectionChange, theOffset); +} + +//! Template function for D3 evaluation of 2D offset curve. +//! Gets D3 and D4 from basis curve and computes offset point and derivatives. +//! Handles singular points where first derivative is nearly zero. +//! +//! @tparam BasisCurveType type of basis curve (must have D3 and DN methods) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theOffset offset distance +//! @param[out] theValue computed offset point +//! @param[out] theD1 computed first derivative +//! @param[out] theD2 computed second derivative +//! @param[out] theD3 computed third derivative +//! @return true if evaluation succeeded, false if computation failed +template +bool EvaluateD3(double theU, + const BasisCurveType& theBasisCurve, + double theOffset, + gp_Pnt2d& theValue, + gp_Vec2d& theD1, + gp_Vec2d& theD2, + gp_Vec2d& theD3) +{ + theBasisCurve->D3(theU, theValue, theD1, theD2, theD3); + gp_Vec2d aD4 = theBasisCurve->DN(theU, 4); + + bool isDirectionChange = false; + if (theD1.SquareMagnitude() <= gp::Resolution()) + { + isDirectionChange = AdjustDerivative(*theBasisCurve, 4, theU, theD1, theD2, theD3, aD4); + } + + return CalculateD3(theValue, theD1, theD2, theD3, aD4, isDirectionChange, theOffset); +} + +//! Template function for DN evaluation of 2D offset curve. +//! Handles derivatives up to order 3 using D1/D2/D3 methods. +//! For derivatives > 3, returns the basis curve derivative directly +//! (offset contribution is negligible for high-order derivatives). +//! +//! @tparam BasisCurveType type of basis curve (must have D1, D2, D3, DN methods) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theOffset offset distance +//! @param[in] theN derivative order (must be >= 1) +//! @param[out] theDN computed N-th derivative +//! @return true if evaluation succeeded, false if computation failed +template +bool EvaluateDN(double theU, + const BasisCurveType& theBasisCurve, + double theOffset, + int theN, + gp_Vec2d& theDN) +{ + gp_Pnt2d aPnt; + gp_Vec2d aDummy; + switch (theN) + { + case 1: + return EvaluateD1(theU, theBasisCurve, theOffset, aPnt, theDN); + case 2: + return EvaluateD2(theU, theBasisCurve, theOffset, aPnt, aDummy, theDN); + case 3: + return EvaluateD3(theU, theBasisCurve, theOffset, aPnt, aDummy, aDummy, theDN); + default: + // For derivatives > 3, return basis curve derivative (no offset contribution) + theDN = theBasisCurve->DN(theU, theN); + return true; + } +} + +} // namespace Geom2d_OffsetCurveUtils + +#endif // _Geom2d_OffsetCurveUtils_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.cxx b/src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.cxx index 8db3a7f080..ac493d331b 100644 --- a/src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.cxx +++ b/src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.cxx @@ -13,9 +13,10 @@ // commercial license or contractual agreement. #include -#include + +#include +#include #include -#include IMPLEMENT_STANDARD_RTTIEXT(Geom2dEvaluator_OffsetCurve, Geom2dEvaluator_Curve) @@ -35,80 +36,136 @@ Geom2dEvaluator_OffsetCurve::Geom2dEvaluator_OffsetCurve(const Handle(Geom2dAdap { } +//================================================================================================== + void Geom2dEvaluator_OffsetCurve::D0(const Standard_Real theU, gp_Pnt2d& theValue) const { - gp_Vec2d aD1; - BaseD1(theU, theValue, aD1); - Geom2dEvaluator::CalculateD0(theValue, aD1, myOffset); + bool isOK = false; + if (!myBaseAdaptor.IsNull()) + { + isOK = Geom2d_OffsetCurveUtils::EvaluateD0(theU, myBaseAdaptor, myOffset, theValue); + } + else + { + isOK = Geom2d_OffsetCurveUtils::EvaluateD0(theU, myBaseCurve, myOffset, theValue); + } + + if (!isOK) + { + throw Geom2d_UndefinedValue("Geom2dEvaluator_OffsetCurve::D0(): Unable to calculate normal"); + } } +//================================================================================================== + void Geom2dEvaluator_OffsetCurve::D1(const Standard_Real theU, gp_Pnt2d& theValue, gp_Vec2d& theD1) const { - gp_Vec2d aD2; - BaseD2(theU, theValue, theD1, aD2); - Geom2dEvaluator::CalculateD1(theValue, theD1, aD2, myOffset); + bool isOK = false; + if (!myBaseAdaptor.IsNull()) + { + isOK = Geom2d_OffsetCurveUtils::EvaluateD1(theU, myBaseAdaptor, myOffset, theValue, theD1); + } + else + { + isOK = Geom2d_OffsetCurveUtils::EvaluateD1(theU, myBaseCurve, myOffset, theValue, theD1); + } + + if (!isOK) + { + throw Geom2d_UndefinedValue("Geom2dEvaluator_OffsetCurve::D1(): Unable to calculate normal"); + } } +//================================================================================================== + void Geom2dEvaluator_OffsetCurve::D2(const Standard_Real theU, gp_Pnt2d& theValue, gp_Vec2d& theD1, gp_Vec2d& theD2) const { - gp_Vec2d aD3; - BaseD3(theU, theValue, theD1, theD2, aD3); - - Standard_Boolean isDirectionChange = Standard_False; - if (theD1.SquareMagnitude() <= gp::Resolution()) + bool isOK = false; + if (!myBaseAdaptor.IsNull()) { - gp_Vec2d aDummyD4; - isDirectionChange = AdjustDerivative(3, theU, theD1, theD2, aD3, aDummyD4); + isOK = + Geom2d_OffsetCurveUtils::EvaluateD2(theU, myBaseAdaptor, myOffset, theValue, theD1, theD2); + } + else + { + isOK = Geom2d_OffsetCurveUtils::EvaluateD2(theU, myBaseCurve, myOffset, theValue, theD1, theD2); } - Geom2dEvaluator::CalculateD2(theValue, theD1, theD2, aD3, isDirectionChange, myOffset); + if (!isOK) + { + throw Geom2d_UndefinedValue("Geom2dEvaluator_OffsetCurve::D2(): Unable to calculate normal"); + } } +//================================================================================================== + void Geom2dEvaluator_OffsetCurve::D3(const Standard_Real theU, gp_Pnt2d& theValue, gp_Vec2d& theD1, gp_Vec2d& theD2, gp_Vec2d& theD3) const { - gp_Vec2d aD4; - BaseD4(theU, theValue, theD1, theD2, theD3, aD4); + bool isOK = false; + if (!myBaseAdaptor.IsNull()) + { + isOK = Geom2d_OffsetCurveUtils::EvaluateD3(theU, + myBaseAdaptor, + myOffset, + theValue, + theD1, + theD2, + theD3); + } + else + { + isOK = Geom2d_OffsetCurveUtils::EvaluateD3(theU, + myBaseCurve, + myOffset, + theValue, + theD1, + theD2, + theD3); + } - Standard_Boolean isDirectionChange = Standard_False; - if (theD1.SquareMagnitude() <= gp::Resolution()) - isDirectionChange = AdjustDerivative(4, theU, theD1, theD2, theD3, aD4); - - Geom2dEvaluator::CalculateD3(theValue, theD1, theD2, theD3, aD4, isDirectionChange, myOffset); + if (!isOK) + { + throw Geom2d_UndefinedValue("Geom2dEvaluator_OffsetCurve::D3(): Unable to calculate normal"); + } } +//================================================================================================== + gp_Vec2d Geom2dEvaluator_OffsetCurve::DN(const Standard_Real theU, const Standard_Integer theDeriv) const { Standard_RangeError_Raise_if(theDeriv < 1, "Geom2dEvaluator_OffsetCurve::DN(): theDeriv < 1"); - gp_Pnt2d aPnt; - gp_Vec2d aDummy, aDN; - switch (theDeriv) + gp_Vec2d aResult; + bool isOK = false; + if (!myBaseAdaptor.IsNull()) { - case 1: - D1(theU, aPnt, aDN); - break; - case 2: - D2(theU, aPnt, aDummy, aDN); - break; - case 3: - D3(theU, aPnt, aDummy, aDummy, aDN); - break; - default: - aDN = BaseDN(theU, theDeriv); + isOK = Geom2d_OffsetCurveUtils::EvaluateDN(theU, myBaseAdaptor, myOffset, theDeriv, aResult); } - return aDN; + else + { + isOK = Geom2d_OffsetCurveUtils::EvaluateDN(theU, myBaseCurve, myOffset, theDeriv, aResult); + } + + if (!isOK) + { + throw Geom2d_UndefinedValue("Geom2dEvaluator_OffsetCurve::DN(): Unable to calculate normal"); + } + + return aResult; } +//================================================================================================== + Handle(Geom2dEvaluator_Curve) Geom2dEvaluator_OffsetCurve::ShallowCopy() const { Handle(Geom2dEvaluator_OffsetCurve) aCopy; @@ -125,138 +182,3 @@ Handle(Geom2dEvaluator_Curve) Geom2dEvaluator_OffsetCurve::ShallowCopy() const return aCopy; } - -void Geom2dEvaluator_OffsetCurve::BaseD0(const Standard_Real theU, gp_Pnt2d& theValue) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D0(theU, theValue); - else - myBaseCurve->D0(theU, theValue); -} - -void Geom2dEvaluator_OffsetCurve::BaseD1(const Standard_Real theU, - gp_Pnt2d& theValue, - gp_Vec2d& theD1) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D1(theU, theValue, theD1); - else - myBaseCurve->D1(theU, theValue, theD1); -} - -void Geom2dEvaluator_OffsetCurve::BaseD2(const Standard_Real theU, - gp_Pnt2d& theValue, - gp_Vec2d& theD1, - gp_Vec2d& theD2) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D2(theU, theValue, theD1, theD2); - else - myBaseCurve->D2(theU, theValue, theD1, theD2); -} - -void Geom2dEvaluator_OffsetCurve::BaseD3(const Standard_Real theU, - gp_Pnt2d& theValue, - gp_Vec2d& theD1, - gp_Vec2d& theD2, - gp_Vec2d& theD3) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D3(theU, theValue, theD1, theD2, theD3); - else - myBaseCurve->D3(theU, theValue, theD1, theD2, theD3); -} - -void Geom2dEvaluator_OffsetCurve::BaseD4(const Standard_Real theU, - gp_Pnt2d& theValue, - gp_Vec2d& theD1, - gp_Vec2d& theD2, - gp_Vec2d& theD3, - gp_Vec2d& theD4) const -{ - if (!myBaseAdaptor.IsNull()) - { - myBaseAdaptor->D3(theU, theValue, theD1, theD2, theD3); - theD4 = myBaseAdaptor->DN(theU, 4); - } - else - { - myBaseCurve->D3(theU, theValue, theD1, theD2, theD3); - theD4 = myBaseCurve->DN(theU, 4); - } -} - -gp_Vec2d Geom2dEvaluator_OffsetCurve::BaseDN(const Standard_Real theU, - const Standard_Integer theDeriv) const -{ - if (!myBaseAdaptor.IsNull()) - return myBaseAdaptor->DN(theU, theDeriv); - return myBaseCurve->DN(theU, theDeriv); -} - -Standard_Boolean Geom2dEvaluator_OffsetCurve::AdjustDerivative( - const Standard_Integer theMaxDerivative, - const Standard_Real theU, - gp_Vec2d& theD1, - gp_Vec2d& theD2, - gp_Vec2d& theD3, - gp_Vec2d& theD4) const -{ - static const Standard_Real aTol = gp::Resolution(); - static const Standard_Real aMinStep = 1e-7; - static const Standard_Integer aMaxDerivOrder = 3; - - Standard_Boolean isDirectionChange = Standard_False; - Standard_Real anUinfium; - Standard_Real anUsupremum; - if (!myBaseAdaptor.IsNull()) - { - anUinfium = myBaseAdaptor->FirstParameter(); - anUsupremum = myBaseAdaptor->LastParameter(); - } - else - { - anUinfium = myBaseCurve->FirstParameter(); - anUsupremum = myBaseCurve->LastParameter(); - } - - static const Standard_Real DivisionFactor = 1.e-3; - Standard_Real du; - if ((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) - du = 0.0; - else - du = anUsupremum - anUinfium; - - const Standard_Real aDelta = std::max(du * DivisionFactor, aMinStep); - - // Derivative is approximated by Taylor-series - Standard_Integer anIndex = 1; // Derivative order - gp_Vec2d V; - - do - { - V = BaseDN(theU, ++anIndex); - } while ((V.SquareMagnitude() <= aTol) && anIndex < aMaxDerivOrder); - - Standard_Real u; - - if (theU - anUinfium < aDelta) - u = theU + aDelta; - else - u = theU - aDelta; - - gp_Pnt2d P1, P2; - BaseD0(std::min(theU, u), P1); - BaseD0(std::max(theU, u), P2); - - gp_Vec2d V1(P1, P2); - isDirectionChange = V.Dot(V1) < 0.0; - Standard_Real aSign = isDirectionChange ? -1.0 : 1.0; - - theD1 = V * aSign; - gp_Vec2d* aDeriv[3] = {&theD2, &theD3, &theD4}; - for (Standard_Integer i = 1; i < theMaxDerivative; i++) - *(aDeriv[i - 1]) = BaseDN(theU, anIndex + i) * aSign; - - return isDirectionChange; -} diff --git a/src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.hxx b/src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.hxx index b6db1cd07b..9bd937e59c 100644 --- a/src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.hxx +++ b/src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.hxx @@ -57,38 +57,6 @@ public: DEFINE_STANDARD_RTTIEXT(Geom2dEvaluator_OffsetCurve, Geom2dEvaluator_Curve) -private: - //! Calculate value of base curve/adaptor - void BaseD0(const Standard_Real theU, gp_Pnt2d& theValue) const; - //! Calculate value and first derivatives of base curve/adaptor - void BaseD1(const Standard_Real theU, gp_Pnt2d& theValue, gp_Vec2d& theD1) const; - //! Calculate value, first and second derivatives of base curve/adaptor - void BaseD2(const Standard_Real theU, gp_Pnt2d& theValue, gp_Vec2d& theD1, gp_Vec2d& theD2) const; - //! Calculate value, first, second and third derivatives of base curve/adaptor - void BaseD3(const Standard_Real theU, - gp_Pnt2d& theValue, - gp_Vec2d& theD1, - gp_Vec2d& theD2, - gp_Vec2d& theD3) const; - //! Calculate value and derivatives till 4th of base curve/adaptor - void BaseD4(const Standard_Real theU, - gp_Pnt2d& theValue, - gp_Vec2d& theD1, - gp_Vec2d& theD2, - gp_Vec2d& theD3, - gp_Vec2d& theD4) const; - //! Calculate N-th derivative of base curve/adaptor - gp_Vec2d BaseDN(const Standard_Real theU, const Standard_Integer theDeriv) const; - - // Recalculate derivatives in the singular point - // Returns true if the direction of derivatives is changed - Standard_Boolean AdjustDerivative(const Standard_Integer theMaxDerivative, - const Standard_Real theU, - gp_Vec2d& theD1, - gp_Vec2d& theD2, - gp_Vec2d& theD3, - gp_Vec2d& theD4) const; - private: Handle(Geom2d_Curve) myBaseCurve; Handle(Geom2dAdaptor_Curve) myBaseAdaptor; diff --git a/src/ModelingData/TKG3d/Geom/FILES.cmake b/src/ModelingData/TKG3d/Geom/FILES.cmake index 044def655d..6aeb9ecdba 100644 --- a/src/ModelingData/TKG3d/Geom/FILES.cmake +++ b/src/ModelingData/TKG3d/Geom/FILES.cmake @@ -37,6 +37,7 @@ set(OCCT_Geom_FILES Geom_Direction.cxx Geom_Direction.hxx Geom_ElementarySurface.cxx + Geom_ExtrusionUtils.pxx Geom_ElementarySurface.hxx Geom_Ellipse.cxx Geom_Ellipse.hxx @@ -49,8 +50,10 @@ set(OCCT_Geom_FILES Geom_Line.hxx Geom_OffsetCurve.cxx Geom_OffsetCurve.hxx + Geom_OffsetCurveUtils.pxx Geom_OffsetSurface.cxx Geom_OffsetSurface.hxx + Geom_OffsetSurfaceUtils.pxx Geom_OsculatingSurface.cxx Geom_OsculatingSurface.hxx Geom_Parabola.cxx @@ -61,6 +64,7 @@ set(OCCT_Geom_FILES Geom_Point.hxx Geom_RectangularTrimmedSurface.cxx Geom_RectangularTrimmedSurface.hxx + Geom_RevolutionUtils.pxx Geom_SequenceOfBSplineSurface.hxx Geom_SphericalSurface.cxx Geom_SphericalSurface.hxx diff --git a/src/ModelingData/TKG3d/Geom/Geom_ExtrusionUtils.pxx b/src/ModelingData/TKG3d/Geom/Geom_ExtrusionUtils.pxx new file mode 100644 index 0000000000..3b9c73ca20 --- /dev/null +++ b/src/ModelingData/TKG3d/Geom/Geom_ExtrusionUtils.pxx @@ -0,0 +1,172 @@ +// Copyright (c) 2025 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 _Geom_ExtrusionUtils_HeaderFile +#define _Geom_ExtrusionUtils_HeaderFile + +#include +#include +#include + +//! @file Geom_ExtrusionUtils.pxx +//! @brief Shared utility functions for extrusion surface evaluation. +//! +//! This file provides template functions for evaluating points and derivatives +//! on linear extrusion surfaces. The functions are templated to work with both +//! Geom_Curve (for Geom_SurfaceOfLinearExtrusion) and Adaptor3d_Curve +//! (for GeomAdaptor_SurfaceOfLinearExtrusion). +//! +//! Extrusion surface: P(U,V) = C(U) + V * Direction + +namespace Geom_ExtrusionUtils +{ + +//! Evaluates point on extrusion surface. +//! @tparam CurveType Type supporting D0(param, point) method +//! @param theU Parameter along the basis curve +//! @param theV Parameter along the extrusion direction +//! @param theBasis Basis curve +//! @param theDir Extrusion direction XYZ (must be normalized) +//! @param theP [out] Evaluated point +template +inline void D0(const double theU, + const double theV, + const CurveType& theBasis, + const gp_XYZ& theDir, + gp_Pnt& theP) +{ + theBasis.D0(theU, theP); + theP.SetXYZ(theP.XYZ() + theV * theDir); +} + +//! Evaluates point and first derivatives on extrusion surface. +//! @tparam CurveType Type supporting D1(param, point, vec) method +//! @param theU Parameter along the basis curve +//! @param theV Parameter along the extrusion direction +//! @param theBasis Basis curve +//! @param theDir Extrusion direction XYZ (must be normalized) +//! @param theP [out] Evaluated point +//! @param theD1U [out] First derivative with respect to U +//! @param theD1V [out] First derivative with respect to V +template +inline void D1(const double theU, + const double theV, + const CurveType& theBasis, + const gp_XYZ& theDir, + gp_Pnt& theP, + gp_Vec& theD1U, + gp_Vec& theD1V) +{ + theBasis.D1(theU, theP, theD1U); + theP.SetXYZ(theP.XYZ() + theV * theDir); + theD1V.SetXYZ(theDir); +} + +//! Evaluates point, first and second derivatives on extrusion surface. +//! @tparam CurveType Type supporting D2(param, point, vec, vec) method +//! @param theU Parameter along the basis curve +//! @param theV Parameter along the extrusion direction +//! @param theBasis Basis curve +//! @param theDir Extrusion direction XYZ (must be normalized) +//! @param theP [out] Evaluated point +//! @param theD1U [out] First derivative with respect to U +//! @param theD1V [out] First derivative with respect to V +//! @param theD2U [out] Second derivative with respect to U +//! @param theD2V [out] Second derivative with respect to V +//! @param theD2UV [out] Mixed second derivative +template +inline void D2(const double theU, + const double theV, + const CurveType& theBasis, + const gp_XYZ& theDir, + gp_Pnt& theP, + gp_Vec& theD1U, + gp_Vec& theD1V, + gp_Vec& theD2U, + gp_Vec& theD2V, + gp_Vec& theD2UV) +{ + theBasis.D2(theU, theP, theD1U, theD2U); + theP.SetXYZ(theP.XYZ() + theV * theDir); + theD1V.SetXYZ(theDir); + theD2V.SetCoord(0.0, 0.0, 0.0); + theD2UV.SetCoord(0.0, 0.0, 0.0); +} + +//! Evaluates point, first, second and third derivatives on extrusion surface. +//! @tparam CurveType Type supporting D3(param, point, vec, vec, vec) method +//! @param theU Parameter along the basis curve +//! @param theV Parameter along the extrusion direction +//! @param theBasis Basis curve +//! @param theDir Extrusion direction XYZ (must be normalized) +//! @param theP [out] Evaluated point +//! @param theD1U [out] First derivative with respect to U +//! @param theD1V [out] First derivative with respect to V +//! @param theD2U [out] Second derivative with respect to U +//! @param theD2V [out] Second derivative with respect to V +//! @param theD2UV [out] Mixed second derivative +//! @param theD3U [out] Third derivative with respect to U +//! @param theD3V [out] Third derivative with respect to V +//! @param theD3UUV [out] Mixed third derivative (UUV) +//! @param theD3UVV [out] Mixed third derivative (UVV) +template +inline void D3(const double theU, + const double theV, + const CurveType& theBasis, + const gp_XYZ& theDir, + gp_Pnt& theP, + gp_Vec& theD1U, + gp_Vec& theD1V, + gp_Vec& theD2U, + gp_Vec& theD2V, + gp_Vec& theD2UV, + gp_Vec& theD3U, + gp_Vec& theD3V, + gp_Vec& theD3UUV, + gp_Vec& theD3UVV) +{ + theBasis.D3(theU, theP, theD1U, theD2U, theD3U); + theP.SetXYZ(theP.XYZ() + theV * theDir); + theD1V.SetXYZ(theDir); + theD2V.SetCoord(0.0, 0.0, 0.0); + theD2UV.SetCoord(0.0, 0.0, 0.0); + theD3V.SetCoord(0.0, 0.0, 0.0); + theD3UUV.SetCoord(0.0, 0.0, 0.0); + theD3UVV.SetCoord(0.0, 0.0, 0.0); +} + +//! Evaluates N-th derivative on extrusion surface. +//! @tparam CurveType Type supporting DN(param, order) method +//! @param theU Parameter along the basis curve +//! @param theBasis Basis curve +//! @param theDir Extrusion direction XYZ (must be normalized) +//! @param theDerU Derivative order with respect to U +//! @param theDerV Derivative order with respect to V +//! @return The derivative vector +template +inline gp_Vec DN(const double theU, + const CurveType& theBasis, + const gp_XYZ& theDir, + const int theDerU, + const int theDerV) +{ + if (theDerV == 0) + return theBasis.DN(theU, theDerU); + else if (theDerU == 0 && theDerV == 1) + return gp_Vec(theDir); + return gp_Vec(0.0, 0.0, 0.0); +} + +} // namespace Geom_ExtrusionUtils + +#endif // _Geom_ExtrusionUtils_HeaderFile diff --git a/src/ModelingData/TKG3d/Geom/Geom_OffsetCurveUtils.pxx b/src/ModelingData/TKG3d/Geom/Geom_OffsetCurveUtils.pxx new file mode 100644 index 0000000000..e0e6c7e142 --- /dev/null +++ b/src/ModelingData/TKG3d/Geom/Geom_OffsetCurveUtils.pxx @@ -0,0 +1,532 @@ +// Copyright (c) 2015-2025 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 _Geom_OffsetCurveUtils_HeaderFile +#define _Geom_OffsetCurveUtils_HeaderFile + +#include +#include +#include +#include +#include +#include + +#include + +//! Internal helper namespace for 3D offset curve calculations. +//! Provides static inline functions to compute offset curve point and derivatives +//! from basis curve derivatives. +//! +//! These functions are used by Geom_OffsetCurve and GeomAdaptor_Curve. +//! +//! Mathematical basis: +//! P(u) = p(u) + Offset * N / ||N|| +//! where N = p'(u) ^ Direction is the local normal direction +namespace Geom_OffsetCurveUtils +{ + +//! Calculates D0 (point) for 3D offset curve. +//! @param[in,out] theValue on input: basis curve point; on output: offset point +//! @param[in] theD1 first derivative of basis curve at the point +//! @param[in] theDir offset reference direction +//! @param[in] theOffset offset distance value +//! @return true if successful, false if normal vector has zero magnitude +inline bool CalculateD0(gp_Pnt& theValue, + const gp_Vec& theD1, + const gp_Dir& theDir, + double theOffset) +{ + gp_XYZ Ndir = (theD1.XYZ()).Crossed(theDir.XYZ()); + double R = Ndir.Modulus(); + if (R <= gp::Resolution()) + { + return false; + } + + Ndir.Multiply(theOffset / R); + theValue.ChangeCoord().Add(Ndir); + return true; +} + +//! Calculates D0 and D1 for 3D offset curve. +//! @param[in,out] theValue on input: basis curve point; on output: offset point +//! @param[in,out] theD1 on input: first derivative of basis; on output: offset curve D1 +//! @param[in] theD2 second derivative of basis curve +//! @param[in] theDir offset reference direction +//! @param[in] theOffset offset distance value +//! @return true if successful, false if computation failed +inline bool CalculateD1(gp_Pnt& theValue, + gp_Vec& theD1, + const gp_Vec& theD2, + const gp_Dir& theDir, + double theOffset) +{ + // P(u) = p(u) + Offset * Ndir / R + // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction) + + // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) + + gp_XYZ Ndir = (theD1.XYZ()).Crossed(theDir.XYZ()); + gp_XYZ DNdir = (theD2.XYZ()).Crossed(theDir.XYZ()); + double R2 = Ndir.SquareModulus(); + double R = std::sqrt(R2); + double R3 = R * R2; + double Dr = Ndir.Dot(DNdir); + if (R3 <= gp::Resolution()) + { + if (R2 <= gp::Resolution()) + { + return false; + } + // We try another computation but the stability is not very good. + DNdir.Multiply(R); + DNdir.Subtract(Ndir.Multiplied(Dr / R)); + DNdir.Multiply(theOffset / R2); + } + else + { + // Same computation as IICURV in EUCLID-IS because the stability is better + DNdir.Multiply(theOffset / R); + DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3)); + } + + Ndir.Multiply(theOffset / R); + // P(u) + theValue.ChangeCoord().Add(Ndir); + // P'(u) + theD1.Add(gp_Vec(DNdir)); + return true; +} + +//! Calculates D0, D1, D2 for 3D offset curve. +//! @param[in,out] theValue on input: basis curve point; on output: offset point +//! @param[in,out] theD1 on input: first derivative of basis; on output: offset curve D1 +//! @param[in,out] theD2 on input: second derivative of basis; on output: offset curve D2 +//! @param[in] theD3 third derivative of basis curve +//! @param[in] theDir offset reference direction +//! @param[in] theOffset offset distance value +//! @param[in] theIsDirChange flag indicating direction change at singular point +//! @return true if successful, false if computation failed +inline bool CalculateD2(gp_Pnt& theValue, + gp_Vec& theD1, + gp_Vec& theD2, + const gp_Vec& theD3, + const gp_Dir& theDir, + double theOffset, + bool theIsDirChange) +{ + // P(u) = p(u) + Offset * Ndir / R + // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction) + + // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) + + // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) + + // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2))) + + gp_XYZ Ndir = (theD1.XYZ()).Crossed(theDir.XYZ()); + gp_XYZ DNdir = (theD2.XYZ()).Crossed(theDir.XYZ()); + gp_XYZ D2Ndir = (theD3.XYZ()).Crossed(theDir.XYZ()); + double R2 = Ndir.SquareModulus(); + double R = std::sqrt(R2); + double R3 = R2 * R; + double R4 = R2 * R2; + double R5 = R3 * R2; + double Dr = Ndir.Dot(DNdir); + double D2r = Ndir.Dot(D2Ndir) + DNdir.Dot(DNdir); + + if (R5 <= gp::Resolution()) + { + if (R4 <= gp::Resolution()) + { + return false; + } + // We try another computation but the stability is not very good + // dixit ISG. + // V2 = P" (U) : + R4 = R2 * R2; + D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R2)); + D2Ndir.Add(Ndir.Multiplied(((3.0 * Dr * Dr) / R4) - (D2r / R2))); + D2Ndir.Multiply(theOffset / R); + + // V1 = P' (U) : + DNdir.Multiply(R); + DNdir.Subtract(Ndir.Multiplied(Dr / R)); + DNdir.Multiply(theOffset / R2); + } + else + { + // Same computation as IICURV in EUCLID-IS because the stability is better. + // V2 = P" (U) : + D2Ndir.Multiply(theOffset / R); + D2Ndir.Subtract(DNdir.Multiplied(2.0 * theOffset * Dr / R3)); + D2Ndir.Add(Ndir.Multiplied(theOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3)))); + + // V1 = P' (U) : + DNdir.Multiply(theOffset / R); + DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3)); + } + + Ndir.Multiply(theOffset / R); + // P(u) + theValue.ChangeCoord().Add(Ndir); + // P'(u) : + theD1.Add(gp_Vec(DNdir)); + // P"(u) : + if (theIsDirChange) + { + theD2.Reverse(); + } + theD2.Add(gp_Vec(D2Ndir)); + return true; +} + +//! Calculates D0, D1, D2, D3 for 3D offset curve. +//! @param[in,out] theValue on input: basis curve point; on output: offset point +//! @param[in,out] theD1 on input: first derivative of basis; on output: offset curve D1 +//! @param[in,out] theD2 on input: second derivative of basis; on output: offset curve D2 +//! @param[in,out] theD3 on input: third derivative of basis; on output: offset curve D3 +//! @param[in] theD4 fourth derivative of basis curve +//! @param[in] theDir offset reference direction +//! @param[in] theOffset offset distance value +//! @param[in] theIsDirChange flag indicating direction change at singular point +//! @return true if successful, false if computation failed +inline bool CalculateD3(gp_Pnt& theValue, + gp_Vec& theD1, + gp_Vec& theD2, + gp_Vec& theD3, + const gp_Vec& theD4, + const gp_Dir& theDir, + double theOffset, + bool theIsDirChange) +{ + // P(u) = p(u) + Offset * Ndir / R + // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction) + + // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) + + // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) + + // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2))) + + // P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir - + // (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir - + // (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir + + // (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir + + gp_XYZ Ndir = (theD1.XYZ()).Crossed(theDir.XYZ()); + gp_XYZ DNdir = (theD2.XYZ()).Crossed(theDir.XYZ()); + gp_XYZ D2Ndir = (theD3.XYZ()).Crossed(theDir.XYZ()); + gp_XYZ D3Ndir = (theD4.XYZ()).Crossed(theDir.XYZ()); + const double R2 = Ndir.SquareModulus(); + const double R = std::sqrt(R2); + const double R3 = R2 * R; + double R4 = R2 * R2; + const double R5 = R3 * R2; + const double R6 = R3 * R3; + const double R7 = R5 * R2; + const double Dr = Ndir.Dot(DNdir); + const double D2r = Ndir.Dot(D2Ndir) + DNdir.Dot(DNdir); + const double D3r = Ndir.Dot(D3Ndir) + 3.0 * DNdir.Dot(D2Ndir); + if (R7 <= gp::Resolution()) + { + if (R6 <= gp::Resolution()) + { + return false; + } + // V3 = P"' (U) : + D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * Dr / R2)); + D3Ndir.Subtract(DNdir.Multiplied(3.0 * ((D2r / R2) + (Dr * Dr / R4)))); + D3Ndir.Add( + Ndir.Multiplied(6.0 * Dr * Dr / R4 + 6.0 * Dr * D2r / R4 - 15.0 * Dr * Dr * Dr / R6 - D3r)); + D3Ndir.Multiply(theOffset / R); + // V2 = P" (U) : + R4 = R2 * R2; + D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R2)); + D2Ndir.Subtract(Ndir.Multiplied((3.0 * Dr * Dr / R4) - (D2r / R2))); + D2Ndir.Multiply(theOffset / R); + // V1 = P' (U) : + DNdir.Multiply(R); + DNdir.Subtract(Ndir.Multiplied(Dr / R)); + DNdir.Multiply(theOffset / R2); + } + else + { + // V3 = P"' (U) : + D3Ndir.Divide(R); + D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * Dr / R3)); + D3Ndir.Subtract(DNdir.Multiplied((3.0 * ((D2r / R3) + (Dr * Dr) / R5)))); + D3Ndir.Add( + Ndir.Multiplied(6.0 * Dr * Dr / R5 + 6.0 * Dr * D2r / R5 - 15.0 * Dr * Dr * Dr / R7 - D3r)); + D3Ndir.Multiply(theOffset); + // V2 = P" (U) : + D2Ndir.Divide(R); + D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R3)); + D2Ndir.Subtract(Ndir.Multiplied((3.0 * Dr * Dr / R5) - (D2r / R3))); + D2Ndir.Multiply(theOffset); + // V1 = P' (U) : + DNdir.Multiply(theOffset / R); + DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3)); + } + + Ndir.Multiply(theOffset / R); + // P(u) + theValue.ChangeCoord().Add(Ndir); + // P'(u) : + theD1.Add(gp_Vec(DNdir)); + // P"(u) + theD2.Add(gp_Vec(D2Ndir)); + // P"'(u) + if (theIsDirChange) + { + theD3.Reverse(); + } + theD3.Add(gp_Vec(D3Ndir)); + return true; +} + +//! Adjusts derivatives at singular points where the first derivative is nearly zero. +//! Uses Taylor series approximation to find a valid tangent direction. +//! @tparam CurveType type supporting D0, DN methods (Geom_Curve or adaptor) +//! @param[in] theCurve basis curve for derivative evaluation +//! @param[in] theMaxDerivative maximum derivative order to compute (3 or 4) +//! @param[in] theU parameter value +//! @param[in,out] theD1 first derivative (will be adjusted) +//! @param[in,out] theD2 second derivative (will be adjusted) +//! @param[in,out] theD3 third derivative (will be adjusted) +//! @param[in,out] theD4 fourth derivative (will be adjusted if theMaxDerivative >= 4) +//! @return true if direction change detected +template +bool AdjustDerivative(const CurveType& theCurve, + int theMaxDerivative, + double theU, + gp_Vec& theD1, + gp_Vec& theD2, + gp_Vec& theD3, + gp_Vec& theD4) +{ + static const double aTol = gp::Resolution(); + static const double aMinStep = 1e-7; + static const int aMaxDerivOrder = 3; + + bool isDirectionChange = false; + const double anUinfium = theCurve.FirstParameter(); + const double anUsupremum = theCurve.LastParameter(); + + static const double DivisionFactor = 1.e-3; + double du; + if ((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + { + du = 0.0; + } + else + { + du = anUsupremum - anUinfium; + } + + const double aDelta = std::max(du * DivisionFactor, aMinStep); + + // Derivative is approximated by Taylor-series + int anIndex = 1; // Derivative order + gp_Vec V; + + do + { + V = theCurve.DN(theU, ++anIndex); + } while ((V.SquareMagnitude() <= aTol) && anIndex < aMaxDerivOrder); + + double u; + + if (theU - anUinfium < aDelta) + { + u = theU + aDelta; + } + else + { + u = theU - aDelta; + } + + gp_Pnt P1, P2; + theCurve.D0(std::min(theU, u), P1); + theCurve.D0(std::max(theU, u), P2); + + gp_Vec V1(P1, P2); + isDirectionChange = V.Dot(V1) < 0.0; + const double aSign = isDirectionChange ? -1.0 : 1.0; + + theD1 = V * aSign; + gp_Vec* aDeriv[3] = {&theD2, &theD3, &theD4}; + for (int i = 1; i < theMaxDerivative; i++) + { + *(aDeriv[i - 1]) = theCurve.DN(theU, anIndex + i) * aSign; + } + + return isDirectionChange; +} + +//! Template function for D0 evaluation of offset curve. +//! Gets D1 from basis curve and computes offset point. +//! +//! @tparam BasisCurveType type of basis curve (must have D1 method) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theDir offset reference direction +//! @param[in] theOffset offset distance +//! @param[out] theValue computed offset point +//! @return true if evaluation succeeded, false if normal is undefined +template +bool EvaluateD0(double theU, + const BasisCurveType& theBasisCurve, + const gp_Dir& theDir, + double theOffset, + gp_Pnt& theValue) +{ + gp_Vec aD1; + theBasisCurve->D1(theU, theValue, aD1); + return CalculateD0(theValue, aD1, theDir, theOffset); +} + +//! Template function for D1 evaluation of offset curve. +//! Gets D2 from basis curve and computes offset point and first derivative. +//! +//! @tparam BasisCurveType type of basis curve (must have D2 method) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theDir offset reference direction +//! @param[in] theOffset offset distance +//! @param[out] theValue computed offset point +//! @param[out] theD1 computed first derivative +//! @return true if evaluation succeeded, false if computation failed +template +bool EvaluateD1(double theU, + const BasisCurveType& theBasisCurve, + const gp_Dir& theDir, + double theOffset, + gp_Pnt& theValue, + gp_Vec& theD1) +{ + gp_Vec aD2; + theBasisCurve->D2(theU, theValue, theD1, aD2); + return CalculateD1(theValue, theD1, aD2, theDir, theOffset); +} + +//! Template function for D2 evaluation of offset curve. +//! Gets D3 from basis curve and computes offset point and derivatives. +//! Handles singular points where first derivative is nearly zero. +//! +//! @tparam BasisCurveType type of basis curve (must have D3 and DN methods) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theDir offset reference direction +//! @param[in] theOffset offset distance +//! @param[out] theValue computed offset point +//! @param[out] theD1 computed first derivative +//! @param[out] theD2 computed second derivative +//! @return true if evaluation succeeded, false if computation failed +template +bool EvaluateD2(double theU, + const BasisCurveType& theBasisCurve, + const gp_Dir& theDir, + double theOffset, + gp_Pnt& theValue, + gp_Vec& theD1, + gp_Vec& theD2) +{ + gp_Vec aD3; + theBasisCurve->D3(theU, theValue, theD1, theD2, aD3); + + bool isDirectionChange = false; + if (theD1.SquareMagnitude() <= gp::Resolution()) + { + gp_Vec aDummyD4; + isDirectionChange = AdjustDerivative(*theBasisCurve, 3, theU, theD1, theD2, aD3, aDummyD4); + } + + return CalculateD2(theValue, theD1, theD2, aD3, theDir, theOffset, isDirectionChange); +} + +//! Template function for D3 evaluation of offset curve. +//! Gets D3 and D4 from basis curve and computes offset point and derivatives. +//! Handles singular points where first derivative is nearly zero. +//! +//! @tparam BasisCurveType type of basis curve (must have D3 and DN methods) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theDir offset reference direction +//! @param[in] theOffset offset distance +//! @param[out] theValue computed offset point +//! @param[out] theD1 computed first derivative +//! @param[out] theD2 computed second derivative +//! @param[out] theD3 computed third derivative +//! @return true if evaluation succeeded, false if computation failed +template +bool EvaluateD3(double theU, + const BasisCurveType& theBasisCurve, + const gp_Dir& theDir, + double theOffset, + gp_Pnt& theValue, + gp_Vec& theD1, + gp_Vec& theD2, + gp_Vec& theD3) +{ + theBasisCurve->D3(theU, theValue, theD1, theD2, theD3); + gp_Vec aD4 = theBasisCurve->DN(theU, 4); + + bool isDirectionChange = false; + if (theD1.SquareMagnitude() <= gp::Resolution()) + { + isDirectionChange = AdjustDerivative(*theBasisCurve, 4, theU, theD1, theD2, theD3, aD4); + } + + return CalculateD3(theValue, theD1, theD2, theD3, aD4, theDir, theOffset, isDirectionChange); +} + +//! Template function for DN evaluation of offset curve. +//! Handles derivatives up to order 3 using D1/D2/D3 methods. +//! For derivatives > 3, returns the basis curve derivative directly +//! (offset contribution is negligible for high-order derivatives). +//! +//! @tparam BasisCurveType type of basis curve (must have D1, D2, D3, DN methods) +//! @param[in] theU parameter value +//! @param[in] theBasisCurve basis curve +//! @param[in] theDir offset reference direction +//! @param[in] theOffset offset distance +//! @param[in] theN derivative order (must be >= 1) +//! @param[out] theDN computed N-th derivative +//! @return true if evaluation succeeded, false if computation failed +template +bool EvaluateDN(double theU, + const BasisCurveType& theBasisCurve, + const gp_Dir& theDir, + double theOffset, + int theN, + gp_Vec& theDN) +{ + gp_Pnt aPnt; + gp_Vec aDummy; + switch (theN) + { + case 1: + return EvaluateD1(theU, theBasisCurve, theDir, theOffset, aPnt, theDN); + case 2: + return EvaluateD2(theU, theBasisCurve, theDir, theOffset, aPnt, aDummy, theDN); + case 3: + return EvaluateD3(theU, theBasisCurve, theDir, theOffset, aPnt, aDummy, aDummy, theDN); + default: + // For derivatives > 3, return basis curve derivative (no offset contribution) + theDN = theBasisCurve->DN(theU, theN); + return true; + } +} + +} // namespace Geom_OffsetCurveUtils + +#endif // _Geom_OffsetCurveUtils_HeaderFile diff --git a/src/ModelingData/TKG3d/Geom/Geom_OffsetSurfaceUtils.pxx b/src/ModelingData/TKG3d/Geom/Geom_OffsetSurfaceUtils.pxx new file mode 100644 index 0000000000..559fe7a0d8 --- /dev/null +++ b/src/ModelingData/TKG3d/Geom/Geom_OffsetSurfaceUtils.pxx @@ -0,0 +1,1500 @@ +// Copyright (c) 2015-2025 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 _Geom_OffsetSurfaceUtils_HeaderFile +#define _Geom_OffsetSurfaceUtils_HeaderFile + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +//! Internal helper namespace for 3D offset surface calculations. +//! Provides static inline functions to compute offset surface point and derivatives +//! from basis surface derivatives. +//! +//! Includes both non-singular (simple cross product normal) and singular +//! (osculating surface) case handling. +//! +//! Mathematical basis: +//! P(u,v) = p(u,v) + Offset * N / ||N|| +//! where N = dP/du ^ dP/dv is the surface normal +namespace Geom_OffsetSurfaceUtils +{ + +//! Default tolerance for normal magnitude check +constexpr double THE_D1_MAGNITUDE_TOL = 1.e-9; + +//! Struct to hold osculating surface query results. +//! Used to abstract osculating surface handling between different classes. +struct OsculatingInfo +{ + bool AlongU = false; //!< True if osculating along U direction + bool AlongV = false; //!< True if osculating along V direction + bool IsOpposite = false; //!< True if normal direction should be reversed + + //! Returns the sign factor for offset calculation + double Sign() const { return ((AlongU || AlongV) && IsOpposite) ? -1.0 : 1.0; } + + //! Returns true if osculating surface is available + bool HasOsculating() const { return AlongU || AlongV; } +}; + +//! Checks if a vector has infinite coordinates. +//! @param[in] theVec vector to check +//! @return true if any coordinate is infinite +inline bool IsInfiniteCoord(const gp_Vec& theVec) +{ + return Precision::IsInfinite(theVec.X()) || Precision::IsInfinite(theVec.Y()) + || Precision::IsInfinite(theVec.Z()); +} + +//! Checks if surface normal is singular (has zero magnitude). +//! @param[in] theD1U first derivative with respect to U +//! @param[in] theD1V first derivative with respect to V +//! @param[in] theTol tolerance for magnitude check +//! @return true if normal magnitude is below tolerance (singular case) +inline bool IsSingular(const gp_Vec& theD1U, + const gp_Vec& theD1V, + double theTol = THE_D1_MAGNITUDE_TOL) +{ + // Normalize derivatives before normal calculation for stability + gp_Vec aD1U(theD1U); + gp_Vec aD1V(theD1V); + double aD1UNorm2 = aD1U.SquareMagnitude(); + double aD1VNorm2 = aD1V.SquareMagnitude(); + if (aD1UNorm2 > 1.0) + aD1U /= std::sqrt(aD1UNorm2); + if (aD1VNorm2 > 1.0) + aD1V /= std::sqrt(aD1VNorm2); + + gp_Vec aNorm = aD1U.Crossed(aD1V); + return aNorm.SquareMagnitude() <= theTol * theTol; +} + +//! Calculates normalized normal vector for non-singular case. +//! @param[in] theD1U first derivative with respect to U +//! @param[in] theD1V first derivative with respect to V +//! @param[out] theNormal computed normalized normal (valid only if return is true) +//! @param[in] theTol tolerance for magnitude check +//! @return true if normal computed successfully, false if singular +inline bool ComputeNormal(const gp_Vec& theD1U, + const gp_Vec& theD1V, + gp_Vec& theNormal, + double theTol = THE_D1_MAGNITUDE_TOL) +{ + // Normalize derivatives before normal calculation for stability + gp_Vec aD1U(theD1U); + gp_Vec aD1V(theD1V); + double aD1UNorm2 = aD1U.SquareMagnitude(); + double aD1VNorm2 = aD1V.SquareMagnitude(); + if (aD1UNorm2 > 1.0) + aD1U /= std::sqrt(aD1UNorm2); + if (aD1VNorm2 > 1.0) + aD1V /= std::sqrt(aD1VNorm2); + + theNormal = aD1U.Crossed(aD1V); + if (theNormal.SquareMagnitude() <= theTol * theTol) + { + return false; + } + theNormal.Normalize(); + return true; +} + +//! Computes dN/du for non-singular offset surface. +//! @param[in] theD1U first derivative with respect to U +//! @param[in] theD1V first derivative with respect to V +//! @param[in] theD2U second derivative d2P/du2 +//! @param[in] theD2UV mixed derivative d2P/dudv +//! @param[in] theNormal unit normal vector +//! @return derivative of normal with respect to U +inline gp_Vec ComputeDNormalU(const gp_Vec& theD1U, + const gp_Vec& theD1V, + const gp_Vec& theD2U, + const gp_Vec& theD2UV, + const gp_Vec& theNormal) +{ + double aScale = (theD1U ^ theD1V).Dot(theNormal); + + gp_Vec aN1U; + aN1U.SetX(theD2U.Y() * theD1V.Z() + theD1U.Y() * theD2UV.Z() - theD2U.Z() * theD1V.Y() + - theD1U.Z() * theD2UV.Y()); + aN1U.SetY(-(theD2U.X() * theD1V.Z() + theD1U.X() * theD2UV.Z() - theD2U.Z() * theD1V.X() + - theD1U.Z() * theD2UV.X())); + aN1U.SetZ(theD2U.X() * theD1V.Y() + theD1U.X() * theD2UV.Y() - theD2U.Y() * theD1V.X() + - theD1U.Y() * theD2UV.X()); + double aScaleU = aN1U.Dot(theNormal); + aN1U.Subtract(aScaleU * theNormal); + aN1U /= aScale; + + return aN1U; +} + +//! Computes dN/dv for non-singular offset surface. +//! @param[in] theD1U first derivative with respect to U +//! @param[in] theD1V first derivative with respect to V +//! @param[in] theD2V second derivative d2P/dv2 +//! @param[in] theD2UV mixed derivative d2P/dudv +//! @param[in] theNormal unit normal vector +//! @return derivative of normal with respect to V +inline gp_Vec ComputeDNormalV(const gp_Vec& theD1U, + const gp_Vec& theD1V, + const gp_Vec& theD2V, + const gp_Vec& theD2UV, + const gp_Vec& theNormal) +{ + double aScale = (theD1U ^ theD1V).Dot(theNormal); + + gp_Vec aN1V; + aN1V.SetX(theD2UV.Y() * theD1V.Z() + theD2V.Z() * theD1U.Y() - theD2UV.Z() * theD1V.Y() + - theD2V.Y() * theD1U.Z()); + aN1V.SetY(-(theD2UV.X() * theD1V.Z() + theD2V.Z() * theD1U.X() - theD2UV.Z() * theD1V.X() + - theD2V.X() * theD1U.Z())); + aN1V.SetZ(theD2UV.X() * theD1V.Y() + theD2V.Y() * theD1U.X() - theD2UV.Y() * theD1V.X() + - theD2V.X() * theD1U.Y()); + double aScaleV = aN1V.Dot(theNormal); + aN1V.Subtract(aScaleV * theNormal); + aN1V /= aScale; + + return aN1V; +} + +//! Calculates D0 (point) for offset surface in non-singular case. +//! @param[in,out] theValue on input: basis surface point; on output: offset point +//! @param[in] theD1U first derivative with respect to U +//! @param[in] theD1V first derivative with respect to V +//! @param[in] theOffset offset distance value +//! @param[in] theSign sign factor (1.0 or -1.0) for offset direction +//! @return false if singular (normal has zero magnitude), true otherwise +inline bool CalculateD0(gp_Pnt& theValue, + const gp_Vec& theD1U, + const gp_Vec& theD1V, + double theOffset, + double theSign = 1.0) +{ + gp_Vec aNorm; + if (!ComputeNormal(theD1U, theD1V, aNorm)) + { + return false; + } + theValue.SetXYZ(theValue.XYZ() + theOffset * theSign * aNorm.XYZ()); + return true; +} + +//! Calculates D0 and D1 for offset surface in non-singular case. +//! @param[in,out] theValue on input: basis surface point; on output: offset point +//! @param[in,out] theD1U on input: basis dP/du; on output: offset surface dP/du +//! @param[in,out] theD1V on input: basis dP/dv; on output: offset surface dP/dv +//! @param[in] theD2U second derivative d2P/du2 of basis surface +//! @param[in] theD2V second derivative d2P/dv2 of basis surface +//! @param[in] theD2UV mixed derivative d2P/dudv of basis surface +//! @param[in] theOffset offset distance value +//! @param[in] theSign sign factor (1.0 or -1.0) for offset direction +//! @return false if singular (normal has zero magnitude), true otherwise +inline bool CalculateD1(gp_Pnt& theValue, + gp_Vec& theD1U, + gp_Vec& theD1V, + const gp_Vec& theD2U, + const gp_Vec& theD2V, + const gp_Vec& theD2UV, + double theOffset, + double theSign = 1.0) +{ + gp_Vec aNorm; + if (!ComputeNormal(theD1U, theD1V, aNorm)) + { + return false; + } + + // Compute offset point + theValue.SetXYZ(theValue.XYZ() + theOffset * theSign * aNorm.XYZ()); + + // Compute normal derivatives + gp_Vec aN1U = ComputeDNormalU(theD1U, theD1V, theD2U, theD2UV, aNorm); + gp_Vec aN1V = ComputeDNormalV(theD1U, theD1V, theD2V, theD2UV, aNorm); + + theD1U += theOffset * theSign * aN1U; + theD1V += theOffset * theSign * aN1V; + + return true; +} + +//! Template function for computing derivatives at singular points. +//! Works with any surface type that provides D1, D2, D3, DN methods. +//! @tparam BasisSurfType type of basis surface (Handle to surface or adaptor) +//! @tparam OscSurfType type of osculating surface (Handle to BSpline surface) +//! @param[in] theMaxOrder maximum derivative order +//! @param[in] theMinOrder minimum derivative order +//! @param[in] theU U parameter +//! @param[in] theV V parameter +//! @param[in] theBasisSurf basis surface +//! @param[in] theNU derivative order in U for output +//! @param[in] theNV derivative order in V for output +//! @param[in] theAlongU true if osculating along U +//! @param[in] theAlongV true if osculating along V +//! @param[in] theOscSurf osculating surface (may be null) +//! @param[out] theDerNUV array of normal derivatives +//! @param[in,out] theDerSurf array of surface derivatives +template +void ComputeDerivatives(int theMaxOrder, + int theMinOrder, + double theU, + double theV, + const BasisSurfType& theBasisSurf, + int theNU, + int theNV, + bool theAlongU, + bool theAlongV, + const OscSurfType& theOscSurf, + TColgp_Array2OfVec& theDerNUV, + TColgp_Array2OfVec& theDerSurf) +{ + gp_Pnt P; + gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V; + + if (theAlongU || theAlongV) + { + theMaxOrder = 0; + // Stack buffer for DerSurfL: max size is (theNU + 2) x (theNV + 2) + // DN can have theNU/theNV up to 6, so max is 8x8=64 + const int aDerSurfLSize = (theMaxOrder + theNU + 2) * (theMaxOrder + theNV + 2); + NCollection_LocalArray aDerSurfLBuffer(aDerSurfLSize); + NCollection_Array2 DerSurfL(aDerSurfLBuffer[0], + 0, + theMaxOrder + theNU + 1, + 0, + theMaxOrder + theNV + 1); + switch (theMinOrder) + { + case 1: + theOscSurf->D1(theU, theV, P, DL1U, DL1V); + DerSurfL.SetValue(1, 0, DL1U); + DerSurfL.SetValue(0, 1, DL1V); + break; + case 2: + theOscSurf->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV); + DerSurfL.SetValue(1, 0, DL1U); + DerSurfL.SetValue(0, 1, DL1V); + DerSurfL.SetValue(1, 1, DL2UV); + DerSurfL.SetValue(2, 0, DL2U); + DerSurfL.SetValue(0, 2, DL2V); + break; + case 3: + theOscSurf->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV); + DerSurfL.SetValue(1, 0, DL1U); + DerSurfL.SetValue(0, 1, DL1V); + DerSurfL.SetValue(1, 1, DL2UV); + DerSurfL.SetValue(2, 0, DL2U); + DerSurfL.SetValue(0, 2, DL2V); + DerSurfL.SetValue(3, 0, DL3U); + DerSurfL.SetValue(2, 1, DL3UUV); + DerSurfL.SetValue(1, 2, DL3UVV); + DerSurfL.SetValue(0, 3, DL3V); + break; + default: + break; + } + + if (theNU <= theNV) + { + for (int i = 0; i <= theMaxOrder + 1 + theNU; i++) + for (int j = i; j <= theMaxOrder + theNV + 1; j++) + if (i + j > theMinOrder) + { + DerSurfL.SetValue(i, j, theOscSurf->DN(theU, theV, i, j)); + theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j)); + if (i != j && j <= theNU + 1) + { + theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i)); + DerSurfL.SetValue(j, i, theOscSurf->DN(theU, theV, j, i)); + } + } + } + else + { + for (int j = 0; j <= theMaxOrder + 1 + theNV; j++) + for (int i = j; i <= theMaxOrder + theNU + 1; i++) + if (i + j > theMinOrder) + { + DerSurfL.SetValue(i, j, theOscSurf->DN(theU, theV, i, j)); + theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j)); + if (i != j && i <= theNV + 1) + { + theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i)); + DerSurfL.SetValue(j, i, theOscSurf->DN(theU, theV, j, i)); + } + } + } + for (int i = 0; i <= theMaxOrder + theNU; i++) + for (int j = 0; j <= theMaxOrder + theNV; j++) + { + if (theAlongU) + theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf)); + if (theAlongV) + theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL)); + } + } + else + { + for (int i = 0; i <= theMaxOrder + theNU + 1; i++) + { + for (int j = i; j <= theMaxOrder + theNV + 1; j++) + { + if (i + j > theMinOrder) + { + theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j)); + if (i != j && j <= theDerSurf.UpperRow() && i <= theDerSurf.UpperCol()) + { + theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i)); + } + } + } + } + for (int i = 0; i <= theMaxOrder + theNU; i++) + for (int j = 0; j <= theMaxOrder + theNV; j++) + theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf)); + } +} + +//! Attempts to replace a zero derivative by stepping away and recomputing. +//! This handles CSLib_InfinityOfSolutions case where one derivative is zero. +//! +//! @tparam BasisSurfType type of basis surface (must have D1 method) +//! @param[in] theU U parameter +//! @param[in] theV V parameter +//! @param[in] theUMin minimum U bound +//! @param[in] theUMax maximum U bound +//! @param[in] theVMin minimum V bound +//! @param[in] theVMax maximum V bound +//! @param[in,out] theDU first derivative with respect to U +//! @param[in,out] theDV first derivative with respect to V +//! @param[in] theSquareTol squared tolerance for zero check +//! @param[in] theBasisSurf basis surface for computing derivatives +//! @return true if derivative was successfully replaced +template +bool ReplaceDerivative(double theU, + double theV, + double theUMin, + double theUMax, + double theVMin, + double theVMax, + gp_Vec& theDU, + gp_Vec& theDV, + double theSquareTol, + const BasisSurfType& theBasisSurf) +{ + bool isReplaceDU = theDU.SquareMagnitude() < theSquareTol; + bool isReplaceDV = theDV.SquareMagnitude() < theSquareTol; + bool isReplaced = false; + + // Only handle case where exactly one derivative is zero + if (isReplaceDU != isReplaceDV) + { + // Calculate step along non-zero derivative + double aStep; + if (isReplaceDV) + { + aStep = Precision::Confusion() * theDU.Magnitude(); + if (aStep > theUMax - theUMin) + aStep = (theUMax - theUMin) / 100.; + } + else + { + aStep = Precision::Confusion() * theDV.Magnitude(); + if (aStep > theVMax - theVMin) + aStep = (theVMax - theVMin) / 100.; + } + + gp_Pnt aP; + gp_Vec aDU, aDV; + + // Step away from current parametric coordinates and calculate derivatives once again. + // Replace zero derivative by the obtained. + for (double aStepSign = -1.0; aStepSign <= 1.0 && !isReplaced; aStepSign += 2.0) + { + double aU = theU; + double aV = theV; + + if (isReplaceDV) + { + aU = theU + aStepSign * aStep; + if (aU < theUMin || aU > theUMax) + continue; + } + else + { + aV = theV + aStepSign * aStep; + if (aV < theVMin || aV > theVMax) + continue; + } + + theBasisSurf->D1(aU, aV, aP, aDU, aDV); + + if (isReplaceDU && aDU.SquareMagnitude() > theSquareTol) + { + theDU = aDU; + isReplaced = true; + } + if (isReplaceDV && aDV.SquareMagnitude() > theSquareTol) + { + theDV = aDV; + isReplaced = true; + } + } + } + return isReplaced; +} + +//! Attempts to shift the evaluation point towards the center of the parametric space. +//! This is used when normal calculation fails at singular points near boundaries. +//! The point is shifted iteratively, each time doubling the distance from the start point. +//! +//! @param[in] theUStart original U parameter (for direction calculation) +//! @param[in] theVStart original V parameter (for direction calculation) +//! @param[in,out] theU current U parameter, modified on success +//! @param[in,out] theV current V parameter, modified on success +//! @param[in] theUMin minimum U bound +//! @param[in] theUMax maximum U bound +//! @param[in] theVMin minimum V bound +//! @param[in] theVMax maximum V bound +//! @param[in] theIsUPeriodic true if surface is U-periodic +//! @param[in] theIsVPeriodic true if surface is V-periodic +//! @param[in] theD1U first derivative with respect to U (for singularity check) +//! @param[in] theD1V first derivative with respect to V (for singularity check) +//! @return true if shift was successful, false if center is reached +inline bool ShiftPoint(double theUStart, + double theVStart, + double& theU, + double& theV, + double theUMin, + double theUMax, + double theVMin, + double theVMax, + bool theIsUPeriodic, + bool theIsVPeriodic, + const gp_Vec& theD1U, + const gp_Vec& theD1V) +{ + // Check if either U or V is singular (normally one of them is) + bool isUSingular = (theD1U.SquareMagnitude() < THE_D1_MAGNITUDE_TOL * THE_D1_MAGNITUDE_TOL); + bool isVSingular = (theD1V.SquareMagnitude() < THE_D1_MAGNITUDE_TOL * THE_D1_MAGNITUDE_TOL); + + // Compute vector to shift from start point to center of the surface; + // if surface is periodic or singular in some direction, take shift in that direction zero + double aDirU = + (theIsUPeriodic || (isUSingular && !isVSingular) ? 0. : 0.5 * (theUMin + theUMax) - theUStart); + double aDirV = + (theIsVPeriodic || (isVSingular && !isUSingular) ? 0. : 0.5 * (theVMin + theVMax) - theVStart); + double aDist = std::sqrt(aDirU * aDirU + aDirV * aDirV); + + // Shift current point from its current position towards center, by value of twice + // current distance from it to start (but not less than Precision::PConfusion()); + // fail if center is overpassed. + double aDU = theU - theUStart; + double aDV = theV - theVStart; + double aStep = std::max(2. * std::sqrt(aDU * aDU + aDV * aDV), Precision::PConfusion()); + if (aStep >= aDist) + { + return false; + } + + aStep /= aDist; + theU += aDirU * aStep; + theV += aDirV * aStep; + + return true; +} + +//! Template function for D0 evaluation with retry mechanism for singular points. +//! When normal calculation fails, attempts to shift the point towards the center +//! and retry the calculation. +//! +//! @tparam BasisSurfType type of basis surface (must have D1 method) +//! @tparam OscSurfQueryType type providing osculating surface query +//! @param[in] theU U parameter +//! @param[in] theV V parameter +//! @param[in] theBasisSurf basis surface adaptor +//! @param[in] theOffset offset distance +//! @param[in] theOscQuery osculating surface query object (may be null) +//! @param[out] theValue computed offset point +//! @return true if calculation succeeded, false if failed at singular point +template +bool EvaluateD0(double theU, + double theV, + const BasisSurfType& theBasisSurf, + double theOffset, + const OscSurfQueryType& theOscQuery, + gp_Pnt& theValue) +{ + const double aUStart = theU; + const double aVStart = theV; + double aUMin, aUMax, aVMin, aVMax; + theBasisSurf->Bounds(aUMin, aUMax, aVMin, aVMax); + const bool isUPer = theBasisSurf->IsUPeriodic(); + const bool isVPer = theBasisSurf->IsVPeriodic(); + + for (;;) + { + gp_Vec aD1U, aD1V; + theBasisSurf->D1(theU, theV, theValue, aD1U, aD1V); + + if (IsInfiniteCoord(aD1U) || IsInfiniteCoord(aD1V)) + { + return false; + } + + // Try non-singular case first + if (CalculateD0(theValue, aD1U, aD1V, theOffset)) + { + return true; + } + + // Singular case - query osculating surface and use higher order derivatives + constexpr int aMaxOrder = 3; + + OsculatingInfo aOscInfo; + Handle(Geom_BSplineSurface) aOscSurf; + if (theOscQuery) + { + aOscInfo.AlongU = theOscQuery->UOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + aOscInfo.AlongV = theOscQuery->VOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + } + + constexpr int aDerNUVSize = (aMaxOrder + 1) * (aMaxOrder + 1); // 16 + constexpr int aDerSurfSize = (aMaxOrder + 2) * (aMaxOrder + 2); // 25 + NCollection_LocalArray aDerNUVBuffer(aDerNUVSize); + NCollection_LocalArray aDerSurfBuffer(aDerSurfSize); + NCollection_Array2 aDerNUV(aDerNUVBuffer[0], 0, aMaxOrder, 0, aMaxOrder); + NCollection_Array2 aDerSurf(aDerSurfBuffer[0], 0, aMaxOrder + 1, 0, aMaxOrder + 1); + + aDerSurf.SetValue(1, 0, aD1U); + aDerSurf.SetValue(0, 1, aD1V); + + // Use ComputeDerivatives which handles osculating surface properly + if (aOscInfo.HasOsculating() && !aOscSurf.IsNull()) + { + ComputeDerivatives(aMaxOrder, + 1, + theU, + theV, + theBasisSurf, + 0, + 0, + aOscInfo.AlongU, + aOscInfo.AlongV, + aOscSurf, + aDerNUV, + aDerSurf); + } + else + { + Handle(Geom_BSplineSurface) aDummy; + ComputeDerivatives(aMaxOrder, + 1, + theU, + theV, + theBasisSurf, + 0, + 0, + false, + false, + aDummy, + aDerNUV, + aDerSurf); + } + + gp_Dir aNormal; + CSLib_NormalStatus aNStatus; + int OrderU, OrderV; + CSLib::Normal(aMaxOrder, + aDerNUV, + THE_D1_MAGNITUDE_TOL, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + aNStatus, + aNormal, + OrderU, + OrderV); + + // Handle CSLib_InfinityOfSolutions by replacing zero derivative + if (aNStatus == CSLib_InfinityOfSolutions) + { + gp_Vec aNewDU = aD1U; + gp_Vec aNewDV = aD1V; + if (ReplaceDerivative(theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + aNewDU, + aNewDV, + THE_D1_MAGNITUDE_TOL * THE_D1_MAGNITUDE_TOL, + theBasisSurf)) + { + CSLib::Normal(aNewDU, aNewDV, THE_D1_MAGNITUDE_TOL, aNStatus, aNormal); + } + } + + if (aNStatus == CSLib_Defined) + { + theValue.SetXYZ(theValue.XYZ() + theOffset * aOscInfo.Sign() * aNormal.XYZ()); + return true; + } + + // Try shifting point towards center - returns false when center is overpassed + if (!ShiftPoint(aUStart, + aVStart, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + isUPer, + isVPer, + aD1U, + aD1V)) + { + return false; + } + } +} + +//! Template function for D1 evaluation with retry mechanism for singular points. +//! +//! @tparam BasisSurfType type of basis surface (must have D2 method) +//! @tparam OscSurfQueryType type providing osculating surface query +//! @param[in] theU U parameter +//! @param[in] theV V parameter +//! @param[in] theBasisSurf basis surface adaptor +//! @param[in] theOffset offset distance +//! @param[in] theOscQuery osculating surface query object (may be null) +//! @param[out] theValue computed offset point +//! @param[out] theD1U computed D1U derivative +//! @param[out] theD1V computed D1V derivative +//! @return true if calculation succeeded, false if failed at singular point +template +bool EvaluateD1(double theU, + double theV, + const BasisSurfType& theBasisSurf, + double theOffset, + const OscSurfQueryType& theOscQuery, + gp_Pnt& theValue, + gp_Vec& theD1U, + gp_Vec& theD1V) +{ + const double aUStart = theU; + const double aVStart = theV; + double aUMin, aUMax, aVMin, aVMax; + theBasisSurf->Bounds(aUMin, aUMax, aVMin, aVMax); + const bool isUPer = theBasisSurf->IsUPeriodic(); + const bool isVPer = theBasisSurf->IsVPeriodic(); + + for (;;) + { + gp_Vec aD2U, aD2V, aD2UV; + theBasisSurf->D2(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV); + + if (IsInfiniteCoord(theD1U) || IsInfiniteCoord(theD1V)) + { + return false; + } + + // Check if singular by normalizing derivatives and computing cross product + gp_Vec aD1U(theD1U); + gp_Vec aD1V(theD1V); + double aD1UNorm2 = aD1U.SquareMagnitude(); + double aD1VNorm2 = aD1V.SquareMagnitude(); + if (aD1UNorm2 > 1.0) + aD1U /= std::sqrt(aD1UNorm2); + if (aD1VNorm2 > 1.0) + aD1V /= std::sqrt(aD1VNorm2); + + bool isSingular = false; + const int aMaxOrder = 3; + gp_Vec aNorm = aD1U.Crossed(aD1V); + + // Query osculating surface only if singular + OsculatingInfo aOscInfo; + Handle(Geom_BSplineSurface) aOscSurf; + if (aNorm.SquareMagnitude() <= THE_D1_MAGNITUDE_TOL * THE_D1_MAGNITUDE_TOL) + { + if (theOscQuery) + { + aOscInfo.AlongU = theOscQuery->UOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + aOscInfo.AlongV = theOscQuery->VOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + } + isSingular = true; + } + + // Compute sign factor + const double aSign = aOscInfo.Sign(); + + // Non-singular case: use direct formulas + if (!isSingular) + { + aNorm.Normalize(); + theValue.SetXYZ(theValue.XYZ() + theOffset * aSign * aNorm.XYZ()); + + // Compute normal derivatives using inline formulas + gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V; + const double aScale = (theD1U ^ theD1V).Dot(aN0); + aN1U.SetX(aD2U.Y() * theD1V.Z() + theD1U.Y() * aD2UV.Z() - aD2U.Z() * theD1V.Y() + - theD1U.Z() * aD2UV.Y()); + aN1U.SetY((aD2U.X() * theD1V.Z() + theD1U.X() * aD2UV.Z() - aD2U.Z() * theD1V.X() + - theD1U.Z() * aD2UV.X()) + * -1.0); + aN1U.SetZ(aD2U.X() * theD1V.Y() + theD1U.X() * aD2UV.Y() - aD2U.Y() * theD1V.X() + - theD1U.Y() * aD2UV.X()); + const double aScaleU = aN1U.Dot(aN0); + aN1U.Subtract(aScaleU * aN0); + aN1U /= aScale; + + aN1V.SetX(aD2UV.Y() * theD1V.Z() + aD2V.Z() * theD1U.Y() - aD2UV.Z() * theD1V.Y() + - aD2V.Y() * theD1U.Z()); + aN1V.SetY((aD2UV.X() * theD1V.Z() + aD2V.Z() * theD1U.X() - aD2UV.Z() * theD1V.X() + - aD2V.X() * theD1U.Z()) + * -1.0); + aN1V.SetZ(aD2UV.X() * theD1V.Y() + aD2V.Y() * theD1U.X() - aD2UV.Y() * theD1V.X() + - aD2V.X() * theD1U.Y()); + const double aScaleV = aN1V.Dot(aN0); + aN1V.Subtract(aScaleV * aN0); + aN1V /= aScale; + + theD1U += theOffset * aSign * aN1U; + theD1V += theOffset * aSign * aN1V; + + return true; + } + + // Singular case - use higher order derivatives + constexpr int aDerNUVSize = (aMaxOrder + 2) * (aMaxOrder + 2); // 25 + constexpr int aDerSurfSize = (aMaxOrder + 3) * (aMaxOrder + 3); // 36 + NCollection_LocalArray aDerNUVBuffer(aDerNUVSize); + NCollection_LocalArray aDerSurfBuffer(aDerSurfSize); + NCollection_Array2 aDerNUV(aDerNUVBuffer[0], 0, aMaxOrder + 1, 0, aMaxOrder + 1); + NCollection_Array2 aDerSurf(aDerSurfBuffer[0], 0, aMaxOrder + 2, 0, aMaxOrder + 2); + + aDerSurf.SetValue(1, 0, theD1U); + aDerSurf.SetValue(0, 1, theD1V); + aDerSurf.SetValue(1, 1, aD2UV); + aDerSurf.SetValue(2, 0, aD2U); + aDerSurf.SetValue(0, 2, aD2V); + + if (aOscInfo.HasOsculating() && !aOscSurf.IsNull()) + { + ComputeDerivatives(aMaxOrder, + 2, + theU, + theV, + theBasisSurf, + 1, + 1, + aOscInfo.AlongU, + aOscInfo.AlongV, + aOscSurf, + aDerNUV, + aDerSurf); + } + else + { + Handle(Geom_BSplineSurface) aDummy; + ComputeDerivatives(aMaxOrder, + 2, + theU, + theV, + theBasisSurf, + 1, + 1, + false, + false, + aDummy, + aDerNUV, + aDerSurf); + } + + gp_Dir aNormal; + CSLib_NormalStatus aNStatus; + int aOrderU, aOrderV; + CSLib::Normal(aMaxOrder, + aDerNUV, + THE_D1_MAGNITUDE_TOL, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + aNStatus, + aNormal, + aOrderU, + aOrderV); + + // Handle CSLib_InfinityOfSolutions by replacing zero derivative + if (aNStatus == CSLib_InfinityOfSolutions) + { + gp_Vec aNewDU = theD1U; + gp_Vec aNewDV = theD1V; + if (ReplaceDerivative(theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + aNewDU, + aNewDV, + THE_D1_MAGNITUDE_TOL * THE_D1_MAGNITUDE_TOL, + theBasisSurf)) + { + // Re-compute with replaced derivatives + aDerSurf.SetValue(1, 0, aNewDU); + aDerSurf.SetValue(0, 1, aNewDV); + if (aOscInfo.HasOsculating() && !aOscSurf.IsNull()) + { + ComputeDerivatives(aMaxOrder, + 2, + theU, + theV, + theBasisSurf, + 1, + 1, + aOscInfo.AlongU, + aOscInfo.AlongV, + aOscSurf, + aDerNUV, + aDerSurf); + } + else + { + Handle(Geom_BSplineSurface) aDummy; + ComputeDerivatives(aMaxOrder, + 2, + theU, + theV, + theBasisSurf, + 1, + 1, + false, + false, + aDummy, + aDerNUV, + aDerSurf); + } + CSLib::Normal(aMaxOrder, + aDerNUV, + THE_D1_MAGNITUDE_TOL, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + aNStatus, + aNormal, + aOrderU, + aOrderV); + } + } + + if (aNStatus == CSLib_Defined) + { + // Compute offset point + theValue.SetXYZ(theValue.XYZ() + theOffset * aSign * aNormal.XYZ()); + // Compute D1 using CSLib + theD1U = CSLib::DNNormal(1, 0, aDerNUV, aOrderU, aOrderV); + theD1V = CSLib::DNNormal(0, 1, aDerNUV, aOrderU, aOrderV); + + theD1U.Multiply(theOffset * aSign); + theD1U.Add(aDerSurf(1, 0)); + theD1V.Multiply(theOffset * aSign); + theD1V.Add(aDerSurf(0, 1)); + return true; + } + + // Try shifting point towards center - returns false when center is overpassed + if (!ShiftPoint(aUStart, + aVStart, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + isUPer, + isVPer, + theD1U, + theD1V)) + { + return false; + } + } +} + +//! Template function for D2 evaluation with retry mechanism for singular points. +//! +//! @tparam BasisSurfType type of basis surface (must have D3 method) +//! @tparam OscSurfQueryType type providing osculating surface query +//! @param[in] theU U parameter +//! @param[in] theV V parameter +//! @param[in] theBasisSurf basis surface adaptor +//! @param[in] theOffset offset distance +//! @param[in] theOscQuery osculating surface query object (may be null) +//! @param[out] theValue computed offset point +//! @param[out] theD1U computed D1U derivative +//! @param[out] theD1V computed D1V derivative +//! @param[out] theD2U computed D2U derivative +//! @param[out] theD2V computed D2V derivative +//! @param[out] theD2UV computed D2UV derivative +//! @return true if calculation succeeded, false if failed at singular point +template +bool EvaluateD2(double theU, + double theV, + const BasisSurfType& theBasisSurf, + double theOffset, + const OscSurfQueryType& theOscQuery, + gp_Pnt& theValue, + gp_Vec& theD1U, + gp_Vec& theD1V, + gp_Vec& theD2U, + gp_Vec& theD2V, + gp_Vec& theD2UV) +{ + const double aUStart = theU; + const double aVStart = theV; + double aUMin, aUMax, aVMin, aVMax; + theBasisSurf->Bounds(aUMin, aUMax, aVMin, aVMax); + const bool isUPer = theBasisSurf->IsUPeriodic(); + const bool isVPer = theBasisSurf->IsVPeriodic(); + + for (;;) + { + gp_Vec aD3U, aD3V, aD3UUV, aD3UVV; + theBasisSurf->D3(theU, + theV, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV, + aD3U, + aD3V, + aD3UUV, + aD3UVV); + + if (IsInfiniteCoord(theD1U) || IsInfiniteCoord(theD1V)) + { + return false; + } + + // Check if singular using CSLib::Normal + gp_Dir aNormal; + CSLib_NormalStatus aNStatus; + CSLib::Normal(theD1U, theD1V, THE_D1_MAGNITUDE_TOL, aNStatus, aNormal); + + // MaxOrder = 0 for non-singular, 3 for singular + const int aMaxOrder = (aNStatus == CSLib_Defined) ? 0 : 3; + + // Get osculating surface info + OsculatingInfo aOscInfo; + Handle(Geom_BSplineSurface) aOscSurf; + if ((aNStatus != CSLib_Defined) && theOscQuery) + { + aOscInfo.AlongU = theOscQuery->UOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + aOscInfo.AlongV = theOscQuery->VOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + } + + // Setup derivative arrays + const int aDerNUVSize = (aMaxOrder + 3) * (aMaxOrder + 3); + const int aDerSurfSize = (aMaxOrder + 4) * (aMaxOrder + 4); + NCollection_LocalArray aDerNUVBuffer(aDerNUVSize); + NCollection_LocalArray aDerSurfBuffer(aDerSurfSize); + NCollection_Array2 aDerNUV(aDerNUVBuffer[0], 0, aMaxOrder + 2, 0, aMaxOrder + 2); + NCollection_Array2 aDerSurf(aDerSurfBuffer[0], 0, aMaxOrder + 3, 0, aMaxOrder + 3); + + aDerSurf.SetValue(1, 0, theD1U); + aDerSurf.SetValue(0, 1, theD1V); + aDerSurf.SetValue(1, 1, theD2UV); + aDerSurf.SetValue(2, 0, theD2U); + aDerSurf.SetValue(0, 2, theD2V); + aDerSurf.SetValue(3, 0, aD3U); + aDerSurf.SetValue(2, 1, aD3UUV); + aDerSurf.SetValue(1, 2, aD3UVV); + aDerSurf.SetValue(0, 3, aD3V); + + // Use ComputeDerivatives to populate DerNUV + if (aOscInfo.HasOsculating() && !aOscSurf.IsNull()) + { + ComputeDerivatives(aMaxOrder, + 3, + theU, + theV, + theBasisSurf, + 2, + 2, + aOscInfo.AlongU, + aOscInfo.AlongV, + aOscSurf, + aDerNUV, + aDerSurf); + } + else + { + Handle(Geom_BSplineSurface) aDummy; + ComputeDerivatives(aMaxOrder, + 3, + theU, + theV, + theBasisSurf, + 2, + 2, + false, + false, + aDummy, + aDerNUV, + aDerSurf); + } + + // Compute normal using CSLib (second call with MaxOrder) + int aOrderU, aOrderV; + CSLib::Normal(aMaxOrder, + aDerNUV, + THE_D1_MAGNITUDE_TOL, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + aNStatus, + aNormal, + aOrderU, + aOrderV); + + if (aNStatus == CSLib_Defined) + { + const double aSign = theOffset * aOscInfo.Sign(); + + // Compute offset point + theValue.SetXYZ(theValue.XYZ() + aSign * aNormal.XYZ()); + + // Compute D1 using CSLib::DNNormal + theD1U = + aDerSurf(1, 0).Added(CSLib::DNNormal(1, 0, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD1V = + aDerSurf(0, 1).Added(CSLib::DNNormal(0, 1, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + + // For D2, re-fetch from basis surface + theD2U = theBasisSurf->DN(theU, theV, 2, 0); + theD2V = theBasisSurf->DN(theU, theV, 0, 2); + theD2UV = theBasisSurf->DN(theU, theV, 1, 1); + + // Add offset corrections using CSLib::DNNormal + theD2U.Add(CSLib::DNNormal(2, 0, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD2V.Add(CSLib::DNNormal(0, 2, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD2UV.Add(CSLib::DNNormal(1, 1, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + return true; + } + + // Try shifting point towards center - returns false when center is overpassed + if (!ShiftPoint(aUStart, + aVStart, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + isUPer, + isVPer, + theD1U, + theD1V)) + { + return false; + } + } +} + +//! Template function for D3 evaluation with retry mechanism for singular points. +//! +//! @tparam BasisSurfType type of basis surface (must have D3 method) +//! @tparam OscSurfQueryType type providing osculating surface query +//! @param[in] theU U parameter +//! @param[in] theV V parameter +//! @param[in] theBasisSurf basis surface adaptor +//! @param[in] theOffset offset distance +//! @param[in] theOscQuery osculating surface query object (may be null) +//! @param[out] theValue computed offset point +//! @param[out] theD1U computed D1U derivative +//! @param[out] theD1V computed D1V derivative +//! @param[out] theD2U computed D2U derivative +//! @param[out] theD2V computed D2V derivative +//! @param[out] theD2UV computed D2UV derivative +//! @param[out] theD3U computed D3U derivative +//! @param[out] theD3V computed D3V derivative +//! @param[out] theD3UUV computed D3UUV derivative +//! @param[out] theD3UVV computed D3UVV derivative +//! @return true if calculation succeeded, false if failed at singular point +template +bool EvaluateD3(double theU, + double theV, + const BasisSurfType& theBasisSurf, + double theOffset, + const OscSurfQueryType& theOscQuery, + gp_Pnt& theValue, + gp_Vec& theD1U, + gp_Vec& theD1V, + gp_Vec& theD2U, + gp_Vec& theD2V, + gp_Vec& theD2UV, + gp_Vec& theD3U, + gp_Vec& theD3V, + gp_Vec& theD3UUV, + gp_Vec& theD3UVV) +{ + const double aUStart = theU; + const double aVStart = theV; + double aUMin, aUMax, aVMin, aVMax; + theBasisSurf->Bounds(aUMin, aUMax, aVMin, aVMax); + const bool isUPer = theBasisSurf->IsUPeriodic(); + const bool isVPer = theBasisSurf->IsVPeriodic(); + + for (;;) + { + theBasisSurf->D3(theU, + theV, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV, + theD3U, + theD3V, + theD3UUV, + theD3UVV); + + if (IsInfiniteCoord(theD1U) || IsInfiniteCoord(theD1V)) + { + return false; + } + + // Check if singular using CSLib::Normal + gp_Dir aNormal; + CSLib_NormalStatus aNStatus; + CSLib::Normal(theD1U, theD1V, THE_D1_MAGNITUDE_TOL, aNStatus, aNormal); + + // MaxOrder = 0 for non-singular, 3 for singular + const int aMaxOrder = (aNStatus == CSLib_Defined) ? 0 : 3; + + // Get osculating surface info + OsculatingInfo aOscInfo; + Handle(Geom_BSplineSurface) aOscSurf; + if ((aNStatus != CSLib_Defined) && theOscQuery) + { + aOscInfo.AlongU = theOscQuery->UOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + aOscInfo.AlongV = theOscQuery->VOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + } + + // Setup derivative arrays + // For D3: DerNUV needs (aMaxOrder + 4), DerSurf needs (aMaxOrder + 5) + const int aDerNUVSize = (aMaxOrder + 4) * (aMaxOrder + 4); + const int aDerSurfSize = (aMaxOrder + 5) * (aMaxOrder + 5); + NCollection_LocalArray aDerNUVBuffer(aDerNUVSize); + NCollection_LocalArray aDerSurfBuffer(aDerSurfSize); + NCollection_Array2 aDerNUV(aDerNUVBuffer[0], 0, aMaxOrder + 3, 0, aMaxOrder + 3); + NCollection_Array2 aDerSurf(aDerSurfBuffer[0], 0, aMaxOrder + 4, 0, aMaxOrder + 4); + + aDerSurf.SetValue(1, 0, theD1U); + aDerSurf.SetValue(0, 1, theD1V); + aDerSurf.SetValue(1, 1, theD2UV); + aDerSurf.SetValue(2, 0, theD2U); + aDerSurf.SetValue(0, 2, theD2V); + aDerSurf.SetValue(3, 0, theD3U); + aDerSurf.SetValue(2, 1, theD3UUV); + aDerSurf.SetValue(1, 2, theD3UVV); + aDerSurf.SetValue(0, 3, theD3V); + + // Use ComputeDerivatives to populate DerNUV + if (aOscInfo.HasOsculating() && !aOscSurf.IsNull()) + { + ComputeDerivatives(aMaxOrder, + 3, + theU, + theV, + theBasisSurf, + 3, + 3, + aOscInfo.AlongU, + aOscInfo.AlongV, + aOscSurf, + aDerNUV, + aDerSurf); + } + else + { + Handle(Geom_BSplineSurface) aDummy; + ComputeDerivatives(aMaxOrder, + 3, + theU, + theV, + theBasisSurf, + 3, + 3, + false, + false, + aDummy, + aDerNUV, + aDerSurf); + } + + // Compute normal using CSLib (second call with MaxOrder) + int aOrderU, aOrderV; + CSLib::Normal(aMaxOrder, + aDerNUV, + THE_D1_MAGNITUDE_TOL, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + aNStatus, + aNormal, + aOrderU, + aOrderV); + + if (aNStatus == CSLib_Defined) + { + const double aSign = theOffset * aOscInfo.Sign(); + + // Compute offset point + theValue.SetXYZ(theValue.XYZ() + aSign * aNormal.XYZ()); + + // Compute D1 using CSLib::DNNormal + theD1U = + aDerSurf(1, 0).Added(CSLib::DNNormal(1, 0, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD1V = + aDerSurf(0, 1).Added(CSLib::DNNormal(0, 1, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + + // For D2 and D3, re-fetch from basis surface + theD2U = theBasisSurf->DN(theU, theV, 2, 0); + theD2V = theBasisSurf->DN(theU, theV, 0, 2); + theD2UV = theBasisSurf->DN(theU, theV, 1, 1); + theD3U = theBasisSurf->DN(theU, theV, 3, 0); + theD3V = theBasisSurf->DN(theU, theV, 0, 3); + theD3UUV = theBasisSurf->DN(theU, theV, 2, 1); + theD3UVV = theBasisSurf->DN(theU, theV, 1, 2); + + // Add offset corrections using CSLib::DNNormal + theD2U.Add(CSLib::DNNormal(2, 0, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD2V.Add(CSLib::DNNormal(0, 2, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD2UV.Add(CSLib::DNNormal(1, 1, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD3U.Add(CSLib::DNNormal(3, 0, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD3V.Add(CSLib::DNNormal(0, 3, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD3UUV.Add(CSLib::DNNormal(2, 1, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + theD3UVV.Add(CSLib::DNNormal(1, 2, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + return true; + } + + // Try shifting point towards center - returns false when center is overpassed + if (!ShiftPoint(aUStart, + aVStart, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + isUPer, + isVPer, + theD1U, + theD1V)) + { + return false; + } + } +} + +//! Template function for DN evaluation with retry mechanism for singular points. +//! Computes arbitrary order derivative of offset surface. +//! +//! @tparam BasisSurfType type of basis surface (must have D1 and DN methods) +//! @tparam OscSurfQueryType type providing osculating surface query +//! @param[in] theU U parameter +//! @param[in] theV V parameter +//! @param[in] theNu derivative order in U +//! @param[in] theNv derivative order in V +//! @param[in] theBasisSurf basis surface adaptor +//! @param[in] theOffset offset distance +//! @param[in] theOscQuery osculating surface query object (may be null) +//! @param[out] theResult computed derivative vector +//! @return true if calculation succeeded, false if failed at singular point +template +bool EvaluateDN(double theU, + double theV, + int theNu, + int theNv, + const BasisSurfType& theBasisSurf, + double theOffset, + const OscSurfQueryType& theOscQuery, + gp_Vec& theResult) +{ + const double aUStart = theU; + const double aVStart = theV; + double aUMin, aUMax, aVMin, aVMax; + theBasisSurf->Bounds(aUMin, aUMax, aVMin, aVMax); + const bool isUPer = theBasisSurf->IsUPeriodic(); + const bool isVPer = theBasisSurf->IsVPeriodic(); + + for (;;) + { + gp_Pnt aP; + gp_Vec aD1U, aD1V; + theBasisSurf->D1(theU, theV, aP, aD1U, aD1V); + + if (IsInfiniteCoord(aD1U) || IsInfiniteCoord(aD1V)) + { + return false; + } + + // Check if singular to determine MaxOrder + gp_Dir aNormal; + CSLib_NormalStatus aNStatus; + CSLib::Normal(aD1U, aD1V, THE_D1_MAGNITUDE_TOL, aNStatus, aNormal); + + const int aMaxOrder = (aNStatus == CSLib_Defined) ? 0 : 3; + int aOrderU, aOrderV; + + // Stack buffers: max size with aMaxOrder=3, theNu=theNv=3 is 7x7=49 and 8x8=64 + const int aDerNUVSize = (aMaxOrder + theNu + 1) * (aMaxOrder + theNv + 1); + const int aDerSurfSize = (aMaxOrder + theNu + 2) * (aMaxOrder + theNv + 2); + NCollection_LocalArray aDerNUVBuffer(aDerNUVSize); + NCollection_LocalArray aDerSurfBuffer(aDerSurfSize); + NCollection_Array2 aDerNUV(aDerNUVBuffer[0], + 0, + aMaxOrder + theNu, + 0, + aMaxOrder + theNv); + NCollection_Array2 aDerSurf(aDerSurfBuffer[0], + 0, + aMaxOrder + theNu + 1, + 0, + aMaxOrder + theNv + 1); + + aDerSurf.SetValue(1, 0, aD1U); + aDerSurf.SetValue(0, 1, aD1V); + + // Check osculating surface only in singular case + OsculatingInfo aOscInfo; + Handle(Geom_BSplineSurface) aOscSurf; + if ((aNStatus != CSLib_Defined) && theOscQuery) + { + aOscInfo.AlongU = theOscQuery->UOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + aOscInfo.AlongV = theOscQuery->VOscSurf(theU, theV, aOscInfo.IsOpposite, aOscSurf); + } + + // Use ComputeDerivatives + if (aOscInfo.HasOsculating() && !aOscSurf.IsNull()) + { + ComputeDerivatives(aMaxOrder, + 1, + theU, + theV, + theBasisSurf, + theNu, + theNv, + aOscInfo.AlongU, + aOscInfo.AlongV, + aOscSurf, + aDerNUV, + aDerSurf); + } + else + { + Handle(Geom_BSplineSurface) aDummy; + ComputeDerivatives(aMaxOrder, + 1, + theU, + theV, + theBasisSurf, + theNu, + theNv, + false, + false, + aDummy, + aDerNUV, + aDerSurf); + } + + // Compute normal with CSLib + CSLib::Normal(aMaxOrder, + aDerNUV, + THE_D1_MAGNITUDE_TOL, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + aNStatus, + aNormal, + aOrderU, + aOrderV); + + if (aNStatus == CSLib_Defined) + { + const double aSign = theOffset * aOscInfo.Sign(); + + // Compute DN result: basis DN + offset * DNNormal + theResult = theBasisSurf->DN(theU, theV, theNu, theNv); + theResult.Add(CSLib::DNNormal(theNu, theNv, aDerNUV, aOrderU, aOrderV).Multiplied(aSign)); + return true; + } + + // Try shifting point towards center - returns false when center is overpassed + if (!ShiftPoint(aUStart, + aVStart, + theU, + theV, + aUMin, + aUMax, + aVMin, + aVMax, + isUPer, + isVPer, + aD1U, + aD1V)) + { + return false; + } + } +} + +} // namespace Geom_OffsetSurfaceUtils + +#endif // _Geom_OffsetSurfaceUtils_HeaderFile diff --git a/src/ModelingData/TKG3d/Geom/Geom_RevolutionUtils.pxx b/src/ModelingData/TKG3d/Geom/Geom_RevolutionUtils.pxx new file mode 100644 index 0000000000..2c2fa7cc32 --- /dev/null +++ b/src/ModelingData/TKG3d/Geom/Geom_RevolutionUtils.pxx @@ -0,0 +1,272 @@ +// Copyright (c) 2025 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 _Geom_RevolutionUtils_HeaderFile +#define _Geom_RevolutionUtils_HeaderFile + +#include +#include +#include +#include +#include + +//! @file Geom_RevolutionUtils.pxx +//! @brief Shared utility functions for surface of revolution evaluation. +//! +//! This file provides template functions for evaluating points and derivatives +//! on surfaces of revolution. The functions are templated to work with both +//! Geom_Curve (for Geom_SurfaceOfRevolution) and Adaptor3d_Curve +//! (for GeomAdaptor_SurfaceOfRevolution). +//! +//! Revolution surface: P(U,V) = Rotation(Axis, U) * BasisCurve(V) +//! where U is the rotation angle and V is the parameter along the basis curve. + +namespace Geom_RevolutionUtils +{ + +//! Evaluates point on surface of revolution. +//! @tparam CurveType Type supporting D0(param, point) method +//! @param theU Rotation angle parameter +//! @param theV Parameter along the basis curve +//! @param theBasis Basis curve +//! @param theAxis Rotation axis +//! @param theP [out] Evaluated point +template +inline void D0(const double theU, + const double theV, + const CurveType& theBasis, + const gp_Ax1& theAxis, + gp_Pnt& theP) +{ + theBasis.D0(theV, theP); + + gp_Trsf aRotation; + aRotation.SetRotation(theAxis, theU); + theP.Transform(aRotation); +} + +//! Evaluates point and first derivatives on surface of revolution. +//! @tparam CurveType Type supporting D1(param, point, vec) method +//! @param theU Rotation angle parameter +//! @param theV Parameter along the basis curve +//! @param theBasis Basis curve +//! @param theAxis Rotation axis +//! @param theP [out] Evaluated point +//! @param theD1U [out] First derivative with respect to U (rotation) +//! @param theD1V [out] First derivative with respect to V (along curve) +template +inline void D1(const double theU, + const double theV, + const CurveType& theBasis, + const gp_Ax1& theAxis, + gp_Pnt& theP, + gp_Vec& theD1U, + gp_Vec& theD1V) +{ + theBasis.D1(theV, theP, theD1V); + + // Vector from center of rotation to the point on rotated curve + gp_XYZ aCQ = theP.XYZ() - theAxis.Location().XYZ(); + theD1U = gp_Vec(theAxis.Direction().XYZ().Crossed(aCQ)); + // If the point is placed on the axis of revolution then derivatives on U are undefined. + // Manually set them to zero. + if (theD1U.SquareMagnitude() < Precision::SquareConfusion()) + { + theD1U.SetCoord(0.0, 0.0, 0.0); + } + + gp_Trsf aRotation; + aRotation.SetRotation(theAxis, theU); + theP.Transform(aRotation); + theD1U.Transform(aRotation); + theD1V.Transform(aRotation); +} + +//! Evaluates point, first and second derivatives on surface of revolution. +//! @tparam CurveType Type supporting D2(param, point, vec, vec) method +//! @param theU Rotation angle parameter +//! @param theV Parameter along the basis curve +//! @param theBasis Basis curve +//! @param theAxis Rotation axis +//! @param theP [out] Evaluated point +//! @param theD1U [out] First derivative with respect to U +//! @param theD1V [out] First derivative with respect to V +//! @param theD2U [out] Second derivative with respect to U +//! @param theD2V [out] Second derivative with respect to V +//! @param theD2UV [out] Mixed second derivative +template +inline void D2(const double theU, + const double theV, + const CurveType& theBasis, + const gp_Ax1& theAxis, + gp_Pnt& theP, + gp_Vec& theD1U, + gp_Vec& theD1V, + gp_Vec& theD2U, + gp_Vec& theD2V, + gp_Vec& theD2UV) +{ + theBasis.D2(theV, theP, theD1V, theD2V); + + // Vector from center of rotation to the point on rotated curve + gp_XYZ aCQ = theP.XYZ() - theAxis.Location().XYZ(); + const gp_XYZ& aDir = theAxis.Direction().XYZ(); + theD1U = gp_Vec(aDir.Crossed(aCQ)); + // If the point is placed on the axis of revolution then derivatives on U are undefined. + // Manually set them to zero. + if (theD1U.SquareMagnitude() < Precision::SquareConfusion()) + { + theD1U.SetCoord(0.0, 0.0, 0.0); + } + theD2U = gp_Vec(aDir.Dot(aCQ) * aDir - aCQ); + theD2UV = gp_Vec(aDir.Crossed(theD1V.XYZ())); + + gp_Trsf aRotation; + aRotation.SetRotation(theAxis, theU); + theP.Transform(aRotation); + theD1U.Transform(aRotation); + theD1V.Transform(aRotation); + theD2U.Transform(aRotation); + theD2V.Transform(aRotation); + theD2UV.Transform(aRotation); +} + +//! Evaluates point and all derivatives up to third order on surface of revolution. +//! @tparam CurveType Type supporting D3(param, point, vec, vec, vec) method +//! @param theU Rotation angle parameter +//! @param theV Parameter along the basis curve +//! @param theBasis Basis curve +//! @param theAxis Rotation axis +//! @param theP [out] Evaluated point +//! @param theD1U [out] First derivative with respect to U +//! @param theD1V [out] First derivative with respect to V +//! @param theD2U [out] Second derivative with respect to U +//! @param theD2V [out] Second derivative with respect to V +//! @param theD2UV [out] Mixed second derivative +//! @param theD3U [out] Third derivative with respect to U +//! @param theD3V [out] Third derivative with respect to V +//! @param theD3UUV [out] Mixed third derivative (UUV) +//! @param theD3UVV [out] Mixed third derivative (UVV) +template +inline void D3(const double theU, + const double theV, + const CurveType& theBasis, + const gp_Ax1& theAxis, + gp_Pnt& theP, + gp_Vec& theD1U, + gp_Vec& theD1V, + gp_Vec& theD2U, + gp_Vec& theD2V, + gp_Vec& theD2UV, + gp_Vec& theD3U, + gp_Vec& theD3V, + gp_Vec& theD3UUV, + gp_Vec& theD3UVV) +{ + theBasis.D3(theV, theP, theD1V, theD2V, theD3V); + + // Vector from center of rotation to the point on rotated curve + gp_XYZ aCQ = theP.XYZ() - theAxis.Location().XYZ(); + const gp_XYZ& aDir = theAxis.Direction().XYZ(); + theD1U = gp_Vec(aDir.Crossed(aCQ)); + // If the point is placed on the axis of revolution then derivatives on U are undefined. + // Manually set them to zero. + if (theD1U.SquareMagnitude() < Precision::SquareConfusion()) + { + theD1U.SetCoord(0.0, 0.0, 0.0); + } + theD2U = gp_Vec(aDir.Dot(aCQ) * aDir - aCQ); + theD2UV = gp_Vec(aDir.Crossed(theD1V.XYZ())); + theD3U = -theD1U; + theD3UUV = gp_Vec(aDir.Dot(theD1V.XYZ()) * aDir - theD1V.XYZ()); + theD3UVV = gp_Vec(aDir.Crossed(theD2V.XYZ())); + + gp_Trsf aRotation; + aRotation.SetRotation(theAxis, theU); + theP.Transform(aRotation); + theD1U.Transform(aRotation); + theD1V.Transform(aRotation); + theD2U.Transform(aRotation); + theD2V.Transform(aRotation); + theD2UV.Transform(aRotation); + theD3U.Transform(aRotation); + theD3V.Transform(aRotation); + theD3UUV.Transform(aRotation); + theD3UVV.Transform(aRotation); +} + +//! Evaluates N-th derivative on surface of revolution. +//! @tparam CurveType Type supporting D0, DN methods +//! @param theU Rotation angle parameter +//! @param theV Parameter along the basis curve +//! @param theBasis Basis curve +//! @param theAxis Rotation axis +//! @param theDerU Derivative order with respect to U +//! @param theDerV Derivative order with respect to V +//! @return The derivative vector +template +inline gp_Vec DN(const double theU, + const double theV, + const CurveType& theBasis, + const gp_Ax1& theAxis, + const int theDerU, + const int theDerV) +{ + gp_Trsf aRotation; + aRotation.SetRotation(theAxis, theU); + + gp_Pnt aP; + gp_Vec aDV; + gp_Vec aResult; + if (theDerU == 0) + { + aResult = theBasis.DN(theV, theDerV); + } + else + { + if (theDerV == 0) + { + theBasis.D0(theV, aP); + aDV = gp_Vec(aP.XYZ() - theAxis.Location().XYZ()); + } + else + { + aDV = theBasis.DN(theV, theDerV); + } + + const gp_XYZ& aDir = theAxis.Direction().XYZ(); + if (theDerU % 4 == 1) + { + aResult = gp_Vec(aDir.Crossed(aDV.XYZ())); + } + else if (theDerU % 4 == 2) + { + aResult = gp_Vec(aDir.Dot(aDV.XYZ()) * aDir - aDV.XYZ()); + } + else if (theDerU % 4 == 3) + { + aResult = gp_Vec(aDir.Crossed(aDV.XYZ())) * (-1.0); + } + else + { + aResult = gp_Vec(aDV.XYZ() - aDir.Dot(aDV.XYZ()) * aDir); + } + } + + aResult.Transform(aRotation); + return aResult; +} + +} // namespace Geom_RevolutionUtils + +#endif // _Geom_RevolutionUtils_HeaderFile diff --git a/src/ModelingData/TKG3d/GeomAdaptor/GeomAdaptor_Surface.hxx b/src/ModelingData/TKG3d/GeomAdaptor/GeomAdaptor_Surface.hxx index 849f4e6760..b92ef116dd 100644 --- a/src/ModelingData/TKG3d/GeomAdaptor/GeomAdaptor_Surface.hxx +++ b/src/ModelingData/TKG3d/GeomAdaptor/GeomAdaptor_Surface.hxx @@ -105,6 +105,22 @@ public: load(theSurf, theUFirst, theULast, theVFirst, theVLast, theTolU, theTolV); } + //! Returns the parametric bounds of the surface. + //! @param[out] theU1 minimum U parameter + //! @param[out] theU2 maximum U parameter + //! @param[out] theV1 minimum V parameter + //! @param[out] theV2 maximum V parameter + void Bounds(Standard_Real& theU1, + Standard_Real& theU2, + Standard_Real& theV1, + Standard_Real& theV2) const + { + theU1 = FirstUParameter(); + theU2 = LastUParameter(); + theV1 = FirstVParameter(); + theV2 = LastVParameter(); + } + const Handle(Geom_Surface)& Surface() const { return mySurface; } virtual Standard_Real FirstUParameter() const Standard_OVERRIDE { return myUFirst; } diff --git a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.cxx b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.cxx index 2e04c94af2..10bdaffd00 100644 --- a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.cxx +++ b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.cxx @@ -15,6 +15,7 @@ #include #include +#include #include IMPLEMENT_STANDARD_RTTIEXT(GeomEvaluator_OffsetCurve, GeomEvaluator_Curve) @@ -39,78 +40,137 @@ GeomEvaluator_OffsetCurve::GeomEvaluator_OffsetCurve(const Handle(GeomAdaptor_Cu { } +//================================================================================================== + void GeomEvaluator_OffsetCurve::D0(const Standard_Real theU, gp_Pnt& theValue) const { - gp_Vec aD1; - BaseD1(theU, theValue, aD1); - CalculateD0(theValue, aD1); + bool isOK = false; + if (!myBaseAdaptor.IsNull()) + isOK = Geom_OffsetCurveUtils::EvaluateD0(theU, myBaseAdaptor, myOffsetDir, myOffset, theValue); + else + isOK = Geom_OffsetCurveUtils::EvaluateD0(theU, myBaseCurve, myOffsetDir, myOffset, theValue); + + if (!isOK) + { + throw Standard_NullValue("GeomEvaluator_OffsetCurve: Undefined normal vector " + "because tangent vector has zero-magnitude!"); + } } +//================================================================================================== + void GeomEvaluator_OffsetCurve::D1(const Standard_Real theU, gp_Pnt& theValue, gp_Vec& theD1) const { - gp_Vec aD2; - BaseD2(theU, theValue, theD1, aD2); - CalculateD1(theValue, theD1, aD2); + bool isOK = false; + if (!myBaseAdaptor.IsNull()) + isOK = Geom_OffsetCurveUtils::EvaluateD1(theU, + myBaseAdaptor, + myOffsetDir, + myOffset, + theValue, + theD1); + else + isOK = + Geom_OffsetCurveUtils::EvaluateD1(theU, myBaseCurve, myOffsetDir, myOffset, theValue, theD1); + + if (!isOK) + { + throw Standard_NullValue("GeomEvaluator_OffsetCurve: Null derivative"); + } } +//================================================================================================== + void GeomEvaluator_OffsetCurve::D2(const Standard_Real theU, gp_Pnt& theValue, gp_Vec& theD1, gp_Vec& theD2) const { - gp_Vec aD3; - BaseD3(theU, theValue, theD1, theD2, aD3); + bool isOK = false; + if (!myBaseAdaptor.IsNull()) + isOK = Geom_OffsetCurveUtils::EvaluateD2(theU, + myBaseAdaptor, + myOffsetDir, + myOffset, + theValue, + theD1, + theD2); + else + isOK = Geom_OffsetCurveUtils::EvaluateD2(theU, + myBaseCurve, + myOffsetDir, + myOffset, + theValue, + theD1, + theD2); - Standard_Boolean isDirectionChange = Standard_False; - if (theD1.SquareMagnitude() <= gp::Resolution()) + if (!isOK) { - gp_Vec aDummyD4; - isDirectionChange = AdjustDerivative(3, theU, theD1, theD2, aD3, aDummyD4); + throw Standard_NullValue("GeomEvaluator_OffsetCurve: Null derivative"); } - - CalculateD2(theValue, theD1, theD2, aD3, isDirectionChange); } +//================================================================================================== + void GeomEvaluator_OffsetCurve::D3(const Standard_Real theU, gp_Pnt& theValue, gp_Vec& theD1, gp_Vec& theD2, gp_Vec& theD3) const { - gp_Vec aD4; - BaseD4(theU, theValue, theD1, theD2, theD3, aD4); + bool isOK = false; + if (!myBaseAdaptor.IsNull()) + isOK = Geom_OffsetCurveUtils::EvaluateD3(theU, + myBaseAdaptor, + myOffsetDir, + myOffset, + theValue, + theD1, + theD2, + theD3); + else + isOK = Geom_OffsetCurveUtils::EvaluateD3(theU, + myBaseCurve, + myOffsetDir, + myOffset, + theValue, + theD1, + theD2, + theD3); - Standard_Boolean isDirectionChange = Standard_False; - if (theD1.SquareMagnitude() <= gp::Resolution()) - isDirectionChange = AdjustDerivative(4, theU, theD1, theD2, theD3, aD4); - - CalculateD3(theValue, theD1, theD2, theD3, aD4, isDirectionChange); + if (!isOK) + { + throw Standard_NullValue("GeomEvaluator_OffsetCurve: Null derivative"); + } } +//================================================================================================== + gp_Vec GeomEvaluator_OffsetCurve::DN(const Standard_Real theU, const Standard_Integer theDeriv) const { Standard_RangeError_Raise_if(theDeriv < 1, "GeomEvaluator_OffsetCurve::DN(): theDeriv < 1"); - gp_Pnt aPnt; - gp_Vec aDummy, aDN; - switch (theDeriv) + gp_Vec aDN; + bool isOK = false; + + if (!myBaseAdaptor.IsNull()) + isOK = + Geom_OffsetCurveUtils::EvaluateDN(theU, myBaseAdaptor, myOffsetDir, myOffset, theDeriv, aDN); + else + isOK = + Geom_OffsetCurveUtils::EvaluateDN(theU, myBaseCurve, myOffsetDir, myOffset, theDeriv, aDN); + + if (!isOK) { - case 1: - D1(theU, aPnt, aDN); - break; - case 2: - D2(theU, aPnt, aDummy, aDN); - break; - case 3: - D3(theU, aPnt, aDummy, aDummy, aDN); - break; - default: - aDN = BaseDN(theU, theDeriv); + throw Standard_NullValue("GeomEvaluator_OffsetCurve: Null derivative"); } + return aDN; } +//================================================================================================== + Handle(GeomEvaluator_Curve) GeomEvaluator_OffsetCurve::ShallowCopy() const { Handle(GeomEvaluator_OffsetCurve) aCopy; @@ -127,340 +187,3 @@ Handle(GeomEvaluator_Curve) GeomEvaluator_OffsetCurve::ShallowCopy() const } return aCopy; } - -void GeomEvaluator_OffsetCurve::BaseD0(const Standard_Real theU, gp_Pnt& theValue) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D0(theU, theValue); - else - myBaseCurve->D0(theU, theValue); -} - -void GeomEvaluator_OffsetCurve::BaseD1(const Standard_Real theU, - gp_Pnt& theValue, - gp_Vec& theD1) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D1(theU, theValue, theD1); - else - myBaseCurve->D1(theU, theValue, theD1); -} - -void GeomEvaluator_OffsetCurve::BaseD2(const Standard_Real theU, - gp_Pnt& theValue, - gp_Vec& theD1, - gp_Vec& theD2) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D2(theU, theValue, theD1, theD2); - else - myBaseCurve->D2(theU, theValue, theD1, theD2); -} - -void GeomEvaluator_OffsetCurve::BaseD3(const Standard_Real theU, - gp_Pnt& theValue, - gp_Vec& theD1, - gp_Vec& theD2, - gp_Vec& theD3) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D3(theU, theValue, theD1, theD2, theD3); - else - myBaseCurve->D3(theU, theValue, theD1, theD2, theD3); -} - -void GeomEvaluator_OffsetCurve::BaseD4(const Standard_Real theU, - gp_Pnt& theValue, - gp_Vec& theD1, - gp_Vec& theD2, - gp_Vec& theD3, - gp_Vec& theD4) const -{ - if (!myBaseAdaptor.IsNull()) - { - myBaseAdaptor->D3(theU, theValue, theD1, theD2, theD3); - theD4 = myBaseAdaptor->DN(theU, 4); - } - else - { - myBaseCurve->D3(theU, theValue, theD1, theD2, theD3); - theD4 = myBaseCurve->DN(theU, 4); - } -} - -gp_Vec GeomEvaluator_OffsetCurve::BaseDN(const Standard_Real theU, - const Standard_Integer theDeriv) const -{ - if (!myBaseAdaptor.IsNull()) - return myBaseAdaptor->DN(theU, theDeriv); - return myBaseCurve->DN(theU, theDeriv); -} - -void GeomEvaluator_OffsetCurve::CalculateD0(gp_Pnt& theValue, const gp_Vec& theD1) const -{ - gp_XYZ Ndir = (theD1.XYZ()).Crossed(myOffsetDir.XYZ()); - Standard_Real R = Ndir.Modulus(); - if (R <= gp::Resolution()) - throw Standard_NullValue("GeomEvaluator_OffsetCurve: Undefined normal vector " - "because tangent vector has zero-magnitude!"); - - Ndir.Multiply(myOffset / R); - theValue.ChangeCoord().Add(Ndir); -} - -void GeomEvaluator_OffsetCurve::CalculateD1(gp_Pnt& theValue, - gp_Vec& theD1, - const gp_Vec& theD2) const -{ - // P(u) = p(u) + Offset * Ndir / R - // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction) - - // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) - - gp_XYZ Ndir = (theD1.XYZ()).Crossed(myOffsetDir.XYZ()); - gp_XYZ DNdir = (theD2.XYZ()).Crossed(myOffsetDir.XYZ()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = std::sqrt(R2); - Standard_Real R3 = R * R2; - Standard_Real Dr = Ndir.Dot(DNdir); - if (R3 <= gp::Resolution()) - { - if (R2 <= gp::Resolution()) - throw Standard_NullValue("GeomEvaluator_OffsetCurve: Null derivative"); - // We try another computation but the stability is not very good. - DNdir.Multiply(R); - DNdir.Subtract(Ndir.Multiplied(Dr / R)); - DNdir.Multiply(myOffset / R2); - } - else - { - // Same computation as IICURV in EUCLID-IS because the stability is better - DNdir.Multiply(myOffset / R); - DNdir.Subtract(Ndir.Multiplied(myOffset * Dr / R3)); - } - - Ndir.Multiply(myOffset / R); - // P(u) - theValue.ChangeCoord().Add(Ndir); - // P'(u) - theD1.Add(gp_Vec(DNdir)); -} - -void GeomEvaluator_OffsetCurve::CalculateD2(gp_Pnt& theValue, - gp_Vec& theD1, - gp_Vec& theD2, - const gp_Vec& theD3, - const Standard_Boolean theIsDirChange) const -{ - // P(u) = p(u) + Offset * Ndir / R - // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction) - - // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) - - // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) + - // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2))) - - gp_XYZ Ndir = (theD1.XYZ()).Crossed(myOffsetDir.XYZ()); - gp_XYZ DNdir = (theD2.XYZ()).Crossed(myOffsetDir.XYZ()); - gp_XYZ D2Ndir = (theD3.XYZ()).Crossed(myOffsetDir.XYZ()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = std::sqrt(R2); - Standard_Real R3 = R2 * R; - Standard_Real R4 = R2 * R2; - Standard_Real R5 = R3 * R2; - Standard_Real Dr = Ndir.Dot(DNdir); - Standard_Real D2r = Ndir.Dot(D2Ndir) + DNdir.Dot(DNdir); - - if (R5 <= gp::Resolution()) - { - if (R4 <= gp::Resolution()) - throw Standard_NullValue("GeomEvaluator_OffsetCurve: Null derivative"); - // We try another computation but the stability is not very good - // dixit ISG. - // V2 = P" (U) : - R4 = R2 * R2; - D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R2)); - D2Ndir.Add(Ndir.Multiplied(((3.0 * Dr * Dr) / R4) - (D2r / R2))); - D2Ndir.Multiply(myOffset / R); - - // V1 = P' (U) : - DNdir.Multiply(R); - DNdir.Subtract(Ndir.Multiplied(Dr / R)); - DNdir.Multiply(myOffset / R2); - } - else - { - // Same computation as IICURV in EUCLID-IS because the stability is better. - // V2 = P" (U) : - D2Ndir.Multiply(myOffset / R); - D2Ndir.Subtract(DNdir.Multiplied(2.0 * myOffset * Dr / R3)); - D2Ndir.Add(Ndir.Multiplied(myOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3)))); - - // V1 = P' (U) : - DNdir.Multiply(myOffset / R); - DNdir.Subtract(Ndir.Multiplied(myOffset * Dr / R3)); - } - - Ndir.Multiply(myOffset / R); - // P(u) - theValue.ChangeCoord().Add(Ndir); - // P'(u) : - theD1.Add(gp_Vec(DNdir)); - // P"(u) : - if (theIsDirChange) - theD2.Reverse(); - theD2.Add(gp_Vec(D2Ndir)); -} - -void GeomEvaluator_OffsetCurve::CalculateD3(gp_Pnt& theValue, - gp_Vec& theD1, - gp_Vec& theD2, - gp_Vec& theD3, - const gp_Vec& theD4, - const Standard_Boolean theIsDirChange) const -{ - // P(u) = p(u) + Offset * Ndir / R - // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction) - - // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) - - // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) + - // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2))) - - // P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir - - // (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir - - // (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir + - // (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir - - gp_XYZ Ndir = (theD1.XYZ()).Crossed(myOffsetDir.XYZ()); - gp_XYZ DNdir = (theD2.XYZ()).Crossed(myOffsetDir.XYZ()); - gp_XYZ D2Ndir = (theD3.XYZ()).Crossed(myOffsetDir.XYZ()); - gp_XYZ D3Ndir = (theD4.XYZ()).Crossed(myOffsetDir.XYZ()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = std::sqrt(R2); - Standard_Real R3 = R2 * R; - Standard_Real R4 = R2 * R2; - Standard_Real R5 = R3 * R2; - Standard_Real R6 = R3 * R3; - Standard_Real R7 = R5 * R2; - Standard_Real Dr = Ndir.Dot(DNdir); - Standard_Real D2r = Ndir.Dot(D2Ndir) + DNdir.Dot(DNdir); - Standard_Real D3r = Ndir.Dot(D3Ndir) + 3.0 * DNdir.Dot(D2Ndir); - if (R7 <= gp::Resolution()) - { - if (R6 <= gp::Resolution()) - throw Standard_NullValue("CSLib_Offset: Null derivative"); - // V3 = P"' (U) : - D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * Dr / R2)); - D3Ndir.Subtract(DNdir.Multiplied(3.0 * ((D2r / R2) + (Dr * Dr / R4)))); - D3Ndir.Add( - Ndir.Multiplied(6.0 * Dr * Dr / R4 + 6.0 * Dr * D2r / R4 - 15.0 * Dr * Dr * Dr / R6 - D3r)); - D3Ndir.Multiply(myOffset / R); - // V2 = P" (U) : - R4 = R2 * R2; - D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R2)); - D2Ndir.Subtract(Ndir.Multiplied((3.0 * Dr * Dr / R4) - (D2r / R2))); - D2Ndir.Multiply(myOffset / R); - // V1 = P' (U) : - DNdir.Multiply(R); - DNdir.Subtract(Ndir.Multiplied(Dr / R)); - DNdir.Multiply(myOffset / R2); - } - else - { - // V3 = P"' (U) : - D3Ndir.Divide(R); - D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * Dr / R3)); - D3Ndir.Subtract(DNdir.Multiplied((3.0 * ((D2r / R3) + (Dr * Dr) / R5)))); - D3Ndir.Add( - Ndir.Multiplied(6.0 * Dr * Dr / R5 + 6.0 * Dr * D2r / R5 - 15.0 * Dr * Dr * Dr / R7 - D3r)); - D3Ndir.Multiply(myOffset); - // V2 = P" (U) : - D2Ndir.Divide(R); - D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R3)); - D2Ndir.Subtract(Ndir.Multiplied((3.0 * Dr * Dr / R5) - (D2r / R3))); - D2Ndir.Multiply(myOffset); - // V1 = P' (U) : - DNdir.Multiply(myOffset / R); - DNdir.Subtract(Ndir.Multiplied(myOffset * Dr / R3)); - } - - Ndir.Multiply(myOffset / R); - // P(u) - theValue.ChangeCoord().Add(Ndir); - // P'(u) : - theD1.Add(gp_Vec(DNdir)); - // P"(u) - theD2.Add(gp_Vec(D2Ndir)); - // P"'(u) - if (theIsDirChange) - theD3.Reverse(); - theD3.Add(gp_Vec(D2Ndir)); -} - -Standard_Boolean GeomEvaluator_OffsetCurve::AdjustDerivative( - const Standard_Integer theMaxDerivative, - const Standard_Real theU, - gp_Vec& theD1, - gp_Vec& theD2, - gp_Vec& theD3, - gp_Vec& theD4) const -{ - static const Standard_Real aTol = gp::Resolution(); - static const Standard_Real aMinStep = 1e-7; - static const Standard_Integer aMaxDerivOrder = 3; - - Standard_Boolean isDirectionChange = Standard_False; - Standard_Real anUinfium; - Standard_Real anUsupremum; - if (!myBaseAdaptor.IsNull()) - { - anUinfium = myBaseAdaptor->FirstParameter(); - anUsupremum = myBaseAdaptor->LastParameter(); - } - else - { - anUinfium = myBaseCurve->FirstParameter(); - anUsupremum = myBaseCurve->LastParameter(); - } - - static const Standard_Real DivisionFactor = 1.e-3; - Standard_Real du; - if ((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) - du = 0.0; - else - du = anUsupremum - anUinfium; - - const Standard_Real aDelta = std::max(du * DivisionFactor, aMinStep); - - // Derivative is approximated by Taylor-series - Standard_Integer anIndex = 1; // Derivative order - gp_Vec V; - - do - { - V = BaseDN(theU, ++anIndex); - } while ((V.SquareMagnitude() <= aTol) && anIndex < aMaxDerivOrder); - - Standard_Real u; - - if (theU - anUinfium < aDelta) - u = theU + aDelta; - else - u = theU - aDelta; - - gp_Pnt P1, P2; - BaseD0(std::min(theU, u), P1); - BaseD0(std::max(theU, u), P2); - - gp_Vec V1(P1, P2); - isDirectionChange = V.Dot(V1) < 0.0; - Standard_Real aSign = isDirectionChange ? -1.0 : 1.0; - - theD1 = V * aSign; - gp_Vec* aDeriv[3] = {&theD2, &theD3, &theD4}; - for (Standard_Integer i = 1; i < theMaxDerivative; i++) - *(aDeriv[i - 1]) = BaseDN(theU, anIndex + i) * aSign; - - return isDirectionChange; -} diff --git a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.hxx b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.hxx index 22a9e13550..152767d3b7 100644 --- a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.hxx +++ b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.hxx @@ -62,62 +62,12 @@ public: DEFINE_STANDARD_RTTIEXT(GeomEvaluator_OffsetCurve, GeomEvaluator_Curve) -private: - //! Recalculate D1 values of base curve into D0 value of offset curve - void CalculateD0(gp_Pnt& theValue, const gp_Vec& theD1) const; - //! Recalculate D2 values of base curve into D1 values of offset curve - void CalculateD1(gp_Pnt& theValue, gp_Vec& theD1, const gp_Vec& theD2) const; - //! Recalculate D3 values of base curve into D2 values of offset curve - void CalculateD2(gp_Pnt& theValue, - gp_Vec& theD1, - gp_Vec& theD2, - const gp_Vec& theD3, - const Standard_Boolean theIsDirChange) const; - //! Recalculate D3 values of base curve into D3 values of offset curve - void CalculateD3(gp_Pnt& theValue, - gp_Vec& theD1, - gp_Vec& theD2, - gp_Vec& theD3, - const gp_Vec& theD4, - const Standard_Boolean theIsDirChange) const; - - //! Calculate value of base curve/adaptor - void BaseD0(const Standard_Real theU, gp_Pnt& theValue) const; - //! Calculate value and first derivatives of base curve/adaptor - void BaseD1(const Standard_Real theU, gp_Pnt& theValue, gp_Vec& theD1) const; - //! Calculate value, first and second derivatives of base curve/adaptor - void BaseD2(const Standard_Real theU, gp_Pnt& theValue, gp_Vec& theD1, gp_Vec& theD2) const; - //! Calculate value, first, second and third derivatives of base curve/adaptor - void BaseD3(const Standard_Real theU, - gp_Pnt& theValue, - gp_Vec& theD1, - gp_Vec& theD2, - gp_Vec& theD3) const; - //! Calculate value and derivatives till 4th of base curve/adaptor - void BaseD4(const Standard_Real theU, - gp_Pnt& theValue, - gp_Vec& theD1, - gp_Vec& theD2, - gp_Vec& theD3, - gp_Vec& theD4) const; - //! Calculate N-th derivative of base curve/adaptor - gp_Vec BaseDN(const Standard_Real theU, const Standard_Integer theDeriv) const; - - // Recalculate derivatives in the singular point - // Returns true if the direction of derivatives is changed - Standard_Boolean AdjustDerivative(const Standard_Integer theMaxDerivative, - const Standard_Real theU, - gp_Vec& theD1, - gp_Vec& theD2, - gp_Vec& theD3, - gp_Vec& theD4) const; - private: Handle(Geom_Curve) myBaseCurve; Handle(GeomAdaptor_Curve) myBaseAdaptor; - Standard_Real myOffset; ///< offset value - gp_Dir myOffsetDir; ///< offset direction + Standard_Real myOffset; //!< offset value + gp_Dir myOffsetDir; //!< offset direction }; DEFINE_STANDARD_HANDLE(GeomEvaluator_OffsetCurve, GeomEvaluator_Curve) diff --git a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx index 954891a84f..7517c94213 100644 --- a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx +++ b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx @@ -14,225 +14,15 @@ #include -#include -#include #include #include +#include #include #include -#include #include -#include IMPLEMENT_STANDARD_RTTIEXT(GeomEvaluator_OffsetSurface, GeomEvaluator_Surface) -namespace -{ - -// tolerance for considering derivative to be null -const Standard_Real the_D1MagTol = 1.e-9; - -// If calculation of normal fails, try shifting the point towards the center -// of the parametric space of the surface, in the hope that derivatives -// are better defined there. -// -// This shift is iterative, starting with Precision::PConfusion() -// and increasing by multiple of 2 on each step. -// -// NB: temporarily this is made as static function and not class method, -// hence code duplications -static Standard_Boolean shiftPoint(const Standard_Real theUStart, - const Standard_Real theVStart, - Standard_Real& theU, - Standard_Real& theV, - const Handle(Geom_Surface)& theSurf, - const Handle(GeomAdaptor_Surface)& theAdaptor, - const gp_Vec& theD1U, - const gp_Vec& theD1V) -{ - // Get parametric bounds and closure status - Standard_Real aUMin, aUMax, aVMin, aVMax; - Standard_Boolean isUPeriodic, isVPeriodic; - if (!theSurf.IsNull()) - { - theSurf->Bounds(aUMin, aUMax, aVMin, aVMax); - isUPeriodic = theSurf->IsUPeriodic(); - isVPeriodic = theSurf->IsVPeriodic(); - } - else - { - aUMin = theAdaptor->FirstUParameter(); - aUMax = theAdaptor->LastUParameter(); - aVMin = theAdaptor->FirstVParameter(); - aVMax = theAdaptor->LastVParameter(); - isUPeriodic = theAdaptor->IsUPeriodic(); - isVPeriodic = theAdaptor->IsVPeriodic(); - } - - // check if either U or V is singular (normally one of them is) - Standard_Boolean isUSingular = (theD1U.SquareMagnitude() < the_D1MagTol * the_D1MagTol); - Standard_Boolean isVSingular = (theD1V.SquareMagnitude() < the_D1MagTol * the_D1MagTol); - - // compute vector to shift from start point to center of the surface; - // if surface is periodic or singular in some direction, take shift in that direction zero - Standard_Real aDirU = - (isUPeriodic || (isUSingular && !isVSingular) ? 0. : 0.5 * (aUMin + aUMax) - theUStart); - Standard_Real aDirV = - (isVPeriodic || (isVSingular && !isUSingular) ? 0. : 0.5 * (aVMin + aVMax) - theVStart); - Standard_Real aDist = std::sqrt(aDirU * aDirU + aDirV * aDirV); - - // shift current point from its current position towards center, by value of twice - // current distance from it to start (but not less than Precision::PConfusion()); - // fail if center is overpassed. - Standard_Real aDU = theU - theUStart; - Standard_Real aDV = theV - theVStart; - Standard_Real aStep = std::max(2. * std::sqrt(aDU * aDU + aDV * aDV), Precision::PConfusion()); - if (aStep >= aDist) - { - return Standard_False; - } - - aStep /= aDist; - theU += aDirU * aStep; - theV += aDirV * aStep; - - // std::cout << "Normal calculation failed at (" << theUStart << ", " << theVStart << "), - // shifting to (" << theU << ", " << theV << ")" << std::endl; - - return Standard_True; -} - -// auxiliary function -template -static void derivatives(Standard_Integer theMaxOrder, - Standard_Integer theMinOrder, - const Standard_Real theU, - const Standard_Real theV, - const SurfOrAdapt& theBasisSurf, - const Standard_Integer theNU, - const Standard_Integer theNV, - const Standard_Boolean theAlongU, - const Standard_Boolean theAlongV, - const Handle(Geom_BSplineSurface)& theL, - TColgp_Array2OfVec& theDerNUV, - TColgp_Array2OfVec& theDerSurf) -{ - Standard_Integer i, j; - gp_Pnt P; - gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V; - - if (theAlongU || theAlongV) - { - theMaxOrder = 0; - TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1); - switch (theMinOrder) - { - case 1: - theL->D1(theU, theV, P, DL1U, DL1V); - DerSurfL.SetValue(1, 0, DL1U); - DerSurfL.SetValue(0, 1, DL1V); - break; - case 2: - theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV); - DerSurfL.SetValue(1, 0, DL1U); - DerSurfL.SetValue(0, 1, DL1V); - DerSurfL.SetValue(1, 1, DL2UV); - DerSurfL.SetValue(2, 0, DL2U); - DerSurfL.SetValue(0, 2, DL2V); - break; - case 3: - theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV); - DerSurfL.SetValue(1, 0, DL1U); - DerSurfL.SetValue(0, 1, DL1V); - DerSurfL.SetValue(1, 1, DL2UV); - DerSurfL.SetValue(2, 0, DL2U); - DerSurfL.SetValue(0, 2, DL2V); - DerSurfL.SetValue(3, 0, DL3U); - DerSurfL.SetValue(2, 1, DL3UUV); - DerSurfL.SetValue(1, 2, DL3UVV); - DerSurfL.SetValue(0, 3, DL3V); - break; - default: - break; - } - - if (theNU <= theNV) - { - for (i = 0; i <= theMaxOrder + 1 + theNU; i++) - for (j = i; j <= theMaxOrder + theNV + 1; j++) - if (i + j > theMinOrder) - { - DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j)); - theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j)); - if (i != j && j <= theNU + 1) - { - theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i)); - DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i)); - } - } - } - else - { - for (j = 0; j <= theMaxOrder + 1 + theNV; j++) - for (i = j; i <= theMaxOrder + theNU + 1; i++) - if (i + j > theMinOrder) - { - DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j)); - theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j)); - if (i != j && i <= theNV + 1) - { - theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i)); - DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i)); - } - } - } - for (i = 0; i <= theMaxOrder + theNU; i++) - for (j = 0; j <= theMaxOrder + theNV; j++) - { - if (theAlongU) - theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf)); - if (theAlongV) - theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL)); - } - } - else - { - for (i = 0; i <= theMaxOrder + theNU + 1; i++) - { - for (j = i; j <= theMaxOrder + theNV + 1; j++) - { - if (i + j > theMinOrder) - { - theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j)); - if (i != j && j <= theDerSurf.UpperRow() && i <= theDerSurf.UpperCol()) - { - theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i)); - } - } - } - } - for (i = 0; i <= theMaxOrder + theNU; i++) - for (j = 0; j <= theMaxOrder + theNV; j++) - theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf)); - } -} - -inline Standard_Boolean IsInfiniteCoord(const gp_Vec& theVec) -{ - return Precision::IsInfinite(theVec.X()) || Precision::IsInfinite(theVec.Y()) - || Precision::IsInfinite(theVec.Z()); -} - -inline void CheckInfinite(const gp_Vec& theVecU, const gp_Vec& theVecV) -{ - if (IsInfiniteCoord(theVecU) || IsInfiniteCoord(theVecV)) - { - throw Standard_NumericError("GeomEvaluator_OffsetSurface: Evaluation of infinite parameters"); - } -} - -} // end of anonymous namespace - GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface( const Handle(Geom_Surface)& theBase, const Standard_Real theOffset, @@ -262,64 +52,78 @@ GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface( { } +//================================================================================================== + void GeomEvaluator_OffsetSurface::D0(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue) const { - Standard_Real aU = theU, aV = theV; - for (;;) + bool isOK = false; + if (!myBaseAdaptor.IsNull()) { - gp_Vec aD1U, aD1V; - BaseD1(aU, aV, theValue, aD1U, aD1V); + isOK = Geom_OffsetSurfaceUtils::EvaluateD0(theU, + theV, + myBaseAdaptor, + myOffset, + myOscSurf.get(), + theValue); + } + else + { + isOK = Geom_OffsetSurfaceUtils::EvaluateD0(theU, + theV, + myBaseSurf, + myOffset, + myOscSurf.get(), + theValue); + } - CheckInfinite(aD1U, aD1V); - - try - { - CalculateD0(aU, aV, theValue, aD1U, aD1V); - break; - } - catch (Geom_UndefinedValue&) - { - // if failed at parametric boundary, try taking derivative at shifted point - if (!shiftPoint(theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, aD1U, aD1V)) - { - throw; - } - } + if (!isOK) + { + throw Geom_UndefinedValue("GeomEvaluator_OffsetSurface::D0(): Unable to calculate normal"); } } +//================================================================================================== + void GeomEvaluator_OffsetSurface::D1(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const { - Standard_Real aU = theU, aV = theV; - for (;;) + bool isOK = false; + if (!myBaseAdaptor.IsNull()) { - gp_Vec aD2U, aD2V, aD2UV; - BaseD2(aU, aV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV); + isOK = Geom_OffsetSurfaceUtils::EvaluateD1(theU, + theV, + myBaseAdaptor, + myOffset, + myOscSurf.get(), + theValue, + theD1U, + theD1V); + } + else + { + isOK = Geom_OffsetSurfaceUtils::EvaluateD1(theU, + theV, + myBaseSurf, + myOffset, + myOscSurf.get(), + theValue, + theD1U, + theD1V); + } - CheckInfinite(theD1U, theD1V); - - try - { - CalculateD1(aU, aV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV); - break; - } - catch (Geom_UndefinedValue&) - { - // if failed at parametric boundary, try taking derivative at shifted point - if (!shiftPoint(theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V)) - { - throw; - } - } + if (!isOK) + { + throw Geom_UndefinedValue("GeomEvaluator_OffsetSurface::D1(): Unable to calculate normal"); } } +//================================================================================================== + void GeomEvaluator_OffsetSurface::D2(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue, @@ -329,41 +133,44 @@ void GeomEvaluator_OffsetSurface::D2(const Standard_Real theU, gp_Vec& theD2V, gp_Vec& theD2UV) const { - Standard_Real aU = theU, aV = theV; - for (;;) + bool isOK = false; + if (!myBaseAdaptor.IsNull()) { - gp_Vec aD3U, aD3V, aD3UUV, aD3UVV; - BaseD3(aU, aV, theValue, theD1U, theD1V, theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV); + isOK = Geom_OffsetSurfaceUtils::EvaluateD2(theU, + theV, + myBaseAdaptor, + myOffset, + myOscSurf.get(), + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV); + } + else + { + isOK = Geom_OffsetSurfaceUtils::EvaluateD2(theU, + theV, + myBaseSurf, + myOffset, + myOscSurf.get(), + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV); + } - CheckInfinite(theD1U, theD1V); - - try - { - CalculateD2(aU, - aV, - theValue, - theD1U, - theD1V, - theD2U, - theD2V, - theD2UV, - aD3U, - aD3V, - aD3UUV, - aD3UVV); - break; - } - catch (Geom_UndefinedValue&) - { - // if failed at parametric boundary, try taking derivative at shifted point - if (!shiftPoint(theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V)) - { - throw; - } - } + if (!isOK) + { + throw Geom_UndefinedValue("GeomEvaluator_OffsetSurface::D2(): Unable to calculate normal"); } } +//================================================================================================== + void GeomEvaluator_OffsetSurface::D3(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue, @@ -377,51 +184,52 @@ void GeomEvaluator_OffsetSurface::D3(const Standard_Real theU, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const { - Standard_Real aU = theU, aV = theV; - for (;;) + bool isOK = false; + if (!myBaseAdaptor.IsNull()) { - BaseD3(aU, - aV, - theValue, - theD1U, - theD1V, - theD2U, - theD2V, - theD2UV, - theD3U, - theD3V, - theD3UUV, - theD3UVV); + isOK = Geom_OffsetSurfaceUtils::EvaluateD3(theU, + theV, + myBaseAdaptor, + myOffset, + myOscSurf.get(), + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV, + theD3U, + theD3V, + theD3UUV, + theD3UVV); + } + else + { + isOK = Geom_OffsetSurfaceUtils::EvaluateD3(theU, + theV, + myBaseSurf, + myOffset, + myOscSurf.get(), + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV, + theD3U, + theD3V, + theD3UUV, + theD3UVV); + } - CheckInfinite(theD1U, theD1V); - - try - { - CalculateD3(aU, - aV, - theValue, - theD1U, - theD1V, - theD2U, - theD2V, - theD2UV, - theD3U, - theD3V, - theD3UUV, - theD3UVV); - break; - } - catch (Geom_UndefinedValue&) - { - // if failed at parametric boundary, try taking derivative at shifted point - if (!shiftPoint(theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V)) - { - throw; - } - } + if (!isOK) + { + throw Geom_UndefinedValue("GeomEvaluator_OffsetSurface::D3(): Unable to calculate normal"); } } +//================================================================================================== + gp_Vec GeomEvaluator_OffsetSurface::DN(const Standard_Real theU, const Standard_Real theV, const Standard_Integer theDerU, @@ -432,30 +240,41 @@ gp_Vec GeomEvaluator_OffsetSurface::DN(const Standard_Real theU, Standard_RangeError_Raise_if(theDerU + theDerV < 1, "GeomEvaluator_OffsetSurface::DN(): theDerU + theDerV < 1"); - Standard_Real aU = theU, aV = theV; - for (;;) + gp_Vec aResult; + bool isOK = false; + if (!myBaseAdaptor.IsNull()) { - gp_Pnt aP; - gp_Vec aD1U, aD1V; - BaseD1(aU, aV, aP, aD1U, aD1V); - - CheckInfinite(aD1U, aD1V); - - try - { - return CalculateDN(aU, aV, theDerU, theDerV, aD1U, aD1V); - } - catch (Geom_UndefinedValue&) - { - // if failed at parametric boundary, try taking derivative at shifted point - if (!shiftPoint(theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, aD1U, aD1V)) - { - throw; - } - } + isOK = Geom_OffsetSurfaceUtils::EvaluateDN(theU, + theV, + theDerU, + theDerV, + myBaseAdaptor, + myOffset, + myOscSurf.get(), + aResult); } + else + { + isOK = Geom_OffsetSurfaceUtils::EvaluateDN(theU, + theV, + theDerU, + theDerV, + myBaseSurf, + myOffset, + myOscSurf.get(), + aResult); + } + + if (!isOK) + { + throw Geom_UndefinedValue("GeomEvaluator_OffsetSurface::DN(): Unable to calculate normal"); + } + + return aResult; } +//================================================================================================== + Handle(GeomEvaluator_Surface) GeomEvaluator_OffsetSurface::ShallowCopy() const { Handle(GeomEvaluator_OffsetSurface) aCopy; @@ -473,699 +292,3 @@ Handle(GeomEvaluator_Surface) GeomEvaluator_OffsetSurface::ShallowCopy() const return aCopy; } - -void GeomEvaluator_OffsetSurface::BaseD0(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D0(theU, theV, theValue); - else - myBaseSurf->D0(theU, theV, theValue); -} - -void GeomEvaluator_OffsetSurface::BaseD1(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D1(theU, theV, theValue, theD1U, theD1V); - else - myBaseSurf->D1(theU, theV, theValue, theD1U, theD1V); -} - -void GeomEvaluator_OffsetSurface::BaseD2(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - gp_Vec& theD2U, - gp_Vec& theD2V, - gp_Vec& theD2UV) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D2(theU, theV, theValue, theD1U, theD1V, theD2U, theD2V, theD2UV); - else - myBaseSurf->D2(theU, theV, theValue, theD1U, theD1V, theD2U, theD2V, theD2UV); -} - -void GeomEvaluator_OffsetSurface::BaseD3(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - gp_Vec& theD2U, - gp_Vec& theD2V, - gp_Vec& theD2UV, - gp_Vec& theD3U, - gp_Vec& theD3V, - gp_Vec& theD3UUV, - gp_Vec& theD3UVV) const -{ - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D3(theU, - theV, - theValue, - theD1U, - theD1V, - theD2U, - theD2V, - theD2UV, - theD3U, - theD3V, - theD3UUV, - theD3UVV); - else - myBaseSurf->D3(theU, - theV, - theValue, - theD1U, - theD1V, - theD2U, - theD2V, - theD2UV, - theD3U, - theD3V, - theD3UUV, - theD3UVV); -} - -void GeomEvaluator_OffsetSurface::CalculateD0(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - const gp_Vec& theD1U, - const gp_Vec& theD1V) const -{ - // Normalize derivatives before normal calculation because it gives more stable result. - // There will be normalized only derivatives greater than 1.0 to avoid differences in last - // significant digit - gp_Vec aD1U(theD1U); - gp_Vec aD1V(theD1V); - Standard_Real aD1UNorm2 = aD1U.SquareMagnitude(); - Standard_Real aD1VNorm2 = aD1V.SquareMagnitude(); - if (aD1UNorm2 > 1.0) - aD1U /= std::sqrt(aD1UNorm2); - if (aD1VNorm2 > 1.0) - aD1V /= std::sqrt(aD1VNorm2); - - gp_Vec aNorm = aD1U.Crossed(aD1V); - if (aNorm.SquareMagnitude() > the_D1MagTol * the_D1MagTol) - { - // Non singular case. Simple computations. - aNorm.Normalize(); - theValue.SetXYZ(theValue.XYZ() + myOffset * aNorm.XYZ()); - } - else - { - const Standard_Integer MaxOrder = 3; - - Handle(Geom_BSplineSurface) L; - Standard_Boolean isOpposite = Standard_False; - Standard_Boolean AlongU = Standard_False; - Standard_Boolean AlongV = Standard_False; - if (!myOscSurf.IsNull()) - { - AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L); - AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L); - } - const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.; - - TColgp_Array2OfVec DerNUV(0, MaxOrder, 0, MaxOrder); - TColgp_Array2OfVec DerSurf(0, MaxOrder + 1, 0, MaxOrder + 1); - Standard_Integer OrderU, OrderV; - Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0; - Bounds(Umin, Umax, Vmin, Vmax); - - DerSurf.SetValue(1, 0, theD1U); - DerSurf.SetValue(0, 1, theD1V); - if (!myBaseSurf.IsNull()) - derivatives(MaxOrder, 1, theU, theV, myBaseSurf, 0, 0, AlongU, AlongV, L, DerNUV, DerSurf); - else - derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, 0, 0, AlongU, AlongV, L, DerNUV, DerSurf); - - gp_Dir Normal; - CSLib_NormalStatus NStatus = CSLib_Singular; - CSLib::Normal(MaxOrder, - DerNUV, - the_D1MagTol, - theU, - theV, - Umin, - Umax, - Vmin, - Vmax, - NStatus, - Normal, - OrderU, - OrderV); - if (NStatus == CSLib_InfinityOfSolutions) - { - // Replace zero derivative and try to calculate normal - gp_Vec aNewDU = theD1U; - gp_Vec aNewDV = theD1V; - if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, the_D1MagTol * the_D1MagTol)) - CSLib::Normal(aNewDU, aNewDV, the_D1MagTol, NStatus, Normal); - } - - if (NStatus != CSLib_Defined) - throw Geom_UndefinedValue( - "GeomEvaluator_OffsetSurface::CalculateD0(): Unable to calculate normal"); - - theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ()); - } -} - -void GeomEvaluator_OffsetSurface::CalculateD1(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - const gp_Vec& theD2U, - const gp_Vec& theD2V, - const gp_Vec& theD2UV) const -{ - // Check offset side. - Handle(Geom_BSplineSurface) L; - Standard_Boolean isOpposite = Standard_False; - Standard_Boolean AlongU = Standard_False; - Standard_Boolean AlongV = Standard_False; - - // Normalize derivatives before normal calculation because it gives more stable result. - // There will be normalized only derivatives greater than 1.0 to avoid differences in last - // significant digit - gp_Vec aD1U(theD1U); - gp_Vec aD1V(theD1V); - Standard_Real aD1UNorm2 = aD1U.SquareMagnitude(); - Standard_Real aD1VNorm2 = aD1V.SquareMagnitude(); - if (aD1UNorm2 > 1.0) - aD1U /= std::sqrt(aD1UNorm2); - if (aD1VNorm2 > 1.0) - aD1V /= std::sqrt(aD1VNorm2); - - Standard_Boolean isSingular = Standard_False; - Standard_Integer MaxOrder = 0; - gp_Vec aNorm = aD1U.Crossed(aD1V); - if (aNorm.SquareMagnitude() <= the_D1MagTol * the_D1MagTol) - { - if (!myOscSurf.IsNull()) - { - AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L); - AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L); - } - - MaxOrder = 3; - isSingular = Standard_True; - } - - const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.; - - if (!isSingular) - { - // AlongU or AlongV leads to more complex D1 computation - // Try to compute D0 and D1 much simpler - aNorm.Normalize(); - theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * aNorm.XYZ()); - - gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V; - Standard_Real aScale = (theD1U ^ theD1V).Dot(aN0); - aN1U.SetX(theD2U.Y() * theD1V.Z() + theD1U.Y() * theD2UV.Z() - theD2U.Z() * theD1V.Y() - - theD1U.Z() * theD2UV.Y()); - aN1U.SetY((theD2U.X() * theD1V.Z() + theD1U.X() * theD2UV.Z() - theD2U.Z() * theD1V.X() - - theD1U.Z() * theD2UV.X()) - * -1.0); - aN1U.SetZ(theD2U.X() * theD1V.Y() + theD1U.X() * theD2UV.Y() - theD2U.Y() * theD1V.X() - - theD1U.Y() * theD2UV.X()); - Standard_Real aScaleU = aN1U.Dot(aN0); - aN1U.Subtract(aScaleU * aN0); - aN1U /= aScale; - - aN1V.SetX(theD2UV.Y() * theD1V.Z() + theD2V.Z() * theD1U.Y() - theD2UV.Z() * theD1V.Y() - - theD2V.Y() * theD1U.Z()); - aN1V.SetY((theD2UV.X() * theD1V.Z() + theD2V.Z() * theD1U.X() - theD2UV.Z() * theD1V.X() - - theD2V.X() * theD1U.Z()) - * -1.0); - aN1V.SetZ(theD2UV.X() * theD1V.Y() + theD2V.Y() * theD1U.X() - theD2UV.Y() * theD1V.X() - - theD2V.X() * theD1U.Y()); - Standard_Real aScaleV = aN1V.Dot(aN0); - aN1V.Subtract(aScaleV * aN0); - aN1V /= aScale; - - theD1U += myOffset * aSign * aN1U; - theD1V += myOffset * aSign * aN1V; - - return; - } - - Standard_Integer OrderU, OrderV; - TColgp_Array2OfVec DerNUV(0, MaxOrder + 1, 0, MaxOrder + 1); - TColgp_Array2OfVec DerSurf(0, MaxOrder + 2, 0, MaxOrder + 2); - Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0; - Bounds(Umin, Umax, Vmin, Vmax); - - DerSurf.SetValue(1, 0, theD1U); - DerSurf.SetValue(0, 1, theD1V); - DerSurf.SetValue(1, 1, theD2UV); - DerSurf.SetValue(2, 0, theD2U); - DerSurf.SetValue(0, 2, theD2V); - if (!myBaseSurf.IsNull()) - derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf); - else - derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf); - - gp_Dir Normal; - CSLib_NormalStatus NStatus; - CSLib::Normal(MaxOrder, - DerNUV, - the_D1MagTol, - theU, - theV, - Umin, - Umax, - Vmin, - Vmax, - NStatus, - Normal, - OrderU, - OrderV); - if (NStatus == CSLib_InfinityOfSolutions) - { - gp_Vec aNewDU = theD1U; - gp_Vec aNewDV = theD1V; - // Replace zero derivative and try to calculate normal - if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, the_D1MagTol * the_D1MagTol)) - { - DerSurf.SetValue(1, 0, aNewDU); - DerSurf.SetValue(0, 1, aNewDV); - if (!myBaseSurf.IsNull()) - derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf); - else - derivatives(MaxOrder, - 2, - theU, - theV, - myBaseAdaptor, - 1, - 1, - AlongU, - AlongV, - L, - DerNUV, - DerSurf); - CSLib::Normal(MaxOrder, - DerNUV, - the_D1MagTol, - theU, - theV, - Umin, - Umax, - Vmin, - Vmax, - NStatus, - Normal, - OrderU, - OrderV); - } - } - - if (NStatus != CSLib_Defined) - throw Geom_UndefinedValue( - "GeomEvaluator_OffsetSurface::CalculateD1(): Unable to calculate normal"); - - theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ()); - - theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV); - theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV); -} - -void GeomEvaluator_OffsetSurface::CalculateD2(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - gp_Vec& theD2U, - gp_Vec& theD2V, - gp_Vec& theD2UV, - const gp_Vec& theD3U, - const gp_Vec& theD3V, - const gp_Vec& theD3UUV, - const gp_Vec& theD3UVV) const -{ - gp_Dir Normal; - CSLib_NormalStatus NStatus; - CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal); - - const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3; - Standard_Integer OrderU, OrderV; - TColgp_Array2OfVec DerNUV(0, MaxOrder + 2, 0, MaxOrder + 2); - TColgp_Array2OfVec DerSurf(0, MaxOrder + 3, 0, MaxOrder + 3); - - Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0; - Bounds(Umin, Umax, Vmin, Vmax); - - DerSurf.SetValue(1, 0, theD1U); - DerSurf.SetValue(0, 1, theD1V); - DerSurf.SetValue(1, 1, theD2UV); - DerSurf.SetValue(2, 0, theD2U); - DerSurf.SetValue(0, 2, theD2V); - DerSurf.SetValue(3, 0, theD3U); - DerSurf.SetValue(2, 1, theD3UUV); - DerSurf.SetValue(1, 2, theD3UVV); - DerSurf.SetValue(0, 3, theD3V); - //********************* - - Handle(Geom_BSplineSurface) L; - Standard_Boolean isOpposite = Standard_False; - Standard_Boolean AlongU = Standard_False; - Standard_Boolean AlongV = Standard_False; - if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull()) - { - AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L); - AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L); - } - const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.; - - if (!myBaseSurf.IsNull()) - derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf); - else - derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf); - - CSLib::Normal(MaxOrder, - DerNUV, - the_D1MagTol, - theU, - theV, - Umin, - Umax, - Vmin, - Vmax, - NStatus, - Normal, - OrderU, - OrderV); - if (NStatus != CSLib_Defined) - throw Geom_UndefinedValue( - "GeomEvaluator_OffsetSurface::CalculateD2(): Unable to calculate normal"); - - theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ()); - - theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV); - theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV); - - if (!myBaseSurf.IsNull()) - { - theD2U = myBaseSurf->DN(theU, theV, 2, 0); - theD2V = myBaseSurf->DN(theU, theV, 0, 2); - theD2UV = myBaseSurf->DN(theU, theV, 1, 1); - } - else - { - theD2U = myBaseAdaptor->DN(theU, theV, 2, 0); - theD2V = myBaseAdaptor->DN(theU, theV, 0, 2); - theD2UV = myBaseAdaptor->DN(theU, theV, 1, 1); - } - - theD2U += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV); - theD2V += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV); - theD2UV += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV); -} - -void GeomEvaluator_OffsetSurface::CalculateD3(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - gp_Vec& theD2U, - gp_Vec& theD2V, - gp_Vec& theD2UV, - gp_Vec& theD3U, - gp_Vec& theD3V, - gp_Vec& theD3UUV, - gp_Vec& theD3UVV) const -{ - gp_Dir Normal; - CSLib_NormalStatus NStatus; - CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal); - const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3; - Standard_Integer OrderU, OrderV; - TColgp_Array2OfVec DerNUV(0, MaxOrder + 3, 0, MaxOrder + 3); - TColgp_Array2OfVec DerSurf(0, MaxOrder + 4, 0, MaxOrder + 4); - Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0; - Bounds(Umin, Umax, Vmin, Vmax); - - DerSurf.SetValue(1, 0, theD1U); - DerSurf.SetValue(0, 1, theD1V); - DerSurf.SetValue(1, 1, theD2UV); - DerSurf.SetValue(2, 0, theD2U); - DerSurf.SetValue(0, 2, theD2V); - DerSurf.SetValue(3, 0, theD3U); - DerSurf.SetValue(2, 1, theD3UUV); - DerSurf.SetValue(1, 2, theD3UVV); - DerSurf.SetValue(0, 3, theD3V); - - //********************* - Handle(Geom_BSplineSurface) L; - Standard_Boolean isOpposite = Standard_False; - Standard_Boolean AlongU = Standard_False; - Standard_Boolean AlongV = Standard_False; - if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull()) - { - AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L); - AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L); - } - const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.; - - if (!myBaseSurf.IsNull()) - derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf); - else - derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf); - - CSLib::Normal(MaxOrder, - DerNUV, - the_D1MagTol, - theU, - theV, - Umin, - Umax, - Vmin, - Vmax, - NStatus, - Normal, - OrderU, - OrderV); - if (NStatus != CSLib_Defined) - throw Geom_UndefinedValue( - "GeomEvaluator_OffsetSurface::CalculateD3(): Unable to calculate normal"); - - theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ()); - - theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV); - theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV); - - if (!myBaseSurf.IsNull()) - { - theD2U = myBaseSurf->DN(theU, theV, 2, 0); - theD2V = myBaseSurf->DN(theU, theV, 0, 2); - theD2UV = myBaseSurf->DN(theU, theV, 1, 1); - theD3U = myBaseSurf->DN(theU, theV, 3, 0); - theD3V = myBaseSurf->DN(theU, theV, 0, 3); - theD3UUV = myBaseSurf->DN(theU, theV, 2, 1); - theD3UVV = myBaseSurf->DN(theU, theV, 1, 2); - } - else - { - theD2U = myBaseAdaptor->DN(theU, theV, 2, 0); - theD2V = myBaseAdaptor->DN(theU, theV, 0, 2); - theD2UV = myBaseAdaptor->DN(theU, theV, 1, 1); - theD3U = myBaseAdaptor->DN(theU, theV, 3, 0); - theD3V = myBaseAdaptor->DN(theU, theV, 0, 3); - theD3UUV = myBaseAdaptor->DN(theU, theV, 2, 1); - theD3UVV = myBaseAdaptor->DN(theU, theV, 1, 2); - } - - theD2U += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV); - theD2V += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV); - theD2UV += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV); - theD3U += aSign * myOffset * CSLib::DNNormal(3, 0, DerNUV, OrderU, OrderV); - theD3V += aSign * myOffset * CSLib::DNNormal(0, 3, DerNUV, OrderU, OrderV); - theD3UUV += aSign * myOffset * CSLib::DNNormal(2, 1, DerNUV, OrderU, OrderV); - theD3UVV += aSign * myOffset * CSLib::DNNormal(1, 2, DerNUV, OrderU, OrderV); -} - -gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(const Standard_Real theU, - const Standard_Real theV, - const Standard_Integer theNu, - const Standard_Integer theNv, - const gp_Vec& theD1U, - const gp_Vec& theD1V) const -{ - gp_Dir Normal; - CSLib_NormalStatus NStatus; - CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal); - const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3; - Standard_Integer OrderU, OrderV; - TColgp_Array2OfVec DerNUV(0, MaxOrder + theNu, 0, MaxOrder + theNv); - TColgp_Array2OfVec DerSurf(0, MaxOrder + theNu + 1, 0, MaxOrder + theNv + 1); - - Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0; - Bounds(Umin, Umax, Vmin, Vmax); - - DerSurf.SetValue(1, 0, theD1U); - DerSurf.SetValue(0, 1, theD1V); - - //********************* - Handle(Geom_BSplineSurface) L; - // Is there any osculatingsurface along U or V; - Standard_Boolean isOpposite = Standard_False; - Standard_Boolean AlongU = Standard_False; - Standard_Boolean AlongV = Standard_False; - if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull()) - { - AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L); - AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L); - } - const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.; - - if (!myBaseSurf.IsNull()) - derivatives(MaxOrder, - 1, - theU, - theV, - myBaseSurf, - theNu, - theNv, - AlongU, - AlongV, - L, - DerNUV, - DerSurf); - else - derivatives(MaxOrder, - 1, - theU, - theV, - myBaseAdaptor, - theNu, - theNv, - AlongU, - AlongV, - L, - DerNUV, - DerSurf); - - CSLib::Normal(MaxOrder, - DerNUV, - the_D1MagTol, - theU, - theV, - Umin, - Umax, - Vmin, - Vmax, - NStatus, - Normal, - OrderU, - OrderV); - if (NStatus != CSLib_Defined) - throw Geom_UndefinedValue( - "GeomEvaluator_OffsetSurface::CalculateDN(): Unable to calculate normal"); - - gp_Vec D; - if (!myBaseSurf.IsNull()) - D = myBaseSurf->DN(theU, theV, theNu, theNv); - else - D = myBaseAdaptor->DN(theU, theV, theNu, theNv); - - D += aSign * myOffset * CSLib::DNNormal(theNu, theNv, DerNUV, OrderU, OrderV); - return D; -} - -void GeomEvaluator_OffsetSurface::Bounds(Standard_Real& theUMin, - Standard_Real& theUMax, - Standard_Real& theVMin, - Standard_Real& theVMax) const -{ - if (!myBaseSurf.IsNull()) - myBaseSurf->Bounds(theUMin, theUMax, theVMin, theVMax); - else - { - theUMin = myBaseAdaptor->FirstUParameter(); - theUMax = myBaseAdaptor->LastUParameter(); - theVMin = myBaseAdaptor->FirstVParameter(); - theVMax = myBaseAdaptor->LastVParameter(); - } -} - -Standard_Boolean GeomEvaluator_OffsetSurface::ReplaceDerivative( - const Standard_Real theU, - const Standard_Real theV, - gp_Vec& theDU, - gp_Vec& theDV, - const Standard_Real theSquareTol) const -{ - Standard_Boolean isReplaceDU = theDU.SquareMagnitude() < theSquareTol; - Standard_Boolean isReplaceDV = theDV.SquareMagnitude() < theSquareTol; - Standard_Boolean isReplaced = Standard_False; - if (isReplaceDU ^ isReplaceDV) - { - Standard_Real aU = theU; - Standard_Real aV = theV; - Standard_Real aUMin = 0, aUMax = 0, aVMin = 0, aVMax = 0; - Bounds(aUMin, aUMax, aVMin, aVMax); - - // Calculate step along non-zero derivative - Standard_Real aStep; - Handle(Adaptor3d_Surface) aSurfAdapt; - if (!myBaseAdaptor.IsNull()) - aSurfAdapt = myBaseAdaptor; - else - aSurfAdapt = new GeomAdaptor_Surface(myBaseSurf); - if (isReplaceDV) - { - aStep = Precision::Confusion() * theDU.Magnitude(); - if (aStep > aUMax - aUMin) - aStep = (aUMax - aUMin) / 100.; - } - else - { - aStep = Precision::Confusion() * theDV.Magnitude(); - if (aStep > aVMax - aVMin) - aStep = (aVMax - aVMin) / 100.; - } - - gp_Pnt aP; - gp_Vec aDU, aDV; - // Step away from current parametric coordinates and calculate derivatives once again. - // Replace zero derivative by the obtained. - for (Standard_Real aStepSign = -1.0; aStepSign <= 1.0 && !isReplaced; aStepSign += 2.0) - { - if (isReplaceDV) - { - aU = theU + aStepSign * aStep; - if (aU < aUMin || aU > aUMax) - continue; - } - else - { - aV = theV + aStepSign * aStep; - if (aV < aVMin || aV > aVMax) - continue; - } - - BaseD1(aU, aV, aP, aDU, aDV); - - if (isReplaceDU && aDU.SquareMagnitude() > theSquareTol) - { - theDU = aDU; - isReplaced = Standard_True; - } - if (isReplaceDV && aDV.SquareMagnitude() > theSquareTol) - { - theDV = aDV; - isReplaced = Standard_True; - } - } - } - return isReplaced; -} diff --git a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.hxx b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.hxx index 0bce0b3793..1b1376b6e1 100644 --- a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.hxx +++ b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.hxx @@ -82,107 +82,12 @@ public: DEFINE_STANDARD_RTTIEXT(GeomEvaluator_OffsetSurface, GeomEvaluator_Surface) -private: - //! Returns bounds of a base surface - void Bounds(Standard_Real& theUMin, - Standard_Real& theUMax, - Standard_Real& theVMin, - Standard_Real& theVMax) const; - - //! Recalculate D1 values of base surface into D0 value of offset surface - void CalculateD0(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - const gp_Vec& theD1U, - const gp_Vec& theD1V) const; - //! Recalculate D2 values of base surface into D1 values of offset surface - void CalculateD1(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - const gp_Vec& theD2U, - const gp_Vec& theD2V, - const gp_Vec& theD2UV) const; - //! Recalculate D3 values of base surface into D2 values of offset surface - void CalculateD2(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - gp_Vec& theD2U, - gp_Vec& theD2V, - gp_Vec& theD2UV, - const gp_Vec& theD3U, - const gp_Vec& theD3V, - const gp_Vec& theD3UUV, - const gp_Vec& theD3UVV) const; - //! Recalculate D3 values of base surface into D3 values of offset surface - void CalculateD3(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - gp_Vec& theD2U, - gp_Vec& theD2V, - gp_Vec& theD2UV, - gp_Vec& theD3U, - gp_Vec& theD3V, - gp_Vec& theD3UUV, - gp_Vec& theD3UVV) const; - //! Calculate DN of offset surface based on derivatives of base surface - gp_Vec CalculateDN(const Standard_Real theU, - const Standard_Real theV, - const Standard_Integer theNu, - const Standard_Integer theNv, - const gp_Vec& theD1U, - const gp_Vec& theD1V) const; - - //! Calculate value of base surface/adaptor - void BaseD0(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue) const; - //! Calculate value and first derivatives of base surface/adaptor - void BaseD1(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V) const; - //! Calculate value, first and second derivatives of base surface/adaptor - void BaseD2(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - gp_Vec& theD2U, - gp_Vec& theD2V, - gp_Vec& theD2UV) const; - //! Calculate value, first, second and third derivatives of base surface/adaptor - void BaseD3(const Standard_Real theU, - const Standard_Real theV, - gp_Pnt& theValue, - gp_Vec& theD1U, - gp_Vec& theD1V, - gp_Vec& theD2U, - gp_Vec& theD2V, - gp_Vec& theD2UV, - gp_Vec& theD3U, - gp_Vec& theD3V, - gp_Vec& theD3UUV, - gp_Vec& theD3UVV) const; - - //! Replace zero derivative by the corresponding derivative in a near point. - //! Return true, if the derivative was replaced. - Standard_Boolean ReplaceDerivative(const Standard_Real theU, - const Standard_Real theV, - gp_Vec& theDU, - gp_Vec& theDV, - const Standard_Real theSquareTol) const; - private: Handle(Geom_Surface) myBaseSurf; Handle(GeomAdaptor_Surface) myBaseAdaptor; - Standard_Real myOffset; ///< offset value - Handle(Geom_OsculatingSurface) myOscSurf; ///< auxiliary osculating surface + Standard_Real myOffset; //!< offset value + Handle(Geom_OsculatingSurface) myOscSurf; //!< auxiliary osculating surface }; DEFINE_STANDARD_HANDLE(GeomEvaluator_OffsetSurface, GeomEvaluator_Surface) diff --git a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.cxx b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.cxx index fa51c2b89d..4d8b7e6d6c 100644 --- a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.cxx +++ b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.cxx @@ -15,6 +15,7 @@ #include #include +#include IMPLEMENT_STANDARD_RTTIEXT(GeomEvaluator_SurfaceOfExtrusion, GeomEvaluator_Surface) @@ -36,33 +37,38 @@ GeomEvaluator_SurfaceOfExtrusion::GeomEvaluator_SurfaceOfExtrusion( { } +//================================================================================================== + void GeomEvaluator_SurfaceOfExtrusion::D0(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue) const { - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D0(theU, theValue); - else - myBaseCurve->D0(theU, theValue); + const gp_XYZ& aDir = myDirection.XYZ(); - Shift(theV, theValue); + if (!myBaseAdaptor.IsNull()) + Geom_ExtrusionUtils::D0(theU, theV, *myBaseAdaptor, aDir, theValue); + else + Geom_ExtrusionUtils::D0(theU, theV, *myBaseCurve, aDir, theValue); } +//================================================================================================== + void GeomEvaluator_SurfaceOfExtrusion::D1(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const { - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D1(theU, theValue, theD1U); - else - myBaseCurve->D1(theU, theValue, theD1U); + const gp_XYZ& aDir = myDirection.XYZ(); - theD1V = myDirection; - Shift(theV, theValue); + if (!myBaseAdaptor.IsNull()) + Geom_ExtrusionUtils::D1(theU, theV, *myBaseAdaptor, aDir, theValue, theD1U, theD1V); + else + Geom_ExtrusionUtils::D1(theU, theV, *myBaseCurve, aDir, theValue, theD1U, theD1V); } +//================================================================================================== + void GeomEvaluator_SurfaceOfExtrusion::D2(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue, @@ -72,18 +78,34 @@ void GeomEvaluator_SurfaceOfExtrusion::D2(const Standard_Real theU, gp_Vec& theD2V, gp_Vec& theD2UV) const { + const gp_XYZ& aDir = myDirection.XYZ(); + if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D2(theU, theValue, theD1U, theD2U); + Geom_ExtrusionUtils::D2(theU, + theV, + *myBaseAdaptor, + aDir, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV); else - myBaseCurve->D2(theU, theValue, theD1U, theD2U); - - theD1V = myDirection; - theD2V.SetCoord(0.0, 0.0, 0.0); - theD2UV.SetCoord(0.0, 0.0, 0.0); - - Shift(theV, theValue); + Geom_ExtrusionUtils::D2(theU, + theV, + *myBaseCurve, + aDir, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV); } +//================================================================================================== + void GeomEvaluator_SurfaceOfExtrusion::D3(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue, @@ -97,21 +119,42 @@ void GeomEvaluator_SurfaceOfExtrusion::D3(const Standard_Real theU, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const { + const gp_XYZ& aDir = myDirection.XYZ(); + if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D3(theU, theValue, theD1U, theD2U, theD3U); + Geom_ExtrusionUtils::D3(theU, + theV, + *myBaseAdaptor, + aDir, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV, + theD3U, + theD3V, + theD3UUV, + theD3UVV); else - myBaseCurve->D3(theU, theValue, theD1U, theD2U, theD3U); - - theD1V = myDirection; - theD2V.SetCoord(0.0, 0.0, 0.0); - theD2UV.SetCoord(0.0, 0.0, 0.0); - theD3V.SetCoord(0.0, 0.0, 0.0); - theD3UUV.SetCoord(0.0, 0.0, 0.0); - theD3UVV.SetCoord(0.0, 0.0, 0.0); - - Shift(theV, theValue); + Geom_ExtrusionUtils::D3(theU, + theV, + *myBaseCurve, + aDir, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV, + theD3U, + theD3V, + theD3UUV, + theD3UVV); } +//================================================================================================== + gp_Vec GeomEvaluator_SurfaceOfExtrusion::DN(const Standard_Real theU, const Standard_Real, const Standard_Integer theDerU, @@ -122,19 +165,16 @@ gp_Vec GeomEvaluator_SurfaceOfExtrusion::DN(const Standard_Real theU, Standard_RangeError_Raise_if(theDerU + theDerV < 1, "GeomEvaluator_SurfaceOfExtrusion::DN(): theDerU + theDerV < 1"); - gp_Vec aResult(0.0, 0.0, 0.0); - if (theDerV == 0) - { - if (!myBaseAdaptor.IsNull()) - aResult = myBaseAdaptor->DN(theU, theDerU); - else - aResult = myBaseCurve->DN(theU, theDerU); - } - else if (theDerU == 0 && theDerV == 1) - aResult = gp_Vec(myDirection); - return aResult; + const gp_XYZ& aDir = myDirection.XYZ(); + + if (!myBaseAdaptor.IsNull()) + return Geom_ExtrusionUtils::DN(theU, *myBaseAdaptor, aDir, theDerU, theDerV); + else + return Geom_ExtrusionUtils::DN(theU, *myBaseCurve, aDir, theDerU, theDerV); } +//================================================================================================== + Handle(GeomEvaluator_Surface) GeomEvaluator_SurfaceOfExtrusion::ShallowCopy() const { Handle(GeomEvaluator_SurfaceOfExtrusion) aCopy; diff --git a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.hxx b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.hxx index 3cf45007b0..105b0c17ed 100644 --- a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.hxx +++ b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.hxx @@ -78,13 +78,6 @@ public: DEFINE_STANDARD_RTTIEXT(GeomEvaluator_SurfaceOfExtrusion, GeomEvaluator_Surface) -private: - //! Shift the point along direction to the given distance (theShift) - void Shift(const Standard_Real theShift, gp_Pnt& thePoint) const - { - thePoint.ChangeCoord() += myDirection.XYZ() * theShift; - } - private: Handle(Geom_Curve) myBaseCurve; Handle(Adaptor3d_Curve) myBaseAdaptor; diff --git a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfRevolution.cxx b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfRevolution.cxx index 1b42d347f8..3d5d47b43e 100644 --- a/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfRevolution.cxx +++ b/src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfRevolution.cxx @@ -15,8 +15,7 @@ #include #include -#include -#include +#include IMPLEMENT_STANDARD_RTTIEXT(GeomEvaluator_SurfaceOfRevolution, GeomEvaluator_Surface) @@ -40,20 +39,20 @@ GeomEvaluator_SurfaceOfRevolution::GeomEvaluator_SurfaceOfRevolution( { } +//================================================================================================== + void GeomEvaluator_SurfaceOfRevolution::D0(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue) const { if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D0(theV, theValue); + Geom_RevolutionUtils::D0(theU, theV, *myBaseAdaptor, myRotAxis, theValue); else - myBaseCurve->D0(theV, theValue); - - gp_Trsf aRotation; - aRotation.SetRotation(myRotAxis, theU); - theValue.Transform(aRotation); + Geom_RevolutionUtils::D0(theU, theV, *myBaseCurve, myRotAxis, theValue); } +//================================================================================================== + void GeomEvaluator_SurfaceOfRevolution::D1(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue, @@ -61,25 +60,13 @@ void GeomEvaluator_SurfaceOfRevolution::D1(const Standard_Real theU, gp_Vec& theD1V) const { if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D1(theV, theValue, theD1V); + Geom_RevolutionUtils::D1(theU, theV, *myBaseAdaptor, myRotAxis, theValue, theD1U, theD1V); else - myBaseCurve->D1(theV, theValue, theD1V); - - // vector from center of rotation to the point on rotated curve - gp_XYZ aCQ = theValue.XYZ() - myRotAxis.Location().XYZ(); - theD1U = gp_Vec(myRotAxis.Direction().XYZ().Crossed(aCQ)); - // If the point is placed on the axis of revolution then derivatives on U are undefined. - // Manually set them to zero. - if (theD1U.SquareMagnitude() < Precision::SquareConfusion()) - theD1U.SetCoord(0.0, 0.0, 0.0); - - gp_Trsf aRotation; - aRotation.SetRotation(myRotAxis, theU); - theValue.Transform(aRotation); - theD1U.Transform(aRotation); - theD1V.Transform(aRotation); + Geom_RevolutionUtils::D1(theU, theV, *myBaseCurve, myRotAxis, theValue, theD1U, theD1V); } +//================================================================================================== + void GeomEvaluator_SurfaceOfRevolution::D2(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue, @@ -90,31 +77,31 @@ void GeomEvaluator_SurfaceOfRevolution::D2(const Standard_Real theU, gp_Vec& theD2UV) const { if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D2(theV, theValue, theD1V, theD2V); + Geom_RevolutionUtils::D2(theU, + theV, + *myBaseAdaptor, + myRotAxis, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV); else - myBaseCurve->D2(theV, theValue, theD1V, theD2V); - - // vector from center of rotation to the point on rotated curve - gp_XYZ aCQ = theValue.XYZ() - myRotAxis.Location().XYZ(); - const gp_XYZ& aDir = myRotAxis.Direction().XYZ(); - theD1U = gp_Vec(aDir.Crossed(aCQ)); - // If the point is placed on the axis of revolution then derivatives on U are undefined. - // Manually set them to zero. - if (theD1U.SquareMagnitude() < Precision::SquareConfusion()) - theD1U.SetCoord(0.0, 0.0, 0.0); - theD2U = gp_Vec(aDir.Dot(aCQ) * aDir - aCQ); - theD2UV = gp_Vec(aDir.Crossed(theD1V.XYZ())); - - gp_Trsf aRotation; - aRotation.SetRotation(myRotAxis, theU); - theValue.Transform(aRotation); - theD1U.Transform(aRotation); - theD1V.Transform(aRotation); - theD2U.Transform(aRotation); - theD2V.Transform(aRotation); - theD2UV.Transform(aRotation); + Geom_RevolutionUtils::D2(theU, + theV, + *myBaseCurve, + myRotAxis, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV); } +//================================================================================================== + void GeomEvaluator_SurfaceOfRevolution::D3(const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue, @@ -129,38 +116,39 @@ void GeomEvaluator_SurfaceOfRevolution::D3(const Standard_Real theU, gp_Vec& theD3UVV) const { if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D3(theV, theValue, theD1V, theD2V, theD3V); + Geom_RevolutionUtils::D3(theU, + theV, + *myBaseAdaptor, + myRotAxis, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV, + theD3U, + theD3V, + theD3UUV, + theD3UVV); else - myBaseCurve->D3(theV, theValue, theD1V, theD2V, theD3V); - - // vector from center of rotation to the point on rotated curve - gp_XYZ aCQ = theValue.XYZ() - myRotAxis.Location().XYZ(); - const gp_XYZ& aDir = myRotAxis.Direction().XYZ(); - theD1U = gp_Vec(aDir.Crossed(aCQ)); - // If the point is placed on the axis of revolution then derivatives on U are undefined. - // Manually set them to zero. - if (theD1U.SquareMagnitude() < Precision::SquareConfusion()) - theD1U.SetCoord(0.0, 0.0, 0.0); - theD2U = gp_Vec(aDir.Dot(aCQ) * aDir - aCQ); - theD2UV = gp_Vec(aDir.Crossed(theD1V.XYZ())); - theD3U = -theD1U; - theD3UUV = gp_Vec(aDir.Dot(theD1V.XYZ()) * aDir - theD1V.XYZ()); - theD3UVV = gp_Vec(aDir.Crossed(theD2V.XYZ())); - - gp_Trsf aRotation; - aRotation.SetRotation(myRotAxis, theU); - theValue.Transform(aRotation); - theD1U.Transform(aRotation); - theD1V.Transform(aRotation); - theD2U.Transform(aRotation); - theD2V.Transform(aRotation); - theD2UV.Transform(aRotation); - theD3U.Transform(aRotation); - theD3V.Transform(aRotation); - theD3UUV.Transform(aRotation); - theD3UVV.Transform(aRotation); + Geom_RevolutionUtils::D3(theU, + theV, + *myBaseCurve, + myRotAxis, + theValue, + theD1U, + theD1V, + theD2U, + theD2V, + theD2UV, + theD3U, + theD3V, + theD3UUV, + theD3UVV); } +//================================================================================================== + gp_Vec GeomEvaluator_SurfaceOfRevolution::DN(const Standard_Real theU, const Standard_Real theV, const Standard_Integer theDerU, @@ -171,52 +159,14 @@ gp_Vec GeomEvaluator_SurfaceOfRevolution::DN(const Standard_Real theU, Standard_RangeError_Raise_if(theDerU + theDerV < 1, "GeomEvaluator_SurfaceOfRevolution::DN(): theDerU + theDerV < 1"); - gp_Trsf aRotation; - aRotation.SetRotation(myRotAxis, theU); - - gp_Pnt aP; - gp_Vec aDV; - gp_Vec aResult; - if (theDerU == 0) - { - if (!myBaseAdaptor.IsNull()) - aResult = myBaseAdaptor->DN(theV, theDerV); - else - aResult = myBaseCurve->DN(theV, theDerV); - } + if (!myBaseAdaptor.IsNull()) + return Geom_RevolutionUtils::DN(theU, theV, *myBaseAdaptor, myRotAxis, theDerU, theDerV); else - { - if (theDerV == 0) - { - if (!myBaseAdaptor.IsNull()) - myBaseAdaptor->D0(theV, aP); - else - myBaseCurve->D0(theV, aP); - aDV = gp_Vec(aP.XYZ() - myRotAxis.Location().XYZ()); - } - else - { - if (!myBaseAdaptor.IsNull()) - aDV = myBaseAdaptor->DN(theV, theDerV); - else - aDV = myBaseCurve->DN(theV, theDerV); - } - - const gp_XYZ& aDir = myRotAxis.Direction().XYZ(); - if (theDerU % 4 == 1) - aResult = gp_Vec(aDir.Crossed(aDV.XYZ())); - else if (theDerU % 4 == 2) - aResult = gp_Vec(aDir.Dot(aDV.XYZ()) * aDir - aDV.XYZ()); - else if (theDerU % 4 == 3) - aResult = gp_Vec(aDir.Crossed(aDV.XYZ())) * (-1.0); - else - aResult = gp_Vec(aDV.XYZ() - aDir.Dot(aDV.XYZ()) * aDir); - } - - aResult.Transform(aRotation); - return aResult; + return Geom_RevolutionUtils::DN(theU, theV, *myBaseCurve, myRotAxis, theDerU, theDerV); } +//================================================================================================== + Handle(GeomEvaluator_Surface) GeomEvaluator_SurfaceOfRevolution::ShallowCopy() const { Handle(GeomEvaluator_SurfaceOfRevolution) aCopy; diff --git a/tests/bugs/modalg_7/bug30679 b/tests/bugs/modalg_7/bug30679 index 0461101753..25418eaff6 100644 --- a/tests/bugs/modalg_7/bug30679 +++ b/tests/bugs/modalg_7/bug30679 @@ -3,7 +3,7 @@ puts "0030679: Attached model hangs most of OCCT common functionality" puts "========" puts "" -puts "REQUIRED ALL: Evaluation of infinite parameters" +puts "REQUIRED ALL: Unable to calculate normal" restore [locate_data_file bug30679_face.brep] a pcurve a