From 170e4c0bf59f1bb7fdebb594d1624222cb8c6623 Mon Sep 17 00:00:00 2001 From: Artem Zhidkov Date: Thu, 19 Mar 2020 14:32:02 +0300 Subject: [PATCH] Fix the wrong calculation of B-spline value during sketch solving (issue #17347) --- src/GeomAPI/GeomAPI_BSpline2d.cpp | 20 ++++++++- .../PlaneGCSSolver_GeoExtensions.cpp | 43 ++++++++++++------- .../PlaneGCSSolver/PlaneGCSSolver_Tools.cpp | 32 ++++++++++++++ .../PlaneGCSSolver/PlaneGCSSolver_Tools.h | 4 ++ .../SketchSolver_ConstraintCoincidence.cpp | 8 ++++ 5 files changed, 90 insertions(+), 17 deletions(-) diff --git a/src/GeomAPI/GeomAPI_BSpline2d.cpp b/src/GeomAPI/GeomAPI_BSpline2d.cpp index f702a278a..5835d9707 100644 --- a/src/GeomAPI/GeomAPI_BSpline2d.cpp +++ b/src/GeomAPI/GeomAPI_BSpline2d.cpp @@ -179,8 +179,24 @@ const bool GeomAPI_BSpline2d::parameter(const std::shared_ptr the const double theTolerance, double& theParameter) const { - return GeomLib_Tool::Parameter(MY_BSPLINE, thePoint->impl(), - theTolerance, theParameter) == Standard_True; + const gp_Pnt2d& aPoint = thePoint->impl(); + bool isOk = GeomLib_Tool::Parameter(MY_BSPLINE, aPoint, + theTolerance, theParameter) == Standard_True; + if (!isOk) { + // Sometimes OCCT's Extrema algorithm cannot find the parameter on B-spline curve + // (usually, if the point is near the curve extremity). + // Workaround: compute distance to each boundary point + isOk = true; + double aDistPS = aPoint.Distance(MY_BSPLINE->Poles().First()); + double aDistPE = aPoint.Distance(MY_BSPLINE->Poles().Last()); + if (aDistPS < aDistPE && aDistPS < theTolerance) + theParameter = MY_BSPLINE->Knots().First(); + else if (aDistPE < aDistPS && aDistPE < theTolerance) + theParameter = MY_BSPLINE->Knots().Last(); + else + isOk = false; + } + return isOk; } void GeomAPI_BSpline2d::D0(const double theU, std::shared_ptr& thePoint) diff --git a/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_GeoExtensions.cpp b/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_GeoExtensions.cpp index 480dbd58d..f69cd48b5 100644 --- a/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_GeoExtensions.cpp +++ b/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_GeoExtensions.cpp @@ -38,7 +38,32 @@ DeriVector2 BSplineImpl::Value(double u, double du, double* derivparam) std::shared_ptr aDeriv; myCurve->D1(u, aValue, aDeriv); - return DeriVector2(aValue->x(), aValue->y(), aDeriv->x() * du, aDeriv->y() * du); + // calculate the derivative on solver's parameter + std::shared_ptr aValueDeriv(new GeomAPI_Pnt2d(0.0, 0.0)); + bool hasParam = false; + std::list aPolesDeriv; + for (GCS::VEC_P::iterator anIt = poles.begin(); anIt != poles.end(); ++anIt) { + double x = 0.0, y = 0.0; + if (anIt->x == derivparam) { + x = 1.0; + hasParam = true; + } + else if (anIt->y == derivparam) { + y = 1.0; + hasParam = true; + } + aPolesDeriv.push_back(GeomPnt2dPtr(new GeomAPI_Pnt2d(x, y))); + } + if (hasParam) { + // use non-periodic curve, because the most of B-spline coefficients are 0, + // thus, it is not necessary to keep original knots and multiplicities to get correct value + std::shared_ptr aCurveDeriv( + new GeomAPI_BSpline2d(degree, aPolesDeriv, myCachedWeights)); + aCurveDeriv->D0(u, aValueDeriv); + } + + return DeriVector2(aValue->x(), aValue->y(), + aValueDeriv->x() + aDeriv->x() * du, aValueDeriv->y() + aDeriv->y() * du); } DeriVector2 BSplineImpl::CalculateNormal(Point &p, double* derivparam) @@ -47,20 +72,8 @@ DeriVector2 BSplineImpl::CalculateNormal(Point &p, double* derivparam) rebuildCache(); double u = 0.0; - if (!myCurve->parameter(GeomPnt2dPtr(new GeomAPI_Pnt2d(*p.x, *p.y)), 1e100, u)) { - // Sometimes OCCT's Extrema algorithm cannot find the parameter on B-spline curve - // (usually, if the point is near the curve extremity). - // Workaround: compute distance to each boundary point - double aDistPS = PlaneGCSSolver_Tools::distance(p, poles.front()); - double aDistPE = PlaneGCSSolver_Tools::distance(p, poles.back()); - static const double THE_TOLERANCE = 1.e-6; - if (aDistPS < aDistPE && aDistPS < THE_TOLERANCE) - u = *knots.front(); - else if (aDistPE < aDistPS && aDistPE < THE_TOLERANCE) - u = *knots.back(); - else - return DeriVector2(); - } + if (!myCurve->parameter(GeomPnt2dPtr(new GeomAPI_Pnt2d(*p.x, *p.y)), 1e100, u)) + return DeriVector2(); std::shared_ptr aValue; std::shared_ptr aDeriv; diff --git a/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_Tools.cpp b/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_Tools.cpp index e29d23018..5cc2e86c1 100644 --- a/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_Tools.cpp +++ b/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_Tools.cpp @@ -65,6 +65,7 @@ #include #include +#include #include #include #include @@ -362,6 +363,37 @@ std::shared_ptr PlaneGCSSolver_Tools::ellipse(EntityWrapperPt new GeomAPI_Ellipse2d(aCenter, anAxis, anEllipse->getRadMaj(), *anEllipse->radmin)); } +std::shared_ptr PlaneGCSSolver_Tools::bspline(EntityWrapperPtr theEntity) +{ + if (theEntity->type() != ENTITY_BSPLINE) + return std::shared_ptr(); + + std::shared_ptr anEntity = + std::dynamic_pointer_cast(theEntity); + std::shared_ptr aSpline = + std::dynamic_pointer_cast(anEntity->entity()); + + std::list aPoles; + for (GCS::VEC_P::iterator anIt = aSpline->poles.begin(); anIt != aSpline->poles.end(); ++anIt) + aPoles.push_back(GeomPnt2dPtr(new GeomAPI_Pnt2d(*anIt->x, *anIt->y))); + + std::list aWeights; + for (GCS::VEC_pD::iterator anIt = aSpline->weights.begin(); + anIt != aSpline->weights.end(); ++anIt) + aWeights.push_back(**anIt); + + std::list aKnots; + for (GCS::VEC_pD::iterator anIt = aSpline->knots.begin(); anIt != aSpline->knots.end(); ++anIt) + aKnots.push_back(**anIt); + + std::list aMultiplicities; + aMultiplicities.assign(aSpline->mult.begin(), aSpline->mult.end()); + + return std::shared_ptr( + new GeomAPI_BSpline2d(aSpline->degree, aPoles, aWeights, + aKnots, aMultiplicities, aSpline->periodic)); +} + void PlaneGCSSolver_Tools::recalculateArcParameters(EntityWrapperPtr theArc) { std::shared_ptr anEdge = diff --git a/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_Tools.h b/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_Tools.h index 2604d1541..aba6c1814 100644 --- a/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_Tools.h +++ b/src/SketchSolver/PlaneGCSSolver/PlaneGCSSolver_Tools.h @@ -24,6 +24,7 @@ #include #include +class GeomAPI_BSpline2d; class GeomAPI_Circ2d; class GeomAPI_Ellipse2d; class GeomAPI_Lin2d; @@ -88,6 +89,9 @@ namespace PlaneGCSSolver_Tools /// \brief Convert entity to ellipse /// \return empty pointer if the entity is not an ellipse std::shared_ptr ellipse(EntityWrapperPtr theEntity); + /// \brief Convert entity to Bs-pline + /// \return empty pointer if the entity is not an ellipse + std::shared_ptr bspline(EntityWrapperPtr theEntity); /// \brief Convert entity to line /// \return empty pointer if the entity is not a line diff --git a/src/SketchSolver/SketchSolver_ConstraintCoincidence.cpp b/src/SketchSolver/SketchSolver_ConstraintCoincidence.cpp index c00dcfce8..cecd7a616 100644 --- a/src/SketchSolver/SketchSolver_ConstraintCoincidence.cpp +++ b/src/SketchSolver/SketchSolver_ConstraintCoincidence.cpp @@ -24,6 +24,9 @@ #include #include +#include +#include + #include #include @@ -270,6 +273,11 @@ void SketchSolver_ConstraintCoincidence::getAttributes( std::shared_ptr aStorage = std::dynamic_pointer_cast(myStorage); myAuxValue.reset(new PlaneGCSSolver_ScalarWrapper(aStorage->createParameter())); + // calculate the parameter of the point on B-spline nearest to the constrained point. + GeomPnt2dPtr aPoint = PlaneGCSSolver_Tools::point(theAttributes[0]); + std::shared_ptr aSpline = PlaneGCSSolver_Tools::bspline(theAttributes[2]); + if (aPoint && aSpline) + aSpline->parameter(aPoint, 1.e100, *myAuxValue->scalar()); } else { // obtain extremity points of the coincident feature for further checking of multi-coincidence -- 2.39.2