lili/OCCT
1
0
Fork 0
mirror of https://github.com/Open-Cascade-SAS/OCCT.git synced 2026-07-18 04:36:41 +08:00
Code Issues Packages Projects Releases Wiki Activity

Modeling Data - Refactor offset curve and surface evaluators (#930)

Browse Source 
Extract derivative calculation logic from GeomEvaluator and
  Geom2dEvaluator classes into dedicated utility headers:
- Geom_OffsetSurfaceUtils.pxx for 3D offset surface derivatives
- Geom_OffsetCurveUtils.pxx for 3D offset curve derivatives  
- Geom2d_OffsetCurveUtils.pxx for 2D offset curve derivatives
- Geom_ExtrusionUtils.pxx for extrusion surface calculations
- Geom_RevolutionUtils.pxx for revolution surface calculations

Utility functions provide unified handling of singular and non-singular
  cases with osculating surface support for higher-order derivatives.

Buffer management uses NCollection_LocalArray for stack-safe allocation
  with heap fallback when sizes exceed expected limits.

Test case expectations updated to reflect improved calculation accuracy.
This commit is contained in:
Pasukhin Dmitry
2025-12-16 21:06:40 +00:00
committed by GitHub
parent e268a55440
commit 4727f7c049
18 changed files with 3506 additions and 1928 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/ModelingData/TKG2d/Geom2d/FILES.cmake
View File

@@ -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

507
src/ModelingData/TKG2d/Geom2d/Geom2d_OffsetCurveUtils.pxx Normal file
View File

@@ -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 <gp.hxx>
#include <gp_Dir2d.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Vec2d.hxx>
#include <gp_XY.hxx>
#include <cmath>
//! 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 <typename CurveType>
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 <typename BasisCurveType>
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 <typename BasisCurveType>
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 <typename BasisCurveType>
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 <typename BasisCurveType>
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 <typename BasisCurveType>
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

268
src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.cxx
View File

@@ -13,9 +13,10 @@
// commercial license or contractual agreement.
#include <Geom2dEvaluator_OffsetCurve.hxx>
#include <Geom2dEvaluator.hxx>
#include <Geom2d_OffsetCurveUtils.pxx>
#include <Geom2d_UndefinedValue.hxx>
#include <Geom2dAdaptor_Curve.hxx>
#include <Standard_NullValue.hxx>
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;
}

32
src/ModelingData/TKG2d/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.hxx
View File

@@ -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;

4
src/ModelingData/TKG3d/Geom/FILES.cmake
View File

@@ -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

172
src/ModelingData/TKG3d/Geom/Geom_ExtrusionUtils.pxx Normal file
View File

@@ -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 <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <gp_XYZ.hxx>
//! @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 <typename CurveType>
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 <typename CurveType>
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 <typename CurveType>
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 <typename CurveType>
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 <typename CurveType>
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

532
src/ModelingData/TKG3d/Geom/Geom_OffsetCurveUtils.pxx Normal file
View File

@@ -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 <gp.hxx>
#include <gp_Dir.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <gp_XYZ.hxx>
#include <Standard_NullValue.hxx>
#include <cmath>
//! 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 <typename CurveType>
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 <typename BasisCurveType>
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 <typename BasisCurveType>
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 <typename BasisCurveType>
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 <typename BasisCurveType>
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 <typename BasisCurveType>
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

1500
src/ModelingData/TKG3d/Geom/Geom_OffsetSurfaceUtils.pxx Normal file
View File

File diff suppressed because it is too large Load Diff

272
src/ModelingData/TKG3d/Geom/Geom_RevolutionUtils.pxx Normal file
View File

@@ -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 <gp_Ax1.hxx>
#include <gp_Pnt.hxx>
#include <gp_Trsf.hxx>
#include <gp_Vec.hxx>
#include <Precision.hxx>
//! @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 <typename CurveType>
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 <typename CurveType>
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 <typename CurveType>
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 <typename CurveType>
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 <typename CurveType>
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

16
src/ModelingData/TKG3d/GeomAdaptor/GeomAdaptor_Surface.hxx
View File

@@ -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; }

467
src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.cxx
View File

@@ -15,6 +15,7 @@
#include <GeomEvaluator_OffsetCurve.hxx>
#include <GeomAdaptor_Curve.hxx>
#include <Geom_OffsetCurveUtils.pxx>
#include <Standard_NullValue.hxx>
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;
}

54
src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetCurve.hxx
View File

@@ -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)

1193
src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx
View File

File diff suppressed because it is too large Load Diff

99
src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_OffsetSurface.hxx
View File

@@ -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)

122
src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.cxx
View File

@@ -15,6 +15,7 @@
#include <GeomEvaluator_SurfaceOfExtrusion.hxx>
#include <GeomAdaptor_Curve.hxx>
#include <Geom_ExtrusionUtils.pxx>
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;

7
src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfExtrusion.hxx
View File

@@ -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;

186
src/ModelingData/TKG3d/GeomEvaluator/GeomEvaluator_SurfaceOfRevolution.cxx
View File

@@ -15,8 +15,7 @@
#include <GeomEvaluator_SurfaceOfRevolution.hxx>
#include <Adaptor3d_Curve.hxx>
#include <gp_Trsf.hxx>
#include <Precision.hxx>
#include <Geom_RevolutionUtils.pxx>
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;

2
tests/bugs/modalg_7/bug30679
View File

@@ -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