00001 #include <math.h>
00002 #include <string.h>
00003 #include <stdlib.h>
00004 #include "newtonraph.h"
00005 #include "nssolver.h"
00006 #include "backupsol.h"
00007 #include "global.h"
00008 #include "globmat.h"
00009 #include "mechprint.h"
00010 #include "gmatrix.h"
00011 #include "gtopology.h"
00012 #include "mechprint.h"
00013 #include "mathem.h"
00014 #include "vector.h"
00015 #include "matrix.h"
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 double gnewton_raphson (long lcid, nonlinman *nlman, double *ra, double *fa, double *fc, double *fp,
00042 long li, double ilambda, answertype outres)
00043 {
00044 long i,j,n,ni,ini,modif;
00045 double lambda;
00046 double blambda;
00047 double dlambda;
00048 double dlambdamin;
00049 double dlambdamax;
00050 double zero;
00051 double ierr;
00052 double norf;
00053 double norfa;
00054 double *rb;
00055 double *dr;
00056 double *fi;
00057 double *fb;
00058 matrix lsm_a(3,3);
00059 vector lsm_r(3), lsm_l(3);
00060 flagsw check_rv_fp = nlman->check_lambdar;
00061
00062
00063
00064 n = Ndofm;
00065
00066 zero = Mp->zero;
00067
00068 lambda = ilambda;
00069
00070 ni = nlman->ninr;
00071
00072 ini = nlman->niilnr;
00073
00074 ierr = nlman->errnr;
00075
00076 dlambda=nlman->incrnr;
00077
00078 dlambdamin=nlman->minincrnr;
00079
00080 dlambdamax=nlman->maxincrnr;
00081
00082
00083 rb = new double [n];
00084 fb = new double [n];
00085 fi = new double [n];
00086 dr = new double [n];
00087 memset (rb,0,n*sizeof(double));
00088 memset (fb,0,n*sizeof(double));
00089 memset (fi,0,n*sizeof(double));
00090 memset (dr,0,n*sizeof(double));
00091
00092 modif=0;
00093
00094
00095
00096
00097
00098
00099
00100 if (lambda == 0.0)
00101 assemble_attained_load_vector(fa, fc, fp, n, lambda, nlman);
00102
00103 for (i=li;i<ni;i++)
00104 {
00105
00106
00107 Mp->strainstate=0;
00108
00109
00110 Mp->stressstate=0;
00111
00112
00113 Mp->otherstate=0;
00114
00115 Mp->istep = i;
00116 Mp->jstep = -1;
00117
00118
00119
00120 copyv(ra, rb, n);
00121
00122
00123 blambda=lambda;
00124 Mp->lambda = lambda;
00125
00126 fprintf (stdout,"\n increment %ld lambda %e, dlambda %e",i,lambda,dlambda);
00127
00128
00129
00130 if ((check_rv_fp == on) && (dlambda+lambda > nlman->lambdar))
00131 {
00132 check_rv_fp = off;
00133 dlambda = nlman->lambdar - lambda;
00134 }
00135
00136
00137 addmultv(fa, fp, dlambda, n);
00138
00139
00140
00141 Mp->lambda += dlambda;
00142 Mp->dlambda = dlambda;
00143
00144
00145 internal_forces(lcid, fi);
00146
00147
00148
00149
00150
00151
00152
00153
00154 subv(fa, fi, fb, n);
00155
00156
00157
00158
00159 assemble_stiffness_matrix(lcid,i,-1,li);
00160
00161
00162 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00163
00164
00165
00166 addv(ra, dr, n);
00167
00168
00169 internal_forces (lcid,fi);
00170
00171
00172
00173 subv(fa, fi, fb, n);
00174
00175
00176 norfa=normv(fa,n);
00177
00178 norf=normv(fb,n);
00179 if (norfa != 0.0)
00180 norf /= norfa;
00181
00182 fprintf (stdout,"\n norm %e",norf);
00183
00184
00185 if (norf<ierr)
00186 {
00187 lambda+=dlambda;
00188 Mm->updateipval();
00189 if (outres == yes)
00190 {
00191 compute_req_val(lcid);
00192 print_step(lcid, i+1, lambda, fa);
00193 print_flush();
00194 }
00195 modif++;
00196 if (modif>1)
00197 {
00198 dlambda*=2.0;
00199 if (dlambda>dlambdamax)
00200 dlambda=dlambdamax;
00201 if (Mespr==1)
00202 fprintf (stdout,"\n increment must be modified (dlambda increase) dlambda=%e",dlambda);
00203 modif=0;
00204 }
00205
00206 if ((nlman->check_lambdar == on) && (check_rv_fp == off))
00207 break;
00208
00209 continue;
00210 }
00211
00212
00213 fillm(0.0, lsm_a);
00214 fillv(0.0, lsm_r);
00215 for (j=0;j<ini;j++)
00216 {
00217 Mp->jstep = j;
00218
00219
00220 assemble_stiffness_matrix(lcid,i,j,li);
00221
00222
00223 Mp->ssle->solve_system(Gtm,Smat,dr,fb,Out);
00224
00225
00226 addv(ra, dr, n);
00227
00228
00229 internal_forces(lcid,fi);
00230
00231
00232
00233 subv(fa, fi, fb, n);
00234
00235
00236 norf=normv(fb,n);
00237 if (norfa != 0.0)
00238 norf /= norfa;
00239
00240 fprintf (stdout,"\n increment %ld inner loop j %ld norf %e dlambda %e",i,j,norf,dlambda);
00241
00242
00243 if (norf<ierr)
00244
00245 break;
00246
00247
00248 if (check_divergency(nlman, lsm_a, lsm_r, lsm_l, j, norf))
00249 break;
00250 }
00251
00252 modif=0;
00253
00254 if (j==ini || norf>ierr)
00255 {
00256
00257
00258 addmultv(fa, fp, -dlambda, n);
00259
00260 copyv(rb, ra, n);
00261 lambda=blambda;
00262 Mp->lambda = blambda;
00263
00264 if (dlambda == dlambdamin)
00265 {
00266 if (Mespr==1) fprintf(stdout,"\n increment cannot be decreased dlambda=%e norf %e",dlambda,norf);
00267 break;
00268 }
00269 dlambda/=2.0;
00270 if (dlambda<dlambdamin)
00271 {
00272 dlambda=dlambdamin;
00273
00274 }
00275 if (Mespr==1)
00276 {
00277 fprintf(stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf:%e",dlambda,norf);
00278
00279 }
00280 }
00281 else
00282 {
00283 lambda+=dlambda;
00284 Mm->updateipval();
00285 if (outres == yes)
00286 {
00287 compute_req_val(lcid);
00288 print_step(lcid, i+1, lambda, fa);
00289 print_flush();
00290 }
00291
00292 if ((nlman->check_lambdar == on) && (check_rv_fp == off))
00293 break;
00294
00295 }
00296
00297
00298
00299 }
00300 fprintf(stdout, "\n The total attained lambda =% le\n", lambda);
00301
00302 if ((i < ni-1) && (norf >= ierr))
00303 {
00304 if (Mespr==1)
00305 {
00306 print_err("the step length was decreased to the minimum required value\n"
00307 " but the equilibrium could not be attained.\n"
00308 " FORCED output of the attained results was performed in this step.",
00309 __FILE__, __LINE__, __func__);
00310 }
00311 Mm->updateipval();
00312 compute_req_val (lcid);
00313 print_step_forced(lcid, i+1, lambda, fb);
00314 print_flush();
00315 }
00316
00317 if (Mp->hdbcont.save_stat())
00318 {
00319
00320 if (Mespr==1)
00321 fprintf (stdout,"\n Creating backup file\n");
00322 solver_save (ra, fp, i, lambda, dlambda, NULL, n);
00323 }
00324
00325 delete [] rb;
00326 delete [] dr;
00327 delete [] fi;
00328 delete [] fb;
00329
00330 return lambda;
00331 }
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357 double gnewton_raphson2 (long lcid, nonlinman *nlman, double *ra, double *fa, double *fc, double *fp,
00358 double *flc, double *flp, long li, double ilambda, answertype outres)
00359 {
00360 long i,j,n,ni,ini,modif;
00361 double lambda;
00362 double blambda;
00363 double dlambda;
00364 double dlambdamin;
00365 double dlambdamax;
00366 double zero;
00367 double dtr;
00368 double ierr;
00369 double norf;
00370 double *rb;
00371 double *dr;
00372 double *fi;
00373 double *fb;
00374 matrix lsm_a(3,3);
00375 vector lsm_r(3), lsm_l(3);
00376 flagsw check_rv_fp = nlman->check_lambdar;
00377
00378
00379
00380 n = Ndofm;
00381
00382 zero = Mp->zero;
00383
00384 lambda = ilambda;
00385
00386 ni = nlman->ninr;
00387
00388 ini = nlman->niilnr;
00389
00390 ierr = nlman->errnr;
00391
00392 dlambda=nlman->incrnr;
00393
00394 dlambdamin=nlman->minincrnr;
00395
00396 dlambdamax=nlman->maxincrnr;
00397
00398
00399 rb = new double [n];
00400 fb = new double [n];
00401 fi = new double [n];
00402 dr = new double [n];
00403 memset (rb,0,n*sizeof(double));
00404 memset (fb,0,n*sizeof(double));
00405 memset (fi,0,n*sizeof(double));
00406 memset (dr,0,n*sizeof(double));
00407
00408 modif=0;
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423 if (lambda == 0.0)
00424 assemble_attained_load_vector(fa, flc, flp, n, lambda, nlman);
00425
00426 for (i=li;i<ni;i++)
00427 {
00428
00429
00430 Mp->strainstate=0;
00431
00432
00433 Mp->stressstate=0;
00434
00435
00436 Mp->otherstate=0;
00437
00438 Mp->istep = i;
00439 Mp->jstep = -1;
00440
00441
00442
00443 copyv(ra, rb, n);
00444
00445
00446 blambda=lambda;
00447 Mp->lambda = lambda;
00448
00449 fprintf (stdout,"\n increment %ld lambda %e, dlambda %e",i,lambda,dlambda);
00450
00451
00452
00453 if ((check_rv_fp == on) && (dlambda+lambda > nlman->lambdar))
00454 {
00455 check_rv_fp = off;
00456 dlambda = nlman->lambdar - lambda;
00457 }
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483 addmultv(fa, flp, dlambda, n);
00484
00485
00486
00487 Mp->lambda += dlambda;
00488 Mp->dlambda = dlambda;
00489
00490
00491
00492 cmulv(dlambda, fp, fb, n);
00493
00494
00495
00496
00497
00498 norf = gnewton_raphson_one_step(lcid, nlman, fa, ra, fb, dr, fi, dtr, i, j, li);
00499
00500
00501 if (norf<ierr)
00502 {
00503 lambda+=dlambda;
00504 Mm->updateipval();
00505 if (outres == yes)
00506 {
00507 compute_req_val(lcid);
00508 print_step(lcid, i+1, lambda, fa);
00509 print_flush();
00510 }
00511
00512
00513 if (j < 0)
00514 modif++;
00515 else
00516 modif=0;
00517
00518
00519 if (modif>1)
00520 {
00521 dlambda*=2.0;
00522 if (dlambda>dlambdamax)
00523 dlambda=dlambdamax;
00524 if (Mespr==1)
00525 fprintf (stdout,"\n increment was modified (dlambda increase) dlambda=%e",dlambda);
00526 modif=0;
00527 }
00528
00529
00530 if ((nlman->check_lambdar == on) && (check_rv_fp == off))
00531 break;
00532
00533 continue;
00534 }
00535
00536
00537
00538 if (norf>ierr)
00539 {
00540
00541
00542 addmultv(fa, fp, -dlambda, n);
00543
00544 copyv(rb, ra, n);
00545 lambda=blambda;
00546 Mp->lambda = blambda;
00547
00548 if (dlambda == dlambdamin)
00549 {
00550 if (Mespr==1)
00551 {
00552 fprintf(stdout,"\n load increment cannot be decreased, dlambda=%e norf %e",dlambda,norf);
00553 print_err("the step length was decreased to the minimum required value\n"
00554 " but the equilibrium could not be attained.\n"
00555 " FORCED output of the attained results was performed in this step.",
00556 __FILE__, __LINE__, __func__);
00557 Mm->updateipval();
00558 compute_req_val (lcid);
00559 print_step_forced(lcid, i+1, lambda, fb);
00560 print_flush();
00561 }
00562 break;
00563 }
00564 if (dtr < 1.0)
00565 dlambda *= dtr;
00566 else
00567 dlambda/=2.0;
00568 if (dlambda<dlambdamin)
00569 {
00570 dlambda=dlambdamin;
00571
00572 }
00573 if (Mespr==1)
00574 {
00575 fprintf(stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf:%e",dlambda,norf);
00576
00577 }
00578 }
00579
00580
00581
00582 }
00583
00584 fprintf(stdout, "\n The total attained lambda =% le\n", lambda);
00585
00586 if (Mp->hdbcont.save_stat())
00587 {
00588
00589 if (Mespr==1)
00590 fprintf (stdout,"\n Creating backup file\n");
00591 solver_save (ra, fa, i, lambda, dlambda, NULL, n);
00592 }
00593
00594 delete [] rb;
00595 delete [] dr;
00596 delete [] fi;
00597 delete [] fb;
00598
00599 return lambda;
00600 }
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628 double gnewton_raphson_one_step(long lcid, nonlinman *nlman, double *fa, double *ra, double *fb, double *dr, double *fi,
00629 double &dtr, long istep, long &j, long li)
00630 {
00631 long n = Ndofm;
00632 long ini;
00633 double ierr;
00634 double norf;
00635 double norfa;
00636 matrix lsm_a(3,3);
00637 vector lsm_r(3), lsm_l(3);
00638
00639
00640 ini = nlman->niilnr;
00641
00642 ierr = nlman->errnr;
00643
00644
00645
00646
00647
00648 assemble_stiffness_matrix(lcid,istep,-1,li);
00649
00650
00651 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00652
00653
00654
00655 addv(ra, dr, n);
00656
00657
00658 internal_forces (lcid,fi);
00659
00660
00661
00662
00663
00664 dtr = Mm->dstep_red_mat();
00665 if (dtr < 1.0)
00666 {
00667 j=-1;
00668 return (2.0*ierr);
00669 }
00670
00671
00672
00673 subv(fa, fi, fb, n);
00674
00675 norfa=normv(fa,n);
00676
00677 norf=normv(fb,n);
00678 if (norfa != 0.0)
00679 norf /= norfa;
00680
00681 fprintf (stdout,"\n norm %e",norf);
00682
00683
00684 if (norf<ierr)
00685 {
00686 j = -1;
00687 return norf;
00688 }
00689
00690
00691 fillm(0.0, lsm_a);
00692 fillv(0.0, lsm_r);
00693 for (j=0; j<ini; j++)
00694 {
00695 Mp->jstep = j;
00696
00697
00698 assemble_stiffness_matrix(lcid,istep,j,li);
00699
00700
00701 Mp->ssle->solve_system(Gtm,Smat,dr,fb,Out);
00702
00703
00704 addv(ra, dr, n);
00705
00706
00707
00708
00709 internal_forces(lcid,fi);
00710
00711
00712
00713
00714
00715 dtr = Mm->dstep_red_mat();
00716 if (dtr < 1.0)
00717 {
00718 j=-1;
00719 return (2.0*ierr);
00720 }
00721
00722
00723
00724 subv(fa, fi, fb, n);
00725
00726
00727 norf=normv(fb,n);
00728 if (norfa != 0.0)
00729 norf /= norfa;
00730
00731 fprintf (stdout,"\n increment %ld inner loop j %ld norf %e",istep,j,norf);
00732
00733
00734
00735
00736
00737 if (norf<ierr)
00738
00739 return norf;
00740
00741
00742 if (check_divergency(nlman, lsm_a, lsm_r, lsm_l, j, norf))
00743 return norf;
00744 }
00745
00746
00747 return norf;
00748 }