Sketcher: solver: Value constraint

(morphed from hyperbola arc range constraint)
(compiles, but untested)

src/Mod/Sketcher/App/planegcs/Constraints.cpp

@@ -1424,81 +1424,83 @@ double ConstraintEllipticalArcRangeToEndPoints::maxStep(MAP_pD_D &dir, double li
 }
 
 // HyperbolicArcRangeToEndPoints
-ConstraintHyperbolicArcRangeToEndPoints::ConstraintHyperbolicArcRangeToEndPoints(Point &p, ArcOfHyperbola &a, double *angle_t)
+ConstraintCurveValue::ConstraintCurveValue(Point &p, double* pcoord, Curve& crv, double *u)
 {
     pvec.push_back(p.x);
     pvec.push_back(p.y);
-    pvec.push_back(angle_t);
-    e = a;
-    e.PushOwnParams(pvec);
+    pvec.push_back(pcoord);
+    pvec.push_back(u);
+    crv.PushOwnParams(pvec);
+    this->crv = crv.Copy();
+    pvecChangedFlag = true;
     origpvec = pvec;
     rescale();
 }
 
-void ConstraintHyperbolicArcRangeToEndPoints::ReconstructGeomPointers()
+ConstraintCurveValue::~ConstraintCurveValue()
+{
+    delete this->crv; this->crv = 0;
+}
+
+void ConstraintCurveValue::ReconstructGeomPointers()
 {
     int i=0;
     p.x=pvec[i]; i++;
     p.y=pvec[i]; i++;
-    i++;//we have an inline function for the angle
-    e.ReconstructOnNewPvec(pvec, i);
+    i++;//we have an inline function for point coordinate
+    i++;//we have an inline function for the parameterU
+    this->crv->ReconstructOnNewPvec(pvec, i);
     pvecChangedFlag = false;
 }
 
-ConstraintType ConstraintHyperbolicArcRangeToEndPoints::getTypeId()
+ConstraintType ConstraintCurveValue::getTypeId()
 {
     return HyperbolicArcRangeToEndPoints;
 }
 
-void ConstraintHyperbolicArcRangeToEndPoints::rescale(double coef)
+void ConstraintCurveValue::rescale(double coef)
 {
     scale = coef * 1;
 }
 
-void ConstraintHyperbolicArcRangeToEndPoints::errorgrad(double *err, double *grad, double *param)
+void ConstraintCurveValue::errorgrad(double *err, double *grad, double *param)
 {
     if (pvecChangedFlag) ReconstructGeomPointers();
-    
-    DeriVector2 c(e.center, param);
-    DeriVector2 f1(e.focus1, param);
-    DeriVector2 emaj = f1.subtr(c).getNormalized();
-    DeriVector2 emin = emaj.rotate90ccw();
-    double b, db;
-    b = *e.radmin; db = e.radmin==param ? 1.0 : 0.0;
-    double a, da;
-    a = e.getRadMaj(c,f1,b,db,da);
-    
-    DeriVector2 multimaj = emaj.multD(b, db);//a vector to muptiply pc by to yield an x for atan2. This is a minor radius drawn along major axis.
-    DeriVector2 multimin = emin.multD(a, da);//to yield y for atan2
-    
-    DeriVector2 pv(p, param);
-    DeriVector2 pc = pv.subtr(c); //point referenced to ellipse's center
-    
-    double x, dx, y, dy;//distorted coordinates of point in ellipse's coordinates, to be fed to atan2 to yiels a t-parameter (called "angle" here)
-    x = pc.scalarProd(multimaj, &dx);
-    y = pc.scalarProd(multimin, &dy);
-    double xylen2 = x*x + y*y ;//square of length of (x,y)
-    
-    double si, co;
-    si = sin(*angle()); co = cos(*angle());
-    
-    double dAngle = param==angle() ? 1.0 : 0.0;
-    
-    if (err)
-        *err = atan2(-si*x+co*y, co*x+si*y);//instead of calculating atan2(y,x) and subtracting angle, we rotate (x,y) by -angle and calculate atan2 of the result. Hopefully, this will not force angles to zero when x,y happen to be zero. Plus, one atan2 is cheaper to compute than two atan2's.
+
+    double u, du;
+    u = *(this->u()); du = ( param == this->u() )   ?   1.0 : 0.0;
+
+    DeriVector2 P_to; //point of curve at parameter value of u, in global coordinates
+    P_to = this->crv->Value(u,du,param);
+
+    DeriVector2 P_from(this->p, param); //point to be constrained
+
+    DeriVector2 err_vec = P_from.subtr(P_to);
+
+    if (this->pcoord() == this->p.x){ //this constraint is for X projection
+        if (err)
+            *err = err_vec.x;
         if (grad)
-            *grad = -dAngle + ( -dx*y / xylen2  +  dy*x / xylen2 );
-        
+            *grad = err_vec.dx;
+    } else if (this->pcoord() == this->p.y) {//this constraint is for Y projection
+        if (err)
+            *err = err_vec.y;
+        if (grad)
+            *grad = err_vec.dy;
+    } else {
+        assert(false/*this constraint is neighter X nor Y. Nothing to do..*/);
+    }
+
 }
 
-double ConstraintHyperbolicArcRangeToEndPoints::error()
+double ConstraintCurveValue::error()
 {
     double err;
     errorgrad(&err,0,0);
     return scale * err;
 }
 
-double ConstraintHyperbolicArcRangeToEndPoints::grad(double *param)
+double ConstraintCurveValue::grad(double *param)
 {
     //first of all, check that we need to compute anything.
     if ( findParamInPvec(param) == -1  ) return 0.0;
@@ -1509,15 +1511,17 @@ double ConstraintHyperbolicArcRangeToEndPoints::grad(double *param)
     return deriv*scale;
 }    
 
-double ConstraintHyperbolicArcRangeToEndPoints::maxStep(MAP_pD_D &dir, double lim)
+double ConstraintCurveValue::maxStep(MAP_pD_D &dir, double lim)
 {
     // step(angle()) <= pi/18 = 10°
-    MAP_pD_D::iterator it = dir.find(angle());
+    /* TODO: curve-dependent parameter change limiting??
+    MAP_pD_D::iterator it = dir.find(this->u());
     if (it != dir.end()) {
         double step = std::abs(it->second);
         if (step > M_PI/18.)
             lim = std::min(lim, (M_PI/18.) / step);
     }
+    */
     return lim;
 }
 

src/Mod/Sketcher/App/planegcs/Constraints.h

@@ -456,16 +456,25 @@ namespace GCS
         virtual double maxStep(MAP_pD_D &dir, double lim=1.);
     };
     
-    class ConstraintHyperbolicArcRangeToEndPoints : public Constraint
+    class ConstraintCurveValue : public Constraint
     {
     private:
-        inline double* angle() { return pvec[2]; }
+        inline double* pcoord() { return pvec[2]; } //defines, which coordinate of point is being constrained by this constraint
+        inline double* u() { return pvec[3]; }
         void errorgrad(double* err, double* grad, double *param); //error and gradient combined. Values are returned through pointers.
         void ReconstructGeomPointers(); //writes pointers in pvec to the parameters of crv1, crv2 and poa
-        Hyperbola e;
+        Curve* crv;
         Point p;
     public:
-        ConstraintHyperbolicArcRangeToEndPoints(Point &p, ArcOfHyperbola &a, double *angle_t);
+        /**
+         * @brief ConstraintCurveValue: solver constraint that ties parameter value with point coordinates, according to curve's parametric equation.
+         * @param p : endpoint to be constrained
+         * @param pcoord : pointer to point coordinate to be constrained. Must be either p.x or p.y
+         * @param crv : the curve (crv->Value() must be functional)
+         * @param u : pointer to u parameter corresponding to the point
+         */
+        ConstraintCurveValue(Point &p, double* pcoord, Curve& crv, double* u);
+        ~ConstraintCurveValue();
         virtual ConstraintType getTypeId();
         virtual void rescale(double coef=1.);
         virtual double error();

src/Mod/Sketcher/App/planegcs/GCS.cpp

@@ -603,7 +603,10 @@ int System::addConstraintArcOfEllipseRules(ArcOfEllipse &a, int tagId)
 
 int System::addConstraintHyperbolicArcRangeToEndPoints(Point &p, ArcOfHyperbola &a, double *angle, int tagId)
 {
-    Constraint *constr = new ConstraintHyperbolicArcRangeToEndPoints(p,a,angle);
+    Constraint *constr = new ConstraintCurveValue(p,p.x,a,angle);
+    constr->setTag(tagId);
+    addConstraint(constr);
+    constr = new ConstraintCurveValue(p,p.y,a,angle);
     constr->setTag(tagId);
     return addConstraint(constr);
 }
@@ -612,10 +615,7 @@ int System::addConstraintHyperbolicArcRangeToEndPoints(Point &p, ArcOfHyperbola
 int System::addConstraintArcOfHyperbolaRules(ArcOfHyperbola &a, int tagId)
 {
     addConstraintHyperbolicArcRangeToEndPoints(a.start,a,a.startAngle, tagId);
-    addConstraintHyperbolicArcRangeToEndPoints(a.end,a,a.endAngle, tagId);
-    
-    addConstraintPointOnHyperbolicArc(a.start, a, tagId);
-    return addConstraintPointOnHyperbolicArc(a.end, a, tagId);
+    return addConstraintHyperbolicArcRangeToEndPoints(a.end,a,a.endAngle, tagId);
 }
 
 int System::addConstraintPointOnArc(Point &p, Arc &a, int tagId)