diff --git a/libsrc/geom2d/splinegeometry.cpp b/libsrc/geom2d/splinegeometry.cpp deleted file mode 100644 index bb9c3104..00000000 --- a/libsrc/geom2d/splinegeometry.cpp +++ /dev/null @@ -1,1287 +0,0 @@ -/* - -2d Spline curve for Mesh generator - -*/ - - -#include -#include -#include - -namespace netgen -{ - - - - template - void SplineGeometry :: LoadDataV2 ( ifstream & infile ) - { - // new parser by Astrid Sinwel - - PrintMessage (1, "Load 2D Geometry V2"); - int nump, leftdom, rightdom; - Point x; - int hi1, hi2, hi3; - double hd; - char buf[50], ch; - int pointnr; - - string keyword; - - Array < GeomPoint > infilepoints (0); - Array 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> 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(x, hd)); - - infilepoints.Append ( GeomPoint(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 * 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(geompoints[hi1-1], - geompoints[hi2-1]); - } - else if (strcmp (buf, "3") == 0) - { // a rational spline - infile >> hi1 >> hi2 >> hi3; - spline = new SplineSeg3 (geompoints[hi1-1], - geompoints[hi2-1], - geompoints[hi3-1]); - } - else if (strcmp (buf, "4") == 0) - { // an arc - infile >> hi1 >> hi2 >> hi3; - spline = new CircleSeg (geompoints[hi1-1], - geompoints[hi2-1], - geompoints[hi3-1]); - break; - } - else if (strcmp (buf, "discretepoints") == 0) - { - int npts; - infile >> npts; - Array< Point > pts(npts); - for (int j = 0; j < npts; j++) - for(int k=0; k> pts[j](k); - - spline = new DiscretePointsSeg (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> 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 - // void SplineGeometry :: 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 - void SplineGeometry :: 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 - SplineGeometry :: ~SplineGeometry() - { - for(int i=0; i - int SplineGeometry :: Load (const Array & 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 > pts(3); - - for(int i=0; i(GeomPoint(pts[0],1), - GeomPoint(pts[1],1)); - //(*testout) << "appending line segment " - // << pts[0] << " -- " << pts[1] << endl; - } - else if (type == 3) - { - splines[i] = new SplineSeg3(GeomPoint(pts[0],1), - GeomPoint(pts[1],1), - GeomPoint(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 - void SplineGeometry :: GetRawData (Array & raw_data) const - { - raw_data.Append(D); - raw_data.Append(elto0); - - - raw_data.Append(splines.Size()); - for(int i=0; iGetRawData(raw_data); - - - } - - /* - template - void SplineGeometry :: CSGLoad (CSGScanner & scan) - { - double hd; - Point x; - int nump, numseg; - - //scan.ReadNext(); - scan >> nump >> ';'; - - hd = 1; - geompoints.SetSize(nump); - for(int i = 0; i> x(0) >> ',' >> x(1) >> ';'; - else if(D==3) - scan >> x(0) >> ',' >> x(1) >> ',' >> x(2) >> ';'; - - geompoints[i] = GeomPoint(x,hd); - } - - scan >> numseg;// >> ';'; - - splines.SetSize(numseg); - - int pnums,pnum1,pnum2,pnum3; - - - for(int i = 0; i> ';' >> pnums >> ','; - if (pnums == 2) - { - scan >> pnum1 >> ',' >> pnum2;// >> ';'; - splines[i] = new LineSeg(geompoints[pnum1-1], - geompoints[pnum2-1]); - } - else if (pnums == 3) - { - scan >> pnum1 >> ',' >> pnum2 >> ',' - >> pnum3;// >> ';'; - splines[i] = new SplineSeg3(geompoints[pnum1-1], - geompoints[pnum2-1], - geompoints[pnum3-1]); - } - else if (pnums == 4) - { - scan >> pnum1 >> ',' >> pnum2 >> ',' - >> pnum3;// >> ';'; - splines[i] = new CircleSeg(geompoints[pnum1-1], - geompoints[pnum2-1], - geompoints[pnum3-1]); - - } - - } - } - */ - - - - template - void SplineGeometry :: Load (const char * filename) - { - - ifstream infile; - Point 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 - void SplineGeometry :: LoadDataNew ( ifstream & infile ) - { - - int nump, numseg, leftdom, rightdom; - Point 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> 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(x, hd)); - geompoints[pointnr-1] = GeomPoint(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 * 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 (geompoints[hi1-1], - geompoints[hi2-1]); - } - else if (strcmp (buf, "3") == 0) - { // a rational spline - infile >> hi1 >> hi2 >> hi3; - spline = new SplineSeg3 (geompoints[hi1-1], - geompoints[hi2-1], - geompoints[hi3-1]); - } - else if (strcmp (buf, "4") == 0) - { // an arc - infile >> hi1 >> hi2 >> hi3; - spline = new CircleSeg (geompoints[hi1-1], - geompoints[hi2-1], - geompoints[hi3-1]); - // break; - } - else if (strcmp (buf, "discretepoints") == 0) - { - int npts; - infile >> npts; - Array< Point > pts(npts); - for (int j = 0; j < npts; j++) - for(int k=0; k> pts[j](k); - - spline = new DiscretePointsSeg (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> 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 - void SplineGeometry :: LoadData ( ifstream & infile ) - { - - int nump, numseg, leftdom, rightdom; - Point 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> 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(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 * 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(geompoints[hi1-1], - geompoints[hi2-1]); - } - else if (strcmp (buf, "3") == 0) - { // a rational spline - infile >> hi1 >> hi2 >> hi3; - spline = new SplineSeg3 (geompoints[hi1-1], - geompoints[hi2-1], - geompoints[hi3-1]); - } - else if (strcmp (buf, "4") == 0) - { // an arc - infile >> hi1 >> hi2 >> hi3; - spline = new CircleSeg (geompoints[hi1-1], - geompoints[hi2-1], - geompoints[hi3-1]); - // break; - } - else if (strcmp (buf, "discretepoints") == 0) - { - int npts; - infile >> npts; - Array< Point > pts(npts); - for (int j = 0; j < npts; j++) - for(int k=0; k> pts[j](k); - - spline = new DiscretePointsSeg (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","") ); - } - } - } - - - - void SplineGeometry2d :: PartitionBoundary (double h, Mesh & mesh2d) - { - enum { D = 2 }; - Box 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;jleftdom : 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 - void SplineGeometry :: CopyEdgeMesh (int from, int to, Mesh & mesh, Point3dTree & searchtree) - { - int i; - - Array mappoints (mesh.GetNP()); - Array 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 newp = splines.Get(to)->GetPoint (param.Get(i)); - Point<3> newp3; - for(int j=0; jbc; - 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 - void SplineGeometry :: GetBoundingBox (Box & box) const - { - if (!splines.Size()) - { - Point auxp = 0.; - box.Set (auxp); - return; - } - - Array > 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 - void SplineGeometry :: SetGrading (const double grading) - { elto0 = grading;} - - /* - template - void SplineGeometry :: AppendPoint (const double x, const double y, const double reffac, const bool hpref) - { - geompoints.Append (GeomPoint(x, y, reffac)); - geompoints.Last().hpref = hpref; - } - */ - - template - void SplineGeometry :: AppendPoint (const Point & p, const double reffac, const bool hpref) - { - geompoints.Append (GeomPoint(p, reffac)); - geompoints.Last().hpref = hpref; - } - - - - - template - void SplineGeometry :: AppendSegment(SplineSeg * 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 - void SplineGeometry :: 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 * spline = new LineSeg(geompoints[n1],geompoints[n2]); - AppendSegment(spline,leftdomain,rightdomain,bc,reffac,hprefleft,hprefright,copyfrom); - } - - template - void SplineGeometry :: 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 * spline = new SplineSeg3(geompoints[n1],geompoints[n2],geompoints[n3]); - AppendSegment(spline,leftdomain,rightdomain,bc,reffac,hprefleft,hprefright,copyfrom); - } - - template - void SplineGeometry :: 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 * spline = new CircleSeg(geompoints[n1],geompoints[n2],geompoints[n3]); - AppendSegment(spline,leftdomain,rightdomain,bc,reffac,hprefleft,hprefright,copyfrom); - } - - template - void SplineGeometry :: AppendDiscretePointsSegment (const Array< Point > & points, const int leftdomain, const int rightdomain, - const int bc, - const double reffac, const bool hprefleft, const bool hprefright, - const int copyfrom) - { - SplineSeg * spline = new DiscretePointsSeg(points); - AppendSegment(spline,leftdomain,rightdomain,bc,reffac,hprefleft,hprefright,copyfrom); - } - - - template - void SplineGeometry :: GetMaterial( const int domnr, char* & material ) - { - if ( materials.Size() >= domnr) - material = materials[domnr-1]; - else material = 0; - return; - } - - - - template - double SplineGeometry :: GetDomainMaxh( const int domnr ) - { - if ( maxh.Size() >= domnr && domnr > 0) - return maxh[domnr-1]; - else - return -1; - } - - - template - string SplineGeometry :: GetBCName( const int bcnr ) const - { - if ( bcnames.Size() >= bcnr) - if (bcnames[bcnr-1] ) - return *bcnames[bcnr-1]; - return "default"; - } - - template - string * SplineGeometry :: BCNamePtr( const int bcnr ) - { - if ( bcnr > bcnames.Size() ) - return 0; - else - return bcnames[bcnr-1]; - } - - SplineGeometry2d :: ~SplineGeometry2d() - { - ; - } - - - extern void MeshFromSpline2D (SplineGeometry2d & geometry, - Mesh *& mesh, - MeshingParameters & mp); - - - int SplineGeometry2d :: GenerateMesh (Mesh*& mesh, MeshingParameters & mparam, - int perfstepsstart, int perfstepsend) - { - MeshFromSpline2D (*this, mesh, mparam); - return 0; - } - - - Refinement & SplineGeometry2d :: GetRefinement () const - { - return * new Refinement2d (*this); - } - - - - template class SplineGeometry<2>; - template class SplineGeometry<3>; -} - - diff --git a/libsrc/geom2d/splinegeometry.hpp b/libsrc/geom2d/splinegeometry.hpp deleted file mode 100644 index 0a67ff89..00000000 --- a/libsrc/geom2d/splinegeometry.hpp +++ /dev/null @@ -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 > geompoints; - Array < SplineSeg* > splines; - double elto0; - Array materials; - Array bcnames; - Array maxh; - Array quadmeshing; - Array tensormeshing; - Array layer; - - void AppendSegment(SplineSeg * 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 & 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 & raw_data) const; - - void CopyEdgeMesh (int from, int to, Mesh & mesh2d, Point3dTree & searchtree); - - const Array*> & GetSplines () const - { return splines; } - - int GetNSplines (void) const { return splines.Size(); } - string GetSplineType (const int i) const { return splines[i]->GetType(); } - SplineSeg & GetSpline (const int i) {return *splines[i];} - const SplineSeg & GetSpline (const int i) const {return *splines[i];} - - void GetBoundingBox (Box & box) const; - Box GetBoundingBox () const - { Box box; GetBoundingBox (box); return box; } - - int GetNP () const { return geompoints.Size(); } - const GeomPoint & 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 & 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 > & 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 diff --git a/libsrc/geom2d/spline.cpp b/libsrc/gprim/spline.cpp similarity index 84% rename from libsrc/geom2d/spline.cpp rename to libsrc/gprim/spline.cpp index 524c3107..f04192bc 100644 --- a/libsrc/geom2d/spline.cpp +++ b/libsrc/gprim/spline.cpp @@ -6,13 +6,11 @@ Spline curve for Mesh generator #include #include -#include -#include +#include +#include "spline.hpp" namespace netgen { -#include "spline.hpp" - // just for testing (JS) template @@ -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 & 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>; + + + + + + + } diff --git a/libsrc/geom2d/spline.hpp b/libsrc/gprim/spline.hpp similarity index 79% rename from libsrc/geom2d/spline.hpp rename to libsrc/gprim/spline.hpp index 44ae3703..b2a5c367 100644 --- a/libsrc/geom2d/spline.hpp +++ b/libsrc/gprim/spline.hpp @@ -11,8 +11,6 @@ namespace netgen { - void CalcPartition (double l, double h, double h1, double h2, - double hcurve, double elto0, Array & 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 & point, Vec & first, Vec & 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 & 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 > & points); + virtual void GetPoints (int n, Array > & 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 double SplineSeg :: Length () const @@ -303,126 +273,8 @@ namespace netgen } - - // partitionizes spline curve template - void SplineSeg :: 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 p, pold, mark, oldmark; - Array 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 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 - void SplineSeg :: GetPoints (int n, Array > & points) + void SplineSeg :: GetPoints (int n, Array > & points) const { points.SetSize (n); if (n >= 2) @@ -873,19 +725,7 @@ namespace netgen return pts[segnr] + rest*Vec(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 diff --git a/libsrc/gprim/splinegeometry.cpp b/libsrc/gprim/splinegeometry.cpp new file mode 100644 index 00000000..580e61ab --- /dev/null +++ b/libsrc/gprim/splinegeometry.cpp @@ -0,0 +1,134 @@ +/* + +2d Spline curve for Mesh generator + +*/ + + +#include +#include +#include +#include "splinegeometry.hpp" + +namespace netgen +{ + + + template + SplineGeometry :: ~SplineGeometry() + { + for(int i = 0; i < splines.Size(); i++) + delete splines[i]; + } + + + template + void SplineGeometry :: GetRawData (Array & raw_data) const + { + raw_data.Append(D); + // raw_data.Append(elto0); + + raw_data.Append(splines.Size()); + for(int i=0; iGetRawData(raw_data); + } + + + + template + int SplineGeometry :: Load (const Array & 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 > pts(3); + + for(int i=0; i(GeomPoint(pts[0],1), + GeomPoint(pts[1],1)); + } + else if (type == 3) + { + splines[i] = new SplineSeg3(GeomPoint(pts[0],1), + GeomPoint(pts[1],1), + GeomPoint(pts[2],1)); + } + else + throw NgException("something wrong with spline raw data"); + + } + return pos; + } + + + + + + + + + + template + void SplineGeometry :: GetBoundingBox (Box & box) const + { + if (!splines.Size()) + { + Point auxp = 0.; + box.Set (auxp); + return; + } + + Array > 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 + void SplineGeometry :: SetGrading (const double grading) + { + elto0 = grading; + } + */ + + template + void SplineGeometry :: AppendPoint (const Point & p, const double reffac, const bool hpref) + { + geompoints.Append (GeomPoint(p, reffac)); + geompoints.Last().hpref = hpref; + } + + + + template class SplineGeometry<2>; + template class SplineGeometry<3>; +} + + diff --git a/libsrc/gprim/splinegeometry.hpp b/libsrc/gprim/splinegeometry.hpp new file mode 100644 index 00000000..f70ae794 --- /dev/null +++ b/libsrc/gprim/splinegeometry.hpp @@ -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 > geompoints; + Array < SplineSeg* > splines; + + ~SplineGeometry(); + + int Load (const Array & raw_data, const int startpos = 0); + + void GetRawData (Array & raw_data) const; + + + const Array*> & GetSplines () const + { return splines; } + + int GetNSplines (void) const { return splines.Size(); } + string GetSplineType (const int i) const { return splines[i]->GetType(); } + SplineSeg & GetSpline (const int i) {return *splines[i];} + const SplineSeg & GetSpline (const int i) const {return *splines[i];} + + void GetBoundingBox (Box & box) const; + Box GetBoundingBox () const + { Box box; GetBoundingBox (box); return box; } + + int GetNP () const { return geompoints.Size(); } + const GeomPoint & GetPoint(int i) const { return geompoints[i]; } + + // void SetGrading (const double grading); + void AppendPoint (const Point & p, const double reffac = 1., const bool hpref = false); + + + void AppendSegment(SplineSeg * spline) + { + splines.Append (spline); + } + }; + + + +} + +#endif // _FILE_SPLINEGEOMETRY