From 0291271c05a0427baa5478c8b2b5826c9558e1d8 Mon Sep 17 00:00:00 2001
From: Joachim Schoeberl <joachim.schoeberl@tuwien.ac.at>
Date: Thu, 1 Sep 2011 16:40:00 +0000
Subject: [PATCH] B-spline curves in 2D geometry

---
 libsrc/geom2d/geometry2d.cpp     |  22 ++++++
 libsrc/gprim/spline.hpp          | 115 ++++++++++++++++++++++++++++++-
 libsrc/occ/Partition_Inter3d.cxx |   4 +-
 libsrc/occ/Partition_Loop2d.cxx  |   2 +
 libsrc/occ/Partition_Spliter.cxx |   3 +-
 5 files changed, 143 insertions(+), 3 deletions(-)

diff --git a/libsrc/geom2d/geometry2d.cpp b/libsrc/geom2d/geometry2d.cpp
index 5017e9c2..e8b10c5e 100644
--- a/libsrc/geom2d/geometry2d.cpp
+++ b/libsrc/geom2d/geometry2d.cpp
@@ -10,6 +10,8 @@
 namespace netgen
 {
 
+
+
   SplineGeometry2d :: ~SplineGeometry2d()
   {
     for ( int i = 0; i < bcnames.Size(); i++ )
@@ -660,6 +662,26 @@ namespace netgen
 		  
 		    spline = new DiscretePointsSeg<D> (pts);
 		  }
+		else if (strcmp (buf, "bsplinepoints") == 0)
+		  {
+		    int npts,order;
+		    infile >> npts;    
+		    infile >> order;
+		    Array< Point<D> > pts(npts);
+		    for (int j = 0; j < npts; j++)
+		      for(int k=0; k<D; k++)
+			infile >> pts[j](k);	    		    
+		    if(order<2)		      
+			cerr<<"Minimum order of 2 is required!!"<<endl;
+		    else if(order==2)
+		      spline = new BSplineSeg<D,2> (pts);
+		      else if(order==3)
+			spline = new BSplineSeg<D,3> (pts);
+		      else if(order==4)
+			spline = new BSplineSeg<D,4> (pts);
+		      else if(order>4)		      
+			cerr<<"Maximum allowed order is 4!!"<<endl;
+		  }
 	      
 		//      infile >> spline->reffak;
 		SplineSegExt * spex = new SplineSegExt (*spline);
diff --git a/libsrc/gprim/spline.hpp b/libsrc/gprim/spline.hpp
index 411c4472..3034fae3 100644
--- a/libsrc/gprim/spline.hpp
+++ b/libsrc/gprim/spline.hpp
@@ -34,7 +34,7 @@ namespace netgen
 
     ///
     GeomPoint (const Point<D> & ap, double aref = 1, bool ahpref=false)
-      : Point<D>(ap), refatpoint(aref), hpref(ahpref) { ; }
+      : Point<D>(ap), refatpoint(aref), hmax(1e99), hpref(ahpref) { ; }
   };
 
 
@@ -252,6 +252,8 @@ namespace netgen
 
 
 
+
+
   // calculates length of spline-curve
   template<int D>
   double SplineSeg<D> :: Length () const
@@ -527,6 +529,117 @@ namespace netgen
     return pts[segnr] + rest*Vec<D>(pts[segnr+1]-pts[segnr]);
   }
 
+
+
+
+
+
+
+  
+
+  // *************************************
+  // Template for B-Splines of order ORDER
+  // *************************************
+  
+  template<int D, int ORDER>
+  class BSplineSeg : public SplineSeg<D>
+  {
+    Array<Point<D> > pts;
+    GeomPoint<D> p1n, p2n;    
+    Array<int> ti; 
+
+  public:
+    ///
+    BSplineSeg (const Array<Point<D> > & apts);
+    ///
+    virtual ~BSplineSeg();
+    ///
+    virtual Point<D> GetPoint (double t) const;
+    ///
+    virtual const GeomPoint<D> & StartPI () const { return p1n; };
+    ///
+    virtual const GeomPoint<D> & EndPI () const { return p2n; }
+    ///
+    virtual void GetCoeff (Vector & coeffs) const {;}
+
+    virtual double MaxCurvature(void) const {return 1;}
+  };
+
+  // Constructor
+  template<int D,int ORDER>
+  BSplineSeg<D,ORDER> :: BSplineSeg (const Array<Point<D> > & apts)
+    : pts (apts)
+  { 
+    /*
+    for(int i=0; i<D; i++)
+      {
+	p1n(i) = apts[0](i);
+	p2n(i) = apts.Last()(i);
+      }
+    */
+    p1n = apts[0];
+    p2n = apts.Last();
+
+    /*
+    p1n.refatpoint = 1;
+    p2n.refatpoint = 1;
+    p1n.hmax = 1e99;
+    p2n.hmax = 1e99;
+    */
+
+    int m=pts.Size()+ORDER;
+    ti.SetSize(m);
+    // b.SetSize(m-1);
+    ti=0;    
+    //    b=0.0;
+    for(int i=ORDER;i<m-ORDER+1;i++)
+      ti[i]=i-ORDER+1;   
+    for(int i=m-ORDER+1;i<m;i++)
+      ti[i]=m-2*ORDER+1;
+  }
+  // Destructor
+  template<int D,int ORDER>
+  BSplineSeg<D, ORDER> :: ~BSplineSeg ()
+  { ; }
+
+
+  // GetPoint Method...(evaluation of BSpline Curve)
+  template<int D,int ORDER>
+  Point<D> BSplineSeg<D,ORDER> :: GetPoint (double t_in) const
+  {    
+    int m=pts.Size()+ORDER;           
+
+    double t = t_in * (m-2*ORDER+1);    
+
+    double b[ORDER];
+    
+    int interval_nr = int(t)+ORDER-1;    
+    if (interval_nr < ORDER-1) interval_nr = ORDER-1;
+    if (interval_nr > m-ORDER-1) interval_nr = m-ORDER-1;
+
+    b[ORDER-1] = 1.0;
+    
+    for(int degree=1;degree<ORDER;degree++)
+      for (int k = 0; k <= degree; k++)
+	{
+	  int j = interval_nr-degree+k;
+	  double bnew = 0;
+
+	  if (k != 0) 
+	    bnew += (t-ti[j]) / ( ti[j+degree]-ti[j] ) * b[k-degree+ORDER-1];
+	  if (k != degree)
+	    bnew += (ti[j+degree+1]-t) / ( ti[j+degree+1]-ti[j+1] ) * b[k-degree+ORDER];
+	  b[k-degree+ORDER-1] = bnew;
+	}
+
+    Point<D> p = 0.0;
+    for(int i=0; i < ORDER; i++) 
+      p += b[i] * Vec<D> (pts[i+interval_nr-ORDER+1]);
+    return p;
+  }
+
+
+
 }
 
 
diff --git a/libsrc/occ/Partition_Inter3d.cxx b/libsrc/occ/Partition_Inter3d.cxx
index 3775e5b3..a47ff2df 100644
--- a/libsrc/occ/Partition_Inter3d.cxx
+++ b/libsrc/occ/Partition_Inter3d.cxx
@@ -243,7 +243,9 @@ static void PutInBounds (const TopoDS_Face&          F,
       Standard_Integer i, nbExt = anExtPS.NbExt();
       Extrema_POnSurf aPOnSurf;
       for (i = 1; i <= nbExt; ++i )
-        if (anExtPS.Value( i ) <= TolE) {
+	if (anExtPS.Value( i ) <= TolE)               // V6.3
+	  // if (anExtPS.SquareDistance( i ) <= TolE)   // V6.5
+	  {
           aPOnSurf = anExtPS.Point( i );
           break;
         }
diff --git a/libsrc/occ/Partition_Loop2d.cxx b/libsrc/occ/Partition_Loop2d.cxx
index 3c1bd5c2..40872bd9 100644
--- a/libsrc/occ/Partition_Loop2d.cxx
+++ b/libsrc/occ/Partition_Loop2d.cxx
@@ -37,6 +37,7 @@
 #include <TopTools_DataMapIteratorOfDataMapOfShapeListOfShape.hxx>
 #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
 #include <TopTools_DataMapOfShapeInteger.hxx>
+// #include <TopTools_DataMapOfShapeReal.hxx>    V6.5
 #include <TopTools_DataMapOfShapeShape.hxx>
 #include <TopTools_IndexedMapOfShape.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
@@ -520,6 +521,7 @@ static void prepareDegen (const TopoDS_Edge&                        DegEdge,
 
   // avoid intersecting twice the same edge
   BRepOffset_DataMapOfShapeReal EUMap ( EdgesList.Extent() );
+  // TopTools_DataMapOfShapeReal EUMap ( EdgesList.Extent() );   // V6.5
 
   Standard_Real U, f, l;
   BRep_Tool::Range (DegEdge, f, l);
diff --git a/libsrc/occ/Partition_Spliter.cxx b/libsrc/occ/Partition_Spliter.cxx
index 3c13cf88..97bd8e33 100644
--- a/libsrc/occ/Partition_Spliter.cxx
+++ b/libsrc/occ/Partition_Spliter.cxx
@@ -1169,7 +1169,8 @@ static void findEqual (const TopTools_ListOfShape& EL1,
           for (; j<=nbj && ok; ++j) {
             if (Extrema.IsMin(j)) {
 	      hasMin = Standard_True;
-              ok = Extrema.Value(j) <= tol;
+	      ok = Extrema.Value(j) <= tol;  // V6.3
+	      // ok = Extrema.SquareDistance(j) <= tol;  // V6.5
 	    }
           }
         }