netgen/libsrc/linalg/polynomial.cpp

199 lines
3.3 KiB
C++
Raw Permalink Normal View History

2009-01-13 04:40:13 +05:00
#include <mystdlib.h>
#include <linalg.hpp>
namespace netgen
{
QuadraticPolynomial1V ::
QuadraticPolynomial1V (double ac, double acx, double acxx)
{
c = ac;
cx = acx;
cxx = acxx;
}
double QuadraticPolynomial1V ::
Value (double x)
{
return c + cx * x + cxx * x * x;
}
double QuadraticPolynomial1V :: MaxUnitInterval ()
{
// inner max
if (cxx < 0 && cx > 0 && cx < -2 * cxx)
{
return c - 0.25 * cx * cx / cxx;
}
if (cx + cxx > 0) // right edge
return c + cx + cxx;
// left end
return c;
}
LinearPolynomial2V ::
LinearPolynomial2V (double ac, double acx, double acy)
{
c = ac;
cx = acx;
cy = acy;
};
QuadraticPolynomial2V ::
QuadraticPolynomial2V ()
{
;
}
QuadraticPolynomial2V ::
QuadraticPolynomial2V (double ac, double acx, double acy,
double acxx, double acxy, double acyy)
{
c = ac;
cx = acx;
cy = acy;
cxx = acxx;
cxy = acxy;
cyy = acyy;
}
void QuadraticPolynomial2V ::
Square (const LinearPolynomial2V & lp)
{
c = lp.c * lp.c;
cx = 2 * lp.c * lp.cx;
cy = 2 * lp.c * lp.cy;
cxx = lp.cx * lp.cx;
cxy = 2 * lp.cx * lp.cy;
cyy = lp.cy * lp.cy;
}
void QuadraticPolynomial2V ::
Add (double lam, const QuadraticPolynomial2V & qp2)
{
c += lam * qp2.c;
cx += lam * qp2.cx;
cy += lam * qp2.cy;
cxx += lam * qp2.cxx;
cxy += lam * qp2.cxy;
cyy += lam * qp2.cyy;
}
double QuadraticPolynomial2V ::
Value (double x, double y)
{
return c + cx * x + cy * y + cxx * x * x + cxy * x * y + cyy * y * y;
}
/*
double QuadraticPolynomial2V ::
MinUnitSquare ()
{
double x, y;
double minv = 1e8;
double val;
for (x = 0; x <= 1; x += 0.1)
for (y = 0; y <= 1; y += 0.1)
{
val = Value (x, y);
if (val < minv)
minv = val;
}
return minv;
};
*/
double QuadraticPolynomial2V ::
MaxUnitSquare ()
{
// find critical point
double maxv = c;
double hv;
double det, x0, y0;
det = 4 * cxx * cyy - cxy * cxy;
if (det > 0)
{
// definite surface
x0 = (-2 * cyy * cx + cxy * cy) / det;
y0 = (cxy * cx -2 * cxx * cy) / det;
if (x0 >= 0 && x0 <= 1 && y0 >= 0 && y0 <= 1)
{
hv = Value (x0, y0);
if (hv > maxv) maxv = hv;
}
}
QuadraticPolynomial1V e1(c, cx, cxx);
QuadraticPolynomial1V e2(c, cy, cyy);
QuadraticPolynomial1V e3(c+cy+cyy, cx+cxy, cxx);
QuadraticPolynomial1V e4(c+cx+cxx, cy+cxy, cyy);
hv = e1.MaxUnitInterval();
if (hv > maxv) maxv = hv;
hv = e2.MaxUnitInterval();
if (hv > maxv) maxv = hv;
hv = e3.MaxUnitInterval();
if (hv > maxv) maxv = hv;
hv = e4.MaxUnitInterval();
if (hv > maxv) maxv = hv;
return maxv;
};
double QuadraticPolynomial2V ::
MaxUnitTriangle ()
{
// find critical point
double maxv = c;
double hv;
double det, x0, y0;
det = 4 * cxx * cyy - cxy * cxy;
if (cxx < 0 && det > 0)
{
// definite surface
x0 = (-2 * cyy * cx + cxy * cy) / det;
y0 = (cxy * cx -2 * cxx * cy) / det;
if (x0 >= 0 && y0 >= 0 && x0+y0 <= 1)
{
return Value (x0, y0);
}
}
QuadraticPolynomial1V e1(c, cx, cxx);
QuadraticPolynomial1V e2(c, cy, cyy);
QuadraticPolynomial1V e3(c+cy+cyy, cx-cy+cxy-2*cyy, cxx-cxy+cyy);
hv = e1.MaxUnitInterval();
if (hv > maxv) maxv = hv;
hv = e2.MaxUnitInterval();
if (hv > maxv) maxv = hv;
hv = e3.MaxUnitInterval();
if (hv > maxv) maxv = hv;
return maxv;
}
}