00001 #include "slopesol.h"
00002 #include "arclength.h"
00003 #include "global.h"
00004 #include "globmat.h"
00005 #include "mechprint.h"
00006 #include "loadcase.h"
00007 #include "dloadcase.h"
00008 #include "intpoints.h"
00009 #include "vector.h"
00010 #include "matrix.h"
00011 #include "mathem.h"
00012 #include "gmatrix.h"
00013 #include "tensor.h"
00014 #include "vecttens.h"
00015 #include "nssolver.h"
00016
00017 #include <math.h>
00018 #include <string.h>
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 double arclengthrv1 (long lcid, double rlambda, double rlerr)
00038 {
00039 long i,j,k,n,ni,ini,stop,modif,li,numr;
00040 double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr;
00041 double lambdao,ss1,ss2,ss3,ss4,ss5;
00042 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi;
00043 const char *mode;
00044 long back_dl = 1;
00045 long check_rv = 1;
00046 long nodiv_dl = 0;
00047 double bdl;
00048
00049
00050
00051
00052 n = Ndofm;
00053
00054 ni = Mp->nlman->nial;
00055
00056 ini = Mp->nlman->niilal;
00057
00058 zero = Mp->zero;
00059
00060 ierr = Mp->nlman->erral;
00061
00062 dl = Mp->nlman->dlal;
00063
00064 dlmax = Mp->nlman->dlmaxal;
00065
00066 dlmin = Mp->nlman->dlminal;
00067
00068 psi = Mp->nlman->psial;
00069
00070
00071 r = new double [n];
00072 ddr = new double [n];
00073 u = new double [n];
00074 v = new double [n];
00075 f = new double [n];
00076 fa = new double [n];
00077 fi = new double [n];
00078
00079
00080 ra = Lsrs->give_lhs (lcid);
00081 fp = Lsrs->give_rhs (lcid*2);
00082 fc = Lsrs->give_rhs (lcid*2+1);
00083
00084 if (Mp->nlman->hdbackupal==2){
00085 arclopen (li,n,lambda,dl,ra,fp);
00086 }
00087 else{
00088 lambda=0.0;
00089 lambdao=0.0;
00090 li=0;
00091 }
00092
00093 if (li == 0)
00094 mode = "wt";
00095 else
00096 mode = "at";
00097
00098
00099 for (j=0;j<n;j++)
00100 {
00101 if (lambda <= rlambda)
00102 fa[j]=lambda*(fc[j]+fp[j]);
00103 else
00104 fa[j]=rlambda*fc[j]+lambda*fp[j];
00105 }
00106
00107 norfp = ss(fp,fp,n);
00108 modif=0;
00109
00110
00111
00112
00113
00114 for (i=li;i<ni;i++) {
00115
00116
00117 copyv (ra, r, n);
00118
00119 blambda=lambda;
00120 if (back_dl)
00121 bdl = dl;
00122
00123 fprintf (stdout,"\n arc-length: increment %ld lambda %e dl %e",i,lambda,dl);
00124
00125
00126
00127
00128
00129
00130 if ((Mp->nlman->stmat==tangent_stiff) || (i == li))
00131 stiffness_matrix (lcid);
00132
00133
00134 if (check_rv)
00135 addv(fc, fp, f, n);
00136 else
00137 copyv (fp, f, n);
00138
00139
00140 Mp->ssle->solve_system (Gtm,Smat,v,f,Out);
00141
00142 norv = displincr (v,v,n,lcid);
00143
00144 dlambda = dl/sqrt(norv+psi*psi*norfp);
00145
00146 for (j=0;j<n;j++){
00147 ddr[j]=dlambda*v[j];
00148 ra[j]+=ddr[j];
00149 if (check_rv)
00150 fa[j]=(lambda+dlambda)*(fc[j]+fp[j]);
00151 else
00152 fa[j]=rlambda*fc[j]+(lambda+dlambda)*fp[j];
00153 }
00154 ddlambda=dlambda;
00155
00156
00157 internal_forces (lcid,fi);
00158 subv(fa, fi, f, n);
00159 norf=normv(f,n);
00160 norfa = normv(fa, n);
00161 norf /= norfa;
00162
00163 if (Mespr==1) fprintf (stdout,"\n %e %e norf %e",lambda,dl,norf);
00164
00165 if (norf<ierr){
00166
00167
00168
00169 modif++;
00170
00171 if (modif>1) {
00172
00173 dl*=2.0;
00174 if (dl>dlmax)
00175 dl=dlmax;
00176 modif=0;
00177 if (Mespr==1){
00178 fprintf (stdout,"\n arc length must be modified (dl increase) dl=%e norf=%e",dl,norf);
00179
00180 }
00181 }
00182 lambda+=dlambda;
00183 if (check_rv && ((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
00184 {
00185 if (back_dl)
00186 {
00187 bdl = dl;
00188 back_dl = 0;
00189 }
00190
00191 dl/=2.0;
00192 if (dl<dlmin)
00193 {
00194 dl=dlmin;
00195 if (nodiv_dl)
00196 break;
00197 else
00198 nodiv_dl = 1;
00199 }
00200 if (Mespr==1)
00201 {
00202 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf);
00203
00204 }
00205
00206 for (j=0;j<n;j++)
00207 ra[j]=r[j];
00208
00209 lambda=blambda;
00210 }
00211 if (fabs(lambda-rlambda) <= rlerr)
00212 {
00213 dl = bdl;
00214 check_rv = 0;
00215 }
00216 Mm->updateipval ();
00217 compute_req_val (lcid);
00218 print_step(lcid, i, lambda, fa);
00219 print_flush();
00220 }
00221 else{
00222
00223
00224
00225 stop=0;
00226 for (j=0;j<ini;j++){
00227
00228
00229 Mp->ssle->solve_system (Gtm,Smat,u,f,Out);
00230 copyv (ddr, f, n);
00231 addv(ddr, u, n);
00232
00233 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2,lcid);
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 numr = solv_polynom_2(a2, a1, a0, l1, l2);
00252 switch (numr)
00253 {
00254 case -1:
00255 fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength");
00256 break;
00257 case 0:
00258 fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength");
00259 break;
00260 case 1:
00261 dlambda = l1;
00262 break;
00263 default:
00264 break;
00265 }
00266 ss1=0.0; ss2=0.0;
00267 ss3=0.0; ss4=0.0; ss5=0.0;
00268 for (k=0;k<n;k++){
00269 ss1+=(ddr[k]+l1*v[k])*f[k];
00270 ss2+=(ddr[k]+l2*v[k])*f[k];
00271 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]);
00272 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]);
00273 ss5+=f[k]*f[k];
00274 }
00275
00276 if (ss1/ss3/ss5>ss2/ss4/ss5)
00277 dlambda=l1;
00278 else
00279 dlambda=l2;
00280
00281 for (k=0;k<n;k++){
00282 ddr[k]+=dlambda*v[k];
00283 ra[k]+=u[k]+dlambda*v[k];
00284 if (check_rv)
00285 fa[k]+=dlambda*(fp[k]+fc[k]);
00286 else
00287 fa[k]+=dlambda*fp[k];
00288 }
00289 ddlambda+=dlambda;
00290
00291
00292 internal_forces (lcid,fi);
00293 subv(fa, fi, f, n);
00294 norf=normv(f,n);
00295 norfa = normv(fa, n);
00296 norf /= norfa;
00297
00298 if (Mespr==1){
00299 fprintf (stdout,"\n increment %ld inner loop %ld norf=%e",i,j,norf);
00300
00301 }
00302 if (norf<ierr){
00303 lambda+=ddlambda;
00304 compute_req_val (lcid);
00305 print_step(lcid, i, lambda, fa);
00306 print_flush();
00307 stop=1;
00308 if (check_rv && ((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
00309 stop=0;
00310 if (check_rv && (fabs(lambda-rlambda) <= rlerr))
00311 {
00312 dl = bdl;
00313 check_rv = 0;
00314 }
00315 if (stop > 0)
00316 Mm->updateipval();
00317 break;
00318 }
00319 }
00320
00321 modif=0;
00322 if (stop==0){
00323
00324 dl/=2.0;
00325 if (dl<dlmin){
00326 dl=dlmin;
00327 if (nodiv_dl)
00328 break;
00329 else
00330 nodiv_dl = 1;
00331 }
00332 if (Mespr==1){
00333 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf);
00334
00335 }
00336
00337 copyv (r, ra, n);
00338
00339 lambda=blambda;
00340 }
00341 }
00342 }
00343
00344
00345
00346
00347
00348 if (Mp->nlman->hdbackupal==1)
00349 arclsave (i,n,lambda,dl,ra,fp);
00350
00351
00352 delete [] r;
00353 delete [] fi;
00354 delete [] fa;
00355 delete [] f;
00356 delete [] v;
00357 delete [] u;
00358 delete [] ddr;
00359
00360 return (i);
00361 }
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378 void newton_raphsonrv1 (long lcid, double rlambda, double rlerr)
00379
00380
00381
00382
00383
00384 {
00385 long i,j,k,n,ni,ini;
00386 double lambda,blambda,dlambda,dlambdamin,dlambdamax,zero,err,norf,norfa;
00387 double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp;
00388 long back_incr = 1;
00389 long check_rv = 1;
00390 double binc;
00391 matrix lsm_a(3,3);
00392 vector lsm_r(3), lsm_l(3);
00393 Mp->nlman->divc_step = new double[10];
00394 Mp->nlman->divc_err = new double[10];
00395
00396
00397
00398 n = Ndofm;
00399
00400 ni = Mp->nlman->ninr;
00401
00402 ini = Mp->nlman->niilnr;
00403
00404 zero = Mp->zero;
00405
00406 err = Mp->nlman->errnr;
00407
00408 dlambda=Mp->nlman->incrnr;
00409
00410 dlambdamin=Mp->nlman->minincrnr;
00411
00412 dlambdamax=Mp->nlman->maxincrnr;
00413
00414 rb = new double [n];
00415 f = new double [n];
00416 fb = new double [n];
00417 fi = new double [n];
00418 dr = new double [n];
00419 memset (rb,0,n*sizeof(double));
00420 memset (f,0,n*sizeof(double));
00421 memset (fi,0,n*sizeof(double));
00422 memset (dr,0,n*sizeof(double));
00423
00424
00425 r = Lsrs->give_lhs (lcid);
00426 fp = Lsrs->give_rhs (lcid*2);
00427 fc = Lsrs->give_rhs (lcid*2+1);
00428 lambda=0.0;
00429
00430
00431
00432
00433 print_init(-1, "wt");
00434 print_step(lcid, 0, 0.0, f);
00435 print_flush();
00436 for (i=0;i<ni;i++){
00437
00438
00439
00440 Mp->strainstate=0;
00441
00442
00443 Mp->stressstate=0;
00444
00445
00446 Mp->otherstate=0;
00447
00448
00449
00450 for (j=0;j<n;j++){
00451 rb[j]=r[j];
00452 }
00453
00454 blambda=lambda;
00455 if (back_incr)
00456 binc = dlambda;
00457
00458 fprintf (stdout,"\n increment %ld lambda %e, dlambda %e",i,lambda,dlambda);
00459
00460
00461 for (j=0;j<n;j++){
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472 if (check_rv)
00473 {
00474 f[j]=lambda*(fc[j]);
00475 fb[j]=dlambda*(fc[j]);
00476 }
00477 else
00478 {
00479 f[j]=rlambda*fc[j]+(lambda-rlambda)*fp[j];
00480 fb[j]=dlambda*fp[j];
00481 }
00482 }
00483
00484
00485 if ((Mp->nlman->stmat==tangent_stiff) || (i == 0))
00486 stiffness_matrix (lcid);
00487
00488
00489
00490 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00491
00492 for (j=0;j<n;j++){
00493 r[j]+=dr[j];
00494 }
00495
00496
00497 internal_forces (lcid,fi);
00498
00499
00500 for (j=0;j<n;j++){
00501
00502
00503
00504
00505 if (check_rv)
00506 fb[j]=f[j]+dlambda*(fc[j]);
00507 else
00508 fb[j]=f[j]+dlambda*fp[j];
00509 }
00510 norfa=normv(fb,n);
00511 for (j=0;j<n;j++){
00512 fb[j] -= fi[j];
00513 }
00514
00515 norf=normv(fb,n);
00516 norf /= norfa;
00517
00518 if (norf<err){
00519 lambda+=dlambda;
00520 if (check_rv && ((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
00521 {
00522 if (back_incr)
00523 {
00524 binc = dlambda;
00525 back_incr = 0;
00526 }
00527
00528 dlambda/=2.0;
00529 if (dlambda<=dlambdamin)
00530 {
00531 if (dlambda == dlambdamin)
00532 break;
00533 dlambda=dlambdamin;
00534 }
00535 if (Mespr==1)
00536 fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf=%e",dlambda,norf);
00537
00538 for (j=0;j<n;j++)
00539 r[j]=rb[j];
00540
00541 lambda=blambda;
00542 }
00543 if (fabs(lambda-rlambda) <= rlerr)
00544 {
00545 dlambda = binc;
00546 check_rv = 0;
00547 }
00548 Mm->updateipval ();
00549 compute_req_val (lcid);
00550 print_step(lcid, i+1, lambda, f);
00551 print_flush();
00552 continue;
00553 }
00554
00555
00556 fillm(0.0, lsm_a);
00557 fillv(0.0, lsm_r);
00558 for (j=0;j<ini;j++) {
00559
00560
00561
00562 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00563
00564 for (k=0;k<n;k++){
00565 r[k]+=dr[k];
00566 }
00567
00568
00569 internal_forces (lcid,fi);
00570
00571
00572 for (k=0;k<n;k++){
00573
00574
00575
00576
00577 if (check_rv)
00578 fb[k]=f[k]+dlambda*(fc[k])-fi[k];
00579 else
00580 fb[k]=f[k]+dlambda*fp[k]-fi[k];
00581 }
00582
00583
00584 norf=normv(fb,n);
00585 norf /= norfa;
00586
00587 fprintf (stdout,"\n increment %ld inner loop j %ld norf=%e dlambda %e",i,j,norf,dlambda);
00588
00589 if (norf<err) break;
00590
00591
00592 if (j > 10)
00593 {
00594 if (lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), Mp->nlman->divc_step[j%10],
00595 Mp->nlman->divc_err[j%10], norf, Mp->zero,1) > 0.0)
00596 {
00597 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
00598 break;
00599 }
00600 Mp->nlman->divc_step[j%10] = double(j+1);
00601 Mp->nlman->divc_err[j%10] = err;
00602 }
00603 else
00604 {
00605 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,0);
00606 Mp->nlman->divc_step[j%10] = double(j+1);
00607 Mp->nlman->divc_err[j%10] = err;
00608 }
00609 }
00610
00611 if (j==ini || norf>err){
00612 if (dlambda == dlambdamin)
00613 {
00614 fprintf (stdout,"\n increment of lambda cannot be decreased, iteration will be stopped\n");
00615 break;
00616 }
00617
00618 dlambda/=2.0;
00619 if (dlambda<dlambdamin) dlambda=dlambdamin;
00620
00621 if (Mespr==1)
00622 fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf=%e",dlambda,norf);
00623
00624 for (j=0;j<n;j++){
00625 r[j]=rb[j];
00626 }
00627 lambda=blambda;
00628 }
00629 else{
00630 lambda+=dlambda;
00631 if (check_rv && ((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
00632 {
00633 if (back_incr)
00634 {
00635 binc = dlambda;
00636 back_incr = 0;
00637 }
00638
00639 dlambda/=2.0;
00640 if (dlambda<=dlambdamin)
00641 {
00642 if (dlambda == dlambdamin)
00643 break;
00644 dlambda=dlambdamin;
00645 }
00646 if (Mespr==1)
00647 fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf=%e",dlambda,norf);
00648
00649 for (j=0;j<n;j++)
00650 r[j]=rb[j];
00651
00652 lambda=blambda;
00653 }
00654 if (fabs(lambda-rlambda) <= rlerr)
00655 {
00656 if (back_incr == 0)
00657 dlambda = binc;
00658 check_rv = 0;
00659 }
00660 Mm->updateipval ();
00661 compute_req_val (lcid);
00662 print_step(lcid, i+1, lambda, f);
00663 print_flush();
00664 }
00665
00666 }
00667
00668 delete [] dr; delete [] fi; delete [] fb; delete [] f;
00669 print_close();
00670 }