move splines

This commit is contained in:
Joachim Schoeberl 2011-02-28 13:13:18 +00:00
parent 5204be0657
commit 901bd944a2
6 changed files with 216 additions and 1679 deletions

File diff suppressed because it is too large Load Diff

View File

@ -1,162 +0,0 @@
/*
JS, Nov 2007
The 2D/3D template-base classes should go into the libsrc/gprim directory
in geom2d only 2D - Geometry classes (with material properties etc.)
*/
#ifndef _FILE_SPLINEGEOMETRY
#define _FILE_SPLINEGEOMETRY
// #include "../csg/csgparser.hpp"
namespace netgen
{
///
extern void LoadBoundarySplines (const char * filename,
Array < GeomPoint<2> > & geompoints,
Array < SplineSeg<2>* > & splines,
double & elto0);
///
extern void PartitionBoundary (const Array < SplineSeg<2>* > & splines,
double h, double elto0,
Mesh & mesh2d);
// allow to turn off messages: cover all couts !!
extern int printmessage_importance;
template < int D >
class SplineGeometry
{
// protected:
public:
Array < GeomPoint<D> > geompoints;
Array < SplineSeg<D>* > splines;
double elto0;
Array<char*> materials;
Array<string*> bcnames;
Array<double> maxh;
Array<bool> quadmeshing;
Array<bool> tensormeshing;
Array<int> layer;
void AppendSegment(SplineSeg<D> * spline, const int leftdomain, const int rightdomain,
const int bc,
const double reffac, const bool hprefleft, const bool hprefright,
const int copyfrom);
public:
~SplineGeometry();
int Load (const Array<double> & raw_data, const int startpos = 0);
void Load (const char * filename);
// void CSGLoad (CSGScanner & scan);
void LoadData( ifstream & infile );
void LoadDataNew ( ifstream & infile );
void LoadDataV2 ( ifstream & infile );
void GetRawData (Array<double> & raw_data) const;
void CopyEdgeMesh (int from, int to, Mesh & mesh2d, Point3dTree & searchtree);
const Array<SplineSeg<D>*> & GetSplines () const
{ return splines; }
int GetNSplines (void) const { return splines.Size(); }
string GetSplineType (const int i) const { return splines[i]->GetType(); }
SplineSeg<D> & GetSpline (const int i) {return *splines[i];}
const SplineSeg<D> & GetSpline (const int i) const {return *splines[i];}
void GetBoundingBox (Box<D> & box) const;
Box<D> GetBoundingBox () const
{ Box<D> box; GetBoundingBox (box); return box; }
int GetNP () const { return geompoints.Size(); }
const GeomPoint<D> & GetPoint(int i) const { return geompoints[i]; }
void SetGrading (const double grading);
// void AppendPoint (const double x, const double y, const double reffac = 1., const bool hpref = false);
void AppendPoint (const Point<D> & p, const double reffac = 1., const bool hpref = false);
void AppendLineSegment (const int n1, const int n2,
const int leftdomain, const int rightdomain, const int bc = -1,
const double reffac = 1.,
const bool hprefleft = false, const bool hprefright = false,
const int copyfrom = -1);
void AppendSplineSegment (const int n1, const int n2, const int n3,
const int leftdomain, const int rightdomain, const int bc = -1,
const double reffac = 1.,
const bool hprefleft = false, const bool hprefright = false,
const int copyfrom = -1);
void AppendCircleSegment (const int n1, const int n2, const int n3,
const int leftdomain, const int rightdomain, const int bc = -1,
const double reffac = 1.,
const bool hprefleft = false, const bool hprefright = false,
const int copyfrom = -1);
void AppendDiscretePointsSegment (const Array< Point<D> > & points,
const int leftdomain, const int rightdomain, const int bc = -1,
const double reffac = 1.,
const bool hprefleft = false, const bool hprefright = false,
const int copyfrom = -1);
void TestComment ( ifstream & infile ) ;
void GetMaterial( const int domnr, char* & material );
double GetDomainMaxh ( const int domnr );
bool GetDomainQuadMeshing ( int domnr )
{
if ( quadmeshing.Size() ) return quadmeshing[domnr-1];
else return false;
}
bool GetDomainTensorMeshing ( int domnr )
{
if ( tensormeshing.Size() ) return tensormeshing[domnr-1];
else return false;
}
int GetDomainLayer ( int domnr )
{
if ( layer.Size() ) return layer[domnr-1];
else return 1;
}
string GetBCName ( const int bcnr ) const;
string * BCNamePtr ( const int bcnr );
};
/*
void MeshFromSpline2D (SplineGeometry<2> & geometry,
Mesh *& mesh,
MeshingParameters & mp);
*/
class SplineGeometry2d : public SplineGeometry<2>, public NetgenGeometry
{
public:
virtual ~SplineGeometry2d();
virtual int GenerateMesh (Mesh*& mesh, MeshingParameters & mparam,
int perfstepsstart, int perfstepsend);
void PartitionBoundary (double h, Mesh & mesh2d);
virtual Refinement & GetRefinement () const;
};
}
#endif // _FILE_SPLINEGEOMETRY

View File

@ -6,13 +6,11 @@ Spline curve for Mesh generator
#include <mystdlib.h>
#include <linalg.hpp>
#include <meshing.hpp>
#include <geometry2d.hpp>
#include <gprim.hpp>
#include "spline.hpp"
namespace netgen
{
#include "spline.hpp"
// just for testing (JS)
template <int D>
@ -236,74 +234,13 @@ namespace netgen
(b2p - b2*fac1) * v2 +
(b3p - b3*fac1) * v3;
/*
second = (b1pp/w - b1p*fac1 - b1*fac2) * v1 +
(b2pp/w - b2p*fac1 - b2*fac2) * v2 +
(b3pp/w - b3p*fac1 - b3*fac2) * v3;
*/
// JS: 2 was missing
second = (b1pp/w - 2*b1p*fac1 - b1*fac2) * v1 +
(b2pp/w - 2*b2p*fac1 - b2*fac2) * v2 +
(b3pp/w - 2*b3p*fac1 - b3*fac2) * v3;
}
void CalcPartition (double l, double h, double h1, double h2,
double hcurve, double elto0, Array<double> & points)
{
// cout << "calcpart, h = " << h << ", h1 = " << h1 << ", h2 = " << h2 << ", hcurve = " << hcurve << endl;
int i, j, n, nel;
double sum, t, dt, fun, fperel, oldf, f;
n = 1000;
points.SetSize (0);
sum = 0;
dt = l / n;
t = 0.5 * dt;
for (i = 1; i <= n; i++)
{
fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2);
sum += dt / fun;
t += dt;
}
nel = int (sum+1);
fperel = sum / nel;
points.Append (0);
i = 1;
oldf = 0;
t = 0.5 * dt;
for (j = 1; j <= n && i < nel; j++)
{
fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2);
f = oldf + dt / fun;
while (f > i * fperel && i < nel)
{
points.Append ( (l/n) * (j-1 + (i * fperel - oldf) / (f - oldf)) );
i++;
}
oldf = f;
t += dt;
}
points.Append (l);
}
template<>
double SplineSeg3<2> :: MaxCurvature(void) const
{
@ -338,4 +275,11 @@ namespace netgen
template class SplineSeg3<2>;
template class SplineSeg3<3>;
}

View File

@ -11,8 +11,6 @@ namespace netgen
{
void CalcPartition (double l, double h, double h1, double h2,
double hcurve, double elto0, Array<double> & points);
/*
Spline curves for 2D mesh generation
@ -47,31 +45,7 @@ namespace netgen
class SplineSeg
{
public:
/// left domain
int leftdom;
/// right domain
int rightdom;
/// refinement at line
double reffak;
/// maximal h;
double hmax;
/// boundary condition number
int bc;
/// copy spline mesh from other spline (-1.. do not copy)
int copyfrom;
/// perfrom anisotropic refinement (hp-refinement) to edge
bool hpref_left;
/// perfrom anisotropic refinement (hp-refinement) to edge
bool hpref_right;
///
int layer;
SplineSeg ()
{
layer = 1;
}
SplineSeg () { ; }
/// calculates length of curve
virtual double Length () const;
/// returns point at curve, 0 <= t <= 1
@ -83,9 +57,8 @@ namespace netgen
Point<D> & point,
Vec<D> & first,
Vec<D> & second) const {;}
/// partitionizes curve
void Partition (double h, double elto0,
Mesh & mesh, Point3dTree & searchtree, int segnr) const;
/// returns initial point on curve
virtual const GeomPoint<D> & StartPI () const = 0;
/// returns terminal point on curve
@ -99,7 +72,7 @@ namespace netgen
virtual void GetCoeff (Vector & coeffs) const = 0;
virtual void GetPoints (int n, Array<Point<D> > & points);
virtual void GetPoints (int n, Array<Point<D> > & points) const;
/** calculates (2D) lineintersections:
for lines $$ a x + b y + c = 0 $$ the interecting points are calculated
@ -279,9 +252,6 @@ namespace netgen
// calculates length of spline-curve
template<int D>
double SplineSeg<D> :: Length () const
@ -303,126 +273,8 @@ namespace netgen
}
// partitionizes spline curve
template<int D>
void SplineSeg<D> :: Partition (double h, double elto0,
Mesh & mesh, Point3dTree & searchtree, int segnr) const
{
int i, j;
double l; // , r1, r2, ra;
double lold, dt, frac;
int n = 100;
Point<D> p, pold, mark, oldmark;
Array<double> curvepoints;
double edgelength, edgelengthold;
l = Length();
double h1 = min (StartPI().hmax, h/StartPI().refatpoint);
double h2 = min (EndPI().hmax, h/EndPI().refatpoint);
double hcurve = min (hmax, h/reffak);
CalcPartition (l, h, h1, h2, hcurve, elto0, curvepoints);
// cout << "curvepoints = " << curvepoints << endl;
dt = 1.0 / n;
l = 0;
j = 1;
pold = GetPoint (0);
lold = 0;
oldmark = pold;
edgelengthold = 0;
Array<int> locsearch;
for (i = 1; i <= n; i++)
{
p = GetPoint (i*dt);
l = lold + Dist (p, pold);
while (j < curvepoints.Size() && (l >= curvepoints[j] || i == n))
{
frac = (curvepoints[j]-lold) / (l-lold);
edgelength = i*dt + (frac-1)*dt;
// mark = pold + frac * (p-pold);
mark = GetPoint (edgelength);
// cout << "mark = " << mark << " =?= " << GetPoint (edgelength) << endl;
{
PointIndex pi1 = -1, pi2 = -1;
Point3d mark3(mark(0), mark(1), 0);
Point3d oldmark3(oldmark(0), oldmark(1), 0);
Vec<3> v (1e-4*h, 1e-4*h, 1e-4*h);
searchtree.GetIntersecting (oldmark3 - v, oldmark3 + v, locsearch);
for (int k = 0; k < locsearch.Size(); k++)
if ( mesh[PointIndex(locsearch[k])].GetLayer() == layer)
pi1 = locsearch[k];
// if (locsearch.Size()) pi1 = locsearch[0];
searchtree.GetIntersecting (mark3 - v, mark3 + v, locsearch);
for (int k = 0; k < locsearch.Size(); k++)
if ( mesh[PointIndex(locsearch[k])].GetLayer() == layer)
pi2 = locsearch[k];
// if (locsearch.Size()) pi2 = locsearch[0];
/*
for (PointIndex pk = PointIndex::BASE;
pk < mesh.GetNP()+PointIndex::BASE; pk++)
{
if (Dist (mesh[pk], oldmark3) < 1e-4 * h) pi1 = pk;
if (Dist (mesh[pk], mark3) < 1e-4 * h) pi2 = pk;
}
*/
// cout << "pi1 = " << pi1 << endl;
// cout << "pi2 = " << pi2 << endl;
if (pi1 == -1)
{
pi1 = mesh.AddPoint(oldmark3, layer);
searchtree.Insert (oldmark3, pi1);
}
if (pi2 == -1)
{
pi2 = mesh.AddPoint(mark3, layer);
searchtree.Insert (mark3, pi2);
}
Segment seg;
seg.edgenr = segnr;
seg.si = bc; // segnr;
seg[0] = pi1;
seg[1] = pi2;
seg.domin = leftdom;
seg.domout = rightdom;
seg.epgeominfo[0].edgenr = segnr;
seg.epgeominfo[0].dist = edgelengthold;
seg.epgeominfo[1].edgenr = segnr;
seg.epgeominfo[1].dist = edgelength;
seg.singedge_left = hpref_left;
seg.singedge_right = hpref_right;
mesh.AddSegment (seg);
}
oldmark = mark;
edgelengthold = edgelength;
j++;
}
pold = p;
lold = l;
}
}
template<int D>
void SplineSeg<D> :: GetPoints (int n, Array<Point<D> > & points)
void SplineSeg<D> :: GetPoints (int n, Array<Point<D> > & points) const
{
points.SetSize (n);
if (n >= 2)
@ -873,19 +725,7 @@ namespace netgen
return pts[segnr] + rest*Vec<D>(pts[segnr+1]-pts[segnr]);
}
typedef GeomPoint<2> GeomPoint2d;
typedef SplineSeg<2> SplineSegment;
typedef LineSeg<2> LineSegment;
typedef SplineSeg3<2> SplineSegment3;
typedef CircleSeg<2> CircleSegment;
typedef DiscretePointsSeg<2> DiscretePointsSegment;
}
#endif

View File

@ -0,0 +1,134 @@
/*
2d Spline curve for Mesh generator
*/
#include <mystdlib.h>
#include <linalg.hpp>
#include <gprim.hpp>
#include "splinegeometry.hpp"
namespace netgen
{
template<int D>
SplineGeometry<D> :: ~SplineGeometry()
{
for(int i = 0; i < splines.Size(); i++)
delete splines[i];
}
template<int D>
void SplineGeometry<D> :: GetRawData (Array<double> & raw_data) const
{
raw_data.Append(D);
// raw_data.Append(elto0);
raw_data.Append(splines.Size());
for(int i=0; i<splines.Size(); i++)
splines[i]->GetRawData(raw_data);
}
template<int D>
int SplineGeometry<D> :: Load (const Array<double> & raw_data, const int startpos)
{
int pos = startpos;
if(raw_data[pos] != D)
throw NgException("wrong dimension of spline raw_data");
pos++;
// elto0 = raw_data[pos]; pos++;
splines.SetSize(int(raw_data[pos]));
pos++;
Array< Point<D> > pts(3);
for(int i=0; i<splines.Size(); i++)
{
int type = int(raw_data[pos]);
pos++;
for(int j=0; j<type; j++)
for(int k=0; k<D; k++)
{
pts[j](k) = raw_data[pos];
pos++;
}
if (type == 2)
{
splines[i] = new LineSeg<D>(GeomPoint<D>(pts[0],1),
GeomPoint<D>(pts[1],1));
}
else if (type == 3)
{
splines[i] = new SplineSeg3<D>(GeomPoint<D>(pts[0],1),
GeomPoint<D>(pts[1],1),
GeomPoint<D>(pts[2],1));
}
else
throw NgException("something wrong with spline raw data");
}
return pos;
}
template<int D>
void SplineGeometry<D> :: GetBoundingBox (Box<D> & box) const
{
if (!splines.Size())
{
Point<D> auxp = 0.;
box.Set (auxp);
return;
}
Array<Point<D> > points;
for (int i = 0; i < splines.Size(); i++)
{
splines[i]->GetPoints (20, points);
if (i == 0) box.Set(points[0]);
for (int j = 0; j < points.Size(); j++)
box.Add (points[j]);
}
}
/*
template<int D>
void SplineGeometry<D> :: SetGrading (const double grading)
{
elto0 = grading;
}
*/
template<int D>
void SplineGeometry<D> :: AppendPoint (const Point<D> & p, const double reffac, const bool hpref)
{
geompoints.Append (GeomPoint<D>(p, reffac));
geompoints.Last().hpref = hpref;
}
template class SplineGeometry<2>;
template class SplineGeometry<3>;
}

View File

@ -0,0 +1,68 @@
/*
JS, Nov 2007
The 2D/3D template-base classes should go into the libsrc/gprim directory
in geom2d only 2D - Geometry classes (with material properties etc.)
*/
#include "spline.hpp"
#ifndef _FILE_SPLINEGEOMETRY
#define _FILE_SPLINEGEOMETRY
namespace netgen
{
template < int D >
class SplineGeometry
{
// protected:
public:
Array < GeomPoint<D> > geompoints;
Array < SplineSeg<D>* > splines;
~SplineGeometry();
int Load (const Array<double> & raw_data, const int startpos = 0);
void GetRawData (Array<double> & raw_data) const;
const Array<SplineSeg<D>*> & GetSplines () const
{ return splines; }
int GetNSplines (void) const { return splines.Size(); }
string GetSplineType (const int i) const { return splines[i]->GetType(); }
SplineSeg<D> & GetSpline (const int i) {return *splines[i];}
const SplineSeg<D> & GetSpline (const int i) const {return *splines[i];}
void GetBoundingBox (Box<D> & box) const;
Box<D> GetBoundingBox () const
{ Box<D> box; GetBoundingBox (box); return box; }
int GetNP () const { return geompoints.Size(); }
const GeomPoint<D> & GetPoint(int i) const { return geompoints[i]; }
// void SetGrading (const double grading);
void AppendPoint (const Point<D> & p, const double reffac = 1., const bool hpref = false);
void AppendSegment(SplineSeg<D> * spline)
{
splines.Append (spline);
}
};
}
#endif // _FILE_SPLINEGEOMETRY