00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #ifndef __MITTELMANNBNDRYCNTRLDIRI3D_27_HPP__
00012 #define __MITTELMANNBNDRYCNTRLDIRI3D_27_HPP__
00013
00014 #include "RegisteredTNLP.hpp"
00015
00016 #ifdef HAVE_CMATH
00017 # include <cmath>
00018 #else
00019 # ifdef HAVE_MATH_H
00020 # include <math.h>
00021 # else
00022 # error "don't have header file for math"
00023 # endif
00024 #endif
00025
00026 using namespace Ipopt;
00027
00037 class MittelmannBndryCntrlDiriBase3D_27 : public RegisteredTNLP
00038 {
00039 public:
00041 MittelmannBndryCntrlDiriBase3D_27();
00042
00044 virtual ~MittelmannBndryCntrlDiriBase3D_27();
00045
00049 virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00050 Index& nnz_h_lag, IndexStyleEnum& index_style);
00051
00053 virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u,
00054 Index m, Number* g_l, Number* g_u);
00055
00057 virtual bool get_starting_point(Index n, bool init_x, Number* x,
00058 bool init_z, Number* z_L, Number* z_U,
00059 Index m, bool init_lambda,
00060 Number* lambda);
00061
00063 virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value);
00064
00066 virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);
00067
00069 virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);
00070
00075 virtual bool eval_jac_g(Index n, const Number* x, bool new_x,
00076 Index m, Index nele_jac, Index* iRow, Index *jCol,
00077 Number* values);
00078
00083 virtual bool eval_h(Index n, const Number* x, bool new_x,
00084 Number obj_factor, Index m, const Number* lambda,
00085 bool new_lambda, Index nele_hess, Index* iRow,
00086 Index* jCol, Number* values);
00087
00089
00091 virtual bool get_scaling_parameters(Number& obj_scaling,
00092 bool& use_x_scaling, Index n,
00093 Number* x_scaling,
00094 bool& use_g_scaling, Index m,
00095 Number* g_scaling);
00096
00101 virtual void finalize_solution(SolverReturn status,
00102 Index n, const Number* x, const Number* z_L, const Number* z_U,
00103 Index m, const Number* g, const Number* lambda,
00104 Number obj_valu,
00105 const IpoptData* ip_data,
00106 IpoptCalculatedQuantities* ip_cq);
00108
00109 protected:
00113 void SetBaseParameters(Index N, Number alpha, Number lb_y,
00114 Number ub_y, Number lb_u, Number ub_u,
00115 Number d_const, Number B, Number C);
00116
00120 virtual Number y_d_cont(Number x1, Number x2, Number x3) const =0;
00122
00123 private:
00135 MittelmannBndryCntrlDiriBase3D_27(const MittelmannBndryCntrlDiriBase3D_27&);
00136 MittelmannBndryCntrlDiriBase3D_27& operator=(const MittelmannBndryCntrlDiriBase3D_27&);
00138
00142 Index N_;
00144 Number h_;
00146 Number hh_;
00148 Number hhh_;
00150 Number lb_y_;
00152 Number ub_y_;
00154 Number lb_u_;
00156 Number ub_u_;
00158 Number d_const_;
00161 Number alpha_;
00163 Number* y_d_;
00165
00170 inline Index y_index(Index i, Index j, Index k) const
00171 {
00172 return k + (N_+2)*j + (N_+2)*(N_+2)*i;
00173 }
00176 inline Index pde_index(Index i, Index j, Index k) const
00177 {
00178 return (k-1) + N_*(j-1) + N_*N_*(i-1);
00179 }
00181 inline Number x1_grid(Index i) const
00182 {
00183 return h_*(Number)i;
00184 }
00186 inline Number x2_grid(Index i) const
00187 {
00188 return h_*(Number)i;
00189 }
00191 inline Number x3_grid(Index i) const
00192 {
00193 return h_*(Number)i;
00194 }
00196 inline Number PenObj(Number t) const
00197 {
00198 if (B_ == 0.) {
00199 return 0.5*t*t;
00200 }
00201 else if (t > B_) {
00202 return B_*B_/2. + C_*(t - B_);
00203 }
00204 else if (t < -B_) {
00205 return B_*B_/2. + C_*(-t - B_);
00206 }
00207 else {
00208 const Number t2 = t*t;
00209 const Number t4 = t2*t2;
00210 const Number t6 = t4*t2;
00211 return PenA_*t2 + PenB_*t4 + PenC_*t6;
00212 }
00213 }
00215 inline Number PenObj_1(Number t) const
00216 {
00217 if (B_ == 0.) {
00218 return t;
00219 }
00220 else if (t > B_) {
00221 return C_;
00222 }
00223 else if (t < -B_) {
00224 return -C_;
00225 }
00226 else {
00227 const Number t2 = t*t;
00228 const Number t3 = t*t2;
00229 const Number t5 = t3*t2;
00230 return 2.*PenA_*t + 4.*PenB_*t3 + 6.*PenC_*t5;
00231 }
00232 }
00234 inline Number PenObj_2(Number t) const
00235 {
00236 if (B_ == 0.) {
00237 return 1.;
00238 }
00239 else if (t > B_) {
00240 return 0.;
00241 }
00242 else if (t < -B_) {
00243 return 0.;
00244 }
00245 else {
00246 const Number t2 = t*t;
00247 const Number t4 = t2*t2;
00248 return 2.*PenA_ + 12.*PenB_*t2 + 30.*PenC_*t4;
00249 }
00250 }
00252
00255 Number B_;
00256 Number C_;
00257 Number PenA_;
00258 Number PenB_;
00259 Number PenC_;
00261 };
00262
00264 class MittelmannBndryCntrlDiri3D_27 : public MittelmannBndryCntrlDiriBase3D_27
00265 {
00266 public:
00267 MittelmannBndryCntrlDiri3D_27()
00268 {}
00269
00270 virtual ~MittelmannBndryCntrlDiri3D_27()
00271 {}
00272
00273 virtual bool InitializeProblem(Index N)
00274 {
00275 if (N<1) {
00276 printf("N has to be at least 1.");
00277 return false;
00278 }
00279 Number alpha = 1e-2;
00280 Number lb_y = -1e20;
00281 Number ub_y = 3.5;
00282 Number lb_u = 0.;
00283 Number ub_u = 10.;
00284 Number d_const = -20.;
00285 Number B = 0.;
00286 Number C = 0.;
00287 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, d_const, B, C);
00288 return true;
00289 }
00290 protected:
00292 virtual Number y_d_cont(Number x1, Number x2, Number x3) const
00293 {
00294 return 3. + 5.*(x1*(x1-1.)*x2*(x2-1.));
00295 }
00296 private:
00299 MittelmannBndryCntrlDiri3D_27(const MittelmannBndryCntrlDiri3D_27&);
00300 MittelmannBndryCntrlDiri3D_27& operator=(const MittelmannBndryCntrlDiri3D_27&);
00302
00303 };
00304
00307 class MittelmannBndryCntrlDiri3D_27BT : public MittelmannBndryCntrlDiriBase3D_27
00308 {
00309 public:
00310 MittelmannBndryCntrlDiri3D_27BT()
00311 {}
00312
00313 virtual ~MittelmannBndryCntrlDiri3D_27BT()
00314 {}
00315
00316 virtual bool InitializeProblem(Index N)
00317 {
00318 if (N<1) {
00319 printf("N has to be at least 1.");
00320 return false;
00321 }
00322 Number alpha = 1e-2;
00323 Number lb_y = -1e20;
00324 Number ub_y = 3.5;
00325 Number lb_u = 0.;
00326 Number ub_u = 10.;
00327 Number d_const = -20.;
00328 Number B = .25;
00329 Number C = 0.01;
00330 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, d_const, B, C);
00331 return true;
00332 }
00333 protected:
00335 virtual Number y_d_cont(Number x1, Number x2, Number x3) const
00336 {
00337 return 3. + 5.*(x1*(x1-1.)*x2*(x2-1.));
00338 }
00339 private:
00342 MittelmannBndryCntrlDiri3D_27BT(const MittelmannBndryCntrlDiri3D_27BT&);
00343 MittelmannBndryCntrlDiri3D_27BT& operator=(const MittelmannBndryCntrlDiri3D_27BT&);
00345
00346 };
00347
00348 #endif