00001 #include "arclength.h"
00002 #include "nssolver.h"
00003 #include "mechprint.h"
00004 #include "mathem.h"
00005 #include "global.h"
00006 #include "globmat.h"
00007
00008 #include "mtsolver.h"
00009
00010 #include <stdlib.h>
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 double garclength(long lcid, nonlinman *nlman, double *ra, double *fa, double *fc, double *fp,
00037 long li, double ilambda, flagsw outres)
00038 {
00039 long i;
00040 long n;
00041 long ni;
00042 long modif;
00043 double dl;
00044 double dlmin;
00045 double dlmax;
00046 double psi;
00047 double psi0;
00048 double lambda;
00049 double norfp;
00050 double tmp;
00051 double *r;
00052 double *ddr;
00053 double *u;
00054 double *v;
00055 double *f;
00056 double *fi;
00057 double totl=0.0;
00058 long modif_psi;
00059 aux_arcl_vec aav;
00060
00061
00062
00063 n = Ndofm;
00064
00065 lambda = ilambda;
00066
00067 ni = nlman->nial;
00068
00069 ini = nlman->niilal;
00070
00071 ierr = nlman->erral;
00072
00073 dl = nlman->dlal;
00074
00075 dlmax = nlman->dlmaxal;
00076
00077 dlmin = nlman->dlminal;
00078
00079 psi = nlman->psial;
00080
00081
00082 r = new double[n];
00083 ddr = new double[n];
00084 u = new double[n];
00085 v = new double[n];
00086 f = new double[n];
00087 fi = new double[n];
00088
00089 aav.r = r;
00090 aav.ddr = ddr;
00091 aav.u = u;
00092 aav.v = v;
00093 aav.f = f;
00094 aav.fi = fi;
00095
00096
00097
00098 dlambda=0.0;
00099
00100
00101
00102
00103 if (lambda == 0.0)
00104 assemble_attained_load_vector(fa, fc, fp, n, lambda, nlman);
00105
00106 norfp = normv(fp, n);
00107
00108 modif=0;
00109 psi0 = psi;
00110 if (psi<0.0){
00111
00112 modif_psi=1;
00113 }else{
00114
00115 modif_psi=0;
00116 }
00117
00118
00119
00120
00121
00122 for (i=li;i<ni;i++)
00123 {
00124
00125
00126 ret = one_step_arcl(lcid, nlman, ra, fa, fc, fp, aav, norfp,
00127 i, li, lambda, totl, psi, psi0, modif_psi, modif, outres);
00128
00129
00130 if (ret == 0)
00131 {
00132
00133 if ((nlman->check_lambdar == on) && (fabs(lambda-nlman->lambdar) < nlman->errl))
00134 break;
00135
00136
00137 if ((nlman->check_tot_al == on) && (totl == nlman->max_tot_al))
00138 break;
00139
00140 if (modif>1)
00141 {
00142
00143 dl*=2.0;
00144 if (dl>dlmax)
00145 dl=dlmax;
00146 modif=0;
00147 if (Mespr==1)
00148 fprintf(stdout,"\n arc-length is be modified (dl increase) dl=%e",dl);
00149 }
00150 }
00151
00152
00153
00154 if (ret == 1)
00155 {
00156 if (dl == dlmin)
00157 {
00158 if (Mespr==1)
00159 {
00160 print_err("the step length was decreased to the minimum required value\n"
00161 " but the equilibrium could not be attained.\n"
00162 " FORCED output of the attained results was performed in this step.",
00163 __FILE__, __LINE__, __func__);
00164 }
00165 compute_req_val (lcid);
00166 print_step_forced(lcid, i+1, lambda, f);
00167 print_flush();
00168 break;
00169 }
00170
00171 dl/=2.0;
00172 if (dl<dlmin)
00173 dl=dlmin;
00174 if (Mespr==1)
00175 {
00176 if (dl == dlmin)
00177 fprintf (stdout,"\n arc length was decreased dl=dlmin=%e",dl);
00178 else
00179 fprintf (stdout,"\n arc length was decreased dl=%e",dl);
00180 }
00181 }
00182
00183
00184
00185 }
00186
00187 if (Mp->hdbcont.save_stat())
00188 {
00189
00190 if (Mespr==1)
00191 fprintf (stdout,"\n Creating backup file\n");
00192 solver_save (ra, fp, i, lambda, dl, NULL, n);
00193 }
00194
00195
00196 delete [] r;
00197 delete [] fi;
00198 delete [] f;
00199 delete [] v;
00200 delete [] u;
00201 delete [] ddr;
00202
00203 return lambda;
00204 }
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 long one_step_arcl(long lcid, nonlinman *nlman, double *ra, double *fa, double *fc, double *fp,
00244 aux_arcl_vec &aav, double norfp, long i, long li, double &lambda, double &totl,
00245 double &psi, double &psi0, long modif_psi, long &modif, flagsw outres)
00246 {
00247 long j;
00248 long n;
00249 long ini;
00250 double dl;
00251 double dlmin;
00252 double dlmax;
00253 double blambda;
00254 double dlambda;
00255 double ddlambda;
00256 double ddlambda0;
00257 double norf;
00258 double norv;
00259 double norfa;
00260 double zero;
00261 double ierr;
00262 double tmp;
00263
00264 double *r;
00265 double *ddr;
00266 double *u;
00267 double *v;
00268 double *f;
00269 double *fi;
00270
00271 double totl = 0.0;
00272 long stopj;
00273 long sm_change = 0;
00274 flagsw check_rv_fp;
00275 flagsw check_max_tot_al;
00276
00277 matrix lsm_a(3,3);
00278 vector lsm_r(3);
00279 vector lsm_l(3);
00280
00281
00282
00283 n = Ndofm;
00284
00285 zero = Mp->zero;
00286
00287 ini = nlman->niilal;
00288
00289 ierr = nlman->erral;
00290
00291 dl = nlman->dlal;
00292
00293 dlmax = nlman->dlmaxal;
00294
00295 dlmin = nlman->dlminal;
00296
00297 psi = nlman->psial;
00298
00299 check_rv_fp = nlman->check_lambdar;
00300
00301 check_max_tot_al = nlman->check_tot_al;
00302
00303
00304 r = aav.r;
00305 ddr = aav.ddr;
00306 u = aav.u;
00307 v = aav.v;
00308 f = aav.f;
00309 fi = aav.fi;
00310
00311
00312
00313 dlambda=0.0;
00314
00315
00316
00317 Mp->strainstate=0;
00318
00319
00320 Mp->stressstate=0;
00321
00322
00323 Mp->otherstate=0;
00324
00325 Mp->istep = i;
00326 Mp->jstep = -1;
00327
00328
00329 copyv(ra, r, n);
00330
00331 blambda=lambda;
00332
00333 if (Mespr==1)
00334 {
00335 fprintf (stdout,"\n increment %ld, lambda=%e dl=%e psi=%e",i,lambda,dl,psi);
00336
00337 }
00338
00339
00340 sm_change = assemble_stiffness_matrix (lcid,i,-1,li);
00341
00342
00343 copyv(fp, f, n);
00344
00345
00346 Mp->ssle->solve_system(Gtm,Smat,v,f,Out);
00347
00348
00349 norv = displincr(v,v,n);
00350
00351 if ((i == 0) && (psi < 0.0) && (norfp > Mp->zero))
00352 psi0 = psi = sqrt(norv/norfp/norfp);
00353
00354
00355 tmp = norv+psi*psi*norfp*norfp;
00356
00357 dlambda = dl/sqrt(tmp);
00358
00359
00360 if ((check_rv_fp == on) && (dlambda+lambda > nlman->lambdar))
00361 {
00362 check_rv_fp = off;
00363 dlambda = nlman->lambdar - lambda;
00364 dl = dlambda*sqrt(tmp);
00365 }
00366
00367 if ((check_max_tot_al == on) && (dl+totl > nlman->max_tot_al))
00368 {
00369 check_max_tot_al = off;
00370 dl = nlman->max_tot_al - totl;
00371 }
00372
00373 if (Mespr == 1)
00374 fprintf(Out, "Step %ld: ddlambda = %e, norv = %e, norfp = %e\n", i, ddlambda, norv, norfp);
00375
00376
00377 copymultv(v, ddr, dlambda, n);
00378
00379
00380 addv(ra, ddr, n);
00381
00382
00383
00384 addmultv(fa, fp, dlambda, n);
00385
00386
00387
00388
00389 ddlambda=dlambda;
00390
00391
00392 internal_forces (lcid,fi);
00393
00394
00395 subv(fa, fi, f, n);
00396
00397
00398 norf=normv(f,n);
00399 norfa = normv(fa, n);
00400
00401
00402
00403 if (fabs(norfa) > Mp->zero)
00404 norf /= norfa;
00405
00406 if (norf<ierr)
00407 {
00408
00409
00410
00411 if (Mespr==1)
00412 fprintf(stdout,"\n increment %ld norf=%e ierr=%e ddlambda=%e dlambda=%e",i,norf,ierr, ddlambda, dlambda);
00413
00414
00415 modif++;
00416 lambda+=dlambda;
00417
00418 if ((nlman->check_tot_al == on) && (check_max_tot_al==off))
00419 totl = nlman->max_tot_al;
00420 else
00421 totl += dl;
00422
00423 if (modif_psi)
00424 adapt_psi_coeff(i, ddlambda, ddlambda0, psi0, psi);
00425
00426 Mm->updateipval();
00427 if (outres == yes)
00428 {
00429 compute_req_val(lcid);
00430 print_step(lcid, i+1, lambda, fa);
00431 print_flush();
00432 }
00433 return 0;
00434 }
00435
00436
00437
00438
00439 if (Mespr==1)
00440 fprintf (stdout,"\n inner loop: 1 norf=%e ierr=%e ddlambda=%e dlambda=%e",norf,ierr,ddlambda, dlambda);
00441 fillm(0.0, lsm_a);
00442 fillv(0.0, lsm_r);
00443 for (j=0;j<ini;j++)
00444 {
00445
00446
00447 Mp->jstep = j;
00448
00449
00450 Mespr = 0;
00451 sm_change = assemble_stiffness_matrix (lcid,i,j,li);
00452 Mespr = 1;
00453
00454
00455 Mp->ssle->solve_system(Gtm,Smat,u,f,Out);
00456 if (sm_change)
00457 {
00458
00459 copyv(fp, f, n);
00460 Mp->ssle->solve_system(Gtm,Smat,v,f,Out);
00461 }
00462
00463
00464
00465 copyv(ddr, f, n);
00466
00467
00468
00469 addv(ddr, u, n);
00470
00471
00472 determine_dlambda (ddr, u, v, f, n, ddlambda,psi,norfp,dl,dlambda,stopj, nlman);
00473 if (stopj)
00474 break;
00475
00476
00477 if ((nlman->check_lambdar == on) && (dlambda+ddlambda+lambda > nlman->lambdar))
00478 dlambda = nlman->lambdar - ddlambda - lambda;
00479
00480
00481
00482 addmultv(ddr, v, dlambda, n);
00483
00484
00485
00486 addv(ra, u, n);
00487 addmultv(ra, v, dlambda, n);
00488
00489
00490 addmultv(fa, fp, dlambda, n);
00491
00492
00493 ddlambda+=dlambda;
00494
00495
00496
00497 internal_forces (lcid,fi);
00498
00499 subv(fa, fi, f, n);
00500
00501 norf=normv(f,n);
00502
00503 norfa = normv(fa, n);
00504 if (fabs(norfa) > zero)
00505 {
00506
00507
00508 norf /= norfa;
00509 }
00510
00511 if (Mespr==1)
00512 {
00513 fprintf (stdout,"\n increment %ld inner loop %ld norf=%e",i,j,norf);
00514
00515 }
00516
00517
00518 if (norf<ierr)
00519 {
00520 lambda+=ddlambda;
00521
00522 if ((nlman->check_tot_al == on) && (check_max_tot_al==off))
00523 totl = nlman->max_tot_al;
00524 else
00525 totl += dl;
00526
00527 if (modif_psi)
00528 adapt_psi_coeff(i, ddlambda, ddlambda0, psi0, psi);
00529
00530 Mm->updateipval();
00531 compute_req_val (lcid);
00532 print_step(lcid, i+1, lambda, fa);
00533 print_flush();
00534
00535 modif = 0;
00536 return 0;
00537 }
00538
00539
00540 if (check_divergency(nlman, lsm_a, lsm_r, lsm_l, j, norf))
00541 break;
00542 }
00543
00544
00545
00546
00547
00548 modif=0;
00549
00550 copyv (r, ra, n);
00551
00552 lambda=blambda;
00553
00554 return 1;
00555 }
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568 void seldofinit ()
00569 {
00570 long i,j,k,l,ndofn;
00571 double x1,x2,y1,y2,z1,z2,length;
00572
00573 switch (Mp->nlman->displnorm){
00574 case alldofs:{ break; }
00575 case seldofs:
00576 case seldofscoord:{
00577 for (i=0;i<Mp->nlman->nsdofal;i++){
00578 j=Mp->nlman->seldofal[i];
00579 Mp->nlman->seldofal[i]=Mt->give_dof (Mp->nlman->selnodal[i],j)-1;
00580 }
00581 break;
00582 }
00583 case selecnodes:{
00584 Mp->nlman->nsdofal=0;
00585 for (i=0;i<Mp->nlman->nsnal;i++){
00586 Mp->nlman->nsdofal+=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00587 }
00588 Mp->nlman->seldofal = new long [Mp->nlman->nsdofal];
00589 k=0;
00590 for (i=0;i<Mp->nlman->nsnal;i++){
00591 ndofn=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00592 for (j=0;j<ndofn;j++){
00593 l=Mt->give_dof (Mp->nlman->selnodal[i],j);
00594 if (l>0){
00595 Mp->nlman->seldofal[k]=l-1;
00596 k++;
00597 }
00598 }
00599 }
00600 Mp->nlman->nsdofal=k;
00601 break;
00602 }
00603 case nodesdistincr:{
00604 Mp->nlman->nsdofal=0;
00605 for (i=0;i<Mp->nlman->nsnal;i++){
00606 Mp->nlman->nsdofal+=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00607 }
00608 Mp->nlman->seldofal = new long [Mp->nlman->nsdofal];
00609 k=0;
00610 for (i=0;i<Mp->nlman->nsnal;i++){
00611 ndofn=Mt->give_ndofn (Mp->nlman->selnodal[i]);
00612 for (j=0;j<ndofn;j++){
00613 Mp->nlman->seldofal[k]=Mt->give_dof (Mp->nlman->selnodal[i],j);
00614 k++;
00615 }
00616 }
00617 x1=Gtm->gnodes[Mp->nlman->selnodal[0]].x; x2=Gtm->gnodes[Mp->nlman->selnodal[1]].x;
00618 y1=Gtm->gnodes[Mp->nlman->selnodal[0]].y; y2=Gtm->gnodes[Mp->nlman->selnodal[1]].y;
00619 z1=Gtm->gnodes[Mp->nlman->selnodal[0]].z; z2=Gtm->gnodes[Mp->nlman->selnodal[1]].z;
00620 length=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
00621 Mp->nlman->nxal=(x2-x1)/length; Mp->nlman->nyal=(y2-y1)/length; Mp->nlman->nzal=(z2-z1)/length;
00622 break;
00623 }
00624 default:{
00625 fprintf (stderr,"\n\n unknown displacement norm is required in function seldofinit (%s, line %d)",__FILE__,__LINE__);
00626 }
00627 }
00628 }
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645 double loadincr (double *fp, long n, double lambda, nonlinman *nlman)
00646 {
00647 long i;
00648 double norm;
00649 double *loadincr = new double[n];
00650
00651 norm = 0.0;
00652
00653 nullv(loadincr, n);
00654 if ((nlman->check_lambdar == off) || (lambda > nlman->lambdar))
00655 addv(loadincr, fp, n);
00656
00657 switch (Mp->nlman->displnorm)
00658 {
00659 case alldofs:
00660 case nodesdistincr:
00661 norm = ss (loadincr,loadincr,n);
00662 break;
00663 case seldofs:
00664 case seldofscoord:
00665 case selecnodes:
00666 for (i=0;i<Mp->nlman->nsdofal;i++)
00667 norm+=loadincr[Mp->nlman->seldofal[i]]*loadincr[Mp->nlman->seldofal[i]];
00668 break;
00669 default:
00670 break;
00671 }
00672 norm = sqrt(norm);
00673 delete [] loadincr;
00674
00675 return norm;
00676 }
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694 double displincr (double *dr1, double *dr2, long n)
00695 {
00696 long i;
00697 double norm,u1,u2,v1,v2,w1,w2;
00698
00699 switch (Mp->nlman->displnorm){
00700 case alldofs:{
00701 norm = scprd(dr1,dr2,n);
00702 break;
00703 }
00704 case seldofs:
00705 case seldofscoord:
00706 case selecnodes:{
00707 norm=0.0;
00708 for (i=0;i<Mp->nlman->nsdofal;i++){
00709 norm+=dr1[Mp->nlman->seldofal[i]]*dr2[Mp->nlman->seldofal[i]];
00710 }
00711 break;
00712 }
00713 case nodesdistincr:{
00714 if (Mp->nlman->probdimal==2){
00715 u1=0.0; u2=0.0; v1=0.0; v2=0.0;
00716 if (Mp->nlman->seldofal[0]>0) u1=dr1[Mp->nlman->seldofal[0]-1];
00717 if (Mp->nlman->seldofal[1]>0) v1=dr1[Mp->nlman->seldofal[1]-1];
00718 if (Mp->nlman->seldofal[2]>0) u2=dr2[Mp->nlman->seldofal[2]-1];
00719 if (Mp->nlman->seldofal[3]>0) v2=dr2[Mp->nlman->seldofal[3]-1];
00720
00721 norm=(u2-u1)*Mp->nlman->nxal + (v2-v1)*Mp->nlman->nyal;
00722 }
00723 if (Mp->nlman->probdimal==3){
00724 u1=0.0; u2=0.0; v1=0.0; v2=0.0; w1=0.0; w2=0.0;
00725 if (Mp->nlman->seldofal[0]>0) u1=dr1[Mp->nlman->seldofal[0]-1];
00726 if (Mp->nlman->seldofal[1]>0) v1=dr1[Mp->nlman->seldofal[1]-1];
00727 if (Mp->nlman->seldofal[2]>0) w1=dr1[Mp->nlman->seldofal[2]-1];
00728 if (Mp->nlman->seldofal[3]>0) u2=dr2[Mp->nlman->seldofal[3]-1];
00729 if (Mp->nlman->seldofal[4]>0) v2=dr2[Mp->nlman->seldofal[4]-1];
00730 if (Mp->nlman->seldofal[5]>0) w2=dr2[Mp->nlman->seldofal[5]-1];
00731
00732 norm=(u2-u1)*Mp->nlman->nxal + (v2-v1)*Mp->nlman->nyal + (w2-w1)*Mp->nlman->nzal;
00733 }
00734 break;
00735 }
00736 default:
00737 print_err("unknown norm of displacement increment is required", __FILE__, __LINE__, __func__);
00738 }
00739 return norm;
00740 }
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763 void quadeqcoeff (double *ddr,double *v,long n,double ddlambda,double psi,double norfp,double dl,
00764 double &a0,double &a1,double &a2)
00765 {
00766 double norddr, norv, norddrv;
00767
00768 switch (Mp->nlman->displnorm){
00769 case alldofs:
00770 case seldofs:
00771 case seldofscoord:
00772 case selecnodes:{
00773 a0 = norddr = displincr(ddr, ddr, n);
00774 a2 = displincr(v, v, n);
00775 norddrv = displincr(ddr, v, n);
00776 a1 = 2.0*norddrv;
00777
00778 a0 += ddlambda*ddlambda*psi*psi*norfp*norfp - dl*dl;
00779 a1 += 2.0*ddlambda*psi*psi*norfp*norfp;
00780 a2 += psi*psi*norfp*norfp;
00781 break;
00782 }
00783 case nodesdistincr:{
00784 norddr = displincr(ddr, ddr, n);
00785 norv = displincr(v, v, n);
00786 norddrv = norddr*norv;
00787 a0 = norddr*norddr + ddlambda*ddlambda*norfp*norfp*psi*psi-dl*dl;
00788 a1 = 2.0*(norddrv + ddlambda*norfp*norfp*psi*psi);
00789 a2 = norv*norv + norfp*norfp*psi*psi;
00790 break;
00791 }
00792 default:
00793 print_err("unknown norm of displacement increment is required", __FILE__, __LINE__, __func__);
00794 }
00795
00796 }
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823 void determine_dlambda (double *ddr, double *u, double *v, double *ddrprev, long n,
00824 double ddlambda, double psi, double norfp, double dl,
00825 double &dlambda, long &stop, nonlinman *nlman)
00826 {
00827 long k,numr;
00828 double l1,l2,aux,a0,a1,a2,ss1,ss2,ss3,ss4,ss5,nom,denom;
00829
00830 switch (Mp->nlman->dlam)
00831 {
00832 case nodetermination:{
00833 print_err("some reasonable type of lambda determination has to be selected",__FILE__,__LINE__,__func__);
00834 break;
00835 }
00836 case minvalue:
00837 case maxvalue:
00838 case minangle:{
00839
00840 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2);
00841
00842 numr = solv_polynom_2(a2, a1, a0, l1, l2);
00843 switch (numr){
00844 case -1:{
00845 print_err("infinite number of solution of constrained condition\n"
00846 "(all coefficients of quadratic equation are equal to zero", __FILE__, __LINE__, __func__);
00847 abort();
00848 break;
00849 }
00850 case 0:{
00851 print_err("nonsolvable constrained condition in function arclength", __FILE__, __LINE__, __func__);
00852 stop=0;
00853 break;
00854 }
00855 case 1:{
00856 dlambda = l1;
00857 break;
00858 }
00859 default:{}
00860 }
00861 break;
00862 }
00863 case linearizedmeth:{
00864 break;
00865 }
00866 case fullmethod:
00867
00868
00869 a0 = displincr(ddrprev, ddrprev, n);
00870 a0 += ddlambda*ddlambda*psi*psi*norfp*norfp - dl*dl;
00871 break;
00872 default:{
00873 print_err("unknown type of detemination of lambda parameter is required",__FILE__,__LINE__,__func__);
00874 }
00875 }
00876
00877
00878 switch (Mp->nlman->dlam){
00879 case nodetermination:{
00880 print_err("some reasonable type of lambda determination has to be selected",__FILE__,__LINE__,__func__);
00881 break;
00882 }
00883 case minvalue:{
00884 if (fabs(l1)>fabs(l2))
00885 dlambda=l2;
00886 else
00887 dlambda=l1;
00888 break;
00889 }
00890 case maxvalue:{
00891 if (fabs(l1)>fabs(l2))
00892 dlambda=l1;
00893 else
00894 dlambda=l2;
00895 break;
00896 }
00897 case minangle:{
00898 ss1=0.0; ss2=0.0;
00899 ss3=0.0; ss4=0.0; ss5=0.0;
00900 for (k=0;k<n;k++){
00901 ss1+=(ddr[k]+l1*v[k])*ddrprev[k];
00902 ss2+=(ddr[k]+l2*v[k])*ddrprev[k];
00903 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]);
00904 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]);
00905 ss5+=ddrprev[k]*ddrprev[k];
00906 }
00907 if (ss1/sqrt(ss3)/sqrt(ss5)>ss2/sqrt(ss4)/sqrt(ss5))
00908 dlambda=l1;
00909 else
00910 dlambda=l2;
00911 break;
00912 }
00913 case linearizedmeth:{
00914 nom=0.0;
00915 denom=0.0;
00916 for (k=0;k<n;k++){
00917 nom-=ddrprev[k]*(ddr[k]-ddrprev[k]);
00918 denom+=ddrprev[k]*v[k];
00919 }
00920 denom+=psi*psi*norfp*norfp*ddlambda;
00921
00922 dlambda=nom/denom;
00923 break;
00924 }
00925 case fullmethod:
00926
00927
00928
00929 aux = displincr(ddrprev, u, n);
00930 dlambda = -a0 - aux;
00931 aux = displincr(ddrprev, v, n);
00932 dlambda /= 2.0*(aux + ddlambda*psi*psi*norfp*norfp);
00933 break;
00934 default:{
00935 print_err("unknown type of detemination of lambda parameter is required",__FILE__,__LINE__,__func__);
00936 }
00937 }
00938 }
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962 void adapt_psi_coeff(long i, double ddlambda, double &ddlambda0, double &psi0, double &psi)
00963 {
00964
00965
00966 if (i==0)
00967 ddlambda0 = ddlambda;
00968
00969
00970
00971 if ((sgn(ddlambda0) != sgn(ddlambda)) && (sgn(ddlambda) != 0.0))
00972 {
00973 psi0= psi;
00974 ddlambda0 = ddlambda;
00975 }
00976
00977
00978
00979 if (ddlambda0 != 0.0)
00980 psi = psi0*fabs(ddlambda/ddlambda0);
00981
00982 return;
00983 }