|
ColPack
|
|
ColPack
|
00001 #include "adolc.h" 00002 00003 // This file implements the dynamic optimization of the simplified SMB 00004 // process developed in Diehl and Walther (2006). The problem discretizes 00005 // the column as a set of mixing tanks. 6 columns are used and the two 00006 // components mix at different rates. 00007 // AMPL model written by L. T. Biegler/ CMU, Jan. 9, 2007 00008 // 00009 // transferred to C by A. Walther / TU Dresden, Jun 5 2007 00010 00011 //#define nel 2 00012 //#define ndis 2 // => 00013 //#define cstr 406 00014 //#define nel 5 00015 //#define ndis 10 // => 00016 //#define cstr 4375 00017 //#define ndis 15 // => 00018 //#define cstr 12925 00019 //#define ndis 20 // => 00020 //#define cstr 8575 00021 //#define ndis 30 // => n = 12780 00022 //#define cstr 12775 00023 //#define ndis 40 // => n = 16980 00024 //#define cstr 16980 00025 //#define ndis 50 // => n = 16980 00026 //#define cstr 21175 00027 //#define ndis 60 // => n = 25380 00028 //#define cstr 25380 00029 00030 //#define nel 10 00031 //#define ndis 10 // => n = 8755 00032 //#define cstr 8750 00033 //#define ndis 15 // => n = 12955 00034 //#define cstr 25845 00035 //#define ndis 20 // => n = 17155 00036 //#define cstr 17150 00037 //#define ndis 30 // => n = 25555 00038 //#define cstr 25550 00039 //#define ndis 40 // => n = 33955 00040 //#define cstr 33950 00041 //#define ndis 50 // => n = 42355 00042 //#define cstr 84645 00043 //#define ndis 60 // => n = 50755 00044 //#define cstr 50750 00045 00046 00047 // ad2008: 00048 #define nel 5 00049 // 1 00050 //#define ndis 10 // => 00051 //#define cstr 4375 00052 // 2 00053 #define ndis 20 // => 00054 #define cstr 8575 00055 // 3 00056 //#define ndis 40 // => n = 16980 00057 //#define cstr 16975 00058 // 4 00059 //#define ndis 60 // => n = 25380 00060 //#define cstr 25375 00061 00062 //#define nel 10 00063 // 11 replaces 3 00064 //#define ndis 20 // => n = 17155 00065 //#define cstr 17150 00066 //12 replaces 4 00067 //#define ndis 30 // => n = 25555 00068 //#define cstr 25550 00069 // 5 00070 //#define ndis 60 // => n = 50755 00071 //#define cstr 50750 00072 00073 // 6 00074 //#define nel 15 00075 //#define ndis 60 // => n = 76115 00076 //#define cstr 76120 00077 00078 // 7 00079 //#define nel 20 00080 //#define ndis 60 // => n = 101505 00081 //#define cstr 101495 00082 00083 // 8 00084 //#define nel 30 00085 //#define ndis 60 // => n = 152255 00086 //#define cstr 152245 00087 00088 // 9 00089 //#define nel 40 00090 //#define ndis 80 // => n = 270205 00091 //#define cstr 270195 00092 00093 // 10 00094 //#define nel 60 00095 //#define ndis 100 // => n = 506105 00096 //#define cstr 506095 00097 00098 #define n1 1 00099 #define n2 2 00100 #define n3 2 00101 #define n4 1 00102 #define ncol 6 00103 00104 #define nex (ndis * n1) 00105 #define nfe ((n1 + n2) * ndis) 00106 #define nra ((n1 + n2 + n3) * ndis) 00107 #define nde ((n1 + n2 + n3 + n4) * ndis) 00108 00109 #define pex 0.95 00110 #define pra 0.95 00111 00112 #define kA (2*ndis) 00113 #define kB ndis 00114 00115 #define cFEA 0.1 00116 #define cFEB 0.1 00117 00118 #define qmax 2 00119 00120 00121 00122 double kk = 0; 00123 00124 double omega[3][3]; 00125 00126 double h[nel]; 00127 00128 //****************************************************************** 00129 // Initialize dimensions 00130 00131 void init_dim(int *n, int *m) 00132 { 00133 00134 // cA, cB,cAdot, cBdot cA0 cB0 m*A,m*B,m*Adot,m*Bdot m*A0,m*B0 00135 *n = 5 + 4*(nde*nel*3) + 2*(nde*nel) + 8*(nel*3) + 4*nel 00136 + 2*nel*3 + nel; 00137 // mfe,mfedot mfe0 00138 *m = cstr-5; 00139 00140 } 00141 00142 //****************************************************************** 00143 // Initialize starting point x 00144 00145 void init_startpoint(double *x, int n) 00146 { 00147 int i,j,l; 00148 int index; 00149 00150 omega[0][0] = 0.19681547722366; 00151 omega[0][1] = 0.39442431473909; 00152 omega[0][2] = 0.37640306270047; 00153 omega[1][0] =-0.06553542585020; 00154 omega[1][1] = 0.29207341166523; 00155 omega[1][2] = 0.51248582618842; 00156 omega[2][0] = 0.02377097434822; 00157 omega[2][1] =-0.04154875212600; 00158 omega[2][2] = 0.11111111111111; 00159 00160 for(i=0;i<nel;i++) 00161 h[i] = 1.0/nel; 00162 00163 x[0] = 2; 00164 x[1] = 0.5; 00165 x[2] = 0.5; 00166 x[3] = 0.5; 00167 x[4] = 1; 00168 00169 index = 5; 00170 for(l=0;l<nde;l++) 00171 for(i=0;i<nel;i++) 00172 { 00173 for(j=0;j<3;j++) 00174 { 00175 x[index] = cFEA; 00176 index++; 00177 x[index] = cFEB; 00178 index++; 00179 x[index] = 1; 00180 index++; 00181 x[index] = 1; 00182 index++; 00183 } 00184 x[index] = cFEA; 00185 index++; 00186 x[index] = cFEB; 00187 index++; 00188 } 00189 00190 for(i=index;i<n;i++) 00191 x[i] = 1; 00192 } 00193 00194 //****************************************************************** 00195 //*************** Function Evaluation ************************* 00196 //*************** Lagrange function of optimization **************** 00197 //****************************************************************** 00198 00199 double feval(double *x, int n) 00200 { 00201 double cA[nde][nel][3],cB[nde][nel][3],cAdot[nde][nel][3],cBdot[nde][nel][3]; 00202 double cA0[nde][nel],cB0[nde][nel]; 00203 double mexA[nel][3],mexB[nel][3],mexAdot[nel][3],mexBdot[nel][3]; 00204 double mexA0[nel],mexB0[nel]; 00205 double mraA[nel][3],mraB[nel][3],mraAdot[nel][3],mraBdot[nel][3]; 00206 double mraA0[nel],mraB0[nel]; 00207 double mfe[nel][3],mfedot[nel][3]; 00208 double mfe0[nel]; 00209 00210 double q1,qde,qfe,qex,time; 00211 double q2,q3,q4,qra; 00212 00213 double res; 00214 double c[cstr],lam[cstr]; 00215 00216 double sum, test; 00217 00218 int index; 00219 int i,j,l,k; 00220 00221 q1 = x[0]; qde = x[1]; qfe = x[2]; qex = x[3]; time = x[4]; 00222 00223 q2 = q1 - qex; 00224 q3 = q1 - qex + qfe; 00225 q4 = q1 - qde; 00226 qra = qde - qex + qfe; 00227 00228 index = 5; 00229 for(l=0;l<nde;l++) 00230 for(i=0;i<nel;i++) 00231 { 00232 for(j=0;j<3;j++) 00233 { 00234 cA[l][i][j] = x[index]; 00235 index++; 00236 cB[l][i][j] = x[index]; 00237 index++; 00238 cAdot[l][i][j] = x[index]; 00239 index++; 00240 cBdot[l][i][j] = x[index]; 00241 index++; 00242 } 00243 cA0[l][i] = x[index]; 00244 index++; 00245 cB0[l][i] = x[index]; 00246 index++; 00247 } 00248 00249 for(i=0;i<nel;i++) 00250 { 00251 for(j=0;j<3;j++) 00252 { 00253 mexA[i][j] = x[index]; 00254 index++; 00255 mexB[i][j] = x[index]; 00256 index++; 00257 mexAdot[i][j] = x[index]; 00258 index++; 00259 mexBdot[i][j] = x[index]; 00260 index++; 00261 mraA[i][j] = x[index]; 00262 index++; 00263 mraB[i][j] = x[index]; 00264 index++; 00265 mraAdot[i][j] = x[index]; 00266 index++; 00267 mraBdot[i][j] = x[index]; 00268 index++; 00269 mfe[i][j] = x[index]; 00270 index++; 00271 mfedot[i][j] = x[index]; 00272 index++; 00273 } 00274 mexA0[i] = x[index]; 00275 index++; 00276 mexB0[i] = x[index]; 00277 index++; 00278 mraA0[i] = x[index]; 00279 index++; 00280 mraB0[i] = x[index]; 00281 index++; 00282 mfe0[i] = x[index]; 00283 index++; 00284 } 00285 00286 // target function 00287 00288 res = -mfe[nel-1][2]/time; 00289 00290 // constraints 00291 00292 index = 0; 00293 00294 for(l=0;l<nde;l++) 00295 for(i=0;i<nel;i++) 00296 { 00297 for(j=0;j<3;j++) 00298 { 00299 sum = 0; 00300 for(k=0;k<3;k++) 00301 sum += omega[k][j]*cAdot[l][i][k]; 00302 c[index] = cA[l][i][j] - cA0[l][i]+time*h[i]*sum; 00303 index++; 00304 sum = 0; 00305 for(k=0;k<3;k++) 00306 sum += omega[k][j]*cBdot[l][i][k]; 00307 c[index] = cB[l][i][j] - cB0[l][i]+time*h[i]*sum; 00308 index++; 00309 } 00310 } 00311 00312 for(l=0;l<nde;l++) 00313 for(i=1;i<nel;i++) 00314 { 00315 sum = 0; 00316 for(j=0;j<3;j++) 00317 sum += omega[j][2]*cAdot[l][i-1][j]; 00318 c[index] = cA0[l][i] - cA0[l][i-1] + time*h[i-1]*sum; 00319 index++; 00320 sum = 0; 00321 for(j=0;j<3;j++) 00322 sum += omega[j][2]*cBdot[l][i-1][j]; 00323 c[index] = cB0[l][i] - cB0[l][i-1] + time*h[i-1]*sum; 00324 index++; 00325 } 00326 00327 for(i=0;i<nel;i++) 00328 { 00329 for(j=0;j<3;j++) 00330 { 00331 sum = 0; 00332 for(k=0;k<3;k++) 00333 sum += omega[k][j]*mexAdot[i][k]; 00334 c[index] = mexA[i][j] - mexA0[i]+time*h[i]*sum; 00335 index++; 00336 sum = 0; 00337 for(k=0;k<3;k++) 00338 sum += omega[k][j]*mexBdot[i][k]; 00339 c[index] = mexB[i][j] - mexB0[i]+time*h[i]*sum; 00340 index++; 00341 sum = 0; 00342 for(k=0;k<3;k++) 00343 sum += omega[k][j]*mraAdot[i][k]; 00344 c[index] = mraA[i][j] - mraA0[i]+time*h[i]*sum; 00345 index++; 00346 sum = 0; 00347 for(k=0;k<3;k++) 00348 sum += omega[k][j]*mraBdot[i][k]; 00349 c[index] = mraB[i][j] - mraB0[i]+time*h[i]*sum; 00350 index++; 00351 sum = 0; 00352 for(k=0;k<3;k++) 00353 sum += omega[k][j]*mfedot[i][k]; 00354 c[index] = mfe[i][j]-mfe0[i]+time*h[i]*sum; 00355 index++; 00356 00357 c[index] = cAdot[0][i][j] - kA*(q4*cA[nde-1][i][j]-q1*cA[0][i][j])* 00358 (2*(1 + kk*cA[0][ i][ j])*(1 + kk*cA[0][ i][ j]))/ (1+(1 +kk*cA[0][i][j])*(1 +kk*cA[0][i][j])); 00359 index++; 00360 c[index] = cBdot[0][i][j] - kB*(q4*cB[nde-1][i][j]-q1*cB[0][i][j])* 00361 (2*(1 + kk*cB[0][ i][ j])*(1 + kk*cB[0][ i][ j]))/ (1+(1 +kk*cB[0][i][j])*(1 +kk*cB[0][i][j])); 00362 index++; 00363 00364 c[index] = cAdot[nfe][i][j] - kA*(q2*cA[nfe-1][i][j] + qfe*cFEA - 00365 q3*cA[nfe][i][j])*(2*(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j]))/ 00366 (1+(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j])); 00367 index++; 00368 c[index] = cBdot[nfe][i][j] - kB*(q2*cB[nfe-1][i][j] + qfe*cFEB - 00369 q3*cB[nfe][i][j])*(2*(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j]))/ 00370 (1+(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j])); 00371 index++; 00372 00373 c[index] = mexAdot[i][j] - qex*cA[nex-1][ i][j]; 00374 index++; 00375 c[index] = mexBdot[i][j] - qex*cB[nex-1][ i][j]; 00376 index++; 00377 c[index] = mraAdot[i][j] - qra*cA[nra-1][ i][j]; 00378 index++; 00379 c[index] = mraBdot[i][j] - qra*cB[nra-1][ i][j]; 00380 index++; 00381 c[index] = mfedot[i][j] - qfe; 00382 index++; 00383 00384 } 00385 } 00386 00387 for(i=1;i<nel;i++) 00388 { 00389 sum = 0; 00390 for(j=0;j<3;j++) 00391 sum += omega[j][2]*mexAdot[i-1][j]; 00392 c[index] = mexA0[i] - mexA0[i-1] + time*h[i-1]*sum; 00393 index++; 00394 sum = 0; 00395 for(j=0;j<3;j++) 00396 sum += omega[j][2]*mexBdot[i-1][j]; 00397 c[index] = mexB0[i] - mexB0[i-1] + time*h[i-1]*sum; 00398 index++; 00399 sum = 0; 00400 for(j=0;j<3;j++) 00401 sum += omega[j][2]*mraAdot[i-1][j]; 00402 c[index] = mraA0[i] - mraA0[i-1] + time*h[i-1]*sum; 00403 index++; 00404 sum = 0; 00405 for(j=0;j<3;j++) 00406 sum += omega[j][2]*mraBdot[i-1][j]; 00407 c[index] = mraB0[i] - mraB0[i-1] + time*h[i-1]*sum; 00408 index++; 00409 sum = 0; 00410 for(j=0;j<3;j++) 00411 sum += omega[j][2]*mfedot[i-1][j]; 00412 c[index] = mfe0[i] - mfe0[i-1] + time*h[i-1]*sum; 00413 index++; 00414 } 00415 00416 for(l=1;l<nex-1;l++) 00417 for(i=0;i<nel;i++) 00418 { 00419 for(j=0;j<3;j++) 00420 { 00421 c[index] = cAdot[l][i][j] - kA*q1*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 00422 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 00423 index++; 00424 c[index] = cBdot[l][i][j] - kB*q1*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 00425 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 00426 index++; 00427 } 00428 } 00429 for(l=nex-1;l<nfe;l++) 00430 for(i=0;i<nel;i++) 00431 { 00432 for(j=0;j<3;j++) 00433 { 00434 c[index] = cAdot[l][i][j] - kA*q2*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 00435 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 00436 index++; 00437 c[index] = cBdot[l][i][j] - kB*q2*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 00438 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 00439 index++; 00440 } 00441 } 00442 00443 for(l=nfe+1;l<nra;l++) 00444 for(i=0;i<nel;i++) 00445 { 00446 for(j=0;j<3;j++) 00447 { 00448 c[index] = cAdot[l][i][j] - kA*q3*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 00449 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 00450 index++; 00451 c[index] = cBdot[l][i][j] - kB*q3*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 00452 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 00453 index++; 00454 } 00455 } 00456 00457 for(l=nra;l<nde;l++) 00458 { 00459 for(i=0;i<nel;i++) 00460 { 00461 for(j=0;j<3;j++) 00462 { 00463 c[index] = cAdot[l][i][j] - kA*q4*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 00464 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 00465 index++; 00466 c[index] = cBdot[l][i][j] - kB*q4*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 00467 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 00468 index++; 00469 } 00470 } 00471 c[index] = cA0[l][0] - cA[l-nra][nel-1][2]; 00472 index++; 00473 c[index] = cB0[l][0] - cB[l-nra][nel-1][2]; 00474 index++; 00475 } 00476 00477 for(l=0;l<nra;l++) 00478 { 00479 c[index] = cA0[l][0] - cA[l+ndis][nel-1][2]; 00480 index++; 00481 c[index] = cB0[l][0] - cB[l+ndis][nel-1][2]; 00482 index++; 00483 } 00484 00485 c[index] = mexA0[0]; 00486 index++; 00487 c[index] = mexB0[0]; 00488 index++; 00489 c[index] = mraA0[0]; 00490 index++; 00491 c[index] = mraB0[0]; 00492 index++; 00493 c[index] = mfe0[0]; 00494 index++; 00495 00496 //printf(" f index %d \n",index); 00497 00498 for(i=0;i<cstr;i++) 00499 lam[i] = 1; 00500 00501 for(i=0;i<cstr;i++) 00502 res += lam[i]*c[i]; 00503 00504 return res; 00505 } 00506 00507 00508 //************** active version ************************ 00509 00510 00511 adouble feval_ad(double *x, int n) 00512 { 00513 adouble cA[nde][nel][3],cB[nde][nel][3],cAdot[nde][nel][3],cBdot[nde][nel][3]; 00514 adouble cA0[nde][nel],cB0[nde][nel]; 00515 adouble mexA[nel][3],mexB[nel][3],mexAdot[nel][3],mexBdot[nel][3]; 00516 adouble mexA0[nel],mexB0[nel]; 00517 adouble mraA[nel][3],mraB[nel][3],mraAdot[nel][3],mraBdot[nel][3]; 00518 adouble mraA0[nel],mraB0[nel]; 00519 adouble mfe[nel][3],mfedot[nel][3]; 00520 adouble mfe0[nel]; 00521 00522 adouble q1,qde,qfe,qex,time; 00523 adouble q2,q3,q4,qra; 00524 00525 adouble res; 00526 adouble c[cstr]; 00527 double lam[cstr]; 00528 adouble sum; 00529 00530 int index; 00531 int i,j,l,k; 00532 00533 q1 <<= x[0]; qde <<= x[1]; qfe <<= x[2]; qex <<= x[3]; time <<= x[4]; 00534 00535 q2 = q1 - qex; 00536 q3 = q1 - qex + qfe; 00537 q4 = q1 - qde; 00538 qra = qde - qex + qfe; 00539 00540 index = 5; 00541 for(l=0;l<nde;l++) 00542 for(i=0;i<nel;i++) 00543 { 00544 for(j=0;j<3;j++) 00545 { 00546 cA[l][i][j] <<= x[index]; 00547 index++; 00548 cB[l][i][j] <<= x[index]; 00549 index++; 00550 cAdot[l][i][j] <<= x[index]; 00551 index++; 00552 cBdot[l][i][j] <<= x[index]; 00553 index++; 00554 } 00555 cA0[l][i] <<= x[index]; 00556 index++; 00557 cB0[l][i] <<= x[index]; 00558 index++; 00559 } 00560 00561 for(i=0;i<nel;i++) 00562 { 00563 for(j=0;j<3;j++) 00564 { 00565 mexA[i][j] <<= x[index]; 00566 index++; 00567 mexB[i][j] <<= x[index]; 00568 index++; 00569 mexAdot[i][j] <<= x[index]; 00570 index++; 00571 mexBdot[i][j] <<= x[index]; 00572 index++; 00573 mraA[i][j] <<= x[index]; 00574 index++; 00575 mraB[i][j] <<= x[index]; 00576 index++; 00577 mraAdot[i][j] <<= x[index]; 00578 index++; 00579 mraBdot[i][j] <<= x[index]; 00580 index++; 00581 mfe[i][j] <<= x[index]; 00582 index++; 00583 mfedot[i][j] <<= x[index]; 00584 index++; 00585 } 00586 mexA0[i] <<= x[index]; 00587 index++; 00588 mexB0[i] <<= x[index]; 00589 index++; 00590 mraA0[i] <<= x[index]; 00591 index++; 00592 mraB0[i] <<= x[index]; 00593 index++; 00594 mfe0[i] <<= x[index]; 00595 index++; 00596 } 00597 00598 res = -mfe[nel-1][2]/time; 00599 00600 index = 0; 00601 00602 for(l=0;l<nde;l++) 00603 for(i=0;i<nel;i++) 00604 { 00605 for(j=0;j<3;j++) 00606 { 00607 sum = 0; 00608 for(k=0;k<3;k++) 00609 sum += omega[k][j]*cAdot[l][i][k]; 00610 c[index] = cA[l][i][j] - cA0[l][i]+time*h[i]*sum; 00611 index++; 00612 sum = 0; 00613 for(k=0;k<3;k++) 00614 sum += omega[k][j]*cBdot[l][i][k]; 00615 c[index] = cB[l][i][j] - cB0[l][i]+time*h[i]*sum; 00616 index++; 00617 } 00618 } 00619 00620 for(l=0;l<nde;l++) 00621 for(i=1;i<nel;i++) 00622 { 00623 sum = 0; 00624 for(j=0;j<3;j++) 00625 sum += omega[j][2]*cAdot[l][i-1][j]; 00626 c[index] = cA0[l][i] - cA0[l][i-1] + time*h[i-1]*sum; 00627 index++; 00628 sum = 0; 00629 for(j=0;j<3;j++) 00630 sum += omega[j][2]*cBdot[l][i-1][j]; 00631 c[index] = cB0[l][i] - cB0[l][i-1] + time*h[i-1]*sum; 00632 index++; 00633 } 00634 00635 for(i=0;i<nel;i++) 00636 { 00637 for(j=0;j<3;j++) 00638 { 00639 sum = 0; 00640 for(k=0;k<3;k++) 00641 sum += omega[k][j]*mexAdot[i][k]; 00642 c[index] = mexA[i][j] - mexA0[i]+time*h[i]*sum; 00643 index++; 00644 sum = 0; 00645 for(k=0;k<3;k++) 00646 sum += omega[k][j]*mexBdot[i][k]; 00647 c[index] = mexB[i][j] - mexB0[i]+time*h[i]*sum; 00648 index++; 00649 sum = 0; 00650 for(k=0;k<3;k++) 00651 sum += omega[k][j]*mraAdot[i][k]; 00652 c[index] = mraA[i][j] - mraA0[i]+time*h[i]*sum; 00653 index++; 00654 sum = 0; 00655 for(k=0;k<3;k++) 00656 sum += omega[k][j]*mraBdot[i][k]; 00657 c[index] = mraB[i][j] - mraB0[i]+time*h[i]*sum; 00658 index++; 00659 sum = 0; 00660 for(k=0;k<3;k++) 00661 sum += omega[k][j]*mfedot[i][k]; 00662 c[index] = mfe[i][j]-mfe0[i]+time*h[i]*sum; 00663 index++; 00664 00665 c[index] = cAdot[0][i][j] - kA*(q4*cA[nde-1][i][j]-q1*cA[0][i][j])* 00666 (2*(1 + kk*cA[0][ i][ j])*(1 + kk*cA[0][ i][ j]))/ (1+(1 +kk*cA[0][i][j])*(1 +kk*cA[0][i][j])); 00667 index++; 00668 c[index] = cBdot[0][i][j] - kB*(q4*cB[nde-1][i][j]-q1*cB[0][i][j])* 00669 (2*(1 + kk*cB[0][ i][ j])*(1 + kk*cB[0][ i][ j]))/ (1+(1 +kk*cB[0][i][j])*(1 +kk*cB[0][i][j])); 00670 index++; 00671 00672 c[index] = cAdot[nfe][i][j] - kA*(q2*cA[nfe-1][i][j] + qfe*cFEA - 00673 q3*cA[nfe][i][j])*(2*(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j]))/ 00674 (1+(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j])); 00675 index++; 00676 c[index] = cBdot[nfe][i][j] - kB*(q2*cB[nfe-1][i][j] + qfe*cFEB - 00677 q3*cB[nfe][i][j])*(2*(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j]))/ 00678 (1+(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j])); 00679 index++; 00680 00681 c[index] = mexAdot[i][j] - qex*cA[nex-1][ i][j]; 00682 index++; 00683 c[index] = mexBdot[i][j] - qex*cB[nex-1][ i][j]; 00684 index++; 00685 c[index] = mraAdot[i][j] - qra*cA[nra-1][ i][j]; 00686 index++; 00687 c[index] = mraBdot[i][j] - qra*cB[nra-1][ i][j]; 00688 index++; 00689 c[index] = mfedot[i][j] - qfe; 00690 index++; 00691 } 00692 } 00693 00694 for(i=1;i<nel;i++) 00695 { 00696 sum = 0; 00697 for(j=0;j<3;j++) 00698 sum += omega[j][2]*mexAdot[i-1][j]; 00699 c[index] = mexA0[i] - mexA0[i-1] + time*h[i-1]*sum; 00700 index++; 00701 sum = 0; 00702 for(j=0;j<3;j++) 00703 sum += omega[j][2]*mexBdot[i-1][j]; 00704 c[index] = mexB0[i] - mexB0[i-1] + time*h[i-1]*sum; 00705 index++; 00706 sum = 0; 00707 for(j=0;j<3;j++) 00708 sum += omega[j][2]*mraAdot[i-1][j]; 00709 c[index] = mraA0[i] - mraA0[i-1] + time*h[i-1]*sum; 00710 index++; 00711 sum = 0; 00712 for(j=0;j<3;j++) 00713 sum += omega[j][2]*mraBdot[i-1][j]; 00714 c[index] = mraB0[i] - mraB0[i-1] + time*h[i-1]*sum; 00715 index++; 00716 sum = 0; 00717 for(j=0;j<3;j++) 00718 sum += omega[j][2]*mfedot[i-1][j]; 00719 c[index] = mfe0[i] - mfe0[i-1] + time*h[i-1]*sum; 00720 index++; 00721 } 00722 00723 for(l=1;l<nex-1;l++) 00724 for(i=0;i<nel;i++) 00725 { 00726 for(j=0;j<3;j++) 00727 { 00728 c[index] = cAdot[l][i][j] - kA*q1*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 00729 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 00730 index++; 00731 c[index] = cBdot[l][i][j] - kB*q1*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 00732 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 00733 index++; 00734 } 00735 } 00736 00737 for(l=nex-1;l<nfe;l++) 00738 for(i=0;i<nel;i++) 00739 { 00740 for(j=0;j<3;j++) 00741 { 00742 c[index] = cAdot[l][i][j] - kA*q2*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 00743 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 00744 index++; 00745 c[index] = cBdot[l][i][j] - kB*q2*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 00746 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 00747 index++; 00748 } 00749 } 00750 00751 for(l=nfe+1;l<nra;l++) 00752 for(i=0;i<nel;i++) 00753 { 00754 for(j=0;j<3;j++) 00755 { 00756 c[index] = cAdot[l][i][j] - kA*q3*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 00757 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 00758 index++; 00759 c[index] = cBdot[l][i][j] - kB*q3*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 00760 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 00761 index++; 00762 } 00763 } 00764 00765 for(l=nra;l<nde;l++) 00766 { 00767 for(i=0;i<nel;i++) 00768 { 00769 for(j=0;j<3;j++) 00770 { 00771 c[index] = cAdot[l][i][j] - kA*q4*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 00772 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 00773 index++; 00774 c[index] = cBdot[l][i][j] - kB*q4*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 00775 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 00776 index++; 00777 } 00778 } 00779 c[index] = cA0[l][0] - cA[l-nra][nel-1][2]; 00780 index++; 00781 c[index] = cB0[l][0] - cB[l-nra][nel-1][2]; 00782 index++; 00783 } 00784 00785 for(l=0;l<nra;l++) 00786 { 00787 c[index] = cA0[l][0] - cA[l+ndis][nel-1][2]; 00788 index++; 00789 c[index] = cB0[l][0] - cB[l+ndis][nel-1][2]; 00790 index++; 00791 } 00792 00793 c[index] = mexA0[0]; 00794 index++; 00795 c[index] = mexB0[0]; 00796 index++; 00797 c[index] = mraA0[0]; 00798 index++; 00799 c[index] = mraB0[0]; 00800 index++; 00801 c[index] = mfe0[0]; 00802 index++; 00803 00804 00805 for(i=0;i<cstr;i++) 00806 { 00807 lam[i] = 1; 00808 } 00809 00810 for(i=0;i<cstr;i++) 00811 { 00812 res += lam[i]*c[i]; 00813 } 00814 00815 return res; 00816 00817 } 00818 00819 //****************************************************************** 00820 //*************** Constraint Evaluation *********************** 00821 //*************** constraints of optimization **************** 00822 //****************************************************************** 00823 00824 void ceval(double *x, double *c, int n) 00825 { 00826 double cA[nde][nel][3],cB[nde][nel][3],cAdot[nde][nel][3],cBdot[nde][nel][3]; 00827 double cA0[nde][nel],cB0[nde][nel]; 00828 double mexA[nel][3],mexB[nel][3],mexAdot[nel][3],mexBdot[nel][3]; 00829 double mexA0[nel],mexB0[nel]; 00830 double mraA[nel][3],mraB[nel][3],mraAdot[nel][3],mraBdot[nel][3]; 00831 double mraA0[nel],mraB0[nel]; 00832 double mfe[nel][3],mfedot[nel][3]; 00833 double mfe0[nel]; 00834 00835 double q1,qde,qfe,qex,time; 00836 double q2,q3,q4,qra; 00837 00838 double sum, test; 00839 00840 int index; 00841 int i,j,l,k; 00842 00843 q1 = x[0]; qde = x[1]; qfe = x[2]; qex = x[3]; time = x[4]; 00844 00845 q2 = q1 - qex; 00846 q3 = q1 - qex + qfe; 00847 q4 = q1 - qde; 00848 qra = qde - qex + qfe; 00849 00850 index = 5; 00851 for(l=0;l<nde;l++) 00852 for(i=0;i<nel;i++) 00853 { 00854 for(j=0;j<3;j++) 00855 { 00856 cA[l][i][j] = x[index]; 00857 index++; 00858 cB[l][i][j] = x[index]; 00859 index++; 00860 cAdot[l][i][j] = x[index]; 00861 index++; 00862 cBdot[l][i][j] = x[index]; 00863 index++; 00864 } 00865 cA0[l][i] = x[index]; 00866 index++; 00867 cB0[l][i] = x[index]; 00868 index++; 00869 } 00870 00871 for(i=0;i<nel;i++) 00872 { 00873 for(j=0;j<3;j++) 00874 { 00875 mexA[i][j] = x[index]; 00876 index++; 00877 mexB[i][j] = x[index]; 00878 index++; 00879 mexAdot[i][j] = x[index]; 00880 index++; 00881 mexBdot[i][j] = x[index]; 00882 index++; 00883 mraA[i][j] = x[index]; 00884 index++; 00885 mraB[i][j] = x[index]; 00886 index++; 00887 mraAdot[i][j] = x[index]; 00888 index++; 00889 mraBdot[i][j] = x[index]; 00890 index++; 00891 mfe[i][j] = x[index]; 00892 index++; 00893 mfedot[i][j] = x[index]; 00894 index++; 00895 } 00896 mexA0[i] = x[index]; 00897 index++; 00898 mexB0[i] = x[index]; 00899 index++; 00900 mraA0[i] = x[index]; 00901 index++; 00902 mraB0[i] = x[index]; 00903 index++; 00904 mfe0[i] = x[index]; 00905 index++; 00906 } 00907 00908 // constraints 00909 00910 index = 0; 00911 00912 for(l=0;l<nde;l++) 00913 for(i=0;i<nel;i++) 00914 { 00915 for(j=0;j<3;j++) 00916 { 00917 sum = 0; 00918 for(k=0;k<3;k++) 00919 sum += omega[k][j]*cAdot[l][i][k]; 00920 c[index] = cA[l][i][j] - cA0[l][i]+time*h[i]*sum; 00921 index++; 00922 sum = 0; 00923 for(k=0;k<3;k++) 00924 sum += omega[k][j]*cBdot[l][i][k]; 00925 c[index] = cB[l][i][j] - cB0[l][i]+time*h[i]*sum; 00926 index++; 00927 } 00928 } 00929 00930 for(l=0;l<nde;l++) 00931 for(i=1;i<nel;i++) 00932 { 00933 sum = 0; 00934 for(j=0;j<3;j++) 00935 sum += omega[j][2]*cAdot[l][i-1][j]; 00936 c[index] = cA0[l][i] - cA0[l][i-1] + time*h[i-1]*sum; 00937 index++; 00938 sum = 0; 00939 for(j=0;j<3;j++) 00940 sum += omega[j][2]*cBdot[l][i-1][j]; 00941 c[index] = cB0[l][i] - cB0[l][i-1] + time*h[i-1]*sum; 00942 index++; 00943 } 00944 00945 for(i=0;i<nel;i++) 00946 { 00947 for(j=0;j<3;j++) 00948 { 00949 sum = 0; 00950 for(k=0;k<3;k++) 00951 sum += omega[k][j]*mexAdot[i][k]; 00952 c[index] = mexA[i][j] - mexA0[i]+time*h[i]*sum; 00953 index++; 00954 sum = 0; 00955 for(k=0;k<3;k++) 00956 sum += omega[k][j]*mexBdot[i][k]; 00957 c[index] = mexB[i][j] - mexB0[i]+time*h[i]*sum; 00958 index++; 00959 sum = 0; 00960 for(k=0;k<3;k++) 00961 sum += omega[k][j]*mraAdot[i][k]; 00962 c[index] = mraA[i][j] - mraA0[i]+time*h[i]*sum; 00963 index++; 00964 sum = 0; 00965 for(k=0;k<3;k++) 00966 sum += omega[k][j]*mraBdot[i][k]; 00967 c[index] = mraB[i][j] - mraB0[i]+time*h[i]*sum; 00968 index++; 00969 sum = 0; 00970 for(k=0;k<3;k++) 00971 sum += omega[k][j]*mfedot[i][k]; 00972 c[index] = mfe[i][j]-mfe0[i]+time*h[i]*sum; 00973 index++; 00974 00975 c[index] = cAdot[0][i][j] - kA*(q4*cA[nde-1][i][j]-q1*cA[0][i][j])* 00976 (2*(1 + kk*cA[0][ i][ j])*(1 + kk*cA[0][ i][ j]))/ (1+(1 +kk*cA[0][i][j])*(1 +kk*cA[0][i][j])); 00977 index++; 00978 c[index] = cBdot[0][i][j] - kB*(q4*cB[nde-1][i][j]-q1*cB[0][i][j])* 00979 (2*(1 + kk*cB[0][ i][ j])*(1 + kk*cB[0][ i][ j]))/ (1+(1 +kk*cB[0][i][j])*(1 +kk*cB[0][i][j])); 00980 index++; 00981 00982 c[index] = cAdot[nfe][i][j] - kA*(q2*cA[nfe-1][i][j] + qfe*cFEA - 00983 q3*cA[nfe][i][j])*(2*(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j]))/ 00984 (1+(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j])); 00985 index++; 00986 c[index] = cBdot[nfe][i][j] - kB*(q2*cB[nfe-1][i][j] + qfe*cFEB - 00987 q3*cB[nfe][i][j])*(2*(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j]))/ 00988 (1+(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j])); 00989 index++; 00990 00991 c[index] = mexAdot[i][j] - qex*cA[nex-1][ i][j]; 00992 index++; 00993 c[index] = mexBdot[i][j] - qex*cB[nex-1][ i][j]; 00994 index++; 00995 c[index] = mraAdot[i][j] - qra*cA[nra-1][ i][j]; 00996 index++; 00997 c[index] = mraBdot[i][j] - qra*cB[nra-1][ i][j]; 00998 index++; 00999 c[index] = mfedot[i][j] - qfe; 01000 index++; 01001 01002 } 01003 } 01004 01005 for(i=1;i<nel;i++) 01006 { 01007 sum = 0; 01008 for(j=0;j<3;j++) 01009 sum += omega[j][2]*mexAdot[i-1][j]; 01010 c[index] = mexA0[i] - mexA0[i-1] + time*h[i-1]*sum; 01011 index++; 01012 sum = 0; 01013 for(j=0;j<3;j++) 01014 sum += omega[j][2]*mexBdot[i-1][j]; 01015 c[index] = mexB0[i] - mexB0[i-1] + time*h[i-1]*sum; 01016 index++; 01017 sum = 0; 01018 for(j=0;j<3;j++) 01019 sum += omega[j][2]*mraAdot[i-1][j]; 01020 c[index] = mraA0[i] - mraA0[i-1] + time*h[i-1]*sum; 01021 index++; 01022 sum = 0; 01023 for(j=0;j<3;j++) 01024 sum += omega[j][2]*mraBdot[i-1][j]; 01025 c[index] = mraB0[i] - mraB0[i-1] + time*h[i-1]*sum; 01026 index++; 01027 sum = 0; 01028 for(j=0;j<3;j++) 01029 sum += omega[j][2]*mfedot[i-1][j]; 01030 c[index] = mfe0[i] - mfe0[i-1] + time*h[i-1]*sum; 01031 index++; 01032 } 01033 01034 for(l=1;l<nex-1;l++) 01035 for(i=0;i<nel;i++) 01036 { 01037 for(j=0;j<3;j++) 01038 { 01039 c[index] = cAdot[l][i][j] - kA*q1*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 01040 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 01041 index++; 01042 c[index] = cBdot[l][i][j] - kB*q1*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 01043 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 01044 index++; 01045 } 01046 } 01047 for(l=nex-1;l<nfe;l++) 01048 for(i=0;i<nel;i++) 01049 { 01050 for(j=0;j<3;j++) 01051 { 01052 c[index] = cAdot[l][i][j] - kA*q2*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 01053 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 01054 index++; 01055 c[index] = cBdot[l][i][j] - kB*q2*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 01056 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 01057 index++; 01058 } 01059 } 01060 01061 for(l=nfe+1;l<nra;l++) 01062 for(i=0;i<nel;i++) 01063 { 01064 for(j=0;j<3;j++) 01065 { 01066 c[index] = cAdot[l][i][j] - kA*q3*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 01067 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 01068 index++; 01069 c[index] = cBdot[l][i][j] - kB*q3*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 01070 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 01071 index++; 01072 } 01073 } 01074 01075 for(l=nra;l<nde;l++) 01076 { 01077 for(i=0;i<nel;i++) 01078 { 01079 for(j=0;j<3;j++) 01080 { 01081 c[index] = cAdot[l][i][j] - kA*q4*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 01082 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 01083 index++; 01084 c[index] = cBdot[l][i][j] - kB*q4*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 01085 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 01086 index++; 01087 } 01088 } 01089 c[index] = cA0[l][0] - cA[l-nra][nel-1][2]; 01090 index++; 01091 c[index] = cB0[l][0] - cB[l-nra][nel-1][2]; 01092 index++; 01093 } 01094 01095 for(l=0;l<nra;l++) 01096 { 01097 c[index] = cA0[l][0] - cA[l+ndis][nel-1][2]; 01098 index++; 01099 c[index] = cB0[l][0] - cB[l+ndis][nel-1][2]; 01100 index++; 01101 } 01102 01103 01104 } 01105 01106 01107 //************** active version ************************ 01108 01109 01110 void ceval_ad(double *x, adouble *c, int n) 01111 { 01112 adouble cA[nde][nel][3],cB[nde][nel][3],cAdot[nde][nel][3],cBdot[nde][nel][3]; 01113 adouble cA0[nde][nel],cB0[nde][nel]; 01114 adouble mexA[nel][3],mexB[nel][3],mexAdot[nel][3],mexBdot[nel][3]; 01115 adouble mexA0[nel],mexB0[nel]; 01116 adouble mraA[nel][3],mraB[nel][3],mraAdot[nel][3],mraBdot[nel][3]; 01117 adouble mraA0[nel],mraB0[nel]; 01118 adouble mfe[nel][3],mfedot[nel][3]; 01119 adouble mfe0[nel]; 01120 01121 adouble q1,qde,qfe,qex,time; 01122 adouble q2,q3,q4,qra; 01123 01124 adouble sum; 01125 01126 int index; 01127 int i,j,l,k; 01128 01129 q1 <<= x[0]; qde <<= x[1]; qfe <<= x[2]; qex <<= x[3]; time <<= x[4]; 01130 01131 q2 = q1 - qex; 01132 q3 = q1 - qex + qfe; 01133 q4 = q1 - qde; 01134 qra = qde - qex + qfe; 01135 01136 index = 5; 01137 for(l=0;l<nde;l++) 01138 for(i=0;i<nel;i++) 01139 { 01140 for(j=0;j<3;j++) 01141 { 01142 cA[l][i][j] <<= x[index]; 01143 index++; 01144 cB[l][i][j] <<= x[index]; 01145 index++; 01146 cAdot[l][i][j] <<= x[index]; 01147 index++; 01148 cBdot[l][i][j] <<= x[index]; 01149 index++; 01150 } 01151 cA0[l][i] <<= x[index]; 01152 index++; 01153 cB0[l][i] <<= x[index]; 01154 index++; 01155 } 01156 01157 for(i=0;i<nel;i++) 01158 { 01159 for(j=0;j<3;j++) 01160 { 01161 mexA[i][j] <<= x[index]; 01162 index++; 01163 mexB[i][j] <<= x[index]; 01164 index++; 01165 mexAdot[i][j] <<= x[index]; 01166 index++; 01167 mexBdot[i][j] <<= x[index]; 01168 index++; 01169 mraA[i][j] <<= x[index]; 01170 index++; 01171 mraB[i][j] <<= x[index]; 01172 index++; 01173 mraAdot[i][j] <<= x[index]; 01174 index++; 01175 mraBdot[i][j] <<= x[index]; 01176 index++; 01177 mfe[i][j] <<= x[index]; 01178 index++; 01179 mfedot[i][j] <<= x[index]; 01180 index++; 01181 } 01182 mexA0[i] <<= x[index]; 01183 index++; 01184 mexB0[i] <<= x[index]; 01185 index++; 01186 mraA0[i] <<= x[index]; 01187 index++; 01188 mraB0[i] <<= x[index]; 01189 index++; 01190 mfe0[i] <<= x[index]; 01191 index++; 01192 } 01193 01194 01195 index = 0; 01196 01197 for(l=0;l<nde;l++) 01198 for(i=0;i<nel;i++) 01199 { 01200 for(j=0;j<3;j++) 01201 { 01202 sum = 0; 01203 for(k=0;k<3;k++) 01204 sum += omega[k][j]*cAdot[l][i][k]; 01205 c[index] = cA[l][i][j] - cA0[l][i]+time*h[i]*sum; 01206 index++; 01207 sum = 0; 01208 for(k=0;k<3;k++) 01209 sum += omega[k][j]*cBdot[l][i][k]; 01210 c[index] = cB[l][i][j] - cB0[l][i]+time*h[i]*sum; 01211 index++; 01212 } 01213 } 01214 01215 for(l=0;l<nde;l++) 01216 for(i=1;i<nel;i++) 01217 { 01218 sum = 0; 01219 for(j=0;j<3;j++) 01220 sum += omega[j][2]*cAdot[l][i-1][j]; 01221 c[index] = cA0[l][i] - cA0[l][i-1] + time*h[i-1]*sum; 01222 index++; 01223 sum = 0; 01224 for(j=0;j<3;j++) 01225 sum += omega[j][2]*cBdot[l][i-1][j]; 01226 c[index] = cB0[l][i] - cB0[l][i-1] + time*h[i-1]*sum; 01227 index++; 01228 } 01229 01230 for(i=0;i<nel;i++) 01231 { 01232 for(j=0;j<3;j++) 01233 { 01234 sum = 0; 01235 for(k=0;k<3;k++) 01236 sum += omega[k][j]*mexAdot[i][k]; 01237 c[index] = mexA[i][j] - mexA0[i]+time*h[i]*sum; 01238 index++; 01239 sum = 0; 01240 for(k=0;k<3;k++) 01241 sum += omega[k][j]*mexBdot[i][k]; 01242 c[index] = mexB[i][j] - mexB0[i]+time*h[i]*sum; 01243 index++; 01244 sum = 0; 01245 for(k=0;k<3;k++) 01246 sum += omega[k][j]*mraAdot[i][k]; 01247 c[index] = mraA[i][j] - mraA0[i]+time*h[i]*sum; 01248 index++; 01249 sum = 0; 01250 for(k=0;k<3;k++) 01251 sum += omega[k][j]*mraBdot[i][k]; 01252 c[index] = mraB[i][j] - mraB0[i]+time*h[i]*sum; 01253 index++; 01254 sum = 0; 01255 for(k=0;k<3;k++) 01256 sum += omega[k][j]*mfedot[i][k]; 01257 c[index] = mfe[i][j]-mfe0[i]+time*h[i]*sum; 01258 index++; 01259 01260 c[index] = cAdot[0][i][j] - kA*(q4*cA[nde-1][i][j]-q1*cA[0][i][j])* 01261 (2*(1 + kk*cA[0][ i][ j])*(1 + kk*cA[0][ i][ j]))/ (1+(1 +kk*cA[0][i][j])*(1 +kk*cA[0][i][j])); 01262 index++; 01263 c[index] = cBdot[0][i][j] - kB*(q4*cB[nde-1][i][j]-q1*cB[0][i][j])* 01264 (2*(1 + kk*cB[0][ i][ j])*(1 + kk*cB[0][ i][ j]))/ (1+(1 +kk*cB[0][i][j])*(1 +kk*cB[0][i][j])); 01265 index++; 01266 01267 c[index] = cAdot[nfe][i][j] - kA*(q2*cA[nfe-1][i][j] + qfe*cFEA - 01268 q3*cA[nfe][i][j])*(2*(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j]))/ 01269 (1+(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j])); 01270 index++; 01271 c[index] = cBdot[nfe][i][j] - kB*(q2*cB[nfe-1][i][j] + qfe*cFEB - 01272 q3*cB[nfe][i][j])*(2*(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j]))/ 01273 (1+(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j])); 01274 index++; 01275 01276 c[index] = mexAdot[i][j] - qex*cA[nex-1][ i][j]; 01277 index++; 01278 c[index] = mexBdot[i][j] - qex*cB[nex-1][ i][j]; 01279 index++; 01280 c[index] = mraAdot[i][j] - qra*cA[nra-1][ i][j]; 01281 index++; 01282 c[index] = mraBdot[i][j] - qra*cB[nra-1][ i][j]; 01283 index++; 01284 c[index] = mfedot[i][j] - qfe; 01285 index++; 01286 } 01287 } 01288 01289 for(i=1;i<nel;i++) 01290 { 01291 sum = 0; 01292 for(j=0;j<3;j++) 01293 sum += omega[j][2]*mexAdot[i-1][j]; 01294 c[index] = mexA0[i] - mexA0[i-1] + time*h[i-1]*sum; 01295 index++; 01296 sum = 0; 01297 for(j=0;j<3;j++) 01298 sum += omega[j][2]*mexBdot[i-1][j]; 01299 c[index] = mexB0[i] - mexB0[i-1] + time*h[i-1]*sum; 01300 index++; 01301 sum = 0; 01302 for(j=0;j<3;j++) 01303 sum += omega[j][2]*mraAdot[i-1][j]; 01304 c[index] = mraA0[i] - mraA0[i-1] + time*h[i-1]*sum; 01305 index++; 01306 sum = 0; 01307 for(j=0;j<3;j++) 01308 sum += omega[j][2]*mraBdot[i-1][j]; 01309 c[index] = mraB0[i] - mraB0[i-1] + time*h[i-1]*sum; 01310 index++; 01311 sum = 0; 01312 for(j=0;j<3;j++) 01313 sum += omega[j][2]*mfedot[i-1][j]; 01314 c[index] = mfe0[i] - mfe0[i-1] + time*h[i-1]*sum; 01315 index++; 01316 } 01317 01318 for(l=1;l<nex-1;l++) 01319 for(i=0;i<nel;i++) 01320 { 01321 for(j=0;j<3;j++) 01322 { 01323 c[index] = cAdot[l][i][j] - kA*q1*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 01324 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 01325 index++; 01326 c[index] = cBdot[l][i][j] - kB*q1*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 01327 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 01328 index++; 01329 } 01330 } 01331 01332 for(l=nex-1;l<nfe;l++) 01333 for(i=0;i<nel;i++) 01334 { 01335 for(j=0;j<3;j++) 01336 { 01337 c[index] = cAdot[l][i][j] - kA*q2*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 01338 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 01339 index++; 01340 c[index] = cBdot[l][i][j] - kB*q2*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 01341 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 01342 index++; 01343 } 01344 } 01345 01346 for(l=nfe+1;l<nra;l++) 01347 for(i=0;i<nel;i++) 01348 { 01349 for(j=0;j<3;j++) 01350 { 01351 c[index] = cAdot[l][i][j] - kA*q3*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 01352 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 01353 index++; 01354 c[index] = cBdot[l][i][j] - kB*q3*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 01355 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 01356 index++; 01357 } 01358 } 01359 01360 for(l=nra;l<nde;l++) 01361 { 01362 for(i=0;i<nel;i++) 01363 { 01364 for(j=0;j<3;j++) 01365 { 01366 c[index] = cAdot[l][i][j] - kA*q4*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ 01367 (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); 01368 index++; 01369 c[index] = cBdot[l][i][j] - kB*q4*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ 01370 (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); 01371 index++; 01372 } 01373 } 01374 c[index] = cA0[l][0] - cA[l-nra][nel-1][2]; 01375 index++; 01376 c[index] = cB0[l][0] - cB[l-nra][nel-1][2]; 01377 index++; 01378 } 01379 01380 for(l=0;l<nra;l++) 01381 { 01382 c[index] = cA0[l][0] - cA[l+ndis][nel-1][2]; 01383 index++; 01384 c[index] = cB0[l][0] - cB[l+ndis][nel-1][2]; 01385 index++; 01386 } 01387 01388 printf(" in ceval_ad index %d \n",index); 01389 /* for(i=0;i<cstr-5;i++) */ 01390 /* { */ 01391 /* printf("%d %e \n",i,c[i].value()); */ 01392 /* } */ 01393 } 01394 01395 01396 /* //\************** active version ************************ */ 01397 01398 01399 /* adouble feval_ad_mod(double *x, int n) */ 01400 /* { */ 01401 /* adouble cA[nde][nel][3],cB[nde][nel][3],cAdot[nde][nel][3],cBdot[nde][nel][3]; */ 01402 /* adouble cA0[nde][nel],cB0[nde][nel]; */ 01403 /* adouble mexA[nel][3],mexB[nel][3],mexAdot[nel][3],mexBdot[nel][3]; */ 01404 /* adouble mexA0[nel],mexB0[nel]; */ 01405 /* adouble mraA[nel][3],mraB[nel][3],mraAdot[nel][3],mraBdot[nel][3]; */ 01406 /* adouble mraA0[nel],mraB0[nel]; */ 01407 /* adouble mfe[nel][3],mfedot[nel][3]; */ 01408 /* adouble mfe0[nel]; */ 01409 01410 /* adouble q1,qde,qfe,qex,time; */ 01411 /* adouble q2,q3,q4,qra; */ 01412 01413 /* adouble res; */ 01414 /* adouble c[cstr]; */ 01415 /* double lam[cstr]; */ 01416 /* adouble sum; */ 01417 01418 /* int index; */ 01419 /* int i,j,l,k; */ 01420 01421 /* q1 = x[0]; qde = x[1]; qfe = x[2]; qex = x[3]; time = x[4]; */ 01422 01423 /* q2 = q1 - qex; */ 01424 /* q3 = q1 - qex + qfe; */ 01425 /* q4 = q1 - qde; */ 01426 /* qra = qde - qex + qfe; */ 01427 01428 /* index = 5; */ 01429 /* for(l=0;l<nde;l++) */ 01430 /* for(i=0;i<nel;i++) */ 01431 /* { */ 01432 /* for(j=0;j<3;j++) */ 01433 /* { */ 01434 /* cA[l][i][j] <<= x[index]; */ 01435 /* index++; */ 01436 /* cB[l][i][j] <<= x[index]; */ 01437 /* index++; */ 01438 /* cAdot[l][i][j] <<= x[index]; */ 01439 /* index++; */ 01440 /* cBdot[l][i][j] <<= x[index]; */ 01441 /* index++; */ 01442 /* } */ 01443 /* cA0[l][i] <<= x[index]; */ 01444 /* index++; */ 01445 /* cB0[l][i] <<= x[index]; */ 01446 /* index++; */ 01447 /* } */ 01448 01449 /* for(i=0;i<nel;i++) */ 01450 /* { */ 01451 /* for(j=0;j<3;j++) */ 01452 /* { */ 01453 /* mexA[i][j] <<= x[index]; */ 01454 /* index++; */ 01455 /* mexB[i][j] <<= x[index]; */ 01456 /* index++; */ 01457 /* mexAdot[i][j] <<= x[index]; */ 01458 /* index++; */ 01459 /* mexBdot[i][j] <<= x[index]; */ 01460 /* index++; */ 01461 /* mraA[i][j] <<= x[index]; */ 01462 /* index++; */ 01463 /* mraB[i][j] <<= x[index]; */ 01464 /* index++; */ 01465 /* mraAdot[i][j] <<= x[index]; */ 01466 /* index++; */ 01467 /* mraBdot[i][j] <<= x[index]; */ 01468 /* index++; */ 01469 /* mfe[i][j] <<= x[index]; */ 01470 /* index++; */ 01471 /* mfedot[i][j] <<= x[index]; */ 01472 /* index++; */ 01473 /* } */ 01474 /* mexA0[i] <<= x[index]; */ 01475 /* index++; */ 01476 /* mexB0[i] <<= x[index]; */ 01477 /* index++; */ 01478 /* mraA0[i] <<= x[index]; */ 01479 /* index++; */ 01480 /* mraB0[i] <<= x[index]; */ 01481 /* index++; */ 01482 /* mfe0[i] <<= x[index]; */ 01483 /* index++; */ 01484 /* } */ 01485 01486 /* res = -mfe[nel-1][2]/time; */ 01487 01488 /* index = 0; */ 01489 01490 /* for(l=0;l<nde;l++) */ 01491 /* for(i=0;i<nel;i++) */ 01492 /* { */ 01493 /* for(j=0;j<3;j++) */ 01494 /* { */ 01495 /* sum = 0; */ 01496 /* for(k=0;k<3;k++) */ 01497 /* sum += omega[k][j]*cAdot[l][i][k]; */ 01498 /* c[index] = cA[l][i][j] - cA0[l][i]+time*h[i]*sum; */ 01499 /* index++; */ 01500 /* sum = 0; */ 01501 /* for(k=0;k<3;k++) */ 01502 /* sum += omega[k][j]*cBdot[l][i][k]; */ 01503 /* c[index] = cB[l][i][j] - cB0[l][i]+time*h[i]*sum; */ 01504 /* index++; */ 01505 /* } */ 01506 /* } */ 01507 01508 /* for(l=0;l<nde;l++) */ 01509 /* for(i=1;i<nel;i++) */ 01510 /* { */ 01511 /* sum = 0; */ 01512 /* for(j=0;j<3;j++) */ 01513 /* sum += omega[j][2]*cAdot[l][i-1][j]; */ 01514 /* c[index] = cA0[l][i] - cA0[l][i-1] + time*h[i-1]*sum; */ 01515 /* index++; */ 01516 /* sum = 0; */ 01517 /* for(j=0;j<3;j++) */ 01518 /* sum += omega[j][2]*cBdot[l][i-1][j]; */ 01519 /* c[index] = cB0[l][i] - cB0[l][i-1] + time*h[i-1]*sum; */ 01520 /* index++; */ 01521 /* } */ 01522 01523 /* for(i=0;i<nel;i++) */ 01524 /* { */ 01525 /* for(j=0;j<3;j++) */ 01526 /* { */ 01527 /* sum = 0; */ 01528 /* for(k=0;k<3;k++) */ 01529 /* sum += omega[k][j]*mexAdot[i][k]; */ 01530 /* c[index] = mexA[i][j] - mexA0[i]+time*h[i]*sum; */ 01531 /* index++; */ 01532 /* sum = 0; */ 01533 /* for(k=0;k<3;k++) */ 01534 /* sum += omega[k][j]*mexBdot[i][k]; */ 01535 /* c[index] = mexB[i][j] - mexB0[i]+time*h[i]*sum; */ 01536 /* index++; */ 01537 /* sum = 0; */ 01538 /* for(k=0;k<3;k++) */ 01539 /* sum += omega[k][j]*mraAdot[i][k]; */ 01540 /* c[index] = mraA[i][j] - mraA0[i]+time*h[i]*sum; */ 01541 /* index++; */ 01542 /* sum = 0; */ 01543 /* for(k=0;k<3;k++) */ 01544 /* sum += omega[k][j]*mraBdot[i][k]; */ 01545 /* c[index] = mraB[i][j] - mraB0[i]+time*h[i]*sum; */ 01546 /* index++; */ 01547 /* sum = 0; */ 01548 /* for(k=0;k<3;k++) */ 01549 /* sum += omega[k][j]*mfedot[i][k]; */ 01550 /* c[index] = mfe[i][j]-mfe0[i]+time*h[i]*sum; */ 01551 /* index++; */ 01552 01553 /* c[index] = cAdot[0][i][j] - kA*(q4*cA[nde-1][i][j]-q1*cA[0][i][j])* */ 01554 /* (2*(1 + kk*cA[0][ i][ j])*(1 + kk*cA[0][ i][ j]))/ (1+(1 +kk*cA[0][i][j])*(1 +kk*cA[0][i][j])); */ 01555 /* index++; */ 01556 /* c[index] = cBdot[0][i][j] - kB*(q4*cB[nde-1][i][j]-q1*cB[0][i][j])* */ 01557 /* (2*(1 + kk*cB[0][ i][ j])*(1 + kk*cB[0][ i][ j]))/ (1+(1 +kk*cB[0][i][j])*(1 +kk*cB[0][i][j])); */ 01558 /* index++; */ 01559 01560 /* c[index] = cAdot[nfe][i][j] - kA*(q2*cA[nfe-1][i][j] + qfe*cFEA - */ 01561 /* q3*cA[nfe][i][j])*(2*(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j]))/ */ 01562 /* (1+(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j])); */ 01563 /* index++; */ 01564 /* c[index] = cBdot[nfe][i][j] - kB*(q2*cB[nfe-1][i][j] + qfe*cFEB - */ 01565 /* q3*cB[nfe][i][j])*(2*(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j]))/ */ 01566 /* (1+(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j])); */ 01567 /* index++; */ 01568 01569 /* c[index] = mexAdot[i][j] - qex*cA[nex-1][ i][j]; */ 01570 /* index++; */ 01571 /* c[index] = mexBdot[i][j] - qex*cB[nex-1][ i][j]; */ 01572 /* index++; */ 01573 /* c[index] = mraAdot[i][j] - qra*cA[nra-1][ i][j]; */ 01574 /* index++; */ 01575 /* c[index] = mraBdot[i][j] - qra*cB[nra-1][ i][j]; */ 01576 /* index++; */ 01577 /* c[index] = mfedot[i][j] - qfe; */ 01578 /* index++; */ 01579 /* } */ 01580 /* } */ 01581 01582 /* for(i=1;i<nel;i++) */ 01583 /* { */ 01584 /* sum = 0; */ 01585 /* for(j=0;j<3;j++) */ 01586 /* sum += omega[j][2]*mexAdot[i-1][j]; */ 01587 /* c[index] = mexA0[i] - mexA0[i-1] + time*h[i-1]*sum; */ 01588 /* index++; */ 01589 /* sum = 0; */ 01590 /* for(j=0;j<3;j++) */ 01591 /* sum += omega[j][2]*mexBdot[i-1][j]; */ 01592 /* c[index] = mexB0[i] - mexB0[i-1] + time*h[i-1]*sum; */ 01593 /* index++; */ 01594 /* sum = 0; */ 01595 /* for(j=0;j<3;j++) */ 01596 /* sum += omega[j][2]*mraAdot[i-1][j]; */ 01597 /* c[index] = mraA0[i] - mraA0[i-1] + time*h[i-1]*sum; */ 01598 /* index++; */ 01599 /* sum = 0; */ 01600 /* for(j=0;j<3;j++) */ 01601 /* sum += omega[j][2]*mraBdot[i-1][j]; */ 01602 /* c[index] = mraB0[i] - mraB0[i-1] + time*h[i-1]*sum; */ 01603 /* index++; */ 01604 /* sum = 0; */ 01605 /* for(j=0;j<3;j++) */ 01606 /* sum += omega[j][2]*mfedot[i-1][j]; */ 01607 /* c[index] = mfe0[i] - mfe0[i-1] + time*h[i-1]*sum; */ 01608 /* index++; */ 01609 /* } */ 01610 01611 /* for(l=1;l<nex-1;l++) */ 01612 /* for(i=0;i<nel;i++) */ 01613 /* { */ 01614 /* for(j=0;j<3;j++) */ 01615 /* { */ 01616 /* c[index] = cAdot[l][i][j] - kA*q1*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ */ 01617 /* (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); */ 01618 /* index++; */ 01619 /* c[index] = cBdot[l][i][j] - kB*q1*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ */ 01620 /* (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); */ 01621 /* index++; */ 01622 /* } */ 01623 /* } */ 01624 01625 /* for(l=nex-1;l<nfe;l++) */ 01626 /* for(i=0;i<nel;i++) */ 01627 /* { */ 01628 /* for(j=0;j<3;j++) */ 01629 /* { */ 01630 /* c[index] = cAdot[l][i][j] - kA*q2*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ */ 01631 /* (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); */ 01632 /* index++; */ 01633 /* c[index] = cBdot[l][i][j] - kB*q2*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ */ 01634 /* (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); */ 01635 /* index++; */ 01636 /* } */ 01637 /* } */ 01638 01639 /* for(l=nfe+1;l<nra;l++) */ 01640 /* for(i=0;i<nel;i++) */ 01641 /* { */ 01642 /* for(j=0;j<3;j++) */ 01643 /* { */ 01644 /* c[index] = cAdot[l][i][j] - kA*q3*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ */ 01645 /* (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); */ 01646 /* index++; */ 01647 /* c[index] = cBdot[l][i][j] - kB*q3*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ */ 01648 /* (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); */ 01649 /* index++; */ 01650 /* } */ 01651 /* } */ 01652 01653 /* for(l=nra;l<nde;l++) */ 01654 /* { */ 01655 /* for(i=0;i<nel;i++) */ 01656 /* { */ 01657 /* for(j=0;j<3;j++) */ 01658 /* { */ 01659 /* c[index] = cAdot[l][i][j] - kA*q4*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ */ 01660 /* (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); */ 01661 /* index++; */ 01662 /* c[index] = cBdot[l][i][j] - kB*q4*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ */ 01663 /* (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); */ 01664 /* index++; */ 01665 /* } */ 01666 /* } */ 01667 /* c[index] = cA0[l][0] - cA[l-nra][nel-1][2]; */ 01668 /* index++; */ 01669 /* c[index] = cB0[l][0] - cB[l-nra][nel-1][2]; */ 01670 /* index++; */ 01671 /* } */ 01672 01673 /* for(l=0;l<nra;l++) */ 01674 /* { */ 01675 /* c[index] = cA0[l][0] - cA[l+ndis][nel-1][2]; */ 01676 /* index++; */ 01677 /* c[index] = cB0[l][0] - cB[l+ndis][nel-1][2]; */ 01678 /* index++; */ 01679 /* } */ 01680 01681 /* c[index] = mexA0[0]; */ 01682 /* index++; */ 01683 /* c[index] = mexB0[0]; */ 01684 /* index++; */ 01685 /* c[index] = mraA0[0]; */ 01686 /* index++; */ 01687 /* c[index] = mraB0[0]; */ 01688 /* index++; */ 01689 /* c[index] = mfe0[0]; */ 01690 /* index++; */ 01691 01692 /* printf(" index = %d\n",index); */ 01693 01694 /* for(i=0;i<cstr;i++) */ 01695 /* lam[i] = 1; */ 01696 01697 /* for(i=0;i<cstr;i++) */ 01698 /* { */ 01699 /* res += lam[i]*c[i]; */ 01700 /* } */ 01701 01702 /* return res; */ 01703 01704 /* } */ 01705 01706 01707 /* //\************** active version ************************ */ 01708 01709 01710 /* adouble feval_ad_modHP(double *x, int n) */ 01711 /* { */ 01712 /* adouble cA[nde][nel][3],cB[nde][nel][3],cAdot[nde][nel][3],cBdot[nde][nel][3]; */ 01713 /* adouble cA0[nde][nel],cB0[nde][nel]; */ 01714 /* adouble mexA[nel][3],mexB[nel][3],mexAdot[nel][3],mexBdot[nel][3]; */ 01715 /* adouble mexA0[nel],mexB0[nel]; */ 01716 /* adouble mraA[nel][3],mraB[nel][3],mraAdot[nel][3],mraBdot[nel][3]; */ 01717 /* adouble mraA0[nel],mraB0[nel]; */ 01718 /* adouble mfe[nel][3],mfedot[nel][3]; */ 01719 /* adouble mfe0[nel]; */ 01720 01721 /* adouble q1,qde,qfe,qex,time; */ 01722 /* adouble q2,q3,q4,qra; */ 01723 01724 /* adouble res; */ 01725 /* adouble c[cstr]; */ 01726 /* double lam[cstr]; */ 01727 /* adouble sum; */ 01728 01729 /* int index; */ 01730 /* int i,j,l,k; */ 01731 01732 /* q1 = x[0]; qde = x[1]; qfe = x[2]; qex = x[3]; time = x[4]; */ 01733 01734 /* q2 = q1 - qex; */ 01735 /* q3 = q1 - qex + qfe; */ 01736 /* q4 = q1 - qde; */ 01737 /* qra = qde - qex + qfe; */ 01738 01739 /* index = 5; */ 01740 /* for(l=0;l<nde;l++) */ 01741 /* for(i=0;i<nel;i++) */ 01742 /* { */ 01743 /* for(j=0;j<3;j++) */ 01744 /* { */ 01745 /* cA[l][i][j] <<= x[index]; */ 01746 /* index++; */ 01747 /* cB[l][i][j] <<= x[index]; */ 01748 /* index++; */ 01749 /* cAdot[l][i][j] <<= x[index]; */ 01750 /* index++; */ 01751 /* cBdot[l][i][j] <<= x[index]; */ 01752 /* index++; */ 01753 /* } */ 01754 /* cA0[l][i] <<= x[index]; */ 01755 /* index++; */ 01756 /* cB0[l][i] <<= x[index]; */ 01757 /* index++; */ 01758 /* } */ 01759 01760 /* for(i=0;i<nel;i++) */ 01761 /* { */ 01762 /* for(j=0;j<3;j++) */ 01763 /* { */ 01764 /* mexA[i][j] <<= x[index]; */ 01765 /* index++; */ 01766 /* mexB[i][j] <<= x[index]; */ 01767 /* index++; */ 01768 /* mexAdot[i][j] <<= x[index]; */ 01769 /* index++; */ 01770 /* mexBdot[i][j] <<= x[index]; */ 01771 /* index++; */ 01772 /* mraA[i][j] <<= x[index]; */ 01773 /* index++; */ 01774 /* mraB[i][j] <<= x[index]; */ 01775 /* index++; */ 01776 /* mraAdot[i][j] <<= x[index]; */ 01777 /* index++; */ 01778 /* mraBdot[i][j] <<= x[index]; */ 01779 /* index++; */ 01780 /* mfe[i][j] <<= x[index]; */ 01781 /* index++; */ 01782 /* mfedot[i][j] <<= x[index]; */ 01783 /* index++; */ 01784 /* } */ 01785 /* mexA0[i] <<= x[index]; */ 01786 /* index++; */ 01787 /* mexB0[i] <<= x[index]; */ 01788 /* index++; */ 01789 /* mraA0[i] <<= x[index]; */ 01790 /* index++; */ 01791 /* mraB0[i] <<= x[index]; */ 01792 /* index++; */ 01793 /* mfe0[i] <<= x[index]; */ 01794 /* index++; */ 01795 /* } */ 01796 01797 /* res = -mfe[nel-1][2]/time; */ 01798 01799 /* index = 0; */ 01800 01801 /* for(l=0;l<nde;l++) */ 01802 /* for(i=0;i<nel;i++) */ 01803 /* { */ 01804 /* for(j=0;j<3;j++) */ 01805 /* { */ 01806 /* sum = 0; */ 01807 /* for(k=0;k<3;k++) */ 01808 /* sum += omega[k][j]*cAdot[l][i][k]; */ 01809 /* c[index] = cA[l][i][j] - cA0[l][i]+time*h[i]*sum; */ 01810 /* index++; */ 01811 /* sum = 0; */ 01812 /* for(k=0;k<3;k++) */ 01813 /* sum += omega[k][j]*cBdot[l][i][k]; */ 01814 /* c[index] = cB[l][i][j] - cB0[l][i]+time*h[i]*sum; */ 01815 /* index++; */ 01816 /* } */ 01817 /* } */ 01818 01819 /* for(l=0;l<nde;l++) */ 01820 /* for(i=1;i<nel;i++) */ 01821 /* { */ 01822 /* sum = 0; */ 01823 /* for(j=0;j<3;j++) */ 01824 /* sum += omega[j][2]*cAdot[l][i-1][j]; */ 01825 /* c[index] = cA0[l][i] - cA0[l][i-1] + time*h[i-1]*sum; */ 01826 /* index++; */ 01827 /* sum = 0; */ 01828 /* for(j=0;j<3;j++) */ 01829 /* sum += omega[j][2]*cBdot[l][i-1][j]; */ 01830 /* c[index] = cB0[l][i] - cB0[l][i-1] + time*h[i-1]*sum; */ 01831 /* index++; */ 01832 /* } */ 01833 01834 /* for(i=0;i<nel;i++) */ 01835 /* { */ 01836 /* for(j=0;j<3;j++) */ 01837 /* { */ 01838 /* sum = 0; */ 01839 /* for(k=0;k<3;k++) */ 01840 /* sum += omega[k][j]*mexAdot[i][k]; */ 01841 /* c[index] = mexA[i][j] - mexA0[i]+time*h[i]*sum; */ 01842 /* index++; */ 01843 /* sum = 0; */ 01844 /* for(k=0;k<3;k++) */ 01845 /* sum += omega[k][j]*mexBdot[i][k]; */ 01846 /* c[index] = mexB[i][j] - mexB0[i]+time*h[i]*sum; */ 01847 /* index++; */ 01848 /* sum = 0; */ 01849 /* for(k=0;k<3;k++) */ 01850 /* sum += omega[k][j]*mraAdot[i][k]; */ 01851 /* c[index] = mraA[i][j] - mraA0[i]+time*h[i]*sum; */ 01852 /* index++; */ 01853 /* sum = 0; */ 01854 /* for(k=0;k<3;k++) */ 01855 /* sum += omega[k][j]*mraBdot[i][k]; */ 01856 /* c[index] = mraB[i][j] - mraB0[i]+time*h[i]*sum; */ 01857 /* index++; */ 01858 /* sum = 0; */ 01859 /* for(k=0;k<3;k++) */ 01860 /* sum += omega[k][j]*mfedot[i][k]; */ 01861 /* c[index] = mfe[i][j]-mfe0[i]+time*h[i]*sum; */ 01862 /* index++; */ 01863 01864 /* c[index] = cAdot[0][i][j] - kA*(q4*cA[nde-1][i][j]-q1*cA[0][i][j])* */ 01865 /* (2*(1 + kk*cA[0][ i][ j])*(1 + kk*cA[0][ i][ j]))/ (1+(1 +kk*cA[0][i][j])*(1 +kk*cA[0][i][j])); */ 01866 /* index++; */ 01867 /* c[index] = cBdot[0][i][j] - kB*(q4*cB[nde-1][i][j]-q1*cB[0][i][j])* */ 01868 /* (2*(1 + kk*cB[0][ i][ j])*(1 + kk*cB[0][ i][ j]))/ (1+(1 +kk*cB[0][i][j])*(1 +kk*cB[0][i][j])); */ 01869 /* index++; */ 01870 01871 /* c[index] = cAdot[nfe][i][j] - kA*(q2*cA[nfe-1][i][j] + qfe*cFEA - */ 01872 /* q3*cA[nfe][i][j])*(2*(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j]))/ */ 01873 /* (1+(1 + kk*cA[nfe][i][j])*(1 + kk*cA[nfe][i][j])); */ 01874 /* index++; */ 01875 /* c[index] = cBdot[nfe][i][j] - kB*(q2*cB[nfe-1][i][j] + qfe*cFEB - */ 01876 /* q3*cB[nfe][i][j])*(2*(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j]))/ */ 01877 /* (1+(1 + kk*cB[nfe][i][j])*(1 + kk*cB[nfe][i][j])); */ 01878 /* index++; */ 01879 01880 /* c[index] = mexAdot[i][j] - qex*cA[nex-1][ i][j]; */ 01881 /* index++; */ 01882 /* c[index] = mexBdot[i][j] - qex*cB[nex-1][ i][j]; */ 01883 /* index++; */ 01884 /* c[index] = mraAdot[i][j] - qra*cA[nra-1][ i][j]; */ 01885 /* index++; */ 01886 /* c[index] = mraBdot[i][j] - qra*cB[nra-1][ i][j]; */ 01887 /* index++; */ 01888 /* c[index] = mfedot[i][j] - qfe; */ 01889 /* index++; */ 01890 /* } */ 01891 /* } */ 01892 01893 /* for(i=1;i<nel;i++) */ 01894 /* { */ 01895 /* sum = 0; */ 01896 /* for(j=0;j<3;j++) */ 01897 /* sum += omega[j][2]*mexAdot[i-1][j]; */ 01898 /* c[index] = mexA0[i] - mexA0[i-1] + time*h[i-1]*sum; */ 01899 /* index++; */ 01900 /* sum = 0; */ 01901 /* for(j=0;j<3;j++) */ 01902 /* sum += omega[j][2]*mexBdot[i-1][j]; */ 01903 /* c[index] = mexB0[i] - mexB0[i-1] + time*h[i-1]*sum; */ 01904 /* index++; */ 01905 /* sum = 0; */ 01906 /* for(j=0;j<3;j++) */ 01907 /* sum += omega[j][2]*mraAdot[i-1][j]; */ 01908 /* c[index] = mraA0[i] - mraA0[i-1] + time*h[i-1]*sum; */ 01909 /* index++; */ 01910 /* sum = 0; */ 01911 /* for(j=0;j<3;j++) */ 01912 /* sum += omega[j][2]*mraBdot[i-1][j]; */ 01913 /* c[index] = mraB0[i] - mraB0[i-1] + time*h[i-1]*sum; */ 01914 /* index++; */ 01915 /* sum = 0; */ 01916 /* for(j=0;j<3;j++) */ 01917 /* sum += omega[j][2]*mfedot[i-1][j]; */ 01918 /* c[index] = mfe0[i] - mfe0[i-1] + time*h[i-1]*sum; */ 01919 /* index++; */ 01920 /* } */ 01921 01922 /* for(l=1;l<nex-1;l++) */ 01923 /* for(i=0;i<nel;i++) */ 01924 /* { */ 01925 /* for(j=0;j<3;j++) */ 01926 /* { */ 01927 /* c[index] = cAdot[l][i][j] - kA*q1*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ */ 01928 /* (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); */ 01929 /* index++; */ 01930 /* c[index] = cBdot[l][i][j] - kB*q1*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ */ 01931 /* (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); */ 01932 /* index++; */ 01933 /* } */ 01934 /* } */ 01935 01936 /* for(l=nex-1;l<nfe;l++) */ 01937 /* for(i=0;i<nel;i++) */ 01938 /* { */ 01939 /* for(j=0;j<3;j++) */ 01940 /* { */ 01941 /* c[index] = cAdot[l][i][j] - kA*q2*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ */ 01942 /* (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); */ 01943 /* index++; */ 01944 /* c[index] = cBdot[l][i][j] - kB*q2*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ */ 01945 /* (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); */ 01946 /* index++; */ 01947 /* } */ 01948 /* } */ 01949 01950 /* for(l=nfe+1;l<nra;l++) */ 01951 /* for(i=0;i<nel;i++) */ 01952 /* { */ 01953 /* for(j=0;j<3;j++) */ 01954 /* { */ 01955 /* c[index] = cAdot[l][i][j] - kA*q3*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ */ 01956 /* (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); */ 01957 /* index++; */ 01958 /* c[index] = cBdot[l][i][j] - kB*q3*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ */ 01959 /* (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); */ 01960 /* index++; */ 01961 /* } */ 01962 /* } */ 01963 01964 /* for(l=nra;l<nde;l++) */ 01965 /* { */ 01966 /* for(i=0;i<nel;i++) */ 01967 /* { */ 01968 /* for(j=0;j<3;j++) */ 01969 /* { */ 01970 /* c[index] = cAdot[l][i][j] - kA*q4*(cA[l-1][i][j] - cA[l][i][j])*(2*(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j]))/ */ 01971 /* (1+(1 + kk*cA[l][i][j])*(1 + kk*cA[l][i][j])); */ 01972 /* index++; */ 01973 /* c[index] = cBdot[l][i][j] - kB*q4*(cB[l-1][i][j] - cB[l][i][j])*(2*(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j]))/ */ 01974 /* (1+(1 + kk*cB[l][i][j])*(1 + kk*cB[l][i][j])); */ 01975 /* index++; */ 01976 /* } */ 01977 /* } */ 01978 /* c[index] = cA0[l][0] - cA[l-nra][nel-1][2]; */ 01979 /* index++; */ 01980 /* c[index] = cB0[l][0] - cB[l-nra][nel-1][2]; */ 01981 /* index++; */ 01982 /* } */ 01983 01984 /* for(l=0;l<nra;l++) */ 01985 /* { */ 01986 /* c[index] = cA0[l][0] - cA[l+ndis][nel-1][2]; */ 01987 /* index++; */ 01988 /* c[index] = cB0[l][0] - cB[l+ndis][nel-1][2]; */ 01989 /* index++; */ 01990 /* } */ 01991 01992 /* c[index] = mexA0[0]; */ 01993 /* index++; */ 01994 /* c[index] = mexB0[0]; */ 01995 /* index++; */ 01996 /* c[index] = mraA0[0]; */ 01997 /* index++; */ 01998 /* c[index] = mraB0[0]; */ 01999 /* index++; */ 02000 /* c[index] = mfe0[0]; */ 02001 /* index++; */ 02002 02003 /* printf(" index = %d\n",index); */ 02004 02005 /* for(i=0;i<cstr;i++) */ 02006 /* lam[i] = 1; */ 02007 02008 /* for(i=0;i<cstr;i++) */ 02009 /* { */ 02010 /* res += lam[i]*c[i]; */ 02011 /* } */ 02012 02013 /* for(l=0;l<nde;l++) */ 02014 /* for(i=0;i<nel;i++) */ 02015 /* { */ 02016 /* for(j=0;j<3;j++) */ 02017 /* { */ 02018 /* res += cA[l][i][j]*cA[l][i][j]; */ 02019 /* res += cB[l][i][j]*cB[l][i][j]; */ 02020 /* res += cAdot[l][i][j]*cAdot[l][i][j]; */ 02021 /* res += cBdot[l][i][j]*cBdot[l][i][j]; */ 02022 /* } */ 02023 /* res += cA0[l][i]*cA0[l][i]; */ 02024 /* res += cA0[l][i]*cB0[l][i]; */ 02025 /* } */ 02026 02027 /* for(i=0;i<nel;i++) */ 02028 /* { */ 02029 /* for(j=0;j<3;j++) */ 02030 /* { */ 02031 /* res += mexA[i][j]*mexA[i][j]; */ 02032 /* res += mexB[i][j]*mexB[i][j]; */ 02033 /* res += mexAdot[i][j]*mexAdot[i][j]; */ 02034 /* res += mexBdot[i][j]*mexBdot[i][j]; */ 02035 /* res += mraA[i][j]*mraA[i][j]; */ 02036 /* res += mraB[i][j]*mraB[i][j]; */ 02037 /* res += mraAdot[i][j]*mraAdot[i][j]; */ 02038 /* res += mraBdot[i][j]*mraBdot[i][j]; */ 02039 /* res += mfe[i][j]*mfe[i][j]; */ 02040 /* res += mfedot[i][j]*mfedot[i][j]; */ 02041 /* } */ 02042 /* res += mexA0[i]*mexA0[i]; */ 02043 /* res += mexB0[i]*mexB0[i]; */ 02044 /* res += mraA0[i]*mraA0[i]; */ 02045 /* res += mraB0[i]*mraB0[i]; */ 02046 /* res += mfe0[i]*mfe0[i]; */ 02047 /* } */ 02048 02049 /* return res; */ 02050 02051 /* } */
1.7.6.1 
