From: azv Date: Fri, 6 Nov 2015 05:49:50 +0000 (+0300) Subject: 0026838: Using GeomEvaluators for calculation of values of curves X-Git-Tag: V7_0_0_beta~87 X-Git-Url: http://git.dev.opencascade.org/gitweb/?p=occt.git;a=commitdiff_plain;h=d660a72aca661bce175a28f57fdd6c9f9b8b2837 0026838: Using GeomEvaluators for calculation of values of curves 1. Implemented evaluators for 2D and 3D offset curves 2. Removed obsolete namespace CSLib_Offset Update of UDLIST (adding no-cdl-pack Geom2dEvaluator) Update TKG2d/CMakeLists.txt after rebase Correction compilation in debug mode --- diff --git a/adm/UDLIST b/adm/UDLIST index 0e596d3014..027f83c73d 100644 --- a/adm/UDLIST +++ b/adm/UDLIST @@ -431,3 +431,4 @@ n IVtkTools t TKIVtk n IVtkDraw t TKIVtkDraw +n Geom2dEvaluator diff --git a/src/CSLib/CSLib_Offset.cxx b/src/CSLib/CSLib_Offset.cxx deleted file mode 100644 index fa8fce47bc..0000000000 --- a/src/CSLib/CSLib_Offset.cxx +++ /dev/null @@ -1,465 +0,0 @@ -// Copyright (c) 2015-... OPEN CASCADE SAS -// -// This file is part of Open CASCADE Technology software library. -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License version 2.1 as published -// by the Free Software Foundation, with special exception defined in the file -// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT -// distribution for complete text of the license and disclaimer of any warranty. -// -// Alternatively, this file may be used under the terms of Open CASCADE -// commercial license or contractual agreement. - -#include - -#include -#include - -#include - - -// ========== Offset values for 2D curves ========== - -void CSLib_Offset::D0(const gp_Pnt2d& theBasePoint, - const gp_Vec2d& theBaseDeriv, - Standard_Real theOffset, - Standard_Boolean , // unused - gp_Pnt2d& theResPoint) -{ - if (theBaseDeriv.SquareMagnitude() <= gp::Resolution()) - Standard_NullValue::Raise("CSLib_Offset: Undefined normal vector " - "because tangent vector has zero-magnitude!"); - - gp_Dir2d aNormal(theBaseDeriv.Y(), -theBaseDeriv.X()); - theResPoint.SetCoord(theBasePoint.X() + aNormal.X() * theOffset, - theBasePoint.Y() + aNormal.Y() * theOffset); -} - -void CSLib_Offset::D1(const gp_Pnt2d& theBasePoint, - const gp_Vec2d& theBaseD1, - const gp_Vec2d& theBaseD2, - Standard_Real theOffset, - Standard_Boolean , // unused - gp_Pnt2d& theResPoint, - gp_Vec2d& theResDeriv) -{ - // 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(theBaseD1.Y(), -theBaseD1.X()); - gp_XY DNdir(theBaseD2.Y(), -theBaseD2.X()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = Sqrt (R2); - Standard_Real R3 = R * R2; - Standard_Real Dr = Ndir.Dot(DNdir); - if (R3 <= gp::Resolution()) - { - if (R2 <= gp::Resolution()) - Standard_NullValue::Raise("CSLib_Offset: Null derivative"); - //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)); - } - - // P(u) - D0(theBasePoint, theBaseD1, theOffset, Standard_False, theResPoint); - // P'(u) - theResDeriv = theBaseD1.Added(gp_Vec2d(DNdir)); -} - -void CSLib_Offset::D2(const gp_Pnt2d& theBasePoint, - const gp_Vec2d& theBaseD1, - const gp_Vec2d& theBaseD2, - const gp_Vec2d& theBaseD3, - Standard_Real theOffset, - Standard_Boolean theIsDirectionChange, - gp_Pnt2d& theResPoint, - gp_Vec2d& theResD1, - gp_Vec2d& theResD2) -{ - // 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(theBaseD1.Y(), -theBaseD1.X()); - gp_XY DNdir(theBaseD2.Y(), -theBaseD2.X()); - gp_XY D2Ndir(theBaseD3.Y(), -theBaseD3.X()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = 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()) - Standard_NullValue::Raise("CSLib_Offset: Null derivative"); - //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)); - } - - // P(u) : - D0(theBasePoint, theBaseD1, theOffset, theIsDirectionChange, theResPoint); - // P'(u) : - theResD1 = theBaseD1.Added(gp_Vec2d(DNdir)); - // P"(u) : - if (theIsDirectionChange) - theResD2 = -theBaseD2; - else - theResD2 = theBaseD2; - theResD2.Add(gp_Vec2d(D2Ndir)); -} - -void CSLib_Offset::D3(const gp_Pnt2d& theBasePoint, - const gp_Vec2d& theBaseD1, - const gp_Vec2d& theBaseD2, - const gp_Vec2d& theBaseD3, - const gp_Vec2d& theBaseD4, - Standard_Real theOffset, - Standard_Boolean theIsDirectionChange, - gp_Pnt2d& theResPoint, - gp_Vec2d& theResD1, - gp_Vec2d& theResD2, - gp_Vec2d& theResD3) -{ - // 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(theBaseD1.Y(), -theBaseD1.X()); - gp_XY DNdir(theBaseD2.Y(), -theBaseD2.X()); - gp_XY D2Ndir(theBaseD3.Y(), -theBaseD3.X()); - gp_XY D3Ndir(theBaseD4.Y(), -theBaseD4.X()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = 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()) - Standard_NullValue::Raise("CSLib_Offset: Null derivative"); - //We try another computation but the stability is not very good dixit ISG. - // V3 = P"' (U) : - D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * theOffset * Dr / R2)); - D3Ndir.Subtract ( - (DNdir.Multiplied ((3.0 * theOffset) * ((D2r/R2) + (Dr*Dr)/R4)))); - D3Ndir.Add (Ndir.Multiplied ( - (theOffset * (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)); - } - - // P(u) - D0(theBasePoint, theBaseD1, theOffset, theIsDirectionChange, theResPoint); - // P'(u) - theResD1 = theBaseD1.Added(gp_Vec2d(DNdir)); - // P"(u) - theResD2 = theBaseD2.Added(gp_Vec2d(D2Ndir)); - // P"'(u) - if (theIsDirectionChange) - theResD3 = -theBaseD3; - else - theResD3 = theBaseD3; - theResD3.Add(gp_Vec2d(D2Ndir)); -} - - -// ========== Offset values for 3D curves ========== - -void CSLib_Offset::D0(const gp_Pnt& theBasePoint, - const gp_Vec& theBaseDeriv, - const gp_Dir& theOffsetDirection, - Standard_Real theOffsetValue, - Standard_Boolean , // unused - gp_Pnt& theResPoint) -{ - gp_XYZ Ndir = (theBaseDeriv.XYZ()).Crossed(theOffsetDirection.XYZ()); - Standard_Real R = Ndir.Modulus(); - if (R <= gp::Resolution()) - Standard_NullValue::Raise("CSLib_Offset: Undefined normal vector " - "because tangent vector has zero-magnitude!"); - - Ndir.Multiply(theOffsetValue / R); - Ndir.Add(theBasePoint.XYZ()); - theResPoint.SetXYZ(Ndir); -} - -void CSLib_Offset::D1(const gp_Pnt& theBasePoint, - const gp_Vec& theBaseD1, - const gp_Vec& theBaseD2, - const gp_Dir& theOffsetDirection, - Standard_Real theOffsetValue, - Standard_Boolean , // unused - gp_Pnt& theResPoint, - gp_Vec& theResDeriv) -{ - // 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 = (theBaseD1.XYZ()).Crossed(theOffsetDirection.XYZ()); - gp_XYZ DNdir = (theBaseD2.XYZ()).Crossed(theOffsetDirection.XYZ()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = Sqrt (R2); - Standard_Real R3 = R * R2; - Standard_Real Dr = Ndir.Dot (DNdir); - if (R3 <= gp::Resolution()) { - if (R2 <= gp::Resolution()) - Standard_NullValue::Raise("CSLib_Offset: Null derivative"); - //We try another computation but the stability is not very good. - DNdir.Multiply(R); - DNdir.Subtract(Ndir.Multiplied(Dr / R)); - DNdir.Multiply(theOffsetValue / R2); - } - else { - // Same computation as IICURV in EUCLID-IS because the stability is - // better - DNdir.Multiply(theOffsetValue / R); - DNdir.Subtract(Ndir.Multiplied(theOffsetValue * Dr / R3)); - } - - // P(u) - D0(theBasePoint, theBaseD1, theOffsetDirection, theOffsetValue, Standard_False, theResPoint); - // P'(u) - theResDeriv = theBaseD1.Added(gp_Vec(DNdir)); -} - -void CSLib_Offset::D2(const gp_Pnt& theBasePoint, - const gp_Vec& theBaseD1, - const gp_Vec& theBaseD2, - const gp_Vec& theBaseD3, - const gp_Dir& theOffsetDirection, - Standard_Real theOffsetValue, - Standard_Boolean theIsDirectionChange, - gp_Pnt& theResPoint, - gp_Vec& theResD1, - gp_Vec& theResD2) -{ - // 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 = (theBaseD1.XYZ()).Crossed(theOffsetDirection.XYZ()); - gp_XYZ DNdir = (theBaseD2.XYZ()).Crossed(theOffsetDirection.XYZ()); - gp_XYZ D2Ndir = (theBaseD3.XYZ()).Crossed(theOffsetDirection.XYZ()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = 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()) - Standard_NullValue::Raise("CSLib_Offset: 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 (theOffsetValue / R); - - // V1 = P' (U) : - DNdir.Multiply(R); - DNdir.Subtract (Ndir.Multiplied (Dr/R)); - DNdir.Multiply (theOffsetValue/R2); - } - else { - // Same computation as IICURV in EUCLID-IS because the stability is - // better. - // V2 = P" (U) : - D2Ndir.Multiply (theOffsetValue/R); - D2Ndir.Subtract (DNdir.Multiplied (2.0 * theOffsetValue * Dr / R3)); - D2Ndir.Add (Ndir.Multiplied (theOffsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3)))); - - // V1 = P' (U) : - DNdir.Multiply (theOffsetValue/R); - DNdir.Subtract (Ndir.Multiplied (theOffsetValue*Dr/R3)); - } - - // P(u) : - D0(theBasePoint, theBaseD1, theOffsetDirection, theOffsetValue, theIsDirectionChange, theResPoint); - // P'(u) : - theResD1 = theBaseD1.Added(gp_Vec(DNdir)); - // P"(u) : - if (theIsDirectionChange) - theResD2 = -theBaseD2; - else - theResD2 = theBaseD2; - theResD2.Add(gp_Vec(D2Ndir)); -} - -void CSLib_Offset::D3(const gp_Pnt& theBasePoint, - const gp_Vec& theBaseD1, - const gp_Vec& theBaseD2, - const gp_Vec& theBaseD3, - const gp_Vec& theBaseD4, - const gp_Dir& theOffsetDirection, - Standard_Real theOffsetValue, - Standard_Boolean theIsDirectionChange, - gp_Pnt& theResPoint, - gp_Vec& theResD1, - gp_Vec& theResD2, - gp_Vec& theResD3) -{ - // 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 = (theBaseD1.XYZ()).Crossed(theOffsetDirection.XYZ()); - gp_XYZ DNdir = (theBaseD2.XYZ()).Crossed(theOffsetDirection.XYZ()); - gp_XYZ D2Ndir = (theBaseD3.XYZ()).Crossed(theOffsetDirection.XYZ()); - gp_XYZ D3Ndir = (theBaseD4.XYZ()).Crossed(theOffsetDirection.XYZ()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = 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()) - Standard_NullValue::Raise("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 (theOffsetValue/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 (theOffsetValue / R); - // V1 = P' (U) : - DNdir.Multiply(R); - DNdir.Subtract (Ndir.Multiplied (Dr/R)); - DNdir.Multiply (theOffsetValue/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 (theOffsetValue); - // 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 (theOffsetValue); - // V1 = P' (U) : - DNdir.Multiply (theOffsetValue/R); - DNdir.Subtract (Ndir.Multiplied (theOffsetValue*Dr/R3)); - } - - // P(u) - D0(theBasePoint, theBaseD1, theOffsetDirection, theOffsetValue, theIsDirectionChange, theResPoint); - // P'(u) - theResD1 = theBaseD1.Added(gp_Vec(DNdir)); - // P"(u) - theResD2 = theBaseD2.Added(gp_Vec(D2Ndir)); - // P"'(u) - if (theIsDirectionChange) - theResD3 = -theBaseD3; - else - theResD3 = theBaseD3; - theResD3.Add(gp_Vec(D2Ndir)); -} - diff --git a/src/CSLib/CSLib_Offset.hxx b/src/CSLib/CSLib_Offset.hxx deleted file mode 100644 index 277dfdc6d8..0000000000 --- a/src/CSLib/CSLib_Offset.hxx +++ /dev/null @@ -1,190 +0,0 @@ -// Copyright (c) 2015-... 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 _CSLib_Offset_Headerfile -#define _CSLib_Offset_Headerfile - -#include -#include -#include -#include -#include -#include - -/** \namespace CSLib_Offset - * \brief Provides a number of static methods to calculate values and derivatives - * of an offset curves and surfaces using values and derivatives of - * a base curve/surface. - */ -namespace CSLib_Offset -{ - /** \brief Calculate value of offset curve in 2D - * \param[in] theBasePoint point on a base curve - * \param[in] theBaseDeriv derivative on a base curve - * \param[in] theOffset size of offset - * \param[in] theIsDirectionChange shows that it is necessary to consider the direction of derivative (not used) - * \param[out] theResPoint point on offset curve - */ - Standard_EXPORT void D0(const gp_Pnt2d& theBasePoint, - const gp_Vec2d& theBaseDeriv, - Standard_Real theOffset, - Standard_Boolean theIsDirectionChange, - gp_Pnt2d& theResPoint); - /** \brief Calculate value of offset curve in 3D - * \param[in] theBasePoint point on a base curve - * \param[in] theBaseDeriv derivative on a base curve - * \param[in] theOffsetDirection direction of the offset - * \param[in] theOffsetValue length of the offset - * \param[in] theIsDirectionChange shows that it is necessary to consider the direction of derivative (not used) - * \param[out] theResPoint point on offset curve - */ - Standard_EXPORT void D0(const gp_Pnt& theBasePoint, - const gp_Vec& theBaseDeriv, - const gp_Dir& theOffsetDirection, - Standard_Real theOffsetValue, - Standard_Boolean theIsDirectionChange, - gp_Pnt& theResPoint); - - - /** \brief Calculate value and the first derivative of offset curve in 2D - * \param[in] theBasePoint point on a base curve - * \param[in] theBaseD1 first derivative on a base curve - * \param[in] theBaseD2 second derivative on a base curve - * \param[in] theOffset size of offset - * \param[in] theIsDirectionChange shows that it is necessary to consider the direction of derivative (not used) - * \param[out] theResPoint point on offset curve - * \param[out] theResDeriv derivative on offset curve - */ - Standard_EXPORT void D1(const gp_Pnt2d& theBasePoint, - const gp_Vec2d& theBaseD1, - const gp_Vec2d& theBaseD2, - Standard_Real theOffset, - Standard_Boolean theIsDirectionChange, - gp_Pnt2d& theResPoint, - gp_Vec2d& theResDeriv); - /** \brief Calculate value and the first derivative of offset curve in 3D - * \param[in] theBasePoint point on a base curve - * \param[in] theBaseD1 first derivative on a base curve - * \param[in] theBaseD2 second derivative on a base curve - * \param[in] theOffsetDirection direction of the offset - * \param[in] theOffsetValue length of the offset - * \param[in] theIsDirectionChange shows that it is necessary to consider the direction of derivative (not used) - * \param[out] theResPoint point on offset curve - * \param[out] theResDeriv derivative on offset curve - */ - Standard_EXPORT void D1(const gp_Pnt& theBasePoint, - const gp_Vec& theBaseD1, - const gp_Vec& theBaseD2, - const gp_Dir& theOffsetDirection, - Standard_Real theOffsetValue, - Standard_Boolean theIsDirectionChange, - gp_Pnt& theResPoint, - gp_Vec& theResDeriv); - - - /** \brief Calculate value and two derivatives of offset curve in 2D - * \param[in] theBasePoint point on a base curve - * \param[in] theBaseD1 first derivative on a base curve - * \param[in] theBaseD2 second derivative on a base curve - * \param[in] theBaseD3 third derivative on a base curve - * \param[in] theOffset size of offset - * \param[in] theIsDirectionChange shows that it is necessary to consider the direction of derivative - * \param[out] theResPoint point on offset curve - * \param[out] theResD1 first derivative on offset curve - * \param[out] theResD2 second derivative on offset curve - */ - Standard_EXPORT void D2(const gp_Pnt2d& theBasePoint, - const gp_Vec2d& theBaseD1, - const gp_Vec2d& theBaseD2, - const gp_Vec2d& theBaseD3, - Standard_Real theOffset, - Standard_Boolean theIsDirectionChange, - gp_Pnt2d& theResPoint, - gp_Vec2d& theResD1, - gp_Vec2d& theResD2); - /** \brief Calculate value and two derivatives of offset curve in 3D - * \param[in] theBasePoint point on a base curve - * \param[in] theBaseD1 first derivative on a base curve - * \param[in] theBaseD2 second derivative on a base curve - * \param[in] theBaseD3 third derivative on a base curve - * \param[in] theOffsetDirection direction of the offset - * \param[in] theOffsetValue length of the offset - * \param[in] theIsDirectionChange shows that it is necessary to consider the direction of derivative - * \param[out] theResPoint point on offset curve - * \param[out] theResD1 first derivative on offset curve - * \param[out] theResD2 second derivative on offset curve - */ - Standard_EXPORT void D2(const gp_Pnt& theBasePoint, - const gp_Vec& theBaseD1, - const gp_Vec& theBaseD2, - const gp_Vec& theBaseD3, - const gp_Dir& theOffsetDirection, - Standard_Real theOffsetValue, - Standard_Boolean theIsDirectionChange, - gp_Pnt& theResPoint, - gp_Vec& theResD1, - gp_Vec& theResD2); - - /** \brief Calculate value and three derivatives of offset curve in 2D - * \param[in] theBasePoint point on a base curve - * \param[in] theBaseD1 first derivative on a base curve - * \param[in] theBaseD2 second derivative on a base curve - * \param[in] theBaseD3 third derivative on a base curve - * \param[in] theBaseD4 fourth derivative on a base curve - * \param[in] theOffset size of offset - * \param[in] theIsDirectionChange shows that it is necessary to consider the direction of derivative - * \param[out] theResPoint point on offset curve - * \param[out] theResD1 first derivative on offset curve - * \param[out] theResD2 second derivative on offset curve - * \param[out] theResD3 third derivative on offset curve - */ - Standard_EXPORT void D3(const gp_Pnt2d& theBasePoint, - const gp_Vec2d& theBaseD1, - const gp_Vec2d& theBaseD2, - const gp_Vec2d& theBaseD3, - const gp_Vec2d& theBaseD4, - Standard_Real theOffset, - Standard_Boolean theIsDirectionChange, - gp_Pnt2d& theResPoint, - gp_Vec2d& theResD1, - gp_Vec2d& theResD2, - gp_Vec2d& theResD3); - /** \brief Calculate value and three derivatives of offset curve in 3D - * \param[in] theBasePoint point on a base curve - * \param[in] theBaseD1 first derivative on a base curve - * \param[in] theBaseD2 second derivative on a base curve - * \param[in] theBaseD3 third derivative on a base curve - * \param[in] theBaseD4 fourth derivative on a base curve - * \param[in] theOffsetDirection direction of the offset - * \param[in] theOffsetValue length of the offset - * \param[in] theIsDirectionChange shows that it is necessary to consider the direction of derivative - * \param[out] theResPoint point on offset curve - * \param[out] theResD1 first derivative on offset curve - * \param[out] theResD2 second derivative on offset curve - * \param[out] theResD3 third derivative on offset curve - */ - Standard_EXPORT void D3(const gp_Pnt& theBasePoint, - const gp_Vec& theBaseD1, - const gp_Vec& theBaseD2, - const gp_Vec& theBaseD3, - const gp_Vec& theBaseD4, - const gp_Dir& theOffsetDirection, - Standard_Real theOffsetValue, - Standard_Boolean theIsDirectionChange, - gp_Pnt& theResPoint, - gp_Vec& theResD1, - gp_Vec& theResD2, - gp_Vec& theResD3); -} - -#endif // _CSLib_Offset_Headerfile diff --git a/src/CSLib/FILES b/src/CSLib/FILES index b96029feda..65d03a194f 100644 --- a/src/CSLib/FILES +++ b/src/CSLib/FILES @@ -6,5 +6,3 @@ CSLib_DerivativeStatus.hxx CSLib_NormalPolyDef.cxx CSLib_NormalPolyDef.hxx CSLib_NormalStatus.hxx -CSLib_Offset.cxx -CSLib_Offset.hxx diff --git a/src/Geom/Geom_OffsetCurve.cxx b/src/Geom/Geom_OffsetCurve.cxx index 86743bcc92..1b1e5071d7 100644 --- a/src/Geom/Geom_OffsetCurve.cxx +++ b/src/Geom/Geom_OffsetCurve.cxx @@ -18,7 +18,6 @@ // avec discernement ! // 19-09-97 : JPI correction derivee seconde -#include #include #include #include @@ -44,28 +43,8 @@ #include #include -typedef Geom_OffsetCurve OffsetCurve; -typedef Geom_Curve Curve; -typedef gp_Dir Dir; -typedef gp_Pnt Pnt; -typedef gp_Trsf Trsf; -typedef gp_Vec Vec; -typedef gp_XYZ XYZ; - -//ordre de derivation maximum pour la recherche de la premiere -//derivee non nulle -static const int maxDerivOrder = 3; -static const Standard_Real MinStep = 1e-7; -static const Standard_Real MyAngularToleranceForG1 = Precision::Angular(); - - -static gp_Vec dummyDerivative; // used as empty value for unused derivatives in AdjustDerivative -// Recalculate derivatives in the singular point -// Returns true if the direction of derivatives is changed -static Standard_Boolean AdjustDerivative( - const Handle(Geom_Curve)& theCurve, Standard_Integer theMaxDerivative, Standard_Real theU, gp_Vec& theD1, - gp_Vec& theD2 = dummyDerivative, gp_Vec& theD3 = dummyDerivative, gp_Vec& theD4 = dummyDerivative); +static const Standard_Real MyAngularToleranceForG1 = Precision::Angular(); //======================================================================= @@ -76,7 +55,7 @@ static Standard_Boolean AdjustDerivative( Handle(Geom_Geometry) Geom_OffsetCurve::Copy () const { Handle(Geom_OffsetCurve) C; - C = new OffsetCurve (basisCurve, offsetValue, direction); + C = new Geom_OffsetCurve (basisCurve, offsetValue, direction); return C; } @@ -107,6 +86,7 @@ void Geom_OffsetCurve::Reverse () { basisCurve->Reverse(); offsetValue = -offsetValue; + myEvaluator->SetOffsetValue(offsetValue); } @@ -133,8 +113,11 @@ const gp_Dir& Geom_OffsetCurve::Direction () const //purpose : //======================================================================= -void Geom_OffsetCurve::SetDirection (const Dir& V) - { direction = V; } +void Geom_OffsetCurve::SetDirection (const gp_Dir& V) +{ + direction = V; + myEvaluator->SetOffsetDirection(direction); +} //======================================================================= //function : SetOffsetValue @@ -142,7 +125,10 @@ void Geom_OffsetCurve::SetDirection (const Dir& V) //======================================================================= void Geom_OffsetCurve::SetOffsetValue (const Standard_Real D) - { offsetValue = D; } +{ + offsetValue = D; + myEvaluator->SetOffsetValue(offsetValue); +} //======================================================================= @@ -242,6 +228,8 @@ void Geom_OffsetCurve::SetBasisCurve (const Handle(Geom_Curve)& C, { basisCurve = aCheckingCurve; } + + myEvaluator = new GeomEvaluator_OffsetCurve(basisCurve, offsetValue, direction); } @@ -283,70 +271,40 @@ GeomAbs_Shape Geom_OffsetCurve::Continuity () const { //purpose : //======================================================================= -void Geom_OffsetCurve::D0 (const Standard_Real U, Pnt& P) const +void Geom_OffsetCurve::D0 (const Standard_Real U, gp_Pnt& P) const { - gp_Pnt PBasis; - gp_Vec VBasis; - D0(U,P,PBasis,VBasis); + myEvaluator->D0(U, P); } - //======================================================================= //function : D1 //purpose : //======================================================================= -void Geom_OffsetCurve::D1 (const Standard_Real U, Pnt& P, Vec& V1) const +void Geom_OffsetCurve::D1 (const Standard_Real U, gp_Pnt& P, gp_Vec& V1) const { - gp_Pnt PBasis; - gp_Vec V1Basis,V2Basis; - D1(U,P,PBasis,V1,V1Basis,V2Basis); + myEvaluator->D1(U, P, V1); } - - //======================================================================= //function : D2 //purpose : //======================================================================= -void Geom_OffsetCurve::D2 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2) const +void Geom_OffsetCurve::D2 (const Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2) const { - gp_Pnt PBasis; - gp_Vec V1Basis,V2Basis,V3Basis; - D2(U,P,PBasis,V1,V2,V1Basis,V2Basis,V3Basis); + myEvaluator->D2(U, P, V1, V2); } - //======================================================================= //function : D3 //purpose : //======================================================================= -void Geom_OffsetCurve::D3 (const Standard_Real theU, Pnt& theP, Vec& theV1, Vec& theV2, Vec& theV3) const +void Geom_OffsetCurve::D3 (const Standard_Real theU, gp_Pnt& theP, + gp_Vec& theV1, gp_Vec& theV2, gp_Vec& theV3) 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 - - Standard_Boolean IsDirectionChange = Standard_False; - - basisCurve->D3 (theU, theP, theV1, theV2, theV3); - Vec aV4 = basisCurve->DN (theU, 4); - if(theV1.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(basisCurve, 4, theU, theV1, theV2, theV3, aV4); - - CSLib_Offset::D3(theP, theV1, theV2, theV3, aV4, direction, offsetValue, - IsDirectionChange, theP, theV1, theV2, theV3); + myEvaluator->D3(theU, theP, theV1, theV2, theV3); } @@ -355,8 +313,8 @@ void Geom_OffsetCurve::D3 (const Standard_Real theU, Pnt& theP, Vec& theV1, Vec& //purpose : //======================================================================= -Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const - { +gp_Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const +{ Standard_RangeError_Raise_if (N < 1, "Exception: " "Geom_OffsetCurve::DN(...). N<1."); @@ -381,83 +339,14 @@ Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const return VN; } -//======================================================================= -//function : D0 -//purpose : -//======================================================================= - -void Geom_OffsetCurve::D0(const Standard_Real theU, gp_Pnt& theP, - gp_Pnt& thePbasis, gp_Vec& theV1basis)const -{ - basisCurve->D1(theU, thePbasis, theV1basis); - Standard_Boolean IsDirectionChange = Standard_False; - if(theV1basis.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(basisCurve, 1, theU, theV1basis); - - CSLib_Offset::D0(thePbasis, theV1basis, direction, offsetValue, IsDirectionChange, theP); -} - -//======================================================================= -//function : D1 -//purpose : -//======================================================================= - -void Geom_OffsetCurve::D1 ( const Standard_Real theU, - Pnt& theP , Pnt& thePBasis , - Vec& theV1, Vec& theV1basis, Vec& theV2basis) 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)) - - basisCurve->D2 (theU, thePBasis, theV1basis, theV2basis); - - Standard_Boolean IsDirectionChange = Standard_False; - if(theV1basis.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(basisCurve, 2, theU, theV1basis, theV2basis); - - CSLib_Offset::D1(thePBasis, theV1basis, theV2basis, direction, offsetValue, IsDirectionChange, theP, theV1); -} - - -//======================================================================= -//function : D2 -//purpose : -//======================================================================= - -void Geom_OffsetCurve::D2 (const Standard_Real theU, - Pnt& theP, Pnt& thePBasis, - Vec& theV1, Vec& theV2, - Vec& theV1basis, Vec& theV2basis, Vec& theV3basis) 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))) - - Standard_Boolean IsDirectionChange = Standard_False; - - basisCurve->D3 (theU, thePBasis, theV1basis, theV2basis, theV3basis); - - if(theV1basis.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(basisCurve, 3, theU, theV1basis, theV2basis, theV3basis); - - CSLib_Offset::D2(thePBasis, theV1basis, theV2basis, theV3basis, direction, offsetValue, - IsDirectionChange, theP, theV1, theV2); -} - //======================================================================= //function : FirstParameter //purpose : //======================================================================= -Standard_Real Geom_OffsetCurve::FirstParameter () const { - +Standard_Real Geom_OffsetCurve::FirstParameter () const +{ return basisCurve->FirstParameter(); } @@ -467,8 +356,8 @@ Standard_Real Geom_OffsetCurve::FirstParameter () const { //purpose : //======================================================================= -Standard_Real Geom_OffsetCurve::LastParameter () const { - +Standard_Real Geom_OffsetCurve::LastParameter () const +{ return basisCurve->LastParameter(); } @@ -478,22 +367,8 @@ Standard_Real Geom_OffsetCurve::LastParameter () const { //purpose : //======================================================================= -Standard_Real Geom_OffsetCurve::Offset () const { return offsetValue; } - -//======================================================================= -//function : Value -//purpose : -//======================================================================= - -void Geom_OffsetCurve::Value (const Standard_Real theU, Pnt& theP, - Pnt& thePbasis, Vec& theV1basis) const -{ - if (myBasisCurveContinuity == GeomAbs_C0) - Geom_UndefinedValue::Raise("Exception: Basis curve is C0 continuity!"); - - basisCurve->D1(theU, thePbasis, theV1basis); - D0(theU,theP); -} +Standard_Real Geom_OffsetCurve::Offset () const +{ return offsetValue; } //======================================================================= @@ -528,11 +403,14 @@ Standard_Boolean Geom_OffsetCurve::IsCN (const Standard_Integer N) const { //purpose : //======================================================================= -void Geom_OffsetCurve::Transform (const Trsf& T) { - +void Geom_OffsetCurve::Transform (const gp_Trsf& T) +{ basisCurve->Transform (T); direction.Transform(T); offsetValue *= T.ScaleFactor(); + + myEvaluator->SetOffsetValue(offsetValue); + myEvaluator->SetOffsetDirection(direction); } //======================================================================= @@ -565,57 +443,3 @@ GeomAbs_Shape Geom_OffsetCurve::GetBasisCurveContinuity() const { return myBasisCurveContinuity; } - - -// ============= Auxiliary functions =================== -Standard_Boolean AdjustDerivative(const Handle(Geom_Curve)& theCurve, Standard_Integer theMaxDerivative, - Standard_Real theU, gp_Vec& theD1, gp_Vec& theD2, - gp_Vec& theD3, gp_Vec& theD4) -{ - static const Standard_Real aTol = gp::Resolution(); - - Standard_Boolean IsDirectionChange = Standard_False; - const Standard_Real anUinfium = theCurve->FirstParameter(); - const Standard_Real anUsupremum = theCurve->LastParameter(); - - 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 = Max(du * DivisionFactor, MinStep); - - //Derivative is approximated by Taylor-series - Standard_Integer anIndex = 1; //Derivative order - gp_Vec V; - - do - { - V = theCurve->DN(theU, ++anIndex); - } - while((V.SquareMagnitude() <= aTol) && anIndex < maxDerivOrder); - - Standard_Real u; - - if(theU-anUinfium < aDelta) - u = theU+aDelta; - else - u = theU-aDelta; - - gp_Pnt P1, P2; - theCurve->D0(Min(theU, u), P1); - theCurve->D0(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]) = theCurve->DN(theU, anIndex + i) * aSign; - - return IsDirectionChange; -} diff --git a/src/Geom/Geom_OffsetCurve.hxx b/src/Geom/Geom_OffsetCurve.hxx index 790ca895f7..38b0fec7cf 100644 --- a/src/Geom/Geom_OffsetCurve.hxx +++ b/src/Geom/Geom_OffsetCurve.hxx @@ -26,6 +26,8 @@ #include #include #include +#include + class Geom_Curve; class Standard_ConstructionError; class Standard_RangeError; @@ -209,35 +211,16 @@ public: //! direction. //! Raised if N < 1. Standard_EXPORT gp_Vec DN (const Standard_Real U, const Standard_Integer N) const Standard_OVERRIDE; - - //! Warning! this should not be called - //! if the basis curve is not at least C1. Nevertheless - //! if used on portion where the curve is C1, it is OK - Standard_EXPORT void Value (const Standard_Real U, gp_Pnt& P, gp_Pnt& Pbasis, gp_Vec& V1basis) const; - - //! Warning! this should not be called - //! if the continuity of the basis curve is not C1. - //! Nevertheless, it's OK to use it on portion - //! where the curve is C1 - Standard_EXPORT void D0 (const Standard_Real U, gp_Pnt& P, gp_Pnt& Pbasis, gp_Vec& V1basis) const; - - //! Warning! this should not be called - //! if the continuity of the basis curve is not C1. - //! Nevertheless, it's OK to use it on portion - //! where the curve is C1 - Standard_EXPORT void D1 (const Standard_Real U, gp_Pnt& P, gp_Pnt& Pbasis, gp_Vec& V1, gp_Vec& V1basis, gp_Vec& V2basis) const; - - //! Warning! this should not be called - //! if the continuity of the basis curve is not C3. - //! Nevertheless, it's OK to use it on portion - //! where the curve is C3 - Standard_EXPORT void D2 (const Standard_Real U, gp_Pnt& P, gp_Pnt& Pbasis, gp_Vec& V1, gp_Vec& V2, gp_Vec& V1basis, gp_Vec& V2basis, gp_Vec& V3basis) const; - - Standard_EXPORT Standard_Real FirstParameter() const Standard_OVERRIDE; - - //! Returns the value of the first or last parameter of this + + //! Returns the value of the first parameter of this //! offset curve. The first parameter corresponds to the - //! start point of the curve. The last parameter + //! start point of the curve. + //! Note: the first and last parameters of this offset curve + //! are also the ones of its basis curve. + Standard_EXPORT Standard_Real FirstParameter() const Standard_OVERRIDE; + + //! Returns the value of the last parameter of this + //! offset curve. The last parameter //! corresponds to the end point. //! Note: the first and last parameters of this offset curve //! are also the ones of its basis curve. @@ -315,7 +298,7 @@ private: gp_Dir direction; Standard_Real offsetValue; GeomAbs_Shape myBasisCurveContinuity; - + Handle(GeomEvaluator_OffsetCurve) myEvaluator; }; diff --git a/src/Geom2d/Geom2d_OffsetCurve.cxx b/src/Geom2d/Geom2d_OffsetCurve.cxx index c05249c329..fb47b715ca 100644 --- a/src/Geom2d/Geom2d_OffsetCurve.cxx +++ b/src/Geom2d/Geom2d_OffsetCurve.cxx @@ -16,7 +16,6 @@ // modified by Edward AGAPOV (eap) Jan 28 2002 --- DN(), occ143(BUC60654) -#include #include #include #include @@ -42,28 +41,9 @@ #include #include -typedef Geom2d_OffsetCurve OffsetCurve; -typedef Geom2d_Curve Curve; -typedef gp_Dir2d Dir2d; -typedef gp_Pnt2d Pnt2d; -typedef gp_Vec2d Vec2d; -typedef gp_Trsf2d Trsf2d; -typedef gp_XY XY; - -//ordre de derivation maximum pour la recherche de la premiere -//derivee non nulle -static const int maxDerivOrder = 3; -static const Standard_Real MinStep = 1e-7; static const Standard_Real MyAngularToleranceForG1 = Precision::Angular(); -static gp_Vec2d dummyDerivative; // used as empty value for unused derivatives in AdjustDerivative - -// Recalculate derivatives in the singular point -// Returns true if the direction of derivatives is changed -static Standard_Boolean AdjustDerivative(const Handle(Geom2d_Curve)& theCurve, Standard_Integer theMaxDerivative, - Standard_Real theU, gp_Vec2d& theD1, gp_Vec2d& theD2 = dummyDerivative, - gp_Vec2d& theD3 = dummyDerivative, gp_Vec2d& theD4 = dummyDerivative); //======================================================================= //function : Copy @@ -73,7 +53,7 @@ static Standard_Boolean AdjustDerivative(const Handle(Geom2d_Curve)& theCurve, S Handle(Geom2d_Geometry) Geom2d_OffsetCurve::Copy () const { Handle(Geom2d_OffsetCurve) C; - C = new OffsetCurve (basisCurve, offsetValue); + C = new Geom2d_OffsetCurve (basisCurve, offsetValue); return C; } @@ -177,6 +157,7 @@ void Geom2d_OffsetCurve::SetBasisCurve (const Handle(Geom2d_Curve)& C, basisCurve = aCheckingCurve; } + myEvaluator = new Geom2dEvaluator_OffsetCurve(basisCurve, offsetValue); } //======================================================================= @@ -184,7 +165,11 @@ void Geom2d_OffsetCurve::SetBasisCurve (const Handle(Geom2d_Curve)& C, //purpose : //======================================================================= -void Geom2d_OffsetCurve::SetOffsetValue (const Standard_Real D) { offsetValue = D; } +void Geom2d_OffsetCurve::SetOffsetValue (const Standard_Real D) +{ + offsetValue = D; + myEvaluator->SetOffsetValue(offsetValue); +} //======================================================================= //function : BasisCurve @@ -223,37 +208,18 @@ GeomAbs_Shape Geom2d_OffsetCurve::Continuity () const //======================================================================= void Geom2d_OffsetCurve::D0 (const Standard_Real theU, - Pnt2d& theP) const - { - Vec2d vD1; - basisCurve->D1 (theU, theP, vD1); - - Standard_Boolean IsDirectionChange = Standard_False; - if(vD1.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(basisCurve, 1, theU, vD1); - - CSLib_Offset::D0(theP, vD1, offsetValue, IsDirectionChange, theP); + gp_Pnt2d& theP) const +{ + myEvaluator->D0(theU, theP); } //======================================================================= //function : D1 //purpose : //======================================================================= -void Geom2d_OffsetCurve::D1 (const Standard_Real theU, Pnt2d& theP, Vec2d& theV1) const +void Geom2d_OffsetCurve::D1 (const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d& theV1) const { - // 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)) - - Vec2d V2; - basisCurve->D2 (theU, theP, theV1, V2); - - Standard_Boolean IsDirectionChange = Standard_False; - if(theV1.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(basisCurve, 2, theU, theV1, V2); - - CSLib_Offset::D1(theP, theV1, V2, offsetValue, IsDirectionChange, theP, theV1); + myEvaluator->D1(theU, theP, theV1); } //======================================================================= @@ -262,25 +228,10 @@ void Geom2d_OffsetCurve::D1 (const Standard_Real theU, Pnt2d& theP, Vec2d& theV1 //======================================================================= void Geom2d_OffsetCurve::D2 (const Standard_Real theU, - Pnt2d& theP, - Vec2d& theV1, Vec2d& theV2) const + gp_Pnt2d& theP, + gp_Vec2d& theV1, gp_Vec2d& theV2) const { - // 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))) - - Vec2d V3; - basisCurve->D3 (theU, theP, theV1, theV2, V3); - - Standard_Boolean IsDirectionChange = Standard_False; - if(theV1.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(basisCurve, 3, theU, theV1, theV2, V3); - - CSLib_Offset::D2(theP, theV1, theV2, V3, offsetValue, IsDirectionChange, theP, theV1, theV2); + myEvaluator->D2(theU, theP, theV1, theV2); } @@ -290,32 +241,10 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real theU, //======================================================================= void Geom2d_OffsetCurve::D3 (const Standard_Real theU, - Pnt2d& theP, - Vec2d& theV1, Vec2d& theV2, Vec2d& theV3) const + gp_Pnt2d& theP, + gp_Vec2d& theV1, gp_Vec2d& theV2, gp_Vec2d& theV3) const { - - // 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 - - basisCurve->D3 (theU, theP, theV1, theV2, theV3); - Vec2d V4 = basisCurve->DN (theU, 4); - - Standard_Boolean IsDirectionChange = Standard_False; - if(theV1.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(basisCurve, 4, theU, theV1, theV2, theV3, V4); - - CSLib_Offset::D3(theP, theV1, theV2, theV3, V4, offsetValue, IsDirectionChange, - theP, theV1, theV2, theV3); + myEvaluator->D3(theU, theP, theV1, theV2, theV3); } //======================================================================= @@ -323,10 +252,10 @@ void Geom2d_OffsetCurve::D3 (const Standard_Real theU, //purpose : //======================================================================= -Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U, - const Standard_Integer N) const +gp_Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U, + const Standard_Integer N) const { -Standard_RangeError_Raise_if (N < 1, "Exception: Geom2d_OffsetCurve::DN(). N<1."); + Standard_RangeError_Raise_if (N < 1, "Exception: Geom2d_OffsetCurve::DN(). N<1."); gp_Vec2d VN, VBidon; gp_Pnt2d PBidon; @@ -342,88 +271,6 @@ Standard_RangeError_Raise_if (N < 1, "Exception: Geom2d_OffsetCurve::DN(). N<1." return VN; } - -//======================================================================= -//function : Value -//purpose : -//======================================================================= - -void Geom2d_OffsetCurve::Value (const Standard_Real theU, - Pnt2d& theP, Pnt2d& thePbasis, - Vec2d& theV1basis ) const -{ - basisCurve->D1(theU, thePbasis, theV1basis); - D0(theU,theP); -} - - -//======================================================================= -//function : D1 -//purpose : -//======================================================================= - -void Geom2d_OffsetCurve::D1 (const Standard_Real U, - Pnt2d& P, Pnt2d& Pbasis, - Vec2d& V1, Vec2d& V1basis, - Vec2d& V2basis ) const -{ - // 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)) - - basisCurve->D2 (U, Pbasis, V1basis, V2basis); - V1 = V1basis; - Vec2d V2 = V2basis; - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= maxDerivOrder) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) { - V2 = basisCurve->DN (U, Index); - } - - CSLib_Offset::D1(P, V1, V2, offsetValue, Standard_False, P, V1); -} - - -//======================================================================= -//function : D2 -//purpose : -//======================================================================= - -void Geom2d_OffsetCurve::D2 (const Standard_Real U, - Pnt2d& P, Pnt2d& Pbasis, - Vec2d& V1, Vec2d& V2, - Vec2d& V1basis, Vec2d& V2basis, - Vec2d& V3basis ) const -{ - // 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))) - - basisCurve->D3 (U, Pbasis, V1basis, V2basis, V3basis); - Standard_Integer Index = 2; - V1 = V1basis; - V2 = V2basis; - Vec2d V3 = V3basis; - while (V1.Magnitude() <= gp::Resolution() && Index <= maxDerivOrder) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) { - V2 = basisCurve->DN (U, Index); - V3 = basisCurve->DN (U, Index + 1); - } - - CSLib_Offset::D2(P, V1, V2, V3, offsetValue, Standard_False, P, V1, V2); -} - //======================================================================= //function : FirstParameter //purpose : @@ -450,7 +297,8 @@ Standard_Real Geom2d_OffsetCurve::LastParameter () const //purpose : //======================================================================= -Standard_Real Geom2d_OffsetCurve::Offset () const { return offsetValue; } +Standard_Real Geom2d_OffsetCurve::Offset () const +{ return offsetValue; } //======================================================================= //function : IsClosed @@ -501,10 +349,12 @@ Standard_Real Geom2d_OffsetCurve::Period() const //purpose : //======================================================================= -void Geom2d_OffsetCurve::Transform (const Trsf2d& T) +void Geom2d_OffsetCurve::Transform (const gp_Trsf2d& T) { basisCurve->Transform (T); offsetValue *= Abs(T.ScaleFactor()); + + myEvaluator->SetOffsetValue(offsetValue); } @@ -537,57 +387,3 @@ GeomAbs_Shape Geom2d_OffsetCurve::GetBasisCurveContinuity() const { return myBasisCurveContinuity; } - - -// ============= Auxiliary functions =================== -Standard_Boolean AdjustDerivative(const Handle(Geom2d_Curve)& theCurve, Standard_Integer theMaxDerivative, - Standard_Real theU, gp_Vec2d& theD1, gp_Vec2d& theD2, - gp_Vec2d& theD3, gp_Vec2d& theD4) -{ - static const Standard_Real aTol = gp::Resolution(); - - Standard_Boolean IsDirectionChange = Standard_False; - const Standard_Real anUinfium = theCurve->FirstParameter(); - const Standard_Real anUsupremum = theCurve->LastParameter(); - - 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 = Max(du * DivisionFactor, MinStep); - - //Derivative is approximated by Taylor-series - Standard_Integer anIndex = 1; //Derivative order - Vec2d V; - - do - { - V = theCurve->DN(theU, ++anIndex); - } - while((V.SquareMagnitude() <= aTol) && anIndex < maxDerivOrder); - - Standard_Real u; - - if(theU-anUinfium < aDelta) - u = theU+aDelta; - else - u = theU-aDelta; - - Pnt2d P1, P2; - theCurve->D0(Min(theU, u),P1); - theCurve->D0(Max(theU, u),P2); - - 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]) = theCurve->DN(theU, anIndex + i) * aSign; - - return IsDirectionChange; -} diff --git a/src/Geom2d/Geom2d_OffsetCurve.hxx b/src/Geom2d/Geom2d_OffsetCurve.hxx index d35ab3dede..8cf9f6c9e4 100644 --- a/src/Geom2d/Geom2d_OffsetCurve.hxx +++ b/src/Geom2d/Geom2d_OffsetCurve.hxx @@ -25,6 +25,8 @@ #include #include #include +#include + class Geom2d_Curve; class Standard_ConstructionError; class Standard_RangeError; @@ -204,28 +206,15 @@ public: //! raised if it is not possible to compute a unique offset direction. Standard_EXPORT gp_Vec2d DN (const Standard_Real U, const Standard_Integer N) const Standard_OVERRIDE; - //! Warning! this should not be called - //! if the basis curve is not at least C1. Nevertheless - //! if used on portion where the curve is C1, it is OK - Standard_EXPORT void Value (const Standard_Real U, gp_Pnt2d& P, gp_Pnt2d& Pbasis, gp_Vec2d& V1basis) const; - - //! Warning! this should not be called - //! if the continuity of the basis curve is not C1. - //! Nevertheless, it's OK to use it on portion - //! where the curve is C1 - Standard_EXPORT void D1 (const Standard_Real U, gp_Pnt2d& P, gp_Pnt2d& Pbasis, gp_Vec2d& V1, gp_Vec2d& V1basis, gp_Vec2d& V2basis) const; - - //! Warning! this should not be called - //! if the continuity of the basis curve is not C3. - //! Nevertheless, it's OK to use it on portion - //! where the curve is C3 - Standard_EXPORT void D2 (const Standard_Real U, gp_Pnt2d& P, gp_Pnt2d& Pbasis, gp_Vec2d& V1, gp_Vec2d& V2, gp_Vec2d& V1basis, gp_Vec2d& V2basis, gp_Vec2d& V3basis) const; - - Standard_EXPORT Standard_Real FirstParameter() const Standard_OVERRIDE; - - //! Returns the value of the first or last parameter of this + //! Returns the value of the first parameter of this //! offset curve. The first parameter corresponds to the - //! start point of the curve. The last parameter + //! start point of the curve. + //! Note: the first and last parameters of this offset curve + //! are also the ones of its basis curve. + Standard_EXPORT Standard_Real FirstParameter() const Standard_OVERRIDE; + + //! Returns the value of the last parameter of this + //! offset curve. The last parameter //! corresponds to the end point. //! Note: the first and last parameters of this offset curve //! are also the ones of its basis curve. @@ -312,7 +301,7 @@ private: Handle(Geom2d_Curve) basisCurve; Standard_Real offsetValue; GeomAbs_Shape myBasisCurveContinuity; - + Handle(Geom2dEvaluator_OffsetCurve) myEvaluator; }; diff --git a/src/Geom2dAdaptor/Geom2dAdaptor_Curve.cxx b/src/Geom2dAdaptor/Geom2dAdaptor_Curve.cxx index 31ad4c1653..432b317a2b 100644 --- a/src/Geom2dAdaptor/Geom2dAdaptor_Curve.cxx +++ b/src/Geom2dAdaptor/Geom2dAdaptor_Curve.cxx @@ -25,7 +25,6 @@ #include #include #include -#include #include #include #include @@ -40,6 +39,7 @@ #include #include #include +#include #include #include #include @@ -64,16 +64,6 @@ //#include static const Standard_Real PosTol = Precision::PConfusion() / 2; -static const Standard_Integer maxDerivOrder = 3; -static const Standard_Real MinStep = 1e-7; - -static gp_Vec2d dummyDerivative; // used as empty value for unused derivatives in AdjustDerivative - -// Recalculate derivatives in the singular point -// Returns true is the direction of derivatives is changed -static Standard_Boolean AdjustDerivative(const Handle(Adaptor2d_HCurve2d)& theAdaptor, Standard_Integer theMaxDerivative, - Standard_Real theU, gp_Vec2d& theD1, gp_Vec2d& theD2 = dummyDerivative, - gp_Vec2d& theD3 = dummyDerivative, gp_Vec2d& theD4 = dummyDerivative); //======================================================================= //function : LocalContinuity @@ -196,6 +186,7 @@ void Geom2dAdaptor_Curve::load(const Handle(Geom2d_Curve)& C, if ( myCurve != C) { myCurve = C; myCurveCache = Handle(BSplCLib_Cache)(); + myNestedEvaluator = Handle(Geom2dEvaluator_Curve)(); Handle(Standard_Type) TheType = C->DynamicType(); if ( TheType == STANDARD_TYPE(Geom2d_TrimmedCurve)) { @@ -236,9 +227,11 @@ void Geom2dAdaptor_Curve::load(const Handle(Geom2d_Curve)& C, else if ( TheType == STANDARD_TYPE(Geom2d_OffsetCurve)) { myTypeCurve = GeomAbs_OffsetCurve; + Handle(Geom2d_OffsetCurve) anOffsetCurve = Handle(Geom2d_OffsetCurve)::DownCast(myCurve); // Create nested adaptor for base curve - Handle(Geom2d_Curve) aBase = Handle(Geom2d_OffsetCurve)::DownCast(myCurve)->BasisCurve(); - myOffsetBaseCurveAdaptor = new Geom2dAdaptor_HCurve(aBase); + Handle(Geom2d_Curve) aBaseCurve = anOffsetCurve->BasisCurve(); + Handle(Geom2dAdaptor_HCurve) aBaseAdaptor = new Geom2dAdaptor_HCurve(aBaseCurve); + myNestedEvaluator = new Geom2dEvaluator_OffsetCurve(aBaseAdaptor, anOffsetCurve->Offset()); } else { myTypeCurve = GeomAbs_OtherCurve; @@ -385,7 +378,8 @@ Standard_Integer Geom2dAdaptor_Curve::NbIntervals(const GeomAbs_Shape S) const case GeomAbs_C2: BaseS = GeomAbs_C3; break; default: BaseS = GeomAbs_CN; } - myNbIntervals = myOffsetBaseCurveAdaptor->NbIntervals(BaseS); + Geom2dAdaptor_Curve anAdaptor( Handle(Geom2d_OffsetCurve)::DownCast(myCurve)->BasisCurve() ); + myNbIntervals = anAdaptor.NbIntervals(BaseS); } return myNbIntervals; @@ -501,8 +495,10 @@ void Geom2dAdaptor_Curve::Intervals(TColStd_Array1OfReal& T, case GeomAbs_C2: BaseS = GeomAbs_C3; break; default: BaseS = GeomAbs_CN; } - myNbIntervals = myOffsetBaseCurveAdaptor->NbIntervals(BaseS); - myOffsetBaseCurveAdaptor->Intervals(T, BaseS); + + Geom2dAdaptor_Curve anAdaptor( Handle(Geom2d_OffsetCurve)::DownCast(myCurve)->BasisCurve() ); + myNbIntervals = anAdaptor.NbIntervals(BaseS); + anAdaptor.Intervals(T, BaseS); } T( T.Lower() ) = myFirst; @@ -589,78 +585,47 @@ void Geom2dAdaptor_Curve::RebuildCache(const Standard_Real theParameter) const } //======================================================================= -//function : Value +//function : IsBoundary //purpose : //======================================================================= - -gp_Pnt2d Geom2dAdaptor_Curve::Value(const Standard_Real U) const +Standard_Boolean Geom2dAdaptor_Curve::IsBoundary(const Standard_Real theU, + Standard_Integer& theSpanStart, + Standard_Integer& theSpanFinish) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - return ValueBSpline(U); - else if (myTypeCurve == GeomAbs_OffsetCurve) - return ValueOffset(U); - else if (myTypeCurve == GeomAbs_BezierCurve) - { // use cached data - gp_Pnt2d aRes; - myCurveCache->D0(U, aRes); - return aRes; - } - - return myCurve->Value(U); -} - -//======================================================================= -//function : ValueBSpline -//purpose : Computes the point of parameter U on the B-spline curve -//======================================================================= -gp_Pnt2d Geom2dAdaptor_Curve::ValueBSpline(const Standard_Real theU) const -{ - if (theU == myFirst || theU == myLast) + Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); + if (!aBspl.IsNull() && (theU == myFirst || theU == myLast)) { - Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + aBspl->LocateU(myFirst, PosTol, theSpanStart, theSpanFinish); + if (theSpanStart < 1) + theSpanStart = 1; + if (theSpanStart >= theSpanFinish) + theSpanFinish = theSpanStart + 1; } - if (theU == myLast) + else if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + aBspl->LocateU(myLast, PosTol, theSpanStart, theSpanFinish); + if (theSpanFinish > aBspl->NbKnots()) + theSpanFinish = aBspl->NbKnots(); + if (theSpanStart >= theSpanFinish) + theSpanStart = theSpanFinish - 1; } - return aBspl->LocalValue(theU, Ideb, Ifin); - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - gp_Pnt2d aRes; - myCurveCache->D0(theU, aRes); - return aRes; + return Standard_True; } - return myCurve->Value(theU); + return Standard_False; } //======================================================================= -//function : ValueOffset -//purpose : Computes the point of parameter U on the offset curve +//function : Value +//purpose : //======================================================================= -gp_Pnt2d Geom2dAdaptor_Curve::ValueOffset(const Standard_Real theU) const + +gp_Pnt2d Geom2dAdaptor_Curve::Value(const Standard_Real U) const { - gp_Pnt2d aP; - gp_Vec2d aD1; - myOffsetBaseCurveAdaptor->D1(theU, aP, aD1); - Standard_Boolean isDirectionChange = Standard_False; - const Standard_Real aTol = gp::Resolution(); - if(aD1.SquareMagnitude() <= aTol) - isDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 1, theU, aD1); - - Standard_Real anOffset = Handle(Geom2d_OffsetCurve)::DownCast(myCurve)->Offset(); - CSLib_Offset::D0(aP, aD1, anOffset, isDirectionChange, aP); - return aP; + gp_Pnt2d aRes; + D0(U, aRes); + return aRes; } //======================================================================= @@ -670,56 +635,35 @@ gp_Pnt2d Geom2dAdaptor_Curve::ValueOffset(const Standard_Real theU) const void Geom2dAdaptor_Curve::D0(const Standard_Real U, gp_Pnt2d& P) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - D0BSpline(U, P); - else if (myTypeCurve == GeomAbs_OffsetCurve) - D0Offset(U, P); - else if (myTypeCurve == GeomAbs_BezierCurve) // use cached data - myCurveCache->D0(U, P); - else - myCurve->D0(U, P); -} - -//======================================================================= -//function : D0BSpline -//purpose : Computes the point of parameter theU on the B-spline curve -//======================================================================= -void Geom2dAdaptor_Curve::D0BSpline(const Standard_Real theU, gp_Pnt2d& theP) const -{ - if (theU == myFirst || theU == myLast) + switch (myTypeCurve) { - Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); + aBspl->LocalD0(U, aStart, aFinish, P); } - if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + else if (!myCurveCache.IsNull()) // use cached data + { + if (!myCurveCache->IsCacheValid(U)) + RebuildCache(U); + myCurveCache->D0(U, P); } - aBspl->LocalD0(theU, Ideb, Ifin, theP); - return; - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - myCurveCache->D0(theU, theP); - return; + else + myCurve->D0(U, P); + break; } - myCurve->D0(theU, theP); -} -//======================================================================= -//function : D0Offset -//purpose : Computes the point of parameter theU on the offset curve -//======================================================================= -void Geom2dAdaptor_Curve::D0Offset(const Standard_Real theU, gp_Pnt2d& theP) const -{ - theP = ValueOffset(theU); + case GeomAbs_OffsetCurve: + myNestedEvaluator->D0(U, P); + break; + + default: + myCurve->D0(U, P); + } } //======================================================================= @@ -730,69 +674,35 @@ void Geom2dAdaptor_Curve::D0Offset(const Standard_Real theU, gp_Pnt2d& theP) con void Geom2dAdaptor_Curve::D1(const Standard_Real U, gp_Pnt2d& P, gp_Vec2d& V) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - D1BSpline(U, P, V); - else if (myTypeCurve == GeomAbs_OffsetCurve) - D1Offset(U, P, V); - else if (myTypeCurve == GeomAbs_BezierCurve) // use cached data - myCurveCache->D1(U, P, V); - else - myCurve->D1(U, P, V); -} - -//======================================================================= -//function : D1BSpline -//purpose : Computes the point of parameter theU on the B-spline curve and its derivative -//======================================================================= -void Geom2dAdaptor_Curve::D1BSpline(const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d& theV) const -{ - if (theU == myFirst || theU == myLast) + switch (myTypeCurve) { - Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); + aBspl->LocalD1(U, aStart, aFinish, P, V); } - if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + else if (!myCurveCache.IsNull()) // use cached data + { + if (!myCurveCache->IsCacheValid(U)) + RebuildCache(U); + myCurveCache->D1(U, P, V); } - aBspl->LocalD1(theU, Ideb, Ifin, theP, theV); - return; - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - myCurveCache->D1(theU, theP, theV); - return; + else + myCurve->D1(U, P, V); + break; } - myCurve->D1(theU, theP, theV); -} - -//======================================================================= -//function : D1Offset -//purpose : Computes the point of parameter theU on the offset curve and its derivative -//======================================================================= -void Geom2dAdaptor_Curve::D1Offset(const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d& theV) const -{ - // 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_Vec2d V2; - myOffsetBaseCurveAdaptor->D2 (theU, theP, theV, V2); - - Standard_Boolean IsDirectionChange = Standard_False; - if(theV.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 2, theU, theV, V2); + case GeomAbs_OffsetCurve: + myNestedEvaluator->D1(U, P, V); + break; - Standard_Real anOffset = Handle(Geom2d_OffsetCurve)::DownCast(myCurve)->Offset(); - CSLib_Offset::D1(theP, theV, V2, anOffset, IsDirectionChange, theP, theV); + default: + myCurve->D1(U, P, V); + } } //======================================================================= @@ -803,73 +713,35 @@ void Geom2dAdaptor_Curve::D1Offset(const Standard_Real theU, gp_Pnt2d& theP, gp_ void Geom2dAdaptor_Curve::D2(const Standard_Real U, gp_Pnt2d& P, gp_Vec2d& V1, gp_Vec2d& V2) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - D2BSpline(U, P, V1, V2); - else if (myTypeCurve == GeomAbs_OffsetCurve) - D2Offset(U, P, V1, V2); - else if (myTypeCurve == GeomAbs_BezierCurve) // use cached data - myCurveCache->D2(U, P, V1, V2); - else - myCurve->D2(U, P, V1, V2); -} - -//======================================================================= -//function : D2BSpline -//purpose : Computes the point of parameter theU on the B-spline curve and its first and second derivatives -//======================================================================= -void Geom2dAdaptor_Curve::D2BSpline(const Standard_Real theU, gp_Pnt2d& theP, - gp_Vec2d& theV1, gp_Vec2d& theV2) const -{ - if (theU == myFirst || theU == myLast) + switch (myTypeCurve) { - Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); + aBspl->LocalD2(U, aStart, aFinish, P, V1, V2); } - if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + else if (!myCurveCache.IsNull()) // use cached data + { + if (!myCurveCache->IsCacheValid(U)) + RebuildCache(U); + myCurveCache->D2(U, P, V1, V2); } - aBspl->LocalD2(theU, Ideb, Ifin, theP, theV1, theV2); - return; - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - myCurveCache->D2(theU, theP, theV1, theV2); - return; + else + myCurve->D2(U, P, V1, V2); + break; } - myCurve->D2(theU, theP, theV1, theV2); -} -//======================================================================= -//function : D2Offset -//purpose : Computes the point of parameter theU on the offset curve and its first and second derivatives -//======================================================================= -void Geom2dAdaptor_Curve::D2Offset(const Standard_Real theU, gp_Pnt2d& theP, - gp_Vec2d& theV1, gp_Vec2d& theV2) const -{ - // 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_Vec2d V3; - myOffsetBaseCurveAdaptor->D3 (theU, theP, theV1, theV2, V3); - - Standard_Boolean IsDirectionChange = Standard_False; - if(theV1.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 3, theU, theV1, theV2, V3); + case GeomAbs_OffsetCurve: + myNestedEvaluator->D2(U, P, V1, V2); + break; - Standard_Real anOffset = Handle(Geom2d_OffsetCurve)::DownCast(myCurve)->Offset(); - CSLib_Offset::D2(theP, theV1, theV2, V3, anOffset, IsDirectionChange, theP, theV1, theV2); + default: + myCurve->D2(U, P, V1, V2); + } } //======================================================================= @@ -881,80 +753,35 @@ void Geom2dAdaptor_Curve::D3(const Standard_Real U, gp_Pnt2d& P, gp_Vec2d& V1, gp_Vec2d& V2, gp_Vec2d& V3) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - D3BSpline(U, P, V1, V2, V3); - else if (myTypeCurve == GeomAbs_OffsetCurve) - D3Offset(U, P, V1, V2, V3); - else if (myTypeCurve == GeomAbs_BezierCurve) // use cached data - myCurveCache->D3(U, P, V1, V2, V3); - else - myCurve->D3(U, P, V1, V2, V3); -} - -//======================================================================= -//function : D3BSpline -//purpose : Computes the point of parameter theU on the B-spline curve and its 1st - 3rd derivatives -//======================================================================= -void Geom2dAdaptor_Curve::D3BSpline(const Standard_Real theU, gp_Pnt2d& theP, - gp_Vec2d& theV1, gp_Vec2d& theV2, gp_Vec2d& theV3) const -{ - if (theU == myFirst || theU == myLast) + switch (myTypeCurve) { - Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); + aBspl->LocalD3(U, aStart, aFinish, P, V1, V2, V3); } - if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + else if (!myCurveCache.IsNull()) // use cached data + { + if (!myCurveCache->IsCacheValid(U)) + RebuildCache(U); + myCurveCache->D3(U, P, V1, V2, V3); } - aBspl->LocalD3(theU, Ideb, Ifin, theP, theV1, theV2, theV3); - return; - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - myCurveCache->D3(theU, theP, theV1, theV2, theV3); - return; + else + myCurve->D3(U, P, V1, V2, V3); + break; } - myCurve->D3(theU, theP, theV1, theV2, theV3); -} -//======================================================================= -//function : D3Offset -//purpose : Computes the point of parameter theU on the offset curve and its 1st - 3rd derivatives -//======================================================================= -void Geom2dAdaptor_Curve::D3Offset(const Standard_Real theU, gp_Pnt2d& theP, - gp_Vec2d& theV1, gp_Vec2d& theV2, gp_Vec2d& theV3) const -{ - // 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 - - Standard_Boolean IsDirectionChange = Standard_False; - - myOffsetBaseCurveAdaptor->D3 (theU, theP, theV1, theV2, theV3); - gp_Vec2d V4 = myOffsetBaseCurveAdaptor->DN (theU, 4); - - if(theV1.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 4, theU, theV1, theV2, theV3, V4); + case GeomAbs_OffsetCurve: + myNestedEvaluator->D3(U, P, V1, V2, V3); + break; - Standard_Real anOffset = Handle(Geom2d_OffsetCurve)::DownCast(myCurve)->Offset(); - CSLib_Offset::D3(theP, theV1, theV2, theV3, V4, anOffset, IsDirectionChange, - theP, theV1, theV2, theV3); + default: + myCurve->D3(U, P, V1, V2, V3); + } } //======================================================================= @@ -965,58 +792,30 @@ void Geom2dAdaptor_Curve::D3Offset(const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d Geom2dAdaptor_Curve::DN(const Standard_Real U, const Standard_Integer N) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - return DNBSpline(U, N); - else if (myTypeCurve == GeomAbs_OffsetCurve) - return DNOffset(U, N); - - return myCurve->DN(U, N); -} - -gp_Vec2d Geom2dAdaptor_Curve::DNBSpline(const Standard_Real U, - const Standard_Integer N) const -{ - if (U==myFirst || U==myLast) + switch (myTypeCurve) { - Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (U==myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom2d_BSplineCurve) aBspl = Handle(Geom2d_BSplineCurve)::DownCast(myCurve); + return aBspl->LocalDN(U, aStart, aFinish, N); } - if (U==myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; - } - return aBspl->LocalDN( U, Ideb, Ifin, N); + else + return myCurve->DN(U, N); + break; } - return myCurve->DN( U, N); -} - -gp_Vec2d Geom2dAdaptor_Curve::DNOffset(const Standard_Real U, - const Standard_Integer N) const -{ - gp_Pnt2d aPnt; - gp_Vec2d aVec, aVN; - - switch (N) - { - case 1: - D1Offset(U, aPnt, aVN); - break; - case 2: - D2Offset(U, aPnt, aVec, aVN); + case GeomAbs_OffsetCurve: + return myNestedEvaluator->DN(U, N); break; - case 3: - D3Offset(U, aPnt, aVec, aVec, aVN); + + default: // to eliminate gcc warning break; - default: - aVN = myCurve->DN(U, N); } - return aVN; + return myCurve->DN(U, N); } //======================================================================= @@ -1233,57 +1032,3 @@ Standard_Integer Geom2dAdaptor_Curve::NbSamples() const { return nbPoints(myCurve); } - - -// ============= Auxiliary functions =================== -Standard_Boolean AdjustDerivative(const Handle(Adaptor2d_HCurve2d)& theAdaptor, Standard_Integer theMaxDerivative, - Standard_Real theU, gp_Vec2d& theD1, gp_Vec2d& theD2, - gp_Vec2d& theD3, gp_Vec2d& theD4) -{ - static const Standard_Real aTol = gp::Resolution(); - - Standard_Boolean IsDirectionChange = Standard_False; - const Standard_Real anUinfium = theAdaptor->FirstParameter(); - const Standard_Real anUsupremum = theAdaptor->LastParameter(); - - 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 = Max(du * DivisionFactor, MinStep); - - //Derivative is approximated by Taylor-series - Standard_Integer anIndex = 1; //Derivative order - gp_Vec2d V; - - do - { - V = theAdaptor->DN(theU, ++anIndex); - } - while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder); - - Standard_Real u; - - if(theU-anUinfium < aDelta) - u = theU+aDelta; - else - u = theU-aDelta; - - gp_Pnt2d P1, P2; - theAdaptor->D0(Min(theU, u),P1); - theAdaptor->D0(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]) = theAdaptor->DN(theU, anIndex + i) * aSign; - - return IsDirectionChange; -} diff --git a/src/Geom2dAdaptor/Geom2dAdaptor_Curve.hxx b/src/Geom2dAdaptor/Geom2dAdaptor_Curve.hxx index ac60255103..89b51cd3e0 100644 --- a/src/Geom2dAdaptor/Geom2dAdaptor_Curve.hxx +++ b/src/Geom2dAdaptor/Geom2dAdaptor_Curve.hxx @@ -29,6 +29,8 @@ #include #include #include +#include + class Geom2d_Curve; class Adaptor2d_HCurve2d; class Standard_NoSuchObject; @@ -174,68 +176,24 @@ protected: private: - - //! Computes the point of parameter U on the B-spline curve - Standard_EXPORT gp_Pnt2d ValueBSpline (const Standard_Real U) const; - - //! Computes the point of parameter U on the offset curve - Standard_EXPORT gp_Pnt2d ValueOffset (const Standard_Real U) const; - - //! Computes the point of parameter U on the B-spline curve - Standard_EXPORT void D0BSpline (const Standard_Real theU, gp_Pnt2d& theP) const; - - //! Computes the point of parameter U on the offset curve - Standard_EXPORT void D0Offset (const Standard_Real theU, gp_Pnt2d& theP) const; - - //! Computes the point of parameter U on the B-spline curve - //! and its derivative - Standard_EXPORT void D1BSpline (const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d& theV) const; - - //! Computes the point of parameter U on the offset curve - //! and its derivative - Standard_EXPORT void D1Offset (const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d& theV) const; - - //! Computes the point of parameter U on the B-spline curve - //! and its first and second derivatives - Standard_EXPORT void D2BSpline (const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d& theV1, gp_Vec2d& theV2) const; - - //! Computes the point of parameter U on the offset curve - //! and its first and second derivatives - Standard_EXPORT void D2Offset (const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d& theV1, gp_Vec2d& theV2) const; - - //! Computes the point of parameter U on the B-spline curve - //! and its first, second and third derivatives - Standard_EXPORT void D3BSpline (const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d& theV1, gp_Vec2d& theV2, gp_Vec2d& theV3) const; - - //! Computes the point of parameter U on the offset curve - //! and its first, second and third derivatives - Standard_EXPORT void D3Offset (const Standard_Real theU, gp_Pnt2d& theP, gp_Vec2d& theV1, gp_Vec2d& theV2, gp_Vec2d& theV3) const; - - - //! The returned vector gives the value of the derivative for the - //! order of derivation N. - Standard_EXPORT gp_Vec2d DNBSpline (const Standard_Real theU, const Standard_Integer N) const; - - - //! The returned vector gives the value of the derivative for the - //! order of derivation N. - Standard_EXPORT gp_Vec2d DNOffset (const Standard_Real theU, const Standard_Integer N) const; - Standard_EXPORT GeomAbs_Shape LocalContinuity (const Standard_Real U1, const Standard_Real U2) const; Standard_EXPORT void load (const Handle(Geom2d_Curve)& C, const Standard_Real UFirst, const Standard_Real ULast); - + + //! Check theU relates to start or finish point of B-spline curve and return indices of span the point is located + Standard_Boolean IsBoundary(const Standard_Real theU, Standard_Integer& theSpanStart, Standard_Integer& theSpanFinish) const; + //! Rebuilds B-spline cache //! \param theParameter the value on the knot axis which identifies the caching span - Standard_EXPORT void RebuildCache (const Standard_Real theParameter) const; + void RebuildCache (const Standard_Real theParameter) const; Handle(Geom2d_Curve) myCurve; GeomAbs_CurveType myTypeCurve; Standard_Real myFirst; Standard_Real myLast; - Handle(BSplCLib_Cache) myCurveCache; - Handle(Adaptor2d_HCurve2d) myOffsetBaseCurveAdaptor; + Handle(BSplCLib_Cache) myCurveCache; ///< Cached data for B-spline or Bezier curve + Handle(Geom2dEvaluator_Curve) myNestedEvaluator; ///< Calculates value of offset curve }; diff --git a/src/Geom2dEvaluator/FILES b/src/Geom2dEvaluator/FILES new file mode 100644 index 0000000000..3eb8ed1b15 --- /dev/null +++ b/src/Geom2dEvaluator/FILES @@ -0,0 +1,3 @@ +Geom2dEvaluator_Curve.hxx +Geom2dEvaluator_OffsetCurve.cxx +Geom2dEvaluator_OffsetCurve.hxx diff --git a/src/Geom2dEvaluator/Geom2dEvaluator_Curve.hxx b/src/Geom2dEvaluator/Geom2dEvaluator_Curve.hxx new file mode 100644 index 0000000000..3d6d0d64b6 --- /dev/null +++ b/src/Geom2dEvaluator/Geom2dEvaluator_Curve.hxx @@ -0,0 +1,53 @@ +// Created on: 2015-09-21 +// Copyright (c) 2015 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 _Geom2dEvaluator_Curve_HeaderFile +#define _Geom2dEvaluator_Curve_HeaderFile + +#include +#include + +class gp_Pnt2d; +class gp_Vec2d; + +//! Interface for calculation of values and derivatives for different kinds of curves in 2D. +//! Works both with adaptors and curves. +class Geom2dEvaluator_Curve : public Standard_Transient +{ +public: + Geom2dEvaluator_Curve() {} + + //! Value of 2D curve + virtual void D0(const Standard_Real theU, + gp_Pnt2d& theValue) const = 0; + //! Value and first derivatives of curve + virtual void D1(const Standard_Real theU, + gp_Pnt2d& theValue, gp_Vec2d& theD1) const = 0; + //! Value, first and second derivatives of curve + virtual void D2(const Standard_Real theU, + gp_Pnt2d& theValue, gp_Vec2d& theD1, gp_Vec2d& theD2) const = 0; + //! Value, first, second and third derivatives of curve + virtual void D3(const Standard_Real theU, + gp_Pnt2d& theValue, gp_Vec2d& theD1, gp_Vec2d& theD2, gp_Vec2d& theD3) const = 0; + //! Calculates N-th derivatives of curve, where N = theDerU. Raises if N < 1 + virtual gp_Vec2d DN(const Standard_Real theU, + const Standard_Integer theDerU) const = 0; + + DEFINE_STANDARD_RTTI(Geom2dEvaluator_Curve, Standard_Transient) +}; + +DEFINE_STANDARD_HANDLE(Geom2dEvaluator_Curve, Standard_Transient) + + +#endif // _Geom2dEvaluator_Curve_HeaderFile diff --git a/src/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.cxx b/src/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.cxx new file mode 100644 index 0000000000..b60cf57c59 --- /dev/null +++ b/src/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.cxx @@ -0,0 +1,450 @@ +// Created on: 2015-09-21 +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include + + +Geom2dEvaluator_OffsetCurve::Geom2dEvaluator_OffsetCurve( + const Handle(Geom2d_Curve)& theBase, + const Standard_Real theOffset) + : Geom2dEvaluator_Curve(), + myBaseCurve(theBase), + myOffset(theOffset) +{ +} + +Geom2dEvaluator_OffsetCurve::Geom2dEvaluator_OffsetCurve( + const Handle(Geom2dAdaptor_HCurve)& theBase, + const Standard_Real theOffset) + : Geom2dEvaluator_Curve(), + myBaseAdaptor(theBase), + myOffset(theOffset) +{ +} + +void Geom2dEvaluator_OffsetCurve::D0(const Standard_Real theU, + gp_Pnt2d& theValue) const +{ + gp_Vec2d aD1; + BaseD1(theU, theValue, aD1); + CalculateD0(theValue, aD1); +} + +void Geom2dEvaluator_OffsetCurve::D1(const Standard_Real theU, + gp_Pnt2d& theValue, + gp_Vec2d& theD1) const +{ + gp_Vec2d aD2; + BaseD2(theU, theValue, theD1, aD2); + CalculateD1(theValue, theD1, aD2); +} + +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()) + { + gp_Vec2d aDummyD4; + isDirectionChange = AdjustDerivative(3, theU, theD1, theD2, aD3, aDummyD4); + } + + CalculateD2(theValue, theD1, theD2, aD3, isDirectionChange); +} + +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); + + 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); +} + +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) + { + 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); + } + return aDN; +} + + +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); +} + + +void Geom2dEvaluator_OffsetCurve::CalculateD0( gp_Pnt2d& theValue, + const gp_Vec2d& theD1) const +{ + if (theD1.SquareMagnitude() <= gp::Resolution()) + Standard_NullValue::Raise("Geom2dEvaluator_OffsetCurve: Undefined normal vector " + "because tangent vector has zero-magnitude!"); + + gp_Dir2d aNormal(theD1.Y(), -theD1.X()); + theValue.ChangeCoord().Add(aNormal.XY() * myOffset); +} + +void Geom2dEvaluator_OffsetCurve::CalculateD1( gp_Pnt2d& theValue, + gp_Vec2d& theD1, + const gp_Vec2d& theD2) const +{ + // 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()); + Standard_Real R2 = Ndir.SquareModulus(); + Standard_Real R = Sqrt(R2); + Standard_Real R3 = R * R2; + Standard_Real Dr = Ndir.Dot(DNdir); + if (R3 <= gp::Resolution()) + { + if (R2 <= gp::Resolution()) + Standard_NullValue::Raise("Geom2dEvaluator_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_Vec2d(DNdir)); +} + +void Geom2dEvaluator_OffsetCurve::CalculateD2( gp_Pnt2d& theValue, + gp_Vec2d& theD1, + gp_Vec2d& theD2, + const gp_Vec2d& theD3, + const Standard_Boolean theIsDirChange) const +{ + // 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()); + Standard_Real R2 = Ndir.SquareModulus(); + Standard_Real R = 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()) + Standard_NullValue::Raise("Geom2dEvaluator_OffsetCurve: Null derivative"); + //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(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_Vec2d(DNdir)); + // P"(u) : + if (theIsDirChange) + theD2.Reverse(); + theD2.Add(gp_Vec2d(D2Ndir)); +} + +void Geom2dEvaluator_OffsetCurve::CalculateD3( gp_Pnt2d& theValue, + gp_Vec2d& theD1, + gp_Vec2d& theD2, + gp_Vec2d& theD3, + const gp_Vec2d& theD4, + const Standard_Boolean theIsDirChange) const +{ + // 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()); + Standard_Real R2 = Ndir.SquareModulus(); + Standard_Real R = 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()) + Standard_NullValue::Raise("Geom2dEvaluator_OffsetCurve: Null derivative"); + //We try another computation but the stability is not very good dixit ISG. + // V3 = P"' (U) : + D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * myOffset * Dr / R2)); + D3Ndir.Subtract( + (DNdir.Multiplied((3.0 * myOffset) * ((D2r / R2) + (Dr*Dr) / R4)))); + D3Ndir.Add(Ndir.Multiplied( + (myOffset * (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 + { + // Same computation as IICURV in EUCLID-IS because the stability is better. + // V3 = P"' (U) : + D3Ndir.Multiply(myOffset / R); + D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * myOffset * Dr / R3)); + D3Ndir.Subtract(DNdir.Multiplied( + ((3.0 * myOffset) * ((D2r / R3) + (Dr*Dr) / R5)))); + D3Ndir.Add(Ndir.Multiplied( + (myOffset * (6.0*Dr*Dr / R5 + 6.0*Dr*D2r / R5 - 15.0*Dr*Dr*Dr / R7 - D3r)))); + // V2 = P" (U) : + D2Ndir.Multiply(myOffset / R); + D2Ndir.Subtract(DNdir.Multiplied(2.0 * myOffset * Dr / R3)); + D2Ndir.Subtract(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_Vec2d(DNdir)); + // P"(u) + theD2.Add(gp_Vec2d(D2Ndir)); + // P"'(u) + if (theIsDirChange) + theD3.Reverse(); + theD3.Add(gp_Vec2d(D2Ndir)); +} + + +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 = 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(Min(theU, u), P1); + BaseD0(Max(theU, u), P2); + + gp_Vec2d V1(P1, P2); + isDirectionChange = V.Dot(V1) < 0.0; + Standard_Real aSign = isDirectionChange ? -1.0 : 1.0; + + theD1 = V * aSign; + gp_Vec2d* aDeriv[3] = { &theD2, &theD3, &theD4 }; + for (Standard_Integer i = 1; i < theMaxDerivative; i++) + *(aDeriv[i - 1]) = BaseDN(theU, anIndex + i) * aSign; + + return isDirectionChange; +} + diff --git a/src/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.hxx b/src/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.hxx new file mode 100644 index 0000000000..b9209da516 --- /dev/null +++ b/src/Geom2dEvaluator/Geom2dEvaluator_OffsetCurve.hxx @@ -0,0 +1,117 @@ +// Created on: 2015-09-21 +// Copyright (c) 2015 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 _Geom2dEvaluator_OffsetCurve_HeaderFile +#define _Geom2dEvaluator_OffsetCurve_HeaderFile + +#include +#include + +class Geom2dAdaptor_HCurve; + +//! Allows to calculate values and derivatives for offset curves in 2D +class Geom2dEvaluator_OffsetCurve : public Geom2dEvaluator_Curve +{ +public: + //! Initialize evaluator by curve + Standard_EXPORT Geom2dEvaluator_OffsetCurve( + const Handle(Geom2d_Curve)& theBase, + const Standard_Real theOffset); + //! Initialize evaluator by curve adaptor + Standard_EXPORT Geom2dEvaluator_OffsetCurve( + const Handle(Geom2dAdaptor_HCurve)& theBase, + const Standard_Real theOffset); + + //! Change the offset value + void SetOffsetValue(Standard_Real theOffset) + { myOffset = theOffset; } + + //! Value of curve + Standard_EXPORT void D0(const Standard_Real theU, + gp_Pnt2d& theValue) const Standard_OVERRIDE; + //! Value and first derivatives of curve + Standard_EXPORT void D1(const Standard_Real theU, + gp_Pnt2d& theValue, gp_Vec2d& theD1) const Standard_OVERRIDE; + //! Value, first and second derivatives of curve + Standard_EXPORT void D2(const Standard_Real theU, + gp_Pnt2d& theValue, gp_Vec2d& theD1, gp_Vec2d& theD2) const Standard_OVERRIDE; + //! Value, first, second and third derivatives of curve + Standard_EXPORT void D3(const Standard_Real theU, + gp_Pnt2d& theValue, gp_Vec2d& theD1, + gp_Vec2d& theD2, gp_Vec2d& theD3) const Standard_OVERRIDE; + //! Calculates N-th derivatives of curve, where N = theDeriv. Raises if N < 1 + Standard_EXPORT gp_Vec2d DN(const Standard_Real theU, + const Standard_Integer theDeriv) const Standard_OVERRIDE; + + DEFINE_STANDARD_RTTI(Geom2dEvaluator_OffsetCurve, Geom2dEvaluator_Curve) + +private: + //! Recalculate D1 values of base curve into D0 value of offset curve + void CalculateD0( gp_Pnt2d& theValue, + const gp_Vec2d& theD1) const; + //! Recalculate D2 values of base curve into D1 values of offset curve + void CalculateD1( gp_Pnt2d& theValue, + gp_Vec2d& theD1, + const gp_Vec2d& theD2) const; + //! Recalculate D3 values of base curve into D2 values of offset curve + void CalculateD2( gp_Pnt2d& theValue, + gp_Vec2d& theD1, + gp_Vec2d& theD2, + const gp_Vec2d& theD3, + const Standard_Boolean theIsDirChange) const; + //! Recalculate D3 values of base curve into D3 values of offset curve + void CalculateD3( gp_Pnt2d& theValue, + gp_Vec2d& theD1, + gp_Vec2d& theD2, + gp_Vec2d& theD3, + const gp_Vec2d& theD4, + const Standard_Boolean theIsDirChange) const; + + //! 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_HCurve) myBaseAdaptor; + + Standard_Real myOffset; ///< offset value +}; + +DEFINE_STANDARD_HANDLE(Geom2dEvaluator_OffsetCurve, Geom2dEvaluator_Curve) + + +#endif // _Geom2dEvaluator_OffsetCurve_HeaderFile diff --git a/src/GeomAdaptor/GeomAdaptor_Curve.cxx b/src/GeomAdaptor/GeomAdaptor_Curve.cxx index e9a3110192..56e8af359a 100644 --- a/src/GeomAdaptor/GeomAdaptor_Curve.cxx +++ b/src/GeomAdaptor/GeomAdaptor_Curve.cxx @@ -25,7 +25,6 @@ #include #include #include -#include #include #include #include @@ -40,6 +39,7 @@ #include #include #include +#include #include #include #include @@ -62,16 +62,6 @@ //#include static const Standard_Real PosTol = Precision::PConfusion() / 2; -static const Standard_Integer maxDerivOrder = 3; -static const Standard_Real MinStep = 1e-7; - -static gp_Vec dummyDerivative; // used as empty value for unused derivatives in AdjustDerivative -// Recalculate derivatives in the singular point -// Returns true if the direction of derivatives is changed -static Standard_Boolean AdjustDerivative( - const Handle(Adaptor3d_HCurve)& theAdaptor, Standard_Integer theMaxDerivative, Standard_Real theU, gp_Vec& theD1, - gp_Vec& theD2 = dummyDerivative, gp_Vec& theD3 = dummyDerivative, gp_Vec& theD4 = dummyDerivative); - //======================================================================= //function : LocalContinuity @@ -153,6 +143,7 @@ void GeomAdaptor_Curve::load(const Handle(Geom_Curve)& C, if ( myCurve != C) { myCurve = C; myCurveCache = Handle(BSplCLib_Cache)(); + myNestedEvaluator = Handle(GeomEvaluator_Curve)(); const Handle(Standard_Type)& TheType = C->DynamicType(); if ( TheType == STANDARD_TYPE(Geom_TrimmedCurve)) { @@ -191,9 +182,12 @@ void GeomAdaptor_Curve::load(const Handle(Geom_Curve)& C, } else if ( TheType == STANDARD_TYPE(Geom_OffsetCurve)) { myTypeCurve = GeomAbs_OffsetCurve; + Handle(Geom_OffsetCurve) anOffsetCurve = Handle(Geom_OffsetCurve)::DownCast(myCurve); // Create nested adaptor for base curve - Handle(Geom_Curve) aBase = Handle(Geom_OffsetCurve)::DownCast(myCurve)->BasisCurve(); - myOffsetBaseCurveAdaptor = new GeomAdaptor_HCurve(aBase); + Handle(Geom_Curve) aBaseCurve = anOffsetCurve->BasisCurve(); + Handle(GeomAdaptor_HCurve) aBaseAdaptor = new GeomAdaptor_HCurve(aBaseCurve); + myNestedEvaluator = new GeomEvaluator_OffsetCurve( + aBaseAdaptor, anOffsetCurve->Offset(), anOffsetCurve->Direction()); } else { myTypeCurve = GeomAbs_OtherCurve; @@ -571,77 +565,47 @@ void GeomAdaptor_Curve::RebuildCache(const Standard_Real theParameter) const } //======================================================================= -//function : Value -//purpose : -//======================================================================= - -gp_Pnt GeomAdaptor_Curve::Value(const Standard_Real U) const -{ - if (myTypeCurve == GeomAbs_BSplineCurve) - return ValueBSpline(U); - else if (myTypeCurve == GeomAbs_OffsetCurve) - return ValueOffset(U); - else if (myTypeCurve == GeomAbs_BezierCurve) - { // use cached data - gp_Pnt aRes; - myCurveCache->D0(U, aRes); - return aRes; - } - return myCurve->Value(U); -} - -//======================================================================= -//function : ValueBSpline +//function : IsBoundary //purpose : //======================================================================= -gp_Pnt GeomAdaptor_Curve::ValueBSpline(const Standard_Real theU) const +Standard_Boolean GeomAdaptor_Curve::IsBoundary(const Standard_Real theU, + Standard_Integer& theSpanStart, + Standard_Integer& theSpanFinish) const { - if (theU == myFirst || theU == myLast) + Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); + if (!aBspl.IsNull() && (theU == myFirst || theU == myLast)) { - Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + if (theU == myFirst) + { + aBspl->LocateU(myFirst, PosTol, theSpanStart, theSpanFinish); + if (theSpanStart < 1) + theSpanStart = 1; + if (theSpanStart >= theSpanFinish) + theSpanFinish = theSpanStart + 1; } - if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + else if (theU == myLast) + { + aBspl->LocateU(myLast, PosTol, theSpanStart, theSpanFinish); + if (theSpanFinish > aBspl->NbKnots()) + theSpanFinish = aBspl->NbKnots(); + if (theSpanStart >= theSpanFinish) + theSpanStart = theSpanFinish - 1; } - return aBspl->LocalValue(theU, Ideb, Ifin); - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - gp_Pnt aRes; - myCurveCache->D0(theU, aRes); - return aRes; + return Standard_True; } - return myCurve->Value(theU); + return Standard_False; } //======================================================================= -//function : ValueOffset +//function : Value //purpose : //======================================================================= -gp_Pnt GeomAdaptor_Curve::ValueOffset(const Standard_Real theU) const + +gp_Pnt GeomAdaptor_Curve::Value(const Standard_Real U) const { - gp_Pnt aP; - gp_Vec aV; - myOffsetBaseCurveAdaptor->D1(theU, aP, aV); - Standard_Boolean IsDirectionChange = Standard_False; - if(aV.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 1, theU, aV); - - Handle(Geom_OffsetCurve) anOffC = Handle(Geom_OffsetCurve)::DownCast(myCurve); - Standard_Real anOffsetVal = anOffC->Offset(); - const gp_Dir& anOffsetDir = anOffC->Direction(); - - CSLib_Offset::D0(aP, aV, anOffsetDir, anOffsetVal, IsDirectionChange, aP); - return aP; + gp_Pnt aValue; + D0(U, aValue); + return aValue; } //======================================================================= @@ -651,56 +615,35 @@ gp_Pnt GeomAdaptor_Curve::ValueOffset(const Standard_Real theU) const void GeomAdaptor_Curve::D0(const Standard_Real U, gp_Pnt& P) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - D0BSpline(U, P); - else if (myTypeCurve == GeomAbs_OffsetCurve) - D0Offset(U, P); - else if (myTypeCurve == GeomAbs_BezierCurve) // use cached data - myCurveCache->D0(U, P); - else - myCurve->D0(U, P); -} - -//======================================================================= -//function : D0BSpline -//purpose : -//======================================================================= -void GeomAdaptor_Curve::D0BSpline(const Standard_Real theU, gp_Pnt& theP) const -{ - if (theU == myFirst || theU == myLast) + switch (myTypeCurve) { - Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); + aBspl->LocalD0(U, aStart, aFinish, P); } - if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + else if (!myCurveCache.IsNull()) // use cached data + { + if (!myCurveCache->IsCacheValid(U)) + RebuildCache(U); + myCurveCache->D0(U, P); } - aBspl->LocalD0(theU, Ideb, Ifin, theP); - return; - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - myCurveCache->D0(theU, theP); - return; + else + myCurve->D0(U, P); + break; } - myCurve->D0(theU, theP); -} -//======================================================================= -//function : D0Offset -//purpose : -//======================================================================= -void GeomAdaptor_Curve::D0Offset(const Standard_Real theU, gp_Pnt& theP) const -{ - theP = ValueOffset(theU); + case GeomAbs_OffsetCurve: + myNestedEvaluator->D0(U, P); + break; + + default: + myCurve->D0(U, P); + } } //======================================================================= @@ -710,69 +653,37 @@ void GeomAdaptor_Curve::D0Offset(const Standard_Real theU, gp_Pnt& theP) const void GeomAdaptor_Curve::D1(const Standard_Real U, gp_Pnt& P, gp_Vec& V) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - D1BSpline(U, P, V); - else if (myTypeCurve == GeomAbs_OffsetCurve) - D1Offset(U, P, V); - else if (myTypeCurve == GeomAbs_BezierCurve) // use cached data - myCurveCache->D1(U, P, V); - else - myCurve->D1(U, P, V); -} - -//======================================================================= -//function : D1BSpline -//purpose : -//======================================================================= -void GeomAdaptor_Curve::D1BSpline(const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV) const -{ - if (theU == myFirst || theU == myLast) + switch (myTypeCurve) { - Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); + aBspl->LocalD1(U, aStart, aFinish, P, V); } - if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + else if (!myCurveCache.IsNull()) // use cached data + { + if (!myCurveCache->IsCacheValid(U)) + RebuildCache(U); + myCurveCache->D1(U, P, V); } - aBspl->LocalD1(theU, Ideb, Ifin, theP, theV); - return; - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - myCurveCache->D1(theU, theP, theV); - return; + else + myCurve->D1(U, P, V); + break; } - myCurve->D1(theU, theP, theV); -} - -//======================================================================= -//function : D1Offset -//purpose : -//======================================================================= -void GeomAdaptor_Curve::D1Offset(const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV) const -{ - gp_Vec aV2; - myOffsetBaseCurveAdaptor->D2 (theU, theP, theV, aV2); - Standard_Boolean IsDirectionChange = Standard_False; - if(theV.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 2, theU, theV, aV2); + case GeomAbs_OffsetCurve: + myNestedEvaluator->D1(U, P, V); + break; - Handle(Geom_OffsetCurve) anOffC = Handle(Geom_OffsetCurve)::DownCast(myCurve); - Standard_Real anOffsetVal = anOffC->Offset(); - const gp_Dir& anOffsetDir = anOffC->Direction(); - CSLib_Offset::D1(theP, theV, aV2, anOffsetDir, anOffsetVal, IsDirectionChange, theP, theV); + default: + myCurve->D1(U, P, V); + } } - //======================================================================= //function : D2 //purpose : @@ -781,68 +692,35 @@ void GeomAdaptor_Curve::D1Offset(const Standard_Real theU, gp_Pnt& theP, gp_Vec& void GeomAdaptor_Curve::D2(const Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - D2BSpline(U, P, V1, V2); - else if (myTypeCurve == GeomAbs_OffsetCurve) - D2Offset(U, P, V1, V2); - else if (myTypeCurve == GeomAbs_BezierCurve) // use cached data - myCurveCache->D2(U, P, V1, V2); - else - myCurve->D2(U, P, V1, V2); -} - -//======================================================================= -//function : D2BSpline -//purpose : -//======================================================================= -void GeomAdaptor_Curve::D2BSpline(const Standard_Real theU, gp_Pnt& theP, - gp_Vec& theV1, gp_Vec& theV2) const -{ - if (theU == myFirst || theU == myLast) + switch (myTypeCurve) { - Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); + aBspl->LocalD2(U, aStart, aFinish, P, V1, V2); } - if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + else if (!myCurveCache.IsNull()) // use cached data + { + if (!myCurveCache->IsCacheValid(U)) + RebuildCache(U); + myCurveCache->D2(U, P, V1, V2); } - aBspl->LocalD2(theU, Ideb, Ifin, theP, theV1, theV2); - return; - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - myCurveCache->D2(theU, theP, theV1, theV2); - return; + else + myCurve->D2(U, P, V1, V2); + break; } - myCurve->D2(theU, theP, theV1, theV2); -} -//======================================================================= -//function : D2Offset -//purpose : -//======================================================================= -void GeomAdaptor_Curve::D2Offset(const Standard_Real theU, gp_Pnt& theP, - gp_Vec& theV1, gp_Vec& theV2) const -{ - gp_Vec V3; - myOffsetBaseCurveAdaptor->D3 (theU, theP, theV1, theV2, V3); - - Standard_Boolean IsDirectionChange = Standard_False; - if(theV1.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 3, theU, theV1, theV2, V3); + case GeomAbs_OffsetCurve: + myNestedEvaluator->D2(U, P, V1, V2); + break; - Handle(Geom_OffsetCurve) anOffC = Handle(Geom_OffsetCurve)::DownCast(myCurve); - Standard_Real anOffsetVal = anOffC->Offset(); - const gp_Dir& anOffsetDir = anOffC->Direction(); - CSLib_Offset::D2(theP, theV1, theV2, V3, anOffsetDir, anOffsetVal, IsDirectionChange, theP, theV1, theV2); + default: + myCurve->D2(U, P, V1, V2); + } } //======================================================================= @@ -854,71 +732,35 @@ void GeomAdaptor_Curve::D3(const Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2, gp_Vec& V3) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - D3BSpline(U, P, V1, V2, V3); - else if (myTypeCurve == GeomAbs_OffsetCurve) - D3Offset(U, P, V1, V2, V3); - else if (myTypeCurve == GeomAbs_BezierCurve) // use cached data - myCurveCache->D3(U, P, V1, V2, V3); - else - myCurve->D3(U, P, V1, V2, V3); -} - -//======================================================================= -//function : D3BSpline -//purpose : -//======================================================================= -void GeomAdaptor_Curve::D3BSpline(const Standard_Real theU, - gp_Pnt& theP, gp_Vec& theV1, - gp_Vec& theV2, gp_Vec& theV3) const -{ - if (theU == myFirst || theU == myLast) + switch (myTypeCurve) { - Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (theU == myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); + aBspl->LocalD3(U, aStart, aFinish, P, V1, V2, V3); } - if (theU == myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; + else if (!myCurveCache.IsNull()) // use cached data + { + if (!myCurveCache->IsCacheValid(U)) + RebuildCache(U); + myCurveCache->D3(U, P, V1, V2, V3); } - aBspl->LocalD3(theU, Ideb, Ifin, theP, theV1, theV2, theV3); - return; - } - else if (!myCurveCache.IsNull()) // use cached B-spline data - { - if (!myCurveCache->IsCacheValid(theU)) - RebuildCache(theU); - myCurveCache->D3(theU, theP, theV1, theV2, theV3); - return; + else + myCurve->D3(U, P, V1, V2, V3); + break; } - myCurve->D3(theU, theP, theV1, theV2, theV3); -} -//======================================================================= -//function : D3Offset -//purpose : -//======================================================================= -void GeomAdaptor_Curve::D3Offset(const Standard_Real theU, - gp_Pnt& theP, gp_Vec& theV1, - gp_Vec& theV2, gp_Vec& theV3) const -{ - myOffsetBaseCurveAdaptor->D3 (theU, theP, theV1, theV2, theV3); - gp_Vec V4 = myOffsetBaseCurveAdaptor->DN(theU, 4); - - Standard_Boolean IsDirectionChange = Standard_False; - if(theV1.SquareMagnitude() <= gp::Resolution()) - IsDirectionChange = AdjustDerivative(myOffsetBaseCurveAdaptor, 4, theU, theV1, theV2, theV3, V4); - - Handle(Geom_OffsetCurve) anOffC = Handle(Geom_OffsetCurve)::DownCast(myCurve); - Standard_Real anOffsetVal = anOffC->Offset(); - const gp_Dir& anOffsetDir = anOffC->Direction(); - CSLib_Offset::D3(theP, theV1, theV2, theV3, V4, anOffsetDir, anOffsetVal, IsDirectionChange, - theP, theV1, theV2, theV3); + case GeomAbs_OffsetCurve: + myNestedEvaluator->D3(U, P, V1, V2, V3); + break; + + default: + myCurve->D3(U, P, V1, V2, V3); + } } //======================================================================= @@ -929,57 +771,30 @@ void GeomAdaptor_Curve::D3Offset(const Standard_Real theU, gp_Vec GeomAdaptor_Curve::DN(const Standard_Real U, const Standard_Integer N) const { - if (myTypeCurve == GeomAbs_BSplineCurve) - return DNBSpline(U, N); - else if (myTypeCurve == GeomAbs_OffsetCurve) - return DNOffset(U, N); - - return myCurve->DN(U, N); -} - -gp_Vec GeomAdaptor_Curve::DNBSpline(const Standard_Real U, - const Standard_Integer N) const -{ - if ((U==myFirst || U==myLast)) + switch (myTypeCurve) { - Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); - Standard_Integer Ideb = 0, Ifin = 0; - if (U==myFirst) { - aBspl->LocateU(myFirst, PosTol, Ideb, Ifin); - if (Ideb<1) Ideb=1; - if (Ideb>=Ifin) Ifin = Ideb+1; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + { + Standard_Integer aStart = 0, aFinish = 0; + if (IsBoundary(U, aStart, aFinish)) + { + Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(myCurve); + return aBspl->LocalDN(U, aStart, aFinish, N); } - if (U==myLast) { - aBspl->LocateU(myLast, PosTol, Ideb, Ifin); - if (Ifin > aBspl->NbKnots()) Ifin = aBspl->NbKnots(); - if (Ideb>=Ifin) Ideb = Ifin-1; - } - return aBspl->LocalDN( U, Ideb, Ifin, N); + else + return myCurve->DN(U, N); + break; } - return myCurve->DN( U, N); -} -gp_Vec GeomAdaptor_Curve::DNOffset(const Standard_Real U, - const Standard_Integer N) const -{ - gp_Pnt aPnt; - gp_Vec aVec, aVN; - - switch (N) - { - case 1: - D1Offset(U, aPnt, aVN); + case GeomAbs_OffsetCurve: + return myNestedEvaluator->DN(U, N); break; - case 2: - D2Offset(U, aPnt, aVec, aVN); - break; - case 3: - D3Offset(U, aPnt, aVec, aVec, aVN); + + default: // to eliminate gcc warning break; - default: - aVN = myCurve->DN(U, N); } - return aVN; + return myCurve->DN(U, N); } //======================================================================= @@ -1164,57 +979,3 @@ Handle(Geom_BSplineCurve) GeomAdaptor_Curve::BSpline() const return Handle(Geom_BSplineCurve)::DownCast (myCurve); } - - -// ============= Auxiliary functions =================== -Standard_Boolean AdjustDerivative(const Handle(Adaptor3d_HCurve)& theAdaptor, Standard_Integer theMaxDerivative, - Standard_Real theU, gp_Vec& theD1, gp_Vec& theD2, - gp_Vec& theD3, gp_Vec& theD4) -{ - static const Standard_Real aTol = gp::Resolution(); - - Standard_Boolean IsDirectionChange = Standard_False; - const Standard_Real anUinfium = theAdaptor->FirstParameter(); - const Standard_Real anUsupremum = theAdaptor->LastParameter(); - - 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 = Max(du * DivisionFactor, MinStep); - - //Derivative is approximated by Taylor-series - Standard_Integer anIndex = 1; //Derivative order - gp_Vec V; - - do - { - V = theAdaptor->DN(theU, ++anIndex); - } - while((V.SquareMagnitude() <= aTol) && anIndex < maxDerivOrder); - - Standard_Real u; - - if(theU-anUinfium < aDelta) - u = theU+aDelta; - else - u = theU-aDelta; - - gp_Pnt P1, P2; - theAdaptor->D0(Min(theU, u), P1); - theAdaptor->D0(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]) = theAdaptor->DN(theU, anIndex + i) * aSign; - - return IsDirectionChange; -} diff --git a/src/GeomAdaptor/GeomAdaptor_Curve.hxx b/src/GeomAdaptor/GeomAdaptor_Curve.hxx index 18360d1417..f74bbdc152 100644 --- a/src/GeomAdaptor/GeomAdaptor_Curve.hxx +++ b/src/GeomAdaptor/GeomAdaptor_Curve.hxx @@ -29,6 +29,8 @@ #include #include #include +#include + class Geom_Curve; class Adaptor3d_HCurve; class Standard_NoSuchObject; @@ -213,68 +215,24 @@ protected: private: - - //! Computes the point of parameter U on the B-spline curve - Standard_EXPORT gp_Pnt ValueBSpline (const Standard_Real U) const; - - //! Computes the point of parameter U on the offset curve - Standard_EXPORT gp_Pnt ValueOffset (const Standard_Real U) const; - - //! Computes the point of parameter U on the B-spline curve - Standard_EXPORT void D0BSpline (const Standard_Real theU, gp_Pnt& theP) const; - - //! Computes the point of parameter U on the offset curve - Standard_EXPORT void D0Offset (const Standard_Real theU, gp_Pnt& theP) const; - - //! Computes the point of parameter U on the B-spline curve - //! and its derivative - Standard_EXPORT void D1BSpline (const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV) const; - - //! Computes the point of parameter U on the offset curve - //! and its derivative - Standard_EXPORT void D1Offset (const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV) const; - - //! Computes the point of parameter U on the B-spline curve - //! and its first and second derivatives - Standard_EXPORT void D2BSpline (const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV1, gp_Vec& theV2) const; - - //! Computes the point of parameter U on the offset curve - //! and its first and second derivatives - Standard_EXPORT void D2Offset (const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV1, gp_Vec& theV2) const; - - //! Computes the point of parameter U on the B-spline curve - //! and its first, second and third derivatives - Standard_EXPORT void D3BSpline (const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV1, gp_Vec& theV2, gp_Vec& theV3) const; - - //! Computes the point of parameter U on the offset curve - //! and its first, second and third derivatives - Standard_EXPORT void D3Offset (const Standard_Real theU, gp_Pnt& theP, gp_Vec& theV1, gp_Vec& theV2, gp_Vec& theV3) const; - - - //! The returned vector gives the value of the derivative for the - //! order of derivation N. - Standard_EXPORT gp_Vec DNBSpline (const Standard_Real theU, const Standard_Integer N) const; - - - //! The returned vector gives the value of the derivative for the - //! order of derivation N. - Standard_EXPORT gp_Vec DNOffset (const Standard_Real theU, const Standard_Integer N) const; - Standard_EXPORT GeomAbs_Shape LocalContinuity (const Standard_Real U1, const Standard_Real U2) const; Standard_EXPORT void load (const Handle(Geom_Curve)& C, const Standard_Real UFirst, const Standard_Real ULast); + //! Check theU relates to start or finish point of B-spline curve and return indices of span the point is located + Standard_Boolean IsBoundary(const Standard_Real theU, Standard_Integer& theSpanStart, Standard_Integer& theSpanFinish) const; + //! Rebuilds B-spline cache //! \param theParameter the value on the knot axis which identifies the caching span - Standard_EXPORT void RebuildCache (const Standard_Real theParameter) const; + void RebuildCache (const Standard_Real theParameter) const; Handle(Geom_Curve) myCurve; GeomAbs_CurveType myTypeCurve; Standard_Real myFirst; Standard_Real myLast; - Handle(BSplCLib_Cache) myCurveCache; - Handle(Adaptor3d_HCurve) myOffsetBaseCurveAdaptor; + Handle(BSplCLib_Cache) myCurveCache; ///< Cached data for B-spline or Bezier curve + Handle(GeomEvaluator_Curve) myNestedEvaluator; ///< Calculates value of offset curve }; diff --git a/src/GeomAdaptor/GeomAdaptor_Surface.hxx b/src/GeomAdaptor/GeomAdaptor_Surface.hxx index 24bfcaf3b7..6851273b8f 100644 --- a/src/GeomAdaptor/GeomAdaptor_Surface.hxx +++ b/src/GeomAdaptor/GeomAdaptor_Surface.hxx @@ -266,7 +266,7 @@ private: Standard_Real myVLast; Standard_Real myTolU; Standard_Real myTolV; - Handle(BSplSLib_Cache) mySurfaceCache; ///< Cached data for B-spline surface + Handle(BSplSLib_Cache) mySurfaceCache; ///< Cached data for B-spline or Bezier surface protected: GeomAbs_SurfaceType mySurfaceType; diff --git a/src/GeomEvaluator/FILES b/src/GeomEvaluator/FILES index f6e1bdc0f9..a042bd5265 100644 --- a/src/GeomEvaluator/FILES +++ b/src/GeomEvaluator/FILES @@ -1,3 +1,6 @@ +GeomEvaluator_Curve.hxx +GeomEvaluator_OffsetCurve.cxx +GeomEvaluator_OffsetCurve.hxx GeomEvaluator_OffsetSurface.cxx GeomEvaluator_OffsetSurface.hxx GeomEvaluator_Surface.hxx diff --git a/src/GeomEvaluator/GeomEvaluator_Curve.hxx b/src/GeomEvaluator/GeomEvaluator_Curve.hxx new file mode 100644 index 0000000000..2c8df0273a --- /dev/null +++ b/src/GeomEvaluator/GeomEvaluator_Curve.hxx @@ -0,0 +1,53 @@ +// Created on: 2015-09-21 +// Copyright (c) 2015 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 _GeomEvaluator_Curve_HeaderFile +#define _GeomEvaluator_Curve_HeaderFile + +#include +#include + +class gp_Pnt; +class gp_Vec; + +//! Interface for calculation of values and derivatives for different kinds of curves in 3D. +//! Works both with adaptors and curves. +class GeomEvaluator_Curve : public Standard_Transient +{ +public: + GeomEvaluator_Curve() {} + + //! Value of 3D curve + virtual void D0(const Standard_Real theU, + gp_Pnt& theValue) const = 0; + //! Value and first derivatives of curve + virtual void D1(const Standard_Real theU, + gp_Pnt& theValue, gp_Vec& theD1) const = 0; + //! Value, first and second derivatives of curve + virtual void D2(const Standard_Real theU, + gp_Pnt& theValue, gp_Vec& theD1, gp_Vec& theD2) const = 0; + //! Value, first, second and third derivatives of curve + virtual void D3(const Standard_Real theU, + gp_Pnt& theValue, gp_Vec& theD1, gp_Vec& theD2, gp_Vec& theD3) const = 0; + //! Calculates N-th derivatives of curve, where N = theDerU. Raises if N < 1 + virtual gp_Vec DN(const Standard_Real theU, + const Standard_Integer theDerU) const = 0; + + DEFINE_STANDARD_RTTI(GeomEvaluator_Curve, Standard_Transient) +}; + +DEFINE_STANDARD_HANDLE(GeomEvaluator_Curve, Standard_Transient) + + +#endif // _GeomEvaluator_Curve_HeaderFile diff --git a/src/GeomEvaluator/GeomEvaluator_OffsetCurve.cxx b/src/GeomEvaluator/GeomEvaluator_OffsetCurve.cxx new file mode 100644 index 0000000000..765662a41d --- /dev/null +++ b/src/GeomEvaluator/GeomEvaluator_OffsetCurve.cxx @@ -0,0 +1,447 @@ +// Created on: 2015-09-21 +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include + + +GeomEvaluator_OffsetCurve::GeomEvaluator_OffsetCurve( + const Handle(Geom_Curve)& theBase, + const Standard_Real theOffset, + const gp_Dir& theDirection) + : GeomEvaluator_Curve(), + myBaseCurve(theBase), + myOffset(theOffset), + myOffsetDir(theDirection) +{ +} + +GeomEvaluator_OffsetCurve::GeomEvaluator_OffsetCurve( + const Handle(GeomAdaptor_HCurve)& theBase, + const Standard_Real theOffset, + const gp_Dir& theDirection) + : GeomEvaluator_Curve(), + myBaseAdaptor(theBase), + myOffset(theOffset), + myOffsetDir(theDirection) +{ +} + +void GeomEvaluator_OffsetCurve::D0(const Standard_Real theU, + gp_Pnt& theValue) const +{ + gp_Vec aD1; + BaseD1(theU, theValue, aD1); + CalculateD0(theValue, aD1); +} + +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); +} + +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); + + Standard_Boolean isDirectionChange = Standard_False; + if (theD1.SquareMagnitude() <= gp::Resolution()) + { + gp_Vec aDummyD4; + isDirectionChange = AdjustDerivative(3, theU, theD1, theD2, aD3, aDummyD4); + } + + 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); + + 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); +} + +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) + { + 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); + } + return aDN; +} + + +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()) + Standard_NullValue::Raise("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 = Sqrt(R2); + Standard_Real R3 = R * R2; + Standard_Real Dr = Ndir.Dot(DNdir); + if (R3 <= gp::Resolution()) { + if (R2 <= gp::Resolution()) + Standard_NullValue::Raise("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 = 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()) + Standard_NullValue::Raise("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 = 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()) + Standard_NullValue::Raise("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 = 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(Min(theU, u), P1); + BaseD0(Max(theU, u), P2); + + gp_Vec V1(P1, P2); + isDirectionChange = V.Dot(V1) < 0.0; + Standard_Real aSign = isDirectionChange ? -1.0 : 1.0; + + theD1 = V * aSign; + gp_Vec* aDeriv[3] = { &theD2, &theD3, &theD4 }; + for (Standard_Integer i = 1; i < theMaxDerivative; i++) + *(aDeriv[i - 1]) = BaseDN(theU, anIndex + i) * aSign; + + return isDirectionChange; +} + diff --git a/src/GeomEvaluator/GeomEvaluator_OffsetCurve.hxx b/src/GeomEvaluator/GeomEvaluator_OffsetCurve.hxx new file mode 100644 index 0000000000..3c74dade7d --- /dev/null +++ b/src/GeomEvaluator/GeomEvaluator_OffsetCurve.hxx @@ -0,0 +1,124 @@ +// Created on: 2015-09-21 +// Copyright (c) 2015 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 _GeomEvaluator_OffsetCurve_HeaderFile +#define _GeomEvaluator_OffsetCurve_HeaderFile + +#include +#include +#include + +class GeomAdaptor_HCurve; + +//! Allows to calculate values and derivatives for offset curves in 3D +class GeomEvaluator_OffsetCurve : public GeomEvaluator_Curve +{ +public: + //! Initialize evaluator by curve + Standard_EXPORT GeomEvaluator_OffsetCurve( + const Handle(Geom_Curve)& theBase, + const Standard_Real theOffset, + const gp_Dir& theDirection); + //! Initialize evaluator by curve adaptor + Standard_EXPORT GeomEvaluator_OffsetCurve( + const Handle(GeomAdaptor_HCurve)& theBase, + const Standard_Real theOffset, + const gp_Dir& theDirection); + + //! Change the offset value + void SetOffsetValue(Standard_Real theOffset) + { myOffset = theOffset; } + + void SetOffsetDirection(const gp_Dir& theDirection) + { myOffsetDir = theDirection; } + + //! Value of curve + Standard_EXPORT void D0(const Standard_Real theU, + gp_Pnt& theValue) const Standard_OVERRIDE; + //! Value and first derivatives of curve + Standard_EXPORT void D1(const Standard_Real theU, + gp_Pnt& theValue, gp_Vec& theD1) const Standard_OVERRIDE; + //! Value, first and second derivatives of curve + Standard_EXPORT void D2(const Standard_Real theU, + gp_Pnt& theValue, gp_Vec& theD1, gp_Vec& theD2) const Standard_OVERRIDE; + //! Value, first, second and third derivatives of curve + Standard_EXPORT void D3(const Standard_Real theU, + gp_Pnt& theValue, gp_Vec& theD1, + gp_Vec& theD2, gp_Vec& theD3) const Standard_OVERRIDE; + //! Calculates N-th derivatives of curve, where N = theDeriv. Raises if N < 1 + Standard_EXPORT gp_Vec DN(const Standard_Real theU, + const Standard_Integer theDeriv) const Standard_OVERRIDE; + + DEFINE_STANDARD_RTTI(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_HCurve) myBaseAdaptor; + + Standard_Real myOffset; ///< offset value + gp_Dir myOffsetDir; ///< offset direction +}; + +DEFINE_STANDARD_HANDLE(GeomEvaluator_OffsetCurve, GeomEvaluator_Curve) + + +#endif // _GeomEvaluator_OffsetCurve_HeaderFile diff --git a/src/TKG2d/PACKAGES b/src/TKG2d/PACKAGES index 0f938b9e0b..617a12fae2 100755 --- a/src/TKG2d/PACKAGES +++ b/src/TKG2d/PACKAGES @@ -5,3 +5,4 @@ Adaptor2d Geom2dLProp Geom2dAdaptor GProp +Geom2dEvaluator