stable evaluation of implizit spline-curve representation

libsrc/csg/revolution.cpp

@@ -156,25 +156,42 @@ namespace netgen
   {
     if(spline_coefficient.Size() == 0)
       spline->GetCoeff(spline_coefficient);
+    if(spline_coefficient_shifted.Size() == 0)
+        spline->GetCoeff(spline_coefficient_shifted, spline->StartPI());
 
     Point<2> p;
     CalcProj(point,p);
 
-    return spline_coefficient(0)*p(0)*p(0) + spline_coefficient(1)*p(1)*p(1)
+    double val = spline_coefficient(0)*p(0)*p(0) + spline_coefficient(1)*p(1)*p(1)
       + spline_coefficient(2)*p(0)*p(1) + spline_coefficient(3)*p(0)
       + spline_coefficient(4)*p(1) + spline_coefficient(5);
+
+    Vec<2> pr = p-spline->StartPI();
+
+
+    // cout << "spline_coefficinet = " << spline_coefficient << endl;
+    // cout << "shifted = " << spline_coefficient_shifted << endl;
+    
+    double val2 = spline_coefficient_shifted(0)*pr(0)*pr(0) + spline_coefficient_shifted(1)*pr(1)*pr(1)
+      + spline_coefficient_shifted(2)*pr(0)*pr(1) + spline_coefficient_shifted(3)*pr(0)
+      + spline_coefficient_shifted(4)*pr(1) + spline_coefficient_shifted(5);
+
+    // cout << "val = " << val << " =?= " << val2 << endl;
+    return val2;
   }
 
   void RevolutionFace :: CalcGradient (const Point<3> & point, Vec<3> & grad) const
   {
     if(spline_coefficient.Size() == 0)
       spline->GetCoeff(spline_coefficient);
+    if(spline_coefficient_shifted.Size() == 0)
+        spline->GetCoeff(spline_coefficient_shifted, spline->StartPI());
 
     Vec<3> point_minus_p0 = point-p0;
 
     Point<2> p;
     CalcProj0(point_minus_p0,p);
-
+    /*
     const double dFdxbar = 2.*spline_coefficient(0)*p(0) + spline_coefficient(2)*p(1) + spline_coefficient(3);
 
     if(fabs(p(1)) > 1e-10)
@@ -193,6 +210,27 @@ namespace netgen
 	grad(2) = dFdxbar*v_axis(2);
 	//(*testout) << "grad2("<<point<<") = " << grad << endl;
       }
+    */
+    Vec<2> pr = p-spline->StartPI();
+    const double dFdxbar = 2.*spline_coefficient_shifted(0)*pr(0) + spline_coefficient_shifted(2)*pr(1) + spline_coefficient_shifted(3);
+
+    if(fabs(p(1)) > 1e-10)
+      {
+	const double dFdybar = 2.*spline_coefficient_shifted(1)*pr(1) + spline_coefficient_shifted(2)*pr(0) + spline_coefficient_shifted(4);
+
+	grad(0) = dFdxbar*v_axis(0) + dFdybar * ( point_minus_p0(0)-v_axis(0)*p(0) )/p(1);
+	grad(1) = dFdxbar*v_axis(1) + dFdybar * ( point_minus_p0(1)-v_axis(1)*p(0) )/p(1);
+	grad(2) = dFdxbar*v_axis(2) + dFdybar * ( point_minus_p0(2)-v_axis(2)*p(0) )/p(1);
+	//(*testout) << "grad1("<<point<<") = " << grad << endl;
+      }
+    else
+      {
+	grad(0) = dFdxbar*v_axis(0);
+	grad(1) = dFdxbar*v_axis(1);
+	grad(2) = dFdxbar*v_axis(2);
+	//(*testout) << "grad2("<<point<<") = " << grad << endl;
+      }
+    
   }
 
   
@@ -559,8 +597,17 @@ namespace netgen
 
     if (!spline -> InConvexHull(p2d, eps)) return false;
 
+    /*
     double val = spline_coefficient(0)*p2d(0)*p2d(0) + spline_coefficient(1)*p2d(1)*p2d(1) + spline_coefficient(2)*p2d(0)*p2d(1) +
       spline_coefficient(3)*p2d(0) + spline_coefficient(4)*p2d(1) + spline_coefficient(5);
+    */
+    Vec<2> pr = p2d - spline->StartPI();
+    double val = spline_coefficient_shifted(0)*pr(0)*pr(0)
+      + spline_coefficient_shifted(1)*pr(1)*pr(1)
+      + spline_coefficient_shifted(2)*pr(0)*pr(1) 
+      + spline_coefficient_shifted(3)*pr(0)
+      + spline_coefficient_shifted(4)*pr(1)
+      + spline_coefficient_shifted(5);
     
     return (fabs(val) < eps);
       /*

libsrc/csg/revolution.hpp

@@ -18,7 +18,9 @@ namespace netgen
 
     int id;
 
+    // coefficient for implicizt polynomial
     mutable Vector spline_coefficient;
+    mutable Vector spline_coefficient_shifted;
 
 
     Array < Vec<2>* > checklines_vec;

libsrc/gprim/spline.cpp

@@ -190,6 +190,51 @@ namespace netgen
   }
 
 
+  template<int D>
+  void SplineSeg3<D> :: GetCoeff (Vector & u, Point<D> pref) const
+  {
+    DenseMatrix a(6, 6);
+    DenseMatrix ata(6, 6);
+    Vector f(6);
+
+    u.SetSize(6);
+
+    //  ata.SetSymmetric(1);
+
+    double t = 0;
+    for (int i = 0; i < 5; i++, t += 0.25)
+      {
+	Vec<D> p = GetPoint (t)-pref;
+	a(i, 0) = p(0) * p(0);
+	a(i, 1) = p(1) * p(1);
+	a(i, 2) = p(0) * p(1);
+	a(i, 3) = p(0);
+	a(i, 4) = p(1);
+	a(i, 5) = 1;
+      }
+    a(5, 0) = 1;
+
+    CalcAtA (a, ata);
+
+    u = 0;
+    u(5) = 1;
+    a.MultTrans (u, f);
+    ata.Solve (f, u);
+
+    // the sign
+    Point<D> p0 = GetPoint(0);
+    Vec<D> ht = GetTangent(0);
+    Vec<2> tang(ht(0), ht(1));
+
+    double gradx = u(3);
+    double grady = u(4);
+    // double gradx = 2.*u(0)*p0(0) + u(2)*p0(1) + u(3);
+    // double grady = 2.*u(1)*p0(1) + u(2)*p0(0) + u(4);
+    Vec<2> gradn (grady, -gradx);
+  
+    if (tang * gradn < 0) u *= -1;
+  }
+  
 
 
   template<int D>

libsrc/gprim/spline.hpp

@@ -85,6 +85,7 @@ namespace netgen
     void PrintCoeff (ostream & ost) const;
 
     virtual void GetCoeff (Vector & coeffs) const = 0;
+    virtual void GetCoeff (Vector & coeffs, Point<D> p0) const { ; } 
 
     virtual void GetPoints (int n, Array<Point<D> > & points) const;
 
@@ -138,6 +139,7 @@ namespace netgen
     virtual const GeomPoint<D> & EndPI () const { return p2; }
     ///
     virtual void GetCoeff (Vector & coeffs) const;
+    virtual void GetCoeff (Vector & coeffs, Point<D> p0) const;
     
     virtual string GetType(void) const {return "line";}
 
@@ -186,6 +188,7 @@ namespace netgen
     virtual const GeomPoint<D> & EndPI () const { return p3; }
     ///
     virtual void GetCoeff (Vector & coeffs) const;
+    virtual void GetCoeff (Vector & coeffs, Point<D> p0) const;
     
     virtual string GetType(void) const {return "spline3";}
 
@@ -394,6 +397,19 @@ namespace netgen
     coeffs[5] = -dx * p1(1) + dy * p1(0);
   }
 
+  template<int D>
+  void LineSeg<D> :: GetCoeff (Vector & coeffs, Point<D> p) const
+  {
+    coeffs.SetSize(6);
+
+    double dx = p2(0) - p1(0);
+    double dy = p2(1) - p1(1);
+
+    coeffs[0] = coeffs[1] = coeffs[2] = 0;
+    coeffs[3] = -dy;
+    coeffs[4] = dx;
+    coeffs[5] = -dx * (p1(1)-p(1)) + dy * (p1(0)-p(0));
+  }
 
 
   template<int D>