Skip to content

Commit

Permalink
GCS: New Constraint for equal size lines, avoiding to use an extra pa…
Browse files Browse the repository at this point in the history 
…rameter

============================================================================

This specific constraint removes the free parameter of the previous implementation. This solves:
https://tracker.freecadweb.org/view.php?id=4501

fixes #4501

However, this implementation of equal size produces zero gradients when coordinates of lines are aligned,
e.g. vertical or horizontal. These zero gradients, which are mathematically right ruin the diagnosis, which
regards corresponding elements as fully constraint (because they are locked from a solver point of view), when
they are simply locked, but are movable and constrainable. For this, when the rightful gradient is small enough
(<1e-10) it is substituted by a surrogate gradient of 1e-10, which solves the problem with the diagnose, which
treats as zero only values under 1e-13 (pivot threshold used in QR decomposition).

This special behaviour fixes the wrong detection here:
https://forum.freecadweb.org/viewtopic.php?f=8&t=53466&start=40#p464168

It also fixes this one:
https://forum.freecadweb.org/viewtopic.php?p=468585#p468587
  • Loading branch information
@abdullahtahiriyo
abdullahtahiriyo committed Jan 23, 2021
1 parent ea49ed9 commit 5609269
Show file tree
Hide file tree
Showing 2 changed files with 160 additions and 45 deletions.
167 changes: 133 additions & 34 deletions src/Mod/Sketcher/App/planegcs/Constraints.cpp
View file
Expand Up @@ -533,18 +533,18 @@ void ConstraintPointOnPerpBisector::errorgrad(double *err, double *grad, double
DeriVector2 p0(Point(p0x(),p0y()), param);
DeriVector2 p1(Point(p1x(),p1y()), param);
DeriVector2 p2(Point(p2x(),p2y()), param);

DeriVector2 d1 = p0.subtr(p1);
DeriVector2 d2 = p0.subtr(p2);
DeriVector2 D = p2.subtr(p1).getNormalized();

double projd1, dprojd1;
projd1 = d1.scalarProd(D, &dprojd1);

double projd2, dprojd2;
projd2 = d2.scalarProd(D, &dprojd2);
if (err)

if (err)
*err = projd1+projd2;
if (grad)
*grad = dprojd1+dprojd2;
Expand Down Expand Up @@ -964,15 +964,15 @@ void ConstraintPointOnEllipse::rescale(double coef)
}

double ConstraintPointOnEllipse::error()
{
{
double X_0 = *p1x();
double Y_0 = *p1y();
double X_c = *cx();
double Y_c = *cy();
double Y_c = *cy();
double X_F1 = *f1x();
double Y_F1 = *f1y();
double b = *rmin();

double err=sqrt(pow(X_0 - X_F1, 2) + pow(Y_0 - Y_F1, 2)) + sqrt(pow(X_0 +
X_F1 - 2*X_c, 2) + pow(Y_0 + Y_F1 - 2*Y_c, 2)) - 2*sqrt(pow(b, 2) +
pow(X_F1 - X_c, 2) + pow(Y_F1 - Y_c, 2));
Expand Down Expand Up @@ -1100,7 +1100,7 @@ double ConstraintEllipseTangentLine::error()
}

double ConstraintEllipseTangentLine::grad(double *param)
{
{
//first of all, check that we need to compute anything.
if ( findParamInPvec(param) == -1 ) return 0.0;

Expand Down Expand Up @@ -1218,15 +1218,15 @@ void ConstraintInternalAlignmentPoint2Ellipse::errorgrad(double *err, double *gr
}

double ConstraintInternalAlignmentPoint2Ellipse::error()
{
{
double err;
errorgrad(&err,0,0);
return scale * err;

}

double ConstraintInternalAlignmentPoint2Ellipse::grad(double *param)
{
{
//first of all, check that we need to compute anything.
if ( findParamInPvec(param) == -1 ) return 0.0;

Expand Down Expand Up @@ -1318,7 +1318,7 @@ void ConstraintInternalAlignmentPoint2Hyperbola::errorgrad(double *err, double *
case HyperbolaNegativeMajorY:
poa = c.sum(emaj.multD(-a, -da));
by_y_not_by_x = AlignmentType == HyperbolaNegativeMajorY;
break;
break;
case HyperbolaPositiveMinorX:
case HyperbolaPositiveMinorY:
{
Expand Down Expand Up @@ -1414,7 +1414,7 @@ void ConstraintEqualMajorAxesConic::errorgrad(double *err, double *grad, double
}

double ConstraintEqualMajorAxesConic::error()
{
{
double err;
errorgrad(&err,0,0);
return scale * err;
Expand Down Expand Up @@ -1490,7 +1490,7 @@ void ConstraintEqualFocalDistance::errorgrad(double *err, double *grad, double *
}

double ConstraintEqualFocalDistance::error()
{
{
double err;
errorgrad(&err,0,0);
return scale * err;
Expand Down Expand Up @@ -1588,12 +1588,12 @@ double ConstraintCurveValue::grad(double *param)
{
//first of all, check that we need to compute anything.
if ( findParamInPvec(param) == -1 ) return 0.0;

double deriv;
errorgrad(0, &deriv, param);

return deriv*scale;
}
}

double ConstraintCurveValue::maxStep(MAP_pD_D &/*dir*/, double lim)
{
Expand Down Expand Up @@ -1647,15 +1647,15 @@ void ConstraintPointOnHyperbola::rescale(double coef)
}

double ConstraintPointOnHyperbola::error()
{
{
double X_0 = *p1x();
double Y_0 = *p1y();
double X_c = *cx();
double Y_c = *cy();
double Y_c = *cy();
double X_F1 = *f1x();
double Y_F1 = *f1y();
double b = *rmin();

// Full sage worksheet at:
// http://forum.freecadweb.org/viewtopic.php?f=10&t=8038&p=110447#p110447
//
Expand All @@ -1682,15 +1682,15 @@ double ConstraintPointOnHyperbola::grad(double *param)
param == f1x() || param == f1y() ||
param == cx() || param == cy() ||
param == rmin()) {

double X_0 = *p1x();
double Y_0 = *p1y();
double X_c = *cx();
double Y_c = *cy();
double X_F1 = *f1x();
double Y_F1 = *f1y();
double b = *rmin();

if (param == p1x())
deriv += -(X_0 - X_F1)/sqrt(pow(X_0 - X_F1, 2) + pow(Y_0 - Y_F1, 2)) +
(X_0 + X_F1 - 2*X_c)/sqrt(pow(X_0 + X_F1 - 2*X_c, 2) + pow(Y_0 + Y_F1 -
Expand Down Expand Up @@ -1778,26 +1778,26 @@ void ConstraintPointOnParabola::errorgrad(double *err, double *grad, double *par
DeriVector2 vertex(this->parab->vertex, param);

DeriVector2 point(this->p, param); //point to be constrained to parabola

DeriVector2 focalvect = focus.subtr(vertex);

DeriVector2 xdir = focalvect.getNormalized();

DeriVector2 point_to_focus = point.subtr(focus);

double focal, dfocal;

focal = focalvect.length(dfocal);

double pf, dpf;

pf = point_to_focus.length(dpf);

double proj, dproj;

proj = point_to_focus.scalarProd(xdir, &dproj);
if (err)

if (err)
*err = pf - 2*focal - proj;
if (grad)
*grad = dpf - 2*dfocal - dproj;
Expand All @@ -1815,10 +1815,10 @@ double ConstraintPointOnParabola::grad(double *param)
{
//first of all, check that we need to compute anything.
if ( findParamInPvec(param) == -1 ) return 0.0;

double deriv;
errorgrad(0, &deriv, param);

return deriv*scale;
}

Expand Down Expand Up @@ -2026,5 +2026,104 @@ double ConstraintSnell::grad(double *param)
return scale * deriv;
}

// ConstraintEqualLineLength
ConstraintEqualLineLength::ConstraintEqualLineLength(Line &l1, Line &l2)
{
this->l1 = l1;
this->l1.PushOwnParams(pvec);

this->l2 = l2;
this->l2.PushOwnParams(pvec);
origpvec = pvec;
pvecChangedFlag = true;
rescale();
}

void ConstraintEqualLineLength::ReconstructGeomPointers()
{
int i=0;
l1.ReconstructOnNewPvec(pvec, i);
l2.ReconstructOnNewPvec(pvec, i);
pvecChangedFlag = false;
}

ConstraintType ConstraintEqualLineLength::getTypeId()
{
return EqualLineLength;
}

void ConstraintEqualLineLength::rescale(double coef)
{
scale = coef * 1;
}

void ConstraintEqualLineLength::errorgrad(double *err, double *grad, double *param)
{
if (pvecChangedFlag) ReconstructGeomPointers();

DeriVector2 p1 (l1.p1, param);
DeriVector2 p2 (l1.p2, param);
DeriVector2 p3 (l2.p1, param);
DeriVector2 p4 (l2.p2, param);

DeriVector2 v1 = p1.subtr(p2);
DeriVector2 v2 = p3.subtr(p4);

double length1, dlength1;
length1 = v1.length(dlength1);

double length2, dlength2;
length2 = v2.length(dlength2);

if (err)
*err = length2 - length1;

if (grad) {
*grad = dlength2 - dlength1;
// if the one of the lines gets vertical or horizontal, the gradients will become zero. this will
// affect the diagnose function and the detection of dependent/independent parameters.
//
// So here we maintain the very small derivative of 1e-10 when the gradient is under such value, such
// that the diagnose function with pivot threshold of 1e-13 treats the value as non-zero and correctly
// detects and can tell appart when a parameter is fully constrained or just locked into a maximum/minimum
if(fabs(*grad) < 1e-10) {
double surrogate = 1e-10;
if( param == l1.p1.x )
*grad = v1.x > 0 ? surrogate : -surrogate;
if( param == l1.p1.y )
*grad = v1.y > 0 ? surrogate : -surrogate;
if( param == l1.p2.x )
*grad = v1.x > 0 ? -surrogate : surrogate;
if( param == l1.p2.y )
*grad = v1.y > 0 ? -surrogate : surrogate;
if( param == l2.p1.x )
*grad = v2.x > 0 ? surrogate : -surrogate;
if( param == l2.p1.y )
*grad = v2.y > 0 ? surrogate : -surrogate;
if( param == l2.p2.x )
*grad = v2.x > 0 ? -surrogate : surrogate;
if( param == l2.p2.y )
*grad = v2.y > 0 ? -surrogate : surrogate;
}
}
}

double ConstraintEqualLineLength::error()
{
double err;
errorgrad(&err,0,0);
return scale * err;
}

double ConstraintEqualLineLength::grad(double *param)
{
if ( findParamInPvec(param) == -1 ) return 0.0;

double deriv;
errorgrad(0, &deriv, param);

return deriv*scale;
}


} //namespace GCS

0 comments on commit 5609269

Please sign in to comment.