[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.
This commit is contained in:
Ajinkya Dahale
2022-11-24 21:07:18 +05:30
committed by Chris Hennes
parent 58b26788be
commit dd64ab8ea0

View File

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