1 #include <mystdlib.h>
2 #include <linalg.hpp>
3
4 namespace netgen
5 {
6
7 QuadraticPolynomial1V ::
QuadraticPolynomial1V(double ac,double acx,double acxx)8 QuadraticPolynomial1V (double ac, double acx, double acxx)
9 {
10 c = ac;
11 cx = acx;
12 cxx = acxx;
13 }
14
15 double QuadraticPolynomial1V ::
Value(double x)16 Value (double x)
17 {
18 return c + cx * x + cxx * x * x;
19 }
20
MaxUnitInterval()21 double QuadraticPolynomial1V :: MaxUnitInterval ()
22 {
23 // inner max
24 if (cxx < 0 && cx > 0 && cx < -2 * cxx)
25 {
26 return c - 0.25 * cx * cx / cxx;
27 }
28
29
30 if (cx + cxx > 0) // right edge
31 return c + cx + cxx;
32
33 // left end
34 return c;
35 }
36
37
38
39
40 LinearPolynomial2V ::
LinearPolynomial2V(double ac,double acx,double acy)41 LinearPolynomial2V (double ac, double acx, double acy)
42 {
43 c = ac;
44 cx = acx;
45 cy = acy;
46 };
47
48
49 QuadraticPolynomial2V ::
QuadraticPolynomial2V()50 QuadraticPolynomial2V ()
51 {
52 ;
53 }
54
55
56 QuadraticPolynomial2V ::
QuadraticPolynomial2V(double ac,double acx,double acy,double acxx,double acxy,double acyy)57 QuadraticPolynomial2V (double ac, double acx, double acy,
58 double acxx, double acxy, double acyy)
59 {
60 c = ac;
61 cx = acx;
62 cy = acy;
63 cxx = acxx;
64 cxy = acxy;
65 cyy = acyy;
66 }
67
68 void QuadraticPolynomial2V ::
Square(const LinearPolynomial2V & lp)69 Square (const LinearPolynomial2V & lp)
70 {
71 c = lp.c * lp.c;
72 cx = 2 * lp.c * lp.cx;
73 cy = 2 * lp.c * lp.cy;
74
75 cxx = lp.cx * lp.cx;
76 cxy = 2 * lp.cx * lp.cy;
77 cyy = lp.cy * lp.cy;
78 }
79
80 void QuadraticPolynomial2V ::
Add(double lam,const QuadraticPolynomial2V & qp2)81 Add (double lam, const QuadraticPolynomial2V & qp2)
82 {
83 c += lam * qp2.c;
84 cx += lam * qp2.cx;
85 cy += lam * qp2.cy;
86 cxx += lam * qp2.cxx;
87 cxy += lam * qp2.cxy;
88 cyy += lam * qp2.cyy;
89 }
90
91 double QuadraticPolynomial2V ::
Value(double x,double y)92 Value (double x, double y)
93 {
94 return c + cx * x + cy * y + cxx * x * x + cxy * x * y + cyy * y * y;
95 }
96
97 /*
98 double QuadraticPolynomial2V ::
99 MinUnitSquare ()
100 {
101 double x, y;
102 double minv = 1e8;
103 double val;
104 for (x = 0; x <= 1; x += 0.1)
105 for (y = 0; y <= 1; y += 0.1)
106 {
107 val = Value (x, y);
108 if (val < minv)
109 minv = val;
110 }
111 return minv;
112 };
113 */
114
115 double QuadraticPolynomial2V ::
MaxUnitSquare()116 MaxUnitSquare ()
117 {
118 // find critical point
119
120 double maxv = c;
121 double hv;
122
123 double det, x0, y0;
124 det = 4 * cxx * cyy - cxy * cxy;
125
126 if (det > 0)
127 {
128 // definite surface
129
130 x0 = (-2 * cyy * cx + cxy * cy) / det;
131 y0 = (cxy * cx -2 * cxx * cy) / det;
132
133 if (x0 >= 0 && x0 <= 1 && y0 >= 0 && y0 <= 1)
134 {
135 hv = Value (x0, y0);
136 if (hv > maxv) maxv = hv;
137 }
138 }
139
140 QuadraticPolynomial1V e1(c, cx, cxx);
141 QuadraticPolynomial1V e2(c, cy, cyy);
142 QuadraticPolynomial1V e3(c+cy+cyy, cx+cxy, cxx);
143 QuadraticPolynomial1V e4(c+cx+cxx, cy+cxy, cyy);
144
145 hv = e1.MaxUnitInterval();
146 if (hv > maxv) maxv = hv;
147 hv = e2.MaxUnitInterval();
148 if (hv > maxv) maxv = hv;
149 hv = e3.MaxUnitInterval();
150 if (hv > maxv) maxv = hv;
151 hv = e4.MaxUnitInterval();
152 if (hv > maxv) maxv = hv;
153
154 return maxv;
155 };
156
157
158
159
160 double QuadraticPolynomial2V ::
MaxUnitTriangle()161 MaxUnitTriangle ()
162 {
163 // find critical point
164
165 double maxv = c;
166 double hv;
167
168 double det, x0, y0;
169 det = 4 * cxx * cyy - cxy * cxy;
170
171 if (cxx < 0 && det > 0)
172 {
173 // definite surface
174
175 x0 = (-2 * cyy * cx + cxy * cy) / det;
176 y0 = (cxy * cx -2 * cxx * cy) / det;
177
178 if (x0 >= 0 && y0 >= 0 && x0+y0 <= 1)
179 {
180 return Value (x0, y0);
181 }
182 }
183
184
185 QuadraticPolynomial1V e1(c, cx, cxx);
186 QuadraticPolynomial1V e2(c, cy, cyy);
187 QuadraticPolynomial1V e3(c+cy+cyy, cx-cy+cxy-2*cyy, cxx-cxy+cyy);
188
189 hv = e1.MaxUnitInterval();
190 if (hv > maxv) maxv = hv;
191 hv = e2.MaxUnitInterval();
192 if (hv > maxv) maxv = hv;
193 hv = e3.MaxUnitInterval();
194 if (hv > maxv) maxv = hv;
195
196 return maxv;
197 }
198 }
199