00001 #include "adolc.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #define nel 2
00012 #define ndis 2 // =>
00013 #define cstr 406
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
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
00130
00131 void init_dim(int *n, int *m)
00132 {
00133
00134
00135 *n = 5 + 4*(nde*nel*3) + 2*(nde*nel) + 8*(nel*3) + 4*nel
00136 + 2*nel*3 + nel;
00137
00138 *m = cstr-5;
00139
00140 }
00141
00142
00143
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
00196
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
00287
00288 res = -mfe[nel-1][2]/time;
00289
00290
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
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
00821
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
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
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
01390
01391
01392
01393 }
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045
02046
02047
02048
02049
02050
02051