00001 #include <math.h>
00002 #include <string.h>
00003 #include <stdlib.h>
00004
00005 #include "nssolver.h"
00006 #include "arclength.h"
00007 #include "newtonraph.h"
00008 #include "slopesol.h"
00009
00010 #include "global.h"
00011 #include "globmat.h"
00012 #include "backupsol.h"
00013 #include "mechprint.h"
00014 #include "gmatrix.h"
00015 #include "gtopology.h"
00016 #include "loadcase.h"
00017 #include "dloadcase.h"
00018 #include "intpoints.h"
00019 #include "mathem.h"
00020 #include "vector.h"
00021
00022 #include "meshtransfer.h"
00023 #include "elemswitch.h"
00024
00025 #include "j2flow.h"
00026
00027
00028
00029
00030
00031
00032
00033
00034 void solve_nonlinear_statics ()
00035 {
00036 long i, li;
00037 double *rhs, *ra, *fa, *fp, *fc, *fl, *flc, *flp;
00038 double ilambda, dl;
00039
00040
00041 fa = new double[Ndofm];
00042
00043 fl = new double[2*Ndofm];
00044
00045 switch (Mp->nlman->tnlinsol)
00046 {
00047 case arcl:
00048 case arclrv1:
00049 case arclrv:
00050 case arclg:
00051 seldofinit();
00052 break;
00053 default:
00054 break;
00055 }
00056
00057 Mp->lambda = 0.0;
00058 rhs=Lsrs->give_rhs (0);
00059
00060 for (i=0;i<Lsrs->nlc;i++){
00061
00062
00063
00064 Mp->lambda = 1.0;
00065
00066 mefel_right_hand_side (i, rhs, fl);
00067
00068
00069 ra = Lsrs->give_lhs(i);
00070
00071 memset(ra, 0, sizeof(*ra)*Ndofm);
00072
00073 memset(fa, 0, sizeof(*fa)*Ndofm);
00074
00075 fp = Lsrs->give_rhs(i*2);
00076
00077 fc = Lsrs->give_rhs(i*2+1);
00078
00079 fp = Lsrs->give_rhs(i*2);
00080
00081 flc = fl+Ndofm;
00082
00083 flp = fl;
00084
00085 if ((Mp->nlman->hdbackupal == 2) || (Mp->hdbcont.restore_stat()))
00086 {
00087 if (Mespr==1)
00088 fprintf (stdout,"\n Reading of backup file\n");
00089
00090 solver_restore(ra, fa, li, ilambda, dl, NULL, Ndofm);
00091 Mp->lambda = ilambda;
00092
00093 Mm->initmaterialmodels(i);
00094 print_init(-1, "at");
00095 }
00096 else
00097 {
00098 Mp->lambda = ilambda = 0.0;
00099 li=0;
00100
00101 Mm->initmaterialmodels(i);
00102 print_init(-1, "wt");
00103 print_step(i, li, Mp->lambda, fa);
00104 print_flush();
00105 }
00106
00107 nonlinear_solver(i, ra, fa, fc, fp, flc, flp, ilambda, li);
00108 print_close();
00109
00110 }
00111 delete [] fa;
00112 delete [] fl;
00113 }
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137 void nonlinear_solver (long lcid, double *ra, double *fa, double *fc, double *fp, double *flc, double *flp,
00138 double ilambda, long li)
00139 {
00140 switch (Mp->nlman->tnlinsol){
00141 case arcl:
00142 case arclg:{
00143 garclength(lcid, Mp->nlman, ra, fa, fc, fp, flc, flp, li, ilambda, yes);
00144 break;
00145 }
00146 case newton:
00147 case newtong:{
00148
00149 gnewton_raphson2(lcid, Mp->nlman, ra, fa, fc, fp, flc, flp, li, ilambda, yes);
00150 break;
00151 }
00152 case arclrv1:{
00153
00154 Mp->nlman->check_lambdar = on;
00155 Mp->nlman->lambdar = 1.1;
00156 Mp->nlman->errl = 1.0e-3;
00157 garclength(lcid, Mp->nlman, ra, fa, fc, fp, flc, flp, li, ilambda, yes);
00158 break;
00159 }
00160 case newtonrv1:{
00161
00162 Mp->nlman->check_lambdar = on;
00163 Mp->nlman->errl = 1.0e-3;
00164 gnewton_raphson2(lcid, Mp->nlman, ra, fa, fc, fp, flc, flp, li, ilambda, yes);
00165 break;
00166 }
00167 case newtonrestart:{
00168
00169 newton_raphson_restart (lcid);
00170 break;
00171 }
00172 case displctrl:{
00173
00174 gnewton_raphson2(lcid, Mp->nlman, ra, fa, fc, fp, flc, flp, li, ilambda, yes);
00175 break;
00176 }
00177 case displctrlrv:{
00178
00179 Mp->nlman->check_lambdar = on;
00180 gnewton_raphson2(lcid, Mp->nlman, ra, fa, fc, fp, flc, flp, li, ilambda, yes);
00181 break;
00182 }
00183 case arclrv:{
00184
00185 Mp->nlman->check_lambdar = on;
00186 garclength(lcid, Mp->nlman, ra, fa, fc, fp, flc, flp, li, ilambda, yes);
00187 break;
00188 }
00189 default:
00190 print_err("unknown solver of nonlinear equations system is required", __FILE__, __LINE__, __func__);
00191 }
00192 }
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215 long assemble_stiffness_matrix (long lcid,long i,long j,long li)
00216 {
00217 long ret = 0;
00218
00219 if (Smat == NULL)
00220 {
00221 stiffness_matrix (lcid);
00222 return 1;
00223 }
00224
00225 switch (Mp->nlman->stmat){
00226 case initial_stiff:{
00227 if (i==li && j<0){
00228
00229 stiffness_matrix (lcid);
00230 ret = 1;
00231 }
00232 break;
00233 }
00234 case tangent_stiff:{
00235 stiffness_matrix (lcid);
00236 ret = 1;
00237 break;
00238 }
00239 case secant_stiff:{
00240 stiffness_matrix (lcid);
00241 ret = 1;
00242 break;
00243 }
00244 case incr_tangent_stiff:{
00245 if (j<0){
00246
00247
00248 stiffness_matrix (lcid);
00249 ret = 1;
00250 }
00251 break;
00252 }
00253 case ijth_tangent_stiff:{
00254 if (j<0){
00255
00256 if ((i%Mp->nlman->ithstiffchange == 0) || (i==li))
00257 {
00258 stiffness_matrix (lcid);
00259 ret = 1;
00260 }
00261 }
00262 else{
00263
00264 if (j%Mp->nlman->jthstiffchange == 0)
00265 {
00266 stiffness_matrix (lcid);
00267 ret = 1;
00268 }
00269 }
00270 break;
00271 }
00272 default:{
00273 print_err("unknown type of stiffness matrix is required",__FILE__,__LINE__,__func__);
00274 }
00275 }
00276 return ret;
00277 }
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299 void assemble_attained_load_vector(double *fa, double *fc, double *fp, long n, double lambda, nonlinman *nlman)
00300 {
00301 double aux = lambda;
00302
00303 if ((nlman->check_lambdar == on) && (lambda > nlman->lambdar))
00304 aux = nlman->lambdar;
00305 addmultv(fc, fp, aux, fa, n);
00306 }
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334 long check_divergency(nonlinman *nlman, matrix &lsm_a, vector &lsm_r, vector &lsm_l, long j, double norf)
00335 {
00336 double ret;
00337 double *divc_step = nlman->divc_step;
00338 double *divc_err = nlman->divc_err;
00339 long min_steps = nlman->div_min_steps;
00340 long i;
00341
00342 if (nlman->check_div == on)
00343 {
00344 if (j >= min_steps-1)
00345 {
00346 i = j%min_steps;
00347 if (j == min_steps-1)
00348 ret = lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero, 1);
00349 else
00350 ret = lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, divc_step[i], divc_err[i], Mp->zero, 2);
00351 divc_step[i] = double(j+1);
00352 divc_err[i] = norf;
00353 i = (j-1)%min_steps;
00354
00355 if ((ret > 0.0) && (divc_err[i] < norf))
00356 {
00357 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
00358 return(1);
00359 }
00360 }
00361 else
00362 {
00363
00364 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero, 0);
00365 i = j%min_steps;
00366 divc_step[i] = double(j+1);
00367 divc_err[i] = norf;
00368 }
00369 }
00370 return(0);
00371 }
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 void newton_raphson_restart (long lcid)
00389
00390
00391
00392
00393
00394 {
00395 long i,j,k,n,ni,ini,nii;
00396 double lambda,blambda,dlambda,dlambdamin,zero,err,norf,norfa;
00397 double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp;
00398
00399
00400
00401 n = Ndofm;
00402
00403 ni = Mp->nlman->ninr;
00404
00405 nii = Mp->nlman->nienr;
00406
00407 ini = Mp->nlman->niilnr;
00408
00409 zero = Mp->zero;
00410
00411 err = Mp->nlman->errnr;
00412
00413 dlambda=Mp->nlman->incrnr;
00414
00415 dlambdamin=Mp->nlman->minincrnr;
00416
00417 rb = new double [n];
00418 f = new double [n];
00419 fb = new double [n];
00420 fi = new double [n];
00421 dr = new double [n];
00422 memset (rb,0,n*sizeof(double));
00423 memset (f,0,n*sizeof(double));
00424 memset (fi,0,n*sizeof(double));
00425 memset (dr,0,n*sizeof(double));
00426
00427
00428 r = Lsrs->give_lhs (lcid);
00429 fp = Lsrs->give_rhs (lcid*2);
00430 fc = Lsrs->give_rhs (lcid*2+1);
00431 lambda=0.0;
00432
00433
00434
00435
00436 long ncompstr;
00437 double tmp;
00438 FILE *aux;
00439
00440
00441 if (Mp->nlman->hdbr)
00442 {
00443 aux = fopen (Mp->nlman->backupfname, "rt");
00444 if (aux == NULL)
00445 {
00446 fprintf(stderr, "\n\nError - unable to open backup file %s\n", Mp->nlman->backupfname);
00447 fprintf(stderr, "file %s, line %d\n", __FILE__, __LINE__);
00448 fprintf(stderr, "Program will be terminated\n");
00449 abort();
00450 }
00451
00452 fprintf(stdout, "\n");
00453 solver_restore(r,fc,i,tmp,tmp,NULL,n);
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471 for (i=0; i<Mm->tnip; i++)
00472 {
00473 if (Mm->ip[i].ncompeqother >= 8)
00474 {
00475 ncompstr = Mm->ip[i].ncompstr;
00476 for (j=0; j<ncompstr; j++)
00477 Mm->ip[i].stress[j] = Mm->ip[i].eqother[ncompstr+j];
00478 }
00479 }
00480 }
00481
00482 Mp->strcomp = 0;
00483 internal_forces(lcid, fc);
00484 for (i=0; i < Ndofm; i++)
00485 fc[i] *= 1.0e+6;
00486 Mp->strcomp = 1;
00487
00488
00489
00490
00491
00492
00493 compute_req_val (lcid);
00494 print_init(-1, "wt");
00495 print_step(lcid, 0, 0, fc);
00496 print_flush();
00497
00498
00499
00500 stiffness_matrix (lcid);
00501
00502 norfa=normv(fc,n);
00503
00504 for (j=0;j<nii;j++)
00505 {
00506
00507 internal_forces (lcid,fi);
00508
00509 for (k=0;k<n;k++)
00510 fb[k]=fc[k]-fi[k];
00511
00512 norf=normv(fb,n);
00513 norf /= norfa;
00514
00515 if (norf<err)
00516 {
00517 Mm->updateipval ();
00518 compute_req_val (lcid);
00519 print_step(lcid, j+1, -0.0000001, fc);
00520 print_flush();
00521 break;
00522 }
00523
00524 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554 for(k=0; k<n; k++)
00555 r[k]+=dr[k];
00556
00557
00558
00559
00560
00561 fprintf (stdout,"\n Equilibrium interation loop j %ld norf=%e", j, norf);
00562 }
00563 if ((nii > 0) && (norf > err))
00564 {
00565 fprintf(stdout, "\n Initial equilibrium could not be reached\n");
00566 fprintf(stdout, " Program will be terminated\n");
00567 print_close();
00568 return;
00569 }
00570 else
00571 fprintf(stdout, "\n\n Starting incremental loading . . .\n");
00572
00573
00574
00575
00576
00577 for (i=0;i<ni;i++){
00578
00579
00580 for (j=0;j<n;j++){
00581 rb[j]=r[j];
00582 }
00583
00584 blambda=lambda;
00585
00586 fprintf (stdout,"\n increment %ld lambda %e, dlambda %e",i,lambda,dlambda);
00587
00588
00589 for (j=0;j<n;j++){
00590 f[j]=fc[j]+lambda*fp[j];
00591 fb[j]=dlambda*fp[j];
00592 }
00593
00594
00595 if ((Mp->nlman->stmat==tangent_stiff) || (i == 0))
00596 stiffness_matrix (lcid);
00597
00598
00599
00600
00601 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00602
00603 for (j=0;j<n;j++){
00604 r[j]+=dr[j];
00605 }
00606
00607 internal_forces (lcid,fi);
00608
00609 for (j=0;j<n;j++){
00610 fb[j]=f[j]+dlambda*fp[j];
00611 }
00612 norfa=normv(fb,n);
00613 for (j=0;j<n;j++){
00614 fb[j] -= fi[j];
00615 }
00616
00617 norf=normv(fb,n);
00618 norf /= norfa;
00619
00620
00621
00622
00623
00624 if (norf<err){
00625 lambda+=dlambda;
00626 Mm->updateipval ();
00627 compute_req_val (lcid);
00628 print_step(lcid, i+1, lambda, f);
00629 print_flush();
00630 continue;
00631 }
00632
00633
00634 for (j=0;j<ini;j++){
00635
00636
00637
00638 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00639
00640 for (k=0;k<n;k++){
00641 r[k]+=dr[k];
00642 }
00643
00644
00645 internal_forces (lcid,fi);
00646
00647
00648 for (k=0;k<n;k++){
00649 fb[k]=f[k]+dlambda*fp[k]-fi[k];
00650 }
00651
00652
00653 norf=normv(fb,n);
00654 norf /= norfa;
00655
00656 fprintf (stdout,"\n increment %ld inner loop j %ld norf=%e dlambda %e",i,j,norf,dlambda);
00657
00658
00659 if (norf<err) break;
00660 }
00661
00662 if (j==ini || norf>err){
00663 dlambda/=2.0;
00664 if (dlambda<dlambdamin) dlambda=dlambdamin;
00665
00666 if (Mespr==1){
00667 fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf=%e",dlambda,norf);
00668
00669 }
00670
00671 for (j=0;j<n;j++){
00672 r[j]=rb[j];
00673 }
00674 lambda=blambda;
00675 }
00676 else{
00677 lambda+=dlambda;
00678 Mm->updateipval ();
00679 compute_req_val (lcid);
00680 print_step(lcid, i+1, lambda, f);
00681 print_flush();
00682 }
00683
00684 }
00685
00686 delete [] dr; delete [] fi; delete [] fb; delete [] f;
00687 print_close();
00688 }
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706 void arclsave (long fi,long n,double blambda,double dl,double *r,double *fp)
00707 {
00708 long i,j;
00709 char *file;
00710 FILE *aux;
00711 file = new char[strlen(Mp->path)+11+1];
00712 sprintf (file,"%sarcl.backup",Mp->path);
00713 aux = fopen (file,"w");
00714
00715
00716
00717
00718 fprintf (aux,"%ld\n",fi);
00719
00720 fprintf (aux,"%ld\n",n);
00721
00722 fprintf (aux,"%e\n",blambda);
00723
00724 fprintf (aux,"%e\n",dl);
00725
00726 fprintf (aux,"%ld\n",Mp->nlman->nsdofal);
00727
00728 for (i=0;i<Mp->nlman->nsdofal;i++){
00729 fprintf (aux,"%ld\n",Mp->nlman->seldofal[i]);
00730 }
00731
00732 for (i=0;i<n;i++){
00733 fprintf (aux,"%e\n",r[i]);
00734 }
00735
00736 for (i=0;i<n;i++){
00737 fprintf (aux,"%e\n",fp[i]);
00738 }
00739
00740 fprintf (aux,"%ld\n",Mm->tnip);
00741
00742 for (i=0;i<Mm->tnip;i++){
00743 fprintf (aux,"%ld %ld\n",Mm->ip[i].ncompstr,Mm->ip[i].ncompother);
00744 for (j=0;j<Mm->ip[i].ncompstr;j++){
00745 fprintf (aux,"%e\n",Mm->ip[i].strain[j]);
00746 }
00747 for (j=0;j<Mm->ip[i].ncompstr;j++){
00748 fprintf (aux,"%e\n",Mm->ip[i].stress[j]);
00749 }
00750 for (j=0;j<Mm->ip[i].ncompother;j++){
00751 fprintf (aux,"%e\n",Mm->ip[i].eqother[j]);
00752 }
00753 }
00754
00755 fclose (aux);
00756 }
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774 void arclopen (long &fi,long &n,double &blambda,double &dl,double *r,double *fp)
00775 {
00776 long i,j;
00777 char *file;
00778 FILE *aux;
00779
00780 file = new char[strlen(Mp->path)+11+1];
00781 sprintf (file,"%sarcl.backup",Mp->path);
00782 aux = fopen (file,"r");
00783
00784
00785 fscanf (aux,"%ld",&fi);
00786
00787 fscanf (aux,"%ld",&n);
00788
00789 fscanf (aux,"%le",&blambda);
00790
00791 fscanf (aux,"%le",&dl);
00792
00793 fscanf (aux,"%ld",&Mp->nlman->nsdofal);
00794
00795 for (i=0;i<Mp->nlman->nsdofal;i++){
00796 fscanf (aux,"%ld",&Mp->nlman->seldofal[i]);
00797 }
00798
00799 for (i=0;i<n;i++){
00800 fscanf (aux,"%le",&r[i]);
00801 }
00802
00803 for (i=0;i<n;i++){
00804 fscanf (aux,"%le",&fp[i]);
00805 }
00806
00807 fscanf (aux,"%ld",&Mm->tnip);
00808
00809 for (i=0;i<Mm->tnip;i++){
00810 fscanf (aux,"%ld %ld",&Mm->ip[i].ncompstr,&Mm->ip[i].ncompother);
00811 for (j=0;j<Mm->ip[i].ncompstr;j++){
00812 fscanf (aux,"%le",&Mm->ip[i].strain[j]);
00813 }
00814 for (j=0;j<Mm->ip[i].ncompstr;j++){
00815 fscanf (aux,"%le",&Mm->ip[i].stress[j]);
00816 }
00817 for (j=0;j<Mm->ip[i].ncompother;j++){
00818 fscanf (aux,"%le",&Mm->ip[i].eqother[j]);
00819 }
00820 }
00821
00822 fclose (aux);
00823 }