optimize CalcPartition() and PartitionBoundary()

This commit is contained in:
Matthias Hochsteger 2020-10-16 14:22:01 +02:00
parent 1a051ec555
commit 95b7720efd

View File

@ -18,29 +18,43 @@ namespace netgen
{ {
double fperel, oldf, f; double fperel, oldf, f;
int n = 10000; int n = 1;
NgArray<Point<2> > xi;
NgArray<Point<2> > xi(n); NgArray<double> hi;
NgArray<double> hi(n);
bool not_fine_enough = true;
for (int i = 0; i < n; i++)
while(not_fine_enough && n < 10000)
{
not_fine_enough = false;
n*=2;
xi.SetSize(n);
hi.SetSize(n);
for (int i = 0; i < n; i++)
{ {
xi[i] = spline.GetPoint ( (i+0.5) / n ); xi[i] = spline.GetPoint ( (i+0.5) / n );
hi[i] = mesh.GetH (Point<3> (xi[i](0), xi[i](1), 0)); hi[i] = mesh.GetH (Point<3> (xi[i](0), xi[i](1), 0));
} }
// limit slope // limit slope
double gradh = min(1/elto0,mp.grading); double gradh = min(1/elto0,mp.grading);
for (int i = 0; i < n-1; i++) for (int i = 0; i < n-1; i++)
{ {
double hnext = hi[i] + gradh * (xi[i+1]-xi[i]).Length(); double hnext = hi[i] + gradh * (xi[i+1]-xi[i]).Length();
hi[i+1] = min(hi[i+1], hnext); if(hnext > 2*hi[i])
not_fine_enough = true;
hi[i+1] = min(hi[i+1], hnext);
} }
for (int i = n-1; i > 1; i--) for (int i = n-1; i > 1; i--)
{ {
double hnext = hi[i] + gradh * (xi[i-1]-xi[i]).Length(); double hnext = hi[i] + gradh * (xi[i-1]-xi[i]).Length();
hi[i-1] = min(hi[i-1], hnext); if(hnext > 2*hi[i])
not_fine_enough = true;
hi[i-1] = min(hi[i-1], hnext);
} }
}
points.SetSize (0); points.SetSize (0);
@ -227,6 +241,10 @@ namespace netgen
mesh2d.RestrictLocalHLine (Point<3>(p1(0),p1(1),0), mesh2d.RestrictLocalHLine (Point<3>(p1(0),p1(1),0),
Point<3>(p2(0),p2(1),0), len/mp.segmentsperedge); Point<3>(p2(0),p2(1),0), len/mp.segmentsperedge);
// skip curvature restrictions for straight lines
if(spline.MaxCurvature()==0)
continue;
double hcurve = min (spline.hmax, h/spline.reffak); double hcurve = min (spline.hmax, h/spline.reffak);
double hl = GetDomainMaxh (spline.leftdom); double hl = GetDomainMaxh (spline.leftdom);
if (hl > 0) hcurve = min2 (hcurve, hl); if (hl > 0) hcurve = min2 (hcurve, hl);