00001 #include <math.h>
00002 #include <string.h>
00003 #include <stdlib.h>
00004 #include "nssolver.h"
00005 #include "slopesol.h"
00006 #include "global.h"
00007 #include "globmat.h"
00008 #include "mechprint.h"
00009 #include "gmatrix.h"
00010 #include "gtopology.h"
00011 #include "loadcase.h"
00012 #include "dloadcase.h"
00013 #include "intpoints.h"
00014 #include "auxdatout.h"
00015 #include "mechprint.h"
00016 #include "j2flow.h"
00017 #include "chen.h"
00018 #include "hardsoft.h"
00019 #include "mathem.h"
00020 #include "vector.h"
00021 #include "mtsolver.h"
00022
00023 #include "adaptivity.h"
00024 #include "meshtransfer.h"
00025
00026
00027
00028
00029
00030
00031 void solve_nonlinear_statics ()
00032 {
00033 long i;
00034 double *rhs;
00035
00036 Mp->lambda = 1.0;
00037 rhs=Lsrs->give_rhs (0);
00038
00039 for (i=0;i<Lsrs->nlc;i++){
00040
00041
00042
00043
00044 mefel_right_hand_side (i,rhs);
00045
00046 if (!Mp->adaptivity)
00047 nonlinear_solver (i);
00048 else
00049 nonlinear_solver_adaptiv (i);
00050
00051 }
00052
00053 }
00054
00055 void print_output(long lcid)
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 void nonlinear_solver (long lcid)
00084 {
00085 switch (Mp->nlman->tnlinsol){
00086 case arcl:{
00087 arclength (lcid,0);
00088 break;
00089 }
00090 case newton:{
00091 newton_raphson (lcid);
00092 break;
00093 }
00094 case arclrv1:{
00095 arclengthrv1 (lcid,1.1,1.0e-3);
00096 break;
00097 }
00098 case newtonrv1:{
00099 newton_raphsonrv1 (lcid,Mp->nlman->rvlam,1.0e-3);
00100 break;
00101 }
00102 case newtonrestart:{
00103 newton_raphson_restart (lcid);
00104 break;
00105 }
00106 case displctrl:{
00107 displ_control (lcid);
00108 break;
00109 }
00110 case adaptram:{
00111 arclengthadapt (lcid,0);
00112 break;
00113 }
00114 default:{
00115 fprintf (stderr,"\n\n unknown solver of nonlinear equation system is required");
00116 fprintf (stderr,"\n in function nonlinear_solver (%s, line %d).\n",__FILE__,__LINE__);
00117 }
00118 }
00119
00120 }
00121
00122 void refresh_plast ()
00123 {
00124 for (long i=0;i<Mm->tnip;i++)
00125 if ( Mm->j2f[Mm->ip[i].idm[0]].give_consparam (i,0) ){
00126 Mm->plast = 1;
00127 return;
00128 }
00129 }
00130
00131
00132
00133
00134
00135
00136
00137
00138 void nonlinear_solver_adaptiv (long lcid)
00139 {
00140 int width;
00141 long i,adaptcontrol;
00142 char file[255];
00143
00144 switch (Mp->nlman->tnlinsol){
00145 case arcl:{
00146 if (Mespr==1){
00147 fprintf (stdout,"\n\n ----------------------------");
00148 fprintf (stdout, "\n *** Arc-length - start ***\n");
00149 }
00150
00151 width = 1 + (long)log10((double)Mp->nlman->nial);
00152 adaptcontrol = 1;
00153
00154 for (i=0; i < Mp->nlman->nial ;){
00155
00156 i = arclength (lcid,adaptcontrol);
00157
00158 if ( i<Mp->nlman->nial || Mp->nlman->nial==1 && (Mp->adaptflag & 16) ){
00159 if (Mespr==1) fprintf (stdout,"\n\n NEW MESH will be generated\n");
00160
00161 newmeshgen (width,i);
00162 sprintf (file,"%s%s.in.%0*ld",Mp->path,Mp->filename,width,i);
00163 newmeshread (file,lcid);
00164
00165 if ( i<Mp->nlman->nial ){
00166 if (Mm->plast == 0)
00167 refresh_plast ();
00168
00169 arclength15 (lcid);
00170 }
00171 else
00172 Ada->run (4,width,i);
00173 }
00174
00175
00176
00177 adaptcontrol++;
00178 }
00179
00180 if (Mespr==1){
00181 fprintf (stdout,"\n\n *** Arc-length is finished ***");
00182 fprintf (stdout, "\n *********************************");
00183 }
00184 break;
00185 }
00186 default:{
00187 fprintf (stderr,"\n\n unknown solver of nonlinear equation system is required");
00188 fprintf (stderr,"\n in function nonlinear_solver_adaptiv (%s, line %d).\n",__FILE__,__LINE__);
00189 }
00190 }
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210 long arclength (long lcid,long adaptcontrol)
00211 {
00212 long i,j,k,n,ni,ini,stop,modif,li,newmesh, numr;
00213 double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr;
00214 double lambdao,ss1,ss2,ss3,ss4,ss5;
00215 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi;
00216
00217
00218
00219
00220 n = Ndofm;
00221
00222 ni = Mp->nlman->nial;
00223
00224 ini = Mp->nlman->niilal;
00225
00226 zero = Mp->zero;
00227
00228 ierr = Mp->nlman->erral;
00229
00230 dl = Mp->nlman->dlal;
00231
00232 dlmax = Mp->nlman->dlmaxal;
00233
00234 dlmin = Mp->nlman->dlminal;
00235
00236 psi = Mp->nlman->psial;
00237
00238 dlambda=0.0;
00239
00240
00241 r = new double [n];
00242 ddr = new double [n];
00243 u = new double [n];
00244 v = new double [n];
00245 f = new double [n];
00246 fa = new double [n];
00247 fi = new double [n];
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260 ra = Lsrs->give_lhs (lcid);
00261 fp = Lsrs->give_rhs (lcid*2);
00262 fc = Lsrs->give_rhs (lcid*2+1);
00263
00264 if (Mp->nlman->hdbackupal == 2){
00265 arclopen (li,n,lambda,dl,ra,fp);
00266 }
00267 else if (adaptcontrol >= 2){
00268 arclopen_adapt (li,lambda,blambda,dl);
00269 }
00270 else{
00271 lambda=0.0;
00272 lambdao=0.0;
00273 li=0;
00274 }
00275
00276
00277
00278
00279
00280
00281 norfp = ss(fp,fp,n);
00282 modif=0;
00283
00284
00285
00286 for (j=0;j<n;j++){
00287 fa[j]=fc[j]+(lambda+dlambda)*fp[j];
00288 }
00289 if (li)
00290 print_init(-1, "at");
00291 else
00292 {
00293 print_init(-1, "wt");
00294 print_step(lcid, li, lambda, fa);
00295 }
00296
00297
00298
00299
00300
00301
00302 for (i=li;i<ni;i++){
00303
00304 fprintf (Out,"\n\n arc-length prirustek %ld",i);
00305
00306
00307
00308
00309 stop=1;
00310
00311 copyv (ra, r, n);
00312
00313 blambda=lambda;
00314
00315
00316
00317
00318
00319
00320
00321
00322 if ((Mp->nlman->stmat==tangent_stiff) || (i == li)){
00323 stiffness_matrix (lcid);
00324 Smat->printmat (Out);
00325 }
00326
00327
00328 copyv(fp, f, n);
00329
00330
00331 Smat->solve_system (v,f);
00332
00333
00334 norv = displincr (v,n);
00335
00336
00337 dlambda = dl/sqrt(norv+psi*psi*norfp);
00338
00339
00340
00341
00342 for (j=0;j<n;j++){
00343 ddr[j]=dlambda*v[j];
00344 ra[j]+=ddr[j];
00345 fa[j]=fc[j]+(lambda+dlambda)*fp[j];
00346 }
00347 ddlambda=dlambda;
00348
00349
00350 internal_forces (lcid,fi);
00351 subv(fa, fi, f, n);
00352 norf=sizev(f,n);
00353 norfa = sizev(fa, n);
00354
00355
00356 norf /= norfa;
00357
00358
00359
00360 if (Mespr==1) fprintf (stdout,"\n increment %ld norv %e norf %e",i,norv,norf);
00361
00362 if (norf<ierr){
00363
00364
00365
00366 modif++;
00367
00368 if (modif>1){
00369
00370 dl*=2.0;
00371 if (dl>dlmax)
00372 dl=dlmax;
00373 modif=0;
00374 if (Mespr==1){
00375 fprintf (stdout,"\n arc-length must be modified (dl increase) dl=%e",dl);
00376
00377 }
00378 }
00379 lambda+=dlambda;
00380 Mm->updateipval ();
00381 compute_req_val (lcid);
00382 print_step(lcid, i+1, lambda, fa);
00383 print_flush();
00384 }
00385 else{
00386
00387
00388
00389 stop=0;
00390 for (j=0;j<ini;j++){
00391 if (Mp->stmat==tangent_stiff)
00392 stiffness_matrix (lcid);
00393
00394 fprintf (Out,"\n arc-length vnitrni smycka %ld %ld",i,j);
00395
00396
00397 Smat->solve_system (u,f);
00398
00399 copyv(ddr, f, n);
00400 addv(ddr, u, n);
00401
00402 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2);
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420 numr = solv_polynom_2(a2, a1, a0, l1, l2);
00421 switch (numr)
00422 {
00423 case -1:
00424 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength");
00425 break;
00426 case 0:
00427 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength");
00428 break;
00429 case 1:
00430 dlambda = l1;
00431 break;
00432 default:
00433 break;
00434 }
00435 ss1=0.0; ss2=0.0;
00436 ss3=0.0; ss4=0.0; ss5=0.0;
00437 for (k=0;k<n;k++){
00438 ss1+=(ddr[k]+l1*v[k])*f[k];
00439 ss2+=(ddr[k]+l2*v[k])*f[k];
00440 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]);
00441 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]);
00442 ss5+=f[k]*f[k];
00443 }
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453 if (fabs(l1)>fabs(l2)) dlambda=l2;
00454 else dlambda=l1;
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478 for (k=0;k<n;k++){
00479 ddr[k]+=dlambda*v[k];
00480 ra[k]+=u[k]+dlambda*v[k];
00481 fa[k]+=dlambda*fp[k];
00482 }
00483 ddlambda+=dlambda;
00484
00485
00486
00487
00488
00489 internal_forces (lcid,fi);
00490 subv(fa, fi, f, n);
00491 norf=sizev(f,n);
00492 norfa = sizev(fa, n);
00493 norf /= norfa;
00494
00495
00496
00497 if (Mespr==1){
00498 fprintf (stdout,"\n norf=%e ierr=%e",norf,ierr);
00499
00500 }
00501
00502 if (norf<ierr){
00503 lambda+=ddlambda;
00504 Mm->updateipval ();
00505 stop=1;
00506 compute_req_val (lcid);
00507 print_step(lcid, i+1, lambda, fa);
00508 print_flush();
00509 break;
00510 }
00511 }
00512 modif=0;
00513 if (stop==0){
00514
00515 dl/=2.0;
00516 if (dl<dlmin){
00517 dl=dlmin;
00518 break;
00519 }
00520 if (Mespr==1){
00521 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e",dl);
00522
00523 }
00524
00525 copyv(r, ra, n);
00526
00527 lambda=blambda;
00528 }
00529 }
00530
00531 fprintf (stdout,"\n increment %ld total lambda %e",i,lambda);
00532
00533 if (stop==0)
00534 continue;
00535
00536 if (Mp->adaptflag & 16 && i%10 && Mm->plast)
00537 continue;
00538
00539 newmesh = 0;
00540 if (Mp->adaptivity && Mp->nlman->nial > (i+1) && !( (Mp->adaptflag & 16) && (Mp->nlman->nial==100 || Mp->nlman->nial==1))){
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550 long width = 1 + (long)log10((double)Mp->nlman->nial);
00551 newmesh = Ada->run (1,width,i);
00552 }
00553
00554 if (newmesh){
00555 i++;
00556 break;
00557 }
00558
00559
00560 }
00561
00562
00563
00564
00565
00566 fprintf (Out,"\n\n blb %lf",Lsrs->lhs[0]);
00567
00568
00569 if (Mp->nlman->hdbackupal==1)
00570 arclsave (i,n,lambda,dl,ra,fp);
00571
00572
00573 if (Mp->adaptivity)
00574 if ((Mp->adaptflag & 16) && (Mp->nlman->nial==100 || Mp->nlman->nial==1))
00575 Ada->run (0,0,0);
00576 else
00577 if (i < ni)
00578 arclsave_adapt (i,lambda,0,dl);
00579
00580
00581
00582 print_close();
00583 delete [] r;
00584 delete [] fi;
00585 delete [] fa;
00586 delete [] f;
00587 delete [] v;
00588 delete [] u;
00589 delete [] ddr;
00590
00591 return (i);
00592 }
00593
00594 long arclength15 (long lcid)
00595 {
00596 long j,k,n,ini,stop,li;
00597 double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,bdlambda,dlambda,ddlambda,norf,norfp,norv,zero,ierr;
00598 double ss1,ss2,ss3,ss4,ss5;
00599 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi;
00600
00601
00602
00603
00604
00605
00606
00607 n = Ndofm;
00608 ini = Mp->nlman->niilal;
00609 zero = Mp->zero;
00610 ierr = Mp->nlman->erral;
00611 dl = Mp->nlman->dlal;
00612 dlmax = Mp->nlman->dlmaxal;
00613 dlmin = Mp->nlman->dlminal;
00614 psi = Mp->nlman->psial;
00615
00616 r = new double [n]; ddr = new double [n]; u = new double [n]; v = new double [n];
00617 f = new double [n]; fa = new double [n]; fi = new double [n];
00618
00619
00620 ra = Lsrs->give_lhs (lcid);
00621 fp = Lsrs->give_rhs (lcid*2);
00622 fc = Lsrs->give_rhs (lcid*2+1);
00623
00624 arclopen_adapt (li,lambda,dlambda,dl);
00625
00626
00627
00628
00629 fprintf (stdout,"\n arc-length15: lambda %e dl %e",lambda,dl);
00630
00631 stiffness_matrix (lcid);
00632 for (j=0;j<n;j++) f[j]=fp[j];
00633
00634 Smat->solve_system (v,f);
00635 norv = displincr (v,n);
00636 norfp = ss (fp,fp,n);
00637
00638 if (0 && !Mm->plast){
00639 for (j=0;j<n;j++)
00640 ra[j]=(lambda+dlambda)*v[j];
00641
00642
00643 }
00644 else{
00645 stop = 0;
00646
00647 for (j=0;j<n;j++){ r[j]=ra[j]; }
00648 blambda=lambda;
00649 bdlambda=dlambda;
00650
00651
00652 for (j=0;j<n;j++){
00653 ddr[j]=dlambda*v[j];
00654 fa[j]=fc[j]+(lambda+dlambda)*fp[j];
00655 }
00656 ddlambda=dlambda;
00657
00658 internal_forces (lcid,fi);
00659
00660 for (k=0;k<n;k++)
00661 f[k]=fa[k]-fi[k];
00662
00663 norf=ss(f,f,n);
00664 norf=sqrt(norf); if (Mespr==1) fprintf (stdout,"\n %e %e norf %e",lambda,dl,norf);
00665
00666 if (norf<ierr){
00667 lambda+=dlambda;
00668 Mm->updateipval ();
00669 stop = 1;
00670 }
00671
00672 while (!stop){
00673
00674
00675
00676
00677
00678
00679 for (j=0;j<10*ini;j++){
00680
00681
00682
00683 Smat->solve_system (u,f);
00684
00685 for (k=0;k<n;k++){
00686 f[k]=ddr[k];
00687 ddr[k]+=u[k];
00688 }
00689
00690
00691
00692
00693
00694
00695 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2);
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739 if (ss1/ss3/ss5>ss2/ss4/ss5) dlambda=l1;
00740 else dlambda=l2;
00741
00742
00743 dlambda = 0;
00744
00745 for (k=0;k<n;k++){
00746 ddr[k]+=dlambda*v[k];
00747 ra[k]+=u[k]+dlambda*v[k];
00748 fa[k]+=dlambda*fp[k];
00749 }
00750 ddlambda+=dlambda;
00751
00752 internal_forces (lcid,fi);
00753
00754 for (k=0;k<n;k++)
00755 f[k]=fa[k]-fi[k];
00756
00757 norf=ss(f,f,n); norf=sqrt(norf);
00758
00759 if (Mespr==1){
00760 fprintf (stdout,"\n increment %ld inner loop %ld norf=%e",li,j,norf);
00761
00762 }
00763
00764 if (norf<0.10*ierr){
00765 lambda+=ddlambda;
00766
00767 Mm->updateipval ();
00768 stop=1; break;
00769 }
00770 }
00771
00772 if (stop==0){
00773
00774 dl/=2.0;
00775 if (dl<dlmin){ dl=dlmin; break; }
00776
00777 if (Mespr==1){
00778 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf);
00779
00780 }
00781
00782 for (j=0;j<n;j++){
00783 ra[j]=r[j];
00784 }
00785
00786 lambda=blambda;
00787 ddlambda=bdlambda;
00788
00789
00790
00791
00792 stop=1;
00793
00794
00795 }
00796 }
00797 }
00798
00799
00800
00801
00802
00803 arclsave_adapt (li,lambda,0.0,dl);
00804
00805 delete [] r;
00806 delete [] fi;
00807 delete [] fa;
00808 delete [] f;
00809 delete [] v;
00810 delete [] u;
00811 delete [] ddr;
00812
00813
00814
00815 return (li);
00816 }
00817
00818 void newton_raphson (long lcid)
00819
00820
00821
00822
00823
00824 {
00825 long i,j,k,n,ni,ini;
00826 double lambda,blambda,dlambda,dlambdamin,zero,err,norf,norfa;
00827 double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp;
00828
00829
00830
00831 n = Ndofm;
00832
00833 ni = Mp->nlman->ninr;
00834
00835 ini = Mp->nlman->niilnr;
00836
00837 zero = Mp->zero;
00838
00839 err = Mp->nlman->errnr;
00840
00841 dlambda=Mp->nlman->incrnr;
00842
00843 dlambdamin=Mp->nlman->minincrnr;
00844
00845 rb = new double [n];
00846 f = new double [n];
00847 fb = new double [n];
00848 fi = new double [n];
00849 dr = new double [n];
00850 memset (rb,0,n*sizeof(double));
00851 memset (f,0,n*sizeof(double));
00852 memset (fi,0,n*sizeof(double));
00853 memset (dr,0,n*sizeof(double));
00854
00855
00856 r = Lsrs->give_lhs (lcid);
00857 fp = Lsrs->give_rhs (lcid*2);
00858 fc = Lsrs->give_rhs (lcid*2+1);
00859 lambda=0.0;
00860
00861
00862
00863
00864 print_init(-1, "wt");
00865 print_step(lcid, 0, 0.0, f);
00866 print_flush();
00867 for (i=0;i<ni;i++){
00868
00869
00870 for (j=0;j<n;j++){
00871 rb[j]=r[j];
00872 }
00873
00874 blambda=lambda;
00875
00876 fprintf (stdout,"\n increment %ld lambda %e, dlambda %e",i,lambda,dlambda);
00877
00878
00879 for (j=0;j<n;j++){
00880 f[j]=fc[j]+lambda*fp[j];
00881 fb[j]=dlambda*fp[j];
00882 }
00883
00884
00885 if ((Mp->stmat==tangent_stiff) || (i == 0))
00886 stiffness_matrix (lcid);
00887
00888
00889 Smat->solve_system (dr,fb);
00890
00891 for (j=0;j<n;j++){
00892 r[j]+=dr[j];
00893 }
00894
00895
00896 internal_forces (lcid,fi);
00897
00898
00899 for (j=0;j<n;j++){
00900 fb[j]=f[j]+dlambda*fp[j];
00901 }
00902 norfa=ss(fb,fb,n);
00903 for (j=0;j<n;j++){
00904 fb[j] -= fi[j];
00905 }
00906
00907 norf=ss(fb,fb,n);
00908 norf /= norfa;
00909
00910
00911
00912 if (norf<err){
00913 lambda+=dlambda;
00914 Mm->updateipval ();
00915 compute_req_val (lcid);
00916 print_step(lcid, i+1, lambda, f);
00917 print_flush();
00918 continue;
00919 }
00920
00921
00922 for (j=0;j<ini;j++){
00923
00924
00925 Smat->solve_system (dr,fb);
00926
00927 for (k=0;k<n;k++){
00928 r[k]+=dr[k];
00929 }
00930
00931
00932 internal_forces (lcid,fi);
00933
00934
00935 for (k=0;k<n;k++){
00936 fb[k]=f[k]+dlambda*fp[k]-fi[k];
00937 }
00938
00939
00940 norf=ss(fb,fb,n);
00941 norf /= norfa;
00942
00943 fprintf (stdout,"\n increment %ld inner loop j %ld norf=%e dlambda %e",i,j,norf,dlambda);
00944
00945
00946 if (norf<err) break;
00947 }
00948
00949 if (j==ini || norf>err){
00950 dlambda/=2.0;
00951 if (dlambda<dlambdamin) dlambda=dlambdamin;
00952
00953 if (Mespr==1){
00954 fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf=%e",dlambda,norf);
00955
00956 }
00957
00958 for (j=0;j<n;j++){
00959 r[j]=rb[j];
00960 }
00961 lambda=blambda;
00962 }
00963 else{
00964 lambda+=dlambda;
00965 Mm->updateipval ();
00966 compute_req_val (lcid);
00967 print_step(lcid, i+1, lambda, f);
00968 print_flush();
00969 }
00970
00971 }
00972
00973 delete [] dr; delete [] fi; delete [] fb; delete [] f;
00974 print_close();
00975 }
00976
00977 void displ_control (long lcid)
00978
00979
00980
00981
00982
00983
00984
00985 {
00986 long i,j,k,n,ni,ini;
00987 double lambda,blambda,dlambda,dlambdamin,zero,err,norf,norfa;
00988 double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp, *fir2;
00989
00990
00991
00992 n = Ndofm;
00993
00994 ni = Mp->nlman->ninr;
00995
00996 ini = Mp->nlman->niilnr;
00997
00998 zero = Mp->zero;
00999
01000 err = Mp->nlman->errnr;
01001
01002 dlambda=Mp->nlman->incrnr;
01003
01004 dlambdamin=Mp->nlman->minincrnr;
01005
01006 rb = new double [n];
01007 f = new double [n];
01008 fb = new double [n];
01009 fi = new double [n];
01010 fir2 = new double[n];
01011 dr = new double [n];
01012 memset (rb,0,n*sizeof(double));
01013 memset (f,0,n*sizeof(double));
01014 memset (fi,0,n*sizeof(double));
01015 memset (dr,0,n*sizeof(double));
01016
01017
01018 r = Lsrs->give_lhs (lcid);
01019 fp = Lsrs->give_rhs (lcid*2);
01020 fc = Lsrs->give_rhs (lcid*2+1);
01021 lambda=0.0;
01022 Mp->lambda = lambda;
01023
01024
01025
01026
01027 print_init(-1, "wt");
01028 print_step(lcid, 0, 0.0, f);
01029 print_flush();
01030 for (i=0;i<ni;i++){
01031
01032
01033 for (j=0;j<n;j++){
01034 rb[j]=r[j];
01035 }
01036
01037 blambda=lambda;
01038 Mp->lambda = lambda;
01039
01040 fprintf (stdout,"\n increment %ld lambda %e, dlambda %e",i,lambda,dlambda);
01041
01042
01043 Mp->lambda += dlambda;
01044 internal_forces(lcid, fir2);
01045
01046 for (j=0;j<n;j++){
01047 f[j]=fc[j]+lambda*fp[j];
01048 fb[j]= f[j] + dlambda*fp[j] - fir2[j];
01049 }
01050
01051
01052 if ((Mp->stmat==tangent_stiff) || (i == 0))
01053 stiffness_matrix (lcid);
01054
01055
01056 Smat->solve_system (dr,fb);
01057
01058 for (j=0;j<n;j++){
01059 r[j]+=dr[j];
01060 }
01061
01062
01063 internal_forces (lcid,fi);
01064
01065
01066 for (j=0;j<n;j++){
01067 fb[j]=f[j]+dlambda*fp[j];
01068 }
01069 norfa=ss(fb,fb,n);
01070 for (j=0;j<n;j++){
01071 fb[j] -= fi[j];
01072 }
01073
01074 norf=ss(fb,fb,n);
01075 if (norfa != 0.0)
01076 norf /= norfa;
01077
01078 if (norf<err){
01079 lambda+=dlambda;
01080 Mm->updateipval ();
01081 compute_req_val (lcid);
01082 print_step(lcid, i+1, lambda, f);
01083 print_flush();
01084 continue;
01085 }
01086
01087
01088 for (j=0;j<ini;j++){
01089
01090
01091 Smat->solve_system (dr,fb);
01092
01093 for (k=0;k<n;k++){
01094 r[k]+=dr[k];
01095 }
01096
01097
01098 internal_forces (lcid,fi);
01099
01100
01101 for (k=0;k<n;k++){
01102 fb[k]=f[k]+dlambda*fp[k]-fi[k];
01103 }
01104
01105
01106 norf=ss(fb,fb,n);
01107 if (norfa != 0.0)
01108 norf /= norfa;
01109
01110 fprintf (stdout,"\n increment %ld inner loop j %ld norf=%e dlambda %e",i,j,norf,dlambda);
01111
01112 if (norf<err) break;
01113 }
01114
01115 if (j==ini || norf>err)
01116 {
01117 if (dlambda == dlambdamin)
01118 {
01119 if (Mespr==1) fprintf (stdout,"\n increment cannot be decreased dlambda=%e norf=%e",dlambda,norf);
01120 break;
01121 }
01122 dlambda/=2.0;
01123 if (dlambda<dlambdamin) dlambda=dlambdamin;
01124 if (Mespr==1) fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf=%e",dlambda,norf);
01125
01126 for (j=0;j<n;j++)
01127 r[j]=rb[j];
01128
01129 lambda=blambda;
01130 }
01131 else
01132 {
01133 lambda+=dlambda;
01134 Mm->updateipval ();
01135 compute_req_val (lcid);
01136 print_step(lcid, i+1, lambda, f);
01137 print_flush();
01138 }
01139 }
01140
01141 delete [] dr; delete [] fi; delete [] fb; delete [] f;
01142 print_close();
01143 }
01144
01145 void newton_raphson_restart (long lcid)
01146
01147
01148
01149
01150
01151 {
01152 long i,j,k,n,ni,ini,nii;
01153 double lambda,blambda,dlambda,dlambdamin,zero,err,norf,norfa;
01154 double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp;
01155
01156
01157
01158 n = Ndofm;
01159
01160 ni = Mp->nlman->ninr;
01161
01162 nii = Mp->nlman->nienr;
01163
01164 ini = Mp->nlman->niilnr;
01165
01166 zero = Mp->zero;
01167
01168 err = Mp->nlman->errnr;
01169
01170 dlambda=Mp->nlman->incrnr;
01171
01172 dlambdamin=Mp->nlman->minincrnr;
01173
01174 rb = new double [n];
01175 f = new double [n];
01176 fb = new double [n];
01177 fi = new double [n];
01178 dr = new double [n];
01179 memset (rb,0,n*sizeof(double));
01180 memset (f,0,n*sizeof(double));
01181 memset (fi,0,n*sizeof(double));
01182 memset (dr,0,n*sizeof(double));
01183
01184
01185 r = Lsrs->give_lhs (lcid);
01186 fp = Lsrs->give_rhs (lcid*2);
01187 fc = Lsrs->give_rhs (lcid*2+1);
01188 lambda=0.0;
01189
01190
01191
01192
01193 long ncompstr;
01194 double tmp;
01195 FILE *aux;
01196
01197
01198 if (Mp->nlman->hdbr)
01199 {
01200 aux = fopen (Mp->nlman->backupfname, "rt");
01201 if (aux == NULL)
01202 {
01203 fprintf(stderr, "\n\nError - unable to open backup file %s\n", Mp->nlman->backupfname);
01204 fprintf(stderr, "file %s, line %d\n", __FILE__, __LINE__);
01205 fprintf(stderr, "Program will be terminated\n");
01206 abort();
01207 }
01208
01209 fprintf(stdout, "\n");
01210 for (k=0; k<Mp->nlman->hdbid; k++)
01211 {
01212 fprintf(stdout, "\r reading backup record %ld", k+1);
01213 fflush(stdout);
01214 mtsolver_restore(aux, r, fc, i, tmp, j, 39, 8);
01215 }
01216 fprintf(stdout, "\n backup time: %e\n", tmp);
01217 if (j != Ndofm)
01218 {
01219 fprintf(stderr, "\n\nError - incorrect number of DOFs is stored in backup file\n");
01220 fprintf(stderr, "file %s, line %d\n", __FILE__, __LINE__);
01221 }
01222
01223
01224
01225 for (i=0; i<Mm->tnip; i++)
01226 {
01227 if (Mm->ip[i].ncompeqother >= 8)
01228 {
01229 ncompstr = Mm->ip[i].ncompstr;
01230 for (j=0; j<ncompstr; j++)
01231 Mm->ip[i].stress[j] = Mm->ip[i].eqother[ncompstr+j];
01232 }
01233 }
01234 }
01235
01236 Mp->strcomp = 0;
01237 internal_forces(lcid, fc);
01238 for (i=0; i < Ndofm; i++)
01239 fc[i] *= 1.0e+6;
01240 Mp->strcomp = 1;
01241
01242
01243
01244
01245
01246
01247 compute_req_val (lcid);
01248 print_init(-1, "wt");
01249 print_step(lcid, 0, 0, fc);
01250 print_flush();
01251
01252
01253
01254 stiffness_matrix (lcid);
01255
01256 norfa=ss(fc,fc,n);
01257
01258 for (j=0;j<nii;j++)
01259 {
01260
01261 internal_forces (lcid,fi);
01262
01263 for (k=0;k<n;k++)
01264 fb[k]=fc[k]-fi[k];
01265
01266 norf=ss(fb,fb,n);
01267 norf /= norfa;
01268
01269 if (norf<err)
01270 {
01271 Mm->updateipval ();
01272 compute_req_val (lcid);
01273 print_step(lcid, j+1, -0.0000001, fc);
01274 print_flush();
01275 break;
01276 }
01277 Smat->solve_system (dr,fb);
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306 for(k=0; k<n; k++)
01307 r[k]+=dr[k];
01308
01309
01310
01311
01312
01313 fprintf (stdout,"\n Equilibrium interation loop j %ld norf=%e", j, norf);
01314 }
01315 if ((nii > 0) && (norf > err))
01316 {
01317 fprintf(stdout, "\n Initial equilibrium could not be reached\n");
01318 fprintf(stdout, " Program will be terminated\n");
01319 print_close();
01320 return;
01321 }
01322 else
01323 fprintf(stdout, "\n\n Starting incremental loading . . .\n");
01324
01325
01326
01327
01328
01329 for (i=0;i<ni;i++){
01330
01331
01332 for (j=0;j<n;j++){
01333 rb[j]=r[j];
01334 }
01335
01336 blambda=lambda;
01337
01338 fprintf (stdout,"\n increment %ld lambda %e, dlambda %e",i,lambda,dlambda);
01339
01340
01341 for (j=0;j<n;j++){
01342 f[j]=fc[j]+lambda*fp[j];
01343 fb[j]=dlambda*fp[j];
01344 }
01345
01346
01347 if ((Mp->stmat==tangent_stiff) || (i == 0))
01348 stiffness_matrix (lcid);
01349
01350
01351
01352 Smat->solve_system (dr,fb);
01353
01354 for (j=0;j<n;j++){
01355 r[j]+=dr[j];
01356 }
01357
01358 internal_forces (lcid,fi);
01359
01360 for (j=0;j<n;j++){
01361 fb[j]=f[j]+dlambda*fp[j];
01362 }
01363 norfa=ss(fb,fb,n);
01364 for (j=0;j<n;j++){
01365 fb[j] -= fi[j];
01366 }
01367
01368 norf=ss(fb,fb,n);
01369 norf /= norfa;
01370
01371
01372
01373
01374
01375 if (norf<err){
01376 lambda+=dlambda;
01377 Mm->updateipval ();
01378 compute_req_val (lcid);
01379 print_step(lcid, i+1, lambda, f);
01380 print_flush();
01381 continue;
01382 }
01383
01384
01385 for (j=0;j<ini;j++){
01386
01387
01388 Smat->solve_system (dr,fb);
01389
01390 for (k=0;k<n;k++){
01391 r[k]+=dr[k];
01392 }
01393
01394
01395 internal_forces (lcid,fi);
01396
01397
01398 for (k=0;k<n;k++){
01399 fb[k]=f[k]+dlambda*fp[k]-fi[k];
01400 }
01401
01402
01403 norf=ss(fb,fb,n);
01404 norf /= norfa;
01405
01406 fprintf (stdout,"\n increment %ld inner loop j %ld norf=%e dlambda %e",i,j,norf,dlambda);
01407
01408
01409 if (norf<err) break;
01410 }
01411
01412 if (j==ini || norf>err){
01413 dlambda/=2.0;
01414 if (dlambda<dlambdamin) dlambda=dlambdamin;
01415
01416 if (Mespr==1){
01417 fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf=%e",dlambda,norf);
01418
01419 }
01420
01421 for (j=0;j<n;j++){
01422 r[j]=rb[j];
01423 }
01424 lambda=blambda;
01425 }
01426 else{
01427 lambda+=dlambda;
01428 Mm->updateipval ();
01429 compute_req_val (lcid);
01430 print_step(lcid, i+1, lambda, f);
01431 print_flush();
01432 }
01433
01434 }
01435
01436 delete [] dr; delete [] fi; delete [] fb; delete [] f;
01437 print_close();
01438 }
01439
01440 void seldofinit ()
01441 {
01442 long i,j,k,l,ndofn;
01443 double x1,x2,y1,y2,z1,z2,length;
01444
01445 switch (Mp->nlman->displnorm){
01446 case alldofs:{ break; }
01447 case seldofs:
01448 case seldofscoord:{
01449 for (i=0;i<Mp->nlman->nsdofal;i++){
01450 j=Mp->nlman->seldofal[i];
01451 Mp->nlman->seldofal[i]=Mt->give_dof (Mp->nlman->selnodal[i],j)-1;
01452 }
01453 break;
01454 }
01455 case selnodes:{
01456 Mp->nlman->nsdofal=0;
01457 for (i=0;i<Mp->nlman->nsnal;i++){
01458 Mp->nlman->nsdofal+=Mt->give_ndofn (Mp->nlman->selnodal[i]);
01459 }
01460 Mp->nlman->seldofal = new long [Mp->nlman->nsdofal];
01461 k=0;
01462 for (i=0;i<Mp->nlman->nsnal;i++){
01463 ndofn=Mt->give_ndofn (Mp->nlman->selnodal[i]);
01464 for (j=0;j<ndofn;j++){
01465 l=Mt->give_dof (Mp->nlman->selnodal[i],j);
01466 if (l>0){
01467 Mp->nlman->seldofal[k]=l-1;
01468 k++;
01469 }
01470 }
01471 }
01472 Mp->nlman->nsdofal=k;
01473 break;
01474 }
01475 case nodesdistincr:{
01476 Mp->nlman->nsdofal=0;
01477 for (i=0;i<Mp->nlman->nsnal;i++){
01478 Mp->nlman->nsdofal+=Mt->give_ndofn (Mp->nlman->selnodal[i]);
01479 }
01480 Mp->nlman->seldofal = new long [Mp->nlman->nsdofal];
01481 k=0;
01482 for (i=0;i<Mp->nlman->nsnal;i++){
01483 ndofn=Mt->give_ndofn (Mp->nlman->selnodal[i]);
01484 for (j=0;j<ndofn;j++){
01485 Mp->nlman->seldofal[k]=Mt->give_dof (Mp->nlman->selnodal[i],j);
01486 k++;
01487 }
01488 }
01489 x1=Gtm->gnodes[Mp->nlman->selnodal[0]].x; x2=Gtm->gnodes[Mp->nlman->selnodal[1]].x;
01490 y1=Gtm->gnodes[Mp->nlman->selnodal[0]].y; y2=Gtm->gnodes[Mp->nlman->selnodal[1]].y;
01491 z1=Gtm->gnodes[Mp->nlman->selnodal[0]].z; z2=Gtm->gnodes[Mp->nlman->selnodal[1]].z;
01492 length=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
01493 Mp->nlman->nxal=(x2-x1)/length; Mp->nlman->nyal=(y2-y1)/length; Mp->nlman->nzal=(z2-z1)/length;
01494 break;
01495 }
01496 default:{
01497 fprintf (stderr,"\n\n unknown displacement norm is required in function seldofinit (%s, line %d)",__FILE__,__LINE__);
01498 }
01499 }
01500 }
01501
01502 double displincr (double *displincr,long n)
01503 {
01504 long i;
01505 double norm,u1,u2,v1,v2,w1,w2;
01506
01507 switch (Mp->nlman->displnorm){
01508 case alldofs:{
01509 norm = ss (displincr,displincr,n);
01510 break;
01511 }
01512 case seldofs:
01513 case seldofscoord:
01514 case selnodes:{
01515 norm=0.0;
01516 for (i=0;i<Mp->nlman->nsdofal;i++){
01517 norm+=displincr[Mp->nlman->seldofal[i]]*displincr[Mp->nlman->seldofal[i]];
01518 }
01519 break;
01520 }
01521 case nodesdistincr:{
01522 if (Mp->nlman->probdimal==2){
01523 u1=0.0; u2=0.0; v1=0.0; v2=0.0;
01524 if (Mp->nlman->seldofal[0]>0) u1=displincr[Mp->nlman->seldofal[0]-1];
01525 if (Mp->nlman->seldofal[1]>0) v1=displincr[Mp->nlman->seldofal[1]-1];
01526 if (Mp->nlman->seldofal[2]>0) u2=displincr[Mp->nlman->seldofal[2]-1];
01527 if (Mp->nlman->seldofal[3]>0) v2=displincr[Mp->nlman->seldofal[3]-1];
01528
01529 norm=(u2-u1)*Mp->nlman->nxal + (v2-v1)*Mp->nlman->nyal;
01530 }
01531 if (Mp->nlman->probdimal==3){
01532 u1=0.0; u2=0.0; v1=0.0; v2=0.0; w1=0.0; w2=0.0;
01533 if (Mp->nlman->seldofal[0]>0) u1=displincr[Mp->nlman->seldofal[0]-1];
01534 if (Mp->nlman->seldofal[1]>0) v1=displincr[Mp->nlman->seldofal[1]-1];
01535 if (Mp->nlman->seldofal[2]>0) w1=displincr[Mp->nlman->seldofal[2]-1];
01536 if (Mp->nlman->seldofal[3]>0) u2=displincr[Mp->nlman->seldofal[3]-1];
01537 if (Mp->nlman->seldofal[4]>0) v2=displincr[Mp->nlman->seldofal[4]-1];
01538 if (Mp->nlman->seldofal[5]>0) w2=displincr[Mp->nlman->seldofal[5]-1];
01539
01540 norm=(u2-u1)*Mp->nlman->nxal + (v2-v1)*Mp->nlman->nyal + (w2-w1)*Mp->nlman->nzal;
01541 }
01542 break;
01543 }
01544 default:{
01545 fprintf (stderr,"\n\n unknown norm of displacement increment is required in function");
01546 fprintf (stderr,"\n displincr (%s, line %d).\n",__FILE__,__LINE__);
01547 }
01548 }
01549 norm=sqrt (norm);
01550 return norm;
01551 }
01552
01553
01554 void quadeqcoeff (double *ddr,double *v,long n,double ddlambda,double psi,double norfp,double dl,
01555 double &a0,double &a1,double &a2)
01556 {
01557 long i;
01558
01559 switch (Mp->nlman->displnorm){
01560 case alldofs:{
01561 a2=ss(v,v,n)+psi*psi*norfp;
01562 a1=2.0*ss(ddr,v,n)+2.0*ddlambda*psi*psi*norfp;
01563 a0=ss(ddr,ddr,n)+ddlambda*ddlambda*psi*psi*norfp-dl*dl;
01564 break;
01565 }
01566 case seldofs:
01567 case seldofscoord:
01568 case selnodes:{
01569 a0=0.0; a1=0.0; a2=0.0;
01570 for (i=0;i<Mp->nlman->nsdofal;i++){
01571 a0+=ddr[Mp->nlman->seldofal[i]]*ddr[Mp->nlman->seldofal[i]];
01572 a1+=ddr[Mp->nlman->seldofal[i]]*v[Mp->nlman->seldofal[i]]*2.0;
01573 a2+=v[Mp->nlman->seldofal[i]]*v[Mp->nlman->seldofal[i]];
01574 }
01575 a0-=dl*dl;
01576 break;
01577 }
01578 case nodesdistincr:{
01579 break;
01580 }
01581 default:{
01582 fprintf (stderr,"\n\n unknown norm of displacement increment is required in function");
01583 fprintf (stderr,"\n displincr (%s, line %d).\n",__FILE__,__LINE__);
01584 }
01585 }
01586
01587 }
01588
01589 void arclsave (long fi,long n,double blambda,double dl,double *r,double *fp)
01590 {
01591 long i,j;
01592 char *file;
01593 FILE *aux;
01594 file = new char[strlen(Mp->path)+11+1];
01595 sprintf (file,"%sarcl.backup",Mp->path);
01596 aux = fopen (file,"w");
01597
01598
01599
01600
01601 fprintf (aux,"%ld\n",fi);
01602
01603 fprintf (aux,"%ld\n",n);
01604
01605 fprintf (aux,"%e\n",blambda);
01606
01607 fprintf (aux,"%e\n",dl);
01608
01609 fprintf (aux,"%ld\n",Mp->nlman->nsdofal);
01610
01611 for (i=0;i<Mp->nlman->nsdofal;i++){
01612 fprintf (aux,"%ld\n",Mp->nlman->seldofal[i]);
01613 }
01614
01615 for (i=0;i<n;i++){
01616 fprintf (aux,"%e\n",r[i]);
01617 }
01618
01619 for (i=0;i<n;i++){
01620 fprintf (aux,"%e\n",fp[i]);
01621 }
01622
01623 fprintf (aux,"%ld\n",Mm->tnip);
01624
01625 for (i=0;i<Mm->tnip;i++){
01626 fprintf (aux,"%ld %ld\n",Mm->ip[i].ncompstr,Mm->ip[i].ncompother);
01627 for (j=0;j<Mm->ip[i].ncompstr;j++){
01628 fprintf (aux,"%e\n",Mm->ip[i].strain[j]);
01629 }
01630 for (j=0;j<Mm->ip[i].ncompstr;j++){
01631 fprintf (aux,"%e\n",Mm->ip[i].stress[j]);
01632 }
01633 for (j=0;j<Mm->ip[i].ncompother;j++){
01634 fprintf (aux,"%e\n",Mm->ip[i].eqother[j]);
01635 }
01636 }
01637
01638 fclose (aux);
01639 }
01640
01641 void arclopen (long &fi,long &n,double &blambda,double &dl,double *r,double *fp)
01642 {
01643 long i,j;
01644 char *file;
01645 FILE *aux;
01646
01647
01648 file = new char[strlen(Mp->path)+11+1];
01649 sprintf (file,"%sarcl.backup",Mp->path);
01650 aux = fopen (file,"r");
01651
01652
01653 fscanf (aux,"%ld",&fi);
01654
01655 fscanf (aux,"%ld",&n);
01656
01657 fscanf (aux,"%le",&blambda);
01658
01659 fscanf (aux,"%le",&dl);
01660
01661 fscanf (aux,"%ld",&Mp->nlman->nsdofal);
01662
01663 for (i=0;i<Mp->nlman->nsdofal;i++){
01664 fscanf (aux,"%ld",&Mp->nlman->seldofal[i]);
01665 }
01666
01667 for (i=0;i<n;i++){
01668 fscanf (aux,"%le",&r[i]);
01669 }
01670
01671 for (i=0;i<n;i++){
01672 fscanf (aux,"%le",&fp[i]);
01673 }
01674
01675 fscanf (aux,"%ld",&Mm->tnip);
01676
01677 for (i=0;i<Mm->tnip;i++){
01678 fscanf (aux,"%ld %ld",&Mm->ip[i].ncompstr,&Mm->ip[i].ncompother);
01679 for (j=0;j<Mm->ip[i].ncompstr;j++){
01680 fscanf (aux,"%le",&Mm->ip[i].strain[j]);
01681 }
01682 for (j=0;j<Mm->ip[i].ncompstr;j++){
01683 fscanf (aux,"%le",&Mm->ip[i].stress[j]);
01684 }
01685 for (j=0;j<Mm->ip[i].ncompother;j++){
01686 fscanf (aux,"%le",&Mm->ip[i].eqother[j]);
01687 }
01688 }
01689
01690 fclose (aux);
01691 }
01692
01693
01694
01695
01696
01697
01698
01699 void arclsave_adapt (long fi,double lambda,double dlambda,double dl)
01700 {
01701 char *file;
01702 FILE *aux;
01703
01704 file = new char[strlen(Mp->path)+19+1];
01705
01706 sprintf (file,"%sarcl.backup.ada",Mp->path);
01707 aux = fopen (file,"w");
01708
01709
01710 fprintf (aux,"%ld\n",fi);
01711
01712
01713 fprintf (aux,"%.16e\n",lambda);
01714
01715
01716 fprintf (aux,"%.16e\n",dlambda);
01717
01718
01719 fprintf (aux,"%.16e\n",dl);
01720
01721 fclose (aux);
01722 }
01723
01724
01725
01726
01727 void arclopen_adapt (long &fi,double &lambda,double &dlambda,double &dl)
01728 {
01729 char *file;
01730 FILE *aux;
01731
01732 file = new char[strlen(Mp->path)+19+1];
01733
01734 sprintf (file,"%sarcl.backup.ada",Mp->path);
01735 aux = fopen (file,"r");
01736
01737
01738 fscanf (aux,"%ld",&fi);
01739
01740
01741 fscanf (aux,"%le",&lambda);
01742
01743
01744 fscanf (aux,"%le",&dlambda);
01745
01746
01747 fscanf (aux,"%le",&dl);
01748
01749 fclose (aux);
01750 }
01751
01752
01753
01754
01755
01756 void newmeshgen (int w,long i)
01757 {
01758 char procname[255];
01759
01760 if (Mespr==1) fprintf (stdout,"\n\n *** T3D generator - start ***\n");
01761
01762 if (Mespr) sprintf (procname,"/home/dr/Bin/T3d -i %sXmesht3d.in -o %smesh.out.%0*ld -m %smesh.bgm.%0*ld -d 100 -p 8 -r 1 -k %ld"
01763 ,Mp->path,Mp->path,w,i,Mp->path,w,i-1,Mt->give_degree(0));
01764 else sprintf (procname,"/home/dr/Bin/T3d -i %sXmesht3d.in -o %smesh.out.%0*ld -m %smesh.bgm.%0*ld -d 100 -p 8 -r 1 -k %ld > /dev/null 2>&1"
01765 ,Mp->path,Mp->path,w,i,Mp->path,w,i-1,Mt->give_degree(0));
01766
01767 if (Mespr==1) fprintf (stdout,"\n\n *** T3D generator - end *****\n");
01768
01769
01770 system (procname);
01771
01772
01773 sprintf (procname,"sed -e s/out.i/out.%0*ld/g -e s/dat.i/dat.%d/g %sprep.pre > %sprep.%0*ld",w,i,0,Mp->path,Mp->path,w,i);
01774 system (procname);
01775
01776 if (Mespr)
01777 sprintf (procname,"./mechprep %sprep.%0*ld %ssifel.in.%0*ld",Mp->path,w,i,Mp->path,w,i);
01778 else
01779 sprintf (procname,"./mechprep %sprep.%0*ld %ssifel.in.%0*ld > /dev/null 2>&1",Mp->path,w,i,Mp->path,w,i);
01780
01781 system (procname);
01782
01783
01784 sprintf (procname,"rm -f %sprep.%0*ld",Mp->path,w,i);
01785 system (procname);
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 long arclengthadapt (long lcid,long adaptcontrol)
01818 {
01819 long i,j,k,n,ni,ini,stop,modif,li,newmesh, numr;
01820 double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr;
01821 double lambdao,ss1,ss2,ss3,ss4,ss5;
01822 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi;
01823
01824
01825
01826
01827 n = Ndofm;
01828
01829 ni = Mp->nlman->nial;
01830
01831 ini = Mp->nlman->niilal;
01832
01833 zero = Mp->zero;
01834
01835 ierr = Mp->nlman->erral;
01836
01837 dl = Mp->nlman->dlal;
01838
01839 dlmax = Mp->nlman->dlmaxal;
01840
01841 dlmin = Mp->nlman->dlminal;
01842
01843 psi = Mp->nlman->psial;
01844
01845 dlambda=0.0;
01846
01847
01848 r = new double [n];
01849 ddr = new double [n];
01850 u = new double [n];
01851 v = new double [n];
01852 f = new double [n];
01853 fa = new double [n];
01854 fi = new double [n];
01855
01856
01857
01858 ra = Lsrs->give_lhs (lcid);
01859 fp = Lsrs->give_rhs (lcid*2);
01860 fc = Lsrs->give_rhs (lcid*2+1);
01861
01862 if (Mp->nlman->hdbackupal == 2){
01863 arclopen (li,n,lambda,dl,ra,fp);
01864 }
01865 else if (adaptcontrol >= 2){
01866 arclopen_adapt (li,lambda,blambda,dl);
01867 }
01868 else{
01869 lambda=0.0;
01870 lambdao=0.0;
01871 li=0;
01872 }
01873
01874
01875 norfp = ss(fp,fp,n);
01876 modif=0;
01877
01878
01879
01880 for (j=0;j<n;j++){
01881 fa[j]=fc[j]+(lambda+dlambda)*fp[j];
01882 f[j]=fp[j];
01883 }
01884 if (li)
01885 print_init(-1, "at");
01886 else
01887 {
01888 print_init(-1, "wt");
01889 print_step(lcid, li, lambda, fa);
01890 }
01891
01892
01893
01894 Mp->phase=1;
01895
01896 stiffness_matrix (lcid);
01897
01898 Smat->solve_system (ra,f);
01899
01900 Mm->computestresses (0);
01901
01902
01903 Mp->phase=2;
01904
01905
01906
01907
01908
01909 for (i=li;i<ni;i++){
01910 stop=1;
01911
01912 copyv (ra, r, n);
01913
01914 blambda=lambda;
01915
01916
01917
01918 if ((Mp->stmat==tangent_stiff) || (i == li))
01919 stiffness_matrix (lcid);
01920
01921
01922 copyv(fp, f, n);
01923
01924
01925 Smat->solve_system (v,f);
01926
01927
01928 norv = displincr (v,n);
01929
01930
01931 dlambda = dl/sqrt(norv+psi*psi*norfp);
01932
01933 for (j=0;j<n;j++){
01934 ddr[j]=dlambda*v[j];
01935 ra[j]+=ddr[j];
01936 fa[j]=fc[j]+(lambda+dlambda)*fp[j];
01937 }
01938 ddlambda=dlambda;
01939
01940
01941
01942 internal_forces (lcid,fi);
01943 subv(fa, fi, f, n);
01944 norf=sizev(f,n);
01945 norfa = sizev(fa, n);
01946 norf /= norfa;
01947
01948
01949
01950 if (Mespr==1) fprintf (stdout,"\n increment %ld norv %e norf %e",i,norv,norf);
01951
01952 if (norf<ierr){
01953
01954
01955
01956 modif++;
01957
01958 if (modif>1){
01959
01960 dl*=2.0;
01961 if (dl>dlmax)
01962 dl=dlmax;
01963 modif=0;
01964 if (Mespr==1){
01965 fprintf (stdout,"\n arc-length must be modified (dl increase) dl=%e",dl);
01966
01967 }
01968 }
01969 lambda+=dlambda;
01970 Mm->updateipval ();
01971 print_step(lcid, i+1, lambda, fa);
01972 print_flush();
01973 }
01974 else{
01975
01976
01977
01978 stop=0;
01979 for (j=0;j<ini;j++){
01980
01981 if (Mp->stmat==tangent_stiff)
01982 stiffness_matrix (lcid);
01983
01984
01985 Smat->solve_system (u,f);
01986
01987 copyv(ddr, f, n);
01988 addv(ddr, u, n);
01989
01990 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2);
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008 numr = solv_polynom_2(a2, a1, a0, l1, l2);
02009 switch (numr)
02010 {
02011 case -1:
02012 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength");
02013 break;
02014 case 0:
02015 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength");
02016 break;
02017 case 1:
02018 dlambda = l1;
02019 break;
02020 default:
02021 break;
02022 }
02023 ss1=0.0; ss2=0.0;
02024 ss3=0.0; ss4=0.0; ss5=0.0;
02025 for (k=0;k<n;k++){
02026 ss1+=(ddr[k]+l1*v[k])*f[k];
02027 ss2+=(ddr[k]+l2*v[k])*f[k];
02028 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]);
02029 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]);
02030 ss5+=f[k]*f[k];
02031 }
02032
02033
02034
02035
02036
02037
02038
02039 if (fabs(l1)>fabs(l2)) dlambda=l2;
02040 else dlambda=l1;
02041
02042
02043
02044
02045
02046 for (k=0;k<n;k++){
02047 ddr[k]+=dlambda*v[k];
02048 ra[k]+=u[k]+dlambda*v[k];
02049 fa[k]+=dlambda*fp[k];
02050 }
02051 ddlambda+=dlambda;
02052
02053
02054
02055
02056
02057 internal_forces (lcid,fi);
02058 subv(fa, fi, f, n);
02059 norf=sizev(f,n);
02060 norfa = sizev(fa, n);
02061 norf /= norfa;
02062
02063 if (Mespr==1){
02064 fprintf (stdout,"\n norf=%e ierr=%e",norf,ierr);
02065
02066 }
02067 if (norf<ierr){
02068 lambda+=ddlambda;
02069 Mm->updateipval ();
02070 stop=1;
02071 print_step(lcid, i+1, lambda, fa);
02072 print_flush();
02073 break;
02074 }
02075 }
02076 modif=0;
02077 if (stop==0){
02078
02079 dl/=2.0;
02080 if (dl<dlmin){
02081 dl=dlmin;
02082 break;
02083 }
02084 if (Mespr==1){
02085 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e",dl);
02086
02087 }
02088
02089 copyv(r, ra, n);
02090
02091 lambda=blambda;
02092 }
02093 }
02094
02095 fprintf (stdout,"\n increment %ld total lambda %e",i,lambda);
02096
02097 if (stop==0)
02098 continue;
02099
02100 if (Mp->adaptflag & 16 && i%10 && Mm->plast)
02101 continue;
02102
02103 newmesh = 0;
02104 if (Mp->adaptivity && Mp->nlman->nial > (i+1) && !( (Mp->adaptflag & 16) && (Mp->nlman->nial==100 || Mp->nlman->nial==1))){
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114 long width = 1 + (long)log10((double)Mp->nlman->nial);
02115 newmesh = Ada->run (1,width,i);
02116 }
02117
02118 if (newmesh){
02119 i++;
02120 break;
02121 }
02122
02123
02124 }
02125
02126
02127
02128
02129
02130 if (Mp->nlman->hdbackupal==1)
02131 arclsave (i,n,lambda,dl,ra,fp);
02132
02133
02134 if (Mp->adaptivity)
02135 if ((Mp->adaptflag & 16) && (Mp->nlman->nial==100 || Mp->nlman->nial==1))
02136 Ada->run (0,0,0);
02137 else
02138 if (i < ni)
02139 arclsave_adapt (i,lambda,0,dl);
02140
02141
02142
02143 print_close();
02144 delete [] r;
02145 delete [] fi;
02146 delete [] fa;
02147 delete [] f;
02148 delete [] v;
02149 delete [] u;
02150 delete [] ddr;
02151
02152 return (i);
02153 }