mirror of
https://github.com/NGSolve/netgen.git
synced 2024-11-14 18:08:33 +05:00
1283 lines
29 KiB
C++
1283 lines
29 KiB
C++
/*
|
|
|
|
2d Spline curve for Mesh generator
|
|
|
|
*/
|
|
|
|
|
|
#include <mystdlib.h>
|
|
#include <csg.hpp>
|
|
#include <geometry2d.hpp>
|
|
#include "meshing.hpp"
|
|
|
|
namespace netgen
|
|
{
|
|
|
|
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: LoadDataV2 ( ifstream & infile )
|
|
{
|
|
// new parser by Astrid Sinwel
|
|
|
|
PrintMessage (1, "Load 2D Geometry V2");
|
|
int nump, leftdom, rightdom;
|
|
Point<D> x;
|
|
int hi1, hi2, hi3;
|
|
double hd;
|
|
char buf[50], ch;
|
|
int pointnr;
|
|
|
|
string keyword;
|
|
|
|
Array < GeomPoint<D> > infilepoints (0);
|
|
Array <int> pointnrs (0);
|
|
nump = 0;
|
|
int numdomains = 0;
|
|
|
|
|
|
TestComment ( infile );
|
|
// refinement factor
|
|
infile >> elto0;
|
|
TestComment ( infile );
|
|
|
|
|
|
// test if next ch is a letter, i.e. new keyword starts
|
|
bool ischar = false;
|
|
|
|
while ( infile.good() )
|
|
{
|
|
infile >> keyword;
|
|
|
|
ischar = false;
|
|
|
|
if ( keyword == "points" )
|
|
{
|
|
PrintMessage (3, "load points");
|
|
infile.get(ch);
|
|
infile.putback(ch);
|
|
|
|
// test if ch is a letter
|
|
if ( int(ch) >= 65 && int(ch) <=90 )
|
|
ischar = true;
|
|
if ( int(ch) >= 97 && int(ch) <= 122 )
|
|
ischar = true;
|
|
|
|
while ( ! ischar )
|
|
{
|
|
TestComment ( infile );
|
|
infile >> pointnr;
|
|
// pointnrs 1-based
|
|
if ( pointnr > nump ) nump = pointnr;
|
|
pointnrs.Append(pointnr);
|
|
|
|
for(int j=0; j<D; j++)
|
|
infile >> x(j);
|
|
// hd is now optional, default 1
|
|
// infile >> hd;
|
|
hd = 1;
|
|
|
|
Flags flags;
|
|
|
|
|
|
// get flags,
|
|
ch = 'a';
|
|
// infile >> ch;
|
|
do
|
|
{
|
|
infile.get (ch);
|
|
// if another int-value, set refinement flag to this value
|
|
// (corresponding to old files)
|
|
if ( int (ch) >= 48 && int(ch) <= 57 )
|
|
{
|
|
infile.putback(ch);
|
|
infile >> hd;
|
|
infile.get(ch);
|
|
}
|
|
}
|
|
while (isspace(ch) && ch != '\n');
|
|
while (ch == '-')
|
|
{
|
|
char flag[100];
|
|
flag[0]='-';
|
|
infile >> (flag+1);
|
|
flags.SetCommandLineFlag (flag);
|
|
ch = 'a';
|
|
do {
|
|
infile.get (ch);
|
|
} while (isspace(ch) && ch != '\n');
|
|
}
|
|
if (infile.good())
|
|
infile.putback (ch);
|
|
|
|
if ( hd == 1 )
|
|
hd = flags.GetNumFlag ( "ref", 1.0);
|
|
// geompoints.Append (GeomPoint<D>(x, hd));
|
|
|
|
infilepoints.Append ( GeomPoint<D>(x, hd) );
|
|
infilepoints.Last().hpref = flags.GetDefineFlag ("hpref");
|
|
infilepoints.Last().hmax = flags.GetNumFlag ("maxh", 1e99);
|
|
|
|
TestComment(infile);
|
|
infile.get(ch);
|
|
infile.putback(ch);
|
|
|
|
// test if letter
|
|
if ( int(ch) >= 65 && int(ch) <=90 )
|
|
ischar = true;
|
|
if ( int(ch) >= 97 && int(ch) <= 122 )
|
|
ischar = true;
|
|
}
|
|
|
|
// infile.putback (ch);
|
|
|
|
geompoints.SetSize(nump);
|
|
for ( int i = 0; i < nump; i++ )
|
|
{
|
|
geompoints[pointnrs[i] - 1] = infilepoints[i];
|
|
geompoints[pointnrs[i] - 1].hpref = infilepoints[i].hpref;
|
|
}
|
|
TestComment(infile);
|
|
}
|
|
|
|
else if ( keyword == "segments" )
|
|
{
|
|
PrintMessage (3, "load segments");
|
|
|
|
bcnames.SetSize(0);
|
|
infile.get(ch);
|
|
infile.putback(ch);
|
|
int i = 0;
|
|
|
|
// test if ch is a letter
|
|
if ( int(ch) >= 65 && int(ch) <=90 )
|
|
ischar = true;
|
|
if ( int(ch) >= 97 && int(ch) <= 122 )
|
|
ischar = true;
|
|
|
|
while ( !ischar ) //ch != 'p' && ch != 'm' )
|
|
{
|
|
i++;
|
|
TestComment ( infile );
|
|
|
|
SplineSeg<D> * spline = 0;
|
|
TestComment ( infile );
|
|
|
|
infile >> leftdom >> rightdom;
|
|
|
|
if ( leftdom > numdomains ) numdomains = leftdom;
|
|
if ( rightdom > numdomains ) numdomains = rightdom;
|
|
|
|
|
|
infile >> buf;
|
|
// type of spline segement
|
|
if (strcmp (buf, "2") == 0)
|
|
{ // a line
|
|
infile >> hi1 >> hi2;
|
|
spline = new LineSeg<D>(geompoints[hi1-1],
|
|
geompoints[hi2-1]);
|
|
}
|
|
else if (strcmp (buf, "3") == 0)
|
|
{ // a rational spline
|
|
infile >> hi1 >> hi2 >> hi3;
|
|
spline = new SplineSeg3<D> (geompoints[hi1-1],
|
|
geompoints[hi2-1],
|
|
geompoints[hi3-1]);
|
|
}
|
|
else if (strcmp (buf, "4") == 0)
|
|
{ // an arc
|
|
infile >> hi1 >> hi2 >> hi3;
|
|
spline = new CircleSeg<D> (geompoints[hi1-1],
|
|
geompoints[hi2-1],
|
|
geompoints[hi3-1]);
|
|
break;
|
|
}
|
|
else if (strcmp (buf, "discretepoints") == 0)
|
|
{
|
|
int npts;
|
|
infile >> npts;
|
|
Array< Point<D> > pts(npts);
|
|
for (int j = 0; j < npts; j++)
|
|
for(int k=0; k<D; k++)
|
|
infile >> pts[j](k);
|
|
|
|
spline = new DiscretePointsSeg<D> (pts);
|
|
}
|
|
|
|
// infile >> spline->reffak;
|
|
spline -> leftdom = leftdom;
|
|
spline -> rightdom = rightdom;
|
|
splines.Append (spline);
|
|
|
|
|
|
// hd is now optional, default 1
|
|
// infile >> hd;
|
|
hd = 1;
|
|
infile >> ch;
|
|
|
|
// get refinement parameter, if it is there
|
|
//infile.get (ch);
|
|
// if another int-value, set refinement flag to this value
|
|
// (corresponding to old files)
|
|
|
|
if ( int (ch) >= 48 && int(ch) <= 57 )
|
|
{
|
|
infile.putback(ch);
|
|
infile >> hd;
|
|
infile >> ch ;
|
|
}
|
|
|
|
// get flags,
|
|
Flags flags;
|
|
while (ch == '-')
|
|
{
|
|
char flag[100];
|
|
flag[0]='-';
|
|
infile >> (flag+1);
|
|
flags.SetCommandLineFlag (flag);
|
|
ch = 'a';
|
|
infile >> ch;
|
|
}
|
|
|
|
if (infile.good())
|
|
infile.putback (ch);
|
|
|
|
splines.Last()->bc = int (flags.GetNumFlag ("bc", i+1));
|
|
splines.Last()->hpref_left = int (flags.GetDefineFlag ("hpref")) ||
|
|
int (flags.GetDefineFlag ("hprefleft"));
|
|
splines.Last()->hpref_right = int (flags.GetDefineFlag ("hpref")) ||
|
|
int (flags.GetDefineFlag ("hprefright"));
|
|
splines.Last()->copyfrom = int (flags.GetNumFlag ("copy", -1));
|
|
splines.Last()->reffak = flags.GetNumFlag ("ref", 1 );
|
|
splines.Last()->hmax = flags.GetNumFlag ("maxh", 1e99 );
|
|
if ( hd != 1 ) splines.Last()->reffak = hd;
|
|
|
|
if ( flags.StringFlagDefined("bcname") )
|
|
{
|
|
int mybc = splines.Last()->bc-1;
|
|
for ( int ii = bcnames.Size(); ii <= mybc; ii++ )
|
|
bcnames.Append ( new string ("default"));
|
|
if ( bcnames[mybc] ) delete bcnames[mybc];
|
|
bcnames[mybc] = new string (flags.GetStringFlag("bcname","") );
|
|
}
|
|
|
|
TestComment(infile);
|
|
infile.get(ch);
|
|
infile.putback(ch);
|
|
|
|
// test if ch is a letter
|
|
if ( int(ch) >= 65 && int(ch) <=90 )
|
|
ischar = true;
|
|
if ( int(ch) >= 97 && int(ch) <= 122 )
|
|
ischar = true;
|
|
|
|
}
|
|
|
|
infile.get(ch);
|
|
infile.putback(ch);
|
|
|
|
|
|
}
|
|
else if ( keyword == "materials" )
|
|
{
|
|
TestComment ( infile );
|
|
int domainnr;
|
|
char material[100];
|
|
|
|
if ( !infile.good() )
|
|
return;
|
|
|
|
materials.SetSize(numdomains) ;
|
|
maxh.SetSize ( numdomains ) ;
|
|
for ( int i = 0; i < numdomains; i++)
|
|
maxh[i] = 1000;
|
|
quadmeshing.SetSize ( numdomains );
|
|
quadmeshing = false;
|
|
tensormeshing.SetSize ( numdomains );
|
|
tensormeshing = false;
|
|
layer.SetSize ( numdomains );
|
|
layer = 1;
|
|
|
|
|
|
TestComment ( infile );
|
|
|
|
for ( int i=0; i<numdomains; i++)
|
|
materials [ i ] = new char[100];
|
|
|
|
for ( int i=0; i<numdomains && infile.good(); i++)
|
|
{
|
|
TestComment ( infile );
|
|
infile >> domainnr;
|
|
infile >> material;
|
|
|
|
strcpy (materials[domainnr-1], material);
|
|
|
|
Flags flags;
|
|
ch = 'a';
|
|
infile >> ch;
|
|
while (ch == '-')
|
|
{
|
|
char flag[100];
|
|
flag[0]='-';
|
|
infile >> (flag+1);
|
|
flags.SetCommandLineFlag (flag);
|
|
ch = 'a';
|
|
infile >> ch;
|
|
}
|
|
|
|
if (infile.good())
|
|
infile.putback (ch);
|
|
|
|
maxh[domainnr-1] = flags.GetNumFlag ( "maxh", 1000);
|
|
if (flags.GetDefineFlag("quad")) quadmeshing[domainnr-1] = true;
|
|
if (flags.GetDefineFlag("tensor")) tensormeshing[domainnr-1] = true;
|
|
layer[domainnr-1] = int(flags.GetNumFlag ("layer", 1));
|
|
}
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
// check if comments in a .in2d file...
|
|
// template <int D>
|
|
// void SplineGeometry<D> :: TestComment ( ifstream & infile )
|
|
// {
|
|
// bool comment = true;
|
|
// char ch;
|
|
// infile.get(ch);
|
|
// infile.putback(ch);
|
|
// int ii = 0;
|
|
// while ( comment == true && ii < 100)
|
|
// {
|
|
// infile.get(ch);
|
|
// if ( ch == '#' )
|
|
// while ( ch != '\n')
|
|
// {
|
|
// infile.get(ch);
|
|
// comment = false;
|
|
// }
|
|
// else if ( ch == '\n' )
|
|
// {
|
|
// comment = true;
|
|
// ii ++;
|
|
// }
|
|
// else
|
|
// {
|
|
// infile.putback(ch);
|
|
// comment = false;
|
|
// }
|
|
//
|
|
// infile.get(ch) ;
|
|
// if ( ch == '\n' || ch == '#' )
|
|
// {
|
|
// comment = true;
|
|
// }
|
|
// infile.putback(ch);
|
|
// if ( !comment ) break;
|
|
// }
|
|
// cerr << "** comment done" << endl;
|
|
// cerr << " * last char was " << ch << endl;
|
|
// return;
|
|
//
|
|
// }
|
|
|
|
// herbert: fixed TestComment
|
|
template <int D>
|
|
void SplineGeometry<D> :: TestComment ( ifstream & infile )
|
|
{
|
|
bool comment = true;
|
|
char ch;
|
|
while ( comment == true && !infile.eof() ) {
|
|
infile.get(ch);
|
|
if ( ch == '#' ) { // skip comments
|
|
while ( ch != '\n' && !infile.eof() ) {
|
|
infile.get(ch);
|
|
}
|
|
}
|
|
else if ( ch == '\n' ) { // skip empty lines
|
|
;
|
|
}
|
|
else if ( isspace(ch) ) { // skip whitespaces
|
|
;
|
|
}
|
|
else { // end of comment
|
|
infile.putback(ch);
|
|
comment = false;
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
|
|
|
|
template<int D>
|
|
SplineGeometry<D> :: ~SplineGeometry()
|
|
{
|
|
for(int i=0; i<splines.Size(); i++)
|
|
delete splines[i];
|
|
splines.DeleteAll();
|
|
geompoints.DeleteAll();
|
|
for (int i=0; i<materials.Size(); i++)
|
|
delete materials[i];
|
|
for ( int i = 0; i < bcnames.Size(); i++ )
|
|
if ( bcnames[i] ) delete bcnames[i];
|
|
}
|
|
|
|
|
|
|
|
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));
|
|
//(*testout) << "appending line segment "
|
|
// << pts[0] << " -- " << pts[1] << endl;
|
|
}
|
|
else if (type == 3)
|
|
{
|
|
splines[i] = new SplineSeg3<D>(GeomPoint<D>(pts[0],1),
|
|
GeomPoint<D>(pts[1],1),
|
|
GeomPoint<D>(pts[2],1));
|
|
//(*testout) << "appending spline segment "
|
|
// << pts[0] << " -- " << pts[1] << " -- " << pts[2] << endl;
|
|
|
|
}
|
|
else
|
|
throw NgException("something wrong with spline raw data");
|
|
|
|
}
|
|
return pos;
|
|
}
|
|
|
|
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>
|
|
void SplineGeometry<D> :: CSGLoad (CSGScanner & scan)
|
|
{
|
|
double hd;
|
|
Point<D> x;
|
|
int nump, numseg;
|
|
|
|
//scan.ReadNext();
|
|
scan >> nump >> ';';
|
|
|
|
hd = 1;
|
|
geompoints.SetSize(nump);
|
|
for(int i = 0; i<nump; i++)
|
|
{
|
|
if(D==2)
|
|
scan >> x(0) >> ',' >> x(1) >> ';';
|
|
else if(D==3)
|
|
scan >> x(0) >> ',' >> x(1) >> ',' >> x(2) >> ';';
|
|
|
|
geompoints[i] = GeomPoint<D>(x,hd);
|
|
}
|
|
|
|
scan >> numseg;// >> ';';
|
|
|
|
splines.SetSize(numseg);
|
|
|
|
int pnums,pnum1,pnum2,pnum3;
|
|
|
|
|
|
for(int i = 0; i<numseg; i++)
|
|
{
|
|
scan >> ';' >> pnums >> ',';
|
|
if (pnums == 2)
|
|
{
|
|
scan >> pnum1 >> ',' >> pnum2;// >> ';';
|
|
splines[i] = new LineSeg<D>(geompoints[pnum1-1],
|
|
geompoints[pnum2-1]);
|
|
}
|
|
else if (pnums == 3)
|
|
{
|
|
scan >> pnum1 >> ',' >> pnum2 >> ','
|
|
>> pnum3;// >> ';';
|
|
splines[i] = new SplineSeg3<D>(geompoints[pnum1-1],
|
|
geompoints[pnum2-1],
|
|
geompoints[pnum3-1]);
|
|
}
|
|
else if (pnums == 4)
|
|
{
|
|
scan >> pnum1 >> ',' >> pnum2 >> ','
|
|
>> pnum3;// >> ';';
|
|
splines[i] = new CircleSeg<D>(geompoints[pnum1-1],
|
|
geompoints[pnum2-1],
|
|
geompoints[pnum3-1]);
|
|
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: Load (const char * filename)
|
|
{
|
|
|
|
ifstream infile;
|
|
Point<D> x;
|
|
char buf[50];
|
|
|
|
|
|
infile.open (filename);
|
|
|
|
if ( ! infile.good() )
|
|
throw NgException(string ("Input file '") +
|
|
string (filename) +
|
|
string ("' not available!"));
|
|
|
|
TestComment ( infile );
|
|
|
|
infile >> buf; // file recognition
|
|
|
|
tensormeshing.SetSize(0);
|
|
quadmeshing.SetSize(0);
|
|
|
|
TestComment ( infile );
|
|
if ( strcmp (buf, "splinecurves2dnew") == 0 )
|
|
{
|
|
LoadDataNew ( infile );
|
|
}
|
|
else if ( strcmp (buf, "splinecurves2dv2") == 0 )
|
|
{
|
|
LoadDataV2 ( infile );
|
|
}
|
|
else
|
|
{
|
|
LoadData(infile );
|
|
}
|
|
infile.close();
|
|
}
|
|
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: LoadDataNew ( ifstream & infile )
|
|
{
|
|
|
|
int nump, numseg, leftdom, rightdom;
|
|
Point<D> x;
|
|
int hi1, hi2, hi3;
|
|
double hd;
|
|
char buf[50], ch;
|
|
int pointnr;
|
|
|
|
|
|
TestComment ( infile );
|
|
infile >> elto0;
|
|
TestComment ( infile );
|
|
|
|
infile >> nump;
|
|
geompoints.SetSize(nump);
|
|
|
|
for (int i = 0; i < nump; i++)
|
|
{
|
|
TestComment ( infile );
|
|
infile >> pointnr;
|
|
if ( pointnr > nump )
|
|
{
|
|
throw NgException(string ("Point number greater than total number of points") );
|
|
}
|
|
for(int j=0; j<D; j++)
|
|
infile >> x(j);
|
|
|
|
|
|
// hd is now optional, default 1
|
|
// infile >> hd;
|
|
hd = 1;
|
|
|
|
Flags flags;
|
|
|
|
|
|
// get flags,
|
|
ch = 'a';
|
|
// infile >> ch;
|
|
do
|
|
{
|
|
|
|
infile.get (ch);
|
|
// if another int-value, set refinement flag to this value
|
|
// (corresponding to old files)
|
|
if ( int (ch) >= 48 && int(ch) <= 57 )
|
|
{
|
|
infile.putback(ch);
|
|
infile >> hd;
|
|
infile.get(ch);
|
|
}
|
|
}
|
|
while (isspace(ch) && ch != '\n');
|
|
while (ch == '-')
|
|
{
|
|
char flag[100];
|
|
flag[0]='-';
|
|
infile >> (flag+1);
|
|
flags.SetCommandLineFlag (flag);
|
|
ch = 'a';
|
|
do {
|
|
infile.get (ch);
|
|
} while (isspace(ch) && ch != '\n');
|
|
}
|
|
|
|
if (infile.good())
|
|
infile.putback (ch);
|
|
|
|
if ( hd == 1 )
|
|
hd = flags.GetNumFlag ( "ref", 1.0);
|
|
// geompoints.Append (GeomPoint<D>(x, hd));
|
|
geompoints[pointnr-1] = GeomPoint<D>(x, hd);
|
|
geompoints[pointnr-1].hpref = flags.GetDefineFlag ("hpref");
|
|
}
|
|
|
|
TestComment ( infile );
|
|
|
|
infile >> numseg;
|
|
bcnames.SetSize(numseg);
|
|
for ( int i = 0; i < numseg; i++ )
|
|
bcnames[i] = 0;//new"default";
|
|
|
|
SplineSeg<D> * spline = 0;
|
|
for (int i = 0; i < numseg; i++)
|
|
{
|
|
TestComment ( infile );
|
|
|
|
infile >> leftdom >> rightdom;
|
|
|
|
// cout << "add spline " << i << ", left = " << leftdom << endl;
|
|
|
|
infile >> buf;
|
|
// type of spline segement
|
|
if (strcmp (buf, "2") == 0)
|
|
{ // a line
|
|
infile >> hi1 >> hi2;
|
|
spline = new LineSeg<D> (geompoints[hi1-1],
|
|
geompoints[hi2-1]);
|
|
}
|
|
else if (strcmp (buf, "3") == 0)
|
|
{ // a rational spline
|
|
infile >> hi1 >> hi2 >> hi3;
|
|
spline = new SplineSeg3<D> (geompoints[hi1-1],
|
|
geompoints[hi2-1],
|
|
geompoints[hi3-1]);
|
|
}
|
|
else if (strcmp (buf, "4") == 0)
|
|
{ // an arc
|
|
infile >> hi1 >> hi2 >> hi3;
|
|
spline = new CircleSeg<D> (geompoints[hi1-1],
|
|
geompoints[hi2-1],
|
|
geompoints[hi3-1]);
|
|
// break;
|
|
}
|
|
else if (strcmp (buf, "discretepoints") == 0)
|
|
{
|
|
int npts;
|
|
infile >> npts;
|
|
Array< Point<D> > pts(npts);
|
|
for (int j = 0; j < npts; j++)
|
|
for(int k=0; k<D; k++)
|
|
infile >> pts[j](k);
|
|
|
|
spline = new DiscretePointsSeg<D> (pts);
|
|
}
|
|
|
|
// infile >> spline->reffak;
|
|
spline -> leftdom = leftdom;
|
|
spline -> rightdom = rightdom;
|
|
splines.Append (spline);
|
|
|
|
// hd is now optional, default 1
|
|
// infile >> hd;
|
|
hd = 1;
|
|
infile >> ch;
|
|
|
|
// get refinement parameter, if it is there
|
|
// infile.get (ch);
|
|
// if another int-value, set refinement flag to this value
|
|
// (corresponding to old files)
|
|
if ( int (ch) >= 48 && int(ch) <= 57 )
|
|
{
|
|
infile.putback(ch);
|
|
infile >> hd;
|
|
infile >> ch ;
|
|
}
|
|
|
|
Flags flags;
|
|
while (ch == '-')
|
|
{
|
|
char flag[100];
|
|
flag[0]='-';
|
|
infile >> (flag+1);
|
|
flags.SetCommandLineFlag (flag);
|
|
ch = 'a';
|
|
infile >> ch;
|
|
}
|
|
|
|
if (infile.good())
|
|
infile.putback (ch);
|
|
|
|
splines.Last()->bc = int (flags.GetNumFlag ("bc", i+1));
|
|
splines.Last()->hpref_left = int (flags.GetDefineFlag ("hpref")) ||
|
|
int (flags.GetDefineFlag ("hprefleft"));
|
|
splines.Last()->hpref_right = int (flags.GetDefineFlag ("hpref")) ||
|
|
int (flags.GetDefineFlag ("hprefright"));
|
|
splines.Last()->copyfrom = int (flags.GetNumFlag ("copy", -1));
|
|
splines.Last()->reffak = flags.GetNumFlag ("ref", 1 );
|
|
splines.Last()->hmax = flags.GetNumFlag ("maxh", 1e99 );
|
|
|
|
if ( flags.StringFlagDefined("bcname") )
|
|
{
|
|
int mybc = splines.Last()->bc-1;
|
|
if ( bcnames[mybc] ) delete bcnames[mybc];
|
|
bcnames[mybc] = new string (flags.GetStringFlag("bcname","") );
|
|
}
|
|
|
|
if ( hd != 1 )
|
|
splines.Last()->reffak = hd;
|
|
}
|
|
if ( !infile.good() )
|
|
return;
|
|
TestComment ( infile );
|
|
int numdomains;
|
|
int domainnr;
|
|
char material[100];
|
|
|
|
if ( !infile.good() )
|
|
return;
|
|
|
|
infile >> numdomains;
|
|
materials.SetSize(numdomains) ;
|
|
maxh.SetSize ( numdomains ) ;
|
|
for ( int i = 0; i < numdomains; i++)
|
|
maxh[i] = 1000;
|
|
|
|
TestComment ( infile );
|
|
|
|
for ( int i=0; i<numdomains; i++)
|
|
materials [ i ] = new char (100);
|
|
|
|
for ( int i=0; i<numdomains && infile.good(); i++)
|
|
{
|
|
TestComment ( infile );
|
|
infile >> domainnr;
|
|
infile >> material;
|
|
strcpy(materials[domainnr-1], material);
|
|
|
|
Flags flags;
|
|
ch = 'a';
|
|
infile >> ch;
|
|
while (ch == '-')
|
|
{
|
|
char flag[100];
|
|
flag[0]='-';
|
|
infile >> (flag+1);
|
|
flags.SetCommandLineFlag (flag);
|
|
ch = 'a';
|
|
infile >> ch;
|
|
}
|
|
|
|
if (infile.good())
|
|
infile.putback (ch);
|
|
|
|
maxh[domainnr-1] = flags.GetNumFlag ( "maxh", 1000);
|
|
}
|
|
return;
|
|
}
|
|
|
|
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: LoadData ( ifstream & infile )
|
|
{
|
|
|
|
int nump, numseg, leftdom, rightdom;
|
|
Point<D> x;
|
|
int hi1, hi2, hi3;
|
|
double hd;
|
|
char buf[50], ch;
|
|
|
|
materials.SetSize(0);
|
|
maxh.SetSize(0);
|
|
infile >> elto0;
|
|
|
|
TestComment ( infile );
|
|
|
|
infile >> nump;
|
|
for (int i = 0; i < nump; i++)
|
|
{
|
|
TestComment ( infile );
|
|
for(int j=0; j<D; j++)
|
|
infile >> x(j);
|
|
infile >> hd;
|
|
|
|
Flags flags;
|
|
|
|
ch = 'a';
|
|
// infile >> ch;
|
|
do {
|
|
infile.get (ch);
|
|
} while (isspace(ch) && ch != '\n');
|
|
while (ch == '-')
|
|
{
|
|
char flag[100];
|
|
flag[0]='-';
|
|
infile >> (flag+1);
|
|
flags.SetCommandLineFlag (flag);
|
|
ch = 'a';
|
|
do {
|
|
infile.get (ch);
|
|
} while (isspace(ch) && ch != '\n');
|
|
}
|
|
|
|
if (infile.good())
|
|
infile.putback (ch);
|
|
|
|
geompoints.Append (GeomPoint<D>(x, hd));
|
|
geompoints.Last().hpref = flags.GetDefineFlag ("hpref");
|
|
geompoints.Last().hmax = 1e99;
|
|
}
|
|
|
|
PrintMessage (3, nump, " points loaded");
|
|
TestComment ( infile );
|
|
|
|
infile >> numseg;
|
|
bcnames.SetSize(numseg);
|
|
for ( int i = 0; i < numseg; i++ )
|
|
bcnames[i] = 0; // "default";
|
|
|
|
SplineSeg<D> * spline = 0;
|
|
|
|
PrintMessage (3, numseg, " segments loaded");
|
|
for (int i = 0; i < numseg; i++)
|
|
{
|
|
TestComment ( infile );
|
|
|
|
infile >> leftdom >> rightdom;
|
|
|
|
// cout << "add spline " << i << ", left = " << leftdom << ", right = " << rightdom << endl;
|
|
|
|
infile >> buf;
|
|
// type of spline segement
|
|
if (strcmp (buf, "2") == 0)
|
|
{ // a line
|
|
infile >> hi1 >> hi2;
|
|
spline = new LineSeg<D>(geompoints[hi1-1],
|
|
geompoints[hi2-1]);
|
|
}
|
|
else if (strcmp (buf, "3") == 0)
|
|
{ // a rational spline
|
|
infile >> hi1 >> hi2 >> hi3;
|
|
spline = new SplineSeg3<D> (geompoints[hi1-1],
|
|
geompoints[hi2-1],
|
|
geompoints[hi3-1]);
|
|
}
|
|
else if (strcmp (buf, "4") == 0)
|
|
{ // an arc
|
|
infile >> hi1 >> hi2 >> hi3;
|
|
spline = new CircleSeg<D> (geompoints[hi1-1],
|
|
geompoints[hi2-1],
|
|
geompoints[hi3-1]);
|
|
// break;
|
|
}
|
|
else if (strcmp (buf, "discretepoints") == 0)
|
|
{
|
|
int npts;
|
|
infile >> npts;
|
|
Array< Point<D> > pts(npts);
|
|
for (int j = 0; j < npts; j++)
|
|
for(int k=0; k<D; k++)
|
|
infile >> pts[j](k);
|
|
|
|
spline = new DiscretePointsSeg<D> (pts);
|
|
}
|
|
|
|
infile >> spline->reffak;
|
|
spline -> leftdom = leftdom;
|
|
spline -> rightdom = rightdom;
|
|
spline -> hmax = 1e99;
|
|
splines.Append (spline);
|
|
|
|
|
|
Flags flags;
|
|
ch = 'a';
|
|
infile >> ch;
|
|
while (ch == '-')
|
|
{
|
|
char flag[100];
|
|
flag[0]='-';
|
|
infile >> (flag+1);
|
|
flags.SetCommandLineFlag (flag);
|
|
ch = 'a';
|
|
infile >> ch;
|
|
}
|
|
|
|
if (infile.good())
|
|
infile.putback (ch);
|
|
|
|
splines.Last()->bc = int (flags.GetNumFlag ("bc", i+1));
|
|
splines.Last()->hpref_left = int (flags.GetDefineFlag ("hpref")) ||
|
|
int (flags.GetDefineFlag ("hprefleft"));
|
|
splines.Last()->hpref_right = int (flags.GetDefineFlag ("hpref")) ||
|
|
int (flags.GetDefineFlag ("hprefright"));
|
|
splines.Last()->copyfrom = int (flags.GetNumFlag ("copy", -1));
|
|
if ( flags.StringFlagDefined("bcname") )
|
|
{
|
|
int mybc = splines.Last()->bc-1;
|
|
if ( bcnames[mybc] ) delete bcnames[mybc];
|
|
bcnames[mybc] = new string (flags.GetStringFlag("bcname","") );
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: PartitionBoundary (double h, Mesh & mesh2d)
|
|
{
|
|
Box<D> bbox;
|
|
GetBoundingBox (bbox);
|
|
double dist = Dist (bbox.PMin(), bbox.PMax());
|
|
Point<3> pmin;
|
|
Point<3> pmax;
|
|
|
|
pmin(2) = -dist; pmax(2) = dist;
|
|
for(int j=0;j<D;j++)
|
|
{
|
|
pmin(j) = bbox.PMin()(j);
|
|
pmax(j) = bbox.PMax()(j);
|
|
}
|
|
|
|
Point3dTree searchtree (pmin, pmax);
|
|
|
|
for (int i = 0; i < splines.Size(); i++)
|
|
for (int side = 0; side <= 1; side++)
|
|
{
|
|
int dom = (side == 0) ? splines[i]->leftdom : splines[i]->rightdom;
|
|
if (dom != 0) splines[i] -> layer = GetDomainLayer (dom);
|
|
}
|
|
|
|
for (int i = 0; i < splines.Size(); i++)
|
|
if (splines[i]->copyfrom == -1)
|
|
{
|
|
// astrid - set boundary meshsize to domain meshsize h
|
|
// if no domain mesh size is given, the max h value from the bounding box is used
|
|
double minimum = min2 ( GetDomainMaxh ( splines[i]->leftdom ), GetDomainMaxh ( splines[i]->rightdom ) );
|
|
double maximum = max2 ( GetDomainMaxh ( splines[i]->leftdom ), GetDomainMaxh ( splines[i]->rightdom ) );
|
|
minimum = min2 ( minimum, h );
|
|
maximum = min2 ( maximum, h);
|
|
if ( minimum > 0 )
|
|
splines[i]->Partition(minimum, elto0, mesh2d, searchtree, i+1);
|
|
else if ( maximum > 0 )
|
|
splines[i]->Partition(maximum, elto0, mesh2d, searchtree, i+1);
|
|
else
|
|
splines[i]->Partition(h, elto0, mesh2d, searchtree, i+1);
|
|
}
|
|
else
|
|
{
|
|
CopyEdgeMesh (splines[i]->copyfrom, i+1, mesh2d, searchtree);
|
|
}
|
|
}
|
|
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: CopyEdgeMesh (int from, int to, Mesh & mesh, Point3dTree & searchtree)
|
|
{
|
|
int i;
|
|
|
|
Array<int, PointIndex::BASE> mappoints (mesh.GetNP());
|
|
Array<double, PointIndex::BASE> param (mesh.GetNP());
|
|
mappoints = -1;
|
|
param = 0;
|
|
|
|
Point3d pmin, pmax;
|
|
mesh.GetBox (pmin, pmax);
|
|
double diam2 = Dist2(pmin, pmax);
|
|
|
|
if (printmessage_importance>0)
|
|
cout << "copy edge, from = " << from << " to " << to << endl;
|
|
|
|
for (i = 1; i <= mesh.GetNSeg(); i++)
|
|
{
|
|
const Segment & seg = mesh.LineSegment(i);
|
|
if (seg.edgenr == from)
|
|
{
|
|
mappoints.Elem(seg[0]) = 1;
|
|
param.Elem(seg[0]) = seg.epgeominfo[0].dist;
|
|
|
|
mappoints.Elem(seg[1]) = 1;
|
|
param.Elem(seg[1]) = seg.epgeominfo[1].dist;
|
|
}
|
|
}
|
|
|
|
bool mapped = false;
|
|
for (i = 1; i <= mappoints.Size(); i++)
|
|
{
|
|
if (mappoints.Get(i) != -1)
|
|
{
|
|
Point<D> newp = splines.Get(to)->GetPoint (param.Get(i));
|
|
Point<3> newp3;
|
|
for(int j=0; j<min2(D,3); j++)
|
|
newp3(j) = newp(j);
|
|
for(int j=min2(D,3); j<3; j++)
|
|
newp3(j) = 0;
|
|
|
|
int npi = -1;
|
|
|
|
for (PointIndex pi = PointIndex::BASE;
|
|
pi < mesh.GetNP()+PointIndex::BASE; pi++)
|
|
if (Dist2 (mesh.Point(pi), newp3) < 1e-12 * diam2)
|
|
npi = pi;
|
|
|
|
if (npi == -1)
|
|
{
|
|
npi = mesh.AddPoint (newp3);
|
|
searchtree.Insert (newp3, npi);
|
|
}
|
|
|
|
mappoints.Elem(i) = npi;
|
|
|
|
mesh.GetIdentifications().Add (i, npi, to);
|
|
mapped = true;
|
|
}
|
|
}
|
|
if(mapped)
|
|
mesh.GetIdentifications().SetType(to,Identifications::PERIODIC);
|
|
|
|
// copy segments
|
|
int oldnseg = mesh.GetNSeg();
|
|
for (i = 1; i <= oldnseg; i++)
|
|
{
|
|
const Segment & seg = mesh.LineSegment(i);
|
|
if (seg.edgenr == from)
|
|
{
|
|
Segment nseg;
|
|
nseg.edgenr = to;
|
|
nseg.si = splines.Get(to)->bc;
|
|
nseg[0] = mappoints.Get(seg[0]);
|
|
nseg[1] = mappoints.Get(seg[1]);
|
|
nseg.domin = splines.Get(to)->leftdom;
|
|
nseg.domout = splines.Get(to)->rightdom;
|
|
|
|
nseg.epgeominfo[0].edgenr = to;
|
|
nseg.epgeominfo[0].dist = param.Get(seg[0]);
|
|
nseg.epgeominfo[1].edgenr = to;
|
|
nseg.epgeominfo[1].dist = param.Get(seg[1]);
|
|
mesh.AddSegment (nseg);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
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 double x, const double y, const double reffac, const bool hpref)
|
|
{
|
|
geompoints.Append (GeomPoint<D>(x, y, reffac));
|
|
geompoints.Last().hpref = hpref;
|
|
}
|
|
*/
|
|
|
|
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<int D>
|
|
void SplineGeometry<D> :: 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)
|
|
{
|
|
spline -> leftdom = leftdomain;
|
|
spline -> rightdom = rightdomain;
|
|
spline -> bc = (bc >= 0) ? bc : (splines.Size()+1);
|
|
spline -> reffak = reffac;
|
|
spline -> hpref_left = hprefleft;
|
|
spline -> hpref_right = hprefright;
|
|
spline -> copyfrom = copyfrom;
|
|
|
|
splines.Append(spline);
|
|
}
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: AppendLineSegment (const int n1, const int n2, const int leftdomain, const int rightdomain,
|
|
const int bc,
|
|
const double reffac, const bool hprefleft, const bool hprefright,
|
|
const int copyfrom)
|
|
{
|
|
SplineSeg<D> * spline = new LineSeg<D>(geompoints[n1],geompoints[n2]);
|
|
AppendSegment(spline,leftdomain,rightdomain,bc,reffac,hprefleft,hprefright,copyfrom);
|
|
}
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: AppendSplineSegment (const int n1, const int n2, const int n3, const int leftdomain, const int rightdomain,
|
|
const int bc,
|
|
const double reffac, const bool hprefleft, const bool hprefright,
|
|
const int copyfrom)
|
|
{
|
|
SplineSeg<D> * spline = new SplineSeg3<D>(geompoints[n1],geompoints[n2],geompoints[n3]);
|
|
AppendSegment(spline,leftdomain,rightdomain,bc,reffac,hprefleft,hprefright,copyfrom);
|
|
}
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: AppendCircleSegment (const int n1, const int n2, const int n3, const int leftdomain, const int rightdomain,
|
|
const int bc,
|
|
const double reffac, const bool hprefleft, const bool hprefright,
|
|
const int copyfrom)
|
|
{
|
|
SplineSeg<D> * spline = new CircleSeg<D>(geompoints[n1],geompoints[n2],geompoints[n3]);
|
|
AppendSegment(spline,leftdomain,rightdomain,bc,reffac,hprefleft,hprefright,copyfrom);
|
|
}
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: AppendDiscretePointsSegment (const Array< Point<D> > & points, const int leftdomain, const int rightdomain,
|
|
const int bc,
|
|
const double reffac, const bool hprefleft, const bool hprefright,
|
|
const int copyfrom)
|
|
{
|
|
SplineSeg<D> * spline = new DiscretePointsSeg<D>(points);
|
|
AppendSegment(spline,leftdomain,rightdomain,bc,reffac,hprefleft,hprefright,copyfrom);
|
|
}
|
|
|
|
|
|
template<int D>
|
|
void SplineGeometry<D> :: GetMaterial( const int domnr, char* & material )
|
|
{
|
|
if ( materials.Size() >= domnr)
|
|
material = materials[domnr-1];
|
|
else material = 0;
|
|
return;
|
|
}
|
|
|
|
|
|
|
|
template<int D>
|
|
double SplineGeometry<D> :: GetDomainMaxh( const int domnr )
|
|
{
|
|
if ( maxh.Size() >= domnr && domnr > 0)
|
|
return maxh[domnr-1];
|
|
else
|
|
return -1;
|
|
}
|
|
|
|
|
|
template<int D>
|
|
string SplineGeometry<D> :: GetBCName( const int bcnr ) const
|
|
{
|
|
if ( bcnames.Size() >= bcnr)
|
|
if (bcnames[bcnr-1] )
|
|
return *bcnames[bcnr-1];
|
|
return "default";
|
|
}
|
|
|
|
template<int D>
|
|
string * SplineGeometry<D> :: BCNamePtr( const int bcnr )
|
|
{
|
|
if ( bcnr > bcnames.Size() )
|
|
return 0;
|
|
else
|
|
return bcnames[bcnr-1];
|
|
}
|
|
|
|
SplineGeometry2d :: ~SplineGeometry2d()
|
|
{
|
|
;
|
|
}
|
|
|
|
|
|
int SplineGeometry2d :: GenerateMesh (Mesh*& mesh,
|
|
int perfstepsstart, int perfstepsend, char* optstring)
|
|
{
|
|
cout << "SplineGeometry2d::GenerateMesh : only a dummy" << endl;
|
|
return 0;
|
|
}
|
|
|
|
|
|
const Refinement & SplineGeometry2d :: GetRefinement () const
|
|
{
|
|
return * new Refinement2d (*this);
|
|
}
|
|
|
|
|
|
|
|
template class SplineGeometry<2>;
|
|
template class SplineGeometry<3>;
|
|
}
|
|
|
|
|