mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-24 21:10:33 +05:00
B-spline curves in 2D geometry
This commit is contained in:
Joachim Schoeberl
parent
ba0319d388
commit
0291271c05
@ -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);
|
||||
|
||||
@ -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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -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;
|
||||
}
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -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
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user