diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 000e5ad6..2849f25e 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -9,12 +9,14 @@ namespace netgen - void CalcPartition (double l, double h, double h1, double h2, - double hcurve, double elto0, Array & points); + void CalcPartition (const SplineSegExt & spline, double l, + MeshingParameters & mp, double h, double h1, double h2, + double hcurve, double elto0, + Array & points); // partitionizes spline curve void Partition (const SplineSegExt & spline, - double h, double elto0, + MeshingParameters & mp, double h, double elto0, Mesh & mesh, Point3dTree & searchtree, int segnr) { enum { D = 2 }; @@ -31,7 +33,7 @@ namespace netgen double h2 = min (spline.EndPI().hmax, h/spline.EndPI().refatpoint); double hcurve = min (spline.hmax, h/spline.reffak); - CalcPartition (l, h, h1, h2, hcurve, elto0, curvepoints); + CalcPartition (spline, l, mp, h, h1, h2, hcurve, elto0, curvepoints); // cout << "curvepoints = " << curvepoints << endl; dt = 1.0 / n; @@ -130,7 +132,7 @@ namespace netgen - void SplineGeometry2d :: PartitionBoundary (double h, Mesh & mesh2d) + void SplineGeometry2d :: PartitionBoundary (MeshingParameters & mp, double h, Mesh & mesh2d) { enum { D = 2 }; Box bbox; @@ -166,13 +168,13 @@ namespace netgen maximum = min2 ( maximum, h); if ( minimum > 0 ) // GetSpline(i).Partition(minimum, elto0, mesh2d, searchtree, i+1); - Partition(GetSpline(i), minimum, elto0, mesh2d, searchtree, i+1); + Partition(GetSpline(i), mp, minimum, elto0, mesh2d, searchtree, i+1); else if ( maximum > 0 ) // GetSpline(i).Partition(maximum, elto0, mesh2d, searchtree, i+1); - Partition(GetSpline(i), maximum, elto0, mesh2d, searchtree, i+1); + Partition(GetSpline(i), mp, maximum, elto0, mesh2d, searchtree, i+1); else // GetSpline(i).Partition(h, elto0, mesh2d, searchtree, i+1); - Partition(GetSpline(i), h, elto0, mesh2d, searchtree, i+1); + Partition(GetSpline(i), mp, h, elto0, mesh2d, searchtree, i+1); } else { @@ -296,7 +298,7 @@ namespace netgen mesh = new Mesh; mesh->SetDimension (2); - geometry.PartitionBoundary (h, *mesh); + geometry.PartitionBoundary (mp, h, *mesh); // marks mesh points for hp-refinement diff --git a/libsrc/geom2d/geometry2d.cpp b/libsrc/geom2d/geometry2d.cpp index e8b10c5e..4e7baf49 100644 --- a/libsrc/geom2d/geometry2d.cpp +++ b/libsrc/geom2d/geometry2d.cpp @@ -827,47 +827,52 @@ namespace netgen - - - - void CalcPartition (double l, double h, double h1, double h2, + /* + void CalcPartition (const SplineSegExt & spline, + 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; + double fperel, oldf, f; - n = 1000; + int n = 1000; points.SetSize (0); - sum = 0; - dt = l / n; - t = 0.5 * dt; - for (i = 1; i <= n; i++) + double dt = l / n; + + double sum = 0; + for (int i = 1; i <= n; i++) { - fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2); + double t = (i-0.5)*dt; + double fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2); + double curv = spline.CalcCurvature (t/l); + cout << "curv = " << curv << endl; + if (curv < 1e-10) curv = 1e-10; + fun = min2 (fun, 0.1/curv); sum += dt / fun; - t += dt; } - nel = int (sum+1); + int nel = int (sum+1); fperel = sum / nel; points.Append (0); - i = 1; + int i = 1; oldf = 0; - t = 0.5 * dt; - for (j = 1; j <= n && i < nel; j++) + // t = 0.5 * dt; + for (int j = 1; j <= n && i < nel; j++) { - fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2); + double t = (j-0.5)*dt; + double fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2); + double curv = spline.CalcCurvature (t/l); + if (curv < 1e-10) curv = 1e-10; + fun = min2 (fun, 0.1/curv); f = oldf + dt / fun; - while (f > i * fperel && i < nel) + while (i * fperel < f && i < nel) { - points.Append ( (l/n) * (j-1 + (i * fperel - oldf) / (f - oldf)) ); + points.Append ( dt * (j-1) + (i * fperel - oldf) * fun); i++; } oldf = f; @@ -875,6 +880,104 @@ namespace netgen } points.Append (l); } + */ + + + + void CalcPartition (const SplineSegExt & spline, double l, + MeshingParameters & mp, double h, double h1, double h2, + double hcurve, double elto0, Array & points) + { + double fperel, oldf, f; + + int n = 10000; + + Array > xi(n); + Array hi(n); + + hi = h; + for (int i = 0; i < n; i++) + { + xi[i] = spline.GetPoint ( (i+0.5) / n ); + double hc = 1.0/mp.curvaturesafety / (1e-20+spline.CalcCurvature ( (i+0.5) / n )); + // *testout << "t = " << (i+0.5)/n << ", xi = " << xi[i] << ", hi = " << hc << endl; + hi[i] = min2(hi[i], hc); + } + hi[0] = min (hi[0], h1); + hi[n-1] = min (hi[n-1], h2); + + if ( (spline.GetPoint(0)-spline.GetPoint(1)).Length() < 1e-12 ) + hi[0] = hi[n-1] = min2 (hi[0], hi[n-1]); + + // limit slope + double gradh = 1/elto0; + for (int i = 0; i < n-1; i++) + { + double hnext = hi[i] + gradh * (xi[i+1]-xi[i]).Length(); + hi[i+1] = min(hi[i+1], hnext); + } + for (int i = n-1; i > 1; i--) + { + double hnext = hi[i] + gradh * (xi[i-1]-xi[i]).Length(); + hi[i-1] = min(hi[i-1], hnext); + } + + + points.SetSize (0); + + double dt = l / n; + + double sum = 0; + for (int i = 1; i <= n; i++) + { + double t = (i-0.5)*dt; + /* + double fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2); + double curv = spline.CalcCurvature (t/l); + cout << "curv = " << curv << endl; + if (curv < 1e-10) curv = 1e-10; + fun = min2 (fun, 0.1/curv); + */ + double fun = hi[i-1]; + sum += dt / fun; + } + + int nel = int (sum+1); + fperel = sum / nel; + + points.Append (0); + + int i = 1; + oldf = 0; + // t = 0.5 * dt; + for (int j = 1; j <= n && i < nel; j++) + { + double t = (j-0.5)*dt; + double fun = hi[j-1]; + /* + double fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2); + double curv = spline.CalcCurvature (t/l); + if (curv < 1e-10) curv = 1e-10; + fun = min2 (fun, 0.1/curv); + */ + f = oldf + dt / fun; + + while (i * fperel < f && i < nel) + { + points.Append ( dt * (j-1) + (i * fperel - oldf) * fun); + i++; + } + oldf = f; + t += dt; + } + points.Append (l); + } + + + + + + diff --git a/libsrc/geom2d/geometry2d.hpp b/libsrc/geom2d/geometry2d.hpp index ef0bd015..26fea07e 100644 --- a/libsrc/geom2d/geometry2d.hpp +++ b/libsrc/geom2d/geometry2d.hpp @@ -97,6 +97,15 @@ namespace netgen return seg.GetType(); } + virtual double CalcCurvature (double t) const + { + Point<2> point; + Vec<2> first, second; + GetDerivatives (t, point, first, second); + double curv = fabs(first(0)*second(1)-first(1)*second(0)) / pow(first.Length(), 3); + return curv; + } + }; @@ -141,7 +150,7 @@ namespace netgen DLL_HEADER virtual int GenerateMesh (Mesh*& mesh, MeshingParameters & mparam, int perfstepsstart, int perfstepsend); - void PartitionBoundary (double h, Mesh & mesh2d); + void PartitionBoundary (MeshingParameters & mp, double h, Mesh & mesh2d); void CopyEdgeMesh (int from, int to, Mesh & mesh2d, Point3dTree & searchtree); diff --git a/libsrc/gprim/spline.hpp b/libsrc/gprim/spline.hpp index c89e15ce..10a6f9a1 100644 --- a/libsrc/gprim/spline.hpp +++ b/libsrc/gprim/spline.hpp @@ -54,11 +54,23 @@ namespace netgen virtual Point GetPoint (double t) const = 0; /// returns a (not necessarily unit-length) tangent vector for 0 <= t <= 1 virtual Vec GetTangent (const double t) const - { cerr << "GetTangent not implemented for spline base-class" << endl; Vec dummy; return dummy;} + { + cerr << "GetTangent not implemented for spline base-class" << endl; + Vec dummy; return dummy; + } + virtual void GetDerivatives (const double t, Point & point, Vec & first, - Vec & second) const {;} + Vec & second) const + { + double eps = 1e-6; + point = GetPoint (t); + Point pl = GetPoint (t-eps); + Point pr = GetPoint (t+eps); + first = 1.0/(2*eps) * (pr-pl); + second = 1.0/sqr(eps) * ( (pr-point)+(pl-point)); + } /// returns initial point on curve