From dd64ab8ea08c9ddf6f74e421cde5e4dabbfe1d1f Mon Sep 17 00:00:00 2001 From: Ajinkya Dahale Date: Thu, 24 Nov 2022 21:07:18 +0530 Subject: [PATCH] [Sketcher][planegcs] Make changes as per review on #7484 Similar to 2715a66ff02a46f94ae3fc6527fd446e666b8e58. Added some comments and removed some TODO's. Return grad values directly rather than doing summation. --- src/Mod/Sketcher/App/planegcs/Constraints.cpp | 57 ++++++++++++------- 1 file changed, 38 insertions(+), 19 deletions(-) diff --git a/src/Mod/Sketcher/App/planegcs/Constraints.cpp b/src/Mod/Sketcher/App/planegcs/Constraints.cpp index 262fa9e778..a583f3bc19 100644 --- a/src/Mod/Sketcher/App/planegcs/Constraints.cpp +++ b/src/Mod/Sketcher/App/planegcs/Constraints.cpp @@ -260,8 +260,9 @@ ConstraintSlopeAtBSplineKnot::ConstraintSlopeAtBSplineKnot(BSpline& b, Line& l, pvec.reserve(3*numpoles + 4); + // `startpole` is the first pole affecting the knot with `knotindex` size_t startpole = 0; - // TODO: Adjust for periodic knot + // See `System::addConstraintInternalAlignmentKnotPoint()` for some elaboration for (size_t j = 1; j <= knotindex; ++j) startpole += b.mult[j]; if (!b.periodic && startpole >= b.poles.size()) @@ -278,7 +279,7 @@ ConstraintSlopeAtBSplineKnot::ConstraintSlopeAtBSplineKnot(BSpline& b, Line& l, pvec.push_back(l.p2.x); pvec.push_back(l.p2.y); - // TODO: Set up factors to get slope at knot point + // Set up factors to get slope at knot point std::vector tempfactors((numpoles + 1), 1.0 / (numpoles + 1)); factors.resize(numpoles); slopefactors.resize(numpoles); @@ -316,7 +317,6 @@ void ConstraintSlopeAtBSplineKnot::rescale(double coef) double ConstraintSlopeAtBSplineKnot::error() { - // TODO: Fill up double xsum = 0., xslopesum = 0.; double ysum = 0., yslopesum = 0.; double wsum = 0., wslopesum = 0.; @@ -333,9 +333,13 @@ double ConstraintSlopeAtBSplineKnot::error() } // This is actually wsum^2 * the respective slopes + // See Eq (19) from: + // https://forum.freecadweb.org/viewtopic.php?f=9&t=71130&start=120#p635538 double slopex = wsum*xslopesum - wslopesum*xsum; double slopey = wsum*yslopesum - wslopesum*ysum; + // Normalizing it ensures that the cross product is not zero just because + // one vector is zero. double linex = *linep2x() - *linep1x(); double liney = *linep2y() - *linep1y(); double dirx = linex / sqrt(linex*linex + liney*liney); @@ -347,15 +351,17 @@ double ConstraintSlopeAtBSplineKnot::error() double ConstraintSlopeAtBSplineKnot::grad(double *param) { + // Equations are from here: + // https://forum.freecadweb.org/viewtopic.php?f=9&t=71130&start=120#p635538 double result = 0.0; double linex = *linep2x() - *linep1x(); double liney = *linep2y() - *linep1y(); double dirx = linex / sqrt(linex*linex + liney*liney); double diry = liney / sqrt(linex*linex + liney*liney); - // TODO: Fill up for (size_t i = 0; i < numpoles; ++i) { if (param == polexat(i)) { + // Eq. (21) double wsum = 0., wslopesum = 0.; for (size_t j = 0; j < numpoles; ++j) { double wcontrib = *weightat(j) * factors[j]; @@ -363,9 +369,11 @@ double ConstraintSlopeAtBSplineKnot::grad(double *param) wsum += wcontrib; wslopesum += wslopecontrib; } - result += (wsum*slopefactors[i] - wslopesum*factors[i]) * diry; + result = (wsum*slopefactors[i] - wslopesum*factors[i]) * diry; + return scale * result; } if (param == poleyat(i)) { + // Eq. (21) double wsum = 0., wslopesum = 0.; for (size_t i = 0; i < numpoles; ++i) { double wcontrib = *weightat(i) * factors[i]; @@ -373,9 +381,11 @@ double ConstraintSlopeAtBSplineKnot::grad(double *param) wsum += wcontrib; wslopesum += wslopecontrib; } - result += - (wsum*slopefactors[i] - wslopesum*factors[i]) * dirx; + result = - (wsum*slopefactors[i] - wslopesum*factors[i]) * dirx; + return scale * result; } if (param == weightat(i)) { + // Eq. (22) double xsum = 0., xslopesum = 0.; double ysum = 0., yslopesum = 0.; for (size_t j = 0; j < numpoles; ++j) { @@ -386,9 +396,10 @@ double ConstraintSlopeAtBSplineKnot::grad(double *param) ysum += wcontrib * (*poleyat(j) - *poleyat(i)); yslopesum += wslopecontrib * (*poleyat(j) - *poleyat(i)); } - result += + result = (factors[i]*xslopesum - slopefactors[i]*xsum) * diry - (factors[i]*yslopesum - slopefactors[i]*ysum) * dirx; + return scale * result; } } @@ -417,27 +428,35 @@ double ConstraintSlopeAtBSplineKnot::grad(double *param) if (param == linep1x()) { getSlopes(); - double dDirxDlinex = (liney*liney) / pow(linex*linex + liney*liney, 1.5); - double dDiryDlinex = -(linex*liney) / pow(linex*linex + liney*liney, 1.5); - result += slopex*(-dDiryDlinex) - slopey*(-dDirxDlinex); + double dDirxDLinex = (liney*liney) / pow(linex*linex + liney*liney, 1.5); + double dDiryDLinex = -(linex*liney) / pow(linex*linex + liney*liney, 1.5); + // NOTE: d(linex)/d(x1) = -1 + result = slopex*(-dDiryDLinex) - slopey*(-dDirxDLinex); + return scale * result; } if (param == linep2x()) { getSlopes(); - double dDirxDlinex = (liney*liney) / pow(linex*linex + liney*liney, 1.5); - double dDiryDlinex = -(linex*liney) / pow(linex*linex + liney*liney, 1.5); - result += slopex*dDiryDlinex - slopey*dDirxDlinex; + double dDirxDLinex = (liney*liney) / pow(linex*linex + liney*liney, 1.5); + double dDiryDLinex = -(linex*liney) / pow(linex*linex + liney*liney, 1.5); + // NOTE: d(linex)/d(x2) = 1 + result = slopex*dDiryDLinex - slopey*dDirxDLinex; + return scale * result; } if (param == linep1y()) { getSlopes(); - double dDirxDliney = -(linex*liney) / pow(linex*linex + liney*liney, 1.5); - double dDiryDliney = (linex*linex) / pow(linex*linex + liney*liney, 1.5); - result += slopex*(-dDiryDliney) - slopey*(-dDirxDliney); + double dDirxDLiney = -(linex*liney) / pow(linex*linex + liney*liney, 1.5); + double dDiryDLiney = (linex*linex) / pow(linex*linex + liney*liney, 1.5); + // NOTE: d(liney)/d(y1) = -1 + result = slopex*(-dDiryDLiney) - slopey*(-dDirxDLiney); + return scale * result; } if (param == linep2y()) { getSlopes(); - double dDirxDliney = -(linex*liney) / pow(linex*linex + liney*liney, 1.5); - double dDiryDliney = (linex*linex) / pow(linex*linex + liney*liney, 1.5); - result += slopex*dDiryDliney - slopey*dDirxDliney; + double dDirxDLiney = -(linex*liney) / pow(linex*linex + liney*liney, 1.5); + double dDiryDLiney = (linex*linex) / pow(linex*linex + liney*liney, 1.5); + // NOTE: d(liney)/d(y2) = 1 + result = slopex*dDiryDLiney - slopey*dDirxDLiney; + return scale * result; } return scale * result;